# Discrete Element Modeling on Mechanical Behavior of Heterogeneous Rock Containing X-Shaped Fissure under Uniaxial Compression

Authors: Ming Chen; Jian Liu; Zhengyong Xie; Jianjun Liu; Xunjian Hu; Bingyang Li; Yue Cen
Journal: Geofluids
 (2020)

---

## Abstract

Based on the experimental results of an intact rock specimen under uniaxial compression, particle flow code (PFC2D) was adopted to carry out a discrete element modeling (DEM) for the mechanical behavior of heterogeneous rocks containing X-shaped fissures (two intersecting symmetric single fissures) under uniaxial compression. The influences of β (the acute angle between two single fissures) and the direction angle α (the acute angle between the bisector of β and perpendicular to the loading direction) on the strength, deformation, energy, crack propagation, and ultimate failure mode were analyzed in detail. Numerical simulated results showed the following: (1) Due to the X-shaped fissures, not only the peak strength, elastic modulus, crack initiation stress, and damage stress were significantly reduced, and the reduced degree of the peak strength was obviously greater than that of the elastic modulus, but also the brittleness and energy were significantly weakened. (2) The peak strength and elastic modulus generally decreased with the increase of β and increased with the increase of α. Moreover, the change trends of crack initiation stress, damage stress, boundary energy, and total strain energy at the peak stress were consistent with the peak strength. (3) Regardless of the changes of α and β, models all firstly initiated wing cracks at the two tips of the single fissure with a larger inclination angle, and the crack initiation angle decreased with the increase of the inclination angle of the single fissure. (4) The fracture was dominated by tensile microcracks, and no microcracks were generated in a certain range of the X-shaped fissure center. The failure mode was mainly split along the axial direction, and the failure surface started from the tips of the fissure and extended to both ends of models. (5) The uniaxial compressive strength and elastic modulus increased exponentially with the increase of the homogeneity factor. When the homogeneity factor was small, the microcracks were more evenly distributed in the models; when the homogeneity factor was large, the microcracks were mainly concentrated at the tips of the fissure in the models. This study can provide some references for the correct understanding of the mechanical properties of rock masses containing X-shaped fissures.


---

## Body

## 1. Introduction

Due to the long-term geological tectonics and human activities, there are often many natural fissures in the rock mass, which will degrade the mechanical properties of rock mass [1–12] and lead to engineering accidents and disasters. Therefore, the research on the mechanical behavior of fissured rocks has always been the focus of the geotechnical engineering community.At present, many scholars at home and abroad have studied the influence of the shape, size, number, and distribution characteristics of the fissures on the mechanical properties, crack propagation, and failure modes of the rock through experiments and numerical simulations [6, 13–25]. For example, Wang YX et al. [13] have carried out uniaxial compression experiments on the granite specimens with a central hole and hole edge flaws, and the results showed that the mechanical parameters and crack types of the specimens distinctly vary with the flaw inclination. Wang CL et al. [15] have studied the propagation of collinear cracks and cracking processes under uniaxial compression by experiment and numerical simulation. Liu XR et al. [16] have studied the failure behavior and fracture mechanism of sandstone containing elliptical holes and fissures through experiments, and the results showed that the ligament angle has a significant impact on rock strength and deformation. Wu JY et al. [17] have studied the length of preexisting fissures effects on the strength and deformation of red sandstone through experiments, and the results showed that as the fissure length increases, the strength and deformation parameters decrease in a negative exponential manner. Dong QQ et al. [18] and Ma GW et al. [19] have shown that the inclination angle and effective curvature of the S-shaped fissure have a great influence on the crack propagation and failure mode of the specimen through experiments and numerical simulations. Yang SQ et al. [20] have conducted experiments and numerical simulations on red sandstone containing a single nonstraight fissure, and the results indicated that as the inclination angle increases, the peak strength, peak strain, and elastic modulus increase as a whole. Moreover, when the inclination angle is smaller or larger, the crack first initiated from the convex point of the prefabricated fissure; when the inclination angle is intermediate, the crack initiation position transferred from the convex point to the tip of the prefabricated fissure. Yang SQ et al. [21] have showed that the oval flaw angle has an effect on the macroscopic strength parameters and crack modes of sandstone through experiments and numerical simulation. Chen LY et al. [22] have analyzed the influence of fissure dips and rock bridge dips on strength, crack propagation, and evolution laws through numerical simulation. Yang SQ et al. [23] have studied the fracture coalescence behavior of red sandstone containing two nonparallel fissures by PFC2D, and the results showed that as the fissure angle increases, the peak strength and elastic modulus first increase and then decrease. Wang GL et al. [25] have analyzed the influence of the inclination angle of the Z-shaped fissure on the strength, deformation, and crack propagation, and the results indicated that the failure modes are mainly divided into tension failure, shear failure, and tensile-shear mixed failure.It can be seen from the above research that the current research mainly focuses on single fissures, double fissures, elliptical fissures, arc-shaped fissures, and Z-shaped fissures. However, the forms of fissures are complex and changeable, and X-shaped fissures are also one of the widely existing fissure forms in actual rock mass engineering, as shown in Figure1. In addition, there are currently few studies on rocks containing X-shaped fissures. There are the following problems in the experiment: (1) It is extremely difficult to prepare a rock sample, especially hard rock. (2) The heterogeneity of the rock leads to large dispersion of experimental results. (3) It is not easy to directly observe the crack propagation inside the rock. However, numerical simulation can solve the above problems. In view of this, this paper intended to use particle flow code (PFC2D) [26] to carry out the numerical simulation for the mechanical behavior and failure characteristics of heterogeneous rock specimens containing X-shaped fissures under uniaxial compression. This study can provide some references for the correct understanding of the mechanical properties of rock masses containing X-shaped fissures.Figure 1
The X-shaped fissures in rock mass.


## 2. Discrete Element Modeling Method

### 2.1. Particle Flow Code (PFC2D)

Based on the discrete element method, Cundall and Strack [27] established the particle flow theory in 1979, and later, the United States Itasca company developed the particle flow program commercial software PFC. In this research, we chose the parallel bond model (PBM) to carry out the numerical simulation, as shown in Figure 2 [28]. The particles are cemented together by parallel bonds, so the relative motion between the particles will cause a force and a moment to develop within the bond material, and this force and moment act on the two bonded particles. The forces and moments can be decomposed into F¯in and M¯in along the normal direction and F¯is and M¯is along the tangential direction. When the normal stress generated in the parallel bond is greater than the normal strength of the bond material, the tensile failure of the parallel bond generates a tensile crack; when the tangential stress generated in the parallel bond is greater than the shear strength of the bond material, the shear failure of the parallel bond generates a shear crack.Figure 2
Mechanical mechanism of parallel bond model [28].


### 2.2. Calibrating Microparameters by Experimental Results of Intact Specimen

