An Anisotropic Peridynamic Model for Simulating Crack Propagation in Isotropic and Anisotropic Rocks

Abstract

In this work, we present a peridynamic-based simulation method for modeling quasi-static fracture propagation in isotropic and anisotropic rock within the framework of peridynamic least square minimization (PDLSM). The original isotropic elastic PDLSM is further extended to investigate fracture propagation in anisotropic materials in this study. The proposed AN-PDLSM model integrates an anisotropic model in fracture mechanics to analyze the failure process of anisotropic rocks. An important advancement in this research lies in the incorporation of the maximum energy release rate criterion (MERR) into the PDLSM framework for the first time. This enhancement enables accurately determining crack propagation and the associated crack angles. The proposed model utilizes the energy release rate calculated through the J-integral method to assess bond breakage, and it employs a mesh-independent, piecewise linear fracture model to describe crack propagation. The proposed method fully combines the merits of traditional fracture mechanics with the unique capabilities of peridynamics. To demonstrate the effectiveness of the proposed model, a simulation of fracture evolution in isotropic plates subjected to semi-circular bending tests is presented and compared with experimental results. It is shown that the proposed model accurately replicates fracture trajectories in isotropic specimens. In the context of anisotropic rock, the effect of a weak coefficient on crack morphology is discussed in order to obtain a suitable value. Additionally, the impact of bedding angles on fracture paths through our proposed model is also explored, revealing excellent agreement with experimental results. The findings in this research demonstrate that the proposed AN-PDLSM model is exceptionally proficient at capturing the intricate, oscillating crack paths observed in anisotropic rock materials.

1. Introduction

Nowadays, much attention has been paid to isotropic rocks, the strength properties of which are not direction dependent. However, sedimentary rocks, such as shale and coal, are very common in nature and used in different types of rock engineering. Sedimentary rocks always compose interior anisotropic structures formed by geological deposition, thus they should be regarded as anisotropic materials rather than isotropic materials. Crack propagation is the inner mechanism in the failure process of brittle rocks and determines the strength and deformation of rocks, thus it has a great influence on the stability of underground excavation, the efficiency of hydraulic fracture and other geotechnical engineering applications [1]. Consequently, in rock engineering, a thorough understanding of the internal mechanism of fracture propagation in anisotropic rocks is essential. However, this problem is more complicated than isotropic rock because the corresponding formulation must reflect the direction dependence of strength properties in anisotropic rock.

In order to analyze fracture processes in anisotropic materials, numerous numerical methods have been extended to simulate fracture propagation in anisotropic materials. The finite element method (FEM) is used to simulate crack propagation in anisotropic rocks [2]. A 2D anisotropic material model based on the Boundary Element Method (BEM) is proposed to simulate crack propagation in rocks [3]. The extended FEM (XFEM) is further developed to simulate crack propagation in 2D anisotropic media [4]. The lattice models are used to simulate the progressive development of rock failure, starting from micro-crack initiation to their coalescence into larger macro-cracks [5,6]. The phase-field method (PFM) is extended to simulate anisotropic material [7]. The combined finite and discrete element method (FDEM) is extended to model the fracture process in anisotropic rock [8]. Nevertheless, these approaches mentioned above are based on classical (local) continuum mechanics and require additional methods to deal with the spatial discontinuities, i.e., cracks.

Peridynamic (PD) theory proposed by Silling [9] is a kind of non-local mechanics that is found to be more suitable for modeling fracture processes than classical continuum mechanics (CCM) because its governing equations are integral equations rather than partial different equations. There are three types of PD theory, bond-based theory [10], ordinary state-based theory [10] and non-ordinary state-based theory [11]. Bond-based PD has a fixed Poisson’s ratio, whereas the other two types have no restriction. Non-ordinary state-based PD has the obvious advantage of inheriting PD’s ability to deal with cracks while enabling convenient implementations of existing failure models in CCM. The traditional non-ordinary state-based PD, however, suffers from numerical oscillation, which may have a significant impact on simulation results. To overcome this limitation, PD least squares minimization (PDLSM) [12] and the revised non-ordinary state-based PD model have been proposed [13]. PDLSM is an improved non-ordinary state-based PD based on the peridynamic differential operator proposed by Madenci (2016) [14] and has the advantage of efficiently providing commonly used information in CCM such as stress and strain. In the field of geotechnical engineering, peridynamic methods have been used to study the fracture processes in rocks, and good agreement with the experimental results has been reached [15,16].

Anisotropic peridynamic models are firstly proposed by Xu et al. [17] for anisotropic media, but there are some limitations in modeling. Ghajari et al. [18] proposed a bond-based model to remove these limitations, which can be used with any discretization of the domain but still has limitation on the ratio of two critical strain energy release rates. Zhang and Qiao [19] used ordinary state-based PD theory to extend this model and overcome the short-comings, but this model suffers from parameter value uncertainty. Hattori [20] simulated crack propagation in anisotropic solids using the anisotropic failure criterion in the framework of non-ordinary state-based PD. There are grid sensitivity problems because the model adopts the stress-based failure criterion. Tian and Zhou [21] extended the inspired PD model to analyze the fracture development in anisotropic rocks. Nevertheless, this method is affected by surface effects, and cannot effectively display stress, stress intensity factors (SIFs) and other information.

The anisotropy of rock strength is the main factor affecting crack propagation in anisotropic rock [22,23] and it is an effective analytical method to treat anisotropic rocks as rock with anisotropic fracture toughness [24]. Therefore, this paper also conducts research based on this model, and integrates an anisotropic model of the maximum energy release rate criterion (MERR) [25] into the PDLSM framework to analyze the fracture of anisotropic rocks. To the best knowledge of the authors, this is the first time that the anisotropic model MERR has been combined with peridynamics. Considering that the stress-based failure criterion may be strongly influenced by the mesh condition [26] in PD, a mesh-independent bond breakage method based on energy release rate G, which is obtained by the J-integral [27], is adopted in this paper to avoid this problem. The proposed AN-PDLSM model suited for the anisotropic rock can be easily reduced to an isotropic model and provides an efficient way to analyze crack propagation in both isotropic and anisotropic rocks. Crack propagation in isotropic and anisotropic plates under semi-circular bending tests is simulated based on the proposed model and shows good agreement with the experimental results.

The remainder of this paper is organized as follows. First, the basic theory of PDLSM is introduced in Section 2. The proposed model AN-PDLSM for anisotropic rocks and a description of the numerical implementation is presented in Section 3. Two numerical examples that validate the proposed model and its ability to simulate fracture propagation in both isotropic and isotropic rocks are included in Section 4. Finally, concluding remarks are summarized in Section 5.

2. A Brief Introduction to Peridynamic Least Square Minimization (PDLSM)

Peridynamic least square minimization (PDLSM) was first developed by Madenci et al. [12], and combines the concepts of peridynamic differential operator [14] and least square minimization (LSM). In contrast to standard PD theory, PDLSM encompasses the concepts of stress and strain, enabling fully utilizing the constitutive models in traditional continuum mechanics (CCM), as well as the stress intensity factors (SIFs) and the J-integral in LEFM to analyze the fracture process. PDLSM is applicable to any arbitrary interaction domain, whereas standard PD theory is limited to sphere interaction domains. As a result, the PDLSM is free of surface effects and volume correction, and it could provide more accurate stress and strain tensor for a wide range of issues compared with standard PD theory.

