# Numerical Test Study of the Microscale Failure Modes and Fractal Analysis of Lower Cambrian Shale Based on Digital Images

Authors: Hengtao Cui; Zhonghu Wu; Liping Li; Jing Wang; Shuguang Song; Shangqu Sun; Yujun Zuo; Hao Liu; Yili Lou
Journal: Advances in Civil Engineering
 (2020)

---

## Abstract

To reveal the effect of confining pressure on the mechanical properties and rupture modes of quartz-bearing shale, the shale core of the no. 3 block of Fenggang, Guizhou Province, China, was analyzed by nuclear magnetic resonance, polarized light microscope thin section observation and identification, and core X-ray whole-rock minerals diffraction analysis to determine the distribution of shale minerals in the Niutitang Formation. On the microscale, based on digital image processing technology, this paper characterizes the nonuniformity of minerals in shale, a numerical model that can reflect the true microstructure of shale. Then, the failure process of shale under different confining pressures was simulated. The results show that when the shale is loaded with vertical displacement under different confining pressures, the compressive strength and elastic modulus of the sample change significantly. The failure mode can be roughly divided into three types: the inverted V-shaped (0 MPa, 2 MPa, and 4 MPa), V-shaped (6 MPa and 8 MPa), and invertedZ-shaped (10 MPa). Since the development of fractal theory provides a new space for studying the damage and fracture of rocks, the damage evolution and failure process of shale can also be regarded as the fractal process of cracks, in which the fractal dimension is the core parameter. The calculation is different under different stress levels. The fractal dimension under the condition of confining pressure shows that the value of the fractal dimension is greatly affected by the effect of confining pressure. When the fractal dimension is higher, the fracture mode is more complicated, and the internal damage degree is more serious. The research results provide important theoretical guidance for shale gas fracturing production.


---

## Body

## 1. Introduction

In recent years, with the continuous consumption of conventional oil and gas, the supply and demand of oil and gas have become a global concern. However, the development of unconventional energy can make up for supply shortage [1, 2]. As a clean and efficient unconventional natural gas resource, shale gas has received increasing attention around the world as it has potential for sustainable development [3–7]. At present, countries in the world are increasing the exploration, exploitation, production, and utilization of shale gas. China’s shale gas resources are very rich, and the potential for the development and utilization of these resources is huge [8]. The recoverable shale gas resources are 3.6 × 1013 m3, slightly higher than the United States and the largest in the world. The prospects for exploitation are very broad [9–13]. The shale reservoir has ultralow porosity and low matrix permeability, which led to great difficulties and challenges to its effective exploitation. To realize commercial shale gas mining, artificial cracks must be created through large-scale hydraulic fracturing technology [14–19]. Therefore, the study of the mechanical properties and fracture modes of shale is of great significance to shale gas fracturing.To better understand the mechanical properties and fracture modes of shale, a large number of scholars have conducted research and obtained many valuable results. Jiang et al. [20] studied uniaxial compression tests, triaxial compression tests, and wave velocity determination tests. The research results show that the confining pressure and burial depth have significant effects on the mechanical properties, failure rules, and physical properties of shale. Mousavi et al. [21] studied the anisotropic mechanical behaviour of shale samples using the Brazilian disk test and proposed a model to simulate the crack propagation, the failure model of heterogeneous layered shale by calculating the tensile strength and fracture toughness of the samples. Wu et al. [22] used the RFPA2D software to perform a uniaxial compression numerical simulation of shale with calcite veins. The experiment showed that the compressive strength and elastic modulus under different directions showed significant anisotropy and that the destruction mode had a large impact. Liu et al. [23] conducted a triaxial mechanical test study of Silurian Longmaxi Formation shale samples under different confining pressures and the same strain rate. They revealed that the confining pressure and strain rate have a significant impact on the elastic modulus and failure modes of the isomechanical properties of shale. Duan et al. [24] obtained CT images under different stress levels by performing uniaxial compression tests on Longmaxi shale, which quantitatively characterized the evolution process of internal cracks in shale. Chuanliang et al. [25] performed an indoor uniaxial compression test, which showed that the mechanical characteristics of shale changed with axial stress and the coring angle along the plane. The research also proposed and verified a strength failure criterion for shale. Heng et al. [26] evaluated the anisotropy of shear strength by conducting direct shear tests with different shear angles on bedding-developed shale samples and used the test data to derive the expression of the shear stress concentration factor.In view of the research on the mechanical properties and failure modes of shale, and although scholars at home and abroad have performed much research and have reported many rich results, researchers have mainly focused on the macroscale failure of shale but have rarely considered the failure problems caused by the nonuniformity of microscale materials and the nonuniformity of stress in shale. The spatial distribution of minerals in shale is irregular. The macromechanical properties cannot fully reflect the influence of shale’s microscale heterogeneity on its microfracture evolution. It is difficult to reveal its microfracture mechanism.Taking the shale gas comprehensive exploration test area in Guizhou, southern China, as an example, the shale cores of the Niutitang Formation in FC-1 well that have been drilled in the area have been tested and analyzed by nuclear magnetic resonance and polarized optical microscopy observations and cores. X-ray whole-rock mineral diffraction analysis found that the shale contains many quartz minerals. As a common brittle mineral, quartz has an important influence on the failure mechanism and brittleness evaluation of shale. Combining digital image processing technology with the rock failure process analysis system RFPA2D, this paper establishes a numerical model of shale containing quartz filling which can reflect the true microstructure of the material to simulate the failure process of shale samples under different loading conditions. We further studied the effect of confining pressure on the strength, deformation, and fracture mode of shale. Combined with the basis of fractal theory, fractal dimension is used to quantitatively analyze the damage of shale microcracks. This study not only provides a reference for revealing the microfracture mechanism of shale, in the process of hydraulic fracturing on marine and continental shale in South China, but also provides an important theoretical support for guiding the exploration and exploitation of deep shale gas resources.


## 2. Geological Structure Characteristics

The no. 3 block of Fenggang is a comprehensive exploration and gas testing area in the northern Guizhou Province, South China, including the Fenggang and Meitan administrative areas, covering an area of 1167.49 km2. The no. 3 block of Fenggang is located at the intersection of the North-South structural belt, the North-Northeast structural belt, and the North-East structural belt of the Sichuan and Guizhou Provinces. The terrain of the no.3 block in the northwest is relatively high. While the terrain in the southeast is relatively low, the terrain is undulating, and the landform type is relatively complex [27] (Figure 1).Figure 1
Regional map of the no. 3 block of Fenggang, northern Guizhou [28].This area is located in the upper Yangtze platform area from the tectonic point of view and is consistent with the evolution of the Yangtze platform from the perspective of evolution characteristics. The area experienced the Xuefeng period, the early-Central Caledonian movement, the late Caledonian period, the Hexi period, the Indosinian period, the Yanshan period, and the Xishan period. Multiphase tectonic movements were superimposed, resulting in relatively complex and changeable structural shapes in the study area. Among them, the Yanshanian movement laid the basic pattern of the current structure, while the Himalayan period transformed and superimposed the Yanshanian period [28, 29].


## 3. Materials and Methods

### 3.1. Basic Information of Tested Samples

The lower Cambrian shale core in the study area is mainly taken from well FC-1, and the structural location lies in the northeast flank of the Fengyan anticline, at the Fenggang north-north-east structural deformation zone of the Zunyi fault-protrusion, located in the North Guizhou platform uplift of the Yangtze quasi-platform. The geological design of the well is a vertical coring well, with a design depth of 2515.00 m and a completed drilling depth of 2585.00 m, including a filling height of 6 m. The target layer for drilling is the lower Cambrian Niutitang Formation, with a stratum thickness of 102 m. The NMR test analysis of the rock samples in the FC-1 well revealed that a large number of nanoscale pores developed in the shale of the Niutitang Formation. The pore throats of the shale are mainly distributed between 0 and 0.1μm, including nanopores, although there are also a few micropores. The pore size distribution of shale is mainly between 0.001 μm∼0.01 μm and 0.01–0.4 μm, indicating that the shale pores in the study area are mainly nanopores, as shown in Figure 2 below.Figure 2
The pore size distribution of shale.The shale core samples of this stratum were selected for observation and identification by optical microslices and X-ray whole-rock mineral diffraction analysis (and clay mineral analysis). First, placing the shale core thin slices under a polarizing microscope for observation and identification allowed clear identification of the microstructural characteristics of the fracture. It was found that the fracture was highly filled and contained a large amount of quartz minerals (Figure3). Second, several groups of shale cores were selected for X-ray diffraction analysis of whole-rock minerals (Figure 4). The instrument model used was X’pert Powder. The analysis results show that the brittle minerals in the reservoir are relatively high, mainly including quartz, feldspar, pyrite, and calcite. Among them, shale brittle minerals are mainly quartz and feldspar, with a content of 45%–89% and an average of 78%. The overall content of calcite and iron ore is low, with an average content of 15%, and clay minerals are relatively low, with an average content of 15% (Figure 5) [27]. The shale with brittle mineral content is prone to form a staggered fracture network during fracturing, and higher shale gas production can be obtained during the mining process. Finally, further analysis of the clay minerals (Figure 6) found mainly illite (with an average content of 87%), which is clay mineral that contains a large number of pores, providing storage space for shale gas.Figure 3
Thin slice of shale core under microscope.Figure 4
X-ray diffraction analysis of whole-rock minerals.Figure 5
Distribution map of shale mineral content.Figure 6
Distribution of clay mineral content.


