# Comparative Study on Mineral-Scale Microcrack Propagation of Shale under Different Loading Methods

Authors: Jianyong Zhang; Zhendong Cui; Weige Han; Kai Si; Yan Zhao
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

Producing a sufficient volume of multiscale crack networks is key to enhancing recovery of shale gas. The formation of crack network largely depends on initiation and propagation of microcracks. To reveal the influence of different loading methods on the propagation of mineral-scale microcracks, this study used the Voronoi tessellation technique to establish a cohesive zone model of shale mineral distribution and applied six different boundary conditions to represent different loading methods. Crack path characteristics, rupture characteristics, continuous crack propagation and turning, and en echelon intermittent crack propagation under different loading methods were compared and analyzed. The essence of different loading methods affecting the length and complexity of cracks was the spreading range of tensile microcracks. The mechanical properties of minerals led to dissimilarities in continuous crack propagation and turning. The formation and propagation of en echelon intermittent fractures of different scales were mainly impacted by the heterogeneity of minerals and mineral aggregates. The spreading direction and connection form of en echelon intermittent fractures were mainly affected by the loading method. Conclusions arising from mineral-scale simulations contribute to understanding the mechanism of microcrack propagation resulting from different loading methods, and these conclusions have a guiding significance to enhanced shale gas recovery.


---

## Body

## 1. Introduction

Shale gas exists in low-porosity and low-permeability shale formations. Producing a sufficient volume of multiscale crack networks is the key to enhancing the recovery of shale gas. At present, academia and industry are more concerned about the macrocrack network, and less research is focused on the propagation of microcracks. However, existing studies show that the formation of crack networks largely depends on the growth of microcracks. Minerals are the basic units of shale. Research on the physical nature of mineral-scale microcracks can reveal their propagation mechanism.Scholars have carried out static experimental observations on the characteristics of mineral-scale microcracks after the failure of rock through the use of optical microscopes, scanning electron microscopes (SEM), computed tomography (CT), and other types of equipment. Fujii et al. [1] found that more than 90% of the tensile cracks in granite after tensile failure exist in the form of intergranular fractures, and the remaining 10% are transgranular fractures. Minerals that form microcracks at different angles to the tensile direction are unusual. Zhong et al. [2] summarized microcrack failure modes of shale in a triaxial compression experiment and categorized them into three types: tensile failure, shear failure, and slip failure. Daigle et al. [3] compared the differences between microcracks in siliceous clay and calcareous clay. Cheng and Wong [4] compared the characteristics of tensile microcracks and shear microcracks produced by the uniaxial compression of marble and analyzed the influence of mineral grains on the types of microcracks. Lan et al. [5] carried out SEM observation of microcracks before and after a triaxial compression test of shale and found that the spatial distribution of microcracks conforms to the power-law distribution, with a defined self-affine and hierarchical structure. This study suggests that the mechanical responses of minerals have different effect on shale fracturing.The above experiments were all static observation and analysis and did not evaluate the growth process of microcracks. SEM with an in situ loading apparatus makes it possible to observe continuous growth of microcracks. Zhao et al. [6] found in a uniaxial compression experiment of Fangshan marble that along with initiation, propagation, and connection of microcracks, the closure of microcracks also occur. Cui et al. [7, 8] observed the dynamic process of microcrack initiation and propagation and discovered microscopic en echelon intermittent fractures for the first time. They summarized the hierarchical distribution characteristics of en echelon intermittent fractures and hypothesized that the cause may be related to the heterogeneity of minerals. Zuo et al. [9] observed the continuous propagation process of main cracks and branch cracks in the microscopic region and performed finite element calculations on their fracture toughness. Dong et al. [10] found that cracks propagate in the form of microcrack connections at the microscopic scale. In other words, the microcracks near the crack tip are always in a state of competition. Microcracks connected to the crack tip become the main crack. After the main crack propagates, microcracks that are not connected to the crack tip are closed owing to stress release. Tang [11] studied the evolution of microdamage in marble under uniaxial compression and proposed a multiscale damage power-law model for cracks. Renard et al. [12] performed CT scans on the entire process of monzonite triaxial compression under 25 MPa confining pressure and used digital image correlation methods to quantitatively analyze the initiation, closure, aggregation, and penetration characteristics of microcracks at different loading stages. These observations of the continuous growth of microcracks found many interesting phenomena that are different from the growth of macroscopic cracks. Unfortunately, none of the observations involved the distribution of minerals in the microregion, so it is impossible to clarify the influence of minerals on the growth of microcracks and thus only hypothetical explanations can be given regarding the mechanism of different observed phenomena.The in situ experimental observation of the continuous growth of microcracks based on mineral distribution is difficult. The strain and stress of the microregion also cannot be obtained easily. Scholars have developed many numerical models based on mineral distribution, such as the universal distinct element code (UDEC) particle boundary model [13, 14], three-dimensional distinct element code (3DEC)-Voronoi model [15, 16], distinct element model (DEM) based on three-dimensional polygons [17], particle flow code (PFC)- grain-based model (GBM) [15, 18, 19], and the Irazu-GBM model [20]. Li et al. [21, 22] used the finite discrete element method based on grains to study factors, including boundary conditions and rock geometry, affecting the fracture growth process in granite under uniaxial compression. Lan et al. [23] carried out a uniaxial compression experiment of a mineral distribution model based on UDEC. The authors found that heterogeneity of grain size affected the distribution of tensile stress, and rocks with such grains are more likely to produce tensile microcracks than rocks with uniform grains. Li et al. [17] studied the influence of grain size on microcrack propagation based on the 3D polycrystalline discrete element method-Voronoi model. The authors reported that in specimens with larger grain sizes, the convenient path for crack propagation is along the loading direction, which results in a series of vertical tensile fractures. When the grain size was small, the microcracks were more scattered. Xu et al. [24] established a uniaxial compression experimental model based on the microstructure of coal. They found that microcracks at the mineral interfaces occur due to heterogeneous strain of minerals. The authors concluded that the generation of tensile microcracks dominated the failure behavior of the rock under uniaxial compression. Zhou et al. [25] established PFC-GBM to study the proportion of shear and tensile cracks during crack propagation. Saadat and Taheri [26] carried out uniaxial compression experiments based on the cohesive contact model of PFC to study the formation and propagation of microcracks in the shear zone.The abovementioned mathematical simulations of microcrack propagation based on mineral distribution models were all set as uniaxial compression tests. The results suggest that tensile fractures dominated the uniaxial compression failure behavior of rock. However, hydraulic fracturing is not a simple uniaxial compression. In hydraulic fracturing, the fractures are first partially opened, and then the proppant enters to propagate the fractures. At the same time, the high pressure of the fracturing fluid pushes both sides of the fractures to move. In this process, the shale is in a complex stress state of tension, shear, and compression [27]. Observations of hydraulic fracturing cracks show that the characteristics of microcracks produced by tensile, shear, tension-shear, and compression-shear are different [28–30]. Other experimental studies focused on tensile and three-point bending loading, and the observed crack network and microcrack propagation were found to be quite different from uniaxial compression [9, 31]. Therefore, it has important significance to study the rupture characteristics of microcracks under different loading methods. On the macroscale, some scholars have compared the peak stress, macrocrack propagation, and fracture characteristics of rocks under different loading methods [32–34]. However, on the mineral scale, it is essential to determine what characteristics of the microcracks are produced under different loading methods such as tension, shear, tension-shear, and compression-shear. It is also essential to determine the difference between microcrack growth and evolutionary form under different loading methods and what types of microcracks dominates the fracture of the rock. These important questions need to be addressed.To reveal the influence of different loading methods on the propagation of mineral-scale microcracks, this study used the Voronoi tessellation technique to establish a cohesive zone model (CZM) based on mineral distribution. Considering the heterogeneity of minerals, various mechanical parameters were set for different minerals. Six loading methods were designed, that is, modeling tension, tension-shear, shear, and compression-shear stresses. We compared and analyzed the crack length, fractal dimension, rupture characteristics, continuous crack propagation and turning, and en echelon intermittent fracture propagation. The comparative analysis uncovered the mechanism of how different loading methods affect the growth of mineral-scale microcracks.


