# Modeling of Shear Crack Propagation in Rock Masses Using Mesh-Free LRPIM

Authors: Qingbo Li; Nengxiong Xu; Weifeng Wan; Yazhe Li
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

The modeling of shear cracks in materials is critical in various engineering applications, such as the safety analysis of concrete structures and stability analysis of rock slopes. Based on the idea of Goodman element, the elastic-plastic constitutive model of the shear cracks is derived, and the elastic-plastic analysis of shear crack propagation is realized in the local radial basis point interpolation method (LRPIM). This method avoids the loss of accuracy caused by the mesh in the analysis of fracture propagation, and the crack propagation of rock brittle material is simulated. The investigation indicates that (1) the LRPIM results are close to the FDM results, which demonstrates that it is feasible to analyze shear cracks in rock masses. (2) Compared with the results of the built-in oblique crack model, when the LRPIM is used to analyze crack propagation, the results are close to the experimental results, showing that the LRPIM can model shear crack propagation in a rock mass. (3) The propagation path using the LRPIM is not sufficiently smooth, which can be explained as the crack tip stress and strain not being sufficiently accurate and still requiring further improvement.


---

## Body

## 1. Introduction

Rock masses are discontinuous media that contain many kinds of defects, such as cracks, bedding, and holes. Existing shear cracks have a substantial influence on the mechanical characteristics of rock masses. A rock mass with many defects is extremely sensitive to the action of forces. Therefore, it is necessary to analyze the influence of shear cracks and their propagation under the action of external forces.Some research results have been obtained through physical experiments. The mechanism of crack propagation in concrete or rock specimens containing cracks under shear loading conditions is studied [1], and ultrasonic waves recorded during direct shear experiments on rock joints were employed to investigate the shear failure processes. [2].Some methods that do not need too many criteria or special elements have been proposed [3, 4], which have great potential in engineering applications. Peridynamics (PD) describes the mechanical behavior of matter by solving the spatial integral equation, which avoids the singularity and complexity of the traditional model based on continuity assumption when facing discontinuity problems. It has been successfully applied to the simulation of various discontinuity problems with different scales and has become a research hotspot in the current international computational mechanics and related fields. In the simulation process of near-field dynamics model, the phenomena of crack initiation, propagation, and bifurcation are naturally generated with the simulation process. In addition to adding a single-bond fracture condition, there is no need for too many judgment conditions or special elements [5, 6].Scholars have proposed various theories and calculation methods to study the discontinuous characteristics of rock masses. The main numerical methods for simulating crack propagation are the finite element method (FEM) [7], extended finite element method (XFEM) [8], and meshless method [9]. The FEM was applied to study the mechanical behavior of rock masses and concrete in references [10, 11]. When the FEM is used to study crack propagation, it is necessary to remesh the elements, which substantially increases the workload and leads to low calculation accuracy and efficiency. Meshing and remeshing has significantly improved in the past years. An alternative crack propagation algorithm has been proposed which effectively circumvents the variable transfer procedure adopted with classical mesh adaptation algorithms [12]. And, neural networks can be used to compute the crack growth in fracture specimens, which finally leads to a prediction of classical fracture toughness parameters [13]. These shortcomings can also be avoided by the XFEM proposed by Belytschko and Moes [8, 9, 14]. The XFEM rationally realizes accurate crack descriptions and has become a popular research topic among numerical methods for crack problems [15].Moes et al. [9] used the XFEM to simulate bond crack propagation with the stress intensity factor as the crack index. Mariani et al. [16] used the XFEM to simulate quasistatic bond crack propagation in quasibrittle materials. Deb and Das [17] combined the XFEM with cohesion to simulate shear crack slip. Zhou and Yang [18] proposed a multiscale XFEM to study crack propagation and rock bridge penetration under multidirectional compression, and the conclusions were in good agreement with the experimental results.Compared with the FEM and XFEM, the meshless method is based on field nodes that can eliminate or partially eliminate the difficulties caused by meshing. The interpolation (approximation) function of the meshless method does not depend on the elements, so it has the advantages of easy higher order approximation formation and realization of local node encryption, bringing convenience to the numerical analysis of cracks.Belytschko [19] first used the EFGM method to analyze crack propagation. Krysl [20] further extended the EFGM method to crack propagation calculations in a three-dimensional model. Belytschko [21] applied the meshless method to simulate dynamic crack propagation in concrete structures, and they simulated the crack propagation direction and velocity in concrete structures very well. Zhuang Xiaoying et al. [22] realized accurate fracture modeling in a two-dimensional model using meshless methods, the visibility criterion, and level sets.The above applications all use the meshless EFGM, but this method requires the background grid for integration and is not a true meshless method. The LRPIM [23] method is a local weak computing method that does not require the background grid and is a real meshless method. Little research has been conducted on crack problems with the LRPIM. Three-dimensional shear crack propagation with the LRPIM has not yet been studied. A three-dimensional shear crack propagation calculation is established based on elastic-plastic theory and fracture mechanics. The reliability of the method is verified through numerical experiments, and the problems that need to be improved are discussed.This paper is organized as follows. Section2 briefly introduces the background of the LRPIM. Section 3 mainly concentrates on the methods of modeling shear cracks with the LRPIM. Section 4 discusses the results. Finally, Section 5 summarizes the conclusions.


## 2. Background: LRPIM

