Research on Shear Behavior and Crack Evolution of Symmetrical Discontinuous Rock Joints Based on FEM-CZM

The discontinuous joints are an essential type of natural joints. The normal force, joint persistency, and distribution exert great influences on the shear resistance of the rock joints. By simulating the uniaxial compression experiment and Brazilian test, the material parameters and the basic size standard for meshing were determined. The symmetrical discontinuous joint distribution of three types were established, the cohesive elements were inserted between the solid elements, and the numerical simulation of the shear test was conducted. The effects of joint distribution, joint continuity, and normal stress on the shear resistance of joint rock were investigated, and the law of crack evolution was analyzed. The results showed that the shear process of discontinuous joints can be divided into four stages: elastic stage, strengthening stage, plastic stage, and residual stress stage. For the scattered joint distribution, the rock bridge can provide more reinforcement for the joints, which enhances the shear resistance of the joints, the stress concentration point at the end of the joint is easy to accumulate more fracture energy, which induces the initiation of the cracks, and under the influence of unbalanced torque, the both-sided joint distribution is more likely to produce tension damage.


Introduction
In nature, when deformed to a certain degree, the rock mass is no longer continuous and intact, with cracks of different sizes initiated inside. If there is no significant displacement of the rock mass on both sides of the fracture surface, the crack surface is called the joint. A lot of engineering projects have shown that the joint features the weak part of the natural rock mass, and the natural rock mass failure starts from the joints [1][2][3][4][5]. The evolution, distribution, lengths, and mechanical properties of joints affect the stability of rock mass. From the perspective of engineering application, the jointed rock mass can be divided into the rock mass with continuous joints and the rock mass with discontinuous joints. Discontinuous joints are commonly seen in nature. The discontinuous joint rock mass contains rock bridge and joints, whose responses under mechanical conditions constitute the overall mechanical similarity with laboratory tests in results. Chang et al. [45] established a two-dimensional Brazilian test model and embedded a zero-thickness cohesive element to investigate the complex crack behaviors in a layered disc with interfacial cracks. Through the comparison between their results with those obtained in laboratory tests, it was found that they are quite consistent in crack development mode and mechanical properties. Besides, they analyzed the mesh sensitivity of cohesive elements. Haddad et al. [46] established a two-dimensional hydraulic fracturing model and analyzed the effects of the initial crack length and mesh refinement on calculation convergence. Wang et al. [47] used the CZM to model the shearing process of bolted jointed rock. The numerical model is verified by the comparison with the experiment results, which demonstrated that it can be effectively applied to the simulation of joint shearing process. Han et al. [48] studied the shear behavior of rock-like materials with cracks by using the FEM-CZM method. The mechanical properties, shear deformation, and cracking behaviors of specimen were respectively discussed.
At present, considering the discontinuity of the jointed rock mass, most scholars prefer to use discrete element software to simulate the shear of discontinuous joints, and it is relatively rare to use the FEM-CZM method. This paper adopted the finite element software ABAQUS to establish shear models with different joint persistence, densities, and spatial distributions, and inserted cohesive elements thoroughly to generate potential crack surfaces inside the model. By simulating the uniaxial compression test and the Brazilian test, and comparing the results with those by laboratory experiments, we determined the material parameters and the basis for mesh division. In addition, we investigated the effect of joint distribution, joint persistence, and normal stress on the shear resistance of rock, and analyzed the crack evolution law during the shear process of discontinuous joints of three different distribution forms.

Initial Linear Elastic Traction-Displacement
Several constitutive relationships of the CZM have been developed on the basis of an effective displacement and an effective traction, and different constitutive models have their own applicable conditions [49]. Cong [50] discovered that the energy-based Separation-Traction linear model (S-T model) delivers better performance in simulating fracture of brittle materials with good applicability. Alfano [51] concluded that the S-T linear model features better simulation accuracy and calculation efficiency through comparative analysis. Therefore, this paper selected the S-T linear model as the criterion for the damage of cohesive elements, whose constitutive relationship is shown in Figure 1.
Symmetry 2020, 12, x FOR PEER REVIEW 3 of 21 established a two-dimensional Brazilian test model and embedded a zero-thickness cohesive element to investigate the complex crack behaviors in a layered disc with interfacial cracks. Through the comparison between their results with those obtained in laboratory tests, it was found that they are quite consistent in crack development mode and mechanical properties. Besides, they analyzed the mesh sensitivity of cohesive elements. Haddad et al. [46] established a two-dimensional hydraulic fracturing model and analyzed the effects of the initial crack length and mesh refinement on calculation convergence. Wang et al. [47] used the CZM to model the shearing process of bolted jointed rock. The numerical model is verified by the comparison with the experiment results, which demonstrated that it can be effectively applied to the simulation of joint shearing process. Han et al. [48] studied the shear behavior of rock-like materials with cracks by using the FEM-CZM method. The mechanical properties, shear deformation, and cracking behaviors of specimen were respectively discussed. At present, considering the discontinuity of the jointed rock mass, most scholars prefer to use discrete element software to simulate the shear of discontinuous joints, and it is relatively rare to use the FEM-CZM method. This paper adopted the finite element software ABAQUS to establish shear models with different joint persistence, densities, and spatial distributions, and inserted cohesive elements thoroughly to generate potential crack surfaces inside the model. By simulating the uniaxial compression test and the Brazilian test, and comparing the results with those by laboratory experiments, we determined the material parameters and the basis for mesh division. In addition, we investigated the effect of joint distribution, joint persistence, and normal stress on the shear resistance of rock, and analyzed the crack evolution law during the shear process of discontinuous joints of three different distribution forms.