## 2. Cohesive Zone Model Based on Mineral Distribution

### 2.1. Voronoi Tessellation Technique

To establish a heterogeneous shale mineral model, we added a Voronoi tessellation technique into Abaqus using a Python script. This technique contains four steps as shown in Figure1. First, a certain number of randomly distributed control points were generated in the model sketch (Figure 1(a)). Then, the adjacent control points were triangulated based on the Delaunay triangulation (Figure 1(b)). When the circumference of the constructed triangle included control points in addition to those at the vertices, the triangle was abandoned, and new triangles were formed until all control points formed vertices. Third, the circumcenters of adjacent triangles were connected (Figure 1(c)). Finally, the control points and triangles were removed and the Voronoi grains were generated (Figure 1(d)).Figure 1
Four steps to generate random Voronoi polygons [35].
(a)(b)(c)(d)


### 2.2. Numerical Calculation Principle of the Cohesion Zone Model

CZM was proposed by Dugdale [36] and Barenblatt [37]. This model helped in solving the research difficulties caused by the crack tip stress singularity and is widely used in the study of interface delamination of composite materials [38]. Strom and Parmigiani [39] used CZM to explain the penetration and deflection behavior of cracks when passing through the bedding surface based on the stress-energy method and verified the applicability of the model.


#### 2.2.1. Constitutive Equation of Cohesive Element

The cohesive element is based on the linear elastic constitutive model of the traction-separation law. Before the cohesive element is damaged, the stress and strain satisfy the following linear elastic relationship:(1)tntstt=KnnKnsKntKnsKssKstKntKstKttεnεsεt.The symbolt represents stress, ε represents strain, the subscript n represents the normal direction, and the subscripts s and t are two tangential directions perpendicular to the normal direction. For uncoupled constitutive relations, only Knn, Kss, and Ktt need to be defined.


#### 2.2.2. Fracture Propagation Criteria

Initial Damage Criterion. The maximum principal stress control criterion was used in the simulation. That is, when the normal tensile stress or tangential stress reaches the corresponding strength, it will break:(2)MAXσnNmax,σsSmax,σtTmax=1,where σn is the normal stress; σs and σt are the shear stresses in two directions; Nmax is the tensile strength; and Smax and Tmax are the shear strengths in the two directions; σn is expressed as(3)σn=0,σn<0,σn=σn,σn≥0.The cohesive element will not produce initial damage in a purely compressed state.Damage Evolution Criterion. In this study, the damage evolution criterion was based on effective displacement, and the damage variable D was defined as(4)D=δmfδmmax−δm0δmmaxδmf−δm0,where δmf is the effective displacement at failure, taken as 0.01 in the model, δm0 is the effective displacement at the initial evolution of damage, δmmax, and is the maximum effective displacement during the loading process.


#### 2.2.3. Principle of Acoustic Emission Simulation

The Abaqus simulation results were extracted using a Python script, and the results were processed by MATLAB programming to simulate the acoustic emission (AE) during the loading process [40]. Through the AE location and characteristic parameters, the crack rupture condition was further described, which provided an effective means for in-depth analysis of the crack propagation process. The method to obtain the AE location can be described as follows. The node coordinates of the damaged element were extracted and the coordinates of the center point of the damaged element were calculated, which were then used as the point of the AE event location. The number of damaged elements was regarded as the number of AE events. Extracting the total energy dissipated per volume of the damaged element gave the AE energy.The parameter MMIXDME, which can be used to determine the type of rupture, is defined as follows:(5)MMIXDME=Gs+GtGT,GT=Gn+Gs+Gt,where Gn is the Type I tensile fracture energy; Gs is the Type II slip fracture energy; Gt is the Type III tear fracture energy.; and GT is the sum of Gn, Gs, and Gt.The value range of MMIXDME is between 0 and 1. When MMIXDME is 0, it indicates a tensile fracture; when the value is 1, it indicates a shear fracture, and when the value is between the 0 and 1, there is a tensile-shear mixed fracture.


## 2.1. Voronoi Tessellation Technique

To establish a heterogeneous shale mineral model, we added a Voronoi tessellation technique into Abaqus using a Python script. This technique contains four steps as shown in Figure1. First, a certain number of randomly distributed control points were generated in the model sketch (Figure 1(a)). Then, the adjacent control points were triangulated based on the Delaunay triangulation (Figure 1(b)). When the circumference of the constructed triangle included control points in addition to those at the vertices, the triangle was abandoned, and new triangles were formed until all control points formed vertices. Third, the circumcenters of adjacent triangles were connected (Figure 1(c)). Finally, the control points and triangles were removed and the Voronoi grains were generated (Figure 1(d)).Figure 1
Four steps to generate random Voronoi polygons [35].
(a)(b)(c)(d)


## 2.2. Numerical Calculation Principle of the Cohesion Zone Model

CZM was proposed by Dugdale [36] and Barenblatt [37]. This model helped in solving the research difficulties caused by the crack tip stress singularity and is widely used in the study of interface delamination of composite materials [38]. Strom and Parmigiani [39] used CZM to explain the penetration and deflection behavior of cracks when passing through the bedding surface based on the stress-energy method and verified the applicability of the model.


### 2.2.1. Constitutive Equation of Cohesive Element

The cohesive element is based on the linear elastic constitutive model of the traction-separation law. Before the cohesive element is damaged, the stress and strain satisfy the following linear elastic relationship:(1)tntstt=KnnKnsKntKnsKssKstKntKstKttεnεsεt.The symbolt represents stress, ε represents strain, the subscript n represents the normal direction, and the subscripts s and t are two tangential directions perpendicular to the normal direction. For uncoupled constitutive relations, only Knn, Kss, and Ktt need to be defined.


### 2.2.2. Fracture Propagation Criteria

Initial Damage Criterion. The maximum principal stress control criterion was used in the simulation. That is, when the normal tensile stress or tangential stress reaches the corresponding strength, it will break:(2)MAXσnNmax,σsSmax,σtTmax=1,where σn is the normal stress; σs and σt are the shear stresses in two directions; Nmax is the tensile strength; and Smax and Tmax are the shear strengths in the two directions; σn is expressed as(3)σn=0,σn<0,σn=σn,σn≥0.The cohesive element will not produce initial damage in a purely compressed state.Damage Evolution Criterion. In this study, the damage evolution criterion was based on effective displacement, and the damage variable D was defined as(4)D=δmfδmmax−δm0δmmaxδmf−δm0,where δmf is the effective displacement at failure, taken as 0.01 in the model, δm0 is the effective displacement at the initial evolution of damage, δmmax, and is the maximum effective displacement during the loading process.


