Numerical Modeling on Crack Propagation Based on a Multi-Grid Bond-Based Dual-Horizon Peridynamics

Abstract

Peridynamics (PD) is a novel nonlocal theory of continuum mechanics capable of describing crack formation and propagation without defining any fracture rules in advance. In this study, a multi-grid bond-based dual-horizon peridynamics (DH-PD) model is presented, which includes varying horizon sizes and can avoid spurious wave reflections. This model incorporates the volume correction, surface correction, and a technique of nonuniformity discretization to improve calculation accuracy and efficiency. Two benchmark problems are simulated to verify the reliability of the proposed model with the effect of the volume correction and surface correction on the computational accuracy confirmed. Two numerical examples, the fracture of an L-shaped concrete specimen and the mixed damage of a double-edged notched specimen, are simulated and analyzed. The simulation results are compared against experimental data, the numerical solution of a traditional PD model, and the output from a finite element model. The comparisons verify the calculation accuracy of the corrected DH-PD model and its advantages over some other models like the traditional PD model.

1. Introduction

Material fracture is a classical problem in solid mechanics [1]. In contrast to experimental research, numerical simulations provide a more efficient approach to understanding material fractures since they eliminate the limitations of complex specimen preparation, test environment conditions, and loading systems. The finite element method (FEM) is one of the commonly used numerical methods for solving fracture problems [2,3]. However, the crack-induced spatial discontinuities are not compatible with the partial differential equations in classical FEM models [4]. To deal with this problem, some remeshing techniques were incorporated into FEM models, which were usually cumbersome and time consuming [5,6]. Recently, the extended finite element method (XFEM) was proposed to solve the intermittent problems of cracks, holes, and inclusions by introducing additional nodal degrees of freedom and local strengthening functions [7,8]. Besides, the element erosion technique [9,10] and the phase field method [11,12] were proposed to simulate the problems of crack propagation with different levels of success. However, these methods are not applicable to simulate complex multi-crack extended penetration problems in three dimensions [13].

In recent years, meshless methods have been developed and applied in various engineering fields. They can simulate complex discontinuity without the remeshing process. For example, Cheng and Li [14] proposed a complex variable meshless method for fracture problems. Gao and Cheng [15] developed a meshless manifold method based on the complex variable moving least-squares approximation and the finite cover theory for fracture problems. Khazal et al. [16] presented an extended element free Galerkin method (XEFGM) for fracture analysis of functionally graded materials. In addition, the discrete element method (DEM) is suitable for simulating crack formation and propagation based on the discontinuity assumption. For instance, Kou et al. [17] adopted the DEM to investigate the dynamic crack propagation and branching in brittle solids. Owayo et al. [18] investigated the formation and propagation of cracks in cement-based materials via DEM simulations. Although some promising achievements have been obtained through the above researches, DEM has limitations when describing the mechanical behavior of the continuous material.

More recently, Silling et al. [19,20] proposed a new nonlocal theory of continuum mechanics called peridynamics (PD) for solving material fracture problems. It describes such a phenomenon by means of the failure of bonds linking material points. This method replaces the differential form of traditional balance equations with an integral form to avoid the singularity problem arising from the derivations at discontinuities. Therefore, PD is appropriate to describe the crack propagation in a unified framework. Owing to these advantages, it has been successfully applied to crack propagation problems [21,22], composite material damage [23,24], fluid mechanics and acoustics [25], and heat conduction analysis [26,27] etc. There are several kinds of peridynamics, such as bond-based PD [19], state-based PD [28], and element-based PD [29] etc., and they can be easily coupled with other numerical methods to deal with the interaction between different materials, such as PD-FEM [30], PD-SPH [31], and PD-DEM [32].