The intact specimens were taken from the tuff in the volcanic reservoir of the Xushen Gas Field in the experiment. First, the intact specimens were processed into a cylindrical standard specimen with a diameter of 25 mm and height of 50 mm and then were loaded with a displacement control method at a rate of 0.005 mm/s.Figure3 showed one intact numerical model generated by PFC2D, and the scale of the numerical model was similarly equal to that of the experimental specimen, i.e., the dimensions of 25 mm in width and 50 mm in height. Each intact numerical model was discretized into 16542 particles and 39379 parallel bonded contacts. By giving the upper and lower walls a certain speed, quasistatic loading was achieved for simulating uniaxial compression experiments.Figure 3
Intact numerical model generated in PFC2D in this research.Due to the differences in the size, shape, mineral types of the rock-forming mineral particles, and diagenesis in the rock, its mechanical properties showed great heterogeneity. Many scholars had adopted the Weibull distribution to characterize the heterogeneity of mechanical parameters in numerical simulations [29, 30]. Therefore, this paper set the values of the parallel bond tensile strength and parallel bond cohesion in the particle flow model to obey the Weibull distribution, and the probability density function of fλ is
(1)fλ=mλ¯λλ¯m−1e−λ/λ¯m,λ≥0,where m is the homogeneity factor (the smaller the value of m, the greater the heterogeneity of the model), λ is the parallel bond tensile strength and parallel bond cohesion, and λ¯ is the characteristic value of parallel bond tensile strength and parallel bond cohesion. For example, Figure 4 shows the density distribution of the parallel bond tensile strength, with the values of λ¯ and m being 78.6 MPa and 5, respectively.Figure 4
The density distribution of the parallel bond tensile strength (λ¯=78.6MPa,m=5).During the calibration process, microparameters were confirmed by using the trial and error method. The macroscopic behavior (i.e., the peak strength, elastic modulus, Poisson’s ratio, the stress-strain curve, and the failure mode) of an intact specimen obtained by experiment was used in this research to calibrate the microparameters. The macroscopic results obtained by numerical simulation after each trial were used to check the microparameters. This process was repeated until the numerical results achieved a good agreement with the experimental results. Table1 lists the microparameters used in the PFC2D model in this research. Table 2 shows the comparison between experimental and numerical mechanical parameters of the intact specimen. In accordance with Table 2, we can conclude that the relative error of mechanical parameters was very small between numerical results and experimental results. Besides, Figure 5 shows the comparison between experimental and numerical stress-strain curves and failure mode of the intact specimen under uniaxial compression. From Figure 5, it can be seen that the stress-strain curve and the failure mode were well in an agreement between numerical results and experimental results. The comparisons shown in Table 2 and Figure 5 calibrated the rightness and reasonability of microparameters used in Table 1. Later, it can be used for the simulation of the numerical models containing X-shaped fissure.Table 1
The microparameters used in the PFC2D model.
Microparameters
Values
Particle density (kg·m-3)
2560
Minimum particle size (mm)
0.11
Particle size ratio
1.66
Particle contact modulus (GPa)
14
Particle stiffness ratio
3.0
Particle friction coefficient
1.2
Homogeneity factor
5
Parallel bonding modulus (GPa)
14
Parallel bonding stiffness ratio
3.0
Characteristic value of parallel bond tensile strength (MPa)
78.6
Characteristic value of parallel bond cohesion (MPa)
138.6
Friction angle of parallel bond (°)
30.0Table 2
Comparison between experimental and numerical mechanical parameters of the intact specimen.
Peak strength (MPa)
Elastic modulus (GPa)
Poisson’s ratio
Numerical results
183.0
20.49
0.26
Experimental results
176.7
20.66
0.25
Relative error
3.57%
0.82%
4.00%Figure 5
Comparison between experimental and numerical stress-strain curves and failure mode of the intact specimen.


### 2.3. Establishment of Numerical Models Containing X-Shaped Fissure

After generating the intact numerical model, we created open X-shaped fissures by deleting particles [23]. The layout of the X-shaped fissure is shown in Figure 6. The intersection point of fissure ① and fissure ② was located in the center of the model. The length and width of fissure ① and fissure ② were both 10 mm and 0.5 mm, respectively. β was the acute angle between fissure ① and fissure ②, and α was the acute angle between bisector of β and perpendicular to the loading direction. In order to investigate in detail the influence of α and β on the mechanical behavior of rock containing X-shaped fissure, we designed the following two simulation schemes: (1) α was constant at 45°, and β was 15°, 30°, 45°, 60°, 75°, and 90°, respectively, and (2) β was constant at 60°, and α was 0°, 15°, 30°, 45°, 60°, 75°, and 90°, respectively.Figure 6
The layout of X-shaped fissure.


## 2.1. Particle Flow Code (PFC2D)

Based on the discrete element method, Cundall and Strack [27] established the particle flow theory in 1979, and later, the United States Itasca company developed the particle flow program commercial software PFC. In this research, we chose the parallel bond model (PBM) to carry out the numerical simulation, as shown in Figure 2 [28]. The particles are cemented together by parallel bonds, so the relative motion between the particles will cause a force and a moment to develop within the bond material, and this force and moment act on the two bonded particles. The forces and moments can be decomposed into F¯in and M¯in along the normal direction and F¯is and M¯is along the tangential direction. When the normal stress generated in the parallel bond is greater than the normal strength of the bond material, the tensile failure of the parallel bond generates a tensile crack; when the tangential stress generated in the parallel bond is greater than the shear strength of the bond material, the shear failure of the parallel bond generates a shear crack.Figure 2
Mechanical mechanism of parallel bond model [28].


## 2.2. Calibrating Microparameters by Experimental Results of Intact Specimen

The intact specimens were taken from the tuff in the volcanic reservoir of the Xushen Gas Field in the experiment. First, the intact specimens were processed into a cylindrical standard specimen with a diameter of 25 mm and height of 50 mm and then were loaded with a displacement control method at a rate of 0.005 mm/s.Figure3 showed one intact numerical model generated by PFC2D, and the scale of the numerical model was similarly equal to that of the experimental specimen, i.e., the dimensions of 25 mm in width and 50 mm in height. Each intact numerical model was discretized into 16542 particles and 39379 parallel bonded contacts. By giving the upper and lower walls a certain speed, quasistatic loading was achieved for simulating uniaxial compression experiments.Figure 3
Intact numerical model generated in PFC2D in this research.Due to the differences in the size, shape, mineral types of the rock-forming mineral particles, and diagenesis in the rock, its mechanical properties showed great heterogeneity. Many scholars had adopted the Weibull distribution to characterize the heterogeneity of mechanical parameters in numerical simulations [29, 30]. Therefore, this paper set the values of the parallel bond tensile strength and parallel bond cohesion in the particle flow model to obey the Weibull distribution, and the probability density function of fλ is
(1)fλ=mλ¯λλ¯m−1e−λ/λ¯m,λ≥0,where m is the homogeneity factor (the smaller the value of m, the greater the heterogeneity of the model), λ is the parallel bond tensile strength and parallel bond cohesion, and λ¯ is the characteristic value of parallel bond tensile strength and parallel bond cohesion. For example, Figure 4 shows the density distribution of the parallel bond tensile strength, with the values of λ¯ and m being 78.6 MPa and 5, respectively.Figure 4
The density distribution of the parallel bond tensile strength (λ¯=78.6MPa,m=5).During the calibration process, microparameters were confirmed by using the trial and error method. The macroscopic behavior (i.e., the peak strength, elastic modulus, Poisson’s ratio, the stress-strain curve, and the failure mode) of an intact specimen obtained by experiment was used in this research to calibrate the microparameters. The macroscopic results obtained by numerical simulation after each trial were used to check the microparameters. This process was repeated until the numerical results achieved a good agreement with the experimental results. Table1 lists the microparameters used in the PFC2D model in this research. Table 2 shows the comparison between experimental and numerical mechanical parameters of the intact specimen. In accordance with Table 2, we can conclude that the relative error of mechanical parameters was very small between numerical results and experimental results. Besides, Figure 5 shows the comparison between experimental and numerical stress-strain curves and failure mode of the intact specimen under uniaxial compression. From Figure 5, it can be seen that the stress-strain curve and the failure mode were well in an agreement between numerical results and experimental results. The comparisons shown in Table 2 and Figure 5 calibrated the rightness and reasonability of microparameters used in Table 1. Later, it can be used for the simulation of the numerical models containing X-shaped fissure.Table 1
The microparameters used in the PFC2D model.
Microparameters
Values
Particle density (kg·m-3)
2560
Minimum particle size (mm)
0.11
Particle size ratio
1.66
Particle contact modulus (GPa)
14
Particle stiffness ratio
3.0
Particle friction coefficient
1.2
Homogeneity factor
5
Parallel bonding modulus (GPa)
14
Parallel bonding stiffness ratio
3.0
Characteristic value of parallel bond tensile strength (MPa)
78.6
Characteristic value of parallel bond cohesion (MPa)
138.6
Friction angle of parallel bond (°)
30.0Table 2
Comparison between experimental and numerical mechanical parameters of the intact specimen.
Peak strength (MPa)
Elastic modulus (GPa)
Poisson’s ratio
Numerical results
183.0
20.49
0.26
Experimental results
176.7
20.66
0.25
Relative error
3.57%
0.82%
4.00%Figure 5
Comparison between experimental and numerical stress-strain curves and failure mode of the intact specimen.