### 2.2.3. Principle of Acoustic Emission Simulation

The Abaqus simulation results were extracted using a Python script, and the results were processed by MATLAB programming to simulate the acoustic emission (AE) during the loading process [40]. Through the AE location and characteristic parameters, the crack rupture condition was further described, which provided an effective means for in-depth analysis of the crack propagation process. The method to obtain the AE location can be described as follows. The node coordinates of the damaged element were extracted and the coordinates of the center point of the damaged element were calculated, which were then used as the point of the AE event location. The number of damaged elements was regarded as the number of AE events. Extracting the total energy dissipated per volume of the damaged element gave the AE energy.The parameter MMIXDME, which can be used to determine the type of rupture, is defined as follows:(5)MMIXDME=Gs+GtGT,GT=Gn+Gs+Gt,where Gn is the Type I tensile fracture energy; Gs is the Type II slip fracture energy; Gt is the Type III tear fracture energy.; and GT is the sum of Gn, Gs, and Gt.The value range of MMIXDME is between 0 and 1. When MMIXDME is 0, it indicates a tensile fracture; when the value is 1, it indicates a shear fracture, and when the value is between the 0 and 1, there is a tensile-shear mixed fracture.


## 2.2.1. Constitutive Equation of Cohesive Element

The cohesive element is based on the linear elastic constitutive model of the traction-separation law. Before the cohesive element is damaged, the stress and strain satisfy the following linear elastic relationship:(1)tntstt=KnnKnsKntKnsKssKstKntKstKttεnεsεt.The symbolt represents stress, ε represents strain, the subscript n represents the normal direction, and the subscripts s and t are two tangential directions perpendicular to the normal direction. For uncoupled constitutive relations, only Knn, Kss, and Ktt need to be defined.


## 2.2.2. Fracture Propagation Criteria

Initial Damage Criterion. The maximum principal stress control criterion was used in the simulation. That is, when the normal tensile stress or tangential stress reaches the corresponding strength, it will break:(2)MAXσnNmax,σsSmax,σtTmax=1,where σn is the normal stress; σs and σt are the shear stresses in two directions; Nmax is the tensile strength; and Smax and Tmax are the shear strengths in the two directions; σn is expressed as(3)σn=0,σn<0,σn=σn,σn≥0.The cohesive element will not produce initial damage in a purely compressed state.Damage Evolution Criterion. In this study, the damage evolution criterion was based on effective displacement, and the damage variable D was defined as(4)D=δmfδmmax−δm0δmmaxδmf−δm0,where δmf is the effective displacement at failure, taken as 0.01 in the model, δm0 is the effective displacement at the initial evolution of damage, δmmax, and is the maximum effective displacement during the loading process.


## 2.2.3. Principle of Acoustic Emission Simulation

The Abaqus simulation results were extracted using a Python script, and the results were processed by MATLAB programming to simulate the acoustic emission (AE) during the loading process [40]. Through the AE location and characteristic parameters, the crack rupture condition was further described, which provided an effective means for in-depth analysis of the crack propagation process. The method to obtain the AE location can be described as follows. The node coordinates of the damaged element were extracted and the coordinates of the center point of the damaged element were calculated, which were then used as the point of the AE event location. The number of damaged elements was regarded as the number of AE events. Extracting the total energy dissipated per volume of the damaged element gave the AE energy.The parameter MMIXDME, which can be used to determine the type of rupture, is defined as follows:(5)MMIXDME=Gs+GtGT,GT=Gn+Gs+Gt,where Gn is the Type I tensile fracture energy; Gs is the Type II slip fracture energy; Gt is the Type III tear fracture energy.; and GT is the sum of Gn, Gs, and Gt.The value range of MMIXDME is between 0 and 1. When MMIXDME is 0, it indicates a tensile fracture; when the value is 1, it indicates a shear fracture, and when the value is between the 0 and 1, there is a tensile-shear mixed fracture.


## 3. Mathematical Model

### 3.1. Model Design

A 2D geometric model (6 mm × 6 mm) was established in Abaqus, as shown in Figure2(a). To create a heterogeneous mineral model of shale, this 2D geometric model was divided into 722 Voronoi units using the Voronoi tessellation technique. Each Voronoi unit represented a mineral grain, and the grain size was assigned according to the study of Ji et al. [41]. Quartz, feldspar, kaolinite, and illite were randomly assigned to Voronoi units for a content of 30%, 30%, 20%, and 20%, respectively, and are represented by gray scales from high to low in Figure 2(b). Finally, the generated Voronoi units were meshed with a 0.08-mm quadrilateral grid. Zero-thickness cohesive elements were globally inserted into grain boundaries and meshes. Each mineral set contains four-node plane strain elements (CPE4) and zero-thickness cohesive elements at the boundaries of the CPE4. The mineral boundary set was created with the remaining cohesive elements at dissimilar mineral boundaries. The parameters of the CPE4 for different minerals are listed in Table 1. The parameters of cohesive elements in minerals (at boundaries of CPE4) and cohesive elements at boundaries of dissimilar minerals are shown in Table 2. It should be noted that mechanical properties of mineral boundaries depend on bond strength, and there is currently no uniform method for characterizing the bond strength of solid clay. Many experimental studies found that compared with clay minerals, mineral boundaries have a higher bulk density and stronger cementation [42–45]. Therefore, the mineral boundary strength parameters were set slightly larger than those of kaolinite but smaller than those of feldspar in this model.Figure 2
(a) 2D geometric model with its dimensions. The notch was 0.3 mm wide and 1.15 mm high. For the ease of reference, the left boundary of the model is marked as L, the lower left boundary is marked as LD, the right boundary is marked as R, and the upper right boundary is marked as RU. (b) Heterogeneous rock mineral model, quartz, feldspar, kaolinite, and illite are represented by gray scales from high to low.Table 1
Mineral composition and parameters of CPE4 of different minerals.
MineralPercentageElastic modulus (GPa)Poisson’s ratioQuartz30400.21Feldspar30350.23Kaolinite20220.25Illite20180.24Table 2
Parameters of cohesive elements in minerals (at boundaries of CPE4) and cohesive elements at boundaries of dissimilar minerals.
Cohesive elementKnn (MPa·m−1)Kss (MPa·m−1)Ktt (MPa·m−1)Tensile strength (MPa)Shear strength (MPa)In quartz120120120612In feldspar110110110510In kaolinite90909024In illite80808012At the boundaries of dissimilar minerals10010010036The total number of elements in the model was 21,242, of which 6,830 were CPE4 and 14,412 were cohesive elements. We set the total loading time to 1 s, the minimum increment size to 1E-7 s, and the maximum increment size to 0.1 s.


### 3.2. Loading Methods

