# Numerical Modeling on Hydraulic Fracturing in Coal-Rock Mass for Enhancing Gas Drainage

Authors: Zhigang Yuan; Yaohua Shao
Journal: Advances in Civil Engineering
 (2018)

---

## Abstract

The mechanism of how hydraulic fracturing influences gas drainage in coal-rock mass is still not clear due to its complex mechanism. In this work, statistical distributions are firstly introduced to describe heterogeneity of coal-rock mass; a novel simultaneously coupled mathematical model, which can describe the fully coupled process including seepage-damage coupling during hydraulic fracturing process and subsequent gas flow during gas drainage process, is established; its numerical implementation procedure is coded into a Matlab program to calculate the damage variables, and it partly uses COMSOL solver to obtain numerical solutions of governing equations with damage-flow coupling; the mathematical model and its implementation are validated for initial damage pressure and mode of a single solid model without considering flow-damage coupling, as well as fracture initiation pressure and influence of heterogeneity on damage evolution of hydraulic fracturing considering flow-damage coupling; and finally, based on an engineering practice of hydraulic fracturing with two boreholes, the mechanism of how hydraulic fracturing influences gas drainage is investigated, numerical simulation results indicate that coal-rock mass pore-fissure structure has been improved, and there would exist a gas migration channel with characteristics of higher porosity and lower stresses, which demonstrates significant effects and mechanism of hydraulic fracturing on improving coal-rock permeability and enhancing gas drainage. The research results provide a guide for operation of hydraulic fracturing and optimal layout of gas drainage boreholes.


---

## Body

## 1. Introduction

Coalbed gas as an unconventional natural gas gradually constitutes a vital part of China’s clean energy structure. Moreover, coalbed gas can also result in gas disaster that seriously affects coal mining safety [1]. For gas recovery and disaster control in underground coal mine, coal-rock mass permeability is an important parameter as it determines the practicability of gas extraction in advance of coal mining. However, most of the permeability of coal-rock mass excavated in China’s coal fields is relatively low, and as its excavation goes deeper, the permeability of coal-rock mass would become much lower, resulting in the difficulty of gas drainage before mining [2]. There are several measures of hydraulics used to improve permeability of coal-rock mass, such as hydraulic punching, water injection, hydraulic cutting, and hydraulic fracturing [3]. Particularly, the hydraulic fracturing measure is generally applied to enhance gas drainage in filed practices of China. However, the mechanism of how hydraulic fracturing affects the gas drainage process in coal-rock mass is still not understood fully due to its complex mechanical process. Related laboratory microscopic experiments denote that the evolution of hydraulic fractures (or cracks) is the damage evolution with initiation, extension of many microfissures due to the impact of external load and hydraulic pressure [3, 4]; additionally, gas drainage in damaged coal-rock mass after fracturing is a process with coal-rock mass deformation and gas flow coupling. Therefore, to investigate the mechanical mechanism of how hydraulic fracturing influences gas drainage process in coal-rock mass, the evolution of hydraulic fractures (or hydraulic cracks) in coal-rock mass and its subsequent coalbed gas seepage in coal-rock mass must be regarded as a simultaneous process with seepage-damage coupling during hydraulic fracturing and gas seepage during gas drainage. However, it is not easy to describe mechanisms of such complex process truly. With the development of numerical methods and modern computers, some models are raised to consider this coupled behavior of coal-rock mass during hydraulic fracturing [5, 6]. In general, there are two types of models used to simulate this process, namely, the direct method and indirect method. In the former method [7–9], the coal-rock mass is discretized into idealized microstructures (e.g., particles or springs) and the hydraulic cracks are denoted by the breakages among these microstructures. On the contrary, in the latter method, coal-rock mass is assumed to be a continuum and the hydraulic cracks are denoted by the nonreversible damage considered in constitutive law of coal-rock mass. Compared to the direct method, models using indirect method are relatively elaborate, as it can easily deal with complex medium containing dispersedly microcosmic cracks or pores and the multifield coupling problems to be consistent with coal-rock mass behavior [10, 11]. Representative capabilities of such models are composed of mesoscopic elements, which would degrade when specific strength criterions are satisfied [10, 11]. However, most existing models assume that the coal-rock mass is not permeable, and the influence of mechanical-hydraulic coupling on characteristic of hydraulic fracturing are usually ignored due to its complex description [11, 12]. Hydromechanical coupling is the essence of hydraulic fracturing, which is the key to describe the hydraulic fracturing process accurately. Zhu et al. [11, 12] developed a RFPA-flow program, in which a model with damage-flow coupling for heterogeneous material was embedded using finite element method. Yuan and Harrison [6] proposed the local degradation method, in which the hydromechanical coupling with evolutions of element’s permeability and dilatancy was considered. More recently, Lyakhovsky and Hamiel [13] simulated the coupling process of crack evolution and fluid flow on the basis of a thermodynamic mathematical model. Lu et al. [14] raised a double-scale mathematical model to describe the actual coal-rock mass comprising microcracks, and the process of hydraulic fracturing with hydromechanical coupling was simulated. In summary, models based on the indirect method have been well-developed and applied by many scholars to study the process of hydraulic fracturing, but very few studies can be found considering the formation of hydraulic fractures and gas seepage within coal-rock mass simultaneously. In the present work, to study the mechanism of how hydraulic fracturing influences gas drainage, inspired from the existing models, a novel simultaneously coupled mathematical model to describe damage of coal-rock mass and gas seepage is built, which can realistically describe the fully coupled process including damage-seepage coupling during hydraulic fracturing process and subsequent gas flow during gas drainage process in damaged coal-rock mass. In the following, the method for heterogeneity assignment of material properties is proposed firstly; then governing equations for the simultaneously coupled mathematical model are introduced secondly; the procedure of numerical implementation is proposed in Section 4; then the mathematical model and numerical procedure are validated and then applied to an engineering case; and lastly, Section 6 lists the conclusions obtained. The research results provide a guide for the operation of hydraulic fracturing and arrangement of gas drainage boreholes.


## 2. Assignment of Material Properties

Coal-rock mass consists of various mineral particles and cementing materials. The dispersedly microcosmic cracks and pores in coal-rock mass exhibit heterogeneity of hydraulic and mechanical properties. To describe such heterogeneity characteristic of coal-rock mass, suppose that it consists of numerous elements, whose mechanical properties satisfy the Weibull distribution with a probability density function (PDF) defined as follows [11]:(1)fz,ω,m=0,z<0,zωm−1mωexp−zωm,z≥0,where z denotes the mechanical property of the element (e.g., the elastic modulus), ω denotes the mean of mechanical properties of all elements, and m denotes homogeneity index for mechanical property of element. On basis of characteristic of Weibull distribution, the mechanical properties of all elements would be closer to ω as homogeneity index m increases, and thus producing a more homogeneous calculation model for coal-rock mass. Preliminary studies show that when m ≥ 1000, the coal-rock mass could be assumed to be homogeneous.Moreover, to describe the heterogeneity of element’s hydraulic property (such as permeability), a log-normal distribution function is introduced [14]:(2)Rξ,σ,μ=12πξσexp−lnμ−ξ22σ2,where ξ represents the hydraulic property of the element, μ and σ represent the average values of hydraulic properties of elements and corresponding deviation, respectively. The larger value of σ, the more the inhomogeneous calculation model for coal-rock mass.According to equations (1) and (2), a Monte Carlo simulation can be used to numerically produce a heterogeneous coal-rock mass model with heterogeneous hydraulic and mechanical properties [14]. The numerically generated heterogeneous model is similar to the realistic coal-rock mass specimen. Therefore, how the heterogeneity of coal-rock mass property affects the mechanical process could be investigated using a numerical simulation method.


## 3. Simultaneously Coupled Mathematical Model

Coal-rock mass is a typical porous medium and has an obvious nonlinear characteristic as plastic deformation occurs. Nevertheless, the elastoplastic model would increase cost and difficulty of calculation with low numerical accuracy. Taking such into account, a linear elastic model combined with damage theory is adopted instead [15]. Two simultaneously coupled processes are involved in the present study, one is the seepage-damage coupling process during hydraulic fracturing, and the other is the gas flow process during gas drainage after fracturing. The former process has a significant influence on the initialization and evolution of the latter. This section briefly summarizes the related equations of simultaneously coupled mathematical model, which are built at element scale.


### 3.1. Governing Equations with Damage-Seepage Coupling during Hydraulic Fracturing

#### 3.1.1. Governing Equations for Deformation Field

Governing equations for deformation field are composed of constitutive relation, displacement-strain relation, and equilibration equation. The indicial notation is adopted to express tensor quantities and vector hereafter. For porous coal-rock mass with single-phase saturated, the effective stressσ′ij (positive for tension), total stress σij, and pore pressure p satisfy the modified principle of Terzaghi effective stress [16].(3)σij′=αpδij+σij.The porous coal-rock mass is assumed to be elastic before damage initiates, then its constitutive relationship conforms to Hooke’s law under isothermal state [17].(4)2νG1−2νεkkδij+2εijG=αpδij+σij,where α and δij are called Biot coefficient and Kronecker delta; G and v are shear modulus and Poisson’s ratio of coal-rock mass; εij is infinitesimal-strain tensor, and εkk = ε11 + ε22 + ε33 with ε11, ε22, ε33 being the principal strain in three directions.The displacement-strain relation for coal-rock mass under infinitesimal strain condition is(5)εij=uj,i+ui,j2,where ui represents displacement in ith direction and the notation [,j] denotes differentiation with respect to the coordinate xj.The equilibration equations for coal-rock mass can be expressed as(6)fi+σji,j=0,where fi represents the volumetric force of ith orientation.On the basis of equations (4)–(6), mathematical governing equations for deformation field can be obtained:(7)fi−G2v−1uj,ji+Gui,jj−αp,i=0.As seen from equation (7), the influence of fluid pressure on mechanical equilibrium of coal-rock mass has been considered, namely, it reflects the unidirectional coupling effect between fluid and solid.


