# Hydraulic Fracturing Mechanism of Rock Mass under Stress-Damage-Seepage Coupling Effect

Authors: Zhengzheng Cao; Yue Wang; Haixiao Lin; Qiang Sun; Xiaogang Wu; Xiaoshuai Yang
Journal: Geofluids
 (2022)

---

## Abstract

According to the damage evolution model of rock mass under stress-seepage coupling effect, the representative element theory is employed to describe the change law of rock mesostructure. Based on the theory of elasticity and Weibull distribution, the statistical damage constitutive model of rock mass and the finite element numerical algorithm are established, by adopting the COMSOL Multiphysics numerical software and MATLAB program. Besides, the validity of the statistical damage constitutive model of rock mass is verified by the triaxial compression test. Besides, the hydraulic fracturing processes of rock mass under equal and unequal in situ stresses are numerically simulated, and the mechanical behavior of rock mass during hydraulic fracturing in complex underground environment is also studied. Under the condition of equal in situ stress, the stress distribution of surrounding rock of circular hole is annular, which is similar to the elastic stress distribution of surrounding rock. Under the condition of unequal in situ stress, the stress distribution tends to be circular with the increase of lateral pressure coefficient, and the stress distribution along the diagonal decreases. The simulation results are in good agreement with the theoretical results, which indicates that the damage mechanical model and the numerical model have correlation and certain accuracy. By analyzing the size and direction of horizontal in situ stress, the shape and extension direction of cracks are judged, which provides an important theoretical basis for water inrush prediction and engineering protection.


---

## Body

## 1. Introduction

The interaction between in situ stress field and groundwater permeability exists in mining engineering, hydropower engineering, and underground tunnel construction engineering. All of them inevitably involve the initiation, expansion, and instability of rock mass cracks under the interaction of in situ stress and pore water pressure [1, 2]. In particular, the seepage instability of rock mass, such as dam instability and water inrush in tunnel surrounding rock, involves the hydraulic fracturing mechanism of rock mass under the stress-damage-seepage coupling effect [3, 4].Hydraulic fracturing is an important aspect of the stress-damage-seepage coupling problem. Many scholars at home and abroad have studied the multiphysical coupling effect of stress and seepage of rock mass damage by the methods of theoretical analysis, laboratory tests, and numerical simulation. For example, Zhu et al. [5] established an isotropic THM damage coupling model of rock, considering the influence of rock damage on temperature, seepage, and stress field distribution. Based on the finite element method, Bao et al. [6] introduced two sets of equations to describe mechanical problems such as rock deformation, fracture propagation, and fluid flow in hydraulic fracturing process. Liu et al. [7] established a mechanical model in the form of fracture, stress, and seepage coupling and further analyzed the basic process of hydraulic fracturing. Wang et al. [8] employed the Galerkin finite element method to simulate the influence of heterogeneity on rock in situ stress distribution, strain energy density, and fluid pressure under fluid-structure coupling. Based on the fluid-structure coupling theory, Li M. et al. [9] established the hydraulic fracturing model of extended permeability in the hydraulic fracturing mode. Li G. et al. [10] established a three-dimensional percolation-stress-damage coupling model that reflects the evolution process of rock mesodamage and simulated the expansion morphology of three-dimensional fractures. Sun et al. [11] conducted shale hydraulic fracturing experiments with a triaxial hydraulic fracturing experimental system and established a three-dimensional hydraulic fracturing calculation model of shale gas reservoir and studied the influence of bedding direction on hydraulic fracturing crack propagation in shale reservoir. Zhang et al. [12] carried out stress-seepage coupling test of sandstone without water and drainage by the automatic triaxial seepage test system and considered that the change of permeability before peak strength of sandstone was related to the development of fracture morphology.The relationship between stress and permeability has been analyzed in the above research results, and the influence of fracture network on the coupling effect has been considered. However, the damage evolution process of rock mass and its coupling effect are not considered systematically. In this paper, the damage variable of rock mass is introduced into the multiphysical coupling fields, by considering the mass conservation. Based on the established mechanical model of rock damage evolution under the coupling interaction of stress-seepage, combined with the theory of seepage mechanics, the evolution equations of porosity and permeability are deduced. Besides, the COMSOL Multiphysics is a finite element software for solving partial differential equations, which is widely employed to analyze the multiple physical fields. In this paper, the COMSOL Multiphysics is employed to simulate the cracking process of rock mass caused by water pressure, and the coupling effect of rock damage and seepage is further discussed.


## 2. The Governing Equation

Based on elastoplastic mechanics and fluid seepage theory, the rock damage evolution model under the coupling action of stress-seepage is established. The model satisfies certain conditional assumptions, and the stress field and seepage field are controlled by the formula equation. The damage is connected with parameters such as porosity and permeability, which builds a multiphysical field coupling model.


### 2.1. The Basic Assumptions

(1)
Rock mass is a kind of heterogeneous material composed of solid skeleton and pores. Besides, the microelement mechanical parameters obey Weibull random distribution, which satisfies the continuum mechanics theory(2)
Ideal elastoplastic damage occurs when rock mass is loaded. The constitutive relation of rock damage is elastoplastic, and the relevant mechanical parameters are functions of damage variables(3)
The seepage of fluid in rock mass follows Biot seepage theory and Darcy’s law(4)
The effect of temperature on the stress-damage-seepage coupling is not considered in this paper