If only the static case is considered, the BVP of three-dimensional solid mechanics can be expressed as(1)σij,j+bi=0,σijnj=ti¯inΓt,ui=u¯iinΓu,ux,t0=u0x,x∈Ω,where u¯i, ti¯, and u0 are the given displacement, the given surface force, and the initial displacement. I, j = (1,2,3) represents x, y, z.Considering the solid mechanics problem defined in domain Ω, the local weighted residual method is used for node I to satisfy the governing equation, and the local weak equation of the node is obtained. The form of the locally weighted residuals is defined on the local integral domainΩqand the corresponding boundary Γq in the following form(2)∫ΩqW^Iσij,j+bidΩ=0,where W^i is a weight function or test function centered at field node I.The follow formula is obtained by integral operation:(3)∫ΩqW^Iσij,jdΩ=∫ΓqW^ItidΓ−∫ΩqW^I,jσijdΩ,where ti is the component of surface force vector on the boundary.We bring formula (10) into formula (6),(4)∫ΓqW^ItiσijdΓ−∫ΩqW^I,jσij−W^IbidΩ=0,∫ΓqiW^ItidΓ+∫ΓquW^ItidΓ+∫ΓqtW^ItidΓ−∫ΩqW^I,jσij−W^IbidΩ=0,where Γqi is the inner boundary of the integral domain; Γqu is the part of the natural boundary that intersects the integral domain; and Γqt is the part of the essential boundary that intersects the integral domainIn the local weak meshless method based on RPIM, the radial base point interpolation shape function is used to obtain the displacement of the calculated point, and the strain and stress vectors are further obtained.(5)ε6∗1=B6∗3nu3n∗1,(6)σ6∗1=D6∗6ε6∗1=D6∗6B6∗3nu3n∗1,B is a differential operator, D is the material constant matrix.The system equation is discretized into matrix form, and the following formula is obtained.(7)∫ΩqV^ITσdΩ−∫ΓqiW^ItdΓ−∫ΓquW^ItdΓ=∫ΓqtW^It¯dΓ+∫ΩqW^ITbdΩ,where W^I is the weight function matrix and V^IT is the derivative matrix of weight function.Bring formulas (5) and (6) into formula (7):(8)∫ΩqV^ITDBudΩ−∫ΓqiW^ITnDBudΓ−∫ΓquW^ITnDBudΓ=∫ΓqtW^It¯dΓ+∫ΩqW^ITbdΩ.The matrix form can be expressed as(9)KI3∗3nu3n∗1=fI3∗1,where KI is the stiffness matrix of the field node I.(10)KI=∫ΩqV^ITDBdΩ−∫ΓqiW^ITnDBdΓ−∫ΓquW^ITnDBdΓ,where fI is the nodal force vector, including the physical force in the problem domain and the surface force on the natural boundary.(11)fI=∫ΓqtW^It¯dΓ+∫ΩqW^ITbdΩ.The standard gauss integral and the rectangular local integral domain are used. The final system equations can be obtained by assembling all these equations based on the system of population numbers.(12)K3N∗3NU3N∗1=F3N∗1.For crack expansion, the relationship between the nodes and Gauss integral points is constantly changing. Therefore, the visual criterion is used to address the range of the influence domain. If the line between the field node and the Gauss integral point (either a solid boundary or crack) intersects, the node is considered to be “masked” After assembling the total rigid matrix, a penalty function is used to apply the displacement boundary conditions.


## 3. Methods

### 3.1. Overview

A three-dimensional elastic-plastic local radial basis point interpolation method for calculating crack propagation is established. The key techniques used in this section are described in detail. (1) According to the principle of the minimum potential energy, the governing equation of the shear crack action is derived based on the idea of the Goodman element, which is commonly used in finite element analysis. (2) The elastic-plastic constitutive model of the shear cracks is derived based on the elastoplastic theory. (3) According to crack mechanics theory, the criterion of shear crack propagation is obtained, which can be used for shear crack propagation.


### 3.2. Derivation of the Governing Equation of the Shear Crack

For the contact problem of shear cracks, the materials on both sides of the crack cannot be embedded into each other. However, to calculate the interaction between cracks, the idea of the Goodman element, which is often used in the FEM, is used to address this problem [19]. It is assumed that there will be a small amount of overlap between the materials on both sides of the crack, and a relatively large normal stiffness is set. When there is a small amount of overlap between the materials above and below the crack, a large normal force will occur. At the same time, the shear stiffness coefficient is set to reflect the tangential force generated by the rubbing of the crack.The crack profile is shown in Figure1. The red line is the crack profile, the hollow dot is the field node, and the solid point is on the crack. It should be pointed out that the points on the crack are not directly involved in the calculation of the model. Taking a point on the crack as the origin, the local coordinate system is established. The normal direction of the crack is along the y′ axis, and the two tangent directions are along the x′ axis and z′ axis. Two points are set up in the positive and negative directions of the y′ axis, and there is a small distance from the origin, which represent the material above and below the crack, respectively.Figure 1
Schematic diagram of the local coordinate system.The energy equation considering crack interaction is shown in equation (13), which represents the work.(13)Πc=12∫ΓΔuTσdΓ,where Δu is the difference in the normal displacement of the materials on both sides of the crack. If Δu is negative, there is an interaction between the two sides of the crack, and it is necessary to calculate the stress-strain effect of this action. Δv and Δw represent the difference in the displacement between the two shear directions in the local coordinate system.(14)Δu=ΔuΔvΔw=uiviwi−ujvjwj,where ui, vi, wi, uj, vj, and wj represent the displacements of points i and j in the local coordinate system. The displacement value can be obtained by interpolation from the displacement information of the field node. In this paper, the radial base point interpolation method is used as follows:(15)uiviwi=Φiu,where u is the displacement vector of the field node in the support domain of the interpolation node, and the vector length is 3n.(16)u=u1iv1iw1iu2iv2iw2i…univniwni,and Φi is the displacement shape function of node i in the following form:(17)Φi=ϕ1i00ϕ2i00…ϕni000ϕ1i00ϕ2i0…0ϕni000ϕ1i00ϕ2i…00ϕni.The shape function and displacement vector should be integrated when solving the problem. Here, we should pay attention to the corresponding point number. We need to expand the form as follows:(18)Φ=ϕ1i00ϕ2i00…ϕni00−ϕ1j00−ϕ2j00…−ϕnj000ϕ1i00ϕ2i0…0ϕni00−ϕ1j00−ϕ2j0…0−ϕnj000ϕ1i00ϕ2i…00ϕni00−ϕ1j00−ϕ2j…0−ϕnj0.Here,σ is set as the node stress vector and D is the material matrix; then, σ(19)σ=τsσnτn=D×Δu,σn is the normal force, and τs and are the component forces in the two tangent directions.kn and ks are the normal and tangential stiffnesses, respectively.(20)D=ks000kn000ks.It is necessary to transform the displacement of the nodes from the global coordinate system to the local coordinate system.R is a transformation matrix that converts global coordinates to local coordinates.(21)R=l1m1n1l2m2n2l3m3n3.The displacement in the local coordinate system is obtained as follows:(22)u′=RΦu.A functional considering the crack energy is established based on the principle of the minimum potential energy. Its form is shown in the following formula:(23)Πc=12∫ΓΔu′TσdΓ,Πc=12∫ΓΔu′TDΔudΓ.By variational calculation,(24)δΠc=δuT∫ΓΦTRTDRΦdΓu=0,KC=∫ΓΦTRTDRΦdΓ,KC is the stiffness matrix reflecting the action of the crack. When there are no repeated calculation points in the field node of the crack, the action of the crack is reflected by the square matrix of 3m+n⋅3m+n. The stiffness matrix can be combined into the stiffness matrix established by the local weak meshless method. In this process, we should pay attention to the corresponding relationship between the nodes.The matrix form of the crack surface energy functional is obtained as follows:(25)KC=∫ΓΦTRTDRΦdΓ=AiΦTRTDRΦ.