### 3.2. Numerical Simulation Method of the Shale Failure Process

#### 3.2.1. Characterization of Shale Microstructure Heterogeneity

In digital image processing, the image as a kind of information resource can effectively reflect the microscale inside the shale through different grey levels and colours. Digital image processing technology selects appropriate colour space for image segmentation according to the characteristics of different images and colours. Generally, two colour spaces are used, namely, the GRB colour space and the HIS colour space, and shale images are mostly grayscale images, and the HIS colour space is chosen for image segmentation [22, 30].First, digital image processing technology was used to characterize the heterogeneity of mineral components in shale and to classify the microscopic medium inside the shale. Figure7 is an image obtained by observing a shale core slice under a polarized optical microscope. The image format is a commonly used 24 BMP bitmap. In Figure 7(a), the black represents the shale matrix, and the white represents the quartz mineral. The segmentation threshold is determined by analyzing the I value change rule on the scan line AA'. The AA′ scan line passes through the mineral medium and the curve in Figure 7(b). The changes are compared. After many experiments and observations, the value of “I” on the shale matrix was found to be below 90, and the “I” value of quartz fluctuates above 90. Through repeated trials, “I” finally chooses “I’”= 90 as the threshold for image segmentation. The “I’” value distinguishes the shale matrix and quartz minerals and determines their geometric shape and spatial location distribution. Second, Figure 7(c) is the shale microstructure characterization image obtained after digital image processing, which is represented by light blue shale matrix, and light yellow indicates quartz. From the figure, it can be seen that the processed digital image, the geometric shape, and spatial distribution of shale and quartz are accurately characterized. Finally, the processed and characterized image was mapped to the finite element grid to establish a numerical model that can reflect its true microstructure.Figure 7
Scanning variation of shale digital image.
(a)(b)(c)


#### 3.2.2. Constitutive Relations of Unit Damage

Using elastic damage mechanics to analyse the mechanical properties of the unit at the microscale, the unit was found to be damaged due to the external force. As the degree of damage increases, the strength and stiffness of the element gradually decrease, and the elastic modulus of the damage can be expressed as [31–33](1)E=1−ωE0,where ω is the damage variable and E0,E are the initial elastic modulus of the material and the elastic modulus of the material after damage, respectively.The microunit behaves elastically in the initial stress state, and its mechanical properties can be reflected by parameters such as the elastic modulus and Poisson’s ratio. Due to the increasing stress on the microunit, when the damage threshold criterion is reached, the unit is damaged. The damage threshold criterion referred to here is divided into two types: the tensile criterion and the Mohr-Coulomb strength criterion. When the maximum tensile stress suffered by the element is greater than the extreme value of the tensile strength of the material, the element begins to undergo tensile damage and failure; when the stress state of the element reaches the Mohr-Coulomb criterion, the element also undergoes shear damage and failure. In this case, the stretching criterion has priority. If the element reaches the tensile criterion, the Mohr-Coulomb strength criterion can no longer be used to determine the damage of the element; only when the tensile strength criterion is not met is the Mohr-Coulomb strength criterion required [34, 35]. The elastic damage constitutive relationship of the element in the uniaxial state in the model is shown in Figures 8 and 9 below.Figure 8
Constitutive relationship of element elastic loss under uniaxial tensile stress.Figure 9
Constitutive relationship of element elastic loss under uniaxial compressive stress.According to Figure8, the green line segment is shown as the stress-strain curve relationship of the unit under uniaxial tension, and the damage constitutive relationship of the microunit under uniaxial tension is [36](2)ω=0,εt0≤ε<0,1−ftrE0ε,εtu≤ε<εt0,1,ε≤εtu,where ftr is the residual tensile strength; ftr=λft0=λE0εt0, where λ is the residual strength coefficient 0<λ≤1; ft0 is the uniaxial tensile strength; εt0 is the tensile strain corresponding to the elastic limit, also known as the threshold strain; that is, εt0=−ft0/E0; and εtu is the tensile strain of the element, and the element has been completely destroyed at this time. Formula (2) is further transformed to obtain the following formula [37, 38]:(3)ω=0,εt0≤ε<0,1−λεt0ε,εtu≤ε<εt0,1,ε≤εtu.The above describes the constitutive relationship of the element under uniaxial tension. When the element is under multiaxial stress, the mechanical properties of the microelement still satisfy the assumption of isotropy and elasticity. In the multiaxial stress state, when the element meets the tensile criterion, the equivalent principal strainε¯ reaches the threshold strain εt0, and the element is still in tension. The equivalent strain ε¯ is defined as [33, 39](4)ε¯=−−ε12+−ε22+−ε32,where ε1,ε2,ε3 are the three principal strains and the expression of the function ∗ is(5)x=x,x≥0,0,x<0.Substituting the equivalent principal strainε¯ for the strain ε in equation (3), the constitutive relationship of the element under the multiaxial stress state can be derived as(6)ω=0,εt0≤ε¯<0,1−λεt0ε¯,εtu<ε¯≤εt0,1,ε¯≤εtu.The above is the constitutive relationship of the elastic damage of the microunit under tension. It can be seen from Figure9 that the red line segment represents the stress-strain curve relationship of the element under uniaxial compression.To study the damage of the element under compression or shear stress, the Mohr-Coulomb strength criterion is selected as the second threshold criterion, and its definition expression is as follows [40, 41]:(7)F=σ1−1+sinφ1−sinφσ3≥fc0,where φ is the internal friction angle of the unit; σ1,σ3 represent large principal stress and small principal stress, respectively; and fc0 is the uniaxial compressive strength. When the microscopic element meets the Mohr-Coulomb strength criterion, the constitutive relationship of the element under uniaxial compression can be defined as follows:(8)ω=0,0≤ε<εc0,1−λεc0ε,ε≥εc0,where λ is the residual strength coefficient of the unit, when the unit is in compression λ=fcr/fc0, and ε0 is the maximum compression principal strain of the element; that is, εc0=fc0/E0.When the element is under the condition of multiaxial stress and the Mohr-Coulomb strength criterion is met, the element is damaged. Meanwhile, considering the influence of other principal stresses on the element damage, the expression of the definition of the element’s maximum compression principal strainεc0 can be obtained as [34](9)εc0=1E0fc0+1+sinφ1−sinφσ3−μσ1+σ2,where μ is Poisson’s ratio and fc0 is the uniaxial compressive strength. To study the evolution of element shear damage, the maximum compression principal strain ε1 of the element is used to replace the strain ε in formula (8), which can be extended to the definition of the constitutive relationship of element damage under multiaxial stress:(10)ω=0,0≤ε1<εc0,1−λεc0ε1,ε1≥εc0.


### 3.3. Establishment of the Numerical Model