Six loading methods were designed as shown in Figure3, which accounted for possible stress states of hydraulic fracturing, such as tension, shear, tension-shear, and compression-shear. The mechanical boundary conditions of the model are listed in Table 3.Figure 3
Six different loading methods. (a) Loading method I: a displacement load ofX = -0.2 mm is applied to the left boundary L and the right boundary R is fixed. (b) Loading method II: a displacement load of X = −0.2 mm is applied to the left boundary L; a displacement load of Y = 0.2 mm is applied to the lower left boundary LD, and the right boundary R is fixed. (c) Loading method III: the left boundary L is applied with X = 0.2 mm, Y = 0.2 mm displacement load, the left lower boundary LD is applied with Y = 0.2 mm displacement load, and the right boundary R is fixed. (d) Loading method IV: Y = 0.2 mm displacement load is applied to the lower left boundary LD, and the right boundary R and the upper right boundary RU are fixed. (e) Loading method V: a pressure load of X = 10 MPa is applied to the left boundary L; a displacement load of Y = 0.2 mm is applied to the lower left boundary LD, and the right boundary R and upper right boundary RU are fixed. (f) Loading method VI: the left boundary L is applied with a pressure load of X = 20 MPa, the left lower boundary LD is applied with a displacement load of Y = 0.2 mm, and the right boundary R and the upper right boundary RU are fixed.
(a)(b)(c)(d)(e)(f)Table 3
Mechanical boundary conditions.
IIIIIIIVVVILX = −0.2 mmX = −0.2 mm(−0.2 mm, 0.2 mm)—X = 10 MPaX = 20 MPaLD—Y = 0.2 mmY = 0.2 mmY = 0.2 mmY = 0.2 mmY = 0.2 mmR(0,0)(0,0)(0,0)(0,0)(0,0)(0,0)RU———(0,0)(0,0)(0,0)


## 3.1. Model Design

A 2D geometric model (6 mm × 6 mm) was established in Abaqus, as shown in Figure2(a). To create a heterogeneous mineral model of shale, this 2D geometric model was divided into 722 Voronoi units using the Voronoi tessellation technique. Each Voronoi unit represented a mineral grain, and the grain size was assigned according to the study of Ji et al. [41]. Quartz, feldspar, kaolinite, and illite were randomly assigned to Voronoi units for a content of 30%, 30%, 20%, and 20%, respectively, and are represented by gray scales from high to low in Figure 2(b). Finally, the generated Voronoi units were meshed with a 0.08-mm quadrilateral grid. Zero-thickness cohesive elements were globally inserted into grain boundaries and meshes. Each mineral set contains four-node plane strain elements (CPE4) and zero-thickness cohesive elements at the boundaries of the CPE4. The mineral boundary set was created with the remaining cohesive elements at dissimilar mineral boundaries. The parameters of the CPE4 for different minerals are listed in Table 1. The parameters of cohesive elements in minerals (at boundaries of CPE4) and cohesive elements at boundaries of dissimilar minerals are shown in Table 2. It should be noted that mechanical properties of mineral boundaries depend on bond strength, and there is currently no uniform method for characterizing the bond strength of solid clay. Many experimental studies found that compared with clay minerals, mineral boundaries have a higher bulk density and stronger cementation [42–45]. Therefore, the mineral boundary strength parameters were set slightly larger than those of kaolinite but smaller than those of feldspar in this model.Figure 2
(a) 2D geometric model with its dimensions. The notch was 0.3 mm wide and 1.15 mm high. For the ease of reference, the left boundary of the model is marked as L, the lower left boundary is marked as LD, the right boundary is marked as R, and the upper right boundary is marked as RU. (b) Heterogeneous rock mineral model, quartz, feldspar, kaolinite, and illite are represented by gray scales from high to low.Table 1
Mineral composition and parameters of CPE4 of different minerals.
MineralPercentageElastic modulus (GPa)Poisson’s ratioQuartz30400.21Feldspar30350.23Kaolinite20220.25Illite20180.24Table 2
Parameters of cohesive elements in minerals (at boundaries of CPE4) and cohesive elements at boundaries of dissimilar minerals.
Cohesive elementKnn (MPa·m−1)Kss (MPa·m−1)Ktt (MPa·m−1)Tensile strength (MPa)Shear strength (MPa)In quartz120120120612In feldspar110110110510In kaolinite90909024In illite80808012At the boundaries of dissimilar minerals10010010036The total number of elements in the model was 21,242, of which 6,830 were CPE4 and 14,412 were cohesive elements. We set the total loading time to 1 s, the minimum increment size to 1E-7 s, and the maximum increment size to 0.1 s.


## 3.2. Loading Methods

Six loading methods were designed as shown in Figure3, which accounted for possible stress states of hydraulic fracturing, such as tension, shear, tension-shear, and compression-shear. The mechanical boundary conditions of the model are listed in Table 3.Figure 3
Six different loading methods. (a) Loading method I: a displacement load ofX = -0.2 mm is applied to the left boundary L and the right boundary R is fixed. (b) Loading method II: a displacement load of X = −0.2 mm is applied to the left boundary L; a displacement load of Y = 0.2 mm is applied to the lower left boundary LD, and the right boundary R is fixed. (c) Loading method III: the left boundary L is applied with X = 0.2 mm, Y = 0.2 mm displacement load, the left lower boundary LD is applied with Y = 0.2 mm displacement load, and the right boundary R is fixed. (d) Loading method IV: Y = 0.2 mm displacement load is applied to the lower left boundary LD, and the right boundary R and the upper right boundary RU are fixed. (e) Loading method V: a pressure load of X = 10 MPa is applied to the left boundary L; a displacement load of Y = 0.2 mm is applied to the lower left boundary LD, and the right boundary R and upper right boundary RU are fixed. (f) Loading method VI: the left boundary L is applied with a pressure load of X = 20 MPa, the left lower boundary LD is applied with a displacement load of Y = 0.2 mm, and the right boundary R and the upper right boundary RU are fixed.
(a)(b)(c)(d)(e)(f)Table 3
Mechanical boundary conditions.
IIIIIIIVVVILX = −0.2 mmX = −0.2 mm(−0.2 mm, 0.2 mm)—X = 10 MPaX = 20 MPaLD—Y = 0.2 mmY = 0.2 mmY = 0.2 mmY = 0.2 mmY = 0.2 mmR(0,0)(0,0)(0,0)(0,0)(0,0)(0,0)RU———(0,0)(0,0)(0,0)


## 4. Results

This study successfully compared the mineral-scale microcrack path, AE, continuous crack propagation and turning, and en echelon intermittent crack propagation under different loading methods.


### 4.1. Crack Path

The longer and more complex the crack produced by fracturing, the more gas-containing holes are opened, and the greater the effect on the improvement of shale gas recovery. The crack length can be measured and the complexity of the crack path can be characterized by the fractal dimension. In this study, the box dimension method was used to calculate the fractal dimension of cracks under six different loading methods. Larger fractal dimensions indicated more complicated cracks.The principle of the box dimension method is as follows. Cover the cracked area with square boxes with side lengthr. Some entire boxes are empty and the rest of the boxes cover a part of the crack. Nr is the number of boxes covering the cracks; when r⟶0, the fractal dimension is as follows:(6)D=−limr⟶0logNrlogr.The geometric features of the fracture paths are shown in Figure4. The crack length exhibited the same trend as the fractal dimension. The crack length of tension and tension-shear was longer, and the fractal dimension was larger. The crack length of shear and compression-shear was smaller, and the fractal dimension was lower. In loading methods I–III, the crack length and fractal dimension gradually increased, whereas in loading methods IV–VI, the crack length and fractal dimension gradually decreased.Figure 4
Geometric features of fracture paths under six different loading methods. (a) Crack length; (b) Fractal dimension.
(a)(b)


### 4.2. Acoustic Emission (AE)