### 3.3. Elastoplastic Theory of the Shear Crack

Limited by the strength of the crack, when the shear stress of the crack exceeds the maximum value, the shear stress of the crack should not continue to increase. It is assumed that the crack is a smooth surface, which accords with the deformation characteristics of ideal plasticity.The stress state of the point on the crack can be decomposed into three directions of stress, as shown in Figure2. A local coordinate system is established at a node on the crack, and the normal stress of the crack is σn. If the resultant force of the tangential force is expressed as τs, then the value of the shear stress in the direction perpendicular to σn and is zero. It is assumed here that the direction of τs is the same as that of the relative displacement of the crack.Figure 2
Local coordinate system.According to the Mohr–Coulomb strength criterion, the general form of the yield function is as follows:(26)f=τs−c−σntanφ,where c and ϕ are the cohesion and internal friction of the crack, τs is the shear stress of the crack, and σn is the normal direction of the crack.(27)∂f∂σ=∂f∂τs∂f∂σn∂f∂τnT=1−tanϕ0T.The elastic matrix of the stress displacement relationship can be written as follows.D=ks000kn000ksIn elastoplastic theory, the deformation matrix includes the elastic modulus and Poisson’s ratio in kPa. The stress is obtained by multiplying the matrix and the strain vector. Here, the stiffness coefficient is in the deformation matrix, and its unit is kPa/m, which corresponds to the displacement difference between the upper and lower nodes of the crack; the unit ism, that is, the strain is replaced by the displacement.The plastic matrix is as follows:(28)Dp=D∂g/∂σ∂f/∂σTDA+∂f/∂σTD∂g/∂σ.Under the ideal elastic-plastic condition,A = 0. An adaptive flow criterion is adopted as g=f. The elastic-plastic deformation matrix of the shear crack surface is as follows:(29)Dp=D∂f/∂σ∂f/∂σTDA+∂f/∂σTD∂f/∂σ=1S0ks2ksS10ksS1S120000,S0=ks+kntan2ϕ,S1=−kntanϕ.The elastic-plastic matrix after yielding is as follows:(30)Dep=D−Dp,Dep=ks−1S0ks2−1S0ksS10−1S0ksS1kn−1S0S12000ks.The deformation matrix under ideal elastoplastic conditions is only related to the stiffness coefficient and internal friction angle.From the above derivation, we can obtain the constitutive relation of the shear cracks under ideal plastic conditions. Before yielding, the deformation matrix isD, and after yielding, the elastic-plastic deformation matrix is Dep. In this way, the yield criterion and elastic-plastic constitutive relation of the cracks are established.In the calculation, when the crack is open, the normal stress is zero, the crack is a free surface, and the contact effect between the cracks is not considered. When the normal stress of the crack is not zero, it indicates that the cracks are in contact with each other and that there is a force effect. In the solution, the stress state of the point on the crack is reflected by modifying the deformation matrix.


### 3.4. Crack Propagation Criterion

In this paper, the maximum circumferential tensile stress theory is selected to analyze crack propagation. For the I-II tension shear composite mode crack [24], under the plane strain state, the critical conditions are as follows:(31)kIC=cosθ02kIcos2θ02−32kIIsinθ0.The expansion step can be set to a fixed value, which is a common setting method. After determining the crack criterion, we can judge whether the crack tip is cracked and determine the direction of the crack. According to previous work, this paper simulates three-dimensional crack propagation by expanding a finite number of crack front points in their respective crack front normal planes.The crack is simulated by a group of small triangles, and a group of line segments are used to simulate the crack front. When the model is balanced, the current stress-strain state of the nodes in the calculation domain can be obtained, and the propagation of the crack can be calculated. As shown in Figure3, the stress states of points B, C, and E in the crack front edge in Figure 3(a) meet the expansion criteria and expand in a certain direction, and a new crack tip endpoint is established in its propagation direction. In diagram 3(b), two new triangular planes are formed by the new tip point and two points A and C, which realize the first step expansion of the crack. The extensions of points C and E are realized by the same method, and the new crack is reconstructed, as shown in Figure 3(d). With the expansion of the crack, the relationship between the calculation point and the field node will constantly change. It is necessary to recalculate the stiffness matrix, solve the equilibrium equation, and then carry out the next step of the expansion calculation.Figure 3
Geometric sketch of crack propagation [25].
(a)(b)(c)(d)


### 3.5. Calculation of the Shear Crack Propagation in a Rock Mass

Combining the solution of the local weak meshless method, crack propagation, and shear crack interaction, the program is compiled according to the discrete rectangular form. The calculation process is as follows:(1)
The initial conditions are set, including the arrangement of the nodes, the definition of the material parameters, the application of the boundary conditions, and the application of the loads.(2)
Applying thenth-order load, the substiffness matrix of the joint and the stiffness matrix of the shear crack are calculated. The stress-strain state of the whole model is obtained by solving the equilibrium convergence of several iterative steps.(3)
After entering the expansion step, we can judge whether the crack is expanding and determine the direction of the crack propagation according to the stress state at the crack tip. If the crack tip meets the crack criterion, crack propagation is realized, and the iterative calculation is carried out again in step 2. If the crack does not expand, then proceed to the next step.(4)
Judge whether level N load has been applied. If so, terminate the calculation; otherwise, return to step ②.The calculation module of shear crack propagation is integrated with the meshless LRPIM program.