Combining the characterized image with the finite element calculation method and using the unit grid mapping method, the image is divided into many small unit grids, and each small square grid is regarded as a pixel of the image. The actual image and pixel size can be represented by the finite element mesh, and then the Monte Carlo method is used to evaluate the parameters of the microelements. After that, we could generate the virtual microstructure of shale matrix and quartz minerals to establish a numerical model reflecting the real situation. At the same time, the nonuniformity coefficient is introduced into the numerical model, assuming that the mechanical properties of the shale matrix and quartz minerals follow a Weibull distribution [42, 43].(11)φa=ma0aa0m−1exp−aa0m,where m is the nonuniformity coefficient of the material and reflects the homogeneity of the shale material. When m is larger, it indicates that the mechanical properties of the material are more uniform. Alternatively, when m is smaller, it indicates that the mechanical properties of the material are more uneven. a represents the mechanical property parameters of the microunits, such as strength, Poisson’s ratio, and elastic modulus. a0 is the average value of the mechanical property parameters of the microunit. φa is a function of the statistical distribution density of the mechanical parameters of the microunits.By magnifying a small part of the numerical model in Figure10(a) and transforming it into a cell grid in Figure 10(b), the whole model is divided into 10×10 cells, where each cell size is 0.33×0.33 mm2. To eliminate the grid effect on crack growth, regular square cells are used here. The brighter the colour in the figure, the greater the elastic modulus; the darker the colour, the smaller the elastic modulus. The material mechanical parameters of the shale matrix and quartz are shown in Table 1 [22, 44].Figure 10
Numerical model and grid after shale transformation.Table 1
Mechanical parameters of shale matrix and quartz minerals.
MaterialElastic modulus (GPa)Compressive strength (MPa)Poisson’s ratioPull ratioInternal friction angle (°)Shale51.61450.221435Quartz963730.081560To study the influence of the shale microstructure on its fracture, based on the combination of digital image processing technology and the rock real fracture process analysis software RFPA2D-DIP, the six groups of quartz-bearing shale with confining pressuresQ of 0 MPa, 2 MPa, 4 MPa, 6 MPa, 8 MPa, and 10 MPa were tested. This problem is converted into a plane strain problem using a displacement-controlled loading method. The initial loading value s is 0.001 mm, and the loading amount per step Δs is 0.002 mm. The test is stopped until the model is completely destroyed. The schematic diagram of the model loading is shown in Figure 11.Figure 11
Model loading diagram.
(a)(b)


## 3.1. Basic Information of Tested Samples

The lower Cambrian shale core in the study area is mainly taken from well FC-1, and the structural location lies in the northeast flank of the Fengyan anticline, at the Fenggang north-north-east structural deformation zone of the Zunyi fault-protrusion, located in the North Guizhou platform uplift of the Yangtze quasi-platform. The geological design of the well is a vertical coring well, with a design depth of 2515.00 m and a completed drilling depth of 2585.00 m, including a filling height of 6 m. The target layer for drilling is the lower Cambrian Niutitang Formation, with a stratum thickness of 102 m. The NMR test analysis of the rock samples in the FC-1 well revealed that a large number of nanoscale pores developed in the shale of the Niutitang Formation. The pore throats of the shale are mainly distributed between 0 and 0.1μm, including nanopores, although there are also a few micropores. The pore size distribution of shale is mainly between 0.001 μm∼0.01 μm and 0.01–0.4 μm, indicating that the shale pores in the study area are mainly nanopores, as shown in Figure 2 below.Figure 2
The pore size distribution of shale.The shale core samples of this stratum were selected for observation and identification by optical microslices and X-ray whole-rock mineral diffraction analysis (and clay mineral analysis). First, placing the shale core thin slices under a polarizing microscope for observation and identification allowed clear identification of the microstructural characteristics of the fracture. It was found that the fracture was highly filled and contained a large amount of quartz minerals (Figure3). Second, several groups of shale cores were selected for X-ray diffraction analysis of whole-rock minerals (Figure 4). The instrument model used was X’pert Powder. The analysis results show that the brittle minerals in the reservoir are relatively high, mainly including quartz, feldspar, pyrite, and calcite. Among them, shale brittle minerals are mainly quartz and feldspar, with a content of 45%–89% and an average of 78%. The overall content of calcite and iron ore is low, with an average content of 15%, and clay minerals are relatively low, with an average content of 15% (Figure 5) [27]. The shale with brittle mineral content is prone to form a staggered fracture network during fracturing, and higher shale gas production can be obtained during the mining process. Finally, further analysis of the clay minerals (Figure 6) found mainly illite (with an average content of 87%), which is clay mineral that contains a large number of pores, providing storage space for shale gas.Figure 3
Thin slice of shale core under microscope.Figure 4
X-ray diffraction analysis of whole-rock minerals.Figure 5
Distribution map of shale mineral content.Figure 6
Distribution of clay mineral content.


## 3.2. Numerical Simulation Method of the Shale Failure Process

### 3.2.1. Characterization of Shale Microstructure Heterogeneity

In digital image processing, the image as a kind of information resource can effectively reflect the microscale inside the shale through different grey levels and colours. Digital image processing technology selects appropriate colour space for image segmentation according to the characteristics of different images and colours. Generally, two colour spaces are used, namely, the GRB colour space and the HIS colour space, and shale images are mostly grayscale images, and the HIS colour space is chosen for image segmentation [22, 30].First, digital image processing technology was used to characterize the heterogeneity of mineral components in shale and to classify the microscopic medium inside the shale. Figure7 is an image obtained by observing a shale core slice under a polarized optical microscope. The image format is a commonly used 24 BMP bitmap. In Figure 7(a), the black represents the shale matrix, and the white represents the quartz mineral. The segmentation threshold is determined by analyzing the I value change rule on the scan line AA'. The AA′ scan line passes through the mineral medium and the curve in Figure 7(b). The changes are compared. After many experiments and observations, the value of “I” on the shale matrix was found to be below 90, and the “I” value of quartz fluctuates above 90. Through repeated trials, “I” finally chooses “I’”= 90 as the threshold for image segmentation. The “I’” value distinguishes the shale matrix and quartz minerals and determines their geometric shape and spatial location distribution. Second, Figure 7(c) is the shale microstructure characterization image obtained after digital image processing, which is represented by light blue shale matrix, and light yellow indicates quartz. From the figure, it can be seen that the processed digital image, the geometric shape, and spatial distribution of shale and quartz are accurately characterized. Finally, the processed and characterized image was mapped to the finite element grid to establish a numerical model that can reflect its true microstructure.Figure 7
Scanning variation of shale digital image.
(a)(b)(c)


### 3.2.2. Constitutive Relations of Unit Damage