The AE varied for crack paths under each of the six loading methods. The type of fracture, AE counts, and AE energy were extracted, as shown in Figure5.Figure 5
AE features under six different loading methods. (a) AE counts; (b) AE energy.
(a)(b)In contrast to the parabolic change of crack length and fractal dimension, for loading methods I–VI, the AE counts, proportion of AE counts, and proportion of AE energy of tensile fractures decreased gradually. Meanwhile, AE counts and AE energy of shear fractures increased. The AE energy of the tensile fractures first increased and then decreased in general. In loading methods I–III, the total AE counts were higher than those in loading methods IV–VI, but the total AE energy for I-II was lower.


### 4.3. Continuous Crack Propagation and Turning

In homogeneous materials, the crack propagation path is controlled only by force; therefore, the crack grows relatively flat and smooth. However, in heterogeneous materials such as shale, the fracture path is often complicated. These different minerals and their boundaries affect continuous propagation and turning of cracks. Figure6 shows the final crack path (in red) under loading methods I–VI. The continuous growth and turning of cracks in different minerals and mineral boundaries are marked with arrows in different colors.Figure 6
Continuous growth and deflection of cracks (in red) under six different loading methods. The yellow arrows refer to crack growth in large scales of kaolinite and illite, the green arrows refer to crack growth in smaller scales of kaolinite and illite surrounded by quartz and feldspar, the blue arrows show crack growth in the boundary between different minerals, and the black arrows refer to crack growth in quartz and feldspar.
(a)(b)(c)(d)(e)(f)


#### 4.3.1. Propagation and Turning of Cracks in Kaolinite and Illite

Among the six loading methods, the cracks were relatively flat and smooth when they propagated in the larger grains of kaolinite and illite. The length of each segment of the crack through these regions was longer, and thus these regions accounted for a higher proportion of the total crack length. When the crack tip was blocked by quartz and feldspar, the crack turned and the turning angle was low.When small grains of kaolinite and illite were surrounded by quartz and feldspar, many small-scale inflections were produced inside the kaolinite and illite. In loading methods I–III, these inflections were relatively small. In loading methods IV–VI, there were more inflections with larger amplitudes.


#### 4.3.2. Propagation and Turning of Cracks at Mineral Boundaries

In loading methods I–III, only a small proportion of the cracks propagated along the mineral boundary. These segments usually propagated for a short distance along the mineral boundary and then turned into a kaolinite or illite grain. The positions of the cracks along the mineral boundary in loading methods IV–VI were similar. In loading method IV, cracks had more inflections, whereas in loading method V, crack inflections were reduced. In loading method VI, the cracks were smooth with almost no inflection, forming a shear arc.


#### 4.3.3. Propagation and Turning of Cracks in Quartz and Feldspar

The cracks rarely extended directly through quartz and feldspar. When it occurred at the end of the crack propagation process, the crack shape was relatively straight. However, when it occurred in the middle of the crack propagation process, there were more inflections.


### 4.4. En Echelon Intermittent Fractures

In addition to continuous propagation and turning, cracks were also generated, extended, and connected in a discontinuous manner, that is, en echelon intermittent cracks were created (Figure7).Figure 7
The initiation and evolution of en echelon intermittent fractures under six different loading methods. Every loading method is shown with five figures that demonstrate the growth of en echelon intermittent fractures over time.
(a)(b)(c)(d)(e)(f)The propagation process of en echelon intermittent fractures is divided into five stages.Stage 1: the small-scale en echelon intermittent fractures initiate and form the first large-scale en echelon intermittent fracture. In loading methods I, V, and VI, the en echelon intermittent fractures spread to the upper left direction. In loading methods II and III, the en echelon intermittent fractures spread to the upper right direction; whereas in loading methods IV, en echelon intermittent fractures were initiated both in the upper left and upper right directions.Stage 2: In loading method I, the second and third large-scale en echelon intermittent fractures were almost formed simultaneously. The direction of the second en echelon intermittent fracture was to the upper left, and the direction of the third en echelon intermittent fracture was to the upper right. The second large-scale en echelon intermittent fracture in loading methods II–V was located in the middle of the sample, whereas it was located at the top of the sample in loading method VI. The directions of the second large-scale en echelon intermittent fractures in loading methods II–V were all to the upper right. From loading method II to loading method V, the length of the second large-scale en echelon intermittent fractures gradually decreased.Stage 3: the third large-scale en echelon intermittent fractures in loading methods II, IV, V, and VI were all formed on the top of the sample. During this stage, the first and second large-scale en echelon intermittent fractures moved toward each other but did not connect. The three large-scale en echelon intermittent fractures in loading method I were also connected. In loading method III, the third en echelon intermittent fracture formed in the middle of the model and connected to the crack at the lower position.Stage 4: in loading methods I, IV, V, and VI, the two large-scale en echelon intermittent fractures at the middle and lower positions connected at this stage. However, in loading methods II and III, the two en echelon intermittent fractures in the middle and upper positions connected.Stage 5: In this stage, the remaining en echelon intermittent fractures connected and the rock was broken.The connection type of the en echelon intermittent fractures of scales under different loading methods was quite different and is summarized as follows:(1)
Wing Crack Connections. After the formation of en echelon intermittent fractures, adjacent en echelon intermittent fractures were connected by generating wing cracks at both ends, and these resembled an S-shape or a reverse S-shape. The wing crack connection was evident in the connection of two large-scale en echelon intermittent fractures and was also produced in the connection of small-scale en echelon intermittent fractures. Wing crack connections appeared in all six loading methods.(2)
Fragment Connection. In loading methods II–VI, there were en echelon intermittent fractures that appeared to connect in the form of fragments. In loading method III, fragments were generated in the small-scale en echelon intermittent fracture connection. The fragments in loading methods II, IV, V, and VI were generated during the large-scale en echelon intermittent fracture connection. The fragment positioning for loading methods IV–VI was the same. With the increase in lateral pressure, the size of the fragments gradually decreased.(3)
X-Shaped Crack Connection. For loading methods I and IV, the en echelon intermittent fractures appeared to create an X-shaped connection. The scale of X-shaped cracks in loading method I was larger than that in loading method IV.(4)
Not Connected or Not Connected at the Crack Tip. The large-scale en echelon intermittent fracture in the lower part of loading method II propagated to the upper right direction, but the crack tips were not connected with other cracks. Instead, fragments were generated in the middle of the crack that connected with the other two penetrating cracks, and the en echelon intermittent fracture that extended to the upper right was abandoned by the main crack. In loading method III, the en echelon intermittent fracture in the lower position directly propagated to the upper right and did not connect with the cracks at the upper position. Consequently, the cracks at the upper position were abandoned by the main crack.


## 4.1. Crack Path

The longer and more complex the crack produced by fracturing, the more gas-containing holes are opened, and the greater the effect on the improvement of shale gas recovery. The crack length can be measured and the complexity of the crack path can be characterized by the fractal dimension. In this study, the box dimension method was used to calculate the fractal dimension of cracks under six different loading methods. Larger fractal dimensions indicated more complicated cracks.The principle of the box dimension method is as follows. Cover the cracked area with square boxes with side lengthr. Some entire boxes are empty and the rest of the boxes cover a part of the crack. Nr is the number of boxes covering the cracks; when r⟶0, the fractal dimension is as follows:(6)D=−limr⟶0logNrlogr.The geometric features of the fracture paths are shown in Figure4. The crack length exhibited the same trend as the fractal dimension. The crack length of tension and tension-shear was longer, and the fractal dimension was larger. The crack length of shear and compression-shear was smaller, and the fractal dimension was lower. In loading methods I–III, the crack length and fractal dimension gradually increased, whereas in loading methods IV–VI, the crack length and fractal dimension gradually decreased.Figure 4
Geometric features of fracture paths under six different loading methods. (a) Crack length; (b) Fractal dimension.
(a)(b)