## 3.1. Overview

A three-dimensional elastic-plastic local radial basis point interpolation method for calculating crack propagation is established. The key techniques used in this section are described in detail. (1) According to the principle of the minimum potential energy, the governing equation of the shear crack action is derived based on the idea of the Goodman element, which is commonly used in finite element analysis. (2) The elastic-plastic constitutive model of the shear cracks is derived based on the elastoplastic theory. (3) According to crack mechanics theory, the criterion of shear crack propagation is obtained, which can be used for shear crack propagation.


## 3.2. Derivation of the Governing Equation of the Shear Crack

For the contact problem of shear cracks, the materials on both sides of the crack cannot be embedded into each other. However, to calculate the interaction between cracks, the idea of the Goodman element, which is often used in the FEM, is used to address this problem [19]. It is assumed that there will be a small amount of overlap between the materials on both sides of the crack, and a relatively large normal stiffness is set. When there is a small amount of overlap between the materials above and below the crack, a large normal force will occur. At the same time, the shear stiffness coefficient is set to reflect the tangential force generated by the rubbing of the crack.The crack profile is shown in Figure1. The red line is the crack profile, the hollow dot is the field node, and the solid point is on the crack. It should be pointed out that the points on the crack are not directly involved in the calculation of the model. Taking a point on the crack as the origin, the local coordinate system is established. The normal direction of the crack is along the y′ axis, and the two tangent directions are along the x′ axis and z′ axis. Two points are set up in the positive and negative directions of the y′ axis, and there is a small distance from the origin, which represent the material above and below the crack, respectively.Figure 1
Schematic diagram of the local coordinate system.The energy equation considering crack interaction is shown in equation (13), which represents the work.(13)Πc=12∫ΓΔuTσdΓ,where Δu is the difference in the normal displacement of the materials on both sides of the crack. If Δu is negative, there is an interaction between the two sides of the crack, and it is necessary to calculate the stress-strain effect of this action. Δv and Δw represent the difference in the displacement between the two shear directions in the local coordinate system.(14)Δu=ΔuΔvΔw=uiviwi−ujvjwj,where ui, vi, wi, uj, vj, and wj represent the displacements of points i and j in the local coordinate system. The displacement value can be obtained by interpolation from the displacement information of the field node. In this paper, the radial base point interpolation method is used as follows:(15)uiviwi=Φiu,where u is the displacement vector of the field node in the support domain of the interpolation node, and the vector length is 3n.(16)u=u1iv1iw1iu2iv2iw2i…univniwni,and Φi is the displacement shape function of node i in the following form:(17)Φi=ϕ1i00ϕ2i00…ϕni000ϕ1i00ϕ2i0…0ϕni000ϕ1i00ϕ2i…00ϕni.The shape function and displacement vector should be integrated when solving the problem. Here, we should pay attention to the corresponding point number. We need to expand the form as follows:(18)Φ=ϕ1i00ϕ2i00…ϕni00−ϕ1j00−ϕ2j00…−ϕnj000ϕ1i00ϕ2i0…0ϕni00−ϕ1j00−ϕ2j0…0−ϕnj000ϕ1i00ϕ2i…00ϕni00−ϕ1j00−ϕ2j…0−ϕnj0.Here,σ is set as the node stress vector and D is the material matrix; then, σ(19)σ=τsσnτn=D×Δu,σn is the normal force, and τs and are the component forces in the two tangent directions.kn and ks are the normal and tangential stiffnesses, respectively.(20)D=ks000kn000ks.It is necessary to transform the displacement of the nodes from the global coordinate system to the local coordinate system.R is a transformation matrix that converts global coordinates to local coordinates.(21)R=l1m1n1l2m2n2l3m3n3.The displacement in the local coordinate system is obtained as follows:(22)u′=RΦu.A functional considering the crack energy is established based on the principle of the minimum potential energy. Its form is shown in the following formula:(23)Πc=12∫ΓΔu′TσdΓ,Πc=12∫ΓΔu′TDΔudΓ.By variational calculation,(24)δΠc=δuT∫ΓΦTRTDRΦdΓu=0,KC=∫ΓΦTRTDRΦdΓ,KC is the stiffness matrix reflecting the action of the crack. When there are no repeated calculation points in the field node of the crack, the action of the crack is reflected by the square matrix of 3m+n⋅3m+n. The stiffness matrix can be combined into the stiffness matrix established by the local weak meshless method. In this process, we should pay attention to the corresponding relationship between the nodes.The matrix form of the crack surface energy functional is obtained as follows:(25)KC=∫ΓΦTRTDRΦdΓ=AiΦTRTDRΦ.


## 3.3. Elastoplastic Theory of the Shear Crack