### 2.2. Deformation Equation

It is assumed that each stress component on the rock microelement satisfies the statics equilibrium condition(1)σij,j+fi=0,where σij,j is the component of the stress tensor and fi is the component of the net body force.The rock microelement is subjected to elastic deformation by external force, and the strain and displacement satisfy the geometric equation(2)εij=12ui,j+uj,i,where εij is the strain tensor and ui,j is the displacement component.The microelement deformation of rock satisfies the generalized Hooke law [13]. Considering the pore water pressure, the modified effective stress formula is employed to derive the stress tensor of rock under the action of external force and fluid
(3)σij=2Gεij+2vG1−2vεVδij−αpwδij,where G=E/21+ is the shear modulus of rock; E is the Young modulus of the rock; v is the Poisson ratio of the rock; εV is the volumetric strain; pw is the pore water pressure; σij is the total stress tensor; α is the Biot coefficient; δij is the Kronecker delta; and δij=1 when i=j, while δij=0 when i≠j.
(4)εV=1−2vEσi+σj+σk.Taking equation (4) into equation (3) to obtain the equilibrium differential equation of rock microelements,
(5)Gui,j+uj,i+2vG1−2vεVδij+fi=αpwδij.


### 2.3. The Seepage Equation

According to the assumption that fluid flows at low speed in compressible porous media [14], the velocity is denoted by Darcy’s Law
(6)u=−kμ∇p.Fluid quality continuity equation is(7)∂ml∂t+∇ρu=Qm,∂ml∂t=n∂ρ∂t+ρ∂n∂t=ρnKw+1−nKs∂p∂t,where u is the seepage velocity; μ is the dynamic viscosity coefficient; k is the permeability of the rock; p is the fluid pressure; Qm is the source origin; n is the porosity; ml=ρn is the mass of fluid in microelement pores; Kw is the bulk modulus of fluid; and Ks is the bulk modulus of the solid constituent.The governing equation of the seepage field is obtained by combining the above three equations(8)nKw+1−nKs∂p∂t+∇⋅−kμ∇p=Qm.


### 2.4. Evolution Model of Porosity and Permeability

Based on the analysis of the coupling effect of rock stress and seepage, the stress field and seepage field are related by porosity and permeability [15]. The relationship between permeability coefficient and porosity is described by Kozeny-Carman cubic law
(9)kk0=nn03,where k0 is the permeability at zero stress and k is the permeability under stress.The porosity evolution equation describing the changes of rock porosity under constant temperature is(10)n=α−α−n0exp−εV−∇pKs.The common form is obtained by the approximate expansion of the exponential form(11)n=α−α−n01−∇p/Ks1+εV,where α is Biot coefficient and n0 is the initial porosity of rock mass.Taking equation (11) into equation (9) to obtain the permeability coefficient evolution model,
(12)k=k0α1+εV−α−n01−∇p/Ksn01+εV3.In consideration of the influence of damage on the permeability of rock mass [16], the equation is obtained
(13)k=k0nn03expαkD,where αk is the damage-permeability effect coefficient, which is 5.0.


### 2.5. Rock Damage Evolution Equation

Since the tensile strength of rock mass is far less than the compressive strength [17], it is necessary to use the theoretical judgment formula of maximum tensile stress (14) firstly to determine whether the microelement has tensile failure, and then, use the Mohr-Coulomb criterion (15) to determine whether the shear failure occurs in the microelement. Based on the test, the real approximation of rock under various stress states accurately reflects the strength conditions of rock mass, and the rock microelement strength is determined:
(14)F1=σ1−ft0,(15)F2=σ1−σ31+sinφ1−sinφ−fc0,where σ1 and σ3 are the first and third principal stresses of rock mass; ft0 and fc0 are uniaxial tensile strength and compressive strength of rock; φ is the internal friction angle of rock; and F1 and F2 are the threshold function of maximum tensile stress theory and Molar Coulomb criterion, respectively.The elastic modulus and damage variable of rock microelement are defined. In specific,F1=0 indicates tensile failure of rock microelements, and dF1>0 means that the element continues to be loaded after failure, and F2=0 means that the element has shear failure, and dF2>0 indicates that the microelement continues to be loaded after shear failure.
(16)E=1−DE0,D=0F1<0,F2<0,1−εt0ε12F1=0,dF1>0,1−εt0ε12F1=0,dF2>0,where E0 is the initial elastic modulus of rock mass; εt0 and εc0 are the tensile strain and compressive strain of rock mass; and ε1 and ε3 are the first and third principal strains of rock mass.


### 2.6. Heterogeneity of Rock Mass

In order to describe the heterogeneity of rock mass, it is necessary to assign rock micromechanical parameters, such as elastic modulus, strength, and permeability [18]. Assuming that the Weibull random distribution is obeyed, the probability density function is
(17)fx,z,m=mzxzm−1exp−xzm,x≥0,where x is the independent variable, representing the mechanical parameters of the rock element; z is the average value of mechanical parameters; and m is the coefficient of heterogeneity.


## 2.1. The Basic Assumptions

(1)
Rock mass is a kind of heterogeneous material composed of solid skeleton and pores. Besides, the microelement mechanical parameters obey Weibull random distribution, which satisfies the continuum mechanics theory(2)
Ideal elastoplastic damage occurs when rock mass is loaded. The constitutive relation of rock damage is elastoplastic, and the relevant mechanical parameters are functions of damage variables(3)
The seepage of fluid in rock mass follows Biot seepage theory and Darcy’s law(4)
The effect of temperature on the stress-damage-seepage coupling is not considered in this paper