Initial Linear Elastic Traction-Displacement
Several constitutive relationships of the CZM have been developed on the basis of an effective displacement and an effective traction, and different constitutive models have their own applicable conditions [49]. Cong [50] discovered that the energy-based Separation-Traction linear model (S-T model) delivers better performance in simulating fracture of brittle materials with good applicability. Alfano [51] concluded that the S-T linear model features better simulation accuracy and calculation efficiency through comparative analysis. Therefore, this paper selected the S-T linear model as the criterion for the damage of cohesive elements, whose constitutive relationship is shown in Figure 1.  The S-T model includes two stages: (1) the linear elastic stage, in which no damage occurs to the cohesive element and the element stiffness remains constant, and (2) the linear damage stage, in which the cohesive element damage evolves, and the element stiffness decreases continuously. In the first stage, we can use a matrix containing nominal strain and nominal stress for description. In order to facilitate parameter calibration, the constitutive thickness of the cohesive element was set to 1, and the geometric thickness was set to 0 by default, which can make the node separation displacement equal to the nominal strain. The elastic stage of the S-T criterion can be expressed as: where t is the nominal traction stress vector consisting of three components, t n , t s , and t t , ε is the nominal strain consisting of ε n , ε s , and ε t , T 0 is initial thickness of the cohesive element, δ is the displacement of the cohesive element, and K is the corresponding stiffness matrix.

Linear Damage Stage of Cohesive Element
It was assumed that there are three traction forces on the crack surface, namely the normal traction force t n and two tangential traction forces t s and t t , all of which gradually decrease as the crack propagates. When the traction force becomes 0, the element will be wholly damaged. This study utilized type I and type II fracture to study the rock properties. The fracture energy can be represented by the triangular area in Figure 1a, which is type I fracture energy and type II fracture energy, respectively.
Describing the damage evolution process of cohesive elements entails at least three of the four parameters, including the initial stiffness, K 0 , the ultimate strength, t 0 , the critical fracture energy, G f , and the final failure displacement, δ f . In this paper, we selected K 0 , t 0 , and G f as the research parameters. The quadratic nominal stress criterion (QUADS) was used for judging the start of damage evolution, and, according to it, the damage starts to evolve when the quadratic interaction function involving the nominal stress ratios (as defined in the expression below) is equal to 1.
where the < > is a Macaulay bracket, which ensures that the damage evolution only occurs under the action of forwarding tension and the compressive stress will not cause damage to the cohesive element. The QUADS criterion ensures that the damage of the cohesive element be generated under the combined action of normal stress and tangential stress.
To describe the degree of damage evolution, we introduced a damage variable, D, the value of which gradually changes from 0 to 1 as the load increases. When D reaches 1, the cohesive element also reaches the ultimate bearing capacity with the cohesive element automatically deleted, and the free edges are formed in the continuously arranged elements, which means the cracks are initiated. The calculation of the damage variable, D, is shown as follows: where δ m f is the corresponding displacement at the moment of complete failure, δ m0 is related to the displacement at damage initiation, and δ m denotes the maximum displacement attained during the loading, which is expressed below: Therefore, the changes of the parameters of the T-S model with the cumulative damage are as follows: Traction:

Cohesive Element Inserting Process
In order to simulate the random development progress of cracks, the cohesive element was inserted into the solid elements to form a potential crack surface. In the simulation model, the solid discrete elements form a continuum through zero-thickness cohesive elements. The top and bottom surfaces of each zero-thickness cohesive element are connected to the solid element, with the nodes shared. The process of embedding cohesive elements is shown in Figure 2. The node renumbering follows the right-hand rule to ensure that each node has a unique number. Since the model has a large number of solid elements, it is infeasible to manually insert cohesive elements. Therefore, a calculation program was developed based on the Python language, which can automatically insert cohesive elements into the solid elements.

Cohesive Element Inserting Process
In order to simulate the random development progress of cracks, the cohesive element was inserted into the solid elements to form a potential crack surface. In the simulation model, the solid discrete elements form a continuum through zero-thickness cohesive elements. The top and bottom surfaces of each zero-thickness cohesive element are connected to the solid element, with the nodes shared. The process of embedding cohesive elements is shown in Figure 2. The node renumbering follows the right-hand rule to ensure that each node has a unique number. Since the model has a large number of solid elements, it is infeasible to manually insert cohesive elements. Therefore, a calculation program was developed based on the Python language, which can automatically insert cohesive elements into the solid elements.

Parameters Determination
To obtain the main mechanical parameters of sandstone, according to the recommendations of the International Society for Rock Mechanics (ISRM), we prepared a cylindrical specimen with the diameter of 50 mm and the height of 100 mm for the uniaxial compression test, as shown in Figure