## 2.3. Establishment of Numerical Models Containing X-Shaped Fissure

After generating the intact numerical model, we created open X-shaped fissures by deleting particles [23]. The layout of the X-shaped fissure is shown in Figure 6. The intersection point of fissure ① and fissure ② was located in the center of the model. The length and width of fissure ① and fissure ② were both 10 mm and 0.5 mm, respectively. β was the acute angle between fissure ① and fissure ②, and α was the acute angle between bisector of β and perpendicular to the loading direction. In order to investigate in detail the influence of α and β on the mechanical behavior of rock containing X-shaped fissure, we designed the following two simulation schemes: (1) α was constant at 45°, and β was 15°, 30°, 45°, 60°, 75°, and 90°, respectively, and (2) β was constant at 60°, and α was 0°, 15°, 30°, 45°, 60°, 75°, and 90°, respectively.Figure 6
The layout of X-shaped fissure.


## 3. Simulation Results

### 3.1. Strength and Deformation Characteristics

When stress drops to 70% of the peak stress, the calculation was stopped in the numerical simulation in this research. From the stress-strain curve in Figure7, the following can be seen: (1) The stress-strain relationship of the intact numerical model and the models containing X-shaped fissure was both nearly linear. Before the peak stress, it was mainly manifested as elastic deformation, with little plastic deformation. After the peak stress, the stress dropped rapidly, showing that the tuff had typical brittle failure characteristics. (2) Compared with the intact numerical model, the elastic stage of numerical models containing X-shaped fissure was significantly shorter, and parts of specimens had obvious stress drops before the peak stress, which was due to local damage caused by new or existing crack propagation and penetration. (3) Compared with the intact numerical model, the numerical models containing X-shaped fissure had a smaller stress drop rate after the peak stress, indicating that the X-shaped fissure weakened the brittleness of the rock and increased the ductility. The above conclusions were consistent with the numerical simulation results of other scholars [25].Figure 7
Stress-strain curves of numerical models.
(a)
α
=
45
°
(b)
β
=
60
°The peak strength and elastic modulus obtained from the stress-strain curves in Figure6 are shown in Table 3. It can be seen that compared with the intact numerical model, the peak strength and elastic modulus of the numerical models containing X-shaped fissure were significantly reduced, indicating that the X-shaped fissure has greatly degraded the mechanical properties of rock. The damage coefficient of ηwas used to quantitatively describe the degree of deterioration, which is determined by the following formula [1]:
(2)η=δ1−δ2δ1,where δ1 and δ2 are the mechanical parameters of the intact numerical model and the models containing X-shaped fissure, respectively. The damage coefficient calculated by the above formula is shown in Table 3. The following can be seen: (1) The peak strength damage coefficient of numerical models containing X-shaped fissure was between 0.31 and 0.53, and the elastic modulus damage coefficient was between 0.08 and 0.21. The peak strength damage coefficient was about 3 times the elastic modulus damage coefficient. It means that compared with the intact numerical model, the X-shaped fissure had a greater degree of deterioration to the peak strength, which was much greater than that to the elastic modulus. This conclusion was consistent with the experimental and numerical simulation results of other scholars [14]. (2) The peak strength and elastic modulus damage coefficients caused by α were between 0.31~0.53 and 0.08~0.21, respectively, while the peak strength and elastic modulus damage coefficients caused by β were between 0.40~0.52 and 0.12~0.18, respectively. Among them, the damage coefficient of the former had a larger variation range than the latter, indicating that α had a greater influence on mechanical parameters than β.Table 3
Mechanical parameters of numerical models.
α (°)
β (°)
Peak strength (MPa)
Elastic modulus (GPa)
Peak strength damage coefficient
Elastic modulus damage coefficient
Intact
183.00
20.49
—
—
45
15
110.28
17.73
0.40
0.13
45
30
104.00
17.94
0.43
0.12
45
45
91.90
17.49
0.50
0.15
45
60
95.08
17.02
0.48
0.17
45
75
96.06
17.54
0.48
0.14
45
90
87.17
16.78
0.52
0.18
0
60
89.48
16.22
0.51
0.21
15
60
86.35
16.83
0.53
0.18
30
60
95.03
16.99
0.48
0.17
45
60
95.08
17.02
0.48
0.17
60
60
101.95
17.92
0.44
0.13
75
60
104.13
18.78
0.43
0.08
90
60
126.95
18.84
0.31
0.08Figure8 shows the curves of mechanical parameters changing with α and β. The following can be seen: (1) With the changes of α and β, the peak strength changed greatly, while the elastic modulus changed less. This was because the elastic modulus was mainly related to the number of fissures; the more the number of fissures, the smaller the elastic modulus [1]. In this paper, the number of fissures was fixed, so the change of elastic modulus was small. (2) With the increase of β, the peak strength and elastic modulus showed a decreasing trend as a whole. Specifically, the peak strength first decreased, then increased and then decreased. When β=15° and β=30°, the peak strength was larger, which was equal to 110.28 MPa and 104.00 MPa, respectively; when β=45° and β=90°, the peak strength was smaller, which was equal to 91.90 MPa and 87.17 MPa, respectively. However, the elastic modulus first increased, then decreased, then increased, and finally decreased. The elastic modulus reached the maximum value of 17.94 GPa at α=30°, while the elastic modulus reached the minimum value of 16.78 GPa at α=90°. (3) With the increase of α, the peak strength and elastic modulus tended to increase as a whole. Specifically, the peak strength increased slowly at α=0°~75° and increased rapidly at α=75°~90°. The peak strength reached the maximum value of 126.95 MPa at α=90°, while the peak strength reached the minimum value of 86.35 MPa at α=15°. However, the elastic modulus growed in an S-shaped curve. The growth was gentle at α=0°~45° and α=75°~90° and rapid at α=45°~75°. When α=90°, the elastic modulus reached the maximum value of 18.84 GPa. When α=0° and α=15°, the elastic modulus was smaller, which was equal to 16.22 GPa and 16.83 GPa, respectively.Figure 8
Change curves of mechanical parameter of numerical models.
(a)
α
=
45
°
(b)
β
=
60
°In summary, the peak strength and elastic modulus decreased with the increase ofβ and increased with the increase of α on the whole. The reason was that as β decreased or α increased, the projection of the fissure length in the direction perpendicular to the loading direction gradually becomes shorter, leading to the gradual reduction of the axial closure effect of the fissure under the axial load, and the effective bearing area of the rock actually bearing the load gradually increased, which caused that the strength and the stiffness gradually increased [25].