## 2.2. Deformation Equation

It is assumed that each stress component on the rock microelement satisfies the statics equilibrium condition(1)σij,j+fi=0,where σij,j is the component of the stress tensor and fi is the component of the net body force.The rock microelement is subjected to elastic deformation by external force, and the strain and displacement satisfy the geometric equation(2)εij=12ui,j+uj,i,where εij is the strain tensor and ui,j is the displacement component.The microelement deformation of rock satisfies the generalized Hooke law [13]. Considering the pore water pressure, the modified effective stress formula is employed to derive the stress tensor of rock under the action of external force and fluid
(3)σij=2Gεij+2vG1−2vεVδij−αpwδij,where G=E/21+ is the shear modulus of rock; E is the Young modulus of the rock; v is the Poisson ratio of the rock; εV is the volumetric strain; pw is the pore water pressure; σij is the total stress tensor; α is the Biot coefficient; δij is the Kronecker delta; and δij=1 when i=j, while δij=0 when i≠j.
(4)εV=1−2vEσi+σj+σk.Taking equation (4) into equation (3) to obtain the equilibrium differential equation of rock microelements,
(5)Gui,j+uj,i+2vG1−2vεVδij+fi=αpwδij.


## 2.3. The Seepage Equation

According to the assumption that fluid flows at low speed in compressible porous media [14], the velocity is denoted by Darcy’s Law
(6)u=−kμ∇p.Fluid quality continuity equation is(7)∂ml∂t+∇ρu=Qm,∂ml∂t=n∂ρ∂t+ρ∂n∂t=ρnKw+1−nKs∂p∂t,where u is the seepage velocity; μ is the dynamic viscosity coefficient; k is the permeability of the rock; p is the fluid pressure; Qm is the source origin; n is the porosity; ml=ρn is the mass of fluid in microelement pores; Kw is the bulk modulus of fluid; and Ks is the bulk modulus of the solid constituent.The governing equation of the seepage field is obtained by combining the above three equations(8)nKw+1−nKs∂p∂t+∇⋅−kμ∇p=Qm.


## 2.4. Evolution Model of Porosity and Permeability

Based on the analysis of the coupling effect of rock stress and seepage, the stress field and seepage field are related by porosity and permeability [15]. The relationship between permeability coefficient and porosity is described by Kozeny-Carman cubic law
(9)kk0=nn03,where k0 is the permeability at zero stress and k is the permeability under stress.The porosity evolution equation describing the changes of rock porosity under constant temperature is(10)n=α−α−n0exp−εV−∇pKs.The common form is obtained by the approximate expansion of the exponential form(11)n=α−α−n01−∇p/Ks1+εV,where α is Biot coefficient and n0 is the initial porosity of rock mass.Taking equation (11) into equation (9) to obtain the permeability coefficient evolution model,
(12)k=k0α1+εV−α−n01−∇p/Ksn01+εV3.In consideration of the influence of damage on the permeability of rock mass [16], the equation is obtained
(13)k=k0nn03expαkD,where αk is the damage-permeability effect coefficient, which is 5.0.


## 2.5. Rock Damage Evolution Equation

Since the tensile strength of rock mass is far less than the compressive strength [17], it is necessary to use the theoretical judgment formula of maximum tensile stress (14) firstly to determine whether the microelement has tensile failure, and then, use the Mohr-Coulomb criterion (15) to determine whether the shear failure occurs in the microelement. Based on the test, the real approximation of rock under various stress states accurately reflects the strength conditions of rock mass, and the rock microelement strength is determined:
(14)F1=σ1−ft0,(15)F2=σ1−σ31+sinφ1−sinφ−fc0,where σ1 and σ3 are the first and third principal stresses of rock mass; ft0 and fc0 are uniaxial tensile strength and compressive strength of rock; φ is the internal friction angle of rock; and F1 and F2 are the threshold function of maximum tensile stress theory and Molar Coulomb criterion, respectively.The elastic modulus and damage variable of rock microelement are defined. In specific,F1=0 indicates tensile failure of rock microelements, and dF1>0 means that the element continues to be loaded after failure, and F2=0 means that the element has shear failure, and dF2>0 indicates that the microelement continues to be loaded after shear failure.
(16)E=1−DE0,D=0F1<0,F2<0,1−εt0ε12F1=0,dF1>0,1−εt0ε12F1=0,dF2>0,where E0 is the initial elastic modulus of rock mass; εt0 and εc0 are the tensile strain and compressive strain of rock mass; and ε1 and ε3 are the first and third principal strains of rock mass.


## 2.6. Heterogeneity of Rock Mass

In order to describe the heterogeneity of rock mass, it is necessary to assign rock micromechanical parameters, such as elastic modulus, strength, and permeability [18]. Assuming that the Weibull random distribution is obeyed, the probability density function is
(17)fx,z,m=mzxzm−1exp−xzm,x≥0,where x is the independent variable, representing the mechanical parameters of the rock element; z is the average value of mechanical parameters; and m is the coefficient of heterogeneity.


## 3. Verification of Damage Constitutive Model of Rock Mass