Parameters Determination
To obtain the main mechanical parameters of sandstone, according to the recommendations of the International Society for Rock Mechanics (ISRM), we prepared a cylindrical specimen with the diameter of 50 mm and the height of 100 mm for the uniaxial compression test, as shown in Figure 3a. The loading rate was set to 0.01 mm/min and we stopped the loading after the test piece was obviously damaged and the loading curve exceeded the peak [52]. During the test, we kept other variables unchanged and did not apply the confining pressure. Then, we established a three-dimensional uniaxial compression numerical model with the same dimensions as the laboratory experiment, which was inserted into the cohesive element by the finite element plug-in. Besides, we set a loading plate, which was given the stiffness much higher than that of the specimen, to ensure that the load can be evenly applied on the top of the test piece, as shown in Figure 3b. When the lower loading plate was fixed, the upper one moved downward at a speed of 0.1 mm/s until the specimen was destroyed. Under ideal conditions, the numerical simulation should set the same loading speed as the laboratory experiment, which, however, requires a smaller time step and longer calculation time. Thus, we amplified the loading speed by ten times. Through the parameter comparison, the parameters of cohesive elements were determined. Multiple sets of parameters were used for simulation, whose results were compared with those obtained by laboratory experiment. If the parameters are similar, the parameters are considered to be finally determined.
Symmetry 2020, 12, x FOR PEER REVIEW 6 of 21 variables unchanged and did not apply the confining pressure. Then, we established a threedimensional uniaxial compression numerical model with the same dimensions as the laboratory experiment, which was inserted into the cohesive element by the finite element plug-in. Besides, we set a loading plate, which was given the stiffness much higher than that of the specimen, to ensure that the load can be evenly applied on the top of the test piece, as shown in Figure 3b. When the lower loading plate was fixed, the upper one moved downward at a speed of 0.1 mm/s until the specimen was destroyed. Under ideal conditions, the numerical simulation should set the same loading speed as the laboratory experiment, which, however, requires a smaller time step and longer calculation time. Thus, we amplified the loading speed by ten times. Through the parameter comparison, the parameters of cohesive elements were determined. Multiple sets of parameters were used for simulation, whose results were compared with those obtained by laboratory experiment. If the parameters are similar, the parameters are considered to be finally determined.
(a) (b) After many attempts, the uniaxial compression numerical model obtained the mechanical properties and crack evolution characteristics of the numerical simulation similar to the laboratory experiment, as shown in Figure 4. It can be found from the results that the peak compressive strength of the laboratory experiment is 43.4 MPa with the corresponding axial strain being 0.41%, while the peak compressive strength of the numerical simulation is 45.6 MPa with the corresponding axial strain being 0.43%, both of which are very close. Both of the stress-strain curves obtained from the numerical simulation and the laboratory experiment show a trend of sudden stress drop after reaching the peak. After many attempts, the uniaxial compression numerical model obtained the mechanical properties and crack evolution characteristics of the numerical simulation similar to the laboratory experiment, as shown in Figure 4. It can be found from the results that the peak compressive strength of the laboratory experiment is 43.4 MPa with the corresponding axial strain being 0.41%, while the peak compressive strength of the numerical simulation is 45.6 MPa with the corresponding axial strain being 0.43%, both of which are very close. Both of the stress-strain curves obtained from the numerical simulation and the laboratory experiment show a trend of sudden stress drop after reaching the peak.
At the same time, through a summary of the crack evolution process in the numerical simulation, we found that the small cracks first appeared in the middle part of the specimen. As the axial load increased, the small cracks further propagated, gradually connected with each other, and finally formed the main fissure penetrated through the specimen. This phenomenon is the same with the uniaxial compression crack growth process observed by Basu et al. [53]. Simultaneously, we monitored SDEG, defined as the scalar stiffness degradation variable of the cohesive element during the uniaxial compression process, which can reflect the damage degree of the cohesive element. It can be found that the damage of the cohesive element derived from the inner part of the specimen, and under the action of tensile stress, extended along the diagonal of the specimen and eventually penetrated the entire specimen. This process corresponds to the crack propagation process. The comparison between the results of the numerical simulation and those of the laboratory experiment reveals the reasonableness of the selection of the parameters. The final parameters of the numerical simulation element are shown in Table 1. At the same time, through a summary of the crack evolution process in the numerical simulation, we found that the small cracks first appeared in the middle part of the specimen. As the axial load increased, the small cracks further propagated, gradually connected with each other, and finally formed the main fissure penetrated through the specimen. This phenomenon is the same with the uniaxial compression crack growth process observed by Basu et al. [53]. Simultaneously, we monitored SDEG, defined as the scalar stiffness degradation variable of the cohesive element during the uniaxial compression process, which can reflect the damage degree of the cohesive element. It can be found that the damage of the cohesive element derived from the inner part of the specimen, and under the action of tensile stress, extended along the diagonal of the specimen and eventually penetrated the entire specimen. This process corresponds to the crack propagation process. The comparison between the results of the numerical simulation and those of the laboratory experiment reveals the reasonableness of the selection of the parameters. The final parameters of the numerical simulation element are shown in Table 1. Since the crack evolution characteristics and mechanical characteristics of the specimen are  Since the crack evolution characteristics and mechanical characteristics of the specimen are closely related to the size of the cohesive element, it is necessary to perform the mesh sensitivity analysis of the specimen. As shown in Figure 5a, a two-dimensional Brazilian test numerical model was established with the disc's diameter of 50 mm and the plane's outer thickness of 40 mm. The sizes of the mesh sizes were set respectively to 0.5, 1, 2, and 3 mm to carry out the Brazilian test simulation, through which we obtained the load-displacement curve in the process of simulation, as shown in Figure 6. With the mesh of 1 mm, when the normal displacement is 0.39 mm, the maximum load along the normal direction is 9.85 kN, while the peak value of the load obtained in the laboratory experiment is 10.0 kN with the corresponding normal displacement of 0.41 mm. Besides, we obtained the final fracture form of the numerical simulation specimen, as shown in Figure 5b, which is similar to the Symmetry 2020, 12, 1314 8 of 20 laboratory experiment result and consistent with the research of Chang et al. [44]. It was found that after the load reached the peak, a vertically penetrated fissure was generated in the middle of the specimen. Therefore, we set the mesh size to 1 mm, with the consideration of the calculation efficiency. simulation, through which we obtained the load-displacement curve in the process of simulation, as shown in Figure 6. With the mesh of 1 mm, when the normal displacement is 0.39 mm, the maximum load along the normal direction is 9.85 kN, while the peak value of the load obtained in the laboratory experiment is 10.0 kN with the corresponding normal displacement of 0.41 mm. Besides, we obtained the final fracture form of the numerical simulation specimen, as shown in Figure 5b, which is similar to the laboratory experiment result and consistent with the research of Chang et al. [44]. It was found that after the load reached the peak, a vertically penetrated fissure was generated in the middle of the specimen. Therefore, we set the mesh size to 1 mm, with the consideration of the calculation efficiency.