### 3.2. Crack Initiation Stress and Damage Stress

The crack initiation stress and damage stress of rock are of great significance for understanding the progressive failure process of rock and establishing its brittle failure mechanism, initiation strength, and long-term strength criterion [31]. Through the method of monitoring microcracks, taking the intact numerical model as an example, Figure 9 shows the change curve of the microcrack number in the intact numerical model. The axial stress corresponding to the crack initiation point was taken as the initiation stress (σci), and the axial stress corresponding to the rapid crack growth point was taken as the damage stress (σcd) [32].Figure 9
Change curve of the microcrack number in intact numerical model.Table4 showed the initiation stress and damage stress of all numerical models determined by the above method. The following can be seen: (1) Compared with the intact numerical model, the crack initiation stress and damage stress of the numerical models containing X-shaped fissure were significantly reduced, but the level of crack initiation stress (the ratio of the initiation stress to the peak stress (σp)) and the level of damage stress (the ratio of damage stress to peak stress) did not change much. (2) The level of crack initiation stress in all numerical models was between 0.36 and 0.5, and the level of damage stress was between 0.76 and 0.84. Studies [33] have shown that the crack initiation stress of brittle rocks is about 0.25 to 0.5 times the peak stress, and the damage stress is about 0.7 to 0.85 times the peak stress.Table 4
The crack initiation stress and damage stress of numerical models.
α (°)
β (°)
σ
c
i (MPa)
σ
c
d (MPa)
σ
p (MPa)
σ
c
i
/
σ
p
σ
c
d
/
σ
p
Intact
74.2
153.8
183.0
0.41
0.84
45
15
44.5
90.3
110.3
0.40
0.82
45
30
43.0
82.7
104.0
0.41
0.80
45
45
36.8
74.4
91.9
0.40
0.81
45
60
45.7
76.5
95.1
0.48
0.80
45
75
44.0
73.4
96.1
0.46
0.76
45
90
43.8
70.8
87.2
0.50
0.81
0
60
37.5
74.4
89.5
0.42
0.83
15
60
34.7
71.3
86.4
0.40
0.83
30
60
34.1
79.5
95.0
0.36
0.84
45
60
45.7
76.5
95.1
0.48
0.80
60
60
49.3
83.6
102.0
0.48
0.82
75
60
43.4
84.2
104.1
0.42
0.81
90
60
55.2
102.5
127.0
0.44
0.81Figure10 shows the curves of the crack initiation stress and damage stress with α and β. From Figure 10(a), the following can be seen: (1) With the increase of β, both the crack initiation stress and the damage stress decreased first, then increased, and then decreased, which was consistent with the peak strength. When β=45°, both the crack initiation stress and the damage stress reached the minimum value. (2) With the increase of β, the level of crack initiation stress and the level of damage stress had no obvious changes, which were between 0.4 and 0.5 and between 0.76 and 0.82, respectively. When β=45°, the level of crack initiation stress reached the minimum value of 0.4, and at this time, cracks were most prone to initiation.Figure 10
Change curves of crack initiation stress and damage stress of numerical models containing X-shaped fissure.
(a)
α
=
45
°
(b)
β
=
60
°From Figure10(b), the following can be seen: (1) With the increase of α, the crack initiation stress and damage stress increased as a whole, which was consistent with the peak strength. (2) With the increase of α, the level of crack initiation stress and the level of damage stress had no obvious changes, which were between 0.36 and 0.48 and between 0.80 and 0.84, respectively. When α=30°, the level of crack initiation stress reached the minimum value of 0.36, and at this time, cracks were most prone to initiation.In comparison, with the changes ofα and β, the range of the level of damage stress was smaller than the level of crack initiation stress, indicating that α and β had a greater impact on the level of crack initiation stress than on the level of damage stress.


### 3.3. Law of Energy Evolution

The essence of rock failure is a state instability phenomenon driven by energy [34], and the study of energy from the perspective of mesomechanics is more in line with actual energy changes. In PFC simulation, the sum of particle bonding energy and strain energy is recorded as total strain energy, and the difference between boundary energy and total strain energy is recorded as dissipation energy [35]. Due to space limitations, the energy evolution law during the loading process was analyzed by taking the numerical model of α=45° and β=60° as an example, as shown in Figure 11, and other numerical models had similar laws.Figure 11
Relationship between energy and stress-strain of numerical model (α=45°, β=60°).From Figure11, the following can be seen: (1) At the initial stage of loading (from point O to point A, which was the crack initiation stress point), there were no microcracks, and the boundary energy and the total strain energy increased nonlinearly and basically coincide, indicating that the input boundary energy was all converted into total strain energy, and the dissipation energy, kinetic energy, and friction energy were almost zero. (2) At the linear elastic stage of loading (from point A to point B, which was the damage stress point), the microcracks initiated and growed steadily. The growth rate of the boundary energy and the total strain energy gradually increased, and the dissipated energy and the friction energy increased slowly. At point B, there was a small separation between the boundary energy and the total strain energy curve, and the strain energy stored in contact during the crack propagation process was gradually transformed into dissipated energy. (3) At the plastic yielding stage of loading (from point B to point C, which was the peak stress point), the microcracks expand unsteadily, and the boundary energy and the total strain energy curve gradually separated. The total strain energy reached the peak at point C, and the growth of dissipation energy gradually increased. (4) After the peak stress (from point C to point D), the boundary energy continued to increase, the stored total strain energy was rapidly released, and the dissipation energy increased sharply until the model was broken. During the entire deformation and failure process, the boundary energy keeps increasing, the total strain energy first increased and then decreased sharply, and the dissipation energy increased slowly and then increased sharply.Table5 shows the energy index of all numerical models at the peak stress. From Table 5, the following can be seen: (1) The total strain energy of all numerical models accounted for a large proportion of the boundary energy at the peak point, indicating that the input energy before the peak stress was mainly stored in the rock in the form of strain energy. (2) Compared with the intact numerical model, the boundary energy, total strain energy, and dissipation energy of the numerical models containing X-shaped fissure were significantly reduced, and the largest decline reached 71.7%, 72.5%, and 71.1%, respectively. It showed that the X-shaped fissure had obvious energy weakening effect on the rock. The total strain energy at the peak stress represents the energy storage limit of the rock and reflects the ability of the rock to resist damage, so the rock containing fissure was more likely to be damaged.Table 5
Energy indexes in the peak stress of numerical models.
α (°)
β (°)
Boundary energy (J)
Total strain energy (J)
Dissipated energy (J)
Intact
1070.8
1035.0
35.7
45
15
439.8
418.9
20.9
45
30
402.1
386.7
15.4
45
45
325.9
308.2
17.7
45
60
356.2
338.3
17.9
45
75
383.5
360.5
23.0
45
90
303.4
284.6
18.8
0
60
371.4
332.6
38.7
15
60
332.8
301.8
30.9
30
60
364.3
342.0
22.3
60
60
391.9
373.6
18.2
75
60
373.3
363.0
10.3
90
60
605.5
567.9
37.6Figure12 shows the curves of the energy indexes at peak stress with α and β. The following can be seen: (1) With the increase of β, the boundary energy and total strain energy showed a decreasing trend as a whole. Specifically, it first decreased, then increased and then decreased, which was consistent with the changing law of peak stress in Figure 8(a). (2) With the increase of α, the boundary energy and the total strain energy showed an overall increasing trend, which was consistent with the change law of the peak stress in Figure 8(b).Figure 12
Change curves of energy indexes of numerical models containing X-shaped fissure in the peak stress.
(a)
α
=
45
°
(b)
β
=
60
°