#### 3.1.2. Governing Equations for Water Seepage Field

The governing equations of water flow in coal-rock mass are composed of water balance equation (continuity equation) and water seepage equation (momentum equation) [17, 18]. Supposing the total bulk volume of coal-rock mass is V, and V = Ve + Vs (Ve represents volume of pore, and Vs represents volume of solid). According to water balance principle, mathematical equation for water balance can be expressed as follows:(8)1V∂V∂t=∂εv∂t=1V∂Vs∂t+∂Ve∂t,where t represents time and εv denotes volumetric strain.The coal-rock mass’s porosity,ϕ, can be given by ϕ=Ve/V. Based on equation (4), equations (9) and (10) are deduced [17]:(9)−ϕβl∂p∂t−∇ql=1V∂Ve∂t,where βl represents the water bulk modulus and ql represents the vector of water seepage velocity; the first and second terms of equation (9) are the change of water quantity due to water pressure and water flux flowing out, respectively.(10)ϕ−1Ks∂p∂t+13Ksδij∂σij′∂t=1V∂Vs∂t,where Ks represents bulk modulus for coal-rock mass matrix; the first and second terms of equation (10) are volume change of solid due to change of water pressure and effective stress, respectively.Substituting equations (9) and (10) in equation (8), one can gain the mathematical governing equations for water seepage field:(11)1V∂V∂t=∂εv∂t=ϕ−1Ks−ϕβl∂p∂t−∇ql+13Ksδij∂σij′∂t.Assuming that the water flow process in coal-rock mass satisfies the Darcy seepage law (water gravity effect is ignored here) [17–19],(12)ql=−kμw⋅∇p,where μw denotes dynamic viscosity coefficient for water and k denotes the coal-rock mass permeability.According to equations (11) and (12), the final mathematical governing equations for seepage field is given by(13)α∂εv∂t−∇⋅kμw∇p−ϕ−1Ks−ϕβl∂p∂t=0.


#### 3.1.3. Evolutions of Damage and Permeability

As mentioned above, the elastic damage constitutive relation is utilized to depict the coal-rock mass behavior in the present study. In other words, the strain-stress curve for coal-rock mass is linear unless the threshold of damage is reached. Two damage modes for coal-rock mass are considered in the numerical procedure. If the shear stress in coal-rock mass meets Mohr–Coulomb criterion (referred as M-C hereafter), it would damage in shear mode, while the coal-rock mass would damage in form of tensile when the maximum principal stress meets maximum tensile stress criterion (referred as M-T hereafter), as illustrated in Figure1. The tensile mode is also prior to the shear mode, namely, the M-T criterion is used first to judge whether the coal-rock mass is damaged in tension or not, and only coal-rock mass that is not damaged in tension would be judged for shear mode by M-C criterion. It is worth noting that equations (14)–(22) operate on effective stress only.Figure 1
Curves of elastic damage constitutive relations under uniaxial states.The M-T criterion and M-C criterion could be written as [20, 21](14)F1=−ft+σ1,(15)F2=−σ3−sinφ+1sinφ−1σ1−fc,where φ represents the internal frictional angle of coal-rock mass, fc and ft represent compressive strength and tensile strength of coal-rock mass, respectively, while F2 and F1 are threshold functions of damage for shear and tensile mode, respectively.As the damage of coal-rock mass evolves, the mechanical index such as elastic modulus would decrease monotonically according to theory of elastic damage [15, 21].(16)E=E0⋅1−D,where E0 is the original elastic modulus for coal-rock mass without damage and D is called the damage variable, whose value is 1.0 for coal-rock mass with full damage and 0.0 for the coal-rock mass without damage.If coal-rock mass is in the uniaxial stress state, the curves of constitutive relationships are shown in Figure1. When F1 > 0, namely, the coal-rock mass is under tensile state, the value of D can be calculated by(17)D=−εt0εn+1,ε≥εt0,0,ε<εt0,where εt0 denotes limit strain before tensile damage occurs.Hooke’s law is given as [17].(18)σi=εi−v2v−1εjjEv+1.According to M-T criterion and equation (18), the following relation can be obtained with η=−ν/v+12v−1.(19)E=ftεv⋅η+ε1/v+1.Based on equations (16) and (19), the expression of D under triaxial stress condition is rewritten as(20)D=−ftηE0εv+E0ε1/v+1+1.Corresponding to its damage evolution under tensile condition, whenF2 > 0, namely, the coal-rock mass is damaged in shear mode as illustrated in Figure 1, the expression of D is calculated by(21)D=1−εc0εn,ε≤εc0,0,ε>εc0.When the coal-rock mass is under triaxial stress condition, the expression ofD can be deduced from equations (15) and (18)(22)D=1−fcE0εvη+ε3/v+1−κεvη+ε1/v+1.When coal-rock mass experiences dilatancy due to its strength or stiffness degradation as damage evolves, its hydraulic property would also change. The porosity,ϕ, strongly associated with rock’s stress state, could be given by [22](23)ϕ=ϕr+ϕ0−ϕrexpσ¯e⋅αϕ,where ϕr denotes the residual porosity under extremely large stress condition, ϕ0 denotes the initial porosity under zero-stress condition, αϕ denotes the coefficient of stress sensitivity, and σ¯e is the average of effective stress calculated as σ¯e=αp+σ3+σ2+σ1/3 with σ3, σ2, and σ1 being the third, the second, and the first total principal stress in three directions.Moreover, as coal-rock mass damage accumulates, its permeability would also change accordingly. According to the related studies of relation between porosity and permeability [22], the damage evolution of permeability k can be described using an exponential function as follows:(24)k=k0ϕϕ0βexpD⋅αD,where k0 represents the initial permeability in zero-stress state, αD is the damage-permeability effect coefficient of 5.0, and β is the index coefficient of 3. The value of D in equation (24) can be calculated by equations (20) and (22). Through some numerical simulation examples given below, equation (24) can well describe the permeability damage evolution of coal-rock mass.


### 3.2. Governing Equation for Gas Flow during Gas Drainage

To describe the gas drainage process including coal-rock mass deformation (the same as that already introduced in Section3.1.1) and gas flow after fracturing, except for the previous hypothesis above, the established governing equations used to describe the gas flow process are based on the following hypotheses: (1) the coal-rock mass is completely saturated by single gas; (2) gas existed in pores is the ideal gas, whose viscosity is kept constant; and (3) the rule of gas flow in coal-rock mass satisfies Darcy’s law and is deemed as an isothermal process.The governing equation of gas mass conservation under isothermal state is given by [23–26](25)Qy=∇⋅ρgqg+∂C∂t,where Qy represents the gas sink or source and ρg is the gas density, which is associated with gas pressure as [26, 27].(26)ρg=MgpTR=pβ.q
g is the vector of Darcy velocity for free-phase gas, which is expressed using Darcy’s law with its gravity effect neglected:(27)qg=−kμg∇p.C is the gas content retained in coal-rock mass, which is given by Langmuir’s formula [27, 28]:(28)C=βa2a1ρsa2p+1+ϕp0p2.In equations (26)–(28), Mg represents gas molecular mass, T represents the Kelvin temperature, p represents gas pressure, R represents ideal gas constant, β represents compressibility coefficient of gas, μg represents gas kinematic viscosity coefficient, k is the coal-rock mass permeability; a2 and a1 are called Langmuir’s coefficients with their units being Pa−1 and m3·kg−1, respectively, ρs is the coal-rock mass density, and p0 represents the standard atmospheric pressure.The ultimate governing equation of gas flow is derived by substituting equations (26)–(28) in equation (25):(29)2pβϕp0−ρspa22a12pa2+12+a2a1ρspa2+1∂p∂t−∇⋅2pβkμg∇p=Qy.In brief, two damage-flow coupling processes are included in the simultaneously coupled mathematical model. equations (7) and (13) are the governing equations with damage-flow coupling during hydraulic fracturing, and equations (29) and (7) are governing equations with gas flow and coal-rock mass deformation coupling during gas drainage. The influence of fluids (including gas and water) seepage on the deformation behavior is considered in equation (7). Besides, the porosity and permeability, strongly associated with its stress state, are considered in equations (23) and (24). Furthermore, the influence of deformation field on the flow behaviors of water and gas are implicit in equations (13) and (29), respectively.


## 3.1. Governing Equations with Damage-Seepage Coupling during Hydraulic Fracturing

### 3.1.1. Governing Equations for Deformation Field