## 4.2. Acoustic Emission (AE)

The AE varied for crack paths under each of the six loading methods. The type of fracture, AE counts, and AE energy were extracted, as shown in Figure5.Figure 5
AE features under six different loading methods. (a) AE counts; (b) AE energy.
(a)(b)In contrast to the parabolic change of crack length and fractal dimension, for loading methods I–VI, the AE counts, proportion of AE counts, and proportion of AE energy of tensile fractures decreased gradually. Meanwhile, AE counts and AE energy of shear fractures increased. The AE energy of the tensile fractures first increased and then decreased in general. In loading methods I–III, the total AE counts were higher than those in loading methods IV–VI, but the total AE energy for I-II was lower.


## 4.3. Continuous Crack Propagation and Turning

In homogeneous materials, the crack propagation path is controlled only by force; therefore, the crack grows relatively flat and smooth. However, in heterogeneous materials such as shale, the fracture path is often complicated. These different minerals and their boundaries affect continuous propagation and turning of cracks. Figure6 shows the final crack path (in red) under loading methods I–VI. The continuous growth and turning of cracks in different minerals and mineral boundaries are marked with arrows in different colors.Figure 6
Continuous growth and deflection of cracks (in red) under six different loading methods. The yellow arrows refer to crack growth in large scales of kaolinite and illite, the green arrows refer to crack growth in smaller scales of kaolinite and illite surrounded by quartz and feldspar, the blue arrows show crack growth in the boundary between different minerals, and the black arrows refer to crack growth in quartz and feldspar.
(a)(b)(c)(d)(e)(f)


### 4.3.1. Propagation and Turning of Cracks in Kaolinite and Illite

Among the six loading methods, the cracks were relatively flat and smooth when they propagated in the larger grains of kaolinite and illite. The length of each segment of the crack through these regions was longer, and thus these regions accounted for a higher proportion of the total crack length. When the crack tip was blocked by quartz and feldspar, the crack turned and the turning angle was low.When small grains of kaolinite and illite were surrounded by quartz and feldspar, many small-scale inflections were produced inside the kaolinite and illite. In loading methods I–III, these inflections were relatively small. In loading methods IV–VI, there were more inflections with larger amplitudes.


### 4.3.2. Propagation and Turning of Cracks at Mineral Boundaries

In loading methods I–III, only a small proportion of the cracks propagated along the mineral boundary. These segments usually propagated for a short distance along the mineral boundary and then turned into a kaolinite or illite grain. The positions of the cracks along the mineral boundary in loading methods IV–VI were similar. In loading method IV, cracks had more inflections, whereas in loading method V, crack inflections were reduced. In loading method VI, the cracks were smooth with almost no inflection, forming a shear arc.


### 4.3.3. Propagation and Turning of Cracks in Quartz and Feldspar

The cracks rarely extended directly through quartz and feldspar. When it occurred at the end of the crack propagation process, the crack shape was relatively straight. However, when it occurred in the middle of the crack propagation process, there were more inflections.


## 4.3.1. Propagation and Turning of Cracks in Kaolinite and Illite

Among the six loading methods, the cracks were relatively flat and smooth when they propagated in the larger grains of kaolinite and illite. The length of each segment of the crack through these regions was longer, and thus these regions accounted for a higher proportion of the total crack length. When the crack tip was blocked by quartz and feldspar, the crack turned and the turning angle was low.When small grains of kaolinite and illite were surrounded by quartz and feldspar, many small-scale inflections were produced inside the kaolinite and illite. In loading methods I–III, these inflections were relatively small. In loading methods IV–VI, there were more inflections with larger amplitudes.


## 4.3.2. Propagation and Turning of Cracks at Mineral Boundaries

In loading methods I–III, only a small proportion of the cracks propagated along the mineral boundary. These segments usually propagated for a short distance along the mineral boundary and then turned into a kaolinite or illite grain. The positions of the cracks along the mineral boundary in loading methods IV–VI were similar. In loading method IV, cracks had more inflections, whereas in loading method V, crack inflections were reduced. In loading method VI, the cracks were smooth with almost no inflection, forming a shear arc.


## 4.3.3. Propagation and Turning of Cracks in Quartz and Feldspar

The cracks rarely extended directly through quartz and feldspar. When it occurred at the end of the crack propagation process, the crack shape was relatively straight. However, when it occurred in the middle of the crack propagation process, there were more inflections.


## 4.4. En Echelon Intermittent Fractures

In addition to continuous propagation and turning, cracks were also generated, extended, and connected in a discontinuous manner, that is, en echelon intermittent cracks were created (Figure7).Figure 7
The initiation and evolution of en echelon intermittent fractures under six different loading methods. Every loading method is shown with five figures that demonstrate the growth of en echelon intermittent fractures over time.
(a)(b)(c)(d)(e)(f)The propagation process of en echelon intermittent fractures is divided into five stages.Stage 1: the small-scale en echelon intermittent fractures initiate and form the first large-scale en echelon intermittent fracture. In loading methods I, V, and VI, the en echelon intermittent fractures spread to the upper left direction. In loading methods II and III, the en echelon intermittent fractures spread to the upper right direction; whereas in loading methods IV, en echelon intermittent fractures were initiated both in the upper left and upper right directions.Stage 2: In loading method I, the second and third large-scale en echelon intermittent fractures were almost formed simultaneously. The direction of the second en echelon intermittent fracture was to the upper left, and the direction of the third en echelon intermittent fracture was to the upper right. The second large-scale en echelon intermittent fracture in loading methods II–V was located in the middle of the sample, whereas it was located at the top of the sample in loading method VI. The directions of the second large-scale en echelon intermittent fractures in loading methods II–V were all to the upper right. From loading method II to loading method V, the length of the second large-scale en echelon intermittent fractures gradually decreased.Stage 3: the third large-scale en echelon intermittent fractures in loading methods II, IV, V, and VI were all formed on the top of the sample. During this stage, the first and second large-scale en echelon intermittent fractures moved toward each other but did not connect. The three large-scale en echelon intermittent fractures in loading method I were also connected. In loading method III, the third en echelon intermittent fracture formed in the middle of the model and connected to the crack at the lower position.Stage 4: in loading methods I, IV, V, and VI, the two large-scale en echelon intermittent fractures at the middle and lower positions connected at this stage. However, in loading methods II and III, the two en echelon intermittent fractures in the middle and upper positions connected.Stage 5: In this stage, the remaining en echelon intermittent fractures connected and the rock was broken.The connection type of the en echelon intermittent fractures of scales under different loading methods was quite different and is summarized as follows:(1)
Wing Crack Connections. After the formation of en echelon intermittent fractures, adjacent en echelon intermittent fractures were connected by generating wing cracks at both ends, and these resembled an S-shape or a reverse S-shape. The wing crack connection was evident in the connection of two large-scale en echelon intermittent fractures and was also produced in the connection of small-scale en echelon intermittent fractures. Wing crack connections appeared in all six loading methods.(2)
Fragment Connection. In loading methods II–VI, there were en echelon intermittent fractures that appeared to connect in the form of fragments. In loading method III, fragments were generated in the small-scale en echelon intermittent fracture connection. The fragments in loading methods II, IV, V, and VI were generated during the large-scale en echelon intermittent fracture connection. The fragment positioning for loading methods IV–VI was the same. With the increase in lateral pressure, the size of the fragments gradually decreased.(3)
X-Shaped Crack Connection. For loading methods I and IV, the en echelon intermittent fractures appeared to create an X-shaped connection. The scale of X-shaped cracks in loading method I was larger than that in loading method IV.(4)
Not Connected or Not Connected at the Crack Tip. The large-scale en echelon intermittent fracture in the lower part of loading method II propagated to the upper right direction, but the crack tips were not connected with other cracks. Instead, fragments were generated in the middle of the crack that connected with the other two penetrating cracks, and the en echelon intermittent fracture that extended to the upper right was abandoned by the main crack. In loading method III, the en echelon intermittent fracture in the lower position directly propagated to the upper right and did not connect with the cracks at the upper position. Consequently, the cracks at the upper position were abandoned by the main crack.