In traditional PD, the problem domain should be discretized into a series of particles with the same size and uniform distribution, thus limiting its application to complex problems such as multiscale analysis, multi-material problems, and adaptive analysis. Dipasquale et al. [33] proposed an adaptive refinement-coarsening method and applied it to a one-dimensional elastic wave propagation simulation. However, the problem of spurious stress waves was not completely solved in this model. Bobaru et al. [34] introduced adaptive refinement algorithms for peridynamics to deal with static and dynamic elasticity problems in 1D. Shojaei et al. [35] proposed an adaptive refinement strategy to use a variable grid size in a peridynamic model, which was verified to be effective and efficient. This approach was then successfully applied to simulate the fracture of a ceramic disk under central thermal shock [36]. Ren et al. [37,38] proposed a dual-horizon peridynamics (DH-PD) model incorporating non-uniform particle distribution and variable horizons to completely solve the problem of spurious stress waves. Wang et al. [39] derived state-based DH-PD governing equations using the Euler-Lagrange description, and they analyzed the dispersion characteristics of the model.

In this study, a particle volume correction and a surface correction are introduced in the DH-PD model to improve the calculation accuracy. In addition, the computational efficiency of this model is increased by incorporating a multi-grid technique. Two benchmark problems are simulated to verify the calculation accuracy of this model. Finally, this model is applied to analyze two numerical examples to validate the model’s reliability and applicability.

2. Peridynamics Theory

2.1. Traditional Peridynamics Theory

In traditional PD, the object occupying the spatial domain R is discretized into a series of closely arranged material points, as shown in Figure 1. The circular (spherical) region around the material point, x, with itself as the center point and δ as the radius is called the horizon, H_δ. The interaction force between the material point, x, with any other neighboring point, xʹ, in the horizon (xʹ ∈ H_δ) is called a pairwise force (bond). At any moment, t, the governing equation for any material point, x, is as follows [19]:

$$\rho \frac{{\partial }^2\mathbf{u}(\mathbf{x},t)}{\partial t^2}={\displaystyle {\int }_{{\mathit{H}}_{\mathbf{x}}}\mathbf{f}}(\mathsf{\eta },\mathsf{\xi })\mathrm{d}{V}_{\mathbf{x}\text{'}}+b(\mathbf{x},t)\hspace{1em}\forall \mathsf{\eta },\mathsf{\xi }.$$

(1)

In Equation (1), f is the pairwise bond force between the material point x and its neighboring point xʹ, ρ is the density of the material, b is the body force vector, and u is the displacement vector. The integration domain, H_x, is the horizon of material point, x, in the reference configuration, mathematically, as ${\mathit{H}}_{\mathbf{x}}=\left\{\mathbf{x}\text{'}\in R\left| ‖ \mathbf{x}\text{'}-\mathbf{x} ‖\le \mathsf{\delta }\right.\right\}$. Additionally, δ is the radius of the nonlocal horizon of the material point, which is constant for each material point in the traditional PD theory. Finally, ξ = xʹ − x and η = uʹ − u denote the relative position vector and the relative displacement vector with neighboring particles in the reference configuration, respectively.

The key to building a PD model is to find an appropriate pairwise force function, f. Theoretically, any function that satisfies the linear admissibility condition and the angular admissibility condition can be a pairwise force function. In a bond-based PD model, the direction of the pairwise force vector is parallel to the displacement vector, ξ + η. For microelastic materials, the pairwise force function in PD model is given by

$$\mathbf{f}(\mathsf{\eta },\mathsf{\xi })=\mu c(\mathsf{\delta })s(\mathsf{\eta },\mathsf{\xi })\frac{\mathsf{\eta }+\mathsf{\xi }}{\Vert \mathsf{\eta }+\mathsf{\xi }\Vert }.$$

(2)

In Equation (2), c(δ) denotes the microelastic modulus, which can be derived from the strain energy equivalence principle, and expressed in terms of E and ν:

$$c\left(\mathsf{\delta }\right)=\left\{\begin{array}{cc}\frac{2E}{A{\mathsf{\delta }}^2}& 1\mathrm{D},\\ \frac{3E}{\pi {\mathsf{\delta }}^3\left(1-\nu \right)}& \mathrm{plane}\mathrm{stress},\\ \frac{3E}{\pi {\mathsf{\delta }}^3\left(1+\nu \right)\left(1-2\nu \right)}& \mathrm{plane}\mathrm{strain},\\ \frac{3E}{\pi {\mathsf{\delta }}^4\left(1-2\nu \right)}& 3\mathrm{D}.\end{array}\right.$$

(3)

Additionally, s represents the bond stretch ratio between two material points, as defined by:

$$s(\mathsf{\eta },\mathsf{\xi })=\frac{\Vert \mathsf{\eta }+\mathsf{\xi }\Vert -\Vert \mathsf{\xi }\Vert }{\Vert \mathsf{\xi }\Vert }.$$

(4)

When s is greater than the critical elongation ratio, s₀, of the material, the bond between these two points is considered to be broken, and the corresponding pairwise force disappears. The critical elongation ratio, s₀, can be determined by the relationship given by Ren et al. [37]:

$$s_0\left(\mathsf{\delta }\right)=\left\{\begin{array}{cc}\sqrt{4\pi G_c/(9E\mathsf{\delta })}& \mathrm{plane}\mathrm{stress},\\ \sqrt{5\pi G_c/(12E\mathsf{\delta })}& \mathrm{plane}\mathrm{strain},\\ \sqrt{5G_c/(6E\mathsf{\delta })}& 3\mathrm{D},\end{array}\right.$$

(5)

where G_c is the critical energy release rate of the material.

The variable μ in Equation (2) is a historical parameter describing whether the bond is broken. In this case, μ is defined as:

$$\mu (t,\mathsf{\eta },\mathsf{\xi })=\left\{\begin{array}{ll}1,\hfill & s<s_0,\hfill \\ 0,\hfill & \mathrm{other}.\hfill \end{array}\right.$$

(6)

According to the above criterion for bond breakage, the degree of damage at the local material point, x, can be expressed as:

$$\phi (\mathbf{x},t)=1-\frac{{\displaystyle {\int }_{{\mathit{H}}_{\mathbf{x}}}\mu }(t,\mathsf{\eta },\mathsf{\xi })\mathrm{d}{V}_{\mathbf{x}\text{'}}}{{\displaystyle {\int }_{{\mathit{H}}_{\mathbf{x}}}\mathrm{d}}{V}_{\mathbf{x}\text{'}}}.$$

(7)

2.2. Dual-Horizon Peridynamics Theory

In traditional PD theory, the force acting on the particle x is the sum of the pairwise forces from all the neighboring particles in the horizon. Therefore, the horizon must be constant for all particles to ensure that the force always appears in pairs. If the horizons are variable, as shown in Figure 2, particle x is subjected to the force f(η, ξ) from the particle xʹ in its horizon. However, xʹ is not subjected to the force from x because its horizon does not include particle x. In this case, the existence of such an unbalanced force destroys the balance of linear momentum and angular momentum and results in ghost forces and spurious wave reflections during the simulation.

Recently, Ren et al. [37] proposed the concept of dual-horizon, which can solve the spurious wave reflection problem. The dual-horizon is a union of points, the horizons of which include x. It is expressed, mathematically, as ${\mathit{H}}_{\mathbf{x}}^{\text{'}}=\left\{\mathbf{x}\text{'}\Vert \mathbf{x}\in {\mathit{H}}_{\mathbf{x}\text{'}},\mathbf{x}\text{'}\in R\right\}$. The superscript prime of ${\mathit{H}}_{\mathbf{x}}^{\text{'}}$ denotes a dual-horizon, and the subscript x denotes the object particle. Therefore, dual-horizon can be interpreted as a set of horizons that belong to particles that can include x in their horizons. In the DH-PD theory, a particle is subjected to two types of forces. The first one is the force from a neighboring particle in the horizon. The other one is the force from the particle in the dual-horizon. Furthermore, these two different types of forces are independent of each other, and the bond and dual bond can break independently in the fracture models.

In Figure 2, H_x is the horizon of particle x, where the bond xxʹ will exert a direct force f_xxʹ on x, and xʹ will be subjected to the reaction force − f_xxʹ according to Newton’s third law. Likewise, the dual bond xʹx will exert a direct force f_xʹx on xʹ, and x will be subject to the reaction force − f_xʹx. The time-dependent governing equations of DH-PD [40] are given as follows:

$$\rho \frac{{\partial }^2\mathbf{u}(\mathbf{x},t)}{\partial t^2}={\displaystyle {\int }_{{\mathit{H}}_{\mathbf{x}}}{\mathbf{f}}_{\mathbf{x}\mathbf{x}\text{'}}}(\mathsf{\eta },\mathsf{\xi })\mathrm{d}{V}_{\mathbf{x}\text{'}}-{\displaystyle {\int }_{{\mathit{H}}_{\mathbf{x}}^{\text{'}}}{\mathbf{f}}_{\mathbf{x}\text{'}\mathbf{x}}}(-\mathsf{\eta },-\mathsf{\xi })\mathrm{d}{V}_{\mathbf{x}\text{'}}+b(\mathbf{x},t),$$

(8)

$$\left\{\begin{array}{l}{\mathbf{f}}_{\mathbf{x}\mathbf{x}\text{'}}=\mu c\left(\mathbf{x}\text{'},{\mathsf{\delta }}_{\mathbf{x}\text{'}}\right)s\frac{\mathsf{\eta }+\mathsf{\xi }}{\Vert \mathsf{\eta }+\mathsf{\xi }\Vert },\hfill \\ {\mathbf{f}}_{\mathbf{x}\text{'}\mathbf{x}}=-\mu c\left(\mathbf{x},{\mathsf{\delta }}_{\mathbf{x}}\right)s\frac{\mathsf{\eta }+\mathsf{\xi }}{\Vert \mathsf{\eta }+\mathsf{\xi }\Vert }.\hfill \end{array}\right.$$

(9)

The micromodulus, c, in Equation (9) depends on the particle spacing of the material point x and is half of the corresponding micromodulus in traditional PD, see Equation (3).

3. DH-PD Numerical Method

Figure 3 shows the program flowchart of the DH-PD model, which contains the pre-processing module, particle generation module, boundary treatment module, and damage degree calculation module. The C++ programming language is used to write the program code in this study.

3.1. Non-Uniform Discretization

The problem domain in the PD model is discretized into a series of particles with the material information in the reference configuration. There are no elements or other geometrical connections between the particles. The size of each particle is uniformly arranged and contains a certain volume. The spatial information of each particle is characterized by the geometric center. Horizons of the particles in the DH-PD model are variable, so, finer discretization with smaller horizons is always used in the concerning region, such as the damage area. Additionally, the variability of the horizons requires coarser particles for the unimportant areas.

The discretization is carried out layer by layer. First, the domain is discretized into a series of coarse particles with the same spacing. Then, the second-level discretization is carried out in the area with high stress concentration. The coordinates of the second-level particles can be determined based on their geometric relationship with the first-level particles. For example, a two-dimensional plane consists of four particles (P₀, P₁, P₂, P₃) with a volume of V and a spacing of Δx, as shown in Figure 4. For the second-level discretization, a coarse particle (P₁) with the coordinates (x₁, y₁) can be divided into four fine particles (P₁₀, P₁₁, P₁₂, P₁₃), with a volume of 0.25 V and a spacing of 0.5 Δx. Therefore, the coordinates of the fine particle P₁₁ are (x₁ + Δx/4, y₁ + Δx/4). In the same way, the second-level particle (P₁₁) can be divided into four third-level particles (P₁₁₀, P₁₁₁, P₁₁₂, P₁₁₃). The volume of the third-level particle is V/16, and the particle spacing is Δx/4. Based on this discretization, the coordinates of these finer particles can be easily determined by their geometrical relationships.

It is worth noting that, when the size difference between adjacent particles is too big, the finer particles cannot include any coarse particles in their horizons, which may result in bond deficiency and lead to unexpected failure near the interface. Therefore, in the presented model, the size difference between adjacent particles is suggested to be smaller than the radius of the finer particle’s horizon.

3.2. Volume Correction and Surface Correction

In PD theory, the problem domain is discretized into a series of particles, and each particle interacts with the neighboring particles in its horizon. However, some particles are located on the boundary and are not completely included in the horizon, as shown in Figure 5. Therefore, the model requires corrections with regards to the volumes of these particles when calculating the pairwise forces. For the multilevel discretization in DH-PD, the horizons of different particles vary with particle size, the particle volume can be corrected as:

$$V_j=\left\{\begin{array}{cc}\Delta {\mathbf{x}}^3\hfill & \hfill (\Vert \mathsf{\xi }\Vert \le \mathsf{\delta }-0.5\Delta \mathbf{x}),\hfill \\ \left[(\mathsf{\delta }+0.5\Delta \mathbf{x}-\Vert \mathsf{\xi }\Vert )/\Delta \mathbf{x}\right]\Delta {\mathbf{x}}^3\hfill & \hfill (\mathsf{\delta }-0.5\Delta \mathbf{x}<\Vert \mathsf{\xi }\Vert \le \mathsf{\delta }),\hfill \\ 0\hfill & \hfill \mathrm{otherwise}.\hfill \end{array}\right.$$

(10)

For the particles near the boundary of the problem domain, the number of neighboring particles in the horizon is lower than other particles, as shown in Figure 6. Therefore, the stiffness of the particles near the boundary is falsely reduced [41]. To deal with this problem, Scabbia et al. [42] presented an effective approach to mitigate the surface effect in state-based Peridynamics. Madenci and Oterkus [43] reconstructed the strain energy density near the boundary through the force density method and energy method to obtain the corresponding correction coefficients. The surface correction factor can be obtained by exerting a simple fictitious load to the object along the direction of the reference coordinate system. Based on this, the ratio of the strain energy densities calculated through the classical continuum mechanics theory and PD theory can be used as a valid surface correction factor.

3.3. Explicit Time Integration

Explicit time integration is the simplest solution for the PD equations in dynamics simulations. In this method, acceleration and velocity are approximatively determined by the finite difference form of the displacement [44]:

$${\ddot{\mathbf{u}}}_{(k)}^n=\frac{{\dot{\mathbf{u}}}_{(k)}^{n+1}-{\dot{\mathbf{u}}}_{(k)}^n}{\Delta t},$$

(11)

$${\dot{\mathbf{u}}}_{(k)}^{n+1}={\ddot{\mathbf{u}}}_{(k)}^n\Delta t+{\dot{\mathbf{u}}}_{(k)}^n.$$

(12)

In Equations (11) and (12), $\mathbf{u}$, $\dot{\mathbf{u}}$, $\ddot{\mathbf{u}}$ define the displacement, velocity, and acceleration vector, respectively. In addition, Δt represents the time increment, k denotes the particle number, and n equals the time step number. With this, displacement can be updated as:

$${\mathbf{u}}_{(k)}^{n+1}={\dot{\mathbf{u}}}_{(k)}^n\Delta t+{\mathbf{u}}_{(k)}^n.$$

(13)

Restricting the time step increment ensures the stability of the time integration, because the computational stability decreases as the time step increases. Silling and Askari [20] proposed a relationship,

$$\Delta t=\sqrt{\frac{2\rho }{{\displaystyle \sum _p\Delta }{V}_p{C}_p}},$$

(14)

to determine the critical time step in a PD model, where p is the number of the neighboring particles in the horizon, ΔV_p is the volume of the neighboring particle, and C_p is the microelastic modulus of the material.

4. Numerical Examples and Discussions

In this section, two benchmark tests are simulated to verify the accuracy and stability of the corrected DH-PD model: longitudinal vibration of a bar and wave reflection in a rectangular plate. The numerical results are compared with the analytical solutions and the finite difference method (FDM) results to validate the computational accuracy and the feasibility of the non-uniform discretization. Then, the corrected DH-PD model simulates two numerical examples: fracture of an L-shaped concrete specimen and mixed damage of a double-edged notched (DEN) specimen. The applicability and unique advantages of this model in crack propagation modeling are discussed.

4.1. Benchmark Problem 1: Longitudinal Vibration of a Bar

Figure 7 shows a bar with a fixed end on the left. The set parameters for the bar include the length of the bar, L = 1 m, the elastic modulus, E = 200 GPa, its Poisson’s ratio, ν = 0.33, and the density, ρ = 7850 kg/m³. The initial velocity of each material point is v_x (x, t) = 0, and the initial displacement gradient is defined as:

$$\frac{\partial {\mathbf{u}}_{\mathbf{x}}}{\partial x}=\epsilon \mathit{H}\left(\Delta t-t\right),$$

(15)

where H is a step function, and ε = 0.001 is the initial strain.

The bar is discretized into 1000 particles with a particle size of Δx = 0.001 m. Three fictitious boundary particles are fixed at the left side of the bar for the enforcement of displacement constraints. “No-fail zone” is commonly used for some applications under extreme loading conditions to avoid unexpected failure between the particles close to the external boundaries. For the loading conditions in the numerical examples presented in this paper, no-fail zone is not used. The computational time step is set as Δt = 0.2 μs, and the horizon associated with each particle is set as 3Δx. The particle with the initial position of x₀ = 0.5 m is selected as the analysis object. Its displacement time history is compared with the analytical solution [44], as shown in Figure 8. The expected longitudinal vibration is successfully captured by the DH-PD model after the volume correction and surface correction. Overall, the numerical results are consistent with the analytical solution, thus verifying the accuracy of the presented DH-PD model.

An error analysis is conducted, below, to assess the effect of the volume correction and surface correction on the calculation accuracy. The parameters are set $\Vert Error\Vert =‖\mathit{u}-{\mathit{u}}_{\mathrm{analytic}}‖$ (shown in Figure 9), and the analysis shows that the introduction of the volume correction and surface correction reduces the error to 0.01 mm.

Figure 10 shows the convergence curves of DH-PD and corrected DH-PD models at t = 4.0 ms. It indicates that the L₂ error in the displacement of the corrected DH-PD model is smaller than that of the DH-PD model, thus verifying that the volume correction and surface correction can improve the computational accuracy to some extent. Figure 10b shows that horizon size has little effect on the computational accuracy for the corrected DH-PH model in this case study. According to Silling and Askari [20], the horizon radius slightly greater than 3∆x usually works well. Therefore, the horizon radius is fixed at 3.015 times the particle size in this work.

4.2. Benchmark Problem 2: Wave Reflection in a Rectangular Plate

Figure 11 shows the geometric size of a rectangular plate in the initial state of: length, L = 10 m, width, W = 1 m, elastic modulus, E = 10 GPa, density, ρ = 2500 kg/m³, and Poisson’s ratio, ν = 0.33. The left end of the plate is fixed, and the right end is assigned a velocity boundary, v(t), based on Equation (16). In Equation (16), H (t₀ − t) is the Heaviside function, and t₀ = 2 ms is the duration of the velocity boundary. The plate is dispersed into 1600 particles, and there is a dense area in the middle section of the plate with the following PD parameters: matter point dispersion spacing (Δx = 0.1 m), particle spacing in the encrypted region (0.5Δx), and the time step (Δt = 1 μs).

$$v(t)=(1-\cos \frac{2\pi }{t_0}t)H(t_0-t).$$

(16)

The numerical results of the velocity wave propagation through the dense zone are shown in Figure 12 (DH-PD above, PD below). When the time is 2.4 ms, the velocity waveform transmitted is about to enter the dense area. Figure 12a shows that both PD and DH-PD simulate this stage accurately. When the time is 3.4 ms, the velocity of the center point in the dense area reaches the theoretical maximum of 2 m/s when using the DH-PD algorithm. When the traditional PD algorithm is used, the phenomenon of false wave reflection to the right will be found, weakening the wave peak of the velocity waveform in the dense area (Figure 12b). This violates the wave propagation theory in the same material and cannot completely simulate the wave propagation process. When the time reaches 4.6 ms, the velocity wave leaves the dense area and enters the non-dense area (Figure 12c).

The material points at the free end of the plate and the middle-end of the plate are monitored and their velocity characteristics are extracted and shown in Figure 13. The FDM simulation curves are almost identical to the analytical solutions, both of which can be used as a reference group. The conventional PD simulation is used to show anomalous spurious reflected waves. At a time of 3.5 ms, the velocity in the dense region is less than 1 m/s, which is half of the real value, and the spurious velocity wave is also received at the free end. DH-PD can simulate the complete propagation process of the wave, and it matches both the analytical solution and the FDM solution.

4.3. Numerical Application 1: Fracture of L-Shaped Concrete Specimen

Figure 14 shows the geometry of an L-shaped concrete specimen [45,46]. A total of 22,728 particles with different particle sizes and varying horizons are used in this case. Specifically, 7504 finer particles (Δx = Δy = 1.25 mm) are used near the corner of the L-shaped specimen, where the stress is relatively concentrated, 4300 coarser particles (Δx = Δy = 5 mm) are distributed at the periphery of the specimen, and 10,924 particles of size Δx = Δy = 2.5 mm are distributed between finer particles and coarser particles. The horizon size is set as 3.015 times the particle spacing. The information of the particles is shown in

Table 1. Values of peridynamics parameters for the DH-PD model.

Level	Δx = Δy	Horizon δ	Particle Numbers
I	5 mm	3.015 Δx	4300
II	2.5 mm	3.015 Δx	10,924
III	1.25 mm	3.015 Δx	7504
Total			22,728

and Figure 14. The other physical parameters include the modulus of elasticity, E = 25.85 GPa, the Poisson’s ratio, ν = 0.18, the energy release rate, G_c = 65 J/m², and the time increment, Δt = 0.1 μs. A displacement boundary, u, is exerted at the left side of the specimen, as shown in Figure 14.

Figure 15 shows the simulated crack propagation in the L-shaped plate. At t = 0.65 ms, a crack appears at the corner point of the L-shaped plate and propagates in the direction of 26° on the horizontal. Then, the direction of crack propagation changes to horizontal at a time close to t = 1 ms. The transverse crack expands to the right boundary of the specimen and cuts across the plate when t = 2.5 ms.

Figure 16 shows the numerical and experimental crack propagation paths. The simulated results essentially coincide with the experimental data recorded in the literature [45]. To quantitatively compare the results, the Pearson correlation coefficient of the experimental and numerical crack paths in Figure 16 is calculated. In this case, 52 points on the crack paths with an equal horizontal distance from the corner point are selected as the analysis objects. Assuming the vector X (x₁, x₂, x₃, ……, x₅₂) represents the calculated vertical coordinates of these points, and Y (y₁, y₂, y₃, ……, y₅₂) represents the measured ones in the experiment, the Pearson correlation coefficient of the two vectors can be calculated as:

$$Corr\left(\mathit{X},\mathit{Y}\right)=\frac{{\displaystyle \sum _{i=1}^n\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}}{\sqrt{{\displaystyle \sum _{i=1}^n{\left(x_i-\overline{x}\right)}^2\times {\displaystyle \sum _{i=1}^n{\left(y_i-\overline{y}\right)}^2}}}}.$$

(17)

In Equation (17), $\overline{\mathbf{x}}$ and $\overline{y}$ are the average values of the two data sets, respectively.

According to Equation (17), the correlation coefficient between the experimental and numerical results is calculated to be 0.988, suggesting that the numerical results have a strong correlation with the experimental data. Therefore, the calculation accuracy of the presented DH-PD model is verified.

To quantitatively compare the numerical and experimental results, the load-displacement curves are obtained based on the numerical results, and are compared with the experimental data in Figure 17. The comparison shows that the varying tendencies of the simulated curves match the experimental data well, and the calculated peak values of the load are close to the test value, thus verifying the validity and applicability of the presented model. In addition, the figure also shows that the corrected DH-PD model has a relatively higher simulation accuracy than DH-PD model.

4.4. Numerical Application 2: Mixed Damage of a DEN Specimen

Figure 18 shows a square thin plate with a length of 200 mm on each side, an elastic modulus of E = 30 GPa, a density of ρ = 2265 kg/m³, a Poisson’s ratio of ν = 0.33, and energy release rates of G_Ic = 110 J/m² and G_IIc = 1100 J/m². There is a prefabricated crack 25 mm in length and 5 mm in width on both sides of the specimen [46,47]. A uniform load, p, with a resultant force of 5 kN is applied in the horizontal direction on the left upper half and the right lower part of the plate. Meanwhile, a velocity load, v = 10 mm/s, is applied at the upper and lower boundaries of the plate, as shown in Figure 18.

The PD and DH-PD models are separately used to simulate this problem with the model parameters shown in

Table 2. Values of peridynamics parameters for the PD and DH-PD models.

Model	Δx = Δy	Horizon δ	Particle Numbers
PD	1.25 mm	3.015 Δx	25,440
DH-PD	2.5 mm, 1.25 mm	3.015 Δx	5040 + 5280 = 10,320

. In the PD model, the square specimen is discretized into 25,440 particles with a spacing of Δx = Δy = 1.25 mm. In the DH-PD model, the specimen is discretized into two kinds of particles: 5280 finer particles with a spacing of Δx = Δy = 1.25 mm in the middle part of the specimen, and 5040 coarser particles with a spacing of Δx = Δy = 2.5 mm near the boundary of the specimen. The uniform force, p, on both sides is transformed into a physical load, p/Δy, which is applied to the particles at the outermost layer of the plate. The velocity load, v, is exerted on the three layers of particles at the upper and lower boundaries of the specimen. Finally, the step length is set to Δt = 0.1 μs.

Figure 19 shows the simulated crack propagation in the DEN specimen based on both the traditional PD model and the DH-PD model proposed in this work. The crack paths and the propagation velocities obtained from these two models match each other. To compare the computational efficiency of the two models,

Table 3. Calculation time and efficiency of the two models.

Model	Particle Numbers	Calculation Time
PD	25,440	390.172 s
DH-PD	10,320	177.866 s
Efficiency Increased		119.4%

shows the calculation times of the same simulation using both models. The DH-PD model drastically improves the use of computational resources largely by increasing efficiency by 119.4% under the same operation mode.

Figure 20 compares the simulated crack paths with experimental data obtained by Nooru-Mohamed et al. [47], those calculated by the FEM model proposed in [48], and the traditional PD model. All the simulated crack paths coincide with the experimental results. To quantitatively analyze the computational accuracy of the corrected DH-PD model,

Table 4. Correlation coefficients between the test and numerical results on the crack paths.

Numerical Models	Corrected DH-PD	DH-PD	PD	FEM
Correlation Coefficients	0.957	0.952	0.952	0.929

lists the Pearson correlation coefficients of the numerical and experimental crack paths shown in Figure 20. The comparison demonstrates that the simulation results obtained by the corrected DH-PD model proposed in this study match the experimental results better than those obtained by other numerical models.

5. Conclusions

In this study, an efficient DH-PD model was developed based on the PD theory and the concept of dual-horizon. In this model, a technique of nonuniformity discretization with variable horizons was incorporated, which significantly improved the computational efficiency. Moreover, volume correction and surface correction were introduced to improve computational accuracy of the model. Finally, the DH-PD model was verified and applied to simulate crack propagation problems, producing valuable results.

The presented DH-PD model was verified through the simulation of benchmark problems. After comparing the numerical results against the analytical solution and FDM results, high goodness of fit was obtained to verify the accuracy of the corrected DH-PD model. No spurious wave reflection was found during the simulation, confirming the advantage of the corrected DH-PD model over the traditional PD method. Additionally, the surface correction and particle volume correction could substantially reduce the calculation error.

In the simulation of the fracture of an L-shaped concrete specimen, the specimen was discretized into particles with three different sizes and varying horizons. The simulated crack propagation produced a crack path consistent with the experimental results. Based on the Pearson correlation coefficient for quantitatively comparing the results, the coefficient value of 0.988 verified the calculation accuracy of the presented DH-PD model.

The corrected DH-PD model simulates the mixed damage of a DEN specimen. The simulated crack paths were compared against the experimental data and the results obtained by the traditional PD model and the FEM model. The comparison showed that the corrected DP-PD model has a higher simulation accuracy than the FEM model and is computationally more efficient than the traditional PD method.