Governing equations for deformation field are composed of constitutive relation, displacement-strain relation, and equilibration equation. The indicial notation is adopted to express tensor quantities and vector hereafter. For porous coal-rock mass with single-phase saturated, the effective stressσ′ij (positive for tension), total stress σij, and pore pressure p satisfy the modified principle of Terzaghi effective stress [16].(3)σij′=αpδij+σij.The porous coal-rock mass is assumed to be elastic before damage initiates, then its constitutive relationship conforms to Hooke’s law under isothermal state [17].(4)2νG1−2νεkkδij+2εijG=αpδij+σij,where α and δij are called Biot coefficient and Kronecker delta; G and v are shear modulus and Poisson’s ratio of coal-rock mass; εij is infinitesimal-strain tensor, and εkk = ε11 + ε22 + ε33 with ε11, ε22, ε33 being the principal strain in three directions.The displacement-strain relation for coal-rock mass under infinitesimal strain condition is(5)εij=uj,i+ui,j2,where ui represents displacement in ith direction and the notation [,j] denotes differentiation with respect to the coordinate xj.The equilibration equations for coal-rock mass can be expressed as(6)fi+σji,j=0,where fi represents the volumetric force of ith orientation.On the basis of equations (4)–(6), mathematical governing equations for deformation field can be obtained:(7)fi−G2v−1uj,ji+Gui,jj−αp,i=0.As seen from equation (7), the influence of fluid pressure on mechanical equilibrium of coal-rock mass has been considered, namely, it reflects the unidirectional coupling effect between fluid and solid.


### 3.1.2. Governing Equations for Water Seepage Field

The governing equations of water flow in coal-rock mass are composed of water balance equation (continuity equation) and water seepage equation (momentum equation) [17, 18]. Supposing the total bulk volume of coal-rock mass is V, and V = Ve + Vs (Ve represents volume of pore, and Vs represents volume of solid). According to water balance principle, mathematical equation for water balance can be expressed as follows:(8)1V∂V∂t=∂εv∂t=1V∂Vs∂t+∂Ve∂t,where t represents time and εv denotes volumetric strain.The coal-rock mass’s porosity,ϕ, can be given by ϕ=Ve/V. Based on equation (4), equations (9) and (10) are deduced [17]:(9)−ϕβl∂p∂t−∇ql=1V∂Ve∂t,where βl represents the water bulk modulus and ql represents the vector of water seepage velocity; the first and second terms of equation (9) are the change of water quantity due to water pressure and water flux flowing out, respectively.(10)ϕ−1Ks∂p∂t+13Ksδij∂σij′∂t=1V∂Vs∂t,where Ks represents bulk modulus for coal-rock mass matrix; the first and second terms of equation (10) are volume change of solid due to change of water pressure and effective stress, respectively.Substituting equations (9) and (10) in equation (8), one can gain the mathematical governing equations for water seepage field:(11)1V∂V∂t=∂εv∂t=ϕ−1Ks−ϕβl∂p∂t−∇ql+13Ksδij∂σij′∂t.Assuming that the water flow process in coal-rock mass satisfies the Darcy seepage law (water gravity effect is ignored here) [17–19],(12)ql=−kμw⋅∇p,where μw denotes dynamic viscosity coefficient for water and k denotes the coal-rock mass permeability.According to equations (11) and (12), the final mathematical governing equations for seepage field is given by(13)α∂εv∂t−∇⋅kμw∇p−ϕ−1Ks−ϕβl∂p∂t=0.


### 3.1.3. Evolutions of Damage and Permeability

As mentioned above, the elastic damage constitutive relation is utilized to depict the coal-rock mass behavior in the present study. In other words, the strain-stress curve for coal-rock mass is linear unless the threshold of damage is reached. Two damage modes for coal-rock mass are considered in the numerical procedure. If the shear stress in coal-rock mass meets Mohr–Coulomb criterion (referred as M-C hereafter), it would damage in shear mode, while the coal-rock mass would damage in form of tensile when the maximum principal stress meets maximum tensile stress criterion (referred as M-T hereafter), as illustrated in Figure1. The tensile mode is also prior to the shear mode, namely, the M-T criterion is used first to judge whether the coal-rock mass is damaged in tension or not, and only coal-rock mass that is not damaged in tension would be judged for shear mode by M-C criterion. It is worth noting that equations (14)–(22) operate on effective stress only.Figure 1
Curves of elastic damage constitutive relations under uniaxial states.The M-T criterion and M-C criterion could be written as [20, 21](14)F1=−ft+σ1,(15)F2=−σ3−sinφ+1sinφ−1σ1−fc,where φ represents the internal frictional angle of coal-rock mass, fc and ft represent compressive strength and tensile strength of coal-rock mass, respectively, while F2 and F1 are threshold functions of damage for shear and tensile mode, respectively.As the damage of coal-rock mass evolves, the mechanical index such as elastic modulus would decrease monotonically according to theory of elastic damage [15, 21].(16)E=E0⋅1−D,where E0 is the original elastic modulus for coal-rock mass without damage and D is called the damage variable, whose value is 1.0 for coal-rock mass with full damage and 0.0 for the coal-rock mass without damage.If coal-rock mass is in the uniaxial stress state, the curves of constitutive relationships are shown in Figure1. When F1 > 0, namely, the coal-rock mass is under tensile state, the value of D can be calculated by(17)D=−εt0εn+1,ε≥εt0,0,ε<εt0,where εt0 denotes limit strain before tensile damage occurs.Hooke’s law is given as [17].(18)σi=εi−v2v−1εjjEv+1.According to M-T criterion and equation (18), the following relation can be obtained with η=−ν/v+12v−1.(19)E=ftεv⋅η+ε1/v+1.Based on equations (16) and (19), the expression of D under triaxial stress condition is rewritten as(20)D=−ftηE0εv+E0ε1/v+1+1.Corresponding to its damage evolution under tensile condition, whenF2 > 0, namely, the coal-rock mass is damaged in shear mode as illustrated in Figure 1, the expression of D is calculated by(21)D=1−εc0εn,ε≤εc0,0,ε>εc0.When the coal-rock mass is under triaxial stress condition, the expression ofD can be deduced from equations (15) and (18)(22)D=1−fcE0εvη+ε3/v+1−κεvη+ε1/v+1.When coal-rock mass experiences dilatancy due to its strength or stiffness degradation as damage evolves, its hydraulic property would also change. The porosity,ϕ, strongly associated with rock’s stress state, could be given by [22](23)ϕ=ϕr+ϕ0−ϕrexpσ¯e⋅αϕ,where ϕr denotes the residual porosity under extremely large stress condition, ϕ0 denotes the initial porosity under zero-stress condition, αϕ denotes the coefficient of stress sensitivity, and σ¯e is the average of effective stress calculated as σ¯e=αp+σ3+σ2+σ1/3 with σ3, σ2, and σ1 being the third, the second, and the first total principal stress in three directions.Moreover, as coal-rock mass damage accumulates, its permeability would also change accordingly. According to the related studies of relation between porosity and permeability [22], the damage evolution of permeability k can be described using an exponential function as follows:(24)k=k0ϕϕ0βexpD⋅αD,where k0 represents the initial permeability in zero-stress state, αD is the damage-permeability effect coefficient of 5.0, and β is the index coefficient of 3. The value of D in equation (24) can be calculated by equations (20) and (22). Through some numerical simulation examples given below, equation (24) can well describe the permeability damage evolution of coal-rock mass.


## 3.1.1. Governing Equations for Deformation Field

Governing equations for deformation field are composed of constitutive relation, displacement-strain relation, and equilibration equation. The indicial notation is adopted to express tensor quantities and vector hereafter. For porous coal-rock mass with single-phase saturated, the effective stressσ′ij (positive for tension), total stress σij, and pore pressure p satisfy the modified principle of Terzaghi effective stress [16].(3)σij′=αpδij+σij.The porous coal-rock mass is assumed to be elastic before damage initiates, then its constitutive relationship conforms to Hooke’s law under isothermal state [17].(4)2νG1−2νεkkδij+2εijG=αpδij+σij,where α and δij are called Biot coefficient and Kronecker delta; G and v are shear modulus and Poisson’s ratio of coal-rock mass; εij is infinitesimal-strain tensor, and εkk = ε11 + ε22 + ε33 with ε11, ε22, ε33 being the principal strain in three directions.The displacement-strain relation for coal-rock mass under infinitesimal strain condition is(5)εij=uj,i+ui,j2,where ui represents displacement in ith direction and the notation [,j] denotes differentiation with respect to the coordinate xj.The equilibration equations for coal-rock mass can be expressed as(6)fi+σji,j=0,where fi represents the volumetric force of ith orientation.On the basis of equations (4)–(6), mathematical governing equations for deformation field can be obtained:(7)fi−G2v−1uj,ji+Gui,jj−αp,i=0.As seen from equation (7), the influence of fluid pressure on mechanical equilibrium of coal-rock mass has been considered, namely, it reflects the unidirectional coupling effect between fluid and solid.


## 3.1.2. Governing Equations for Water Seepage Field