## 5. Discussion

In rock damage theory, microdamages converge to form microcracks, and then microcracks extend or connect with other microcracks to form macrocracks [46]. Based on this theory, the strain contour and AE location before the formation of macroscopic cracks in loading methods I and IV were derived as shown in Figure 8. The AE characteristics of microcracks before macrocrack formation were extracted by programming in Python, and the AE location map was drawn using MATLAB software. In tensile loading (loading method I), the deformations of different minerals were nonuniform. Owing to the lower tensile strength, kaolinite, and illite have greater strain compared to quartz and feldspar. Because of the widespread distribution of such minerals, a large number of widely distributed microcracks were generated in the rock. Most of them were tensile microcracks, whereas very few were shear microcracks. The widely distributed tensile microcracks significantly increased the potential spread range of the cracks. The cracks could extend freely or could connect to the remaining microcracks. This led to the formation of longer length and more complicated cracks. In loading method IV (shear loading), strain first occurred at both ends of the maximum shear stress line. The strain in the rest of the rock did not reach failure displacement, even when the strain of the quartz and feldspar minerals around the notch did not reach the failure displacement. There were no widely distributed microcracks in the front of the crack propagation path; only the right-inclined and parallel tensile microcracks were generated on the right side. When the mineral at the crack tip could not be directly sheared, the crack connected with the nearest tensile microcrack. This considerably reduced the potential spread range of cracks.Figure 8
In loading method I (tensile loading) and loading method IV (shear loading), the strain and microfracture before the formation of macro-cracks. (a, b) Strain contour before macro-crack formation in loading methods I and IV, respectively. (c, d) In loading methods I and IV, the AE location of the microcracks distributed in the rock before the macrocrack formed. Red dots are the AE location of tensile fractures, purple dots are the AE location of shear fractures, and the colors between the two are the AE location of tensile-shear fractures. The size of the dot represents the level of AE energy. The AE energy of microfracture before the macrocrack formed is relatively small, so the AE energy magnification of this figure is larger than the others in this paper.
(a)(b)(c)(d)Compared with loading method IV, the lateral pressure in loading methods V and VI was equivalent to directly increasing the horizontal tensile strength of the rock. However, the increment in shear strength was the normal stress multiplied by the friction coefficient that was less than 1; therefore, the increased amplitude of shear strength was much smaller than that of tensile strength. The minerals at the crack tip could be sheared directly, but the surrounding minerals had not reached its tensile strength. Lateral pressure inhibited the generation of tensile microcracks around the crack tip. As a result, the potential spread range of cracks was further reduced, and the length and complexity of the cracks were also reduced. However, the changing trend of the crack length and fractal dimension of loading methods I–III was opposite to the trend of the amount and proportion of tensile fractures. This can be explained with the following analysis. In loading method I, the direction of the maximum tensile force is horizontal. The tensile microcracks spread and tended to propagate in the vertical direction. In loading methods II and III (tensile-shear loading), the shear and tensile forces were superimposed, causing the maximum tensile force to rotate clockwise from the horizontal direction. The tensile microcracks tended to spread in the upper right direction; as a result, the main crack connected more with the tensile microcracks on the upper right of the rock. This led to an increase in the length and complexity of the cracks. This study showed that the length and complexity of the crack cannot be simply determined by measuring the number of tensile cracks, and the spread range of the tensile microcracks is the fundamental reason.In geology, it is believed that tortuous cracks occur due to tension, and the smooth and straight cracks occur due to shear [2]. This study demonstrated that different loading methods affecting the complexity of cracks were in the spread range of tensile microcracks. Tensile loading can not only cause bedding and fissures to fracture but can also produce tensile microcracks in widely distributed weak minerals such as kaolinite and illite. Therefore, the crack has a high degree of complexity and a large spread range. Shear loading only produces parallel tensile microcracks on one side of the shear line and can only rupture the cracks on one side. The crack propagates along the direction of the maximum shear force; therefore, the tensile microcrack on one side cannot be extended easily, and the possibility of the main crack and the tensile microcrack connecting is reduced. As a result, the shear crack has a low degree of complexity and a small spread range.The differences in the continuous propagation and turning of cracks in different minerals are mainly caused by different mechanical responses of the minerals. Figure9 shows the superposition of the AE location and crack path. Owing to the low tensile and shear strengths of weak minerals such as kaolinite and illite, microcracks easily produced tensile or shear fractures, so there was no need to connect with the surrounding tensile microcracks. Therefore, in most cases, the crack propagation was relatively straight and the steering was relatively smooth. When weak minerals were surrounded by hard minerals such as quartz and feldspar, there were more external barriers to crack propagation, causing many types of mineral fractures. As a result, more crack inflections occurred. The tensile and shear strengths at the mineral boundary were higher than that of kaolinite and illite. Therefore, in loading methods I–III, the crack only extended at the mineral boundary for a short period of time and then went back into the weak mineral. The fracture types were all tensile fractures. Loading method IV mainly produced shear fractures at the mineral boundary. When the crack deviated from the maximum shear stress line (centerline), it connected to the tensile fracture occurred in the adjacent kaolinite and illite and then returned to the shear centerline. As the lateral pressure in loading methods V and VI increased, the tensile strength of kaolinite and illite increased, and it was not easy to generate tensile fracture. Therefore, the crack continued to extend along the mineral boundary in the form of a shear arc. The tensile and shear strengths of quartz and feldspar were considerably high. The tensile strength was lower than the shear strength. Therefore, when cracks passed through the quartz and feldspar, they were all tensile fractures. In loading methods I and III, the situation when cracks passed through the quartz and feldspar was similar to that of breaking the last piece of strong supporting mineral, whereas in loading methods II and VI, the cracks broke the hard mineral after several inflections.Figure 9
AE location on the fracture path under six different load methods. The yellow arrows refer to the crack growth in a wide region of kaolinite and illite, the green arrows refer to the crack growth in a small region of kaolinite and illite that was surrounded by quartz and feldspar, the blue arrows refer to the crack growth in the boundary of different minerals, and the black arrows refer to the crack growth in quartz and feldspar. Red dots are the AE location of tensile fractures, purple dots are the AE location of shear fractures, and the colors between the two are the AE location of tensile-shear fractures. The size of the dot represents the level of AE energy.
(a)(b)(c)(d)(e)(f)En echelon intermittent cracks were first proposed in macroscopic fractures. Cui and Han [31] observed discontinuous microcracks in experiments for the first time and explained their hierarchical structure. Figure 10 shows the AE location of small-scale en echelon intermittent crack initiation. The small-scale en echelon intermittent cracks were all blocked by strong minerals such as quartz and feldspar, and no cracks were generated in the strong minerals. The heterogeneous strain on strong and weak minerals caused the weak minerals to reach failure displacement first, and thus, small-scale en echelon intermittent cracks were generated.Figure 10
The initiation of small-scale en echelon intermittent fractures under six different load methods. Red dots are the AE location of tensile fractures, purple dots are the AE location of shear fractures, and the colors between the two are the AE location of tensile-shear fractures. The size of the dot represents the level of AE energy.
(a)(b)(c)(d)(e)(f)Based on this, the occurrence of large-scale en echelon intermittent cracks was also related to the barrier of strong mineral aggregates. The AE location of the second large-scale en echelon intermittent crack formation is shown in Figure11, and there was no AE location in the strong mineral aggregates. The heterogeneous strain between the strong mineral aggregates and the weak minerals caused the weak minerals to be blocked by the strong mineral aggregates from reaching the failure displacement. Cui and Han [31] defined the position between the en echelon intermittent cracks as rock bridges. This study found that, regardless of the formation of large-scale or small-scale en echelon intermittent cracks, the rock bridges were composed of strong minerals or strong mineral aggregates.Figure 11
The AE location when the second large-scale en echelon intermittent fractures were formed under six different loading methods. The red arrows refer to the strong mineral aggregate. Red dots are the AE location of tensile fractures, purple dots are the AE location of shear fractures, and the colors between the two are the AE location of tensile-shear fractures. The size of the dot represents the level of AE energy.
(a)(b)(c)(d)(e)(f)Regardless of the scale of the en echelon intermittent crack, the spreading of cracks is related to the loading method. In loading methods I–III, en echelon intermittent cracks were mostly tensile fractures, and these were generated perpendicular to the direction of maximum tensile force. However, in loading methods IV–VI, en echelon intermittent cracks were mostly shear fractures, which were generated in the direction close to the shear centerline, and the greater the lateral pressure, the more obvious this tendency. It is worth noting that during the formation of large-scale en echelon intermittent cracks in loading methods IV and V, the shear strain was concentrated at both ends of the centerline; therefore, the middle part of the rock did not reach displacement due to shear failure but reached displacement from tensile failure. As the lateral pressure increased, the length of the crack in the middle of the rock decreased. In loading method VI, no cracks occurred in the middle part of the rock, and shear fracture occurred in the upper part. Consistent with the previous analysis, lateral pressure suppressed the occurrence of tensile fracture.Wing cracks were the most common connection types for en echelon intermittent cracks. When there were strong minerals or strong mineral aggregates in the en echelon intermittent cracks, the wing cracks bypassed the strong minerals or strong mineral aggregates to produce fragmented connections. The AE location of the fragment formation process for loading method IV is shown in Figure12. As the shear force could not cut the strong mineral aggregates (shown in Figure 8), tensile microcracks were formed on the right side of the shear crack. The cracks connected to the tensile microcracks around the boundary of the strong mineral aggregates and eventually formed fragments. In loading methods II and III, the tensile microcracks had a large distribution range, so there were many small fragments. The lateral pressure in loading methods V and VI restrained the generation range of tensile microcracks. As lateral pressure increased, the size of the fragments decreased.Figure 12
The formation process and AE location of fragment in loading method IV. Red dots are the AE location of tensile fractures, purple dots are the AE location of shear fractures, and the colors between the two are the AE location of tensile-shear fractures. The size of the dot represents the level of AE energy.
(a)(b)(c)The AE location of the X-shaped crack formation process for loading methods I and IV is shown in Figure13. In loading method I, tensile microcracks were widely generated, and en echelon intermittent cracks with different directions were formed. When the en echelon intermittent cracks extended to the center, they met in the middle to form an X-shaped crack. However, in loading method IV, the en echelon intermittent crack in the middle of the rock was blocked by strong minerals and could not continue to grow. Shear microcracks that extended to the upper left formed in the middle of the crack under shear loading. The crack in the lower part propagated and connected to it, forming an X-shaped crack. The spreading range of tensile microcracks is larger than that of shear microcracks. Therefore, the size of the X crack in loading method I is larger than that in loading method IV.Figure 13
X fracture formation and AE location for loading methods I ((a), (b), (c)) and IV ((d), (e), (f)). Red dots are the AE location of tensile fractures, purple dots are the AE location of shear fractures, and the colors between the two are the AE location of tensile-shear fractures. The size of the dot represents the level of AE energy.
(a)(b)(c)(d)(e)(f)Following this research, in-depth studies can still be conducted in many aspects. In the future, we could set different degrees of mineral heterogeneity by ratio and size in each loading method to further confirm the role of heterogeneity and the loading method. For different loading methods, further study is needed on the similarities and differences in the en echelon intermittent crack connections when compared with connections between naturally occurring cracks, pores, and fissures. Previous studies used fractal theory to summarize the multiscale and hierarchical distribution of cracks [11], and it was found that the en echelon intermittent cracks also show the characteristics of multiscale and hierarchical distribution. The fractal study of en echelon intermittent cracks should consider the effects of loading methods, mineral heterogeneity, pores, natural fissures, and bedding. Besides, how to make the en echelon intermittent cracks more fully connected is also an important issue for the improvement of shale gas recovery. This study provides a significant reference for future research for the loading method to increase the complexity of fractures. The combination of fractal dimension and crack length per unit area demonstrated here is a promising way to evaluate the effect of fracturing in future studies.