Using elastic damage mechanics to analyse the mechanical properties of the unit at the microscale, the unit was found to be damaged due to the external force. As the degree of damage increases, the strength and stiffness of the element gradually decrease, and the elastic modulus of the damage can be expressed as [31–33](1)E=1−ωE0,where ω is the damage variable and E0,E are the initial elastic modulus of the material and the elastic modulus of the material after damage, respectively.The microunit behaves elastically in the initial stress state, and its mechanical properties can be reflected by parameters such as the elastic modulus and Poisson’s ratio. Due to the increasing stress on the microunit, when the damage threshold criterion is reached, the unit is damaged. The damage threshold criterion referred to here is divided into two types: the tensile criterion and the Mohr-Coulomb strength criterion. When the maximum tensile stress suffered by the element is greater than the extreme value of the tensile strength of the material, the element begins to undergo tensile damage and failure; when the stress state of the element reaches the Mohr-Coulomb criterion, the element also undergoes shear damage and failure. In this case, the stretching criterion has priority. If the element reaches the tensile criterion, the Mohr-Coulomb strength criterion can no longer be used to determine the damage of the element; only when the tensile strength criterion is not met is the Mohr-Coulomb strength criterion required [34, 35]. The elastic damage constitutive relationship of the element in the uniaxial state in the model is shown in Figures 8 and 9 below.Figure 8
Constitutive relationship of element elastic loss under uniaxial tensile stress.Figure 9
Constitutive relationship of element elastic loss under uniaxial compressive stress.According to Figure8, the green line segment is shown as the stress-strain curve relationship of the unit under uniaxial tension, and the damage constitutive relationship of the microunit under uniaxial tension is [36](2)ω=0,εt0≤ε<0,1−ftrE0ε,εtu≤ε<εt0,1,ε≤εtu,where ftr is the residual tensile strength; ftr=λft0=λE0εt0, where λ is the residual strength coefficient 0<λ≤1; ft0 is the uniaxial tensile strength; εt0 is the tensile strain corresponding to the elastic limit, also known as the threshold strain; that is, εt0=−ft0/E0; and εtu is the tensile strain of the element, and the element has been completely destroyed at this time. Formula (2) is further transformed to obtain the following formula [37, 38]:(3)ω=0,εt0≤ε<0,1−λεt0ε,εtu≤ε<εt0,1,ε≤εtu.The above describes the constitutive relationship of the element under uniaxial tension. When the element is under multiaxial stress, the mechanical properties of the microelement still satisfy the assumption of isotropy and elasticity. In the multiaxial stress state, when the element meets the tensile criterion, the equivalent principal strainε¯ reaches the threshold strain εt0, and the element is still in tension. The equivalent strain ε¯ is defined as [33, 39](4)ε¯=−−ε12+−ε22+−ε32,where ε1,ε2,ε3 are the three principal strains and the expression of the function ∗ is(5)x=x,x≥0,0,x<0.Substituting the equivalent principal strainε¯ for the strain ε in equation (3), the constitutive relationship of the element under the multiaxial stress state can be derived as(6)ω=0,εt0≤ε¯<0,1−λεt0ε¯,εtu<ε¯≤εt0,1,ε¯≤εtu.The above is the constitutive relationship of the elastic damage of the microunit under tension. It can be seen from Figure9 that the red line segment represents the stress-strain curve relationship of the element under uniaxial compression.To study the damage of the element under compression or shear stress, the Mohr-Coulomb strength criterion is selected as the second threshold criterion, and its definition expression is as follows [40, 41]:(7)F=σ1−1+sinφ1−sinφσ3≥fc0,where φ is the internal friction angle of the unit; σ1,σ3 represent large principal stress and small principal stress, respectively; and fc0 is the uniaxial compressive strength. When the microscopic element meets the Mohr-Coulomb strength criterion, the constitutive relationship of the element under uniaxial compression can be defined as follows:(8)ω=0,0≤ε<εc0,1−λεc0ε,ε≥εc0,where λ is the residual strength coefficient of the unit, when the unit is in compression λ=fcr/fc0, and ε0 is the maximum compression principal strain of the element; that is, εc0=fc0/E0.When the element is under the condition of multiaxial stress and the Mohr-Coulomb strength criterion is met, the element is damaged. Meanwhile, considering the influence of other principal stresses on the element damage, the expression of the definition of the element’s maximum compression principal strainεc0 can be obtained as [34](9)εc0=1E0fc0+1+sinφ1−sinφσ3−μσ1+σ2,where μ is Poisson’s ratio and fc0 is the uniaxial compressive strength. To study the evolution of element shear damage, the maximum compression principal strain ε1 of the element is used to replace the strain ε in formula (8), which can be extended to the definition of the constitutive relationship of element damage under multiaxial stress:(10)ω=0,0≤ε1<εc0,1−λεc0ε1,ε1≥εc0.


## 3.2.1. Characterization of Shale Microstructure Heterogeneity

In digital image processing, the image as a kind of information resource can effectively reflect the microscale inside the shale through different grey levels and colours. Digital image processing technology selects appropriate colour space for image segmentation according to the characteristics of different images and colours. Generally, two colour spaces are used, namely, the GRB colour space and the HIS colour space, and shale images are mostly grayscale images, and the HIS colour space is chosen for image segmentation [22, 30].First, digital image processing technology was used to characterize the heterogeneity of mineral components in shale and to classify the microscopic medium inside the shale. Figure7 is an image obtained by observing a shale core slice under a polarized optical microscope. The image format is a commonly used 24 BMP bitmap. In Figure 7(a), the black represents the shale matrix, and the white represents the quartz mineral. The segmentation threshold is determined by analyzing the I value change rule on the scan line AA'. The AA′ scan line passes through the mineral medium and the curve in Figure 7(b). The changes are compared. After many experiments and observations, the value of “I” on the shale matrix was found to be below 90, and the “I” value of quartz fluctuates above 90. Through repeated trials, “I” finally chooses “I’”= 90 as the threshold for image segmentation. The “I’” value distinguishes the shale matrix and quartz minerals and determines their geometric shape and spatial location distribution. Second, Figure 7(c) is the shale microstructure characterization image obtained after digital image processing, which is represented by light blue shale matrix, and light yellow indicates quartz. From the figure, it can be seen that the processed digital image, the geometric shape, and spatial distribution of shale and quartz are accurately characterized. Finally, the processed and characterized image was mapped to the finite element grid to establish a numerical model that can reflect its true microstructure.Figure 7
Scanning variation of shale digital image.
(a)(b)(c)


## 3.2.2. Constitutive Relations of Unit Damage

Using elastic damage mechanics to analyse the mechanical properties of the unit at the microscale, the unit was found to be damaged due to the external force. As the degree of damage increases, the strength and stiffness of the element gradually decrease, and the elastic modulus of the damage can be expressed as [31–33](1)E=1−ωE0,where ω is the damage variable and E0,E are the initial elastic modulus of the material and the elastic modulus of the material after damage, respectively.The microunit behaves elastically in the initial stress state, and its mechanical properties can be reflected by parameters such as the elastic modulus and Poisson’s ratio. Due to the increasing stress on the microunit, when the damage threshold criterion is reached, the unit is damaged. The damage threshold criterion referred to here is divided into two types: the tensile criterion and the Mohr-Coulomb strength criterion. When the maximum tensile stress suffered by the element is greater than the extreme value of the tensile strength of the material, the element begins to undergo tensile damage and failure; when the stress state of the element reaches the Mohr-Coulomb criterion, the element also undergoes shear damage and failure. In this case, the stretching criterion has priority. If the element reaches the tensile criterion, the Mohr-Coulomb strength criterion can no longer be used to determine the damage of the element; only when the tensile strength criterion is not met is the Mohr-Coulomb strength criterion required [34, 35]. The elastic damage constitutive relationship of the element in the uniaxial state in the model is shown in Figures 8 and 9 below.Figure 8
Constitutive relationship of element elastic loss under uniaxial tensile stress.Figure 9
Constitutive relationship of element elastic loss under uniaxial compressive stress.According to Figure8, the green line segment is shown as the stress-strain curve relationship of the unit under uniaxial tension, and the damage constitutive relationship of the microunit under uniaxial tension is [36](2)ω=0,εt0≤ε<0,1−ftrE0ε,εtu≤ε<εt0,1,ε≤εtu,where ftr is the residual tensile strength; ftr=λft0=λE0εt0, where λ is the residual strength coefficient 0<λ≤1; ft0 is the uniaxial tensile strength; εt0 is the tensile strain corresponding to the elastic limit, also known as the threshold strain; that is, εt0=−ft0/E0; and εtu is the tensile strain of the element, and the element has been completely destroyed at this time. Formula (2) is further transformed to obtain the following formula [37, 38]:(3)ω=0,εt0≤ε<0,1−λεt0ε,εtu≤ε<εt0,1,ε≤εtu.The above describes the constitutive relationship of the element under uniaxial tension. When the element is under multiaxial stress, the mechanical properties of the microelement still satisfy the assumption of isotropy and elasticity. In the multiaxial stress state, when the element meets the tensile criterion, the equivalent principal strainε¯ reaches the threshold strain εt0, and the element is still in tension. The equivalent strain ε¯ is defined as [33, 39](4)ε¯=−−ε12+−ε22+−ε32,where ε1,ε2,ε3 are the three principal strains and the expression of the function ∗ is(5)x=x,x≥0,0,x<0.Substituting the equivalent principal strainε¯ for the strain ε in equation (3), the constitutive relationship of the element under the multiaxial stress state can be derived as(6)ω=0,εt0≤ε¯<0,1−λεt0ε¯,εtu<ε¯≤εt0,1,ε¯≤εtu.The above is the constitutive relationship of the elastic damage of the microunit under tension. It can be seen from Figure9 that the red line segment represents the stress-strain curve relationship of the element under uniaxial compression.To study the damage of the element under compression or shear stress, the Mohr-Coulomb strength criterion is selected as the second threshold criterion, and its definition expression is as follows [40, 41]:(7)F=σ1−1+sinφ1−sinφσ3≥fc0,where φ is the internal friction angle of the unit; σ1,σ3 represent large principal stress and small principal stress, respectively; and fc0 is the uniaxial compressive strength. When the microscopic element meets the Mohr-Coulomb strength criterion, the constitutive relationship of the element under uniaxial compression can be defined as follows:(8)ω=0,0≤ε<εc0,1−λεc0ε,ε≥εc0,where λ is the residual strength coefficient of the unit, when the unit is in compression λ=fcr/fc0, and ε0 is the maximum compression principal strain of the element; that is, εc0=fc0/E0.When the element is under the condition of multiaxial stress and the Mohr-Coulomb strength criterion is met, the element is damaged. Meanwhile, considering the influence of other principal stresses on the element damage, the expression of the definition of the element’s maximum compression principal strainεc0 can be obtained as [34](9)εc0=1E0fc0+1+sinφ1−sinφσ3−μσ1+σ2,where μ is Poisson’s ratio and fc0 is the uniaxial compressive strength. To study the evolution of element shear damage, the maximum compression principal strain ε1 of the element is used to replace the strain ε in formula (8), which can be extended to the definition of the constitutive relationship of element damage under multiaxial stress:(10)ω=0,0≤ε1<εc0,1−λεc0ε1,ε1≥εc0.