The rock mass for laboratory test is taken from a field project and processed into standard rock samples with a diameter of 50 mm and a height of 100 mm. The conventional triaxial compression test is carried out by the RMT-150 rock mechanics test system. The confining pressures of 5 MPa, 10 MPa, 15 MPa, 20 MPa, and 25 MPa are set for 5 samples, respectively, and gradually applied at a rate of 0.500 MPa/s. The vertical displacement is applied after the confining pressure stabilized. The conventional triaxial compression test is shown in Figure1, and the broken samples in the test are shown in Figure 2.Figure 1
Conventional triaxial compression test.Figure 2
Broken samples.In order to verify the accuracy of rock damage statistical constitutive model,m and z are determined by processing test data.
(18)∂σ∂ε1ε1=εb=0,m=1lnEε1b/σ1b−2μσ3,z=1m−1/m,where Fb is the peak strength of rock microelement; σ1b is the peak stress; ε1b is the peak strain; and σ3 represents the confining pressure.According to the test data, the basic mechanical parameters and model parametersm and z are calculated in Table 1, and the theoretical curve of rock stress-strain under different confining pressures is obtained in Figure 3.Table 1
Rock parameters of damage constitutive model.
Sampleσ3/MPaE/GPaμmF0/MPaP14.9823.860.251.20212.75P29.9824.140.252.08225.63P314.9823.100.241.76214.56P419.9722.050.231.78277.51P524.7335.440.281.91360.72Figure 3
Theoretical curves of rock stress-strain under different confining pressures.Regarding the peak strength, the test results are in consistent with the theoretical results. The test results show that the stress decreases rapidly after reaching the peak, while the stress in the theoretical curve decreases slowly after the peak; the postpeak stage of P1 and P5 is shorter, which is the same as the case of the rock yielding failure immediately after reaching the bearing limit in the test results. Therefore, the statistical constitutive model of rock damage established in this paper is reasonable and correct.


## 4. Numerical Simulation of Hydraulic Fracturing of Rock Mass

### 4.1. Establishment of Numerical Model

In order to study the cracking process of rock mass with a circular hole caused by different water pressure and in situ stress, a numerical calculation model of hydraulic fracturing of rock mass with a circular hole is established as shown in Figure4. Assuming that there is no obvious deformation along the length of rock, it is regarded as a plane strain problem [19]. The square rock model is 50 mm long on each side and has a circular hole with a radius of 3 mm in the center, and the water pressure pw is applied step by step to the edge of the hole, taking the original pore pressure p0 into account. According to the representative element theory, the Weibull random distribution function is employed to assign random numbers of strength and elastic modulus of rock mass.Figure 4
Numerical model of hydraulic fracturing of rock mass.By controlling the relevant variables, three kinds of rock damage under different water pressure, equal in situ stress, and unequal in situ stress are analyzed in this numerical simulation. According to the theory of elasticity, the original in situ stress and water pressure, taking a pointM outside the circular hole on the rock plane model as the research object, the radial stress, circumferential stress, and shear stress are obtained.
(19)σθ=σx+σy21−r2R2−σx−σy2⋅cos2θ1+3R4r4−pwr2R2,σr=σx+σy21−r2R2+σx−σy2⋅cos2θ1−4R2r2+3R4r4−pwr2R2,where R is the distance from the point M to the center of the circle and θ is the angle between the line between the point M and the center of the circle and σx.The rock parameters in the simulation model are shown in Table2. Considering the influence of in situ stress, roll support is set on the left and lower boundaries of the model. The upper boundary is the maximum horizontal principal stress σy, and the right boundary is the minimum horizontal principal stress σx. Besides, γ=σx/σy is taken as the lateral pressure coefficient.Table 2
The basic parameters.
The parameter nameValueUniaxial compressive strength of rock,fc100 MPaUniaxial tensile strength of rock,ft6.04 MPaThe elastic modulus of rock,E35 GPaThe initial porosity of rock,n00.05Poisson’s ratio of rock,v0.25Internal friction angle of rock,φ34°Permeability,k1×10−12m2Biot coefficient,α0.1Fluid density,ρw1020kg/m3Dynamic viscosity coefficient,μ0.0018Pa⋅sInitial pore pressure,p00 MPaDensity of rock,ρs2500kg/m3


### 4.2. Result Analysis

#### 4.2.1. Damage Evolution of Rock Mass under Different Pore Pressures