The governing equations of water flow in coal-rock mass are composed of water balance equation (continuity equation) and water seepage equation (momentum equation) [17, 18]. Supposing the total bulk volume of coal-rock mass is V, and V = Ve + Vs (Ve represents volume of pore, and Vs represents volume of solid). According to water balance principle, mathematical equation for water balance can be expressed as follows:(8)1V∂V∂t=∂εv∂t=1V∂Vs∂t+∂Ve∂t,where t represents time and εv denotes volumetric strain.The coal-rock mass’s porosity,ϕ, can be given by ϕ=Ve/V. Based on equation (4), equations (9) and (10) are deduced [17]:(9)−ϕβl∂p∂t−∇ql=1V∂Ve∂t,where βl represents the water bulk modulus and ql represents the vector of water seepage velocity; the first and second terms of equation (9) are the change of water quantity due to water pressure and water flux flowing out, respectively.(10)ϕ−1Ks∂p∂t+13Ksδij∂σij′∂t=1V∂Vs∂t,where Ks represents bulk modulus for coal-rock mass matrix; the first and second terms of equation (10) are volume change of solid due to change of water pressure and effective stress, respectively.Substituting equations (9) and (10) in equation (8), one can gain the mathematical governing equations for water seepage field:(11)1V∂V∂t=∂εv∂t=ϕ−1Ks−ϕβl∂p∂t−∇ql+13Ksδij∂σij′∂t.Assuming that the water flow process in coal-rock mass satisfies the Darcy seepage law (water gravity effect is ignored here) [17–19],(12)ql=−kμw⋅∇p,where μw denotes dynamic viscosity coefficient for water and k denotes the coal-rock mass permeability.According to equations (11) and (12), the final mathematical governing equations for seepage field is given by(13)α∂εv∂t−∇⋅kμw∇p−ϕ−1Ks−ϕβl∂p∂t=0.


## 3.1.3. Evolutions of Damage and Permeability

As mentioned above, the elastic damage constitutive relation is utilized to depict the coal-rock mass behavior in the present study. In other words, the strain-stress curve for coal-rock mass is linear unless the threshold of damage is reached. Two damage modes for coal-rock mass are considered in the numerical procedure. If the shear stress in coal-rock mass meets Mohr–Coulomb criterion (referred as M-C hereafter), it would damage in shear mode, while the coal-rock mass would damage in form of tensile when the maximum principal stress meets maximum tensile stress criterion (referred as M-T hereafter), as illustrated in Figure1. The tensile mode is also prior to the shear mode, namely, the M-T criterion is used first to judge whether the coal-rock mass is damaged in tension or not, and only coal-rock mass that is not damaged in tension would be judged for shear mode by M-C criterion. It is worth noting that equations (14)–(22) operate on effective stress only.Figure 1
Curves of elastic damage constitutive relations under uniaxial states.The M-T criterion and M-C criterion could be written as [20, 21](14)F1=−ft+σ1,(15)F2=−σ3−sinφ+1sinφ−1σ1−fc,where φ represents the internal frictional angle of coal-rock mass, fc and ft represent compressive strength and tensile strength of coal-rock mass, respectively, while F2 and F1 are threshold functions of damage for shear and tensile mode, respectively.As the damage of coal-rock mass evolves, the mechanical index such as elastic modulus would decrease monotonically according to theory of elastic damage [15, 21].(16)E=E0⋅1−D,where E0 is the original elastic modulus for coal-rock mass without damage and D is called the damage variable, whose value is 1.0 for coal-rock mass with full damage and 0.0 for the coal-rock mass without damage.If coal-rock mass is in the uniaxial stress state, the curves of constitutive relationships are shown in Figure1. When F1 > 0, namely, the coal-rock mass is under tensile state, the value of D can be calculated by(17)D=−εt0εn+1,ε≥εt0,0,ε<εt0,where εt0 denotes limit strain before tensile damage occurs.Hooke’s law is given as [17].(18)σi=εi−v2v−1εjjEv+1.According to M-T criterion and equation (18), the following relation can be obtained with η=−ν/v+12v−1.(19)E=ftεv⋅η+ε1/v+1.Based on equations (16) and (19), the expression of D under triaxial stress condition is rewritten as(20)D=−ftηE0εv+E0ε1/v+1+1.Corresponding to its damage evolution under tensile condition, whenF2 > 0, namely, the coal-rock mass is damaged in shear mode as illustrated in Figure 1, the expression of D is calculated by(21)D=1−εc0εn,ε≤εc0,0,ε>εc0.When the coal-rock mass is under triaxial stress condition, the expression ofD can be deduced from equations (15) and (18)(22)D=1−fcE0εvη+ε3/v+1−κεvη+ε1/v+1.When coal-rock mass experiences dilatancy due to its strength or stiffness degradation as damage evolves, its hydraulic property would also change. The porosity,ϕ, strongly associated with rock’s stress state, could be given by [22](23)ϕ=ϕr+ϕ0−ϕrexpσ¯e⋅αϕ,where ϕr denotes the residual porosity under extremely large stress condition, ϕ0 denotes the initial porosity under zero-stress condition, αϕ denotes the coefficient of stress sensitivity, and σ¯e is the average of effective stress calculated as σ¯e=αp+σ3+σ2+σ1/3 with σ3, σ2, and σ1 being the third, the second, and the first total principal stress in three directions.Moreover, as coal-rock mass damage accumulates, its permeability would also change accordingly. According to the related studies of relation between porosity and permeability [22], the damage evolution of permeability k can be described using an exponential function as follows:(24)k=k0ϕϕ0βexpD⋅αD,where k0 represents the initial permeability in zero-stress state, αD is the damage-permeability effect coefficient of 5.0, and β is the index coefficient of 3. The value of D in equation (24) can be calculated by equations (20) and (22). Through some numerical simulation examples given below, equation (24) can well describe the permeability damage evolution of coal-rock mass.


## 3.2. Governing Equation for Gas Flow during Gas Drainage

To describe the gas drainage process including coal-rock mass deformation (the same as that already introduced in Section3.1.1) and gas flow after fracturing, except for the previous hypothesis above, the established governing equations used to describe the gas flow process are based on the following hypotheses: (1) the coal-rock mass is completely saturated by single gas; (2) gas existed in pores is the ideal gas, whose viscosity is kept constant; and (3) the rule of gas flow in coal-rock mass satisfies Darcy’s law and is deemed as an isothermal process.The governing equation of gas mass conservation under isothermal state is given by [23–26](25)Qy=∇⋅ρgqg+∂C∂t,where Qy represents the gas sink or source and ρg is the gas density, which is associated with gas pressure as [26, 27].(26)ρg=MgpTR=pβ.q
g is the vector of Darcy velocity for free-phase gas, which is expressed using Darcy’s law with its gravity effect neglected:(27)qg=−kμg∇p.C is the gas content retained in coal-rock mass, which is given by Langmuir’s formula [27, 28]:(28)C=βa2a1ρsa2p+1+ϕp0p2.In equations (26)–(28), Mg represents gas molecular mass, T represents the Kelvin temperature, p represents gas pressure, R represents ideal gas constant, β represents compressibility coefficient of gas, μg represents gas kinematic viscosity coefficient, k is the coal-rock mass permeability; a2 and a1 are called Langmuir’s coefficients with their units being Pa−1 and m3·kg−1, respectively, ρs is the coal-rock mass density, and p0 represents the standard atmospheric pressure.The ultimate governing equation of gas flow is derived by substituting equations (26)–(28) in equation (25):(29)2pβϕp0−ρspa22a12pa2+12+a2a1ρspa2+1∂p∂t−∇⋅2pβkμg∇p=Qy.In brief, two damage-flow coupling processes are included in the simultaneously coupled mathematical model. equations (7) and (13) are the governing equations with damage-flow coupling during hydraulic fracturing, and equations (29) and (7) are governing equations with gas flow and coal-rock mass deformation coupling during gas drainage. The influence of fluids (including gas and water) seepage on the deformation behavior is considered in equation (7). Besides, the porosity and permeability, strongly associated with its stress state, are considered in equations (23) and (24). Furthermore, the influence of deformation field on the flow behaviors of water and gas are implicit in equations (13) and (29), respectively.


## 4. Numerical Implementation of Mathematical Model

The aforementioned equations, especially equations (7), (13), and (29), are nonlinear partial differential equations, and their analytical solution is difficult to solve in theory. Therefore, a numerical method is used instead. The flow chart of the numerical implementation procedure of mathematical model is presented in Figure 2, and its basic procedure can be divided into four steps:Figure 2
Flow chart of the numerical implementation procedure of mathematical model.Step 1.
After the model geometry is determined, it is then discretized into numerous mesoscopic elements (such as triangular or quadrilateral elements). To ensure that the distribution of elements’ properties conforms to specific distributions, a Monte Carlo simulation, coded by Matlab program, is adopted to assign hydraulic or mechanical property of coal-rock mass [14]. Meanwhile, to guarantee the static response of numerical model, prescribed boundary stress is divided into discrete stress increment, which is gradually loaded on the boundary of numerical model.Step 2.
For each stress increment, a fully coupled solution is executed using the Comsol solver [29] to solve governing equations (7) and (13), and the solution results of each element’s strain and effective stress are stored for the next step.Step 3.
The damage state of each element is judged by substituting the stored effective stress into equations (14) and (15). If it is in damage state, the value of D on each element node is calculated using equations (17)–(22). Furthermore, the elastic modulus E and permeability k with damage on each element node are calculated by equations (16) and (24), thus forming a new finite element model with updated mechanical and hydraulic properties. Then, steps (2) and (3) are repeated to check the stress condition of each element. A new stress increment is loaded unless the given convergence criterion is reached, namely, ||(Di+1 − Di)/Di+1||∞ < 1 × 10−5, where D is the matrix of damage variable and i denotes the computation step of numerical model.Step 4.
According to the characteristic of damage and hydraulic and mechanical properties after hydraulic fracturing, another fully coupled solution is executed using Comsol solver according to governing equations (7) and (29), and then, the change of pore pressure with time can be obtained. As illustrated in Figure 2, only the dotted portion represents the computation procedure for the gas flow during gas drainage after hydraulic fracturing.
In brief, the proposed procedure is coded into a Matlab program to calculate the related damage variables in constitutive relation and partly uses Comsol solver to obtain solutions of governing equations with damage-flow coupling.