Model Establishment
Since the shear-induced cracks can cause collision and misalignment between the elements, we adopted global contact to define the contact behaviors. The global contact algorithm, however, can simulation, through which we obtained the load-displacement curve in the process of simulation, as shown in Figure 6. With the mesh of 1 mm, when the normal displacement is 0.39 mm, the maximum load along the normal direction is 9.85 kN, while the peak value of the load obtained in the laboratory experiment is 10.0 kN with the corresponding normal displacement of 0.41 mm. Besides, we obtained the final fracture form of the numerical simulation specimen, as shown in Figure 5b, which is similar to the laboratory experiment result and consistent with the research of Chang et al. [44]. It was found that after the load reached the peak, a vertically penetrated fissure was generated in the middle of the specimen. Therefore, we set the mesh size to 1 mm, with the consideration of the calculation efficiency.

Model Establishment
Since the shear-induced cracks can cause collision and misalignment between the elements, we adopted global contact to define the contact behaviors. The global contact algorithm, however, can

Model Establishment
Since the shear-induced cracks can cause collision and misalignment between the elements, we adopted global contact to define the contact behaviors. The global contact algorithm, however, can only be applied to the three-dimensional surface contact, thus we utilized a three-dimensional model for analysis. A model with the dimensions of 200 × 100 × 1 mm was established, and the zero-thickness cohesive elements were inserted into the specimens, as shown in Figure 7. During the numerical simulation, we simulated the loading process of the shear box through setting appropriate boundary conditions. During the shearing process, we applied a constant normal load, σ, to the top and an appropriate shear speed, τ, to the upper part, with the lower part remaining static. only be applied to the three-dimensional surface contact, thus we utilized a three-dimensional model for analysis. A model with the dimensions of 200 × 100 × 1 mm was established, and the zero-thickness cohesive elements were inserted into the specimens, as shown in Figure 7. During the numerical simulation, we simulated the loading process of the shear box through setting appropriate boundary conditions. During the shearing process, we applied a constant normal load, σ, to the top and an appropriate shear speed, τ, to the upper part, with the lower part remaining static. To accurately describe the characteristics of the discontinuous joints, we introduced the joint persistence, , as shown below: where is the length of the rock bridge, and is the length of the joint. In this paper, 0 < < 1, which means that all the specimens are composed of joints and rock bridges.
To analyze the effect of the discontinuous joint distribution, joint persistence, and normal stress on the shear resistance of rock, we established joint distribution models of three types, namely, both sided (type-I), scattered (type-II), and central (type-III), as shown in Figure 8. Among them, only the type-I joints are open, and there are two discontinuous joints in the type-II and type-III, both of which share the same length (joint persistence) and are symmetrical about the central axis. We set three normal stresses of 0.5, 1.0, and 1.5 MPa respectively, and set four different joint persistences of 0.2, 0.4, 0.6, and 0.8 respectively, with the corresponding rock bridge lengths of 160, 120, 80, and 40 mm. In this paper, a total of 36 discontinuous joint shear models under different conditions were established. It should be stated that this paper did not take consideration of the effect of shear rate on shear performance. The displacement method was used to control the application of the shear load. The shear rate was set to 0.06 mm/min to ensure that the shear be performed under quasi-static conditions. Since the model owns a large number of calculation elements, in order to improve the calculation efficiency, we set a dynamic explicit step and employed the mass scaling technology, which can change the discrete mass matrix by setting mass scaling factor, leading to a corresponding improvement in the calculation efficiency. Mass scaling mainly affects the inertial force of the element; especially, the generated virtual inertial force can easily affect the calculation accuracy and convergence. For quasi-static simulations, it is usually necessary to ensure that the ratio of kinetic energy to potential energy do not exceed 0.1. The load application method adopted in this paper can satisfy this condition. To accurately describe the characteristics of the discontinuous joints, we introduced the joint persistence, p j , as shown below: where L r is the length of the rock bridge, and L j is the length of the joint. In this paper, 0 < p j < 1, which means that all the specimens are composed of joints and rock bridges.
To analyze the effect of the discontinuous joint distribution, joint persistence, and normal stress on the shear resistance of rock, we established joint distribution models of three types, namely, both sided (type-I), scattered (type-II), and central (type-III), as shown in Figure 8. Among them, only the type-I joints are open, and there are two discontinuous joints in the type-II and type-III, both of which share the same length (joint persistence) and are symmetrical about the central axis. We set three normal stresses of 0.5, 1.0, and 1.5 MPa respectively, and set four different joint persistences of 0.2, 0.4, 0.6, and 0.8 respectively, with the corresponding rock bridge lengths of 160, 120, 80, and 40 mm. In this paper, a total of 36 discontinuous joint shear models under different conditions were established. It should be stated that this paper did not take consideration of the effect of shear rate on shear performance. The displacement method was used to control the application of the shear load. The shear rate was set to 0.06 mm/min to ensure that the shear be performed under quasi-static conditions. Since the model owns a large number of calculation elements, in order to improve the calculation efficiency, we set a dynamic explicit step and employed the mass scaling technology, which can change the discrete mass matrix by setting mass scaling factor, leading to a corresponding improvement in the calculation efficiency. Mass scaling mainly affects the inertial force of the element; especially, the generated virtual inertial force can easily affect the calculation accuracy and convergence. For quasi-static simulations, it is usually necessary to ensure that the ratio of kinetic energy to potential energy do not exceed 0.1. The load application method adopted in this paper can satisfy this condition.