### 3.4. Crack Propagation Process and Failure Mode

Due to space limitations, the crack propagation process during the loading process was analyzed by taking the numerical model ofα=45° and β=60° as an example, as shown in Figure 13. Figures 14 and 15, respectively, show the failure modes of numerical models with different α and β. The red and black in Figures 13 and 14 represented tensile microcracks and shear microcracks, respectively, and the marked numbers indicated the sequence of macrocracks.Figure 13
Crack propagation process during the loading process of the numerical model containing X-shaped fissure (α=45°, β=60°). (Red represents tensile microcracks. Black represents shear microcracks. Numbers indicate the order of macrocrack generation. The angle value indicates the crack initiation angle).
(a)
σ
=
56.3
MPa
,
ε
=
3.3
×
10
−
3
(b)
σ
=
75.8
MPa
,
ε
=
4.4
×
10
−
3
(c)
σ
=
84.8
MPa
,
ε
=
5.2
×
10
−
3
(d)
σ
=
94.5
MPa
,
ε
=
5.6
×
10
−
3
(e)
σ
=
88.4
MPa
,
ε
=
5.8
×
10
−
3
(f)
σ
=
66.6
MPa
,
ε
=
6.0
×
10
−
3Figure 14
Failure modes of numerical models with different values ofβ (α=45°). (Red represents tensile microcracks. Black represents shear microcracks. Numbers indicate the order of macrocrack generation. The angle values indicate the crack initiation angle).
(a)
β
=
15
°
(b)
β
=
30
°
(c)
β
=
45
°
(d)
β
=
60
°
(e)
β
=
75
°
(f)
β
=
90
°Figure 15
Failure modes of numerical models with different values ofα (β=60°). (Red represents tensile microcracks. Black represents shear microcracks. Numbers indicate the order of macrocrack generation. The angle values indicate the crack initiation angle).
(a)
α
=
0
°
(b)
α
=
15
°
(c)
α
=
30
°
(d)
α
=
45
°
(e)
α
=
60
°
(f)
α
=
75
°
(g)
α
=
90
°From Figure13, it can be seen that when the load reached a certain stress level, the model first initiated wing cracks 1a and 1b at the tip of fissure ② along the extension direction of fissure ② at an inclination angle of 30°, which was the crack initiation angle [36]. Continuing to load, the wing cracks 1a and 1b propagated to the upper and lower ends of the model along the loading direction, and macrocracks 2, 3, and 4 were generated at the tip of fissure ①, which eventually led to the overall failure of the model.From Figure14, the following can be seen: (1) Wing cracks 1a and 1b were generated at the two tips of fissure ②. When β=15°, the crack initiation angle was 90°, and when β=90°, the crack initiation angle was 0°. Moreover, with the increase of β, the crack initiation angle decreased. (2) Macrocrack 2 was produced at the two tips of fissure ①. When β=15° and 60°, macrocrack 2 was produced at the left tip of fissure ①, and when β=30°, 45°, 75°, and 90°, macrocracks were generated at the right tip of fissure ①. (3) Macrocracks 3 and 4 were both produced at the two tips of fissure ① or at the ends of existing macrocracks. (4) When β=15°, 30°, 60°, and 75°, the four tips of the X-shaped fissure all had macrocracks and the distribution of microcracks was relatively dispersed, while when β=45° and 90°, only three tips had macrocracks and microcracks were concentrated at the fissure tips. Since the stress concentration of the former was not as high as that of the latter, the latter will fail at a lower stress level, which was also the reason why the peak strength of the model was small at β=45° and 90°.From Figure15, the following can be seen: (1) Except for α=0° and 90°, the numerical models first initiated wing cracks 1a and 1b at the two tips of fissure ②. When α=0°, the crack initiation angle was 90°; when α=90°, the crack initiation angle was 67°; and when α=60°, the crack initiation reached the minimum value of 0°. Moreover, with the increase of α, the crack initiation angle first decreased and then increased. (2) Macrocracks 2, 3, and 4 were all generated at the tips of fissure ① and fissure ② or at the end of the macrocracks that had been generated. (3) When α=0° and 15°, cracks were formed at the tips of fissure ① and fissure ②, so the peak strength was small. When α=90°, the four tips of the X-shaped fissure all had macrocracks and the distribution of microcracks was relatively scattered, so the peak strength was large.From Figures14 and 15, the following can be seen: (1) Tensile microcracks were the main cracks in the numerical models, and shear microcracks were almost rare. However, there were no microcracks within a certain range of the X-shaped fissure center. (2) The failure mode was mainly split along the axial direction, and the failure surface started from the tips of the fissure and extended to both ends of the numerical models.In summary, with the changes ofα and β, the numerical models first initiated wing cracks 1a and 1b at the two tips of fissure ② with a larger dip, and the crack initiation angle decreased with the increase of the inclination angle of fissure ②. From the theory of crack initiation, the relationship between the crack initiation angle of θ0 of a single-opening fissure tip and the inclination angle of α under uniaxial compression is as follows [36]:
(3)θ0=2arctan1+1+8sin2α/1+cos2α24sin2α/1+cos2α.The above formula shows that the crack initiation angle ofθ0 decreases with the increase of the inclination angle of α, which was consistent with the conclusion of this paper.


## 3.1. Strength and Deformation Characteristics