## 5. Numerical Model Validation and a Case Study

### 5.1. Validation of Mathematical Model and Procedures

In order to verify the numerical model, two numerical examples are presented. The first example is a single solid model without considering flow-damage coupling. The second example analyzes the damage evolution of hydraulic fracturing considering flow-damage coupling.Example 1.
Figure3 presents the numerical model geometry and its boundary conditions for the first example. The related coal-rock mass parameters used are shown in Table 1. The upper boundary stress σv and horizontal boundary stress σh applied on model are 15.45 MPa and 10.3 MPa, respectively. The initial borehole wall pressure σn is 0 MPa, and the pressurization rate is set to 0.1 MPa/step to analyze the initial damage evolution.
The distribution diagram of initial elastic modulus with homogeneity indexes of 1000 is shown in Figure4 (only the mechanical heterogeneity is examined here). As seen from Figure 4, the values of elastic modulus in entire computation domain range from 5.94 GPa to 6.02 GPa, and their deviations from the average value of 6 GPa are relatively small; thus, it can be assumed as a homogeneous material.
According to elastic mechanics, the circumferential stress surrounding borehole can be defined [30]:(30)σθ=3R4r4+1σv−σh2cos2θ+R2r2+1σv+σh2−R2r2σn,where R is the distance between point location and borehole center, r is borehole radius, and θ is the intersection angle from horizontal stress, as shown in Figure 4.
As the value ofσv is greater than σh, it can be concluded that σθA < σθB. Thus, the Point A is the position where tensile damage could occur first, and its critical pressure is calculated as σn=3σh−σv+σt. The numerical simulation results indicates that a initial damage occurs at Point A when the borehole wall pressure σn reaches 21.40 MPa, which agrees well with the theoretical solution of 21.45 MPa. Additionally, the tensile yielded elements using Phase2 software are also consistent with the damaged elements using the method in this paper, as shown in Figure 5. The aforementioned results validate the effectiveness of the mathematical model and its numerical procedure presented in this paper.Figure 3
Geometry and boundary conditions of numerical model.Table 1
Parameters of coal-rock mass for numerical calculation.
Parameters
Values
Parameters
Values
Poisson’s ratiov
0.25
Homogeneity indexm
1000
Elastic modulusE
6 GPa
Internal friction angleφ
30°
Tensile strengthft
6 MPa
Compressive strengthfc
60 MPaFigure 4
Distribution diagram of initial elastic modulus (m = 1000).Figure 5
Tensile yielded elements or damaged elements surrounding borehole. (a) Yielded elements using Phase2 software. (b) Damaged elements using method in this paper.
(a)
(b)Example 2
The geometry and definite conditions of numerical model for the second example are shown in Figure6. Two models, whose homogeneity indexes are 20 and 1000, are established in this example. Related coal-rock mass parameters used are listed in Table 2. The upper boundary stress σv and horizontal boundary stress σh applied are 15.45 MPa and 10.3 MPa, respectively. The original water pressure p is 0 MPa, and the pressurization rate of borehole wall is set to 0.5 MPa/step to analyze damage evolution of hydraulic fracturing considering flow-damage coupling.
For the homogeneous model with homogeneity index ofm = 1000, the simulation results show that if hydraulic pressure p reaches 13.5 MPa, two nearly straight cracks (represented by damaged elements) occur, whose propagation directions are approximately perpendicular to horizontal stress direction. There are two classical equations currently used to estimate the breakdown pressure under far field stress condition. One is the Hubbert–Willis (H-W) equation, which is applied in impermeable coal-rock mass where the injecting fluid cannot penetrate into the borehole wall, while the other is the Haimson–Fairhurst (H-F) equation, which is applied in permeable coal-rock mass where the injecting fluid can permeate into the borehole wall [14]. The two different equations are as follows:(31)pHF=−ft−3σh+σH+2p02η−2+p0,(32)pHW=ft−σH+3σh−p0,where pHF is the breakdown pressure calculated by the H-F equation, pHW is the breakdown pressure calculated by the H-W equation, σh denotes minor horizontal stress, σH denotes maximum horizontal stress, p0 denotes original water pressure existed in coal-rock mass, and η=α2v−1/2v−2.
According to equations (31) and (32), the theoretical solutions for the permeable and impermeable coal-rock mass are 11.6 MPa and 18.45 MPa, respectively. It is worth noting that the breakdown pressure, calculated by equation (31) or (32), is on the basis of a tensile strength criterion. Namely, when the circumferential stress on borehole equals its coal-rock mass tensile strength, a breakdown that leads to full failure of coal-rock mass would occur. However, engineering practices indicate that the propagation of hydraulic cracks would keep stable unless the fluid pressure applied on borehole wall is increased to another critical pressure [11, 12]. Thus, the pressure estimated using equation (31) or (32) can be regarded as the fracture initiation pressure instead of the breakdown pressure of coal-rock mass. Moreover, as H-F equation and H-W equation represent the theoretical solution of fracture initiation pressure for permeable coal-rock mass and impermeable coal-rock mass, they should be regarded as the lower and upper limit values for fracture initiation pressure, respectively. Simulated fracture initiation pressure of 13.5 MPa is between the lower limit (11.6 MPa) and upper limit (18.45 MPa), and the simulated fracture initiation pressure approaches the lower H-F solution due to the reason that both the models have considered the fluid seepage influence on stress, which also demonstrates the excellent performance of the numerical model and its numerical procedure.
Figure7 presents the distribution diagram of initial elastic modulus with homogeneity index of 20.
As illustrated in Figure7, the values of elastic modulus in the entire computation domain range from 2.00 GPa to 3.88 GPa, and their deviations from the average value of 3.35 GPa are relatively large; thus, it can be assumed as a heterogenous material. The result shows that the Monte Carlo simulation can be used to numerically produce a heterogeneous coal-rock mass model.
To analyze the effect of homogeneity index on stress distribution within coal-rock mass, Figure8 presents the distribution diagram of the first principal stress with homogeneity indexes of 1000 and 20. As illustrated in Figure 8, the stress distribution of the heterogenous model would exhibit discontinuity and roughness with respect to the homogeneous model under the same condition. The reason for this is that grains and defects are randomly distributed in coal-rock mass with lower homogeneity index; when the coal-rock mass is loaded, the distributions of inner stress would be nonuniform due to the behavior differences of deformations and force transfer of rock constitutes; when local stress concentrations occur, it would cause the propagation of local microcracks, thus changing the failure mode of coal-rock mass.
For a heterogeneous model with homogeneity index ofm = 20, simulation results show that when hydraulic pressure p in borehole reaches 12 MPa, two nearly straight cracks occur (represented by damaged elements), whose propagation directions are approximately perpendicular to the horizontal stress direction, as shown in Figure 9.
As shown in Figure9, a loading process is represented by notation of “Step 12_25,” in which “12” denotes the hydraulic pressure of 12 MPa applied at borehole wall and “25” represents the iteration number in a load increment under current hydraulic pressure. It can be understood that when the initial crack occurs under specific hydraulic pressure, the water flows into the crack, the equilibrium state is disturbed, the mechanical and hydraulic properties need to be updated, and then the coupling calculation should be calculated again to obtain a new convergence state. The numerical simulation results show that fracture initiation pressure of 12 MPa with homogeneity index of 20 is lower than that of 13.5 MPa with homogeneity index of 1000. As mentioned above, with the increase of m, the material would become more homogeneous, thus leading to a higher strength of material. In short, to accurately depict the crack initiation and evolution during coal-rock mass fracturing, it is necessary to consider realistic heterogeneity of coal-rock mass for engineering practice of hydraulic fracturing.
The whole evolution of hydraulic cracks (represented by damaged elements) in coal-rock mass with homogeneity index of 20 is shown in Figure10. As illustrated in Figure 10, the whole evolution process can be divided into three steps, namely, the initiation process (Figure 10(a)), propagation process (Figures 10(b) and 10(c)), and coalescence process (Figure 10(d)).
As can be seen from Figures9 and 10, the initial cracks occur when the pressure on borehole wall reaches the fracture initiation pressure, but hydraulic cracks do not keep propagating only if the pressure on borehole wall is increased. When the pressure on borehole wall increases to an extent, hydraulic cracks propagate unsteadily until the failure of coal-rock mass reaches. Thus, it can be inferred that for heterogenous coal-rock mass, even if hydraulic cracks resulting from increasing pressure on borehole wall are initiated, they do not immediately result in full failure of coal-rock mass, which also demonstrates that the breakdown pressure is not consistent with fracture initiation pressure. As the evolution of hydraulic fractures, the hydraulic pressure is increased. In other words, the breakdown pressure mass is greater than its fracture initiation pressure, which is meaningful to comprehend the instability mechanism of practical coal-rock mass engineering induced by fracturing [12, 31]. Based on the aforementioned results, it can be seen that the mathematical model and its procedure proposed are able to accurately describe the mechanical processes during hydraulic fracturing.Figure 6
Definite and geometry conditions for numerical simulation.Table 2
Parameters of coal-rock mass and water for numerical calculation.
Parameters
Values
Parameters
Values
Poisson’s ratiov
0.33
Initial permeabilityk0
3 × 10−18 m2
Elastic modulusE
3.35 GPa
Initial porosityϕ0
0.05
Internal friction angleφ
30°
Residual porosityϕr
0.01
Compressive strengthfc
40 MPa
Water dynamic viscosityμw
3.55 × 10−4 Pa·s
Tensile strengthft
3 MPa
Bulk modulus of waterβl
2.5 × 109 Pa
Homogeneity indexm
1000, 20
Biot coefficientα
0.8Figure 7
Distribution diagram of initial elastic modulus (m = 20).Figure 8
Distribution diagram of the first principal stress. (a)m = 1000. (b) m = 20.
(a)
(b)Figure 9
Distribution of initial cracks around borehole (p = 12 MPa). (a) Step 12_01. (b) Step 12_25. (c) Step 12_50. (d) Step 12_99.
(a)
(b)
(c)
(d)Figure 10
Distribution diagram of cracks evolution during hydraulic fracturing process (m = 20). (a) Step 12_99. (b) Step 14_55. (c) Step 15.5_28. (d) Step 18.5_9.
(a)
(b)
(c)
(d)