## 3.3. Establishment of the Numerical Model

Combining the characterized image with the finite element calculation method and using the unit grid mapping method, the image is divided into many small unit grids, and each small square grid is regarded as a pixel of the image. The actual image and pixel size can be represented by the finite element mesh, and then the Monte Carlo method is used to evaluate the parameters of the microelements. After that, we could generate the virtual microstructure of shale matrix and quartz minerals to establish a numerical model reflecting the real situation. At the same time, the nonuniformity coefficient is introduced into the numerical model, assuming that the mechanical properties of the shale matrix and quartz minerals follow a Weibull distribution [42, 43].(11)φa=ma0aa0m−1exp−aa0m,where m is the nonuniformity coefficient of the material and reflects the homogeneity of the shale material. When m is larger, it indicates that the mechanical properties of the material are more uniform. Alternatively, when m is smaller, it indicates that the mechanical properties of the material are more uneven. a represents the mechanical property parameters of the microunits, such as strength, Poisson’s ratio, and elastic modulus. a0 is the average value of the mechanical property parameters of the microunit. φa is a function of the statistical distribution density of the mechanical parameters of the microunits.By magnifying a small part of the numerical model in Figure10(a) and transforming it into a cell grid in Figure 10(b), the whole model is divided into 10×10 cells, where each cell size is 0.33×0.33 mm2. To eliminate the grid effect on crack growth, regular square cells are used here. The brighter the colour in the figure, the greater the elastic modulus; the darker the colour, the smaller the elastic modulus. The material mechanical parameters of the shale matrix and quartz are shown in Table 1 [22, 44].Figure 10
Numerical model and grid after shale transformation.Table 1
Mechanical parameters of shale matrix and quartz minerals.
MaterialElastic modulus (GPa)Compressive strength (MPa)Poisson’s ratioPull ratioInternal friction angle (°)Shale51.61450.221435Quartz963730.081560To study the influence of the shale microstructure on its fracture, based on the combination of digital image processing technology and the rock real fracture process analysis software RFPA2D-DIP, the six groups of quartz-bearing shale with confining pressuresQ of 0 MPa, 2 MPa, 4 MPa, 6 MPa, 8 MPa, and 10 MPa were tested. This problem is converted into a plane strain problem using a displacement-controlled loading method. The initial loading value s is 0.001 mm, and the loading amount per step Δs is 0.002 mm. The test is stopped until the model is completely destroyed. The schematic diagram of the model loading is shown in Figure 11.Figure 11
Model loading diagram.
(a)(b)


## 4. Test Results

### 4.1. Mechanical Properties of Shale under the Confining Pressure Effect

The stress-strain curves of shale samples under different confining pressures are shown in Figure12. The abscissa in the figure represents the strain, and the ordinate represents the stress. It can be clearly seen that, with the increase of confining pressure, the peak strength and residual strength also increase. Reflecting the significant confining pressure effect, which is consistent with the results of the general test, shale containing quartz minerals mainly exhibits elastic brittle failure. The stress value of the sample undergoes changes in different stages with the increase of strain, which can be roughly divided into three stages: linear growth stage, postpeak softening stage, and gentle stage. First, there was no obvious compaction stage because, in the numerical test, when the boundary conditions were set, the load in Y direction was loaded by displacement and the initial value was 0.001 mm, which led to a curve that did not pass the origin at the beginning, different from the conventional stress-strain curve. Second, the nonlinear performance before the peak is not obvious, the curve shows a linear growth stage, and the confining pressure curves basically overlap. When the peak strength is reached, with the continuous loading of displacement, the specimen is unstable and damaged, and the internal stress is released. However, the load-bearing capacity is not completely lost, and the curve drops sharply, which is the postpeak softening stage. As the confining pressure increases, the drop becomes smaller, which is mainly caused by the confining pressure effect. Finally, the curve shows a gentle stage, and the stress release tends to be stable. It is worth noting that, under a low confining pressure (0 MPa and 2 MPa), when the peak strength is reached, the curve has the largest steep drop and the curve experiences multiple stress drops, showing the characteristics of gradual energy release. According to classical elastoplastic mechanics, the calculation formulas for the compressive strength and elastic modulus of the shale specimens are as follows:(12)σc=PcA,where Pc is the maximum compressive load when the shale specimen fails and A is the cross-sectional area of the specimen.(13)E=σcε.Figure 12
Stress-strain curves under different confining pressures.According to formulas (12) and (13), the summary results of the shale compressive strength and elastic modulus can be obtained in Table 2.Table 2
Statistics of numerical calculation results of elastic modulus and compressive strength of shale.
Confining pressure (MPa)Elastic modulus (GPa)Compressive strength (MPa)053.2485.82252.06117.13452.32129.22651.84139.45852.44146.291052.80153.66Figure13 is a graph showing the change trend of the compressive strength and elastic modulus of shale under different confining pressures. In the figure, the abscissa indicates the confining pressure, the left side of the ordinate indicates the compressive strength, and the right side of the ordinate indicates the elastic modulus. It can be seen from the figure that the confining pressure effect shows obvious changes in compressive strength and deformability. When the confining pressure is equal to 0 MPa, the compressive strength is the smallest, 85.82 MPa. When the confining pressure reaches 10 MPa, the compressive strength is the largest, 153.66 MPa, and the difference between the maximum value and the minimum value is 67.84 MPa. When the confining pressure is equal to 6 MPa, the elastic modulus is the smallest, 51.84 GPa. When there is no confining pressure, the elastic modulus is the largest, 53.24 GPa, and the difference between the maximum elastic modulus and minimum elastic modulus is 1.4 GPa. Among them, due to the increase in the confining pressure, the compressive strength increases approximately linearly. The confining pressure increased from 0 MPa to 2 MPa. At this time, the slope of the curve is the largest and is approximately 45° from the horizontal axis, but the elastic modulus decreases. From 2 MPa to 10 MPa, the slope is basically the same, but the elastic modulus tends to increase first, then decrease, and, finally, gradually increase. The reason for this is mainly due to the heterogeneity of the internal microstructure in shale, which results in different stress concentrations under different confining pressure loading. As a result, the inducing microscopic and small element damage then gradually transitions to the macroscopic damage, which fully shows that the mechanical parameters of shale are greatly affected by the confining pressure effect.Figure 13
Variation trend of compressive strength and elastic modulus of shale under different confining pressures.


### 4.2. Fracture Damage Evolution Process of Shale under the Confining Pressure Effect