Influence of Joint Distribution on Shear Resistance
The stress-strain curves of the specimens with the persistence of 0.2 are shown in Figure 9. It can be divided into the elastic stage, strengthening stage, plastic stage, and residual stress stage (I-IV). By comparing the curves with different normal stresses, it can be found that the elastic modulus does not change with the variation of the normal stress (I), as proved by the fact that the change of normal stress has little effect on the elastic stage. During this stage, the cohesive elements remained intact with no crack appearing on the specimen. As the shear load was applied, the shear curve enters the strengthening stage (II), and the curve slope increases with the increase of normal stress. During this stage, part of the cohesive elements began to enter the linear damage stage (0 < D < 1), with the damaged cohesive elements (D = 1) deleted and initial cracks appearing. It is worth noting that the local failure of the specimen induced the sudden reduction of the shear resistance, which is shown as a sudden variation of the stress-strain curve in Figure 9a and c. As the shear load increases further, the curve enters the plastic stage (III). During this stage, the cracks further propagated, the shear resistance of the specimen began to decrease, and the peak shear stress showed a positive correlation with the normal stress. As the crack propagated further to form a continuous fissure, the curve enters the residual stage (IV). During this stage, the shear strength is provided by the mechanical bite force and the friction of the test piece.

Influence of Joint Distribution on Shear Resistance
The stress-strain curves of the specimens with the persistence of 0.2 are shown in Figure 9. It can be divided into the elastic stage, strengthening stage, plastic stage, and residual stress stage (I-IV). By comparing the curves with different normal stresses, it can be found that the elastic modulus does not change with the variation of the normal stress (I), as proved by the fact that the change of normal stress has little effect on the elastic stage. During this stage, the cohesive elements remained intact with no crack appearing on the specimen. As the shear load was applied, the shear curve enters the strengthening stage (II), and the curve slope increases with the increase of normal stress. During this stage, part of the cohesive elements began to enter the linear damage stage (0 < D < 1), with the damaged cohesive elements (D = 1) deleted and initial cracks appearing. It is worth noting that the local failure of the specimen induced the sudden reduction of the shear resistance, which is shown as a sudden variation of the stress-strain curve in Figure 9a,c. As the shear load increases further, the curve enters the plastic stage (III). During this stage, the cracks further propagated, the shear resistance of the specimen began to decrease, and the peak shear stress showed a positive correlation with the normal stress. As the crack propagated further to form a continuous fissure, the curve enters the residual stage (IV). During this stage, the shear strength is provided by the mechanical bite force and the friction of the test piece.
Comparing the peak strength of joints with different distributions, as shown in Figure 10a, it was found that the shear strengths of the type-II and type-III are the highest, with that of the type-I being the lowest. Stress concentration tends to occur at the ends of joints, from which the cracks are generated and propagated. Compared with the type-II and type-III, the type-I is more likely to produce tensile fractures at the joint ends, resulting in tensile cracks. As the cracks gradually extend through the rock bridge, the specimen enters the residual strength stage. Though the cracks of the type-II and type-III also start from the joint end, they can result in shear failure. Since the rock is more resistant to shear than the tension, the type-I cracks grow with the highest speed, also the specimen enters the residual strength stage. The difference between the strength peak and the residual strength was analyzed, as shown in Figure 10b. It was found that the type-II intensity decreased the most, followed by the type-III, and the type-I decreased the least. In the stress-strain curve, the type-II residual strength stage has a higher slope. It can be seen that type-II features the highest brittleness, while type-I features the highest plastic properties. This is due to more dispersed distribution of the type-II joints and more joint reinforcement provided by the rock bridge, thereby restraining the deformation of the specimen under the shear load. When the shear stress reaches the peak, the penetrated crack is formed in the rock bridge and the accumulated fracture energy is quickly released, making the shear strength decrease more rapidly. Comparing the peak strength of joints with different distributions, as shown in Figure 10a, it was found that the shear strengths of the type-II and type-III are the highest, with that of the type-I being the lowest. Stress concentration tends to occur at the ends of joints, from which the cracks are generated and propagated. Compared with the type-II and type-III, the type-I is more likely to produce tensile fractures at the joint ends, resulting in tensile cracks. As the cracks gradually extend through the rock bridge, the specimen enters the residual strength stage. Though the cracks of the type-II and type-III also start from the joint end, they can result in shear failure. Since the rock is more resistant to shear than the tension, the type-I cracks grow with the highest speed, also the specimen enters the residual strength stage. The difference between the strength peak and the residual strength was analyzed, as shown in Figure 10b. It was found that the type-Ⅱ intensity decreased the most, followed by the type-III, and the type-I decreased the least. In the stress-strain curve, the type-II residual strength stage has a higher slope. It can be seen that type-II features the highest brittleness, while type-I features the highest plastic properties. This is due to more dispersed distribution of the type-II joints and more joint reinforcement provided by the rock bridge, thereby restraining the deformation of the specimen under the shear load. When the shear stress reaches the peak, the penetrated crack is formed in the rock bridge and the accumulated fracture energy is quickly released, making the shear strength decrease more rapidly. Specimens with the persistence of 0.6 and the normal stress of 1.0 MPa were selected to obtain the vertical displacement field, U2. Figure 11a is the displacement field with no shear load applied, and Figure 11b is the displacement field when the specimen finally fails. The largest vertical displacement always occurs above the joint, while the displacement at the rock bridge is the smallest. During the final failure stage, the vertical displacement of type-I is the largest, while that of type-II is the smallest, which can reflect the degree of the specimen affected by the dilatancy effect. For type-I, since the joints are open, the upper part of the specimen is similar to a cantilever beam, and under Specimens with the persistence of 0.6 and the normal stress of 1.0 MPa were selected to obtain the vertical displacement field, U2. Figure 11a is the displacement field with no shear load applied, and Figure 11b is the displacement field when the specimen finally fails. The largest vertical displacement always occurs above the joint, while the displacement at the rock bridge is the smallest. During the final failure stage, the vertical displacement of type-I is the largest, while that of type-II is the smallest, which can reflect the degree of the specimen affected by the dilatancy effect. For type-I, since the joints are open, the upper part of the specimen is similar to a cantilever beam, and under normal stress, a negative vertical displacement was generated. During the final failure stage, under the action of the shear stress, a clockwise moment around the centroid of the specimen was generated. Under the effect of the moment, the two ends of the specimen produced vertical displacements in opposite directions. Since the rock bridge can provide a certain degree of reinforcement, the vertical displacement on the side away from the loading site is smaller, and the overall displacement field of type-I is not center-symmetric. For type-II and type-III, during the initial stage, the stress in type-II is more dispersed, making the vertical displacement smaller. During the final failure stage, the whole upper part of the specimen has a positive vertical displacement. In type-II, the three-segmented rock bridge provides more reinforcement for the joint, and the specimen has a smaller deformation. In summary, the type-II has the smallest dilatancy effect during the shearing process, and the specimen has higher rigidity and stability. In comparison, type-I is more likely to produce unbalanced moments, which leads to the instability of the specimen.