Limited by the strength of the crack, when the shear stress of the crack exceeds the maximum value, the shear stress of the crack should not continue to increase. It is assumed that the crack is a smooth surface, which accords with the deformation characteristics of ideal plasticity.The stress state of the point on the crack can be decomposed into three directions of stress, as shown in Figure2. A local coordinate system is established at a node on the crack, and the normal stress of the crack is σn. If the resultant force of the tangential force is expressed as τs, then the value of the shear stress in the direction perpendicular to σn and is zero. It is assumed here that the direction of τs is the same as that of the relative displacement of the crack.Figure 2
Local coordinate system.According to the Mohr–Coulomb strength criterion, the general form of the yield function is as follows:(26)f=τs−c−σntanφ,where c and ϕ are the cohesion and internal friction of the crack, τs is the shear stress of the crack, and σn is the normal direction of the crack.(27)∂f∂σ=∂f∂τs∂f∂σn∂f∂τnT=1−tanϕ0T.The elastic matrix of the stress displacement relationship can be written as follows.D=ks000kn000ksIn elastoplastic theory, the deformation matrix includes the elastic modulus and Poisson’s ratio in kPa. The stress is obtained by multiplying the matrix and the strain vector. Here, the stiffness coefficient is in the deformation matrix, and its unit is kPa/m, which corresponds to the displacement difference between the upper and lower nodes of the crack; the unit ism, that is, the strain is replaced by the displacement.The plastic matrix is as follows:(28)Dp=D∂g/∂σ∂f/∂σTDA+∂f/∂σTD∂g/∂σ.Under the ideal elastic-plastic condition,A = 0. An adaptive flow criterion is adopted as g=f. The elastic-plastic deformation matrix of the shear crack surface is as follows:(29)Dp=D∂f/∂σ∂f/∂σTDA+∂f/∂σTD∂f/∂σ=1S0ks2ksS10ksS1S120000,S0=ks+kntan2ϕ,S1=−kntanϕ.The elastic-plastic matrix after yielding is as follows:(30)Dep=D−Dp,Dep=ks−1S0ks2−1S0ksS10−1S0ksS1kn−1S0S12000ks.The deformation matrix under ideal elastoplastic conditions is only related to the stiffness coefficient and internal friction angle.From the above derivation, we can obtain the constitutive relation of the shear cracks under ideal plastic conditions. Before yielding, the deformation matrix isD, and after yielding, the elastic-plastic deformation matrix is Dep. In this way, the yield criterion and elastic-plastic constitutive relation of the cracks are established.In the calculation, when the crack is open, the normal stress is zero, the crack is a free surface, and the contact effect between the cracks is not considered. When the normal stress of the crack is not zero, it indicates that the cracks are in contact with each other and that there is a force effect. In the solution, the stress state of the point on the crack is reflected by modifying the deformation matrix.


## 3.4. Crack Propagation Criterion

In this paper, the maximum circumferential tensile stress theory is selected to analyze crack propagation. For the I-II tension shear composite mode crack [24], under the plane strain state, the critical conditions are as follows:(31)kIC=cosθ02kIcos2θ02−32kIIsinθ0.The expansion step can be set to a fixed value, which is a common setting method. After determining the crack criterion, we can judge whether the crack tip is cracked and determine the direction of the crack. According to previous work, this paper simulates three-dimensional crack propagation by expanding a finite number of crack front points in their respective crack front normal planes.The crack is simulated by a group of small triangles, and a group of line segments are used to simulate the crack front. When the model is balanced, the current stress-strain state of the nodes in the calculation domain can be obtained, and the propagation of the crack can be calculated. As shown in Figure3, the stress states of points B, C, and E in the crack front edge in Figure 3(a) meet the expansion criteria and expand in a certain direction, and a new crack tip endpoint is established in its propagation direction. In diagram 3(b), two new triangular planes are formed by the new tip point and two points A and C, which realize the first step expansion of the crack. The extensions of points C and E are realized by the same method, and the new crack is reconstructed, as shown in Figure 3(d). With the expansion of the crack, the relationship between the calculation point and the field node will constantly change. It is necessary to recalculate the stiffness matrix, solve the equilibrium equation, and then carry out the next step of the expansion calculation.Figure 3
Geometric sketch of crack propagation [25].
(a)(b)(c)(d)


## 3.5. Calculation of the Shear Crack Propagation in a Rock Mass

Combining the solution of the local weak meshless method, crack propagation, and shear crack interaction, the program is compiled according to the discrete rectangular form. The calculation process is as follows:(1)
The initial conditions are set, including the arrangement of the nodes, the definition of the material parameters, the application of the boundary conditions, and the application of the loads.(2)
Applying thenth-order load, the substiffness matrix of the joint and the stiffness matrix of the shear crack are calculated. The stress-strain state of the whole model is obtained by solving the equilibrium convergence of several iterative steps.(3)
After entering the expansion step, we can judge whether the crack is expanding and determine the direction of the crack propagation according to the stress state at the crack tip. If the crack tip meets the crack criterion, crack propagation is realized, and the iterative calculation is carried out again in step 2. If the crack does not expand, then proceed to the next step.(4)
Judge whether level N load has been applied. If so, terminate the calculation; otherwise, return to step ②.The calculation module of shear crack propagation is integrated with the meshless LRPIM program.


## 4. Results: Case Studies

The correctness and stability of the calculation method are verified by the following three examples. Example 1 verifies the calculation method and realizes the force transfer effect of cracks, and case 2 verifies the influence of oblique shear cracks on the deformation characteristics of the model under different tangential stiffnesses to verify the correctness of the calculation results. In case 3, the propagation path of the built-in inclined crack under axial compression is obtained.


### 4.1. Model Information

The information for the example model is as follows: the dimensions are50mm×50mm×100mm, the deformation parameters include an elastic modulus of 1.6 GPa, and a Poisson’s ratio of 0.25; the upper and lower surfaces of the model are fixed, and a vertical displacement is applied on the upper surface; and the displacement value is 1 mm.According to whether the interaction between the cracks is considered, the experimental calculation is carried out in two groups. The first experiment does not use the developed calculation module, and the other test uses the developed calculation module. The normal stiffness is set to 20 times the elastic modulus, and the tangential stiffness is set to 0.1 times the elastic modulus. The node distribution of the model is shown in Figure4.Figure 4
Model and field node distribution diagram.If the interaction between cracks is not considered, the calculated displacement nephogram is shown in Figure5(a). The upper part of the model has a consistent vertical displacement value of 1 mm. Because the interaction between the cracks is not considered, the upper part of the model cannot transfer the force to the lower part, so the displacement value of the lower part is 0 mm because there is no force on the lower part. This calculation result is consistent with the expected results.Figure 5
Vertical displacement nephogram (m).
(a)(b)If the interaction of cracks is considered, the vertical displacement nephogram of the model is obtained, as shown in Figure5(b). The displacement of the model presents a regular ladder shape, and the force between the cracks is effectively transferred. This shows that the local weak meshless method considering the interaction between crack surfaces is effective, and the established calculation method can effectively calculate the interaction between the two sides of the crack.


### 4.2. Oblique Crack Model

#### 4.2.1. Elastic Crack Model