### 5.2. Numerical Simulation Study on a Case

A numerical model of hydraulic fracturing with two boreholes for enhancing gas drainage is built according to an engineering practice of Nantong mine district in China [32]. The parameters used for numerical calculation are shown in Table 3. The definite conditions and geometry for the numerical model are presented in Figure 11.Table 3
Parameters for the numerical calculation.
Material parameters
Values
Material parameters
Values
Elastic modulusE
3.35 GPa
Initial permeabilityk0
3 × 10−18 m2
Poisson’s ratiov
0.33
Initial porosityϕ0
0.05
Uniaxial compressive strengthfc
40 MPa
Residual porosityϕr
0.01
Internal friction angleφ
30°
Water dynamic viscosityμw
3.55 × 10−4 Pa·s
Uniaxial tensile strengthft
3 MPa
Bulk modulus of waterβl
2.5 GPa
Homogeneity indexm
1.5
Biot coefficientα
0.8
Langmuir coefficienta1
4 × 10−2 m−3·kg−1
Langmuir coefficienta2
8 × 103 Pa−1
Coal densityρs
1.25 × 103 kg·m−3
Compressibility coefficient of gasβ
1.18 × 10−5 kgPa−1·m−3
Gas dynamical viscosityμg
1.84 × 10−5 Pa·s
Initial gas pressurepw
2.1 MPaFigure 11
Geometry and definite conditions of numerical model.The borehole 1 or borehole 2 is called the fracturing borehole as it is used for fracturing, while other boreholes are called control boreholes (as illustrated in Figure11), whose aim is to make the hydraulic fracture easily propagate to the directions of them [33–35]. It is understood that strengths of surrounding coal-rock mass of control boreholes would decline after they have been excavated; thus, they can be regarded as weak parts in coal-rock mass during fracturing. Two simultaneously coupled processes are involved in this case, one is the seepage-damage coupling process during hydraulic fracturing, and the other is gas flow process in damaged coal-rock mass during gas drainage. Therefore, the numerical simulation results and corresponding analysis are divided into two parts. Part 1 is the simulation of hydraulic fracturing with two boreholes, and part 2 is the simulation of gas flow during gas drainage based on the studies of part 1.Part 1.
The initial hydraulic pressure applied at walls of fracturing boreholes are both 7 MPa, while the hydraulic pressures of control boreholes are 0 MPa, and both the vertical and horizontal distances between two boreholes are 3 m, as shown in Figure11. The pressurization rates of fracturing boreholes are set to 0.5 MPa/step (for transient analysis, 1 step represents 0.5 minutes in this case) to analyze the damage evolution. The simulation results show that when hydraulic pressure p reaches 8 MPa, two straight initial cracks occur (represented by damaged elements), whose propagation directions are perpendicular to the horizontal stress direction, as shown in Figure 12(a). Moreover, the whole initiation, propagation, and coalescence of hydraulic cracks (represented by damaged elements) and their corresponding pore pressure evolution during hydraulic fracturing process are presented in Figure 12.
In the present case, the lateral coefficientσh/σv is set to 0.5. Most of the cracks propagate parallel to the orientation of vertical stress during primary process of hydraulic fracturing (Figure 12(a)); no larger deflection occurs, which demonstrates that extension direction of hydraulic cracks is dominated by initial ground stress of coal-rock mass but not the control holes [33]. As the computation step increases (namely, increase of hydraulic pressure and time), the effective stresses increase, and the cracks further propagate under the water wedge with high pressure (Figure 12(b)); the damage extent nearby the zones of crack tip is increased largely, which could greatly increase the permeability of coal-rock mass, thus providing a high flow channel for gas seepage after fracturing. When the distance of hydraulic cracks surrounding borehole 1 and 2 reduces to an extent (Figure 12(c)), the pore pressures on the crack propagation paths of two fracturing boreholes would be larger than that of other directions, then the cracks exhibit deflection tendency to each other due to the effect of local pore pressure, which would also increase the coalesce velocity of two cracks. When the hydraulic pressure reaches 24 MPa (namely, 17 minutes) in Figure 12(d), the two fracturing boreholes are connected through hydraulic cracks. It can be seen from the aforementioned simulation results that the initial ground stress of coal-rock mass and local pore pressure gradient combine to influence the crack propagation direction. Because the dominated directions of initial ground stress and local pore pressure gradient are almost identical in this case, the control boreholes seem to have no effect in this study. It could be inferred that if a rational layout of multiple boreholes is designed, then hydraulic pressures resulting from injected fluid in boreholes can be applied to generate an appropriate gradient of pore pressure to control the trajectory of fracture propagation.Figure 12
Evolution process of pore pressure (figure on left side) and cracks (figure on right side). (a) Step 9.5_45 (t = 2.5 min). (b) Step 14.5_3 (t = 7.5 min). (c) Step 19.5_23 (t = 12.5 min). (d) Step 24_9 (t = 17 min).
(a)
(b)
(c)
(d)Part 2.
After fracturing, based on the simulation results of “Step 14.5_3” in the first part, a fully coupling simulation based on governing equations (7) and (29) is conducted. The original gas pressure in solution domain is 2.1 MPa, the negative pressures in gas drainage boreholes is 20 KPa, and related parameters of this numerical computation are presented as shown in Table 3. In China, gas pressure of coal-rock mass is used to evaluate gas outburst risk. Therefore, compared with gas drainage in coal-rock mass without hydraulic fracturing, if the gas pressure can be reduced to a larger extent through gas drainage after hydraulic fracturing, then measure of hydraulic fracturing is deemed to be effective. The distribution diagram of gas pressure under varying gas drainage times is presented in Figure 13.
As seen in Figure13, due to the initial nonuniform distribution of hydraulic cracks after fracturing, the initial permeability of coal-rock mass would also exhibit nonuniform distribution, thus resulting in the nonuniform distribution of gas pressure.
The curves for gas pressure under varying gas drainage times on the I-I survey line and II-II survey line (as shown in Figure11) are presented in Figure 14.
Moreover, the curves for gas pressure distribution on the vertical fracture path during gas drainage in borehole 1 are presented in Figure15.
It can be seen from Figures(13)–(15) that as gas drainage time increases, gas pressure surrounding fracturing boreholes would decline gradually, and the influence regions of gas drainage increase gradually; the decrease degree of gas pressure on the path of fracture propagation is greater than that of other directions; the distributions of gas pressure surrounding fracturing boreholes are extremely nonuniform, which indicate that the gas pressure of damaged zones declines significantly compared with the undamaged zones; and the gas pressure with a distance of 0.5 m from center of fracturing borehole 1 is 1.5 MPa, while the gas pressure with the same distance from center of control borehole 5 is 2.05 MPa, indicating that the reduction rate of gas pressure has increased by 27% due to hydraulic fracturing. In summary, after hydraulic fracturing, pore-fissure structures within coal-rock mass can be improved, there would exist a gas migration channel with characteristics of lower stress and higher permeability, which demonstrate that hydraulic fracturing has significant influences on improving the coal-rock permeability and thus enhancing effectiveness of gas drainage.Figure 13
Distribution diagram of gas pressure under varying gas drainage times. (a) Two days. (b) Six days. (c) Ten days.
(a)
(b)
(c)Figure 14
Curves for gas pressure under varying gas drainage times on survey lines. (a) Horizontal survey line. (b) Vertical survey line.
(a)
(b)Figure 15
Curves for gas pressure distribution on the vertical fracture path.


## 5.1. Validation of Mathematical Model and Procedures