Effect of Joint Persistence on Shear Resistance
Taking the specimen with the normal stress of 1.0 MPa as an example, we investigated the effect of different joint persistence on the joint shear performance. The stress-strain curves of different types of joints distributed under 1 MPa normal stress are shown in Figure 12a-c. With the increase of the joint persistency, the shear stress and shear modulus are gradually decreasing. This is attributed to the fact that the rock bridge provides less reinforcement during the shearing process, with the stiffness of the specimen becoming lower. At the same time, the ratio of peak stress to residual stress gradually decreases, the curves become smoother, and the plasticity of the specimen is higher. The lower the joint persistence is, the steeper the stress-strain curve is, and the more brittle the specimen will be, which shows that the numerical simulation is closer to the direct shear test. The lower joint persistence can result in larger contact area of the specimen after the formation of the penetrated joint, which can provide more friction and mechanical bite force, thus the residual strength of the test piece is higher. The peak shear strengths of these twelve specimens were summarized, as shown in Figure 12d. It can be found that for the specimens with the three joint distribution forms, the peak shear strength and the joint continuity are negatively correlated. The type-II features the highest shear strength with different joint persistence, followed by type-III, and the type-I is the lowest. This shows that with the same overall length of the rock bridge, the scattered rock bridge can play a more significant role in strengthening the joint. In engineering practices, this finding may provide new ideas for strengthening joints by grouting.

Effect of Joint Persistence on Shear Resistance
Taking the specimen with the normal stress of 1.0 MPa as an example, we investigated the effect of different joint persistence on the joint shear performance. The stress-strain curves of different types of joints distributed under 1 MPa normal stress are shown in Figure 12a-c. With the increase of the joint persistency, the shear stress and shear modulus are gradually decreasing. This is attributed to the fact that the rock bridge provides less reinforcement during the shearing process, with the stiffness of the specimen becoming lower. At the same time, the ratio of peak stress to residual stress gradually decreases, the curves become smoother, and the plasticity of the specimen is higher. The lower the joint persistence is, the steeper the stress-strain curve is, and the more brittle the specimen will be, which shows that the numerical simulation is closer to the direct shear test. The lower joint persistence can result in larger contact area of the specimen after the formation of the penetrated joint, which can provide more friction and mechanical bite force, thus the residual strength of the test piece is higher. The peak shear strengths of these twelve specimens were summarized, as shown in Figure  12d. It can be found that for the specimens with the three joint distribution forms, the peak shear strength and the joint continuity are negatively correlated. The type-II features the highest shear strength with different joint persistence, followed by type-III, and the type-I is the lowest. This shows that with the same overall length of the rock bridge, the scattered rock bridge can play a more significant role in strengthening the joint. In engineering practices, this finding may provide new ideas for strengthening joints by grouting. We summarized the relationship of the peak shear strength of all specimens with the joint continuity and normal stress, as shown in Figure 13. It can be found that the shear strength of the specimen rises with the increase of the normal stress and drops with the increase of joint continuity. Compared with type-II and type-III, the type-I always has the lowest shear strength. This is because the rock bridge in the type-I cannot provide effective reinforcement for the joints, resulting in the specimen being more prone to fail from the joint tip, especially in the form of tensile failure. We summarized the relationship of the peak shear strength of all specimens with the joint continuity and normal stress, as shown in Figure 13. It can be found that the shear strength of the specimen rises with the increase of the normal stress and drops with the increase of joint continuity. Compared with type-II and type-III, the type-I always has the lowest shear strength. This is because the rock bridge in the type-I cannot provide effective reinforcement for the joints, resulting in the specimen being more prone to fail from the joint tip, especially in the form of tensile failure.