In the evolution process of shale fracture damage under different confining pressure loading, when the stress intensity reaches the peak, the sample fractures. The fracture mode and corresponding acoustic emission diagram are shown in Figure14. It can be seen from the figure that when the confining pressure is 0 MPa, the cracks start from the relatively concentrated position of the upper quartz minerals, approximately at an angle of 45°. After, it extending to both ends, and the cracks gradually expand until the instability shows an inverted V-shaped failure. When the confining pressure is 2 MPa, the cracks rupture along the angle of 60° from the left side of the shale, where the quartz minerals are concentrated. Soon after, sample results in obvious cracks, which continue to extend to both ends, bifurcating to produce many small cracks, and show an inverted V-shaped failure finally. When the confining pressure is 4 MPa, a small crack will be generated from the middle part of the quartz mineral concentration area on the left along the angle of 45°, which will continue to expand to the two ends, then penetrate the entire sample, and finally show an inverted V-shaped failure. When the confining pressure is 6 MPa, the cracks start from the quartz-concentrated area in the middle of the shale. At the same time, small cracks also occur at the upper left end, and the cracks propagate downwards until the instability failure presents a V-shaped failure. When the confining pressure is 8 MPa, the cracks first cracked from the quartz concentration area at the bottom of the shale, producing two small cracks, each extending along an angle of 45° at both ends, and finally the damage appeared in an inverted V shape. When the confining pressure is 10 MPa, the cracks also start to crack along the bottom, then continue to expand to both ends, and finally show a Z-shaped failure. In summary, due to the different stratum depths of shale, the confining pressure effect exists, which effectively limits the development of small microcracks. In addition, shale mineral components usually contain quartz, which is a brittle mineral with a weak structure and uneven distribution. It is affected by the mechanical parameters of the shale sample, and cracks occur first during the loading process. Therefore, the failure modes of shale under different confining pressures can be summarized into three types: inverted V-shaped (0 MPa, 2 MPa, and 4 MPa), V-shaped (6 MPa and 8 MPa), and inverted Z-shaped (10 MPa).Figure 14
Fracture modes diagram of shale under different confining pressures and corresponding acoustic emission diagram. (a) 0 MPa. (b) 2 MPa. (c) 4 MPa. (d) 6 MPa. (e) 8 MPa. (f) 10 MPa.It can be seen from Figure14 that different fracture modes are also different from the corresponding acoustic emission. In the figure, red indicates that the unit has tensile damage at the current loading step, yellow indicates that the unit has compression shear damage at the current loading step, and black indicates that the unit has been completely damaged at the current loading step. From the microscopic point of view, at low confining pressures (0 MPa and 2 MPa), the failure of the specimen is mainly due to the coupling of tensile and shear failure. The unit failure is mainly due to the tensile damage (red), which also reflects that the final macrocracks of the specimen are mainly caused by tensile failure. Because quartz is a brittle mineral, the concentrated area of the quartz mineral distribution is damaged first when loading. In the presence of confining pressure, the unit failure is mainly compression shear damage (yellow), which is also accompanied by tensile damage. As shale is easy to compress and not easy to tense, during the loading process, the tensile stress will first reach the tensile strength of rock, which meets the tensile criterion. When the tensile stress meets the tensile criterion, tensile damage occurs, and then compression shear failure occurs again. From the perspective of acoustic emission localization, the development of microcracks within the specimen determines the final macroscopic failure mode. The nonuniformity of the shale microstructure has an important influence on its mechanical properties and fracture mode.


### 4.3. Fractal Characteristics of the Shale Micro Failure Process

Fractal theory, as a very popular and active new theory and discipline, belongs to a branch of nonlinear science. Fractal theory was first proposed by the American French mathematician Roth and can be used to describe many irregular and self-similarities in nature. The geometric shape of features has been widely used in geology, computer science, medicine, and many other fields [45–48].As we all know, the damage evolution process of rock can also be regarded as the fractal process of cracks. In 1996, the academician Xie Heping applied fractal theory to the field of rock engineering and created the “fractal-rock mechanics” theory, which unified the macroscopic, mesoscopic, and microscopic theoretical studies of rock in a new theoretical system. Research on the application of fractal theory in geotechnical engineering is open. Rock rupture is mainly due to the initiation of microcracks in rock. The continuous expansion of a crack is accompanied by the rapid propagation of elastic sound waves and the release of energy in the rock, namely, acoustic emissions. Rock has fractal characteristics, including microfractures, macrofractures, and acoustic emissions, which form the most basic physical quantity of rock failure. A microunit fracture is generated for each acoustic emission event, so acoustic emissions can be used to analyze its fractal characteristics.The fractal dimension is the core parameter describing the fractal theory. The fractal dimension of the microstructure can be used as a quantitative parameter for the growth of shale microscopic cracks. The box dimension in the fractal dimension can more accurately represent the self-similarity of objects. The calculation principle is relatively simple, and the fractal dimension of irregular geometry can be accurately determined. The box dimension is used in this article, and the calculation formula is as follows:(14)Ds=−limr⟶0logNrlogr,where Ds is the self-similar fractal dimension of the destruction area, r is the side length of the square box, and Nr is the number of boxes required to cover the damaged area in the entire figure with a square box with side length r.The calculation of the fractal dimension of the shale compression test under different confining pressures is performed on the MATLAB software platform by writing a fractal dimension-related calculation program. The main steps of the algorithm are as follows: the first step is to process an image of acoustic emissions as a512×512 pixel image and read the processed image and then to preprocess the image grayscale, to process the image binary, and to store the relevant data. The second step is to cover the processed image with a square box with a side length of r, where r=2n,n=1,2,…,512, to count the number of square boxes in the failure area of the rock sample Nr, and to save the relevant data. Last, the double logarithmic coordinates of r and Nr are established for linear regression fitting, and the fractal dimension values under different confining pressures are calculated.Figure15 is the fitting curve of the fractal characteristics of the shale failure area when the stress level σ/σmax is 100% and the confining pressure is 10 MPa. The correlation coefficient R2=0.9968 can reflect the fact that the fractal dimension has high credibility, and the fractal dimension is Ds=1.9210.Figure 15
Fractal characteristics of the failure area when the stress level is 100% and the confining pressure is 10 MPa.Figures16–18 reflect various changes in stress levels, AE energy, and fractal dimensions under different confining pressures, respectively. Table 3 lists the values of the AE energy and the fractal dimension under six different confining pressures and different stress levels.Figure 16
Relationship between stress levels and AE Energy under different confining pressures.Figure 17
Relationship between stress levels and fractal dimensions under different confining pressures.Figure 18
Relationship between stress level, confining pressure, and fractal dimension.Table 3
AE energy and fractal dimension values under different stress levels and different confining pressures.
Stress levelConfining pressure(MPa)10%20%30%40%50%60%70%80%90%100%0AE energy00000.00010.00130.01120.05620.303540.3338Ds0000.47520.66170.83301.20601.18801.36701.51302AE energy0000.00010.0020.01960.15321.20846.028664.1970Ds0000.53290.7721.00201.26101.50001.59201.67704AE energy00000.00120.02950.24131.518419.9760108.7755Ds0000.41970.64551.00301.26401.5061.70101.78406AE energy00000.00200.03320.30252.153624.1428154.0042Ds0000.40760.63600.99521.28201.52301.7251.83508AE energy000000.0190.21061.797133.5573279.7562Ds00000.41970.88181.21101.50401.71501.893010AE energy00000.00060.01750.20691.761545.5702377.1735Ds00000.50120.82741.16801.48501.71301.9210The relationship between the stress level and the AE energy under different confining pressures is shown in Figure16. In the figure, the abscissa indicates different stress level values, and the ordinate indicates the energy value. When the stress level is less than 80%, the AE energy at each confining pressure is zero, and there is no trend. When the stress level is between 80% and 90%, the increase of the AE energy curve is relatively low, and the trend of change is relatively gentle and in a calm state. When the stress level exceeds 90% and reaches 100%, the change trend is more significant and is in a surge state. When the confining pressure is 10 MPa, the AE energy curve increases rapidly to the maximum, and when the confining pressure is 0 MPa, the AE energy value is the minimum. The results show that the energy released from the damage and fracture of the specimen is the largest when the confining pressure is 10 MPa. When the damage degree is the most severe, the internal damage is also the most serious, the fracture mode is relatively complex, and the confining pressure effect on the acoustic emission energy release is significant.From the relationship between the stress level and the fractal dimension under different confining pressures, as shown in Figure17, It can be seen that as the damage inside the specimen continues to accumulate, the fractal dimension value also increases with the increase of the confining pressure. When the stress level is less than 30%, the fractal dimension of each confining pressure state is 0. When the stress level reaches 30%, the fractal dimension values with confining pressures of 0 MPa, 2 MPa, 4 MPa, and 6 MPa change, and the curve change trend is relatively slow. When the confining pressure is 0 MPa, the fractal dimension value is the smallest, which is 1.5120, and the stress level of 8 MPa and 10 MPa reaches 40% and the fractal dimension value changes. However, the curve increases sharply, and the change trend is significant. When the confining pressure is 10 MPa, the fractal dimension reaches the maximum value of 1.9210. Further analysis shows that the fracture crack of the specimen with a confining pressure of 0 MPa is an inverted V-shaped, and the microfracture starts from the concentrated position of the quartz minerals in the upper part of the shale. At this time, the fractal dimension value is also the smallest. With increasing confining pressure, the fractal dimension increases gradually. When the confining pressure is 10 MPa, the fracture crack of the specimen is an inverted Z shape, the microcrack starts to fracture from the bottom, and the fractal dimension reaches the maximum value. The more complex the failure, the more the internal damage of the specimen. Therefore, when the fractal dimension value is larger, the fracture mode is more complicated, and the fractal dimension value gradually increases with crack propagation. The fractal dimension of the microstructure can be used as the quantitative parameter of the shale microcrack propagation.It can be clearly seen from Figure18 that the relationship between the confining pressure, stress level, and fractal dimension changes. Here, the fractal perspective is used to illustrate that the confining pressure effect has a significant influence on the shale failure mode.


