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.
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. Section 2 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 aswhere , , and 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 and the corresponding boundary in the following formwhere is a weight function or test function centered at field node I.
The follow formula is obtained by integral operation:where is the component of surface force vector on the boundary.
We bring formula (10) into formula (6),where is the inner boundary of the integral domain; is the part of the natural boundary that intersects the integral domain; and is the part of the essential boundary that intersects the integral domain
In 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.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.where is the weight function matrix and is the derivative matrix of weight function.
Bring formulas (5) and (6) into formula (7):
The matrix form can be expressed aswhere is the stiffness matrix of the field node I.where is the nodal force vector, including the physical force in the problem domain and the surface force on the natural boundary.
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.
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 Figure 1. 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 axis, and the two tangent directions are along the axis and axis. Two points are set up in the positive and negative directions of the axis, and there is a small distance from the origin, which represent the material above and below the crack, respectively.
The energy equation considering crack interaction is shown in equation (13), which represents the work.where is the difference in the normal displacement of the materials on both sides of the crack. If 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. and represent the difference in the displacement between the two shear directions in the local coordinate system.where , , , , , and 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:where is the displacement vector of the field node in the support domain of the interpolation node, and the vector length is .and is the displacement shape function of node i in the following form:
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:
Here, is set as the node stress vector and is the material matrix; then, σn is the normal force, and and are the component forces in the two tangent directions.
kn and ks are the normal and tangential stiffnesses, respectively.
It is necessary to transform the displacement of the nodes from the global coordinate system to the local coordinate system. is a transformation matrix that converts global coordinates to local coordinates.
The displacement in the local coordinate system is obtained as follows:
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:
By variational calculation, 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 . 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:
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 Figure 2. A local coordinate system is established at a node on the crack, and the normal stress of the crack is . If the resultant force of the tangential force is expressed as , then the value of the shear stress in the direction perpendicular to and is zero. It is assumed here that the direction of is the same as that of the relative displacement of the crack.
According to the Mohr–Coulomb strength criterion, the general form of the yield function is as follows:where and are the cohesion and internal friction of the crack, is the shear stress of the crack, and is the normal direction of the crack.
The elastic matrix of the stress displacement relationship can be written as follows.
In 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 is m, that is, the strain is replaced by the displacement.
The plastic matrix is as follows:
Under the ideal elastic-plastic condition, A = 0. An adaptive flow criterion is adopted as . The elastic-plastic deformation matrix of the shear crack surface is as follows:
The elastic-plastic matrix after yielding is as follows:
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 is , and after yielding, the elastic-plastic deformation matrix is . 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:
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 Figure 3, 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.
(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 the nth-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 are , 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 Figure 4.
If the interaction between cracks is not considered, the calculated displacement nephogram is shown in Figure 5(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.
(a)
(b)
If the interaction of cracks is considered, the vertical displacement nephogram of the model is obtained, as shown in Figure 5(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 Figure 6, 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.
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 Table 1.
The horizontal displacement nephograms of the five test results are shown in Figure 7. 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.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
As shown in Table 1, 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 Figure 8. 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%.
(a)
(b)
4.3. Built-In Inclined Crack Model
For shear cracks, the stress intensity factor is because there is no relative displacement in the normal direction. According to the maximum circumferential stress theory, .
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 Figure 9 [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 Figure 10.
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. Figure 11 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.
(a)
(b)
Figure 12 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.
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.
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.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the Natural Key R&D Program of China (Grant no. 2018YFC1508703) and level I R&D project for Yellow River Engineering Consulting Co., Ltd. (Grant no. 2017-ky(02)).