In PDLSM, the equation of motion for linear isotropic material is given by

(1)$\begin{array}{cc}\rho \left(\mathit{x}\right)\ddot{\mathit{u}}\left(\mathit{x},t\right)& ={\mathit{L}}^{\mathrm{p}\mathrm{d}}+\mathit{b}\left(\mathit{x},\mathit{t}\right)\\ & ={\int }_{{\mathrm{H}}_x}^{}\omega \left(\left|\mathit{\xi }\right|\right)\mathit{G}\cdot \mathit{\eta }\mathrm{d}{V}_{{\mathit{x}}^{\prime }}+\mathit{b}(\mathit{x},t)\end{array}$

where ρ is the mass density, $\ddot{\mathit{u}}$ is the acceleration of node x, b is the external force vector, and ξ = x′ − x and η = u(x′) − u(x) are the relative position and relative displacement of nodes x′ and x, respectively, ω(|ξ|) is the weight function, and the integral part is the internal force vector ${\mathit{L}}^{\mathrm{p}\mathrm{d}}$, which is calculated by the integration of bond force between node x and x′ over the interaction domain H_x of node x, as shown in Figure 1. The matrix G connects linear elastic Hooke’s law and peridynamic internal force vector ${\mathit{L}}^{\mathrm{p}\mathrm{d}}$, according to Liu et al. [26] the matrix G in the 2D case can be defined for plane stress problems as

(2)$\left[\mathit{G}\right]=\left[\begin{array}{cc}\frac{E}{2\left(1-\nu \right)}{d}_1+S\left(d_1+d_2\right)& \frac{E}{2\left(1-\nu \right)}{d}_3\\ \frac{E}{2\left(1-\nu \right)}{d}_3& \frac{E}{2\left(1-\nu \right)}{d}_2+S\left(d_1+d_2\right)\end{array}\right]$

and for plane strain problems

(3)$\left[\mathit{G}\right]=\left[\begin{array}{cc}\left(\lambda +S\right)d_1+S\left(d_1+d_2\right)& \left(\lambda +S\right)d_3\\ \left(\lambda +S\right)d_3& \left(\lambda +S\right)d_2+S\left(d_1+d_2\right)\end{array}\right]$

where λ is the Lame’s material constant, S is the shear modulus, ν is Poisson’s ratio, E is Young’s modulus, and d₁, d₂, and d₃ are peridynamic differential operators, defined as

(4)$[\frac{\partial }{\partial x} \frac{\partial }{\partial y} \frac{{\partial }^2}{\partial x^2} \frac{{\partial }^2}{\partial y^2} \frac{{\partial }^2}{\partial xy}]\mathit{u}(\mathit{x},t)={\int }_{{\mathit{H}}_{\mathit{x}}}^{}\omega \left(\left|\mathit{\xi }\right|\right)[d_1 d_2 { d}_3 g_1 g_2]\cdot \mathit{\eta }\mathrm{d}{V}_{x^{\prime }}$

The non-local displacement gradient can be expressed according to the peridynamic differential operator theory,

(5)$𝛁⨂\mathit{u}(\mathit{x},t)={\int }_{H_x}^{}\omega \left(\left|\mathit{\xi }\right|\right)\mathit{\eta }⨂\mathit{g}d{V}_{\mathit{x}\prime }$

where ⨂ is the inner product and g = [g₁ g₂] is the peridynamic differential operator vector, as explained by Liu et al. [27] and Madenci et al. [12].

According to the CCM, the stain tensor ε in the framework of infinitesimal deformation is expressed as

(6)$\mathit{\epsilon }=\frac{1}{2}(𝛁⨂\mathit{u}+\mathit{u}⨂𝛁)$

Thus, the stress tensor σ could also be determined directly based on the strain tensor ε and the stiffness tensor D, given as

(7)$\mathit{\sigma }=\mathit{D}:\mathit{\epsilon }$

3. AN-PDLSM Model for Analyzing Fracture in Anisotropic Rock

In this section, we propose AN-PDLSM, a new method for simulating quasi-static fracture propagation in anisotropic materials, particularly brittle rock, based on the framework of PDLSM and linear elastic fracture mechanics (LEFM). The displacement, strain tensor, and stress tensor are calculated by PDLSM. In this paper, the maximum energy release rate (MERR) criterion proposed by Gao et al. [23] and Wei et al. [22] is adopted to analyze the fracture behavior for anisotropic rocks. In the framework of PDLSM, the MERR criterion is evaluated from stress intensity factors (SIFs), which are calculated from a path-independent J-integral. Thus, the proposed method AN-PDLSM makes full use of the existing theories in LEFM to analyze the fracture behavior of anisotropic rocks within the framework of PDLSM, thereby broadening PDLSM’s scope of application.

3.1. The Crack Propagation Criterion for Anisotropic Rock

An effective method for studying the fracture behavior of anisotropic rocks is to consider anisotropic rock as rock with anisotropic fracture toughness, which has been widely used in the field of rock mechanics [22,23]. This paper adopts this useful method and combines it with peridynamics in the framework of PDLSM for the first time. In this method, the fracture toughness G_c(θ) is assumed to vary with the orientation angle θ in anisotropic materials, while remaining constant in isotropic materials. According to the study of Hakim et al. [28], the MERR criterion is extended to materials with anisotropic fracture energy as

(8)$\underset{\theta }{\max }\frac{G\left(\theta \right)}{G_c\left(\theta \right)}>1$

where G(θ) is the energy release rate in different directions.

According to Amestoy and Leblond (1992), the energy release rate G(β) in a specific direction β is defined as

(9)$G\left(\beta \right)=\alpha \left[K_I^2\left(\beta \right)+K_{II}^2\left(\beta \right)\right]$

Here, K_I(β) and K_II(β) are kinked stress intensity factors (SIFs) in mode I and mode II with a deflection angle β of the old crack propagation direction, respectively. Parameter α is denoted by