The damage condition in the rock mass is characterized by number change of microelement damage in rock mass, and the number of damage points under each loading step is obtained by statistical analysis. The pore pressure of 10 MPa, 20 MPa, and 30 MPa is applied at the boundary of circular hole, respectively, when the in situ stressσy is 10 MPa, and the γ is 0.1. The number of damage points changes with the loading step as shown in Figure 5.Figure 5
The number of damage points changes with the loading step.It is obvious that the rock damage occurs at the second step in three cases. The number of damage points shows a fluctuating and increasing trend on the whole, with the gradual increase of water pressure. However, the growth velocity is different in three cases, and the number of damage points increases faster when the water pressure is higher. When the water pressure increases, the stress on the boundary of the circular hole decreases and gradually approaches the tensile strength of the rock element. The greater the increase of water pressure at each step, the more obvious the effect of water pressure on rock fracture. Therefore, the pressure of circular hole is one of the important factors affecting rock damage.The rock damage conditions at the last step under different water pressures are shown in Figure6. When the in situ stress σy is larger than σx, the tension cracks occur firstly at the upper and lower boundaries in the circular hole. Besides, the greater the pressure applied, the longer and wider the crack propagation. When the water pressure is 30 MPa, the boundary of circular hole is damaged by the interaction of rock elements. The tension cracks occur, without compression shear crack, which is related with the in situ stress condition. Meanwhile, tensile cracks extend along the direction of the maximum principal stress in the far field, consistent with the theoretical result of crack occurrence, which verifies the accuracy of the hydraulic fracturing numerical model.Figure 6
The rock damage conditions at the last step. (a) 10MPa. (b) 20MPa. (c) 30MPa.Rock mass is continuously damaged by water pressure, followed by the damage in the internal structure, and the porosity is a physical quantity that characterizes the internal structure of rock materials. Taking BC section as the monitoring line, the porosity curve along thex axis is shown in Figure 7.Figure 7
The rock porosity in the monitoring line BC.The porosity increases with the increase of the distance from the center. The porosity changes obviously in the damage area of surrounding rock circular hole and then rises to a certain value, fluctuating up and down. And the greater the water pressure, the greater the porosity. The result indicates that porosity is a constantly changing physical quantity when water pressure is applied to rock mass, which is related to the deformation of rock skeleton caused by pore pressure.


#### 4.2.2. Hydraulic Fracturing Characteristics of Rock Mass under Equal In Situ Stress

In order to compare the theoretical solutions with simulated results, six equal confining pressures (σy is equal to σx) are set, with confining pressures of 1 MPa, 2 MPa, 4 MPa, 6 MPa, 8 MPa, and 10 MPa, respectively. The water pressure is applied to the boundary of the circular hole by increasing the water pressure with 1 MPa at each step. The Mises stress distribution under equal confining pressure is shown in Figure 8.Figure 8
The Mises stress distribution under equal confining pressure.
(a)
1MPa(b)
2MPa(c)
4MPa(d)
6MPa(e)
8MPa(f)
10MPaWhen the in situ stressσy is equal to σx, the failure point occurs around the circular hole in rock mass, and the maximum stress is in red and yellow areas, corresponding to the failure area of surrounding rock. Besides, the lower stress is in cyan area, corresponding to the plastic zone of surrounding rock, and the blue area is the elastic zone of surrounding rock. With the increase of confining pressure, the damage on the boundary of circular hole is gradually connected, and the failure zone and plastic zone tend to expand outwards.According to the number of damage points, the loading step of rock initial damage under different confining pressures is known, and then, the crack initiation pressure is also obtained. The linear relationship between crack initiation pressure and confining pressure is shown in Figure9.Figure 9
The linear relationship between crack initiation pressure and confining pressure.Specifically,y1 is the numerical result, and y2 is the theoretical result, and it is obvious that the two lines are quite closely parallel. The physical meaning of straight intercept is the tensile strength of rock mass, and the simulated tensile strength is 6.77 MPa. The relative error is 12%, compared with the original rock tensile strength parameter 6.04 MPa. Therefore, it is obvious that the numerical simulation is close to the result of theoretical analysis.


#### 4.2.3. Hydraulic Fracturing Characteristics of Rock Mass under Unequal In Situ Stress

The in situ stress field includes gravity stress, tectonic stress, pore fluid pressure, and thermal stress. In specific, the tectonic stress caused by the horizontal tectonic movement has the greatest influence on the in situ stress [20]. In order to analyze the influence of in situ stress further, five unequal confining pressures (σy is unequal to σx) are set, with lateral pressure coefficient γ of 0.1, 0.2, 0.4, 0.6, and 0.8, respectively. Besides, one equal confining pressure (lateral pressure coefficient is 1) is set. The Mises stress distribution under unequal confining pressure is obtained in Figure 10.Figure 10
The Mises stress distribution under unequal confining pressure.
(a)
γ =0.1(b)
γ =0.2(c)
γ =0.4(d)
γ =0.6(e)
γ =0.8(f)
γ =1.0Whenγ≤0.4, the stress distribution of circular hole rock shows along the diagonal, and the stress on the left and right boundary of the circular hole is larger; when γ=0.4, the stress of the upper and lower boundary of circular hole is larger; when γ≥0.6, the upper and lower boundary stress of circular hole increases; when γ=0.8, the stress appears obvious circular distribution. The stress distribution along the diagonal decreases with the increase of γ. On the whole, with the increase of the lateral pressure coefficient, the Mises stress tends to a circular distribution, which is the process of the reconciliation between the minimum horizontal principal stress σx and the maximum horizontal principal stress σy. In practical engineering, the size and direction of horizontal in situ stress are analyzed to judge the shape and extension direction of cracks, which provides an important basis for water inrush prediction and engineering protection.The change of crack initiation pressure with lateral pressure coefficient is shown in Figure11. The green line is the linear relationship derived from fitting the data, and the red line is the relationship derived from the formula. When the γ is small, the simulated value ps has a certain error compared with the theoretical value pt. With the increase of γ value, the error becomes smaller. On the whole, the influence of unequal confining pressures simulated in this section on rock hydraulic fracturing is consistent with the results of theoretical analysis, which indicates that the results of numerical simulation is reliable.Figure 11
The relationship between crack initiation pressure and lateral pressure coefficient.


## 4.1. Establishment of Numerical Model