## 4.1. Mechanical Properties of Shale under the Confining Pressure Effect

The stress-strain curves of shale samples under different confining pressures are shown in Figure12. The abscissa in the figure represents the strain, and the ordinate represents the stress. It can be clearly seen that, with the increase of confining pressure, the peak strength and residual strength also increase. Reflecting the significant confining pressure effect, which is consistent with the results of the general test, shale containing quartz minerals mainly exhibits elastic brittle failure. The stress value of the sample undergoes changes in different stages with the increase of strain, which can be roughly divided into three stages: linear growth stage, postpeak softening stage, and gentle stage. First, there was no obvious compaction stage because, in the numerical test, when the boundary conditions were set, the load in Y direction was loaded by displacement and the initial value was 0.001 mm, which led to a curve that did not pass the origin at the beginning, different from the conventional stress-strain curve. Second, the nonlinear performance before the peak is not obvious, the curve shows a linear growth stage, and the confining pressure curves basically overlap. When the peak strength is reached, with the continuous loading of displacement, the specimen is unstable and damaged, and the internal stress is released. However, the load-bearing capacity is not completely lost, and the curve drops sharply, which is the postpeak softening stage. As the confining pressure increases, the drop becomes smaller, which is mainly caused by the confining pressure effect. Finally, the curve shows a gentle stage, and the stress release tends to be stable. It is worth noting that, under a low confining pressure (0 MPa and 2 MPa), when the peak strength is reached, the curve has the largest steep drop and the curve experiences multiple stress drops, showing the characteristics of gradual energy release. According to classical elastoplastic mechanics, the calculation formulas for the compressive strength and elastic modulus of the shale specimens are as follows:(12)σc=PcA,where Pc is the maximum compressive load when the shale specimen fails and A is the cross-sectional area of the specimen.(13)E=σcε.Figure 12
Stress-strain curves under different confining pressures.According to formulas (12) and (13), the summary results of the shale compressive strength and elastic modulus can be obtained in Table 2.Table 2
Statistics of numerical calculation results of elastic modulus and compressive strength of shale.
Confining pressure (MPa)Elastic modulus (GPa)Compressive strength (MPa)053.2485.82252.06117.13452.32129.22651.84139.45852.44146.291052.80153.66Figure13 is a graph showing the change trend of the compressive strength and elastic modulus of shale under different confining pressures. In the figure, the abscissa indicates the confining pressure, the left side of the ordinate indicates the compressive strength, and the right side of the ordinate indicates the elastic modulus. It can be seen from the figure that the confining pressure effect shows obvious changes in compressive strength and deformability. When the confining pressure is equal to 0 MPa, the compressive strength is the smallest, 85.82 MPa. When the confining pressure reaches 10 MPa, the compressive strength is the largest, 153.66 MPa, and the difference between the maximum value and the minimum value is 67.84 MPa. When the confining pressure is equal to 6 MPa, the elastic modulus is the smallest, 51.84 GPa. When there is no confining pressure, the elastic modulus is the largest, 53.24 GPa, and the difference between the maximum elastic modulus and minimum elastic modulus is 1.4 GPa. Among them, due to the increase in the confining pressure, the compressive strength increases approximately linearly. The confining pressure increased from 0 MPa to 2 MPa. At this time, the slope of the curve is the largest and is approximately 45° from the horizontal axis, but the elastic modulus decreases. From 2 MPa to 10 MPa, the slope is basically the same, but the elastic modulus tends to increase first, then decrease, and, finally, gradually increase. The reason for this is mainly due to the heterogeneity of the internal microstructure in shale, which results in different stress concentrations under different confining pressure loading. As a result, the inducing microscopic and small element damage then gradually transitions to the macroscopic damage, which fully shows that the mechanical parameters of shale are greatly affected by the confining pressure effect.Figure 13
Variation trend of compressive strength and elastic modulus of shale under different confining pressures.


## 4.2. Fracture Damage Evolution Process of Shale under the Confining Pressure Effect

In the evolution process of shale fracture damage under different confining pressure loading, when the stress intensity reaches the peak, the sample fractures. The fracture mode and corresponding acoustic emission diagram are shown in Figure14. It can be seen from the figure that when the confining pressure is 0 MPa, the cracks start from the relatively concentrated position of the upper quartz minerals, approximately at an angle of 45°. After, it extending to both ends, and the cracks gradually expand until the instability shows an inverted V-shaped failure. When the confining pressure is 2 MPa, the cracks rupture along the angle of 60° from the left side of the shale, where the quartz minerals are concentrated. Soon after, sample results in obvious cracks, which continue to extend to both ends, bifurcating to produce many small cracks, and show an inverted V-shaped failure finally. When the confining pressure is 4 MPa, a small crack will be generated from the middle part of the quartz mineral concentration area on the left along the angle of 45°, which will continue to expand to the two ends, then penetrate the entire sample, and finally show an inverted V-shaped failure. When the confining pressure is 6 MPa, the cracks start from the quartz-concentrated area in the middle of the shale. At the same time, small cracks also occur at the upper left end, and the cracks propagate downwards until the instability failure presents a V-shaped failure. When the confining pressure is 8 MPa, the cracks first cracked from the quartz concentration area at the bottom of the shale, producing two small cracks, each extending along an angle of 45° at both ends, and finally the damage appeared in an inverted V shape. When the confining pressure is 10 MPa, the cracks also start to crack along the bottom, then continue to expand to both ends, and finally show a Z-shaped failure. In summary, due to the different stratum depths of shale, the confining pressure effect exists, which effectively limits the development of small microcracks. In addition, shale mineral components usually contain quartz, which is a brittle mineral with a weak structure and uneven distribution. It is affected by the mechanical parameters of the shale sample, and cracks occur first during the loading process. Therefore, the failure modes of shale under different confining pressures can be summarized into three types: inverted V-shaped (0 MPa, 2 MPa, and 4 MPa), V-shaped (6 MPa and 8 MPa), and inverted Z-shaped (10 MPa).Figure 14
Fracture modes diagram of shale under different confining pressures and corresponding acoustic emission diagram. (a) 0 MPa. (b) 2 MPa. (c) 4 MPa. (d) 6 MPa. (e) 8 MPa. (f) 10 MPa.It can be seen from Figure14 that different fracture modes are also different from the corresponding acoustic emission. In the figure, red indicates that the unit has tensile damage at the current loading step, yellow indicates that the unit has compression shear damage at the current loading step, and black indicates that the unit has been completely damaged at the current loading step. From the microscopic point of view, at low confining pressures (0 MPa and 2 MPa), the failure of the specimen is mainly due to the coupling of tensile and shear failure. The unit failure is mainly due to the tensile damage (red), which also reflects that the final macrocracks of the specimen are mainly caused by tensile failure. Because quartz is a brittle mineral, the concentrated area of the quartz mineral distribution is damaged first when loading. In the presence of confining pressure, the unit failure is mainly compression shear damage (yellow), which is also accompanied by tensile damage. As shale is easy to compress and not easy to tense, during the loading process, the tensile stress will first reach the tensile strength of rock, which meets the tensile criterion. When the tensile stress meets the tensile criterion, tensile damage occurs, and then compression shear failure occurs again. From the perspective of acoustic emission localization, the development of microcracks within the specimen determines the final macroscopic failure mode. The nonuniformity of the shale microstructure has an important influence on its mechanical properties and fracture mode.


## 4.3. Fractal Characteristics of the Shale Micro Failure Process