In order to verify the numerical model, two numerical examples are presented. The first example is a single solid model without considering flow-damage coupling. The second example analyzes the damage evolution of hydraulic fracturing considering flow-damage coupling.Example 1.
Figure3 presents the numerical model geometry and its boundary conditions for the first example. The related coal-rock mass parameters used are shown in Table 1. The upper boundary stress σv and horizontal boundary stress σh applied on model are 15.45 MPa and 10.3 MPa, respectively. The initial borehole wall pressure σn is 0 MPa, and the pressurization rate is set to 0.1 MPa/step to analyze the initial damage evolution.
The distribution diagram of initial elastic modulus with homogeneity indexes of 1000 is shown in Figure4 (only the mechanical heterogeneity is examined here). As seen from Figure 4, the values of elastic modulus in entire computation domain range from 5.94 GPa to 6.02 GPa, and their deviations from the average value of 6 GPa are relatively small; thus, it can be assumed as a homogeneous material.
According to elastic mechanics, the circumferential stress surrounding borehole can be defined [30]:(30)σθ=3R4r4+1σv−σh2cos2θ+R2r2+1σv+σh2−R2r2σn,where R is the distance between point location and borehole center, r is borehole radius, and θ is the intersection angle from horizontal stress, as shown in Figure 4.
As the value ofσv is greater than σh, it can be concluded that σθA < σθB. Thus, the Point A is the position where tensile damage could occur first, and its critical pressure is calculated as σn=3σh−σv+σt. The numerical simulation results indicates that a initial damage occurs at Point A when the borehole wall pressure σn reaches 21.40 MPa, which agrees well with the theoretical solution of 21.45 MPa. Additionally, the tensile yielded elements using Phase2 software are also consistent with the damaged elements using the method in this paper, as shown in Figure 5. The aforementioned results validate the effectiveness of the mathematical model and its numerical procedure presented in this paper.Figure 3
Geometry and boundary conditions of numerical model.Table 1
Parameters of coal-rock mass for numerical calculation.
Parameters
Values
Parameters
Values
Poisson’s ratiov
0.25
Homogeneity indexm
1000
Elastic modulusE
6 GPa
Internal friction angleφ
30°
Tensile strengthft
6 MPa
Compressive strengthfc
60 MPaFigure 4
Distribution diagram of initial elastic modulus (m = 1000).Figure 5
Tensile yielded elements or damaged elements surrounding borehole. (a) Yielded elements using Phase2 software. (b) Damaged elements using method in this paper.
(a)
(b)Example 2
The geometry and definite conditions of numerical model for the second example are shown in Figure6. Two models, whose homogeneity indexes are 20 and 1000, are established in this example. Related coal-rock mass parameters used are listed in Table 2. The upper boundary stress σv and horizontal boundary stress σh applied are 15.45 MPa and 10.3 MPa, respectively. The original water pressure p is 0 MPa, and the pressurization rate of borehole wall is set to 0.5 MPa/step to analyze damage evolution of hydraulic fracturing considering flow-damage coupling.
For the homogeneous model with homogeneity index ofm = 1000, the simulation results show that if hydraulic pressure p reaches 13.5 MPa, two nearly straight cracks (represented by damaged elements) occur, whose propagation directions are approximately perpendicular to horizontal stress direction. There are two classical equations currently used to estimate the breakdown pressure under far field stress condition. One is the Hubbert–Willis (H-W) equation, which is applied in impermeable coal-rock mass where the injecting fluid cannot penetrate into the borehole wall, while the other is the Haimson–Fairhurst (H-F) equation, which is applied in permeable coal-rock mass where the injecting fluid can permeate into the borehole wall [14]. The two different equations are as follows:(31)pHF=−ft−3σh+σH+2p02η−2+p0,(32)pHW=ft−σH+3σh−p0,where pHF is the breakdown pressure calculated by the H-F equation, pHW is the breakdown pressure calculated by the H-W equation, σh denotes minor horizontal stress, σH denotes maximum horizontal stress, p0 denotes original water pressure existed in coal-rock mass, and η=α2v−1/2v−2.
According to equations (31) and (32), the theoretical solutions for the permeable and impermeable coal-rock mass are 11.6 MPa and 18.45 MPa, respectively. It is worth noting that the breakdown pressure, calculated by equation (31) or (32), is on the basis of a tensile strength criterion. Namely, when the circumferential stress on borehole equals its coal-rock mass tensile strength, a breakdown that leads to full failure of coal-rock mass would occur. However, engineering practices indicate that the propagation of hydraulic cracks would keep stable unless the fluid pressure applied on borehole wall is increased to another critical pressure [11, 12]. Thus, the pressure estimated using equation (31) or (32) can be regarded as the fracture initiation pressure instead of the breakdown pressure of coal-rock mass. Moreover, as H-F equation and H-W equation represent the theoretical solution of fracture initiation pressure for permeable coal-rock mass and impermeable coal-rock mass, they should be regarded as the lower and upper limit values for fracture initiation pressure, respectively. Simulated fracture initiation pressure of 13.5 MPa is between the lower limit (11.6 MPa) and upper limit (18.45 MPa), and the simulated fracture initiation pressure approaches the lower H-F solution due to the reason that both the models have considered the fluid seepage influence on stress, which also demonstrates the excellent performance of the numerical model and its numerical procedure.
Figure7 presents the distribution diagram of initial elastic modulus with homogeneity index of 20.
As illustrated in Figure7, the values of elastic modulus in the entire computation domain range from 2.00 GPa to 3.88 GPa, and their deviations from the average value of 3.35 GPa are relatively large; thus, it can be assumed as a heterogenous material. The result shows that the Monte Carlo simulation can be used to numerically produce a heterogeneous coal-rock mass model.
To analyze the effect of homogeneity index on stress distribution within coal-rock mass, Figure8 presents the distribution diagram of the first principal stress with homogeneity indexes of 1000 and 20. As illustrated in Figure 8, the stress distribution of the heterogenous model would exhibit discontinuity and roughness with respect to the homogeneous model under the same condition. The reason for this is that grains and defects are randomly distributed in coal-rock mass with lower homogeneity index; when the coal-rock mass is loaded, the distributions of inner stress would be nonuniform due to the behavior differences of deformations and force transfer of rock constitutes; when local stress concentrations occur, it would cause the propagation of local microcracks, thus changing the failure mode of coal-rock mass.
For a heterogeneous model with homogeneity index ofm = 20, simulation results show that when hydraulic pressure p in borehole reaches 12 MPa, two nearly straight cracks occur (represented by damaged elements), whose propagation directions are approximately perpendicular to the horizontal stress direction, as shown in Figure 9.
As shown in Figure9, a loading process is represented by notation of “Step 12_25,” in which “12” denotes the hydraulic pressure of 12 MPa applied at borehole wall and “25” represents the iteration number in a load increment under current hydraulic pressure. It can be understood that when the initial crack occurs under specific hydraulic pressure, the water flows into the crack, the equilibrium state is disturbed, the mechanical and hydraulic properties need to be updated, and then the coupling calculation should be calculated again to obtain a new convergence state. The numerical simulation results show that fracture initiation pressure of 12 MPa with homogeneity index of 20 is lower than that of 13.5 MPa with homogeneity index of 1000. As mentioned above, with the increase of m, the material would become more homogeneous, thus leading to a higher strength of material. In short, to accurately depict the crack initiation and evolution during coal-rock mass fracturing, it is necessary to consider realistic heterogeneity of coal-rock mass for engineering practice of hydraulic fracturing.
The whole evolution of hydraulic cracks (represented by damaged elements) in coal-rock mass with homogeneity index of 20 is shown in Figure10. As illustrated in Figure 10, the whole evolution process can be divided into three steps, namely, the initiation process (Figure 10(a)), propagation process (Figures 10(b) and 10(c)), and coalescence process (Figure 10(d)).
As can be seen from Figures9 and 10, the initial cracks occur when the pressure on borehole wall reaches the fracture initiation pressure, but hydraulic cracks do not keep propagating only if the pressure on borehole wall is increased. When the pressure on borehole wall increases to an extent, hydraulic cracks propagate unsteadily until the failure of coal-rock mass reaches. Thus, it can be inferred that for heterogenous coal-rock mass, even if hydraulic cracks resulting from increasing pressure on borehole wall are initiated, they do not immediately result in full failure of coal-rock mass, which also demonstrates that the breakdown pressure is not consistent with fracture initiation pressure. As the evolution of hydraulic fractures, the hydraulic pressure is increased. In other words, the breakdown pressure mass is greater than its fracture initiation pressure, which is meaningful to comprehend the instability mechanism of practical coal-rock mass engineering induced by fracturing [12, 31]. Based on the aforementioned results, it can be seen that the mathematical model and its procedure proposed are able to accurately describe the mechanical processes during hydraulic fracturing.Figure 6
Definite and geometry conditions for numerical simulation.Table 2
Parameters of coal-rock mass and water for numerical calculation.
Parameters
Values
Parameters
Values
Poisson’s ratiov
0.33
Initial permeabilityk0
3 × 10−18 m2
Elastic modulusE
3.35 GPa
Initial porosityϕ0
0.05
Internal friction angleφ
30°
Residual porosityϕr
0.01
Compressive strengthfc
40 MPa
Water dynamic viscosityμw
3.55 × 10−4 Pa·s
Tensile strengthft
3 MPa
Bulk modulus of waterβl
2.5 × 109 Pa
Homogeneity indexm
1000, 20
Biot coefficientα
0.8Figure 7
Distribution diagram of initial elastic modulus (m = 20).Figure 8
Distribution diagram of the first principal stress. (a)m = 1000. (b) m = 20.
(a)
(b)Figure 9
Distribution of initial cracks around borehole (p = 12 MPa). (a) Step 12_01. (b) Step 12_25. (c) Step 12_50. (d) Step 12_99.
(a)
(b)
(c)
(d)Figure 10
Distribution diagram of cracks evolution during hydraulic fracturing process (m = 20). (a) Step 12_99. (b) Step 14_55. (c) Step 15.5_28. (d) Step 18.5_9.
(a)
(b)
(c)
(d)