(10)$\alpha =\left\{\begin{array}{c}\frac{(1-{\upsilon }^2)}{E} (plane strain)\\ \frac{1}{E} (plane stress)\end{array}\right.$

The relationship between the kinked SIFs {K_I(β), K_II(β)} and the SIFs before deflection {K_I, K_II} has been studied by Niusmer [29], and Amestoy and Leblond [30]. In this study, we decide to use the simplest but useful model proposed by Niusmer [29], which is given as

(11)$$\begin{gathered} K_I(\beta )=\frac{1}{2}\mathrm{c}\mathrm{o}\mathrm{s}\frac{\beta }{2}[K_I\left(1+\mathrm{c}\mathrm{o}\mathrm{s}\beta \right)-3K_{II}\mathrm{s}\mathrm{i}\mathrm{n}\beta ] \\ {K}_{II}(\beta )=\frac{1}{2}\mathrm{c}\mathrm{o}\mathrm{s}\frac{\beta }{2}[K_I\mathrm{s}\mathrm{i}\mathrm{n}\beta +K_{II}\left(3\mathrm{c}\mathrm{o}\mathrm{s}\beta -1\right)] \end{gathered}$$

Because the presence of weak planes in anisotropic rocks causes anisotropy in fracture toughness, Gao et al. [23], Zeng and Wei [22] assumed that the rocks have an intact fracture toughness G_ci in all directions except the direction of weak plane, and a weak planar fracture toughness G_cw in the direction of the weak plane, as shown in Figure 2. The weak plane model is utilized in this paper to describe the anisotropic rock, which is given by

(12)$G_c(\theta )=\left\{\begin{array}{c}{G}_{cw} (\theta =\gamma or \gamma +\pi )\\ G_{ci} (\mathrm{other} \mathrm{directions})\end{array}\right.$

where G_c(θ) is the fracture energy in the direction of θ and γ is the direction of weak plane. It should be noted that this weak plane model can be reduced to the maximum energy release rate criterion for isotropic material by setting G_cw equals to G_ci.

Determining the crack propagation angle θ_c is an important part in analyzing fracture propagation, which is defined as the angle between the direction of old crack and the direction of new crack. The value of θ_c in the isotropic case is denoted by θ_ci, which is calculated by considering the maximum value of the energy release rate G(θ) in Equation (9), and is found to be

(13)${\theta }_c={\theta }_{ci}=\left\{\begin{array}{c}2{\mathrm{t}\mathrm{a}\mathrm{n}}^{-1}\left(\frac{K_I}{4K_{II}}-\frac{1}{4}\sqrt{{\left(\frac{K_I}{K_{II}}\right)}^2+8} \right) (if K_{II}>0)\\ 2{\mathrm{t}\mathrm{a}\mathrm{n}}^{-1}\left(\frac{K_I}{4K_{II}}+\frac{1}{4}\sqrt{{\left(\frac{K_I}{K_{II}}\right)}^2+8} \right) (if K_{II}<0)\end{array}\right.$

However, the value of θ_c in the anisotropic case is obtained in a more complicated way than in the isotropic case because of the existence of the weak planes. And the propagation angle θ_c is calculated according to the MERR criterion for material with anisotropic fracture energy, which is given by

(14)${\theta }_c=\left\{\begin{array}{c}{ \theta }_{ci} ( if \frac{G\left({\theta }_{ci}\right)}{G_{ci}}\ge \frac{G\left(\beta \right)}{G_{cw}} )\\ \beta ( if \frac{G\left({\theta }_{ci}\right)}{G_{ci}}<\frac{G\left(\beta \right)}{G_{cw}} )\end{array}\right.$

Once the crack propagation angle θ_c is determined using Equation (14), the condition for determining the crack propagation could be established. If G(θ_c) > G_c(θ_c), crack propagate will occur in the direction of θ_c.

According to the author’s knowledge, there are many other anisotropic criteria for the occurrence of cracks for anisotropy solid, such as the anisotropic maximum circumferential stress criterion [31] and anisotropic maximum strain theory [32]. The MERR model used in this article is an anisotropic criterion based on fracture energy. When a crack reaches the critical fracture energy, the crack will propagate. Therefore, this MERR model can effectively consider the occurrence of cracks. Of course, other types of fracture criteria could be considered in future work.

3.2. Fracture Energy and the J-Integral

In order to accomplish the MERR within the framework of PDLSM, stress intensity factors (SIFs) and energy release rate G must be precisely calculated. Considering that the J-integral model is the most common and robust method to obtain the SIFs and is widely used in XFEM [23], the J-integral is used to calculate the SIFs in this paper. And the J-integral is defined as

(15)$J={\int }_{\Gamma }^{}\left(\frac{1}{2}{\sigma }_{ij}{\epsilon }_{ij}{\delta }_{1k}-{\sigma }_{ik}\frac{\partial u_i}{\partial x_1}\right)n_{k}d\Gamma$

in which δ is the Kronecker delta, σ_ij is the stress tensor, ε_ij is the strain tensor, u_k is the displacement vector, and n_k is the unit outward normal vector on the integration contour Γ.

In the field of linear elastic fracture mechanics (LEFM), the J-integral is equal to the energy release rate G, and is related to the SIFs K_I and K_II by the well-known relationships,

(16)$J=G=\alpha {(K}_I^2+K_{II}^2)$

where the SIFs K_I and K_II can be directly determined by solving two equations based on the interaction integral (I-integral) theory proposed by Yau et al. [33],

(17)$\left\{\begin{array}{c}{M}^{\left(\mathrm{1,2}a\right)}={2\alpha K}_I\\ M^{\left(\mathrm{1,2}b\right)}={2\alpha K}_{II}\end{array}\right.$

where M^(1,2a) and M^(1,2b) are two interaction integrals defined as

(18)$\left\{\begin{array}{c}{M}^{\left(\mathrm{1,2}a\right)}={\int }_{\mathsf{\Gamma }}^{}\left({\sigma }_{ij}^{\left(1\right)}{\epsilon }_{ij}^{\left(2a\right)}{\delta }_{1j}-{\sigma }_{ij}^{\left(1\right)}\frac{\partial u_i^{\left(2a\right)}}{\partial x_1}-{\sigma }_{ij}^{\left(2a\right)}\frac{\partial u_i^{\left(1\right)}}{\partial x_1}\right)n_{j}d\Gamma \\ M^{\left(\mathrm{1,2}b\right)}={\int }_{\mathsf{\Gamma }}^{}\left({\sigma }_{ij}^{\left(1\right)}{\epsilon }_{ij}^{\left(2b\right)}{\delta }_{1j}-{\sigma }_{ij}^{\left(1\right)}\frac{\partial u_i^{\left(2b\right)}}{\partial x_1}-{\sigma }_{ij}^{\left(2b\right)}\frac{\partial u_i^{\left(1\right)}}{\partial x_1}\right)n_{j}d\Gamma \end{array}\right.$

{${\sigma }_{ij}^{\left(1\right)}$, $u_i^{\left(1\right)}$} denote the actual solution of stress tensor and displacement vector calculated by PDLSM; {${\epsilon }_{ij}^{\left(2a\right)}$, $u_i^{\left(2a\right)}$, ${\sigma }_{ij}^{\left(2a\right)}$} are solutions of strain tensor, displacement vector, stress tensor under the assumption of K_I = 1 and K_II = 0; and {${\epsilon }_{ij}^{\left(2b\right)}$, $u_i^{\left(2b\right)}$, ${\sigma }_{ij}^{\left(2b\right)}$} are solutions of strain tensor, displacement vector, stress tensor under the assumption of K_I = 0 and K_II = 1. The detailed description of those parts could be found in Yau’s [33] paper.

From a numerical perspective, an integration contour path should be specified to obtain the J-integral. The numerical integration method proposed by Liu et al. [26] is employed in this study, as shown in Figure 3. The selected PD nodes and integral path for calculating the J-integral could be de determined if the integral radius r is specified, which should be a suitable value, usually 6 fold the minimal node space. With the integral path from the selected PD nodes, the J-integral could be calculated using Equation (16).

3.3. A Mesh-Independent Piecewise Linear Model for Bond Breakage

The traditional XFEM, CZM and FDEM methods usually limit fracture propagation along the element boundaries and could only deflect at the element boundaries. Obviously, those methods are sensitive to the mesh condition and require fine mesh to reduce this effect, while fine mesh will significantly increase the computational burden.

Peridynamic bonds are assumed to be broken when they reach a certain state, and the fracture process is accomplished by breaking the bonds that connects two PD nodes. Existing studies have shown that this bond breakage method also exhibits grid dependencies [34]. In address this shortcoming, a mesh-independent piecewise linear model for bond breakage is adopted in this paper to simulate the fracture process of rocks, as shown in Figure 4. Once the new crack is determined using the MERR criterion described in Section 3.1, the bonds intersected with the new crack are set to break and the interaction force on the bond vanishes. This bond breaking model differs from the original bond breaking manner. As mentioned by Bazant et al. [35], the original bond breakage method may cause an unreal phenomenon since the PD nodes, on opposite sides of the crack, may have ghost interaction in some cases. On the contrary, this mesh-independent piecewise linear model could eliminate all interactions, making it more physical, and is adopted by some researchers like Liu et al. [26] and Ni et al. [36].

In the proposed model, the new crack could propagate in any direction, not limited to the element boundaries. In numerical simulation, the fracture propagation direction can be determined by the LEFM exactly, whereas the fracture growth length is chosen in a somewhat arbitrary way. The fracture growth length is set to be proportional to the sqrt of minimum area of PD nodes according to Liu et al. [26].

(19)${∆}_{min}=k\sqrt{A_{min}}$

where A_min is the minimum area of PD nodes and k is a fracture growth length parameter, which is set to be 1.0 in this study.

A scalar variable μ(ξ,t) is used in the AN-PDLSM model to record the damage state of the bonds varying with time. Once the bond reaches a critical state, the bond is set to be broken, and the force on the bond vanishes. The scalar variable μ(ξ,t) is defined as

(20)$\mu (\xi ,t)=\left\{\begin{array}{c}1 (\mathrm{if} \mathrm{the} \mathrm{bond} \mathrm{is} \mathrm{intact})\\ 0 (\mathrm{if} \mathrm{the} \mathrm{bond} \mathrm{is} \mathrm{broken})\end{array}\right.$

A damage variable D is fined on the node level rather than the bond level in order to conveniently display the damage of the entire body. The damage D(x,t) of a node x is calculated by the proportion of broken bonds and intact bonds in the corresponding interaction domain H_x, which is defined as

(21)$D\left(\mathit{x},t\right)=1-\frac{{\int }_{H_x}^{}\mu (\xi ,t)\mathrm{d}{V}_{{\mathit{x}}^{\prime }}}{{\int }_{H_x}^{}\mathrm{d}{V}_{{\mathit{x}}^{\prime }}}$

3.4. Implicit Solver of AN-PDLSM

The quasi-static fracture problem is widely studied in the fields of rock mechanics and rock engineering, and this paper aims to model a quasi-static fracture process rather than a dynamic fracture process. The explicit local damping method [37] and the adaptive dynamic relaxation method (ADR) [38] are widely used and preferred for quasi-static problems. However, compared to those explicit methods, the implicit algorithm is better suited for quasi-static problems in some cases because it does not require satisfying the C-F-L condition and is unconditionally convergent [36].

Considering the foregoing, this paper employs the implicit method to solve the governing equation, i.e., the equilibrium equations. The weak form of AN-PDLSM is obtained by applying the weighted residual method (WRM) on the equation of motion and the boundary conditions (B.C.), which is given as

(22)${\int }_{V_{pd}}^{}\delta \mathit{u}\cdot ({\mathit{L}}_{pd}+\mathit{b}-\rho \ddot{\mathit{u}})\mathrm{d}V-{\int }_{S_{pd}}^{}\delta \mathit{u}\cdot ({\mathit{\sigma }}_{pd}\cdot \mathit{n}-\mathit{T})\mathrm{d}S=0$

where T is the external traction on the boundary, n is the unit normal vector of the boundary, and δu is a virtual displacement. V_pd means the volume of the entire body, and S_pd is the boundary of the body.

According to the study of Liu et al. (2021b), the matrix form of the governing Equation (22) can be written in the same form as FEM, which is given as

(23)${\mathbf{M}}^{pd}\ddot{\mathit{u}} +{\mathit{K}}^{pd}\mathit{u}={\mathit{F}}^{pd}$

and the global stiffness matrix K^pd is defined as

(24)${\mathit{K}}^{pd}={\mathit{K}}_{\mathit{b}}^{pd}+{\mathit{K}}_{\mathit{s}}^{pd}$

where u is the nodal displacement vector, F^PD is the global equivalent node force, and M^PD is the global equivalent nodal mass. In the AN-PDLSM model, the global stiffness matrix K^PD consists of two parts, K_b^PD associated with PD domain and K_s^PD associated with PD boundary. The stiffness matrix K^PD varies with the bond breakage state. The matrix will be updated if any bond is set to be broken, and the corresponding stiffness of the broken bond vanishes. The detailed description of the global stiffness matrix K^PD could refer to the study of Liu [26] and Medanci [12].

The global displacement vector u can be expressed as

(25)$\mathit{u}=\left[\begin{array}{c}{\mathit{u}}_u\\ {\mathit{u}}_p\end{array}\right]$

Here, u_u is the set of all unknown displacement components and u_p is the set of prescribed displacement components. For static or quasi-static problems, the acceleration vector $\ddot{\mathit{u}}$ is negligibly small, the dynamic term ${\mathbf{M}}^{pd}\ddot{\mathit{u}} \mathrm{i}\mathrm{s} \mathrm{v}\mathrm{a}\mathrm{n}\mathrm{i}\mathrm{s}\mathrm{h}\mathrm{e}\mathrm{d}$. Thus, governing Equation (23) can be reduced to

(26)$\left[\begin{array}{cc}{\mathbf{K}}_{uu}& {\mathbf{K}}_{up}\\ {\mathbf{K}}_{pu}& {\mathbf{K}}_{pp}\end{array}\right]\left[\begin{array}{c}{\mathit{u}}_u\\ {\mathit{u}}_p\end{array}\right]=\left[\begin{array}{c}{\mathit{F}}_u\\ {\mathit{F}}_p\end{array}\right]$

Since u_p is the prescribed displacement constraint, the governing Equation (26) for quasi-static problems reduces to

(27)${\mathbf{K}}_{uu}{\mathit{u}}_u={\mathit{F}}_u-{\mathbf{K}}_{up}{\mathit{u}}_p$

By solving the above large linear systems of equations, the unknown displacement components u_u can be solved.

In this paper, these large linear systems of equations are solved by the Parallel Direct Sparse Solver of Intel Math Kernel Library (MKL). Compared to the explicit dynamic relaxation method, directly and implicitly solving Equation (27) avoids a large number of iterations and is highly computationally efficient. The direct solving method based on LU decomposition is selected to solve the large linear systems of equations. The large sparse matrix is stored in a compressed sparse row (CSR) format, which is one of the inputs of the Parallel Direct Sparse Solver. The above numerical calculation process is implemented using code written in C++.

Generally, the explicit dynamic relaxation method requires a large number of iterations for each time step to reach a convergent state. Compared to the explicit ADR method, the implicit method does not need a large number of iterations for solutions, and convergent is guaranteed in each times step. Although the implicit method requires to solve the inversion of matrix K^PD, the stiffness matrix K^PD is symmetric and sparse and its inversion is relatively easier and fast computationally using the Parallel Direct Sparse Solver.

In the traditional PD theory, applying boundary conditions is much more complex than in FEM. In those PD theories, a fictitious layer with non-zero volume is required for applying boundary conditions. However, within the framework of PDLSM, the constraint and traction boundary conditions can be treated the same way as the FEM method, making it more convenient than traditional PD theories. And those feathers are the advantages of the AN-PDLSM method compared with classic PD theory.

The implicit numerical implementation is divided into seven steps, step 1 is to build the computing mesh model; Step 2 is to input calculation parameters; Step 3 is to apply boundary conditions; Step 4 is to integrate the stiffness matrix K^pd and node force vector F^pd; Step 5 is to solve linear algebraic equations for displacement, strain and stress; Step 6 is to calculate SIFs and the energy release rate G according to the J-integral; Step 7 is to calculate the crack propagation angle θ_c and the corresponding energy release rate G(θ_c); Step 8 is to determine the bond breakage condition and break the bond according to the MERR.

In Step 1 of establishing a computational mesh, the traditional PD often uses uniform mesh, but will greatly increase the amount of computation because it divides the mesh of unimportant areas as small as the areas of interest. On the contrary, the proposed AN-PDLSM can use non-uniform mesh, with fine mesh in the regions of interest and coarse mesh in the less important regions. Previous studies have shown that PDLSM has grid-independent properties, and the proposed AN-PDLSM also inherits this advantage, which is superior to the traditional PD method [12]. All the computational meshes in this paper are generated by commercial software Abaqus 2020, but only nodal information is used, while information of elements is not required, because AN-PDLSM is still a meshless method in essence.

4. Numerical Verification and Simulation

Two plane stress examples are performed to demonstrate the capability of the proposed AN-PDLSM model to simulate fracture propagation in isotropic and anisotropic rocks. The first example verifies the fracture propagation of an isotropic plate under a semi-circular bending test, while the second example simulates the fracture process of a anisotropic semi-circular bending specimen and discusses the influence of a weak coefficient and a weak planar direction on the crack trajectory.

4.1. Fracture Analysis of Isotropic Semi-Circular Bending Tests

As a benchmark problem, an isotropic semi-circular bending specimen made of polymethylmethacrylate (PMMA) with a prefabricated notch is simulated to validate the capability of the AN-PDLSM model to capture crack evolution in isotropic materials, as shown in Figure 5. If the material is isotropic, the AN-PDLSM model is reduced to the PDLSM model, as described in Section 3. The isotropic semi-circular bending plate is subjected to displacement loading u, which increases from zero with incremental Δu =1.6 × 10⁻⁵ m on the top of the numerical model. The displacement of the left support pin in the x and y directions is fixed, while the y direction of the right pin is fixed and the x direction is free. The angle formed by the prefabricated notch and the vertical line is denoted by α (Figure 5a). The dimensions of the specimen are determined as radius R = 50 mm, crack length a = 15 mm, support distance 2S = 43 mm, and thickness t = 5 mm. The mechanical properties of PMMA are specified as Young’s modulus E = 3.75 GPa, Poisson’s ratio ν = 0.31, and fracture toughness K_IC = 2.13 MPa·m^1/2.

As shown in Figure 5b, the entire domain is discretized with 10,512 PD nodes with a non-uniform mesh rather than uniform mesh, with the mesh of the middle part of the specimen being finer than the edge part. The size of the horizon is specified as δ = 3d, where d is the PD node space. The corresponding numerical simulation process is carried out as described in Section 3.4. Since the sample is isotropic, the intact fracture toughness G_ci is equal to the weak planar fracture toughness G_cw, and the crack propagation angle θ_c is calculated using Equation (13).

Figure 6 illustrate the numerical results of five specimens with different angles α = 0°, 10°, 20°, 30°, 40°. When α = 0°, the fracture propagates along the direction of the prefabricated notch, resulting in pure model I fracture pattern. When α is set to other angles, the fracture propagate propagates as a wing crack from the notch tip and eventually aligns with the horizontal direction. The fracture trajectories of the five cases clearly agree with the experiment results of Ayatollahi et al. [39]. Consequently, the proposed AN-PDLSM model can be reduce to the isotropic case and effectively simulate the fracture process of an isotropic rock under a semi-circular bending test.

In order to analyze the robustness of the proposed model and examine the influence of mesh size on the calculation results. Semi-circular bending tests with different mesh sizes are simulated, respectively, and the notch angle α of the specimens is set to 30°. Because non-uniform meshes are used in this paper, the number of nodes is used to represent the size of the mesh. Four cases with different mesh size are set, as shown in Table 1. To ensure a single variable, the time steps of the four numerical cases are set to 350 steps. Table 1 shows the calculation time (CPU time) for different cases, and Figure 7 plots the relationship between the calculation time and the number of nodes. The results show that the calculation time has an approximately linear relationship with the number of nodes.

Figure 8 shows the crack shapes corresponding to different mesh size. It can be found that the crack shapes of the four cases are generally consistent. The results suggest that the mesh size does not have a significant impact on the numerical results. Even coarser meshes simulate crack shapes well and achieve good results. According to the author’s experience, to obtain the best simulation results, the grid distance dx should be smaller than 1/100 of the length of the sample.

4.2. Fracture Analysis of Anisotropic Semi-Circular Bending Tests

4.2.1. Influence of a Bedding Weak Coefficient

In this section, the proposed AN-PDLSM model described in Section 3 is used to analyze the crack propagation of anisotropic rock subjected to a semi-circular bending test. An anisotropic semi-circular plate is subjected to displacement loading u increasing from zero with an incremental Δu = 1.6 × 10⁻⁵ m on the top of the specimen. The dimensions and the boundary conditions of this specimen are determined in the same way as the isotropic model described in Section 4.1. The mechanical properties of the specimen are specified according to the study of Wang and Teng. [40] as E = 2326.31 MPa, ν = 0.31, and intact fracture toughness G_ci = 20 J. The entire domain is discretized with 10,512 PD nodes, the same as the numerical model in Section 4.1. The size of the horizon is specified as δ = 3d, where d is the PD node space.

In the anisotropic model, the fracture toughness G_c is no longer constant, but contains two different values G_cw and G_ci. For the convenience of analysis, the ratio of G_cw and G_ci is defined as the bedding weak coefficient η = G_cw/G_ci, which reflects the degree of softening of the bedding relative to the intact rock. The bedding weak coefficient η = 1 indicates that the bedding is consistent with the intact rock, and η close to 0 indicates that the bedding is very weak. In order to quantitatively evaluate the influence of a weak coefficient η on the fracture pattern of rock under a semi-circular bending test, several simulations were carried out in the case of η = 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, and other variables remained unchanged. The Angle of prefabricated notch was 0°, crack length a is equal to 15 mm and the direction of the bedding is π/6.

Figure 9 illustrate the fracture trajectory of specimens with a different weak coefficient η. The numerical results show that η has a significantly effect on the morphology of fracture propagation in a semi-circular bending test. When η = 0.5, the bedding strength is much lower than the intact rock, and the fracture extends completely along the bedding direction. When η = 1.0, the bedding strength is equal to that of intact rock, i.e., there is no bedding, and the fracture propagates vertically upward completely, which is consistent with the isotropic case in Section 4.1. When η is between 0.5 and 1.0, the fracture presents a typical two-stage expansion, with the crack first expanding along the bedding, which is referred to as bedding crack in this paper, and then deviating from the bedding direction to the vertical line, which is referred to as wing crack in this paper. Bedding crack represents crack along the bedding, while wing crack represents fracture that occurs in intact rock.

Figure 10 depicts a comparison of the crack morphology and experiment results for various weak coefficients η. The turning point is the point at which the crack transitions from bedding crack to wing crack. As η decreases from 1.0 to 0.5, the turning point appears closer to the top part of the specimen; The proportion of bedding cracks increases while the proportion of wing crack decreases. Compared with the experiment results by Wang and Teng. [40], it is clear that when η equals to 0.6, the numerical fracture path is basically consistent with the experimental results. Thus, η = 0.6 is selected to describe anisotropic rock in the proposed AN-PDLSM model in the following part.

Unlike the bond-based peridynamic model, which lacks the definition of stress tensor, AN-PDLSM contains the concept of stress, enabling it to analyze the evolution of stress during the failure process. Because horizontal tensile failure will occur in the rock sample a under semi-circular bending test, we provide the contour of horizontal stress σ_x corresponding to the failure process of plate under loading when η = 0.6, as shown in Figure 11. The tensile stress σ_x at the crack tip presents a symmetrical distribution when the fracture begins to propagate, as shown in Figure 11a; When the fracture reaches the turning point, the stress distribution at the crack tip changes, the direction of the stress deflates, and the crack type changes from bedding crack to wing crack. The stress symmetry axis of the wing crack tip is generally perpendicular to the crack.

In order to analyze the mechanism of the turning point, the SIFs {K_I, K_II} and energy release rate G of each time step have been recorded, as shown in Figure 12. The loading stage can be divided into three different stages based on the variation in energy release rate G during the loading process. The first stage is the loading stage without crack propagation. The numerical results show that at this stage, the K_I increases linearly under loading, whereas the K_II remains zero. Thus, the fracture pattern in stage one is type I crack mode in this test. The energy release rate G increases in the form of a quadratic function, which is consistent with Formula (16). The second stage is the bedding crack propagation stage. Figure 12b shows that the energy release rate G in this stage is fixed around G_cw, which is equal to η*G_ci, indicating the accuracy of the calculated results. In the second stage, the SIFs K_I and K_II decreases almost linearly, respectively. The third stage corresponds to the wing crack propagation stage. In this stage, most of the energy release rate G equals to G_ci, while only a small part equals to G_cw. It implies that the predominant crack at this stage is the wing crack, with the bedding crack making up a much smaller portion.

When compared to the isotropic case in Section 4.1, the crack trajectory of anisotropic model is clearly more complicated than that of the isotropic model due to the presence of the oscillation cracks and transition zones, which have been reported in anisotropic rocks [41]. According to the experimental results of Wang and Teng [40], the fracture starts from the notch tip and develops along the bedding angle before deflecting in the opposite direction in the transition zone, producing the intricate oscillation cracks depicted in Figure 10b. The proposed AN-PDLSM model captures this complex crack path well, demonstrating the method’s effectiveness in simulating crack propagation in anisotropic rock under a semi-circular bending test.

4.2.2. Influence of Bedding Plane Angle

On the basis of Section 4.2.1, the proposed AN-PDLSM model is used to analyze the influence of the bedding angle on the fracture trajectory. The specimen with 5 different bedding angles 0°, 22.5°, 45°, 67.5°, and 90° are simulated, respectively. Except for the bedding angle, the other parameters are consistent with Section 4.2.1. It should be noted that in those simulations, the weak coefficient η is set to 0.6.

All of those numerical results are compared with the experimental results, respectively, as shown in Figure 13. The numerical results are found to be in good agreement with the experimental results. The crack propagates vertically upward when the bedding angle θ is equal to 0° and 90°, respectively. When the bedding angle θ is equal to 22.5°, 45° and 67.5°, the crack path appears as a complex oscillating crack; In the experimental results, the oscillating crack first propagates along the bedding before converting to a wing crack after reaching the turning point. Furthermore, as the bedding angle θ increases from 22.5° to 67.5°, the length of the bedding crack increases; while the length of wing crack decreases.

Figure 14 demonstrates a detailed comparison of the crack length between the numerical results and the experimental results. The length of the simulated bedding crack is consistent with the experiment, and the crack propagation angles of bedding crack and wing crack are in good agreement, respectively. There is a slight discrepancy between numerical results and experimental results in the case of bedding angel equal to 45°, since there occurs a vertical tension crack at the crack tip before entering the bedding crack. Regardless, the numerical simulation precisely captured the bedding crack and wing crack following the tensile crack. The proposed AN-PDLSM model can accurately simulate this complex oscillating crack and the length of the bedding crack that varies with the bedding angle.

5. Conclusions

In this study, peridynamic least square minimization (PDLSM) has been extended to anisotropic materials named AN-PDLSM model, based on the weighted residual method. The proposed AN-PDLSM model integrates an anisotropic model of fracture mechanics to analyze the failure process of anisotropic rocks. The rock matrix and bedding part are assumed to possess different levels of fracture toughness, and the improved maximum energy release rate criterion (MERR) is used to determine crack propagation condition and the angle of crack. In addition, the proposed AN-PDLSM model can be easily reduced to the isotropic PDLSM model. The proposed AN-PDLSM could apply boundary conditions the same as FEM and does not require surface correction and volume correction, thus improving computation accuracy. To simulate crack propagation, bonds are interrupted by a grid-independent piecewise linear model. In terms of numerical accomplishments, the proposed model uses an implicit method for quasi-static problems.

Crack propagation in isotropic plates during semi-circular bending testing is simulated based on the suggested model in order to confirm its accuracy. The numerical results show good agreement with the experimental results, confirming the effectiveness of the proposed model for analyzing fracture in isotropic rock and its ability to be reduced to the isotropic PDLSM model.

For anisotropic rocks, crack propagation in anisotropic plates under semi-circular bending tests is simulated, and the influence of the bedding weak coefficient η on the crack morphology is analyzed. The value of η has an impact on the location of a turning point. The proportion of bedding cracks increases with the decrease in the value of η, while the proportion of wing crack decreases, and vice versa. The highest agreement between the simulated and experimental results occurs when η = 0.6, which is a suitable value to characterize the weak degree of the bedding of the anisotropic rock analyzed in this paper. By analyzing the stress distribution of the failure process, it is found to be symmetrical before reaching the turning point, the stress symmetry axis will turn after the turning point. The mechanism of bedding crack and wing crack are examined by analyzing the stress intensity factor SIFs and energy release rate G, respectively. The rationality of the proposed numerical model is further verified by those analyses.

Additionally, a simulation of rock failure with various bedding angles is performed and the results are in good agreement with the experimental findings. Simple vertical cracks form at bedding angles of 0° and 90°; complex oscillating cracks appear at bedding angles of 22.5°, 45°, and 67.5°. For oscillating cracks, the crack morphology, bedding crack length, and wing crack length simulated by the model are quite similar to the experimental results. The results demonstrate the proposed AN-PDLSM model’s capacity to capture the complex crack path. In conclusion, the proposed AN-PDLSM model can accurately simulate fracture propagation in isotropic and isotropic rocks.

Author Contributions

Conceptualization, Z.Z. and Q.S.; Methodology, K.T.; Software, K.T.; Validation, K.T.; Formal analysis, N.T.; Investigation, K.T.; Data curation, K.T. and N.T.; Writing—original draft, N.T.; Writing—review & editing, Z.Z.; Visualization, N.T.; Supervision, Z.Z. and Q.S. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data are contained within the article.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. Schematics of peridynamic.

Figure 2. Anisotropic rock with weak plane and corresponding anisotropic fracture toughness Gc. (a) anisotropic rock, (b) enlarged weak planes, (c) anisotropic fracture toughness.

Figure 3. A schematic diagram of the J-integral in PDLSM. (Modified form [26]).

Figure 4. A schematic diagram of mesh-independent piecewise linear fracture growth.

Figure 5. (a) Schematic of a semi-circular bending specimen with a prefabricated notch; (b) non-uniform PD mesh.

Figure 6. Comparison of numerical results and experimental results [39] of isotropic specimens with different prefabricated notch angle α under semi-circular bending tests.

Figure 7. Computation time (CPU time) for simulation with different mesh sizes.

Figure 8. Comparison of numerical results with different mesh size. (a) mesh with 26,138 nodes, (b) mesh with 10,512 nodes, (c) mesh with 4218 nodes, (d) mesh with 1490 nodes.

Figure 9. The fracture process of a semi-circular bending specimen with different G0/Gw ratios, that is the weak coefficient η.

Figure 10. (a) The crack trajectory of a plate under a semi-circular bending test with prefabricated notch angle α = 90°, bedding angle θ = 67.5°, The weak coefficients η are, respectively, 1.0, 0.9, 0.8, 0.7, 0.6, 0.5. (b) Experimental results [40].

Figure 11. The distribution of stress component σx of different failure stage in the case of η = 0.6. (a) crack initiation, (b) reach turning point, (c) wing crack development.

Figure 12. Three typical loading stages and corresponding (a) stress intensity factor, SIFs; (b) energy release rate, G.

Figure 13. Comparison of numerical crack trajectory and experimental crack path (Wang et al., 2022) of an anisotropic semi-circular bending specimen with prefabricated notch angle α = 0° and different bedding angle θ.

Figure 14. Comparison of numerical crack length and experimental results [40] of an anisotropic semi-circular bending specimen with prefabricated notch angle α = 0° (a) bedding angle θ = 22.5° (b) bedding angle θ = 45° (c) bedding angle θ = 67.5°.

Table 1

The CPU time for simulation with different mesh sizes.

	Case 1	Case 2	Case 3	Case 4
Number of nodes	26,138	10,512	4218	1490
CPU time (s)	1230.91	411.93	147.14	50.95

References

1. Xie, Y.; Cao, P.; Jin, J.; Wang, M. Mixed mode fracture analysis of semi-circular bend (scb) specimen: A numerical study based on extended finite element method. Comput. Geotech.; 2017; 82, pp. 157-172. [DOI: https://dx.doi.org/10.1016/j.compgeo.2016.10.012]

2. Boone, T.J.; Wawrzynek, P.A.; Ingraffea, A.R. Finite element modelling of fracture propagation in orthotropic materials. Eng. Fract. Mech.; 1987; 26, pp. 185-201. [DOI: https://dx.doi.org/10.1016/0013-7944(87)90196-2]

3. Aliabadi, M.; Sollero, P. Crack growth analysis in homogeneous orthotropic laminates. Compos. Sci. Technol.; 1998; 58, pp. 1697-1703. [DOI: https://dx.doi.org/10.1016/S0266-3538(97)00240-6]

4. Motamedi, D.; Mohammadi, S. Dynamic crack propagation analysis of orthotropic media by the extended finite element method. Int. J. Fract.; 2010; 161, pp. 21-39. [DOI: https://dx.doi.org/10.1007/s10704-009-9423-7]

5. Nikolic, M.; Ibrahimbegovic, A. Rock mechanics model capable of representing initial heterogeneities and full set of 3D failure mechanisms. Comput. Methods Appl. Mech. Eng.; 2015; 290, pp. 209-227. [DOI: https://dx.doi.org/10.1016/j.cma.2015.02.024]

6. Nikolić, M.; Karavelić, E.; Ibrahimbegovic, A.; Miščević, P. Lattice Element Models and Their Peculiarities. Arch. Comput. Methods Eng. State Art Rev.; 2018; 25, pp. 753-784. [DOI: https://dx.doi.org/10.1007/s11831-017-9210-y]

7. Wang, L.; Zhou, X. Phase field model for simulating the fracture behaviors of some disc-type specimens. Eng. Fract. Mech.; 2020; 226, 106870. [DOI: https://dx.doi.org/10.1016/j.engfracmech.2020.106870]

8. Lisjak, A.; Tatone, B.S.A.; Grasselli, G.; Vietor, T. Numerical modelling of the anisotropic mechanical behaviour of opalinus clay at the laboratory-scale using FEM/DEM. Rock Mech. Rock Eng.; 2012; 47, pp. 187-206. [DOI: https://dx.doi.org/10.1007/s00603-012-0354-7]

9. Silling, S. Reformulation of elasticity theory for discontinuities and long-range forces. J. Mech. Phys. Solids; 2000; 48, pp. 175-209. [DOI: https://dx.doi.org/10.1016/S0022-5096(99)00029-0]

10. Silling, S.; Epton, M.; Weckner, O.; Xu, J.; Askari, E. Peridynamic states and constitutive modeling. J. Elast.; 2007; 88, pp. 151-184. [DOI: https://dx.doi.org/10.1007/s10659-007-9125-1]

11. Warren, T.; Silling, S.; Askari, A.; Weckner, O.; Epton, M.; Xu, J. A non-ordinary state-based peridynamic method to model solid material deformation and fracture. Int. J. Solids Struct.; 2009; 46, pp. 1186-1195. [DOI: https://dx.doi.org/10.1016/j.ijsolstr.2008.10.029]

12. Madenci, E.; Dorduncu, M.; Gu, X. Peridynamic least squares minimization. Comput. Methods Appl. Mech. Eng.; 2019; 348, pp. 846-874. [DOI: https://dx.doi.org/10.1016/j.cma.2019.01.032]

13. Liu, Q.; Xin, X. Revised non-ordinary state-based peridynamics and a new framework for coupling with finite element method. Eng. Fract. Mech.; 2021; 242, 107483. [DOI: https://dx.doi.org/10.1016/j.engfracmech.2020.107483]

14. Madenci, E.; Barut, A.; Futch, M. Peridynamic differential operator and its applications. Comput. Methods Appl. Mech. Eng.; 2016; 304, pp. 408-451. [DOI: https://dx.doi.org/10.1016/j.cma.2016.02.028]

15. Zhu, F.; Zhao, J. Peridynamic modelling of blasting induced rock fractures. J. Mech. Phys. Solids; 2021; 153, 104469. [DOI: https://dx.doi.org/10.1016/j.jmps.2021.104469]

16. Zhu, F.; Zhao, J. A Peridynamic Investigation on Crushing of Sand Particles. Géotechnique; 2018; 6, pp. 526-540. [DOI: https://dx.doi.org/10.1680/jgeot.17.P.274]

17. Xu, J.; Askari, A.; Weckner, O.; Silling, S. Peridynamic analysis of impact damage in composite laminates. J. Aerosp. Eng.; 2008; 21, pp. 187-194. [DOI: https://dx.doi.org/10.1061/(ASCE)0893-1321(2008)21:3(187)]

18. Ghajari, M.; Iannucci, L.; Curtis, P. A peridynamic material model for the analysis of dynamic crack propagation in orthotropic media. Comput. Methods Appl. Mech. Eng.; 2014; 276, pp. 431-452. [DOI: https://dx.doi.org/10.1016/j.cma.2014.04.002]

19. Zhang, H.; Qiao, P. A state-based peridynamic model for quantitative elastic and fracture analysis of orthotropic materials. Eng. Fract. Mech.; 2019; 206, pp. 147-171. [DOI: https://dx.doi.org/10.1016/j.engfracmech.2018.10.003]

20. Hattori, G.; Trevelyan, J.; Coombs, W. A non-ordinary state-based peridynamics framework for anisotropic materials. Comput. Methods Appl. Mech. Eng.; 2018; 339, pp. 416-442. [DOI: https://dx.doi.org/10.1016/j.cma.2018.05.007]

21. Tian, D.; Zhou, X. A continuum-kinematics-inspired peridynamic model of anisotropic continua: Elasticity, damage, and fracture. Int. J. Mech. Sci.; 2021; 199, 106413. [DOI: https://dx.doi.org/10.1016/j.ijmecsci.2021.106413]

22. Zeng, X.; Wei, Y. Crack deflection in brittle media with heterogeneous interfaces and its application in shale fracking. J. Mech. Phys. Solids; 2017; 101, pp. 235-249. [DOI: https://dx.doi.org/10.1016/j.jmps.2016.12.012]

23. Gao, Y.; Liu, Z.; Zeng, Q.; Wang, T.; Zhuang, Z.; Hwang, K.-C. Theoretical and numerical prediction of crack path in the material with anisotropic fracture toughness. Eng. Fract. Mech.; 2017; 180, pp. 330-347. [DOI: https://dx.doi.org/10.1016/j.engfracmech.2017.06.013]

24. Hoek, E. Fracture of anisotropic rock. J. S. Afr. Inst. Min. Metall.; 1964; 64, pp. 501-518.

25. Hussain, M.A.; Pu, S.L.; Underwood, J. Strain-Energy-Release Rate for a Crack under Combined Mode I and Mode II; ASTM: West Conshohocken, PA, USA, 1974; Volume 560, pp. 2-28.

26. Liu, Q.; Xin, X.J.; Ma, J.; Wang, Y. Simulating quasi-static crack propagation by coupled peridynamics least square minimization with finite element method. Eng. Fract. Mech.; 2021; 252, 107862. [DOI: https://dx.doi.org/10.1016/j.engfracmech.2021.107862]

27. Rice, J.R. A path independent integral and the approximate analysis of strain concentration by notches and cracks. J. Appl. Mech.; 1968; 35, pp. 379-386. [DOI: https://dx.doi.org/10.1115/1.3601206]

28. Hakim, V.; Karma, A. Crack path prediction in anisotropic brittle materials. Phys. Rev. Lett.; 2005; 95, 235501. [DOI: https://dx.doi.org/10.1103/PhysRevLett.95.235501] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16384316]

29. Nuismer, R.J. An energy release rate criterion for mixed mode fracture. Int. J. Fract.; 1975; 11, pp. 245-250. [DOI: https://dx.doi.org/10.1007/BF00038891]

30. Amestoy, M.; Leblond, J. Crack paths in plane situations—II. Detailed form of the expansion of the stress intensity factors. Int. J. Solids Struct.; 1992; 29, pp. 465-501. [DOI: https://dx.doi.org/10.1016/0020-7683(92)90210-K]

31. Saouma, V.E.; Ayari, M.L.; Leavell, D.A. Mixed mode crack propagation in homogeneous anisotropic solids. Eng. Fract. Mech.; 1987; 27, pp. 171-184. [DOI: https://dx.doi.org/10.1016/0013-7944(87)90166-4]

32. Ayari, M.L.; Ye, Z. Maximum strain theory for mixed mode crack propagation in anisotropic solids. Eng. Fract. Mech.; 1995; 52, pp. 389-400. [DOI: https://dx.doi.org/10.1016/0013-7944(95)00040-3]

33. Yau, J.; Wang, S.; Corten, H. A Mixed-Mode Crack Analysis of Isotropic Solids Using Conservation Laws of Elasticity. J. Appl. Mech. Trans. ASME; 1980; 47, pp. 335-341. [DOI: https://dx.doi.org/10.1115/1.3153665]

34. Dipasquale, D.; Sarego, G.; Zaccariotto, M.; Galvanetto, U. Dependence of crack paths on the orientation of regular 2D peridynamic grids. Eng. Fract. Mech.; 2016; 160, pp. 248-263. [DOI: https://dx.doi.org/10.1016/j.engfracmech.2016.03.022]

35. Bažant, Z.P.; Nguyen, H.T.; Donmez, A.A. Critical Comparison of Phase-Field, Peridynamics, and Crack Band Model M7 in Light of Gap Test and Classical Fracture Tests. J. Appl. Mech.; 2022; 89, 061008. [DOI: https://dx.doi.org/10.1115/1.4054221]

36. Ni, T.; Zaccariotto, M.; Zhu, Q.-Z.; Galvanetto, U. Static solution of crack propagation problems in Peridynamics. Comput. Methods Appl. Mech. Eng.; 2019; 346, pp. 126-151. [DOI: https://dx.doi.org/10.1016/j.cma.2018.11.028]

37. Huang, D.; Lu, G.; Wang, C.; Qiao, P. An extended peridynamic approach for deformation and fracture analysis. Eng. Fract. Mech.; 2015; 141, pp. 196-211. [DOI: https://dx.doi.org/10.1016/j.engfracmech.2015.04.036]

38. Kilic, B.; Madenci, E. An adaptive dynamic relaxation method for quasi-static simulations using the peridynamic theory. Theor. Appl. Fract. Mech.; 2010; 53, pp. 194-204. [DOI: https://dx.doi.org/10.1016/j.tafmec.2010.08.001]

39. Ayatollahi, M.R.; Aliha MR, M.; Hasani, M.M. Mixed mode brittle fracture in PMMA—An experimental study using SCB specimens. Mater. Sci. Eng. A; 2006; 417, pp. 348-356. [DOI: https://dx.doi.org/10.1016/j.msea.2005.11.002]

40. Wang, W.; Teng, T. Experimental study on anisotropic fracture characteristics of coal using notched semi-circular bend specimen. Theor. Appl. Fract. Mech.; 2022; 122, 103559. [DOI: https://dx.doi.org/10.1016/j.tafmec.2022.103559]

41. Ju, M.; Li, J.; Li, J.; Zhao, J. Loading rate effects on anisotropy and crack propagation of weak bedding plane-rich rocks. Eng. Fract. Mech.; 2020; 230, 106983. [DOI: https://dx.doi.org/10.1016/j.engfracmech.2020.106983]