There is an oblique through crack in the middle of the model in Figure6, which divides the model into two parts. The parameters and boundary conditions of the model are consistent with those of test 1, but the crack is oblique. The normal stiffness is 20 times the elastic modulus, and the tangential stiffness changes. The correctness of the algorithm is verified by setting different tangential stiffness coefficients. Five test models are set, and the tangential stiffness is set to 0 times, 0.05 times, 0.1 times, 0.15 times, and 0.2 times the elastic modulus.Figure 6
Schematic diagram of oblique through crack.It should be noted that the crack is implicit in this algorithm. The results do not include crack information, only field node information. The field node at the crack of the displacement nephogram is not on the crack surface. In addition, the value shown in the legend is not the maximum and minimum value of the cloud map but the segmented value. The maximum values of the calculated results are compared in Table1.Table 1
Horizontal displacement of the model.
Tangential stiffness (Pa)0.0E0.05E0.1E0.15E0.2EMaximum horizontal displacement (mm)LRPIM7.069e − 26.829e − 26.604e − 26.393e − 26.195e − 2FDM6.997e − 26.887e − 26.782e − 26.682e − 26.585e − 2Error (%)1.00.82.64.35.9The horizontal displacement nephograms of the five test results are shown in Figure7. The maximum horizontal displacements are approximately 7.069e-2 mm, 6.829e-2 mm, 6.604e-2 mm, 6.393e-2 mm, and 6.195e-2 mm. It can be understood that, due to the increasing tangential stiffness, the deformation of the model on both sides of the crack is constrained, resulting in the decrease in the lateral deformation, which is consistent with the expected calculation results.Figure 7
Horizontal displacement nephogram (m).
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)As shown in Table1, the maximum horizontal displacement decreased with increasing tangential stiffness, and the calculation results obtained by the two algorithms were very close.The above derivation and calculation results assume that, under elastic conditions, the material and crack have no yield criterion or constitutive relation under elastic-plastic conditions. The model considering the plastic deformation of cracks will be calculated and analyzed in the following section.


#### 4.2.2. Elastoplastic Crack Model

The interface elements in FLAC3D also adopt the Mohr–Coulomb strength criterion and ideal plastic constitutive model. The contact surface parameters include the normal stiffness, tangential stiffness, cohesion, internal friction angle, tensile strength, and expansion angle. The tensile strength and expansion angle of the contact surface are not considered in this paper.The crack is set as elastic-plastic, the cohesion is 0 kPa, and the friction coefficient is 0.01. A normal stiffness of 20 and tangential stiffness of 0.1 times the elastic modulus are set.The calculation results of the LRPIM and FLAC3D are shown in Figure8. The results of the two numerical examples have the same displacement nephogram. As shown in Table 2, under ideal plastic conditions, the maximum horizontal displacement of the model using the LRPIM is 7.051e − 2 mm. The results of the two methods are close, and the error is only 0.8%.Figure 8
Horizontal displacement nephogram (m). (a) FLAC3D. (b) LRPIM.
(a)(b)Table 2
Horizontal displacement comparison.
MethodLRPIMFLAC3DMaximum horizontal displacement (mm)7.051e − 26.995e − 2Error (%)0.8


### 4.3. Built-In Inclined Crack Model

For shear cracks, the stress intensity factor isKI=0 because there is no relative displacement in the normal direction. According to the maximum circumferential stress theory, θ=70.5°.Under the action of axial pressure, the cracks of brittle rock with built-in inclined cracks will expand, and the crack propagation pattern is shown in Figure9 [26].Figure 9
Crack propagation under uniaxial compression [26].The model parameters are the same as the parameters of common marble. The model information is set as follows: dimensions of 62 mm × 25 mm × 110 mm, elastic modulus of 17 GPa, Poisson’s ratio of 0.25, internal friction angle of 30°, and fixed displacement of the upper and lower surfaces of the model. The angle between the crack and the horizontal plane is 45°, and the crack length is 20 mm. The cohesive force of the model material is 20.9 kPa, the cohesive force of the crack is 2 kPa, the friction coefficient is 0.01, the normal stiffness is twice the elastic modulus, and the tangential stiffness is 0.01 times the elastic modulus. The force diagram of the model is shown in Figure10.Figure 10
Schematic diagram of cracks subjected to uniaxial compression.A vertical displacement of 0.2 mm is applied on the upper surface of the model to simulate the stress state of compression and promote the expansion of cracks. Figure11 shows the vertical displacement nephogram and stress nephogram in the initial state. There is discontinuity in the model, while the vertical stress at the crack tip is obviously larger.Figure 11
Nephogram of the model subjected to unidirectional compression. (a) Vertical displacement (m). (b) Vertical stress (Pa).
(a)(b)Figure12 shows the crack morphology calculated by the LRPIM, and the two ends of the crack are curved and extend upward, which demonstrates the propagation process of the wing crack.Figure 12
Crack propagation under uniaxial compression.


## 4.1. Model Information

The information for the example model is as follows: the dimensions are50mm×50mm×100mm, the deformation parameters include an elastic modulus of 1.6 GPa, and a Poisson’s ratio of 0.25; the upper and lower surfaces of the model are fixed, and a vertical displacement is applied on the upper surface; and the displacement value is 1 mm.According to whether the interaction between the cracks is considered, the experimental calculation is carried out in two groups. The first experiment does not use the developed calculation module, and the other test uses the developed calculation module. The normal stiffness is set to 20 times the elastic modulus, and the tangential stiffness is set to 0.1 times the elastic modulus. The node distribution of the model is shown in Figure4.Figure 4
Model and field node distribution diagram.If the interaction between cracks is not considered, the calculated displacement nephogram is shown in Figure5(a). The upper part of the model has a consistent vertical displacement value of 1 mm. Because the interaction between the cracks is not considered, the upper part of the model cannot transfer the force to the lower part, so the displacement value of the lower part is 0 mm because there is no force on the lower part. This calculation result is consistent with the expected results.Figure 5
Vertical displacement nephogram (m).
(a)(b)If the interaction of cracks is considered, the vertical displacement nephogram of the model is obtained, as shown in Figure5(b). The displacement of the model presents a regular ladder shape, and the force between the cracks is effectively transferred. This shows that the local weak meshless method considering the interaction between crack surfaces is effective, and the established calculation method can effectively calculate the interaction between the two sides of the crack.


## 4.2. Oblique Crack Model

### 4.2.1. Elastic Crack Model