Crack Evolution Analysis
Under the action of shear load, the cohesive elements fail, and cracks appear on the specimen. The process of crack evolution can reflect the internal mechanical characteristics of the specimen. By modifying the model file, the ratio parameter of the mixed mode of initial damage MMIXDMI is defined. Its initial value was set to -1.0, and it can be expressed as follows: where is the sum of three types of fracture energy, and , , and are the ratios of type I, II, and III fracture energy to total fracture energy, respectively.
Cohesive elements can effectively simulate the three forms of failure, including open cracks, slip cracks, and tear cracks. The last two forms of failure can be collectively referred to as a shear failure. Therefore, we can analyze the driving force for cracking by calculating . With the application of shear load, the cohesive elements begin to bear shear stress, tensile stress, or mixed stresses. At this time, changes from 0 to 0.5, the fracture energy accumulated by shearing accounts for most of the total fracture energy of cohesive elements, corresponding to the MMIXDMI value of 1~0.5, and the cohesive elements fail mainly in the form of shear failure. When changes from 0.5 to 1, the fracture energy accumulated by traction accounts for most of the total fracture energy of cohesive

Crack Evolution Analysis
Under the action of shear load, the cohesive elements fail, and cracks appear on the specimen. The process of crack evolution can reflect the internal mechanical characteristics of the specimen. By modifying the model file, the ratio parameter of the mixed mode of initial damage MMIXDMI is defined. Its initial value was set to -1.0, and it can be expressed as follows: where G T is the sum of three types of fracture energy, and m 1 , m 2 , and m 3 are the ratios of type I, II, and III fracture energy to total fracture energy, respectively. Cohesive elements can effectively simulate the three forms of failure, including open cracks, slip cracks, and tear cracks. The last two forms of failure can be collectively referred to as a shear failure. Therefore, we can analyze the driving force for cracking by calculating m 1 . With the application of shear load, the cohesive elements begin to bear shear stress, tensile stress, or mixed stresses. At this time, m 1 changes from 0 to 0.5, the fracture energy accumulated by shearing accounts for most of the total fracture energy of cohesive elements, corresponding to the MMIXDMI value of 1~0.5, and the cohesive elements fail mainly in the form of shear failure. When m 1 changes from 0.5 to 1, the fracture energy accumulated by traction accounts for most of the total fracture energy of cohesive elements, corresponding to the MMIXDMI value of 0.5~0, and the cohesive elements fail mainly in the form of tensile failure. In particular, when the value of m 1 is 1, 0.5, and 0, the corresponding failure modes of cohesive elements are pure shear failure, tensile-shear failure, and pure tensile failure.
When the joint persistence is 0.6 and the normal load is 1.0 MPa, the stress field S11 of the specimen along the shear direction and the MMIXDMI of the cohesive elements are shown in Figure 14. The crack evolution process of the specimen can be divided into three stages, namely the crack initiation stage (a), crack evolution stage (b), and final failure stage (c).
Symmetry 2020, 12, x FOR PEER REVIEW 16 of 20 elements, corresponding to the MMIXDMI value of 0.5~0, and the cohesive elements fail mainly in the form of tensile failure. In particular, when the value of is 1, 0.5, and 0, the corresponding failure modes of cohesive elements are pure shear failure, tensile-shear failure, and pure tensile failure.
When the joint persistence is 0.6 and the normal load is 1.0 MPa, the stress field S11 of the specimen along the shear direction and the MMIXDMI of the cohesive elements are shown in Figure  14. The crack evolution process of the specimen can be divided into three stages, namely the crack initiation stage (a), crack evolution stage (b), and final failure stage (c). Firstly, when the shear load is just applied, the cohesive elements begin to accumulate fracture energy. At this time, the cohesive elements are in the linear elastic stage, and the traction force does not reach the initial damage stress, . No cracks appear on the specimen, which corresponds to the elastic stage in the stress-strain curve.
As the shear load continues to be applied, the cohesive elements accumulate more fracture energy. Due to the presence of joints inside the specimen, the stress concentration tends to occur at the end of the joint. As shown in Figure 15a, the stress gradually decreases from the joint end to other directions. The cohesive elements near the stress concentration point tend to accumulate more fracture energy, and thus reach the stage of crack evolution, as shown in Figure 15b. The stiffness degradation rate SDEG of the cohesive element at the joint end reaches 1.0 as the first, which is deleted to generate the initial crack, and the specimen enters the crack initiation stage. For the three different types of joint distribution, the initial cracks all appear at the joint ends. As the load continues to be applied, the tip of the initial crack becomes a new stress concentration point, and the cohesive element at the crack tip quickly reaches the stage of damage evolution, which is thus deleted. The crack continues to grow, and the specimen enters the crack evolution stage. For type-I, the joint ends are more prone to stretching, while for type-II and type-III, the joint ends are more subjected to shearing, thus the crack evolution of type-I joints is faster. Furthermore, under load, the cracks continue to propagate and eventually form the main fissure through the rock bridge, and the specimen enters the final failure stage. At this stage, the type-I joint produces the largest shear deformation, followed by the type-III joint, and the type-II joint has the smallest. Only type-I Firstly, when the shear load is just applied, the cohesive elements begin to accumulate fracture energy. At this time, the cohesive elements are in the linear elastic stage, and the traction force does not reach the initial damage stress, T 0 . No cracks appear on the specimen, which corresponds to the elastic stage in the stress-strain curve.
As the shear load continues to be applied, the cohesive elements accumulate more fracture energy. Due to the presence of joints inside the specimen, the stress concentration tends to occur at the end of the joint. As shown in Figure 15a, the stress gradually decreases from the joint end to other directions. The cohesive elements near the stress concentration point tend to accumulate more fracture energy, and thus reach the stage of crack evolution, as shown in Figure 15b. The stiffness degradation rate SDEG of the cohesive element at the joint end reaches 1.0 as the first, which is deleted to generate the initial crack, and the specimen enters the crack initiation stage. For the three different types of joint distribution, the initial cracks all appear at the joint ends. Firstly, when the shear load is just applied, the cohesive elements begin to accumulate fracture energy. At this time, the cohesive elements are in the linear elastic stage, and the traction force does not reach the initial damage stress, . No cracks appear on the specimen, which corresponds to the elastic stage in the stress-strain curve.
As the shear load continues to be applied, the cohesive elements accumulate more fracture energy. Due to the presence of joints inside the specimen, the stress concentration tends to occur at the end of the joint. As shown in Figure 15a, the stress gradually decreases from the joint end to other directions. The cohesive elements near the stress concentration point tend to accumulate more fracture energy, and thus reach the stage of crack evolution, as shown in Figure 15b. The stiffness degradation rate SDEG of the cohesive element at the joint end reaches 1.0 as the first, which is deleted to generate the initial crack, and the specimen enters the crack initiation stage. For the three different types of joint distribution, the initial cracks all appear at the joint ends. As the load continues to be applied, the tip of the initial crack becomes a new stress concentration point, and the cohesive element at the crack tip quickly reaches the stage of damage evolution, which is thus deleted. The crack continues to grow, and the specimen enters the crack evolution stage. For type-I, the joint ends are more prone to stretching, while for type-II and type-III, the joint ends are more subjected to shearing, thus the crack evolution of type-I joints is faster. Furthermore, under load, the cracks continue to propagate and eventually form the main fissure through the rock bridge, and the specimen enters the final failure stage. At this stage, the type-I joint produces the largest shear deformation, followed by the type-III joint, and the type-II joint has the smallest. Only type-I produces tensile through cracks, which is mainly attributed to the clockwise moment around the centroid of the specimen. For type-II and type-III, the penetrated crack extends along the rock bridge and eventually merges with the joint to form a new joint through the specimen. At this time, the shear resistance of the test piece consists of two parts, one is the mechanical engagement force from the As the load continues to be applied, the tip of the initial crack becomes a new stress concentration point, and the cohesive element at the crack tip quickly reaches the stage of damage evolution, which is thus deleted. The crack continues to grow, and the specimen enters the crack evolution stage. For type-I, the joint ends are more prone to stretching, while for type-II and type-III, the joint ends are more subjected to shearing, thus the crack evolution of type-I joints is faster. Furthermore, under load, the cracks continue to propagate and eventually form the main fissure through the rock bridge, and the specimen enters the final failure stage. At this stage, the type-I joint produces the largest shear deformation, followed by the type-III joint, and the type-II joint has the smallest. Only type-I produces tensile through cracks, which is mainly attributed to the clockwise moment around the centroid of the specimen. For type-II and type-III, the penetrated crack extends along the rock bridge and eventually merges with the joint to form a new joint through the specimen. At this time, the shear resistance of the test piece consists of two parts, one is the mechanical engagement force from the irregular protrusions on both sides of the crack, and the other is the friction resistance on the fracture surface. Through the comparison between the type-II and type-III, it can be found that more cracks are produced on the rock bridge of the type-II joint, and there are even completely free solid elements in the rock bridge on the left, from which it can be inferred that the fracture surface of the type-II rock bridge is rougher.