When stress drops to 70% of the peak stress, the calculation was stopped in the numerical simulation in this research. From the stress-strain curve in Figure7, the following can be seen: (1) The stress-strain relationship of the intact numerical model and the models containing X-shaped fissure was both nearly linear. Before the peak stress, it was mainly manifested as elastic deformation, with little plastic deformation. After the peak stress, the stress dropped rapidly, showing that the tuff had typical brittle failure characteristics. (2) Compared with the intact numerical model, the elastic stage of numerical models containing X-shaped fissure was significantly shorter, and parts of specimens had obvious stress drops before the peak stress, which was due to local damage caused by new or existing crack propagation and penetration. (3) Compared with the intact numerical model, the numerical models containing X-shaped fissure had a smaller stress drop rate after the peak stress, indicating that the X-shaped fissure weakened the brittleness of the rock and increased the ductility. The above conclusions were consistent with the numerical simulation results of other scholars [25].Figure 7
Stress-strain curves of numerical models.
(a)
α
=
45
°
(b)
β
=
60
°The peak strength and elastic modulus obtained from the stress-strain curves in Figure6 are shown in Table 3. It can be seen that compared with the intact numerical model, the peak strength and elastic modulus of the numerical models containing X-shaped fissure were significantly reduced, indicating that the X-shaped fissure has greatly degraded the mechanical properties of rock. The damage coefficient of ηwas used to quantitatively describe the degree of deterioration, which is determined by the following formula [1]:
(2)η=δ1−δ2δ1,where δ1 and δ2 are the mechanical parameters of the intact numerical model and the models containing X-shaped fissure, respectively. The damage coefficient calculated by the above formula is shown in Table 3. The following can be seen: (1) The peak strength damage coefficient of numerical models containing X-shaped fissure was between 0.31 and 0.53, and the elastic modulus damage coefficient was between 0.08 and 0.21. The peak strength damage coefficient was about 3 times the elastic modulus damage coefficient. It means that compared with the intact numerical model, the X-shaped fissure had a greater degree of deterioration to the peak strength, which was much greater than that to the elastic modulus. This conclusion was consistent with the experimental and numerical simulation results of other scholars [14]. (2) The peak strength and elastic modulus damage coefficients caused by α were between 0.31~0.53 and 0.08~0.21, respectively, while the peak strength and elastic modulus damage coefficients caused by β were between 0.40~0.52 and 0.12~0.18, respectively. Among them, the damage coefficient of the former had a larger variation range than the latter, indicating that α had a greater influence on mechanical parameters than β.Table 3
Mechanical parameters of numerical models.
α (°)
β (°)
Peak strength (MPa)
Elastic modulus (GPa)
Peak strength damage coefficient
Elastic modulus damage coefficient
Intact
183.00
20.49
—
—
45
15
110.28
17.73
0.40
0.13
45
30
104.00
17.94
0.43
0.12
45
45
91.90
17.49
0.50
0.15
45
60
95.08
17.02
0.48
0.17
45
75
96.06
17.54
0.48
0.14
45
90
87.17
16.78
0.52
0.18
0
60
89.48
16.22
0.51
0.21
15
60
86.35
16.83
0.53
0.18
30
60
95.03
16.99
0.48
0.17
45
60
95.08
17.02
0.48
0.17
60
60
101.95
17.92
0.44
0.13
75
60
104.13
18.78
0.43
0.08
90
60
126.95
18.84
0.31
0.08Figure8 shows the curves of mechanical parameters changing with α and β. The following can be seen: (1) With the changes of α and β, the peak strength changed greatly, while the elastic modulus changed less. This was because the elastic modulus was mainly related to the number of fissures; the more the number of fissures, the smaller the elastic modulus [1]. In this paper, the number of fissures was fixed, so the change of elastic modulus was small. (2) With the increase of β, the peak strength and elastic modulus showed a decreasing trend as a whole. Specifically, the peak strength first decreased, then increased and then decreased. When β=15° and β=30°, the peak strength was larger, which was equal to 110.28 MPa and 104.00 MPa, respectively; when β=45° and β=90°, the peak strength was smaller, which was equal to 91.90 MPa and 87.17 MPa, respectively. However, the elastic modulus first increased, then decreased, then increased, and finally decreased. The elastic modulus reached the maximum value of 17.94 GPa at α=30°, while the elastic modulus reached the minimum value of 16.78 GPa at α=90°. (3) With the increase of α, the peak strength and elastic modulus tended to increase as a whole. Specifically, the peak strength increased slowly at α=0°~75° and increased rapidly at α=75°~90°. The peak strength reached the maximum value of 126.95 MPa at α=90°, while the peak strength reached the minimum value of 86.35 MPa at α=15°. However, the elastic modulus growed in an S-shaped curve. The growth was gentle at α=0°~45° and α=75°~90° and rapid at α=45°~75°. When α=90°, the elastic modulus reached the maximum value of 18.84 GPa. When α=0° and α=15°, the elastic modulus was smaller, which was equal to 16.22 GPa and 16.83 GPa, respectively.Figure 8
Change curves of mechanical parameter of numerical models.
(a)
α
=
45
°
(b)
β
=
60
°In summary, the peak strength and elastic modulus decreased with the increase ofβ and increased with the increase of α on the whole. The reason was that as β decreased or α increased, the projection of the fissure length in the direction perpendicular to the loading direction gradually becomes shorter, leading to the gradual reduction of the axial closure effect of the fissure under the axial load, and the effective bearing area of the rock actually bearing the load gradually increased, which caused that the strength and the stiffness gradually increased [25].


## 3.2. Crack Initiation Stress and Damage Stress

The crack initiation stress and damage stress of rock are of great significance for understanding the progressive failure process of rock and establishing its brittle failure mechanism, initiation strength, and long-term strength criterion [31]. Through the method of monitoring microcracks, taking the intact numerical model as an example, Figure 9 shows the change curve of the microcrack number in the intact numerical model. The axial stress corresponding to the crack initiation point was taken as the initiation stress (σci), and the axial stress corresponding to the rapid crack growth point was taken as the damage stress (σcd) [32].Figure 9
Change curve of the microcrack number in intact numerical model.Table4 showed the initiation stress and damage stress of all numerical models determined by the above method. The following can be seen: (1) Compared with the intact numerical model, the crack initiation stress and damage stress of the numerical models containing X-shaped fissure were significantly reduced, but the level of crack initiation stress (the ratio of the initiation stress to the peak stress (σp)) and the level of damage stress (the ratio of damage stress to peak stress) did not change much. (2) The level of crack initiation stress in all numerical models was between 0.36 and 0.5, and the level of damage stress was between 0.76 and 0.84. Studies [33] have shown that the crack initiation stress of brittle rocks is about 0.25 to 0.5 times the peak stress, and the damage stress is about 0.7 to 0.85 times the peak stress.Table 4
The crack initiation stress and damage stress of numerical models.
α (°)
β (°)
σ
c
i (MPa)
σ
c
d (MPa)
σ
p (MPa)
σ
c
i
/
σ
p
σ
c
d
/
σ
p
Intact
74.2
153.8
183.0
0.41
0.84
45
15
44.5
90.3
110.3
0.40
0.82
45
30
43.0
82.7
104.0
0.41
0.80
45
45
36.8
74.4
91.9
0.40
0.81
45
60
45.7
76.5
95.1
0.48
0.80
45
75
44.0
73.4
96.1
0.46
0.76
45
90
43.8
70.8
87.2
0.50
0.81
0
60
37.5
74.4
89.5
0.42
0.83
15
60
34.7
71.3
86.4
0.40
0.83
30
60
34.1
79.5
95.0
0.36
0.84
45
60
45.7
76.5
95.1
0.48
0.80
60
60
49.3
83.6
102.0
0.48
0.82
75
60
43.4
84.2
104.1
0.42
0.81
90
60
55.2
102.5
127.0
0.44
0.81Figure10 shows the curves of the crack initiation stress and damage stress with α and β. From Figure 10(a), the following can be seen: (1) With the increase of β, both the crack initiation stress and the damage stress decreased first, then increased, and then decreased, which was consistent with the peak strength. When β=45°, both the crack initiation stress and the damage stress reached the minimum value. (2) With the increase of β, the level of crack initiation stress and the level of damage stress had no obvious changes, which were between 0.4 and 0.5 and between 0.76 and 0.82, respectively. When β=45°, the level of crack initiation stress reached the minimum value of 0.4, and at this time, cracks were most prone to initiation.Figure 10
Change curves of crack initiation stress and damage stress of numerical models containing X-shaped fissure.
(a)
α
=
45
°
(b)
β
=
60
°From Figure10(b), the following can be seen: (1) With the increase of α, the crack initiation stress and damage stress increased as a whole, which was consistent with the peak strength. (2) With the increase of α, the level of crack initiation stress and the level of damage stress had no obvious changes, which were between 0.36 and 0.48 and between 0.80 and 0.84, respectively. When α=30°, the level of crack initiation stress reached the minimum value of 0.36, and at this time, cracks were most prone to initiation.In comparison, with the changes ofα and β, the range of the level of damage stress was smaller than the level of crack initiation stress, indicating that α and β had a greater impact on the level of crack initiation stress than on the level of damage stress.