There is an oblique through crack in the middle of the model in Figure6, which divides the model into two parts. The parameters and boundary conditions of the model are consistent with those of test 1, but the crack is oblique. The normal stiffness is 20 times the elastic modulus, and the tangential stiffness changes. The correctness of the algorithm is verified by setting different tangential stiffness coefficients. Five test models are set, and the tangential stiffness is set to 0 times, 0.05 times, 0.1 times, 0.15 times, and 0.2 times the elastic modulus.Figure 6
Schematic diagram of oblique through crack.It should be noted that the crack is implicit in this algorithm. The results do not include crack information, only field node information. The field node at the crack of the displacement nephogram is not on the crack surface. In addition, the value shown in the legend is not the maximum and minimum value of the cloud map but the segmented value. The maximum values of the calculated results are compared in Table1.Table 1
Horizontal displacement of the model.
Tangential stiffness (Pa)0.0E0.05E0.1E0.15E0.2EMaximum horizontal displacement (mm)LRPIM7.069e − 26.829e − 26.604e − 26.393e − 26.195e − 2FDM6.997e − 26.887e − 26.782e − 26.682e − 26.585e − 2Error (%)1.00.82.64.35.9The horizontal displacement nephograms of the five test results are shown in Figure7. The maximum horizontal displacements are approximately 7.069e-2 mm, 6.829e-2 mm, 6.604e-2 mm, 6.393e-2 mm, and 6.195e-2 mm. It can be understood that, due to the increasing tangential stiffness, the deformation of the model on both sides of the crack is constrained, resulting in the decrease in the lateral deformation, which is consistent with the expected calculation results.Figure 7
Horizontal displacement nephogram (m).
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)As shown in Table1, the maximum horizontal displacement decreased with increasing tangential stiffness, and the calculation results obtained by the two algorithms were very close.The above derivation and calculation results assume that, under elastic conditions, the material and crack have no yield criterion or constitutive relation under elastic-plastic conditions. The model considering the plastic deformation of cracks will be calculated and analyzed in the following section.


### 4.2.2. Elastoplastic Crack Model

The interface elements in FLAC3D also adopt the Mohr–Coulomb strength criterion and ideal plastic constitutive model. The contact surface parameters include the normal stiffness, tangential stiffness, cohesion, internal friction angle, tensile strength, and expansion angle. The tensile strength and expansion angle of the contact surface are not considered in this paper.The crack is set as elastic-plastic, the cohesion is 0 kPa, and the friction coefficient is 0.01. A normal stiffness of 20 and tangential stiffness of 0.1 times the elastic modulus are set.The calculation results of the LRPIM and FLAC3D are shown in Figure8. The results of the two numerical examples have the same displacement nephogram. As shown in Table 2, under ideal plastic conditions, the maximum horizontal displacement of the model using the LRPIM is 7.051e − 2 mm. The results of the two methods are close, and the error is only 0.8%.Figure 8
Horizontal displacement nephogram (m). (a) FLAC3D. (b) LRPIM.
(a)(b)Table 2
Horizontal displacement comparison.
MethodLRPIMFLAC3DMaximum horizontal displacement (mm)7.051e − 26.995e − 2Error (%)0.8


## 4.2.1. Elastic Crack Model

There is an oblique through crack in the middle of the model in Figure6, which divides the model into two parts. The parameters and boundary conditions of the model are consistent with those of test 1, but the crack is oblique. The normal stiffness is 20 times the elastic modulus, and the tangential stiffness changes. The correctness of the algorithm is verified by setting different tangential stiffness coefficients. Five test models are set, and the tangential stiffness is set to 0 times, 0.05 times, 0.1 times, 0.15 times, and 0.2 times the elastic modulus.Figure 6
Schematic diagram of oblique through crack.It should be noted that the crack is implicit in this algorithm. The results do not include crack information, only field node information. The field node at the crack of the displacement nephogram is not on the crack surface. In addition, the value shown in the legend is not the maximum and minimum value of the cloud map but the segmented value. The maximum values of the calculated results are compared in Table1.Table 1
Horizontal displacement of the model.
Tangential stiffness (Pa)0.0E0.05E0.1E0.15E0.2EMaximum horizontal displacement (mm)LRPIM7.069e − 26.829e − 26.604e − 26.393e − 26.195e − 2FDM6.997e − 26.887e − 26.782e − 26.682e − 26.585e − 2Error (%)1.00.82.64.35.9The horizontal displacement nephograms of the five test results are shown in Figure7. The maximum horizontal displacements are approximately 7.069e-2 mm, 6.829e-2 mm, 6.604e-2 mm, 6.393e-2 mm, and 6.195e-2 mm. It can be understood that, due to the increasing tangential stiffness, the deformation of the model on both sides of the crack is constrained, resulting in the decrease in the lateral deformation, which is consistent with the expected calculation results.Figure 7
Horizontal displacement nephogram (m).
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)As shown in Table1, the maximum horizontal displacement decreased with increasing tangential stiffness, and the calculation results obtained by the two algorithms were very close.The above derivation and calculation results assume that, under elastic conditions, the material and crack have no yield criterion or constitutive relation under elastic-plastic conditions. The model considering the plastic deformation of cracks will be calculated and analyzed in the following section.


## 4.2.2. Elastoplastic Crack Model

The interface elements in FLAC3D also adopt the Mohr–Coulomb strength criterion and ideal plastic constitutive model. The contact surface parameters include the normal stiffness, tangential stiffness, cohesion, internal friction angle, tensile strength, and expansion angle. The tensile strength and expansion angle of the contact surface are not considered in this paper.The crack is set as elastic-plastic, the cohesion is 0 kPa, and the friction coefficient is 0.01. A normal stiffness of 20 and tangential stiffness of 0.1 times the elastic modulus are set.The calculation results of the LRPIM and FLAC3D are shown in Figure8. The results of the two numerical examples have the same displacement nephogram. As shown in Table 2, under ideal plastic conditions, the maximum horizontal displacement of the model using the LRPIM is 7.051e − 2 mm. The results of the two methods are close, and the error is only 0.8%.Figure 8
Horizontal displacement nephogram (m). (a) FLAC3D. (b) LRPIM.
(a)(b)Table 2
Horizontal displacement comparison.
MethodLRPIMFLAC3DMaximum horizontal displacement (mm)7.051e − 26.995e − 2Error (%)0.8