Conclusions
This study used the FEM-CZM method, established the symmetrical discontinuous joint shear model, and inserted cohesive elements into the solid elements to form potential crack surfaces. This method can not only show the mechanical deformation of the element but also have unique advantages in the simulation of complex cracks. Through the uniaxial compression test and the Brazilian tests, we determined the material parameters and the basis for meshing. Three types of discontinuous joints with different distribution forms were established, including both-sided, scattered, and central. The influence of the distribution form, joint persistence, and normal stress on the shear resistance of rock was investigated. The law of crack propagation during the shear process of the joints with three different distribution forms was analyzed. Based on this method, the following main conclusions can be obtained: (1) The shearing process can be divided into four stages: elastic stage, strengthening stage, plastic stage, and residual stress stage. In the stress-strain curve, the residual stress stage of type-II has a higher slope, which shows that the type-II has more brittleness and the type-I has higher plasticity. At the same time, with the decrease of joint persistence, the specimen turns more brittle, and the joint shear is closer to the direct shear test of intact rock. (2) Under the same conditions, the specimen in the type-I is more likely to produce an unbalanced moment around the centroid of the specimen, which causes the opposite ends of the specimen to produce vertical displacements in opposite directions. Due to the strengthening effect of the rock bridge, the vertical displacement on the side away from the loading site is smaller, and as a result, the overall displacement field of the specimen is not simply center-symmetric. At the same time, the distribution of rock bridges in type-II is more dispersed, providing more reinforcement for the joints, and the specimens are subjected to lower dilatancy. (3) The crack propagation process can be divided into three stages: crack initiation stage, crack evolution stage, and final failure stage. Under the load, the joint tips are prone to stress concentration, where the cohesive elements accumulate more fracture energy and reach the damage evolution stage faster. Initial cracks always start from the joint ends. Affected by the unbalanced moment, more tensile cracks are generated in type-I, while the type-II and type-III penetrated cracks are all shear cracks, and the cracks mainly propagate along the rock bridge.