In order to study the cracking process of rock mass with a circular hole caused by different water pressure and in situ stress, a numerical calculation model of hydraulic fracturing of rock mass with a circular hole is established as shown in Figure4. Assuming that there is no obvious deformation along the length of rock, it is regarded as a plane strain problem [19]. The square rock model is 50 mm long on each side and has a circular hole with a radius of 3 mm in the center, and the water pressure pw is applied step by step to the edge of the hole, taking the original pore pressure p0 into account. According to the representative element theory, the Weibull random distribution function is employed to assign random numbers of strength and elastic modulus of rock mass.Figure 4
Numerical model of hydraulic fracturing of rock mass.By controlling the relevant variables, three kinds of rock damage under different water pressure, equal in situ stress, and unequal in situ stress are analyzed in this numerical simulation. According to the theory of elasticity, the original in situ stress and water pressure, taking a pointM outside the circular hole on the rock plane model as the research object, the radial stress, circumferential stress, and shear stress are obtained.
(19)σθ=σx+σy21−r2R2−σx−σy2⋅cos2θ1+3R4r4−pwr2R2,σr=σx+σy21−r2R2+σx−σy2⋅cos2θ1−4R2r2+3R4r4−pwr2R2,where R is the distance from the point M to the center of the circle and θ is the angle between the line between the point M and the center of the circle and σx.The rock parameters in the simulation model are shown in Table2. Considering the influence of in situ stress, roll support is set on the left and lower boundaries of the model. The upper boundary is the maximum horizontal principal stress σy, and the right boundary is the minimum horizontal principal stress σx. Besides, γ=σx/σy is taken as the lateral pressure coefficient.Table 2
The basic parameters.
The parameter nameValueUniaxial compressive strength of rock,fc100 MPaUniaxial tensile strength of rock,ft6.04 MPaThe elastic modulus of rock,E35 GPaThe initial porosity of rock,n00.05Poisson’s ratio of rock,v0.25Internal friction angle of rock,φ34°Permeability,k1×10−12m2Biot coefficient,α0.1Fluid density,ρw1020kg/m3Dynamic viscosity coefficient,μ0.0018Pa⋅sInitial pore pressure,p00 MPaDensity of rock,ρs2500kg/m3


## 4.2. Result Analysis

### 4.2.1. Damage Evolution of Rock Mass under Different Pore Pressures

The damage condition in the rock mass is characterized by number change of microelement damage in rock mass, and the number of damage points under each loading step is obtained by statistical analysis. The pore pressure of 10 MPa, 20 MPa, and 30 MPa is applied at the boundary of circular hole, respectively, when the in situ stressσy is 10 MPa, and the γ is 0.1. The number of damage points changes with the loading step as shown in Figure 5.Figure 5
The number of damage points changes with the loading step.It is obvious that the rock damage occurs at the second step in three cases. The number of damage points shows a fluctuating and increasing trend on the whole, with the gradual increase of water pressure. However, the growth velocity is different in three cases, and the number of damage points increases faster when the water pressure is higher. When the water pressure increases, the stress on the boundary of the circular hole decreases and gradually approaches the tensile strength of the rock element. The greater the increase of water pressure at each step, the more obvious the effect of water pressure on rock fracture. Therefore, the pressure of circular hole is one of the important factors affecting rock damage.The rock damage conditions at the last step under different water pressures are shown in Figure6. When the in situ stress σy is larger than σx, the tension cracks occur firstly at the upper and lower boundaries in the circular hole. Besides, the greater the pressure applied, the longer and wider the crack propagation. When the water pressure is 30 MPa, the boundary of circular hole is damaged by the interaction of rock elements. The tension cracks occur, without compression shear crack, which is related with the in situ stress condition. Meanwhile, tensile cracks extend along the direction of the maximum principal stress in the far field, consistent with the theoretical result of crack occurrence, which verifies the accuracy of the hydraulic fracturing numerical model.Figure 6
The rock damage conditions at the last step. (a) 10MPa. (b) 20MPa. (c) 30MPa.Rock mass is continuously damaged by water pressure, followed by the damage in the internal structure, and the porosity is a physical quantity that characterizes the internal structure of rock materials. Taking BC section as the monitoring line, the porosity curve along thex axis is shown in Figure 7.Figure 7
The rock porosity in the monitoring line BC.The porosity increases with the increase of the distance from the center. The porosity changes obviously in the damage area of surrounding rock circular hole and then rises to a certain value, fluctuating up and down. And the greater the water pressure, the greater the porosity. The result indicates that porosity is a constantly changing physical quantity when water pressure is applied to rock mass, which is related to the deformation of rock skeleton caused by pore pressure.


### 4.2.2. Hydraulic Fracturing Characteristics of Rock Mass under Equal In Situ Stress