## 5.2. Numerical Simulation Study on a Case

A numerical model of hydraulic fracturing with two boreholes for enhancing gas drainage is built according to an engineering practice of Nantong mine district in China [32]. The parameters used for numerical calculation are shown in Table 3. The definite conditions and geometry for the numerical model are presented in Figure 11.Table 3
Parameters for the numerical calculation.
Material parameters
Values
Material parameters
Values
Elastic modulusE
3.35 GPa
Initial permeabilityk0
3 × 10−18 m2
Poisson’s ratiov
0.33
Initial porosityϕ0
0.05
Uniaxial compressive strengthfc
40 MPa
Residual porosityϕr
0.01
Internal friction angleφ
30°
Water dynamic viscosityμw
3.55 × 10−4 Pa·s
Uniaxial tensile strengthft
3 MPa
Bulk modulus of waterβl
2.5 GPa
Homogeneity indexm
1.5
Biot coefficientα
0.8
Langmuir coefficienta1
4 × 10−2 m−3·kg−1
Langmuir coefficienta2
8 × 103 Pa−1
Coal densityρs
1.25 × 103 kg·m−3
Compressibility coefficient of gasβ
1.18 × 10−5 kgPa−1·m−3
Gas dynamical viscosityμg
1.84 × 10−5 Pa·s
Initial gas pressurepw
2.1 MPaFigure 11
Geometry and definite conditions of numerical model.The borehole 1 or borehole 2 is called the fracturing borehole as it is used for fracturing, while other boreholes are called control boreholes (as illustrated in Figure11), whose aim is to make the hydraulic fracture easily propagate to the directions of them [33–35]. It is understood that strengths of surrounding coal-rock mass of control boreholes would decline after they have been excavated; thus, they can be regarded as weak parts in coal-rock mass during fracturing. Two simultaneously coupled processes are involved in this case, one is the seepage-damage coupling process during hydraulic fracturing, and the other is gas flow process in damaged coal-rock mass during gas drainage. Therefore, the numerical simulation results and corresponding analysis are divided into two parts. Part 1 is the simulation of hydraulic fracturing with two boreholes, and part 2 is the simulation of gas flow during gas drainage based on the studies of part 1.Part 1.
The initial hydraulic pressure applied at walls of fracturing boreholes are both 7 MPa, while the hydraulic pressures of control boreholes are 0 MPa, and both the vertical and horizontal distances between two boreholes are 3 m, as shown in Figure11. The pressurization rates of fracturing boreholes are set to 0.5 MPa/step (for transient analysis, 1 step represents 0.5 minutes in this case) to analyze the damage evolution. The simulation results show that when hydraulic pressure p reaches 8 MPa, two straight initial cracks occur (represented by damaged elements), whose propagation directions are perpendicular to the horizontal stress direction, as shown in Figure 12(a). Moreover, the whole initiation, propagation, and coalescence of hydraulic cracks (represented by damaged elements) and their corresponding pore pressure evolution during hydraulic fracturing process are presented in Figure 12.
In the present case, the lateral coefficientσh/σv is set to 0.5. Most of the cracks propagate parallel to the orientation of vertical stress during primary process of hydraulic fracturing (Figure 12(a)); no larger deflection occurs, which demonstrates that extension direction of hydraulic cracks is dominated by initial ground stress of coal-rock mass but not the control holes [33]. As the computation step increases (namely, increase of hydraulic pressure and time), the effective stresses increase, and the cracks further propagate under the water wedge with high pressure (Figure 12(b)); the damage extent nearby the zones of crack tip is increased largely, which could greatly increase the permeability of coal-rock mass, thus providing a high flow channel for gas seepage after fracturing. When the distance of hydraulic cracks surrounding borehole 1 and 2 reduces to an extent (Figure 12(c)), the pore pressures on the crack propagation paths of two fracturing boreholes would be larger than that of other directions, then the cracks exhibit deflection tendency to each other due to the effect of local pore pressure, which would also increase the coalesce velocity of two cracks. When the hydraulic pressure reaches 24 MPa (namely, 17 minutes) in Figure 12(d), the two fracturing boreholes are connected through hydraulic cracks. It can be seen from the aforementioned simulation results that the initial ground stress of coal-rock mass and local pore pressure gradient combine to influence the crack propagation direction. Because the dominated directions of initial ground stress and local pore pressure gradient are almost identical in this case, the control boreholes seem to have no effect in this study. It could be inferred that if a rational layout of multiple boreholes is designed, then hydraulic pressures resulting from injected fluid in boreholes can be applied to generate an appropriate gradient of pore pressure to control the trajectory of fracture propagation.Figure 12
Evolution process of pore pressure (figure on left side) and cracks (figure on right side). (a) Step 9.5_45 (t = 2.5 min). (b) Step 14.5_3 (t = 7.5 min). (c) Step 19.5_23 (t = 12.5 min). (d) Step 24_9 (t = 17 min).
(a)
(b)
(c)
(d)Part 2.
After fracturing, based on the simulation results of “Step 14.5_3” in the first part, a fully coupling simulation based on governing equations (7) and (29) is conducted. The original gas pressure in solution domain is 2.1 MPa, the negative pressures in gas drainage boreholes is 20 KPa, and related parameters of this numerical computation are presented as shown in Table 3. In China, gas pressure of coal-rock mass is used to evaluate gas outburst risk. Therefore, compared with gas drainage in coal-rock mass without hydraulic fracturing, if the gas pressure can be reduced to a larger extent through gas drainage after hydraulic fracturing, then measure of hydraulic fracturing is deemed to be effective. The distribution diagram of gas pressure under varying gas drainage times is presented in Figure 13.
As seen in Figure13, due to the initial nonuniform distribution of hydraulic cracks after fracturing, the initial permeability of coal-rock mass would also exhibit nonuniform distribution, thus resulting in the nonuniform distribution of gas pressure.
The curves for gas pressure under varying gas drainage times on the I-I survey line and II-II survey line (as shown in Figure11) are presented in Figure 14.
Moreover, the curves for gas pressure distribution on the vertical fracture path during gas drainage in borehole 1 are presented in Figure15.
It can be seen from Figures(13)–(15) that as gas drainage time increases, gas pressure surrounding fracturing boreholes would decline gradually, and the influence regions of gas drainage increase gradually; the decrease degree of gas pressure on the path of fracture propagation is greater than that of other directions; the distributions of gas pressure surrounding fracturing boreholes are extremely nonuniform, which indicate that the gas pressure of damaged zones declines significantly compared with the undamaged zones; and the gas pressure with a distance of 0.5 m from center of fracturing borehole 1 is 1.5 MPa, while the gas pressure with the same distance from center of control borehole 5 is 2.05 MPa, indicating that the reduction rate of gas pressure has increased by 27% due to hydraulic fracturing. In summary, after hydraulic fracturing, pore-fissure structures within coal-rock mass can be improved, there would exist a gas migration channel with characteristics of lower stress and higher permeability, which demonstrate that hydraulic fracturing has significant influences on improving the coal-rock permeability and thus enhancing effectiveness of gas drainage.Figure 13
Distribution diagram of gas pressure under varying gas drainage times. (a) Two days. (b) Six days. (c) Ten days.
(a)
(b)
(c)Figure 14
Curves for gas pressure under varying gas drainage times on survey lines. (a) Horizontal survey line. (b) Vertical survey line.
(a)
(b)Figure 15
Curves for gas pressure distribution on the vertical fracture path.


## 6. Conclusions

It is important to investigate mechanism of how hydraulic fracturing influences gas drainage in coal-rock mass. In present work, the method for heterogeneity assignment of material properties is proposed; a novel simultaneously coupled mathematical model with flow-damage coupling is established and its numerical implementation procedure is also proposed; then two examples are presented to validate their effectiveness; and finally, numerical simulation of an engineering case on hydraulic fracturing in coal-rock mass to improve gas extraction is presented. The obtained conclusions are as follows:(1)
Monte Carlo simulations can be used to numerically produce a heterogeneous coal-rock mass model including heterogeneity of hydraulic and mechanical properties. To accurately depict the crack initiation and evolution during hydraulic fracturing, it is necessary to consider realistic heterogeneity of coal-rock mass for engineering practice of hydraulic fracturing.(2)
Simulation results of two numerical examples indicate that the simultaneously coupled mathematical model and its numerical implementation procedure proposed can be used to describe the fully coupled hydraulic fracturing processes including seepage-damage coupling.(3)
Simulated fracture initiation pressure is between H-F equation and H-W equation, and it approaches lower H-F equation as two models both have considered seepage influence on stress. The higher the homogeneity index, the greater the fracture initiation pressure. For heterogenous coal-rock mass, breakdown pressure is greater than fracture initiation pressure, which is meaningful to understand the hydraulic fracturing mechanism in practical coal-rock mass engineering.(4)
The engineering case of hydraulic fracturing with two boreholes for improving gas drainage shows that decreased degree of gas pressure on fracture propagation path is greater than that of other directions; the gas pressure of damaged zones declines significantly compared with the undamaged zones; the reduction rate of gas pressure surrounding borehole has increased after fracturing; the pore-fissure structure of coal-rock mass is improved, and there exists a gas migration channel with characteristics of lower stress and higher porosity, which is the main mechanism of hydraulic fracturing on increasing coal-rock permeability and enhancing gas drainage.



---

*Source: 1485672-2018-12-12.xml*