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.
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 where is the component of the stress tensor and 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 where is the strain tensor and 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 where is the shear modulus of rock; is the Young modulus of the rock; is the Poisson ratio of the rock; is the volumetric strain; is the pore water pressure; is the total stress tensor; is the Biot coefficient; is the Kronecker delta; and when , while when .
Taking equation (4) into equation (3) to obtain the equilibrium differential equation of rock microelements,
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
Fluid quality continuity equation is where is the seepage velocity; is the dynamic viscosity coefficient; is the permeability of the rock; is the fluid pressure; is the source origin; is the porosity; is the mass of fluid in microelement pores; is the bulk modulus of fluid; and is the bulk modulus of the solid constituent.
The governing equation of the seepage field is obtained by combining the above three equations
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 where is the permeability at zero stress and is the permeability under stress.
The porosity evolution equation describing the changes of rock porosity under constant temperature is
The common form is obtained by the approximate expansion of the exponential form where is Biot coefficient and is the initial porosity of rock mass.
Taking equation (11) into equation (9) to obtain the permeability coefficient evolution model,
In consideration of the influence of damage on the permeability of rock mass [16], the equation is obtained where 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: where and are the first and third principal stresses of rock mass; and are uniaxial tensile strength and compressive strength of rock; is the internal friction angle of rock; and and 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, indicates tensile failure of rock microelements, and means that the element continues to be loaded after failure, and means that the element has shear failure, and indicates that the microelement continues to be loaded after shear failure. where is the initial elastic modulus of rock mass; and are the tensile strain and compressive strain of rock mass; and and 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 where is the independent variable, representing the mechanical parameters of the rock element; is the average value of mechanical parameters; and 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 Figure 1, and the broken samples in the test are shown in Figure 2.
In order to verify the accuracy of rock damage statistical constitutive model, and are determined by processing test data. where is the peak strength of rock microelement; is the peak stress; is the peak strain; and represents the confining pressure.
According to the test data, the basic mechanical parameters and model parameters and are calculated in Table 1, and the theoretical curve of rock stress-strain under different confining pressures is obtained in Figure 3.
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 Figure 4. 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 is applied step by step to the edge of the hole, taking the original pore pressure 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.
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 point outside the circular hole on the rock plane model as the research object, the radial stress, circumferential stress, and shear stress are obtained. where is the distance from the point to the center of the circle and is the angle between the line between the point and the center of the circle and .
The rock parameters in the simulation model are shown in Table 2. 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 , and the right boundary is the minimum horizontal principal stress . Besides, is taken as the lateral pressure coefficient.
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 is 10 MPa, and the is 0.1. The number of damage points changes with the loading step as shown in Figure 5.
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 Figure 6. When the in situ stress is larger than , 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.
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 the axis is shown in Figure 7.
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 ( is equal to ) 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.
(a) 1MPa
(b) 2MPa
(c) 4MPa
(d) 6MPa
(e) 8MPa
(f) 10MPa
When the in situ stress is equal to , 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 Figure 9.
Specifically, is the numerical result, and 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 ( is unequal to ) 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.
(a) γ =0.1
(b) γ =0.2
(c) γ =0.4
(d) γ =0.6
(e) γ =0.8
(f) γ =1.0
When , 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 , the stress of the upper and lower boundary of circular hole is larger; when , the upper and lower boundary stress of circular hole increases; when , 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 and the maximum horizontal principal stress . 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 Figure 11. 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 has a certain error compared with the theoretical value . 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.
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
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (52004082), the Science and Technology Project of Henan Province (222102320381), the Foundation for Higher Education Key Research Project by Henan Province (22A130002), the Project of Henan Key Laboratory of Underground Engineering and Disaster Prevention (Henan Polytechnic University), the Strategic Consulting Research Project of Henan Institute for China Engineering Science and Technology Development Strategy (2021HENZT03), the Ph.D. Programs Foundation of Henan Polytechnic University (B2021-57), the Virtual Simulation Experiment Teaching Project of Henan Province, the Research and Practice Project of Educational and Teaching Reformation of Henan Polytechnic University (2021JG020 and 2019JG074), and the Postdoctoral Research Projects of Henan Province.