## 3.3. Law of Energy Evolution

The essence of rock failure is a state instability phenomenon driven by energy [34], and the study of energy from the perspective of mesomechanics is more in line with actual energy changes. In PFC simulation, the sum of particle bonding energy and strain energy is recorded as total strain energy, and the difference between boundary energy and total strain energy is recorded as dissipation energy [35]. Due to space limitations, the energy evolution law during the loading process was analyzed by taking the numerical model of α=45° and β=60° as an example, as shown in Figure 11, and other numerical models had similar laws.Figure 11
Relationship between energy and stress-strain of numerical model (α=45°, β=60°).From Figure11, the following can be seen: (1) At the initial stage of loading (from point O to point A, which was the crack initiation stress point), there were no microcracks, and the boundary energy and the total strain energy increased nonlinearly and basically coincide, indicating that the input boundary energy was all converted into total strain energy, and the dissipation energy, kinetic energy, and friction energy were almost zero. (2) At the linear elastic stage of loading (from point A to point B, which was the damage stress point), the microcracks initiated and growed steadily. The growth rate of the boundary energy and the total strain energy gradually increased, and the dissipated energy and the friction energy increased slowly. At point B, there was a small separation between the boundary energy and the total strain energy curve, and the strain energy stored in contact during the crack propagation process was gradually transformed into dissipated energy. (3) At the plastic yielding stage of loading (from point B to point C, which was the peak stress point), the microcracks expand unsteadily, and the boundary energy and the total strain energy curve gradually separated. The total strain energy reached the peak at point C, and the growth of dissipation energy gradually increased. (4) After the peak stress (from point C to point D), the boundary energy continued to increase, the stored total strain energy was rapidly released, and the dissipation energy increased sharply until the model was broken. During the entire deformation and failure process, the boundary energy keeps increasing, the total strain energy first increased and then decreased sharply, and the dissipation energy increased slowly and then increased sharply.Table5 shows the energy index of all numerical models at the peak stress. From Table 5, the following can be seen: (1) The total strain energy of all numerical models accounted for a large proportion of the boundary energy at the peak point, indicating that the input energy before the peak stress was mainly stored in the rock in the form of strain energy. (2) Compared with the intact numerical model, the boundary energy, total strain energy, and dissipation energy of the numerical models containing X-shaped fissure were significantly reduced, and the largest decline reached 71.7%, 72.5%, and 71.1%, respectively. It showed that the X-shaped fissure had obvious energy weakening effect on the rock. The total strain energy at the peak stress represents the energy storage limit of the rock and reflects the ability of the rock to resist damage, so the rock containing fissure was more likely to be damaged.Table 5
Energy indexes in the peak stress of numerical models.
α (°)
β (°)
Boundary energy (J)
Total strain energy (J)
Dissipated energy (J)
Intact
1070.8
1035.0
35.7
45
15
439.8
418.9
20.9
45
30
402.1
386.7
15.4
45
45
325.9
308.2
17.7
45
60
356.2
338.3
17.9
45
75
383.5
360.5
23.0
45
90
303.4
284.6
18.8
0
60
371.4
332.6
38.7
15
60
332.8
301.8
30.9
30
60
364.3
342.0
22.3
60
60
391.9
373.6
18.2
75
60
373.3
363.0
10.3
90
60
605.5
567.9
37.6Figure12 shows the curves of the energy indexes at peak stress with α and β. The following can be seen: (1) With the increase of β, the boundary energy and total strain energy showed a decreasing trend as a whole. Specifically, it first decreased, then increased and then decreased, which was consistent with the changing law of peak stress in Figure 8(a). (2) With the increase of α, the boundary energy and the total strain energy showed an overall increasing trend, which was consistent with the change law of the peak stress in Figure 8(b).Figure 12
Change curves of energy indexes of numerical models containing X-shaped fissure in the peak stress.
(a)
α
=
45
°
(b)
β
=
60
°


## 3.4. Crack Propagation Process and Failure Mode