In order to compare the theoretical solutions with simulated results, six equal confining pressures (σy is equal to σx) are set, with confining pressures of 1 MPa, 2 MPa, 4 MPa, 6 MPa, 8 MPa, and 10 MPa, respectively. The water pressure is applied to the boundary of the circular hole by increasing the water pressure with 1 MPa at each step. The Mises stress distribution under equal confining pressure is shown in Figure 8.Figure 8
The Mises stress distribution under equal confining pressure.
(a)
1MPa(b)
2MPa(c)
4MPa(d)
6MPa(e)
8MPa(f)
10MPaWhen the in situ stressσy is equal to σx, the failure point occurs around the circular hole in rock mass, and the maximum stress is in red and yellow areas, corresponding to the failure area of surrounding rock. Besides, the lower stress is in cyan area, corresponding to the plastic zone of surrounding rock, and the blue area is the elastic zone of surrounding rock. With the increase of confining pressure, the damage on the boundary of circular hole is gradually connected, and the failure zone and plastic zone tend to expand outwards.According to the number of damage points, the loading step of rock initial damage under different confining pressures is known, and then, the crack initiation pressure is also obtained. The linear relationship between crack initiation pressure and confining pressure is shown in Figure9.Figure 9
The linear relationship between crack initiation pressure and confining pressure.Specifically,y1 is the numerical result, and y2 is the theoretical result, and it is obvious that the two lines are quite closely parallel. The physical meaning of straight intercept is the tensile strength of rock mass, and the simulated tensile strength is 6.77 MPa. The relative error is 12%, compared with the original rock tensile strength parameter 6.04 MPa. Therefore, it is obvious that the numerical simulation is close to the result of theoretical analysis.


### 4.2.3. Hydraulic Fracturing Characteristics of Rock Mass under Unequal In Situ Stress

The in situ stress field includes gravity stress, tectonic stress, pore fluid pressure, and thermal stress. In specific, the tectonic stress caused by the horizontal tectonic movement has the greatest influence on the in situ stress [20]. In order to analyze the influence of in situ stress further, five unequal confining pressures (σy is unequal to σx) are set, with lateral pressure coefficient γ of 0.1, 0.2, 0.4, 0.6, and 0.8, respectively. Besides, one equal confining pressure (lateral pressure coefficient is 1) is set. The Mises stress distribution under unequal confining pressure is obtained in Figure 10.Figure 10
The Mises stress distribution under unequal confining pressure.
(a)
γ =0.1(b)
γ =0.2(c)
γ =0.4(d)
γ =0.6(e)
γ =0.8(f)
γ =1.0Whenγ≤0.4, the stress distribution of circular hole rock shows along the diagonal, and the stress on the left and right boundary of the circular hole is larger; when γ=0.4, the stress of the upper and lower boundary of circular hole is larger; when γ≥0.6, the upper and lower boundary stress of circular hole increases; when γ=0.8, the stress appears obvious circular distribution. The stress distribution along the diagonal decreases with the increase of γ. On the whole, with the increase of the lateral pressure coefficient, the Mises stress tends to a circular distribution, which is the process of the reconciliation between the minimum horizontal principal stress σx and the maximum horizontal principal stress σy. In practical engineering, the size and direction of horizontal in situ stress are analyzed to judge the shape and extension direction of cracks, which provides an important basis for water inrush prediction and engineering protection.The change of crack initiation pressure with lateral pressure coefficient is shown in Figure11. The green line is the linear relationship derived from fitting the data, and the red line is the relationship derived from the formula. When the γ is small, the simulated value ps has a certain error compared with the theoretical value pt. With the increase of γ value, the error becomes smaller. On the whole, the influence of unequal confining pressures simulated in this section on rock hydraulic fracturing is consistent with the results of theoretical analysis, which indicates that the results of numerical simulation is reliable.Figure 11
The relationship between crack initiation pressure and lateral pressure coefficient.


## 4.2.1. Damage Evolution of Rock Mass under Different Pore Pressures

The damage condition in the rock mass is characterized by number change of microelement damage in rock mass, and the number of damage points under each loading step is obtained by statistical analysis. The pore pressure of 10 MPa, 20 MPa, and 30 MPa is applied at the boundary of circular hole, respectively, when the in situ stressσy is 10 MPa, and the γ is 0.1. The number of damage points changes with the loading step as shown in Figure 5.Figure 5
The number of damage points changes with the loading step.It is obvious that the rock damage occurs at the second step in three cases. The number of damage points shows a fluctuating and increasing trend on the whole, with the gradual increase of water pressure. However, the growth velocity is different in three cases, and the number of damage points increases faster when the water pressure is higher. When the water pressure increases, the stress on the boundary of the circular hole decreases and gradually approaches the tensile strength of the rock element. The greater the increase of water pressure at each step, the more obvious the effect of water pressure on rock fracture. Therefore, the pressure of circular hole is one of the important factors affecting rock damage.The rock damage conditions at the last step under different water pressures are shown in Figure6. When the in situ stress σy is larger than σx, the tension cracks occur firstly at the upper and lower boundaries in the circular hole. Besides, the greater the pressure applied, the longer and wider the crack propagation. When the water pressure is 30 MPa, the boundary of circular hole is damaged by the interaction of rock elements. The tension cracks occur, without compression shear crack, which is related with the in situ stress condition. Meanwhile, tensile cracks extend along the direction of the maximum principal stress in the far field, consistent with the theoretical result of crack occurrence, which verifies the accuracy of the hydraulic fracturing numerical model.Figure 6
The rock damage conditions at the last step. (a) 10MPa. (b) 20MPa. (c) 30MPa.Rock mass is continuously damaged by water pressure, followed by the damage in the internal structure, and the porosity is a physical quantity that characterizes the internal structure of rock materials. Taking BC section as the monitoring line, the porosity curve along thex axis is shown in Figure 7.Figure 7
The rock porosity in the monitoring line BC.The porosity increases with the increase of the distance from the center. The porosity changes obviously in the damage area of surrounding rock circular hole and then rises to a certain value, fluctuating up and down. And the greater the water pressure, the greater the porosity. The result indicates that porosity is a constantly changing physical quantity when water pressure is applied to rock mass, which is related to the deformation of rock skeleton caused by pore pressure.