Fractal theory, as a very popular and active new theory and discipline, belongs to a branch of nonlinear science. Fractal theory was first proposed by the American French mathematician Roth and can be used to describe many irregular and self-similarities in nature. The geometric shape of features has been widely used in geology, computer science, medicine, and many other fields [45–48].As we all know, the damage evolution process of rock can also be regarded as the fractal process of cracks. In 1996, the academician Xie Heping applied fractal theory to the field of rock engineering and created the “fractal-rock mechanics” theory, which unified the macroscopic, mesoscopic, and microscopic theoretical studies of rock in a new theoretical system. Research on the application of fractal theory in geotechnical engineering is open. Rock rupture is mainly due to the initiation of microcracks in rock. The continuous expansion of a crack is accompanied by the rapid propagation of elastic sound waves and the release of energy in the rock, namely, acoustic emissions. Rock has fractal characteristics, including microfractures, macrofractures, and acoustic emissions, which form the most basic physical quantity of rock failure. A microunit fracture is generated for each acoustic emission event, so acoustic emissions can be used to analyze its fractal characteristics.The fractal dimension is the core parameter describing the fractal theory. The fractal dimension of the microstructure can be used as a quantitative parameter for the growth of shale microscopic cracks. The box dimension in the fractal dimension can more accurately represent the self-similarity of objects. The calculation principle is relatively simple, and the fractal dimension of irregular geometry can be accurately determined. The box dimension is used in this article, and the calculation formula is as follows:(14)Ds=−limr⟶0logNrlogr,where Ds is the self-similar fractal dimension of the destruction area, r is the side length of the square box, and Nr is the number of boxes required to cover the damaged area in the entire figure with a square box with side length r.The calculation of the fractal dimension of the shale compression test under different confining pressures is performed on the MATLAB software platform by writing a fractal dimension-related calculation program. The main steps of the algorithm are as follows: the first step is to process an image of acoustic emissions as a512×512 pixel image and read the processed image and then to preprocess the image grayscale, to process the image binary, and to store the relevant data. The second step is to cover the processed image with a square box with a side length of r, where r=2n,n=1,2,…,512, to count the number of square boxes in the failure area of the rock sample Nr, and to save the relevant data. Last, the double logarithmic coordinates of r and Nr are established for linear regression fitting, and the fractal dimension values under different confining pressures are calculated.Figure15 is the fitting curve of the fractal characteristics of the shale failure area when the stress level σ/σmax is 100% and the confining pressure is 10 MPa. The correlation coefficient R2=0.9968 can reflect the fact that the fractal dimension has high credibility, and the fractal dimension is Ds=1.9210.Figure 15
Fractal characteristics of the failure area when the stress level is 100% and the confining pressure is 10 MPa.Figures16–18 reflect various changes in stress levels, AE energy, and fractal dimensions under different confining pressures, respectively. Table 3 lists the values of the AE energy and the fractal dimension under six different confining pressures and different stress levels.Figure 16
Relationship between stress levels and AE Energy under different confining pressures.Figure 17
Relationship between stress levels and fractal dimensions under different confining pressures.Figure 18
Relationship between stress level, confining pressure, and fractal dimension.Table 3
AE energy and fractal dimension values under different stress levels and different confining pressures.
Stress levelConfining pressure(MPa)10%20%30%40%50%60%70%80%90%100%0AE energy00000.00010.00130.01120.05620.303540.3338Ds0000.47520.66170.83301.20601.18801.36701.51302AE energy0000.00010.0020.01960.15321.20846.028664.1970Ds0000.53290.7721.00201.26101.50001.59201.67704AE energy00000.00120.02950.24131.518419.9760108.7755Ds0000.41970.64551.00301.26401.5061.70101.78406AE energy00000.00200.03320.30252.153624.1428154.0042Ds0000.40760.63600.99521.28201.52301.7251.83508AE energy000000.0190.21061.797133.5573279.7562Ds00000.41970.88181.21101.50401.71501.893010AE energy00000.00060.01750.20691.761545.5702377.1735Ds00000.50120.82741.16801.48501.71301.9210The relationship between the stress level and the AE energy under different confining pressures is shown in Figure16. In the figure, the abscissa indicates different stress level values, and the ordinate indicates the energy value. When the stress level is less than 80%, the AE energy at each confining pressure is zero, and there is no trend. When the stress level is between 80% and 90%, the increase of the AE energy curve is relatively low, and the trend of change is relatively gentle and in a calm state. When the stress level exceeds 90% and reaches 100%, the change trend is more significant and is in a surge state. When the confining pressure is 10 MPa, the AE energy curve increases rapidly to the maximum, and when the confining pressure is 0 MPa, the AE energy value is the minimum. The results show that the energy released from the damage and fracture of the specimen is the largest when the confining pressure is 10 MPa. When the damage degree is the most severe, the internal damage is also the most serious, the fracture mode is relatively complex, and the confining pressure effect on the acoustic emission energy release is significant.From the relationship between the stress level and the fractal dimension under different confining pressures, as shown in Figure17, It can be seen that as the damage inside the specimen continues to accumulate, the fractal dimension value also increases with the increase of the confining pressure. When the stress level is less than 30%, the fractal dimension of each confining pressure state is 0. When the stress level reaches 30%, the fractal dimension values with confining pressures of 0 MPa, 2 MPa, 4 MPa, and 6 MPa change, and the curve change trend is relatively slow. When the confining pressure is 0 MPa, the fractal dimension value is the smallest, which is 1.5120, and the stress level of 8 MPa and 10 MPa reaches 40% and the fractal dimension value changes. However, the curve increases sharply, and the change trend is significant. When the confining pressure is 10 MPa, the fractal dimension reaches the maximum value of 1.9210. Further analysis shows that the fracture crack of the specimen with a confining pressure of 0 MPa is an inverted V-shaped, and the microfracture starts from the concentrated position of the quartz minerals in the upper part of the shale. At this time, the fractal dimension value is also the smallest. With increasing confining pressure, the fractal dimension increases gradually. When the confining pressure is 10 MPa, the fracture crack of the specimen is an inverted Z shape, the microcrack starts to fracture from the bottom, and the fractal dimension reaches the maximum value. The more complex the failure, the more the internal damage of the specimen. Therefore, when the fractal dimension value is larger, the fracture mode is more complicated, and the fractal dimension value gradually increases with crack propagation. The fractal dimension of the microstructure can be used as the quantitative parameter of the shale microcrack propagation.It can be clearly seen from Figure18 that the relationship between the confining pressure, stress level, and fractal dimension changes. Here, the fractal perspective is used to illustrate that the confining pressure effect has a significant influence on the shale failure mode.


## 5. Conclusion

(1)
Through NMR test analysis, polarized optical microscopy observation, and core X-ray whole-rock mineral diffraction analysis of the Niutitang Formation shale core, the results show that a large number of nanoscale pores and reservoirs are developed in the Niutitang Formation shale. The brittle minerals are relatively high in composition, mainly quartz, feldspar, pyrite, calcite, and so forth. These brittle minerals are easy to form cracks under the action of external force, which is more conducive to the exploitation of shale gas.(2)
The effect of the confining pressure on the strength and deformation of the shale samples from the lower Cambrian Niutitang Formation is significant. With increasing confining pressure, the compressive strength of the samples increases gradually, which is approximately linear. Regarding the elastic modulus, it shows the trend of first decreasing and then increasing. In the process of modelling, the displacement loading method is adopted, and the curve does not pass through the origin, which is quite different from the conventional triaxial test. The stress-strain curve of the specimen under loading can be roughly divided into three stages: the linear growth stage, postpeak softening stage, and gentle stage.(3)
The heterogeneity of the shale microstructure has a significant impact on the mechanical properties and fracture mode of shale. The failure modes of samples under different confining pressures are also quite different and can be summarized into three types: invertedV-shaped (0 MPa, 2 MPa, and 4 MPa), V-shaped (6 MPa and 8 MPa), and inverted Z-shaped (10 MPa). The three failure modes are caused by the confining pressure effect on the one hand and the uneven distribution and irregular geometry of quartz brittle minerals on the other hand.(4)
Due to the fractal characteristics of acoustic emission, the fractal dimension is greatly affected by the confining pressure effect and stress level. When the fractal dimension is larger, the fracture mode is more complex, and the internal damage of shale is more serious. The fractal dimension increases with the crack growth, so the fractal dimension of the microstructure can be used as the quantitative parameter of shale microcrack growth.(5)
In the later stage, the further research work is to consider the relationship and evolution mechanism between the micro and macro failure processes of shale under triaxial compression. According to the test results, the evolution model of fractal dimension, in spatial coordinates, is established. We do several qualitatively analyses of the fractal characteristics of each component in shale samples.



---

*Source: 8819711-2020-12-28.xml*