## 6. Conclusions

In this study, the Voronoi tessellation technique was used to establish a shale CZM based on mineral distribution, and six different loading methods were applied. Based on the rock damage theory, combined with the crack path, strain contour, and AE, the initiation and propagation of microcracks under different loading methods were compared and analyzed. The following are the main conclusions drawn from this study.(1)
The essence of different loading methods affecting crack length and complexity was the spread range of tensile microcracks. Tensile loading formed a large number of widely distributed tensile microcracks in the rock, so the crack length was longer and the twists and turns were more complicated. In contrast, shear loading only produced parallel tensile microcracks on the right side of the maximum shear line, and the cracks tended to extend along the direction of the maximum shear force, resulting in a shorter crack length and low complexity.(2)
The mechanical properties of various minerals and mineral boundaries are different; therefore, the mechanical responses to different loading methods are different. This leads to differences in continuous crack propagation and steering under different loading methods.(3)
The formation and propagation of en echelon intermittent cracks of different scales were mainly affected by the heterogeneity of minerals and mineral aggregates. The differences in displacements at failure of different minerals led to the generation of small-scale en echelon intermittent cracks separated by strong minerals. Similarly, the differences in displacements at failure of different mineral aggregates led to the generation of large-scale en echelon intermittent cracks separated by strong mineral aggregates.(4)
The spreading direction and connection form of en echelon intermittent cracks were mainly affected by the loading method. The en echelon intermittent cracks propagated perpendicular to the direction of maximum tensile stress during tensile and tensile-shear loading, with various connection methods. In shear and compression-shear loading, the en echelon intermittent cracks spread near the maximum shear stress line. As the lateral pressure increased, the cracks tended to concentrate more on the centerline, and the connection method became simpler.



---

*Source: 6621347-2021-01-29.xml*