## 4.3. Built-In Inclined Crack Model

For shear cracks, the stress intensity factor isKI=0 because there is no relative displacement in the normal direction. According to the maximum circumferential stress theory, θ=70.5°.Under the action of axial pressure, the cracks of brittle rock with built-in inclined cracks will expand, and the crack propagation pattern is shown in Figure9 [26].Figure 9
Crack propagation under uniaxial compression [26].The model parameters are the same as the parameters of common marble. The model information is set as follows: dimensions of 62 mm × 25 mm × 110 mm, elastic modulus of 17 GPa, Poisson’s ratio of 0.25, internal friction angle of 30°, and fixed displacement of the upper and lower surfaces of the model. The angle between the crack and the horizontal plane is 45°, and the crack length is 20 mm. The cohesive force of the model material is 20.9 kPa, the cohesive force of the crack is 2 kPa, the friction coefficient is 0.01, the normal stiffness is twice the elastic modulus, and the tangential stiffness is 0.01 times the elastic modulus. The force diagram of the model is shown in Figure10.Figure 10
Schematic diagram of cracks subjected to uniaxial compression.A vertical displacement of 0.2 mm is applied on the upper surface of the model to simulate the stress state of compression and promote the expansion of cracks. Figure11 shows the vertical displacement nephogram and stress nephogram in the initial state. There is discontinuity in the model, while the vertical stress at the crack tip is obviously larger.Figure 11
Nephogram of the model subjected to unidirectional compression. (a) Vertical displacement (m). (b) Vertical stress (Pa).
(a)(b)Figure12 shows the crack morphology calculated by the LRPIM, and the two ends of the crack are curved and extend upward, which demonstrates the propagation process of the wing crack.Figure 12
Crack propagation under uniaxial compression.


## 5. Discussion

### 5.1. Simulation Analysis of Shear Cracks Using the LRPIM

Under certain assumptions, the governing equation of the shear cracks under three-dimensional conditions is derived and discretized into a programmable mode to realize the solution and expansion of shear crack propagation in the LRPIM.The constitutive relation of the shear cracks is established according to elastic-plastic theory. The elastic-plastic stress-strain relationship and yield criterion of the cracks are derived according to the Mohr–Coulomb strength criterion, which is commonly used in rock mechanics. Based on the above theory, the incremental iteration method is used to solve the problem.The results of several numerical examples are basically consistent with the expected results, which verifies the correctness of the relevant algorithms. The horizontal through crack model verifies that the crack realizes force transfer. A model of the oblique through cracks is established, and the deformation of the cracks in the elastic and elastic-plastic conditions is calculated. Compared with the FDM results, the source of the difference is analyzed. Compared with the calculation results of the classical model with built-in inclined cracks, the ability of this method to simulate the propagation of shear composite cracks is verified.The results of the numerical model calculation show that the developed method can effectively simulate the shear crack and its propagation in the rock mass and can accurately calculate the force between the two sides of the crack.


### 5.2. Comparative Analysis

In the simulation of crack propagation, compared with the results of classical physical experiments, the crack propagation is not sufficiently smooth, which indicates that the calculation accuracy needs to be improved. In meshless computing, there are many factors that affect the accuracy of the model, and the influence law is complex. Different distributions of model nodes may result in different calculation results, which may affect the crack propagation path of the model. When simulating the crack propagation of a rock mass, the crack representation methods used are all implicit. However, this representation is more flexible than other methods. Because there are no nodes on the crack, the stress-strain state of the crack cannot be accurately obtained. In the future, the display method can be used to characterize cracks.


## 5.1. Simulation Analysis of Shear Cracks Using the LRPIM

Under certain assumptions, the governing equation of the shear cracks under three-dimensional conditions is derived and discretized into a programmable mode to realize the solution and expansion of shear crack propagation in the LRPIM.The constitutive relation of the shear cracks is established according to elastic-plastic theory. The elastic-plastic stress-strain relationship and yield criterion of the cracks are derived according to the Mohr–Coulomb strength criterion, which is commonly used in rock mechanics. Based on the above theory, the incremental iteration method is used to solve the problem.The results of several numerical examples are basically consistent with the expected results, which verifies the correctness of the relevant algorithms. The horizontal through crack model verifies that the crack realizes force transfer. A model of the oblique through cracks is established, and the deformation of the cracks in the elastic and elastic-plastic conditions is calculated. Compared with the FDM results, the source of the difference is analyzed. Compared with the calculation results of the classical model with built-in inclined cracks, the ability of this method to simulate the propagation of shear composite cracks is verified.The results of the numerical model calculation show that the developed method can effectively simulate the shear crack and its propagation in the rock mass and can accurately calculate the force between the two sides of the crack.


## 5.2. Comparative Analysis

In the simulation of crack propagation, compared with the results of classical physical experiments, the crack propagation is not sufficiently smooth, which indicates that the calculation accuracy needs to be improved. In meshless computing, there are many factors that affect the accuracy of the model, and the influence law is complex. Different distributions of model nodes may result in different calculation results, which may affect the crack propagation path of the model. When simulating the crack propagation of a rock mass, the crack representation methods used are all implicit. However, this representation is more flexible than other methods. Because there are no nodes on the crack, the stress-strain state of the crack cannot be accurately obtained. In the future, the display method can be used to characterize cracks.


## 6. Conclusion

In this paper, the governing equation of the shear cracks under three-dimensional conditions is derived. The constitutive relation of the shear cracks is established according to elastic-plastic theory. The stress-strain relationship of the model material and crack is unified, and the solution is solved by the incremental iteration method, which realizes the solution and expansion of shear crack propagation in the LRPIM. The second development in the LRPIM program is carried out to establish the related computing ability. Several numerical examples are used to verify the reliability and accuracy of the algorithm. Compared with the calculation results, the source of the difference is analyzed.The research results show that the LRPIM, as a real meshless method, can simulate the shear crack of a rock mass and avoid the inconvenience of meshing, so it has great application potential. The experimental results show that further work is needed to improve the calculation accuracy.



---

*Source: 6654812-2021-03-25.xml*