## 4.2.2. Hydraulic Fracturing Characteristics of Rock Mass under Equal In Situ Stress

In order to compare the theoretical solutions with simulated results, six equal confining pressures (σy is equal to σx) are set, with confining pressures of 1 MPa, 2 MPa, 4 MPa, 6 MPa, 8 MPa, and 10 MPa, respectively. The water pressure is applied to the boundary of the circular hole by increasing the water pressure with 1 MPa at each step. The Mises stress distribution under equal confining pressure is shown in Figure 8.Figure 8
The Mises stress distribution under equal confining pressure.
(a)
1MPa(b)
2MPa(c)
4MPa(d)
6MPa(e)
8MPa(f)
10MPaWhen the in situ stressσy is equal to σx, the failure point occurs around the circular hole in rock mass, and the maximum stress is in red and yellow areas, corresponding to the failure area of surrounding rock. Besides, the lower stress is in cyan area, corresponding to the plastic zone of surrounding rock, and the blue area is the elastic zone of surrounding rock. With the increase of confining pressure, the damage on the boundary of circular hole is gradually connected, and the failure zone and plastic zone tend to expand outwards.According to the number of damage points, the loading step of rock initial damage under different confining pressures is known, and then, the crack initiation pressure is also obtained. The linear relationship between crack initiation pressure and confining pressure is shown in Figure9.Figure 9
The linear relationship between crack initiation pressure and confining pressure.Specifically,y1 is the numerical result, and y2 is the theoretical result, and it is obvious that the two lines are quite closely parallel. The physical meaning of straight intercept is the tensile strength of rock mass, and the simulated tensile strength is 6.77 MPa. The relative error is 12%, compared with the original rock tensile strength parameter 6.04 MPa. Therefore, it is obvious that the numerical simulation is close to the result of theoretical analysis.


## 4.2.3. Hydraulic Fracturing Characteristics of Rock Mass under Unequal In Situ Stress

The in situ stress field includes gravity stress, tectonic stress, pore fluid pressure, and thermal stress. In specific, the tectonic stress caused by the horizontal tectonic movement has the greatest influence on the in situ stress [20]. In order to analyze the influence of in situ stress further, five unequal confining pressures (σy is unequal to σx) are set, with lateral pressure coefficient γ of 0.1, 0.2, 0.4, 0.6, and 0.8, respectively. Besides, one equal confining pressure (lateral pressure coefficient is 1) is set. The Mises stress distribution under unequal confining pressure is obtained in Figure 10.Figure 10
The Mises stress distribution under unequal confining pressure.
(a)
γ =0.1(b)
γ =0.2(c)
γ =0.4(d)
γ =0.6(e)
γ =0.8(f)
γ =1.0Whenγ≤0.4, the stress distribution of circular hole rock shows along the diagonal, and the stress on the left and right boundary of the circular hole is larger; when γ=0.4, the stress of the upper and lower boundary of circular hole is larger; when γ≥0.6, the upper and lower boundary stress of circular hole increases; when γ=0.8, the stress appears obvious circular distribution. The stress distribution along the diagonal decreases with the increase of γ. On the whole, with the increase of the lateral pressure coefficient, the Mises stress tends to a circular distribution, which is the process of the reconciliation between the minimum horizontal principal stress σx and the maximum horizontal principal stress σy. In practical engineering, the size and direction of horizontal in situ stress are analyzed to judge the shape and extension direction of cracks, which provides an important basis for water inrush prediction and engineering protection.The change of crack initiation pressure with lateral pressure coefficient is shown in Figure11. The green line is the linear relationship derived from fitting the data, and the red line is the relationship derived from the formula. When the γ is small, the simulated value ps has a certain error compared with the theoretical value pt. With the increase of γ value, the error becomes smaller. On the whole, the influence of unequal confining pressures simulated in this section on rock hydraulic fracturing is consistent with the results of theoretical analysis, which indicates that the results of numerical simulation is reliable.Figure 11
The relationship between crack initiation pressure and lateral pressure coefficient.


## 5. Conclusions

According to the mechanical model of rock damage evolution under the stress and seepage coupling effect, a numerical model of rock mass under water pressure is established. By controlling the conditions of homogeneity, loading mode, confining pressure and water pressure, the mechanism, and process of rock damage evolution are deeply explored, and the following conclusions are obtained.(1)
The conventional triaxial compression test is carried out by the RMT-150 rock mechanics test system. The theoretical curve is consistent with the test results of the immediate yield failure of rock mass, when it reaches the bearing limit, which indicates that the statistical constitutive model of rock damage established in this paper is reasonable and correct(2)
Under the condition of equal in situ stress, the stress distribution of circular hole surrounding rock is annular, which is similar to the elastic stress distribution of surrounding rock. Under the condition of unequal in situ stress, the stress distribution tends to be circular with the increase of lateral pressure coefficient, and the stress distribution along the diagonal decreases(3)
The numerical simulation results are in good agreement with the theoretical results, which indicates that the rock mechanical damage model and the numerical model have correlation and are reasonable. In practical engineering, the size and direction of horizontal in situ stress are analyzed to judge the shape and extension direction of cracks, which provides an important basis for water inrush prediction and engineering protection



---

*Source: 5241708-2022-04-20.xml*