Due to space limitations, the crack propagation process during the loading process was analyzed by taking the numerical model ofα=45° and β=60° as an example, as shown in Figure 13. Figures 14 and 15, respectively, show the failure modes of numerical models with different α and β. The red and black in Figures 13 and 14 represented tensile microcracks and shear microcracks, respectively, and the marked numbers indicated the sequence of macrocracks.Figure 13
Crack propagation process during the loading process of the numerical model containing X-shaped fissure (α=45°, β=60°). (Red represents tensile microcracks. Black represents shear microcracks. Numbers indicate the order of macrocrack generation. The angle value indicates the crack initiation angle).
(a)
σ
=
56.3
MPa
,
ε
=
3.3
×
10
−
3
(b)
σ
=
75.8
MPa
,
ε
=
4.4
×
10
−
3
(c)
σ
=
84.8
MPa
,
ε
=
5.2
×
10
−
3
(d)
σ
=
94.5
MPa
,
ε
=
5.6
×
10
−
3
(e)
σ
=
88.4
MPa
,
ε
=
5.8
×
10
−
3
(f)
σ
=
66.6
MPa
,
ε
=
6.0
×
10
−
3Figure 14
Failure modes of numerical models with different values ofβ (α=45°). (Red represents tensile microcracks. Black represents shear microcracks. Numbers indicate the order of macrocrack generation. The angle values indicate the crack initiation angle).
(a)
β
=
15
°
(b)
β
=
30
°
(c)
β
=
45
°
(d)
β
=
60
°
(e)
β
=
75
°
(f)
β
=
90
°Figure 15
Failure modes of numerical models with different values ofα (β=60°). (Red represents tensile microcracks. Black represents shear microcracks. Numbers indicate the order of macrocrack generation. The angle values indicate the crack initiation angle).
(a)
α
=
0
°
(b)
α
=
15
°
(c)
α
=
30
°
(d)
α
=
45
°
(e)
α
=
60
°
(f)
α
=
75
°
(g)
α
=
90
°From Figure13, it can be seen that when the load reached a certain stress level, the model first initiated wing cracks 1a and 1b at the tip of fissure ② along the extension direction of fissure ② at an inclination angle of 30°, which was the crack initiation angle [36]. Continuing to load, the wing cracks 1a and 1b propagated to the upper and lower ends of the model along the loading direction, and macrocracks 2, 3, and 4 were generated at the tip of fissure ①, which eventually led to the overall failure of the model.From Figure14, the following can be seen: (1) Wing cracks 1a and 1b were generated at the two tips of fissure ②. When β=15°, the crack initiation angle was 90°, and when β=90°, the crack initiation angle was 0°. Moreover, with the increase of β, the crack initiation angle decreased. (2) Macrocrack 2 was produced at the two tips of fissure ①. When β=15° and 60°, macrocrack 2 was produced at the left tip of fissure ①, and when β=30°, 45°, 75°, and 90°, macrocracks were generated at the right tip of fissure ①. (3) Macrocracks 3 and 4 were both produced at the two tips of fissure ① or at the ends of existing macrocracks. (4) When β=15°, 30°, 60°, and 75°, the four tips of the X-shaped fissure all had macrocracks and the distribution of microcracks was relatively dispersed, while when β=45° and 90°, only three tips had macrocracks and microcracks were concentrated at the fissure tips. Since the stress concentration of the former was not as high as that of the latter, the latter will fail at a lower stress level, which was also the reason why the peak strength of the model was small at β=45° and 90°.From Figure15, the following can be seen: (1) Except for α=0° and 90°, the numerical models first initiated wing cracks 1a and 1b at the two tips of fissure ②. When α=0°, the crack initiation angle was 90°; when α=90°, the crack initiation angle was 67°; and when α=60°, the crack initiation reached the minimum value of 0°. Moreover, with the increase of α, the crack initiation angle first decreased and then increased. (2) Macrocracks 2, 3, and 4 were all generated at the tips of fissure ① and fissure ② or at the end of the macrocracks that had been generated. (3) When α=0° and 15°, cracks were formed at the tips of fissure ① and fissure ②, so the peak strength was small. When α=90°, the four tips of the X-shaped fissure all had macrocracks and the distribution of microcracks was relatively scattered, so the peak strength was large.From Figures14 and 15, the following can be seen: (1) Tensile microcracks were the main cracks in the numerical models, and shear microcracks were almost rare. However, there were no microcracks within a certain range of the X-shaped fissure center. (2) The failure mode was mainly split along the axial direction, and the failure surface started from the tips of the fissure and extended to both ends of the numerical models.In summary, with the changes ofα and β, the numerical models first initiated wing cracks 1a and 1b at the two tips of fissure ② with a larger dip, and the crack initiation angle decreased with the increase of the inclination angle of fissure ②. From the theory of crack initiation, the relationship between the crack initiation angle of θ0 of a single-opening fissure tip and the inclination angle of α under uniaxial compression is as follows [36]:
(3)θ0=2arctan1+1+8sin2α/1+cos2α24sin2α/1+cos2α.The above formula shows that the crack initiation angle ofθ0 decreases with the increase of the inclination angle of α, which was consistent with the conclusion of this paper.


## 4. Heterogeneity Analysis

In order to evaluate the influence of the heterogeneity of the mechanical parameters on the macroscopic mechanical properties and failure modes of the rock containing X-shaped fissures under uniaxial compression, the numerical model ofα=45°,β=60° was taken as an example, and the homogeneity factor of m was 1, 2, 3, 5, 10, 20, 30, and 40, respectively. Figure 16 shows the relationship between the stress-strain curves and the macroscopic mechanical parameters of the numerical models under different homogeneity factors of m.Figure 16
Relationship between macromechanical properties and homogeneity factor ofm.
(a)
The stress-strain curves
(b)
The uniaxial compressive strength and elastic modulusFigure16 shows the following: (1) When the homogeneity factor of m increased from 1 to 40, the uniaxial compressive strength and elastic modulus increased from 55.6 MPa and 15.7 GPa to 103.3 MPa and 18.0 GPa, respectively, and increased by 85.8% and 14.6%, respectively. It can be seen that the heterogeneity of mechanical parameters had an impact on the uniaxial compressive strength and elastic modulus, and the influence on the uniaxial compressive strength was obviously greater than the elastic modulus. The reason for this phenomenon was that when the homogeneity factor was large, the strength of the mesoelement was relatively uniform, all the elements shared the external load, and the bearing capacity was large [20]. (2) When the homogeneity factor of m increased from 1 to 10, the uniaxial compressive strength and elastic modulus increased rapidly; while when the homogeneity factor of m increased from 10 to 40, the uniaxial compressive strength and elastic modulus increased slowly. (3) The relationship of uniaxial compressive strength and elastic modulus with the homogeneity factor of m conformed to the exponential function relationship (R2 was 0.99 and 0.80, respectively). These conclusions were consistent with the numerical simulation results of other scholars [37].Figure17 shows the failure modes of numerical models containing X-shaped fissure under different homogeneity factors of m. The red and black in the figure represented tensile microcracks and shear microcracks, respectively. From Figure 17, it can be seen that the heterogeneity of mechanical parameters had an obvious impact on the failure mode. When the homogeneity factor of m was small, the distribution of microcracks was relatively scattered, while when the homogeneity factor of m was large, the microcracks were mainly distributed in the tip of the X-shaped fissures. These conclusions were consistent with the numerical simulation results of other scholars [29].Figure 17
Failure modes of numerical models containing X-shaped fissure with different homogeneity factor values ofm. (Red represents tensile microcracks. Black represents shear microcracks).
(a)
m
=
1
(b)
m
=
2
(c)
m
=
3
(d)
m
=
5
(e)
m
=
10
(f)
m
=
20
(g)
m
=
30
(h)
m
=
40


## 5. Conclusion

In this paper, the particle flow program PFC2D was used to study the influence of X-shaped fissure and heterogeneity of mechanical parameters on the mechanical behavior of rock under uniaxial compression. The following conclusions were drawn:
(1)
Compared with intact rock, the X-shaped fissure not only weakened the brittleness of the rock but also deteriorated the mechanical properties. The deterioration of the peak strength was obviously greater than the elastic modulus. Moreover, the peak strength and elastic modulus generally decreased asβ increases and increased as α increases(2)
Compared with the intact numerical model, the initiation stress and damage stress of the models containing X-shaped fissure were significantly reduced, but the crack initiation stress level and damage stress level did not change much. Moreover, the change of the crack initiation stress and damage stress withα and β was consistent with peak strength(3)
Compared with the numerical intact model, the X-shaped fissure had obvious energy weakening effect on the rock, so the rock containing fissure was more likely to be damaged. The change of boundary energy and total strain energy withα and β at peak stress was consistent with peak strength(4)
The fracture was dominated by tensile microcracks, and no microcracks were generated in a certain range of the X-shaped fissure center. The failure mode was mainly split along the axial direction, and the failure surface started from the tips of the fissure and extended to both ends of the model(5)
Regardless of the changes ofα and β, the models all firstly initiated wing cracks 1a and 1b at the two tips of fissure ② with a larger dip, and the crack initiation angle decreased with the increase of the inclination angle of fissure ②(6)
The uniaxial compressive strength and elastic modulus increased exponentially with the increase of the homogeneity factor ofm. When the homogeneity factor of m was small, the microcracks were more evenly distributed in the model, and when the homogeneity factor of m was large, the microcracks mainly concentrated at the tips of the fissure



---

*Source: 8828199-2020-10-17.xml*