# An Improved Coupled Modeling Method for Coalbed Methane Extraction after Hydraulic Fracturing

Authors: Hao-Hao Zhang; Jun-Xing Wang; Chao-Jun Fan; Hui-Jie Bi; Li Chen; Chen-Yu Xu
Journal: Geofluids
 (2022)

---

## Abstract

This paper proposes an improved coupled modeling method for coalbed methane extraction after hydraulic fracturing. Based on the damage mechanics theory, we introduce a state damage variable to record the damage changes at different times and construct the coupling governing equations of the seepage field and solid mechanics field in the hydraulic fracturing process. Then, the results of hydraulic fracturing simulation are used as initial values to perform coupling modeling for coalbed methane extraction. The method is implemented by the combination of different modules in COMSOL to solve the state damage variable, displacement, and pressure. Additionally, two 2D numerical examples are used to verify the precision of this method, and a 3D coalbed methane extraction simulation is performed to analyze the extraction effects and the evolution characteristics of the permeability before and after hydraulic fracturing. Numerical examples show that the proposed coupling modeling method can not only simplify the calculation steps of the CwM method and reduce the numerical errors caused, but also characterize the damage, porosity, permeability, and extraction effect caused by hydraulic fracturing in coal seams with permeability anisotropy.


---

## Body

## 1. Introduction

In recent years, with the vigorous exploitation of coalbed methane resources and the improvement of extraction technology, to a certain extent, the global energy shortage situation has been alleviated [1, 2], and the incidence of coal mine gas accidents has also been greatly reduced [3–6]. Hydraulic fracturing is often used to enhance the permeability of coal seams and promote the effective exploitation of coalbed methane. However, it is controversial because of the leakage of fracturing fluid or gas and the pollution of groundwater caused by hydraulic fracturing [7]. Therefore, under the background of huge economic benefits and corresponding risks brought by hydraulic fracturing, many numerical methods for predicting the effect of hydraulic fracturing in porous media have emerged.Numerical methods of hydraulic fracturing in porous media such as coal seam or shale rock can be generally divided into classified into continuous and discontinuous methods [8]. It is common for continuous methods to simulate hydraulic fracturing by using special elements and technologies [9–15]. Discontinuous methods emphasize the aspects of discontinuous behavior and large displacement. They are developed based on many discontinuities in rock masses [16, 17], so they are more popular for simulating the damage and failure behavior. But it is difficult to simulate the evolution of porosity and permeability and the effect of coalbed methane extraction after hydraulic fracturing. The primary pores and fractures in the coal seam are penetrated and a large number of secondary fractures are generated after hydraulic fracturing. The porosity and permeability of the coal seam increase significantly, and the coalbed methane is easier to be pumped out. Some scholars [18–21] have analyzed the impact of hydraulic fracturing on coalbed methane extraction based on the continuous method. Among them, a method based on COMSOL with MATLAB (CwM) is widely used. This method can not only simulate coal seam hydraulic fracturing but also analyze the coupling effect of different factors on the effect of coalbed methane extraction. Hydraulic fracturing simulation based on the CwM method is mainly used to extract stress, strain, and other data from COMSOL files in unit time step in the form of matrix data points through MATLAB, and to determine damage in MATLAB. The damage variables, elastic modulus, and other parameters in this step are taken as the initial values of the COMSOL files in the next time step. The precision of this method depends mainly on the size of the unit time step and the density of matrix data points. A reasonable unit time step and matrix data point density can make the simulation results more realistic. However, since the damage judgment of this method is based on matrix data points instead of discretized grids, the implementation process of this method becomes very cumbersome when the geometric model is irregular. In addition, when a long-term hydraulic fracturing simulation is performed, a large number of COMSOL files need to be stored, which easily leads to problems such as insufficient memory and difficulty in data postprocessing.The coal seam is regarded as a typical dual-pore system composed of microporous matrix and fracture, which contains a large amount of adsorbed methane and free methane. It is not only the source of coalbed methane but also the reservoir. Coalbed methane extraction is a typical process of coupling multiple physical fields, which is affected by a variety of factors such as in situ stress, methane adsorption/desorption, groundwater, porosity, and permeability. Xia et al. [22] analyzed the coupling mechanism of coal deformation, methane transportation, and airflow in the coal seam and proposed a coupling composition mode. Hu et al. [23] and Wang et al. [24, 25] consider the influence of coal deformation caused by methane adsorption or desorption to establish a gas-solid coupling model. Li et al. [26, 27] and Yuan et al. [28] considered the influence of coal deformation caused by methane adsorption and desorption, coal seam temperature, and groundwater and established a thermo-fluid-solid coupling model. Fan et al. [21] analyzed the influence of non-Darcy law seepage on coalbed methane extraction in the case of fracturing damage based on the thermal-fluid-solid coupling model. However, because the distribution of pores and fractures in the coal seams presents obvious randomness, the coalbed methane migration process in the extraction process shows anisotropy characteristics. Permeability, as one of the key factors of coalbed methane extraction, mainly controls the seepage velocity of coalbed methane in fractures, and fractures are the main migration channel of coalbed methane and also the main factor affecting permeability [29–31]. Therefore, the anisotropy characteristics of fractures in coal seams can be attributed to the anisotropy of permeability in the process of coalbed methane extraction. Gu et al. [32], Wu et al. [33], and Liu et al. [34] analyzed the influence of permeability anisotropy on methane migration in coal seams and established a permeability anisotropy coupling model.In this paper, based on the CwM method and migration characteristics of coalbed methane in coal seams, we used the state variable to determine the damage in COMSOL and combined it with its multifield coupling characteristics to carry out the coupled modeling for coalbed methane extraction after hydraulic fracturing. And the effect of coalbed methane extraction before and after hydraulic fracturing was simulated and analyzed. The paper is organized as follows. We begin with a short introduction to hydraulic fracturing modeling and coupling modeling for coalbed methane extraction in Section2. In Section 3, we present the numerical implementation method and some 2D numerical examples for coalbed methane extraction after hydraulic fracturing in COMSOL Multiphysics. Finally, we take a 3D numerical example to analyze the effect of coalbed methane extraction after hydraulic fracturing in Section 4.


## 2. Coupled Modeling Governing Equations

### 2.1. Hydraulic Fracturing Modeling

#### 2.1.1. Damage Evolution Model

Coal rock contains mineral particles, cement, and a large number of pores and fractures. It is one of the complex solid materials formed in nature after hundreds of millions of years of geological evolution and tectonic movement. It has obvious heterogeneity characteristics. The heterogeneity of coal rock under external load has a significant influence on its destructive behavior. It can be assumed to be composed of a large number of representative element volumes (REVs), and its heterogeneity can be characterized by the Weibull distribution of the material parameters of different REVs in the material space [35]. A large number of fractures were developed in the coal seam due to the combined action of high-pressure water disturbance and geological tectonic stress, resulting in damage. The elastic modulus of REVs in the coal seam can be described as
(1)E=E01−Hd,where E0 and E are the initial elastic modulus of REVs before and after damage, MPa, and Hd is the state damage variable, which is a historical state variable related to time and strain tensor. The relationship between Hd and the damage-driven evolution function Dd satisfies the Kuhn Tucker condition, which can be described as,
(2)Hdt,εijt=Hd0,εij0t=0maxDdt,εijt,Hdt−1,εijt−1t≥1,where t is the time, s, and εijt is the strain tensor at time t.It can be seen from Equation (2) that when t=0, Hd=Hd (0, εij0), it is used to characterize the initial damage, which is used to describe the damage caused by drilling or other measures before hydraulic fracturing. If the initial damage is ignored, Hd (0, εij0) = 0. When t≥1, Hd is the larger value between the damage-driven evolution function Dd at time t and Hd at time t−1. And the damage-driven evolution function Dd is defined as
(3)Ddt,εijt=0Ft<0,Fc<01−εt0ε12Ft=0,dFt>01−εc0ε32Fc=0,dFc>0,where εt0 and εc0 are the maximum tensile and compressive principal strain, εt0=σt/E0, εc0=σc/E0; σt and σc are the ultimate tensile strength and compressive strength, Pa; ε1 and ε3 are the first and third principal strains; and Ft and Fc are the damage threshold of tensile and compression in REVs.During hydraulic fracturing, when the stress state of the REV satisfies the maximum tensile stress criterion, the tensile damage of the coal seam is produced. Similarly, when the stress state meets the Drucker-Prager criterion, shear damage is produced. They can be described as follows:(4)Ft=σ1−σt=0Fc=2sinθ33−sinθI1+J2−6ccosθ33−sinθ=0,in which I1=σ1+σ2+σ3, J2=1/6σ1−σ22+σ2−σ32+σ3−σ12, where σ1, σ2, and σ3 are the first, second, and third principal stress, Pa; θ is the internal friction angle, °; and c is the cohesion, Pa.


#### 2.1.2. Governing Equations

During hydraulic fracturing, high-pressure water flows into the coal seam through the fracturing borehole and causes damage. Pores and fractures of the coal seam are further developed; the permeability of the coal seam is also significantly improved. It is the result of the coupling of the solid mechanics field and the seepage field. Therefore, the solid mechanics field and the seepage field in hydraulic fracturing can be expressed by the following equations:(5)ρsd2udt2−∇⋅σ+αwpw+ρsg=0,(6)ρwφfχf∂pw∂t+∇⋅ρwvw=Qm,where ρs is the material density, kg/m3; u is the variable of the displacement field, m; t is the time, s; σ is the Cauchy stress tensor, Pa; αw is the water Biot coefficient; pw is the water pressure, Pa; g is the acceleration of gravity, m/s2; φf is the fracture porosity; χf is the compressibility factor, χf=1×10−10; vw is the water flow velocity, m/s; and Qm is the source term.Equation (5) is the solid mechanics field; the third term in the equation is used to characterize the effect of water pressure on the solid mechanics field, which is applied in the form of volume load. Equation (6) is the governing equation of the poroelastic storage water model; the flow velocity vw follows Darcy’s law, and it is defined by the following:
(7)vw=−kfμw∇pw+ρwg,where μw is the water dynamic viscosity, Pa.s, and kf is the fracture permeability, m2.Hydraulic fracturing makes the damage accumulate, and the porosity and permeability of the coal seam increase with the development of the damaged area, which makes it easier for water to enter the coal seam, thereby promoting the effect of hydraulic fracturing. The evolution of porosity and permeability during hydraulic fracturing can be expressed as follows:(8)φf=φf01+ξ1Hdkf=kf01+ξ2Hd,where φf0 and kf0 are the initial fracture porosity and permeability; ξ1 is the increased coefficient of porosity after damage; and ξ2 is the increased coefficient of permeability after damage.Substituting Equation (7) and Equation (8) into Equation (6), Equation (6) can be written as
(9)1+ξ1Hdρwφf0χf∂pw∂t+∇⋅−ρwkf01+ξ2Hdμw∇pw+ρwg=Qm.Equation (5) and Equation (9) are the governing equations of the mechanical field and the seepage field in hydraulic fracturing.


### 2.2. Coupling Modeling for Coalbed Methane Extraction

#### 2.2.1. Permeability Anisotropy Evolution Model

According to literature [26, 27], coalbed methane exists and migrates in the matrix in adsorption and free state, while only free methane exists in fracture. It is assumed that methane in the matrix is mainly transported to fractures by diffusion, and a small amount of methane is transported to fractures by seepage (Figure 1). In this process, methane diffusion follows Fick’s law and seepage follows Darcy’s law. Methane in fractures migrates to gas wells in the way of seepage, and the seepage rate is mainly determined by the fracture opening and gas pressure difference. However, the permeability of coal seam fracture has obvious anisotropy due to the different fracture opening bi (i=x,y,z) in different directions.Figure 1
Methane migration modes from matrix to fracture.As the cuboid fracture model in Figure1, the initial fracture porosity and permeability can be obtained by the following (Equation (19))
(10)ϕf0=∑ibiakfi0=∑i≠jbj312awhere a is the matrix size, m; kfi0 is the initial fracture permeability in i direction, m2; i, j = x, y, z, represent the three directions of the Cartesian coordinate; and bi and bj are the opening degree of the fracture perpendicular to the directions iand j, m.Takingx direction as an example, the relationship between stress and strain per unit volume is
(11)εx=13εLΔpmpL+Δpm−byaLfΔσfx−1KΔσmx,where εx is the volumetric strain per unit volume in x direction; εL and pL are the Langmuir strain constant and pressure constant; Δpm is the matrix methane pressure increment, Pa; Lf is the improved fracture stiffness, Pa; K is the bulk modulus, Pa; and Δσfx and Δσmx are the stress increments of coal matrix and fracture, respectively.The effective stress increment of the matrix and fracture is equal to the total incrementΔσex; that is, Δσmx=Δσfx=Δσex [33], and then, Equation (11) can be modified to
(12)εx=13εLΔpmpL+Δpm−byaLf+1KΔσex.In coal seam, the fracture opening degree incrementΔby is closely related to the effective stress increment and the initial fracture opening degree, and its value can be obtained by
(13)Δby=−byLfΔσex.Then, the fracture porosityφf in the coal seam is
(14)ϕfϕf0=1+∑iΔbi∑ibi.According to the permeability cubic law, the fracture permeability inx direction is
(15)kfxkfx0=121+Δbyby3+121+Δbzbz3.Similarly, the permeability inx and y directions is
(16)kfykfy0=121+Δbxbx3+121+Δbzbz3kfzkfz0=121+Δbxbx3+121+Δbyby3,where kfi represents the fracture permeability in i direction, m2, and bj and Δbj are the fracture opening degree and the opening degree increments perpendicular to j direction, respectively, m.The porosity of the coal matrix is defined as(17)ϕm=11+Sϕm01+S0+αmS‐S0where
(18)S=εv+pm/Ks‐εLpm/pm+pLS0=εv0+pm0/Ks‐εLpm0/pm0+pLαm=1−K/Lfwhere φm0 is the initial matrix porosity; αm is the matrix methane Biot coefficient; Ks is the matrix bulk modulus, Pa; v is the Poisson’s ratio; εv is the volumetric strain; and pm is the matrix methane pressure, Pa.The permeability of the coal matrix can be described as(19)kmkm0=11+S1+S0+αmS‐S0ϕm03where km0 is the initial matrix permeability, m2.Equation (14) is the fracture porosity equation, and Equation (15) and Equation (16) are fracture permeability equations considering anisotropy. Equation (17) and Equation (19) are porosity and permeability equations in matrix permeability.


#### 2.2.2. Governing Equations

Based on the assumption that coal has dual porosity and dual permeability and the above evolution equations of porosity and permeability in matrix and fracture, according to the mass conservation equation, the migration equations of methane in matrix and fracture can be described as(20)∂∂tVLpmpL+pmρsMgRTsps+∂∂tφmMgRTpm+∇⋅−MgkmμgRTpm∇pm=−MgτRTpm−pfg∂∂tφfMgRTpfg+∇⋅−Mgpfg+hRTkfμg∇pfg=1−φfMgτRTpm−pfg,where VL is the Langmuir volume constant, m3/kg; μg is the dynamic viscosity of methane, Pa.s; Ts and ps are the temperature and the methane pressure under standard condition; T is the temperature of the coal seam; τ is the methane desorption time, s; Mg is the methane molar mass, kg/mol; R is the molar constant; pfg is the pressure methane in fracture, Pa; and h is the Klinkenberg factor, Pa.Similarly, based on the Navier equation, the stress field governing equation for coalbed methane migration can be modified to(21)Gu+G1−2vu−αmpm−αfpfg−KεLΔpmpL+Δpm+F=0,where G is the shear modulus, Pa; αf is the fracture methane Biot coefficient; and F is the volume force, Pa.Equation (20) and Equation (21) are the governing equations for coalbed methane extraction considering the permeability anisotropy of coal seam.


## 2.1. Hydraulic Fracturing Modeling

### 2.1.1. Damage Evolution Model

Coal rock contains mineral particles, cement, and a large number of pores and fractures. It is one of the complex solid materials formed in nature after hundreds of millions of years of geological evolution and tectonic movement. It has obvious heterogeneity characteristics. The heterogeneity of coal rock under external load has a significant influence on its destructive behavior. It can be assumed to be composed of a large number of representative element volumes (REVs), and its heterogeneity can be characterized by the Weibull distribution of the material parameters of different REVs in the material space [35]. A large number of fractures were developed in the coal seam due to the combined action of high-pressure water disturbance and geological tectonic stress, resulting in damage. The elastic modulus of REVs in the coal seam can be described as
(1)E=E01−Hd,where E0 and E are the initial elastic modulus of REVs before and after damage, MPa, and Hd is the state damage variable, which is a historical state variable related to time and strain tensor. The relationship between Hd and the damage-driven evolution function Dd satisfies the Kuhn Tucker condition, which can be described as,
(2)Hdt,εijt=Hd0,εij0t=0maxDdt,εijt,Hdt−1,εijt−1t≥1,where t is the time, s, and εijt is the strain tensor at time t.It can be seen from Equation (2) that when t=0, Hd=Hd (0, εij0), it is used to characterize the initial damage, which is used to describe the damage caused by drilling or other measures before hydraulic fracturing. If the initial damage is ignored, Hd (0, εij0) = 0. When t≥1, Hd is the larger value between the damage-driven evolution function Dd at time t and Hd at time t−1. And the damage-driven evolution function Dd is defined as
(3)Ddt,εijt=0Ft<0,Fc<01−εt0ε12Ft=0,dFt>01−εc0ε32Fc=0,dFc>0,where εt0 and εc0 are the maximum tensile and compressive principal strain, εt0=σt/E0, εc0=σc/E0; σt and σc are the ultimate tensile strength and compressive strength, Pa; ε1 and ε3 are the first and third principal strains; and Ft and Fc are the damage threshold of tensile and compression in REVs.During hydraulic fracturing, when the stress state of the REV satisfies the maximum tensile stress criterion, the tensile damage of the coal seam is produced. Similarly, when the stress state meets the Drucker-Prager criterion, shear damage is produced. They can be described as follows:(4)Ft=σ1−σt=0Fc=2sinθ33−sinθI1+J2−6ccosθ33−sinθ=0,in which I1=σ1+σ2+σ3, J2=1/6σ1−σ22+σ2−σ32+σ3−σ12, where σ1, σ2, and σ3 are the first, second, and third principal stress, Pa; θ is the internal friction angle, °; and c is the cohesion, Pa.


### 2.1.2. Governing Equations

During hydraulic fracturing, high-pressure water flows into the coal seam through the fracturing borehole and causes damage. Pores and fractures of the coal seam are further developed; the permeability of the coal seam is also significantly improved. It is the result of the coupling of the solid mechanics field and the seepage field. Therefore, the solid mechanics field and the seepage field in hydraulic fracturing can be expressed by the following equations:(5)ρsd2udt2−∇⋅σ+αwpw+ρsg=0,(6)ρwφfχf∂pw∂t+∇⋅ρwvw=Qm,where ρs is the material density, kg/m3; u is the variable of the displacement field, m; t is the time, s; σ is the Cauchy stress tensor, Pa; αw is the water Biot coefficient; pw is the water pressure, Pa; g is the acceleration of gravity, m/s2; φf is the fracture porosity; χf is the compressibility factor, χf=1×10−10; vw is the water flow velocity, m/s; and Qm is the source term.Equation (5) is the solid mechanics field; the third term in the equation is used to characterize the effect of water pressure on the solid mechanics field, which is applied in the form of volume load. Equation (6) is the governing equation of the poroelastic storage water model; the flow velocity vw follows Darcy’s law, and it is defined by the following:
(7)vw=−kfμw∇pw+ρwg,where μw is the water dynamic viscosity, Pa.s, and kf is the fracture permeability, m2.Hydraulic fracturing makes the damage accumulate, and the porosity and permeability of the coal seam increase with the development of the damaged area, which makes it easier for water to enter the coal seam, thereby promoting the effect of hydraulic fracturing. The evolution of porosity and permeability during hydraulic fracturing can be expressed as follows:(8)φf=φf01+ξ1Hdkf=kf01+ξ2Hd,where φf0 and kf0 are the initial fracture porosity and permeability; ξ1 is the increased coefficient of porosity after damage; and ξ2 is the increased coefficient of permeability after damage.Substituting Equation (7) and Equation (8) into Equation (6), Equation (6) can be written as
(9)1+ξ1Hdρwφf0χf∂pw∂t+∇⋅−ρwkf01+ξ2Hdμw∇pw+ρwg=Qm.Equation (5) and Equation (9) are the governing equations of the mechanical field and the seepage field in hydraulic fracturing.


## 2.1.1. Damage Evolution Model

Coal rock contains mineral particles, cement, and a large number of pores and fractures. It is one of the complex solid materials formed in nature after hundreds of millions of years of geological evolution and tectonic movement. It has obvious heterogeneity characteristics. The heterogeneity of coal rock under external load has a significant influence on its destructive behavior. It can be assumed to be composed of a large number of representative element volumes (REVs), and its heterogeneity can be characterized by the Weibull distribution of the material parameters of different REVs in the material space [35]. A large number of fractures were developed in the coal seam due to the combined action of high-pressure water disturbance and geological tectonic stress, resulting in damage. The elastic modulus of REVs in the coal seam can be described as
(1)E=E01−Hd,where E0 and E are the initial elastic modulus of REVs before and after damage, MPa, and Hd is the state damage variable, which is a historical state variable related to time and strain tensor. The relationship between Hd and the damage-driven evolution function Dd satisfies the Kuhn Tucker condition, which can be described as,
(2)Hdt,εijt=Hd0,εij0t=0maxDdt,εijt,Hdt−1,εijt−1t≥1,where t is the time, s, and εijt is the strain tensor at time t.It can be seen from Equation (2) that when t=0, Hd=Hd (0, εij0), it is used to characterize the initial damage, which is used to describe the damage caused by drilling or other measures before hydraulic fracturing. If the initial damage is ignored, Hd (0, εij0) = 0. When t≥1, Hd is the larger value between the damage-driven evolution function Dd at time t and Hd at time t−1. And the damage-driven evolution function Dd is defined as
(3)Ddt,εijt=0Ft<0,Fc<01−εt0ε12Ft=0,dFt>01−εc0ε32Fc=0,dFc>0,where εt0 and εc0 are the maximum tensile and compressive principal strain, εt0=σt/E0, εc0=σc/E0; σt and σc are the ultimate tensile strength and compressive strength, Pa; ε1 and ε3 are the first and third principal strains; and Ft and Fc are the damage threshold of tensile and compression in REVs.During hydraulic fracturing, when the stress state of the REV satisfies the maximum tensile stress criterion, the tensile damage of the coal seam is produced. Similarly, when the stress state meets the Drucker-Prager criterion, shear damage is produced. They can be described as follows:(4)Ft=σ1−σt=0Fc=2sinθ33−sinθI1+J2−6ccosθ33−sinθ=0,in which I1=σ1+σ2+σ3, J2=1/6σ1−σ22+σ2−σ32+σ3−σ12, where σ1, σ2, and σ3 are the first, second, and third principal stress, Pa; θ is the internal friction angle, °; and c is the cohesion, Pa.


## 2.1.2. Governing Equations

During hydraulic fracturing, high-pressure water flows into the coal seam through the fracturing borehole and causes damage. Pores and fractures of the coal seam are further developed; the permeability of the coal seam is also significantly improved. It is the result of the coupling of the solid mechanics field and the seepage field. Therefore, the solid mechanics field and the seepage field in hydraulic fracturing can be expressed by the following equations:(5)ρsd2udt2−∇⋅σ+αwpw+ρsg=0,(6)ρwφfχf∂pw∂t+∇⋅ρwvw=Qm,where ρs is the material density, kg/m3; u is the variable of the displacement field, m; t is the time, s; σ is the Cauchy stress tensor, Pa; αw is the water Biot coefficient; pw is the water pressure, Pa; g is the acceleration of gravity, m/s2; φf is the fracture porosity; χf is the compressibility factor, χf=1×10−10; vw is the water flow velocity, m/s; and Qm is the source term.Equation (5) is the solid mechanics field; the third term in the equation is used to characterize the effect of water pressure on the solid mechanics field, which is applied in the form of volume load. Equation (6) is the governing equation of the poroelastic storage water model; the flow velocity vw follows Darcy’s law, and it is defined by the following:
(7)vw=−kfμw∇pw+ρwg,where μw is the water dynamic viscosity, Pa.s, and kf is the fracture permeability, m2.Hydraulic fracturing makes the damage accumulate, and the porosity and permeability of the coal seam increase with the development of the damaged area, which makes it easier for water to enter the coal seam, thereby promoting the effect of hydraulic fracturing. The evolution of porosity and permeability during hydraulic fracturing can be expressed as follows:(8)φf=φf01+ξ1Hdkf=kf01+ξ2Hd,where φf0 and kf0 are the initial fracture porosity and permeability; ξ1 is the increased coefficient of porosity after damage; and ξ2 is the increased coefficient of permeability after damage.Substituting Equation (7) and Equation (8) into Equation (6), Equation (6) can be written as
(9)1+ξ1Hdρwφf0χf∂pw∂t+∇⋅−ρwkf01+ξ2Hdμw∇pw+ρwg=Qm.Equation (5) and Equation (9) are the governing equations of the mechanical field and the seepage field in hydraulic fracturing.


## 2.2. Coupling Modeling for Coalbed Methane Extraction

### 2.2.1. Permeability Anisotropy Evolution Model

According to literature [26, 27], coalbed methane exists and migrates in the matrix in adsorption and free state, while only free methane exists in fracture. It is assumed that methane in the matrix is mainly transported to fractures by diffusion, and a small amount of methane is transported to fractures by seepage (Figure 1). In this process, methane diffusion follows Fick’s law and seepage follows Darcy’s law. Methane in fractures migrates to gas wells in the way of seepage, and the seepage rate is mainly determined by the fracture opening and gas pressure difference. However, the permeability of coal seam fracture has obvious anisotropy due to the different fracture opening bi (i=x,y,z) in different directions.Figure 1
Methane migration modes from matrix to fracture.As the cuboid fracture model in Figure1, the initial fracture porosity and permeability can be obtained by the following (Equation (19))
(10)ϕf0=∑ibiakfi0=∑i≠jbj312awhere a is the matrix size, m; kfi0 is the initial fracture permeability in i direction, m2; i, j = x, y, z, represent the three directions of the Cartesian coordinate; and bi and bj are the opening degree of the fracture perpendicular to the directions iand j, m.Takingx direction as an example, the relationship between stress and strain per unit volume is
(11)εx=13εLΔpmpL+Δpm−byaLfΔσfx−1KΔσmx,where εx is the volumetric strain per unit volume in x direction; εL and pL are the Langmuir strain constant and pressure constant; Δpm is the matrix methane pressure increment, Pa; Lf is the improved fracture stiffness, Pa; K is the bulk modulus, Pa; and Δσfx and Δσmx are the stress increments of coal matrix and fracture, respectively.The effective stress increment of the matrix and fracture is equal to the total incrementΔσex; that is, Δσmx=Δσfx=Δσex [33], and then, Equation (11) can be modified to
(12)εx=13εLΔpmpL+Δpm−byaLf+1KΔσex.In coal seam, the fracture opening degree incrementΔby is closely related to the effective stress increment and the initial fracture opening degree, and its value can be obtained by
(13)Δby=−byLfΔσex.Then, the fracture porosityφf in the coal seam is
(14)ϕfϕf0=1+∑iΔbi∑ibi.According to the permeability cubic law, the fracture permeability inx direction is
(15)kfxkfx0=121+Δbyby3+121+Δbzbz3.Similarly, the permeability inx and y directions is
(16)kfykfy0=121+Δbxbx3+121+Δbzbz3kfzkfz0=121+Δbxbx3+121+Δbyby3,where kfi represents the fracture permeability in i direction, m2, and bj and Δbj are the fracture opening degree and the opening degree increments perpendicular to j direction, respectively, m.The porosity of the coal matrix is defined as(17)ϕm=11+Sϕm01+S0+αmS‐S0where
(18)S=εv+pm/Ks‐εLpm/pm+pLS0=εv0+pm0/Ks‐εLpm0/pm0+pLαm=1−K/Lfwhere φm0 is the initial matrix porosity; αm is the matrix methane Biot coefficient; Ks is the matrix bulk modulus, Pa; v is the Poisson’s ratio; εv is the volumetric strain; and pm is the matrix methane pressure, Pa.The permeability of the coal matrix can be described as(19)kmkm0=11+S1+S0+αmS‐S0ϕm03where km0 is the initial matrix permeability, m2.Equation (14) is the fracture porosity equation, and Equation (15) and Equation (16) are fracture permeability equations considering anisotropy. Equation (17) and Equation (19) are porosity and permeability equations in matrix permeability.


### 2.2.2. Governing Equations

Based on the assumption that coal has dual porosity and dual permeability and the above evolution equations of porosity and permeability in matrix and fracture, according to the mass conservation equation, the migration equations of methane in matrix and fracture can be described as(20)∂∂tVLpmpL+pmρsMgRTsps+∂∂tφmMgRTpm+∇⋅−MgkmμgRTpm∇pm=−MgτRTpm−pfg∂∂tφfMgRTpfg+∇⋅−Mgpfg+hRTkfμg∇pfg=1−φfMgτRTpm−pfg,where VL is the Langmuir volume constant, m3/kg; μg is the dynamic viscosity of methane, Pa.s; Ts and ps are the temperature and the methane pressure under standard condition; T is the temperature of the coal seam; τ is the methane desorption time, s; Mg is the methane molar mass, kg/mol; R is the molar constant; pfg is the pressure methane in fracture, Pa; and h is the Klinkenberg factor, Pa.Similarly, based on the Navier equation, the stress field governing equation for coalbed methane migration can be modified to(21)Gu+G1−2vu−αmpm−αfpfg−KεLΔpmpL+Δpm+F=0,where G is the shear modulus, Pa; αf is the fracture methane Biot coefficient; and F is the volume force, Pa.Equation (20) and Equation (21) are the governing equations for coalbed methane extraction considering the permeability anisotropy of coal seam.


## 2.2.1. Permeability Anisotropy Evolution Model

According to literature [26, 27], coalbed methane exists and migrates in the matrix in adsorption and free state, while only free methane exists in fracture. It is assumed that methane in the matrix is mainly transported to fractures by diffusion, and a small amount of methane is transported to fractures by seepage (Figure 1). In this process, methane diffusion follows Fick’s law and seepage follows Darcy’s law. Methane in fractures migrates to gas wells in the way of seepage, and the seepage rate is mainly determined by the fracture opening and gas pressure difference. However, the permeability of coal seam fracture has obvious anisotropy due to the different fracture opening bi (i=x,y,z) in different directions.Figure 1
Methane migration modes from matrix to fracture.As the cuboid fracture model in Figure1, the initial fracture porosity and permeability can be obtained by the following (Equation (19))
(10)ϕf0=∑ibiakfi0=∑i≠jbj312awhere a is the matrix size, m; kfi0 is the initial fracture permeability in i direction, m2; i, j = x, y, z, represent the three directions of the Cartesian coordinate; and bi and bj are the opening degree of the fracture perpendicular to the directions iand j, m.Takingx direction as an example, the relationship between stress and strain per unit volume is
(11)εx=13εLΔpmpL+Δpm−byaLfΔσfx−1KΔσmx,where εx is the volumetric strain per unit volume in x direction; εL and pL are the Langmuir strain constant and pressure constant; Δpm is the matrix methane pressure increment, Pa; Lf is the improved fracture stiffness, Pa; K is the bulk modulus, Pa; and Δσfx and Δσmx are the stress increments of coal matrix and fracture, respectively.The effective stress increment of the matrix and fracture is equal to the total incrementΔσex; that is, Δσmx=Δσfx=Δσex [33], and then, Equation (11) can be modified to
(12)εx=13εLΔpmpL+Δpm−byaLf+1KΔσex.In coal seam, the fracture opening degree incrementΔby is closely related to the effective stress increment and the initial fracture opening degree, and its value can be obtained by
(13)Δby=−byLfΔσex.Then, the fracture porosityφf in the coal seam is
(14)ϕfϕf0=1+∑iΔbi∑ibi.According to the permeability cubic law, the fracture permeability inx direction is
(15)kfxkfx0=121+Δbyby3+121+Δbzbz3.Similarly, the permeability inx and y directions is
(16)kfykfy0=121+Δbxbx3+121+Δbzbz3kfzkfz0=121+Δbxbx3+121+Δbyby3,where kfi represents the fracture permeability in i direction, m2, and bj and Δbj are the fracture opening degree and the opening degree increments perpendicular to j direction, respectively, m.The porosity of the coal matrix is defined as(17)ϕm=11+Sϕm01+S0+αmS‐S0where
(18)S=εv+pm/Ks‐εLpm/pm+pLS0=εv0+pm0/Ks‐εLpm0/pm0+pLαm=1−K/Lfwhere φm0 is the initial matrix porosity; αm is the matrix methane Biot coefficient; Ks is the matrix bulk modulus, Pa; v is the Poisson’s ratio; εv is the volumetric strain; and pm is the matrix methane pressure, Pa.The permeability of the coal matrix can be described as(19)kmkm0=11+S1+S0+αmS‐S0ϕm03where km0 is the initial matrix permeability, m2.Equation (14) is the fracture porosity equation, and Equation (15) and Equation (16) are fracture permeability equations considering anisotropy. Equation (17) and Equation (19) are porosity and permeability equations in matrix permeability.


## 2.2.2. Governing Equations

Based on the assumption that coal has dual porosity and dual permeability and the above evolution equations of porosity and permeability in matrix and fracture, according to the mass conservation equation, the migration equations of methane in matrix and fracture can be described as(20)∂∂tVLpmpL+pmρsMgRTsps+∂∂tφmMgRTpm+∇⋅−MgkmμgRTpm∇pm=−MgτRTpm−pfg∂∂tφfMgRTpfg+∇⋅−Mgpfg+hRTkfμg∇pfg=1−φfMgτRTpm−pfg,where VL is the Langmuir volume constant, m3/kg; μg is the dynamic viscosity of methane, Pa.s; Ts and ps are the temperature and the methane pressure under standard condition; T is the temperature of the coal seam; τ is the methane desorption time, s; Mg is the methane molar mass, kg/mol; R is the molar constant; pfg is the pressure methane in fracture, Pa; and h is the Klinkenberg factor, Pa.Similarly, based on the Navier equation, the stress field governing equation for coalbed methane migration can be modified to(21)Gu+G1−2vu−αmpm−αfpfg−KεLΔpmpL+Δpm+F=0,where G is the shear modulus, Pa; αf is the fracture methane Biot coefficient; and F is the volume force, Pa.Equation (20) and Equation (21) are the governing equations for coalbed methane extraction considering the permeability anisotropy of coal seam.


## 3. Implementation Method and Model Validation

### 3.1. Implementation Method

Coupled modeling governing equations in Section2 can be implemented in COMSOL Multiphysics software to simulate the coalbed methane extraction after hydraulic fracturing. The mechanical field and the water seepage field of the hydraulic fracturing model adopt the modified solid mechanics module and Darcy’s law module. The mechanical field, the matrix methane seepage field, and the fracture methane seepage field for coalbed methane extraction adopt the modified solid mechanics module and coefficient form PDE modules. The relationship between all modules used can be described in Figure 2, which shows that coalbed methane extraction simulation after hydraulic fracturing can be carried out in the following two stages:Figure 2
Relationship between all the modules adopted.


#### 3.1.1. Stage I: Hydraulic Fracturing Stage

In this stage, the water pressurepw of the seepage field is applied to the mechanical field in the form of body load. Under the combined action of water pressure pw and other loads, when the REV satisfies the damage criterion (Equation (4)), damage occurs, and the damage takes the form of the damage-driven evolution function Dd. Meanwhile, the state damage variable Hd records the damage at this moment and transmits it to the elastic modulus E of the REV. As the elastic modulus of some element nodes weakened, the stress field will be adjusted. And the porosity φf and the permeability kfi of the damaged area are transferred to the water seepage field, realizing the coupling effect of the mechanical field and the water seepage field.


#### 3.1.2. Stage II: Coalbed Methane Extraction Stage

The state damage variableHd, elastic modulus E, porosity φf, and permeability kfi of the last time step of stage I are used as the initial variables for coalbed methane extraction simulation. This stage mainly involves the coupling relationship of solid mechanics field, matrix methane seepage field, and fractured methane seepage field. It can be seen from Equation (20) that the matrix methane pressure and the fracture methane pressure are one of the variables in the original terms of each other, and there is a material exchange between the two fields. Both and the adsorbed methane strain act on the evolution of the stress and strain of the solid mechanics field. The stress-strain evolution of solid mechanics, in turn, acts on the matrix methane seepage field and the fracture methane seepage field in the form of porosity and permeability.Figure3 is the flowchart of our implementation method for coalbed methane extraction after hydraulic fracturing in COMSOL Multiphysics. The details are as follows: (1) Set parameters for all fields in the global definition node; (2) create two variable nodes to input damage variables and extraction variables, respectively, and create a state variable node to input state damage variable Hd; (3) create 2D or 3D geometry; (4) set up a solid mechanics module, a Darcy’s law module, a coefficient form PDE (pm), and a coefficient form PDE (pfg), and input parameters, variables, initial and boundary conditions, and take standard Lagrangian elements to discretize the space domain for all fields; (5) select the solid mechanics module and the Darcy’s law module, disable the coefficient form PDE modules, and create and optimize the geometry meshes; and (6) select the solver and set the tolerance to solve and save data. It should be noted here that the iterative convergence speed is slow, and it is necessary to select a separate solver and use the Anderson acceleration method to accelerate convergence [36]. (7) Update variables such as elastic modulus E, state damage variable Hd, porosity φf, and permeability kfi; (8) select the solid mechanics module and the coefficient form PDE modules, disable the Darcy’s law module, and create and optimize the geometry meshes; and (9) repeat step (6).Figure 3
COMSOL implementation of coalbed methane extraction after hydraulic fracturing.


### 3.2. Model Validation

In this section, we have performed a 2D hydraulic fracturing simulation and a 2D coalbed methane extraction simulation after hydraulic fracturing. In the hydraulic fracturing simulation, we tested the influence of different water injection pressures on the coalbed fracturing effect. In the coalbed methane extraction simulation, a permeability anisotropy coefficient is introduced to test the influence of permeability anisotropy on the effect of hydraulic fracturing and coalbed methane extraction. The impact of fracturing and extracting of a single vertical well on the coal seam mainly occurs near the well wall of the gas well. If the simulated coal seam range is too large, it will consume a lot of calculation resources when taking measures such as heterogeneity of elastic modulus and strength and grid subdivision, and the accuracy of numerical simulation is difficult to ensure. Therefore, both simulations are carried out on the geometry as shown in Figure4. The geometric size is 10 m ×10 m, there is a circular hole in the middle at the geometric center, and the hole diameter is 0.1 m.Figure 4
2D geometry used in both simulations.


#### 3.2.1. Hydraulic Fracturing Simulation

The boundary around the geometry is set as the fixed boundary of the solid mechanical field and the impermeable boundary of the water seepage field, and the central hole is set as the water injection hole. And the following material parameters are adopted:E0=8.47GPa, v=0.23, ρs=1650kg/m3, σt=0.10MPa, σc=16.80MPa, ξ1=5, ξ2=5,φf0=0.018, and kf0=0.01mD. The water injection lasts for 120 minutes, and we select the water injection pressure pw0=1MPa, pw0=3MPa, pw0=5MPa, and pw0=10MPa for simulation analysis. Grid generation of the domain, including 8408 elements, is done using the free triangular mesh method.The state damage variableHd is used to record the damaged area of the coal seam, so the surface integral of Hd in the simulation domain could obtain the distribution of the damaged area. Figure 5 is the distribution of damaged area in coal seam under different water injection pressure at 5 min. As the material of the simulation domain follows the Weibull distribution, and the permeability is set as isotropic, the distribution of the damaged area is radially outward from the central hole. It can be seen from Figure 6 that damaged area size increases with the water pressure increases, showing a trend of gradually increasing first and then stabilizing, and the greater the water injection pressure, the longer it takes for the damaged area to develop.Figure 5
Distribution of the damaged area under different water injection pressures at 5 min. (a)pw0=1MPa, (b) pw0=3MPa, (c) pw0=5MPa, and (d) pw0=10MPa.
(a)(b)(c)(d)Figure 6
Damaged area size under different water injection pressures.Based on the geometry in Figure4, the proposed method is compared with the CwM method [18–21] to verify its precision. Figures 7 and 8 are the damaged area distribution cloud diagram and damaged area size variation curves, respectively, of the CwM method at different time steps Δt and the proposed method when the injection pressure pw0=5MPa. It can be seen from the damage area distribution that there are certain differences in the damage distribution of the two methods; at the same time, only the time step Δt=0.1min of the CwM method is similar to the proposed method. And the damaged area size is shown as CwM, Δt=0.1min<CwM, Δt=0.5min<Proposedmethod, Hd<CwM, Δt=1min. The damage area size of the CwM method increases with the time step Δt. When using the CwM method to simulate hydraulic fracturing, the total damage is obtained by the accumulation of the damage generated by the simulation within a unit time step. A longer time step will result in a larger cumulative error in the simulation. When the time step is smaller, the error can be reduced, but due to the time difference between the update of the damage variable and the variables of the stress field and seepage field, it is easy to cause incomplete development of the damaged area and underestimate the size. In addition, because the principle of this method is based on the mutual conversion between the matrix data points of MATLAB and the geometric grid in COMSOL, because the matrix data points and the grid are not in a one-to-one correspondence, there is a little error. However, the proposed method is carried out on the geometric grid in COMSOL, and the state damage variable Hd is updated in real time with the dependent variables such as water pressure and displacement when solving, which could avoid these problems to some extent.Figure 7
Distribution of damaged area in the CwM method and proposed method at 5 min. (a) CwM,Δt=1min; (b) CwM, Δt=0.5min; (c) CwM, Δt=0.1min; and (d) proposed method, Hd.
(a)(b)(c)(d)Figure 8
Damaged area size of the CwM method and proposed method.


#### 3.2.2. Coalbed Methane Extraction Simulation after Hydraulic Fracturing

According to Section3.1, the coalbed methane extraction simulation after hydraulic fracturing can be divided into two stages. In stage I, the boundary conditions and parameters are similar to the hydraulic fracturing simulation, and the water injection pressure is set as pw0=5MPa. In stage II, the boundary around the geometry in Figure 4 is set as the fixed boundary of the solid mechanical field and the impermeable boundary of the methane seepage field, and the central hole is set as the drainage hole. The relevant parameters of coalbed methane extraction are as follows: pfg0=pm0=1.5MPa, τ=153360s, φm0=0.035, PL=0.93MPa, VL=37.04m3/t, km0=0.002mD, and kfx0=0.01mD. And the point M in the geometry is the permeability monitoring point. To analyze the effect of permeability anisotropy, we introduce an anisotropy coefficient λ to describe the evolution characteristics of coal seam fracture permeability, λ=kfy0/kfx0, λ∈(0,1). The anisotropy coefficients λ=0.1, λ=0.2, λ=0.3, and λ=0.5 are adopted for simulation analysis. The coalbed methane extraction simulation lasts for 180 days.It can be seen from Figure9 that damaged area size mainly develops in x direction of coal seam under different anisotropy coefficients. With increasing anisotropy coefficient λ, the damaged area size in xdirection decreases gradually, while the damaged area size in y direction increases gradually. And the total damaged area size increases with the increase of the anisotropy coefficient λ (Figure 10(a)). The distributions of gas pressure in Figure 11 show obvious anisotropic characteristics. The smaller the anisotropy coefficient λ is, the larger the difference of gas pressure distribution in x direction andydirection is. Figure 10 shows the damaged area size and gas production rate under different anisotropy coefficients. The development of the damaged area is affected by the anisotropy of the permeability. When λ=0.1, the development of the damaged area is slower, and the development time of the damage lasts about 100 minutes. And when λ=0.5, the development of the damaged area is faster and lasts about 100 minutes. This is because the permeability in ydirection in the model is small, and x direction is the main seepage direction of the injected water. The water pressure of the damaged area in a short time is not enough to cause the failure of the undamaged coal body. It requires continuous water injection to increase the water pressure of the damaged area further fractures. The coal seam gas production rate also shows a gradual decrease as time increases (Figure 10(b)). The permeability of the damaged area increases, and the gas is easier to be pumped out, while the gas in the undamaged area has a lower permeability and the gas is not easy to be pumped out.Figure 9
Distribution of damaged area under different anisotropy coefficients at 120 min. (a)λ=0.1, (b) λ=0.2, (c) λ=0.3, and (d) λ=0.5.
(a)(b)(c)(d)Figure 10
Damaged area size and gas production rate under different anisotropy coefficients. (a) Damaged area size; (b) Gas production rate.
(a)(b)Figure 11
Distribution of gas pressure under different anisotropy coefficients at 180 days. (a)λ=0.1, (b) λ=0.2, (c) λ=0.3, and (d) λ=0.5.Figure12 shows the evolution law of the matrix and fracture permeability at point M in stage II. The overall trend of the ratio of matrix and fracture permeability is polynomial increase with the time increases, and the increase of matrix permeability is greater than fracture permeability, and the increase of the fracture permeability ratio in x direction is greater than that in ydirection. As the anisotropy coefficient increases, the ratio of matrix and fracture permeability both increase significantly.Figure 12
Matrix permeability ratio and fracture permeability ratio under different anisotropy coefficients: (a) matrix permeability ratio and (b) fracture permeability ratio.
(a)(b)


## 3.1. Implementation Method

Coupled modeling governing equations in Section2 can be implemented in COMSOL Multiphysics software to simulate the coalbed methane extraction after hydraulic fracturing. The mechanical field and the water seepage field of the hydraulic fracturing model adopt the modified solid mechanics module and Darcy’s law module. The mechanical field, the matrix methane seepage field, and the fracture methane seepage field for coalbed methane extraction adopt the modified solid mechanics module and coefficient form PDE modules. The relationship between all modules used can be described in Figure 2, which shows that coalbed methane extraction simulation after hydraulic fracturing can be carried out in the following two stages:Figure 2
Relationship between all the modules adopted.


### 3.1.1. Stage I: Hydraulic Fracturing Stage

In this stage, the water pressurepw of the seepage field is applied to the mechanical field in the form of body load. Under the combined action of water pressure pw and other loads, when the REV satisfies the damage criterion (Equation (4)), damage occurs, and the damage takes the form of the damage-driven evolution function Dd. Meanwhile, the state damage variable Hd records the damage at this moment and transmits it to the elastic modulus E of the REV. As the elastic modulus of some element nodes weakened, the stress field will be adjusted. And the porosity φf and the permeability kfi of the damaged area are transferred to the water seepage field, realizing the coupling effect of the mechanical field and the water seepage field.


### 3.1.2. Stage II: Coalbed Methane Extraction Stage

The state damage variableHd, elastic modulus E, porosity φf, and permeability kfi of the last time step of stage I are used as the initial variables for coalbed methane extraction simulation. This stage mainly involves the coupling relationship of solid mechanics field, matrix methane seepage field, and fractured methane seepage field. It can be seen from Equation (20) that the matrix methane pressure and the fracture methane pressure are one of the variables in the original terms of each other, and there is a material exchange between the two fields. Both and the adsorbed methane strain act on the evolution of the stress and strain of the solid mechanics field. The stress-strain evolution of solid mechanics, in turn, acts on the matrix methane seepage field and the fracture methane seepage field in the form of porosity and permeability.Figure3 is the flowchart of our implementation method for coalbed methane extraction after hydraulic fracturing in COMSOL Multiphysics. The details are as follows: (1) Set parameters for all fields in the global definition node; (2) create two variable nodes to input damage variables and extraction variables, respectively, and create a state variable node to input state damage variable Hd; (3) create 2D or 3D geometry; (4) set up a solid mechanics module, a Darcy’s law module, a coefficient form PDE (pm), and a coefficient form PDE (pfg), and input parameters, variables, initial and boundary conditions, and take standard Lagrangian elements to discretize the space domain for all fields; (5) select the solid mechanics module and the Darcy’s law module, disable the coefficient form PDE modules, and create and optimize the geometry meshes; and (6) select the solver and set the tolerance to solve and save data. It should be noted here that the iterative convergence speed is slow, and it is necessary to select a separate solver and use the Anderson acceleration method to accelerate convergence [36]. (7) Update variables such as elastic modulus E, state damage variable Hd, porosity φf, and permeability kfi; (8) select the solid mechanics module and the coefficient form PDE modules, disable the Darcy’s law module, and create and optimize the geometry meshes; and (9) repeat step (6).Figure 3
COMSOL implementation of coalbed methane extraction after hydraulic fracturing.


## 3.1.1. Stage I: Hydraulic Fracturing Stage

In this stage, the water pressurepw of the seepage field is applied to the mechanical field in the form of body load. Under the combined action of water pressure pw and other loads, when the REV satisfies the damage criterion (Equation (4)), damage occurs, and the damage takes the form of the damage-driven evolution function Dd. Meanwhile, the state damage variable Hd records the damage at this moment and transmits it to the elastic modulus E of the REV. As the elastic modulus of some element nodes weakened, the stress field will be adjusted. And the porosity φf and the permeability kfi of the damaged area are transferred to the water seepage field, realizing the coupling effect of the mechanical field and the water seepage field.


## 3.1.2. Stage II: Coalbed Methane Extraction Stage

The state damage variableHd, elastic modulus E, porosity φf, and permeability kfi of the last time step of stage I are used as the initial variables for coalbed methane extraction simulation. This stage mainly involves the coupling relationship of solid mechanics field, matrix methane seepage field, and fractured methane seepage field. It can be seen from Equation (20) that the matrix methane pressure and the fracture methane pressure are one of the variables in the original terms of each other, and there is a material exchange between the two fields. Both and the adsorbed methane strain act on the evolution of the stress and strain of the solid mechanics field. The stress-strain evolution of solid mechanics, in turn, acts on the matrix methane seepage field and the fracture methane seepage field in the form of porosity and permeability.Figure3 is the flowchart of our implementation method for coalbed methane extraction after hydraulic fracturing in COMSOL Multiphysics. The details are as follows: (1) Set parameters for all fields in the global definition node; (2) create two variable nodes to input damage variables and extraction variables, respectively, and create a state variable node to input state damage variable Hd; (3) create 2D or 3D geometry; (4) set up a solid mechanics module, a Darcy’s law module, a coefficient form PDE (pm), and a coefficient form PDE (pfg), and input parameters, variables, initial and boundary conditions, and take standard Lagrangian elements to discretize the space domain for all fields; (5) select the solid mechanics module and the Darcy’s law module, disable the coefficient form PDE modules, and create and optimize the geometry meshes; and (6) select the solver and set the tolerance to solve and save data. It should be noted here that the iterative convergence speed is slow, and it is necessary to select a separate solver and use the Anderson acceleration method to accelerate convergence [36]. (7) Update variables such as elastic modulus E, state damage variable Hd, porosity φf, and permeability kfi; (8) select the solid mechanics module and the coefficient form PDE modules, disable the Darcy’s law module, and create and optimize the geometry meshes; and (9) repeat step (6).Figure 3
COMSOL implementation of coalbed methane extraction after hydraulic fracturing.


## 3.2. Model Validation

In this section, we have performed a 2D hydraulic fracturing simulation and a 2D coalbed methane extraction simulation after hydraulic fracturing. In the hydraulic fracturing simulation, we tested the influence of different water injection pressures on the coalbed fracturing effect. In the coalbed methane extraction simulation, a permeability anisotropy coefficient is introduced to test the influence of permeability anisotropy on the effect of hydraulic fracturing and coalbed methane extraction. The impact of fracturing and extracting of a single vertical well on the coal seam mainly occurs near the well wall of the gas well. If the simulated coal seam range is too large, it will consume a lot of calculation resources when taking measures such as heterogeneity of elastic modulus and strength and grid subdivision, and the accuracy of numerical simulation is difficult to ensure. Therefore, both simulations are carried out on the geometry as shown in Figure4. The geometric size is 10 m ×10 m, there is a circular hole in the middle at the geometric center, and the hole diameter is 0.1 m.Figure 4
2D geometry used in both simulations.


### 3.2.1. Hydraulic Fracturing Simulation

The boundary around the geometry is set as the fixed boundary of the solid mechanical field and the impermeable boundary of the water seepage field, and the central hole is set as the water injection hole. And the following material parameters are adopted:E0=8.47GPa, v=0.23, ρs=1650kg/m3, σt=0.10MPa, σc=16.80MPa, ξ1=5, ξ2=5,φf0=0.018, and kf0=0.01mD. The water injection lasts for 120 minutes, and we select the water injection pressure pw0=1MPa, pw0=3MPa, pw0=5MPa, and pw0=10MPa for simulation analysis. Grid generation of the domain, including 8408 elements, is done using the free triangular mesh method.The state damage variableHd is used to record the damaged area of the coal seam, so the surface integral of Hd in the simulation domain could obtain the distribution of the damaged area. Figure 5 is the distribution of damaged area in coal seam under different water injection pressure at 5 min. As the material of the simulation domain follows the Weibull distribution, and the permeability is set as isotropic, the distribution of the damaged area is radially outward from the central hole. It can be seen from Figure 6 that damaged area size increases with the water pressure increases, showing a trend of gradually increasing first and then stabilizing, and the greater the water injection pressure, the longer it takes for the damaged area to develop.Figure 5
Distribution of the damaged area under different water injection pressures at 5 min. (a)pw0=1MPa, (b) pw0=3MPa, (c) pw0=5MPa, and (d) pw0=10MPa.
(a)(b)(c)(d)Figure 6
Damaged area size under different water injection pressures.Based on the geometry in Figure4, the proposed method is compared with the CwM method [18–21] to verify its precision. Figures 7 and 8 are the damaged area distribution cloud diagram and damaged area size variation curves, respectively, of the CwM method at different time steps Δt and the proposed method when the injection pressure pw0=5MPa. It can be seen from the damage area distribution that there are certain differences in the damage distribution of the two methods; at the same time, only the time step Δt=0.1min of the CwM method is similar to the proposed method. And the damaged area size is shown as CwM, Δt=0.1min<CwM, Δt=0.5min<Proposedmethod, Hd<CwM, Δt=1min. The damage area size of the CwM method increases with the time step Δt. When using the CwM method to simulate hydraulic fracturing, the total damage is obtained by the accumulation of the damage generated by the simulation within a unit time step. A longer time step will result in a larger cumulative error in the simulation. When the time step is smaller, the error can be reduced, but due to the time difference between the update of the damage variable and the variables of the stress field and seepage field, it is easy to cause incomplete development of the damaged area and underestimate the size. In addition, because the principle of this method is based on the mutual conversion between the matrix data points of MATLAB and the geometric grid in COMSOL, because the matrix data points and the grid are not in a one-to-one correspondence, there is a little error. However, the proposed method is carried out on the geometric grid in COMSOL, and the state damage variable Hd is updated in real time with the dependent variables such as water pressure and displacement when solving, which could avoid these problems to some extent.Figure 7
Distribution of damaged area in the CwM method and proposed method at 5 min. (a) CwM,Δt=1min; (b) CwM, Δt=0.5min; (c) CwM, Δt=0.1min; and (d) proposed method, Hd.
(a)(b)(c)(d)Figure 8
Damaged area size of the CwM method and proposed method.


### 3.2.2. Coalbed Methane Extraction Simulation after Hydraulic Fracturing

According to Section3.1, the coalbed methane extraction simulation after hydraulic fracturing can be divided into two stages. In stage I, the boundary conditions and parameters are similar to the hydraulic fracturing simulation, and the water injection pressure is set as pw0=5MPa. In stage II, the boundary around the geometry in Figure 4 is set as the fixed boundary of the solid mechanical field and the impermeable boundary of the methane seepage field, and the central hole is set as the drainage hole. The relevant parameters of coalbed methane extraction are as follows: pfg0=pm0=1.5MPa, τ=153360s, φm0=0.035, PL=0.93MPa, VL=37.04m3/t, km0=0.002mD, and kfx0=0.01mD. And the point M in the geometry is the permeability monitoring point. To analyze the effect of permeability anisotropy, we introduce an anisotropy coefficient λ to describe the evolution characteristics of coal seam fracture permeability, λ=kfy0/kfx0, λ∈(0,1). The anisotropy coefficients λ=0.1, λ=0.2, λ=0.3, and λ=0.5 are adopted for simulation analysis. The coalbed methane extraction simulation lasts for 180 days.It can be seen from Figure9 that damaged area size mainly develops in x direction of coal seam under different anisotropy coefficients. With increasing anisotropy coefficient λ, the damaged area size in xdirection decreases gradually, while the damaged area size in y direction increases gradually. And the total damaged area size increases with the increase of the anisotropy coefficient λ (Figure 10(a)). The distributions of gas pressure in Figure 11 show obvious anisotropic characteristics. The smaller the anisotropy coefficient λ is, the larger the difference of gas pressure distribution in x direction andydirection is. Figure 10 shows the damaged area size and gas production rate under different anisotropy coefficients. The development of the damaged area is affected by the anisotropy of the permeability. When λ=0.1, the development of the damaged area is slower, and the development time of the damage lasts about 100 minutes. And when λ=0.5, the development of the damaged area is faster and lasts about 100 minutes. This is because the permeability in ydirection in the model is small, and x direction is the main seepage direction of the injected water. The water pressure of the damaged area in a short time is not enough to cause the failure of the undamaged coal body. It requires continuous water injection to increase the water pressure of the damaged area further fractures. The coal seam gas production rate also shows a gradual decrease as time increases (Figure 10(b)). The permeability of the damaged area increases, and the gas is easier to be pumped out, while the gas in the undamaged area has a lower permeability and the gas is not easy to be pumped out.Figure 9
Distribution of damaged area under different anisotropy coefficients at 120 min. (a)λ=0.1, (b) λ=0.2, (c) λ=0.3, and (d) λ=0.5.
(a)(b)(c)(d)Figure 10
Damaged area size and gas production rate under different anisotropy coefficients. (a) Damaged area size; (b) Gas production rate.
(a)(b)Figure 11
Distribution of gas pressure under different anisotropy coefficients at 180 days. (a)λ=0.1, (b) λ=0.2, (c) λ=0.3, and (d) λ=0.5.Figure12 shows the evolution law of the matrix and fracture permeability at point M in stage II. The overall trend of the ratio of matrix and fracture permeability is polynomial increase with the time increases, and the increase of matrix permeability is greater than fracture permeability, and the increase of the fracture permeability ratio in x direction is greater than that in ydirection. As the anisotropy coefficient increases, the ratio of matrix and fracture permeability both increase significantly.Figure 12
Matrix permeability ratio and fracture permeability ratio under different anisotropy coefficients: (a) matrix permeability ratio and (b) fracture permeability ratio.
(a)(b)


## 3.2.1. Hydraulic Fracturing Simulation

The boundary around the geometry is set as the fixed boundary of the solid mechanical field and the impermeable boundary of the water seepage field, and the central hole is set as the water injection hole. And the following material parameters are adopted:E0=8.47GPa, v=0.23, ρs=1650kg/m3, σt=0.10MPa, σc=16.80MPa, ξ1=5, ξ2=5,φf0=0.018, and kf0=0.01mD. The water injection lasts for 120 minutes, and we select the water injection pressure pw0=1MPa, pw0=3MPa, pw0=5MPa, and pw0=10MPa for simulation analysis. Grid generation of the domain, including 8408 elements, is done using the free triangular mesh method.The state damage variableHd is used to record the damaged area of the coal seam, so the surface integral of Hd in the simulation domain could obtain the distribution of the damaged area. Figure 5 is the distribution of damaged area in coal seam under different water injection pressure at 5 min. As the material of the simulation domain follows the Weibull distribution, and the permeability is set as isotropic, the distribution of the damaged area is radially outward from the central hole. It can be seen from Figure 6 that damaged area size increases with the water pressure increases, showing a trend of gradually increasing first and then stabilizing, and the greater the water injection pressure, the longer it takes for the damaged area to develop.Figure 5
Distribution of the damaged area under different water injection pressures at 5 min. (a)pw0=1MPa, (b) pw0=3MPa, (c) pw0=5MPa, and (d) pw0=10MPa.
(a)(b)(c)(d)Figure 6
Damaged area size under different water injection pressures.Based on the geometry in Figure4, the proposed method is compared with the CwM method [18–21] to verify its precision. Figures 7 and 8 are the damaged area distribution cloud diagram and damaged area size variation curves, respectively, of the CwM method at different time steps Δt and the proposed method when the injection pressure pw0=5MPa. It can be seen from the damage area distribution that there are certain differences in the damage distribution of the two methods; at the same time, only the time step Δt=0.1min of the CwM method is similar to the proposed method. And the damaged area size is shown as CwM, Δt=0.1min<CwM, Δt=0.5min<Proposedmethod, Hd<CwM, Δt=1min. The damage area size of the CwM method increases with the time step Δt. When using the CwM method to simulate hydraulic fracturing, the total damage is obtained by the accumulation of the damage generated by the simulation within a unit time step. A longer time step will result in a larger cumulative error in the simulation. When the time step is smaller, the error can be reduced, but due to the time difference between the update of the damage variable and the variables of the stress field and seepage field, it is easy to cause incomplete development of the damaged area and underestimate the size. In addition, because the principle of this method is based on the mutual conversion between the matrix data points of MATLAB and the geometric grid in COMSOL, because the matrix data points and the grid are not in a one-to-one correspondence, there is a little error. However, the proposed method is carried out on the geometric grid in COMSOL, and the state damage variable Hd is updated in real time with the dependent variables such as water pressure and displacement when solving, which could avoid these problems to some extent.Figure 7
Distribution of damaged area in the CwM method and proposed method at 5 min. (a) CwM,Δt=1min; (b) CwM, Δt=0.5min; (c) CwM, Δt=0.1min; and (d) proposed method, Hd.
(a)(b)(c)(d)Figure 8
Damaged area size of the CwM method and proposed method.


## 3.2.2. Coalbed Methane Extraction Simulation after Hydraulic Fracturing

According to Section3.1, the coalbed methane extraction simulation after hydraulic fracturing can be divided into two stages. In stage I, the boundary conditions and parameters are similar to the hydraulic fracturing simulation, and the water injection pressure is set as pw0=5MPa. In stage II, the boundary around the geometry in Figure 4 is set as the fixed boundary of the solid mechanical field and the impermeable boundary of the methane seepage field, and the central hole is set as the drainage hole. The relevant parameters of coalbed methane extraction are as follows: pfg0=pm0=1.5MPa, τ=153360s, φm0=0.035, PL=0.93MPa, VL=37.04m3/t, km0=0.002mD, and kfx0=0.01mD. And the point M in the geometry is the permeability monitoring point. To analyze the effect of permeability anisotropy, we introduce an anisotropy coefficient λ to describe the evolution characteristics of coal seam fracture permeability, λ=kfy0/kfx0, λ∈(0,1). The anisotropy coefficients λ=0.1, λ=0.2, λ=0.3, and λ=0.5 are adopted for simulation analysis. The coalbed methane extraction simulation lasts for 180 days.It can be seen from Figure9 that damaged area size mainly develops in x direction of coal seam under different anisotropy coefficients. With increasing anisotropy coefficient λ, the damaged area size in xdirection decreases gradually, while the damaged area size in y direction increases gradually. And the total damaged area size increases with the increase of the anisotropy coefficient λ (Figure 10(a)). The distributions of gas pressure in Figure 11 show obvious anisotropic characteristics. The smaller the anisotropy coefficient λ is, the larger the difference of gas pressure distribution in x direction andydirection is. Figure 10 shows the damaged area size and gas production rate under different anisotropy coefficients. The development of the damaged area is affected by the anisotropy of the permeability. When λ=0.1, the development of the damaged area is slower, and the development time of the damage lasts about 100 minutes. And when λ=0.5, the development of the damaged area is faster and lasts about 100 minutes. This is because the permeability in ydirection in the model is small, and x direction is the main seepage direction of the injected water. The water pressure of the damaged area in a short time is not enough to cause the failure of the undamaged coal body. It requires continuous water injection to increase the water pressure of the damaged area further fractures. The coal seam gas production rate also shows a gradual decrease as time increases (Figure 10(b)). The permeability of the damaged area increases, and the gas is easier to be pumped out, while the gas in the undamaged area has a lower permeability and the gas is not easy to be pumped out.Figure 9
Distribution of damaged area under different anisotropy coefficients at 120 min. (a)λ=0.1, (b) λ=0.2, (c) λ=0.3, and (d) λ=0.5.
(a)(b)(c)(d)Figure 10
Damaged area size and gas production rate under different anisotropy coefficients. (a) Damaged area size; (b) Gas production rate.
(a)(b)Figure 11
Distribution of gas pressure under different anisotropy coefficients at 180 days. (a)λ=0.1, (b) λ=0.2, (c) λ=0.3, and (d) λ=0.5.Figure12 shows the evolution law of the matrix and fracture permeability at point M in stage II. The overall trend of the ratio of matrix and fracture permeability is polynomial increase with the time increases, and the increase of matrix permeability is greater than fracture permeability, and the increase of the fracture permeability ratio in x direction is greater than that in ydirection. As the anisotropy coefficient increases, the ratio of matrix and fracture permeability both increase significantly.Figure 12
Matrix permeability ratio and fracture permeability ratio under different anisotropy coefficients: (a) matrix permeability ratio and (b) fracture permeability ratio.
(a)(b)


## 4. Numerical Examples

### 4.1. Numerical Simulation Scheme

In order to further study the feasibility of the proposed method, this section carried out a 3D simulation of coalbed methane extraction after hydraulic fracturing based on the actual engineering background. Taking a single branch of a multibranch well as an example, a geometric model as shown in Figure13 is established to study the effect of coalbed methane extraction after hydraulic fracturing. The size of the physical model is 50 m ×50 m ×5.6 m, the vertical well with a height of 5.6 m is located in the middle of the model xy plane, and the horizontal branch well with a length of 20 m is located at 2.8 m from the coal seam floor, which is connected with the vertical well. The diameter of the gas well is 0.1 m, the buried depth of the coal seam is 320 m, the thickness is 5 m, and the gas pressure is 1.32 MPa. And other numerical simulation parameters are shown in Table 1.Figure 13
Geometric model for numerical simulation.Table 1
Parameters for numerical simulation.
VariableParametersValueVariableParametersValueEYoung’s modulus/GPa8.47ξ2Increased coefficient of permeability10νPoisson’s ratio0.23pfg0Initial methane pressure in fracture/MPa1.32ρsCoal density/kg/m31650pm0Initial methane pressure in matrix/MPa1.32σtUltimate tensile strength/MPa0.10φm0Initial porosity for coal matrix0.035σcUltimate compressive strength/MPa16.8φf0Initial porosity for fracture0.018θInternal friction angle/°32aCoal matrix size/m0.01p0Water injection pressure/MPa10TTemperature of coal seam/K300kfx0Initial fracture permeability inxdirection/m22.25×10−17μgMethane dynamic viscosity/pa·s1.84×10−5kfy0Initial fracture permeability inydirection/m21.125×10−17τAdsorption time/s153360kfz0Initial fracture permeability inz direction/m22.25×10−18pLLangmuir pressure constant/MPa0.93km0Initial matrix permeability/mD0.002VLLangmuir volume constant/m3/kg0.0371pcExtraction pressure/kPa20εLLangmuir strain constant0.012ξ1Increased coefficient of porosity10hKlinkenberg factor/MPa0.76The initial boundary conditions are as follows: the top of the model is affected by the stress of the overlying strata, the bottom is a fixed boundary, and the four sides are sliding boundaries. And the outer boundary of the physical model is the impermeable boundary of coalbed methane and water. The maximum permeability is inx direction. The AB line is the air pressure and water pressure monitoring line, which is in the same plane as the central axis of the horizontal well, and point C is a point on the AB line, which is used to monitor the evolution of permeability. The simulation scheme is divided into two stages. In stage I, the boundary of the horizontal well is set as the water injection pressure boundary, the injection pressure is pw0=10MPa, and the water injection lasts for 480 minutes. In stage II, there are two schemes of coalbed methane extraction before and after hydraulic fracturing. The boundaries of the horizontal well and the vertical well are set as the drainage boundary, and the extraction time is set as 1000 days.


### 4.2. Simulation Results of Hydraulic Fracturing

To intuitively describe the effect of coal seam hydraulic fracturing, we made two cross sections along the central axis of the horizontal well on thexy plane and yz plane of the model and plotted the damage area distribution and water pressure distribution cloud diagram as shown in Figure 14. It shows the distribution damaged area, and water pressure at fracturing time is 5 min, 50 min, and 480 min. In the initial stage of fracturing, the damaged area and water pressure are mainly distributed near the horizontal well. As time increases, the damaged area expands slowly, and the water pressure distribution has expanded from the horizontal well to the boundary of the model.Figure 14
Distribution of the damaged area and water pressure at different fracturing times: (a) distribution of the damaged area and (b) distribution of water pressure.In combination with Figure15(a), it can be seen that the expansion of the damaged area occurs mainly in the initial stage of fracturing. After about 120 minutes, the damaged area size tends to be generally stable, although there may be a small increase. The water pressure at point C increases with time, but its increasing amplitude gradually decreases. Figure 15(b) shows the water pressure distribution curve of line AB at different times. The 25 m distance is the junction of the horizontal well and the vertical well, where the water pressure is 10 MPa. Affected by hydraulic fracturing, the coal seam near here has been damaged, so that the water pressure at different times is higher near here than far away from here, and the water pressure decreases with the distance increasing. As the fracturing time increases, the water pressure on line AB increases significantly.Figure 15
Damaged area size and water pressure: (a) damaged area size and water pressure at point C and (b) water pressure of line AB.
(a)(b)


### 4.3. Simulation Results of Coalbed Methane Extraction

Based on the cloud diagram drawing method in Figure14, we added a section perpendicular to the horizontal well along the xz plane to analyze the gas pressure distribution of coalbed methane extraction before and after hydraulic fracturing, as shown in Figure 16. It can be seen from Figure 16 that after hydraulic fracturing, the gas pressure near the gas well is significantly lower than before fracturing, and as the extraction time increases, the low-pressure gas pressure area gradually extends outward, and the extension range is significantly larger than that before fracturing. In the early stage of coalbed methane extraction, the low-pressure area of the gas pressure is mainly distributed around the gas well. Affected by the anisotropy of permeability, the gas pressure distribution before and after fracturing extends along the direction of maximum permeability with increasing extraction time.Figure 16
Distribution of gas pressure before fracturing and after fracturing at different times: (a) distribution of gas pressure before fracturing and (b) distribution of gas pressure after fracturing.It can be seen from Figure17(a) that the gas pressure after hydraulic fracturing at different extraction times is significantly lower than before fracturing, and the closer the gas well is, the lower the gas pressure. When the extraction time is 1000 days, the gas pressure before fracturing increases sharply and then gradually decreases with the increase of the distance from the gas well. There is an inflection point at about 13 m and 37 m; the influence of coalbed methane extraction is gradually reduced. After fracturing, the gas pressure increases with increasing distance from the gas well, forming an obvious pressure drop funnel, and the range of extraction influence is significantly greater than before fracturing. Figure 17(b) shows that the cumulative production of coalbed methane after fracturing is significantly greater than that before fracturing. When the extraction time is about 700 days, the cumulative production after fracturing has a small inflection point. It can be seen that the impact of hydraulic fracturing is gradually reduced, and the gas production rate is gradually reduced.Figure 17
Gas pressure and cumulative production before and after fracturing: (a) gas pressure of line AB and (b) cumulative production of coalbed methane.
(a)(b)After hydraulic fracturing, the porosity and permeability of the damaged area increase significantly, which promotes the increase of coalbed methane gas production. For point C of the undamaged area before and after fracturing, it can be seen from Equations (15)–(16) and Equation (19) that the permeability of the matrix and fractures is mainly controlled by the initial stress and the adsorption/desorption stress of methane. Due to the increase in gas production after fracturing, the desorption of adsorbed methane increases, which indirectly leads to an increase in permeability. As can be seen from Figure 18, the overall trend of the permeability ratio variation curve is as follows: with the increase of extraction time, both matrix and fracture permeability show a polynomial increase, and the increase of matrix permeability is greater than that of fracture permeability. The permeability of matrix and fracture increases greatly after fracturing, which is the result of hydraulic fracturing and coalbed methane extraction. The variation trend of fracture permeability ratio in xand y directions is basically the same before hydraulic fracturing, indicating that permeability anisotropy has a similar influence on them, while there is a difference between them after hydraulic fracturing. It can be seen that hydraulic fracturing can not only increase coal seam permeability through fracturing damage but also indirectly increase coal seam permeability by increasing the production of coalbed methane.Figure 18
Matrix permeability ratio and fracture permeability ratio before fracturing and after fracturing: (a) matrix permeability ratio and (b) fracture permeability ratio.
(a)(b)


## 4.1. Numerical Simulation Scheme

In order to further study the feasibility of the proposed method, this section carried out a 3D simulation of coalbed methane extraction after hydraulic fracturing based on the actual engineering background. Taking a single branch of a multibranch well as an example, a geometric model as shown in Figure13 is established to study the effect of coalbed methane extraction after hydraulic fracturing. The size of the physical model is 50 m ×50 m ×5.6 m, the vertical well with a height of 5.6 m is located in the middle of the model xy plane, and the horizontal branch well with a length of 20 m is located at 2.8 m from the coal seam floor, which is connected with the vertical well. The diameter of the gas well is 0.1 m, the buried depth of the coal seam is 320 m, the thickness is 5 m, and the gas pressure is 1.32 MPa. And other numerical simulation parameters are shown in Table 1.Figure 13
Geometric model for numerical simulation.Table 1
Parameters for numerical simulation.
VariableParametersValueVariableParametersValueEYoung’s modulus/GPa8.47ξ2Increased coefficient of permeability10νPoisson’s ratio0.23pfg0Initial methane pressure in fracture/MPa1.32ρsCoal density/kg/m31650pm0Initial methane pressure in matrix/MPa1.32σtUltimate tensile strength/MPa0.10φm0Initial porosity for coal matrix0.035σcUltimate compressive strength/MPa16.8φf0Initial porosity for fracture0.018θInternal friction angle/°32aCoal matrix size/m0.01p0Water injection pressure/MPa10TTemperature of coal seam/K300kfx0Initial fracture permeability inxdirection/m22.25×10−17μgMethane dynamic viscosity/pa·s1.84×10−5kfy0Initial fracture permeability inydirection/m21.125×10−17τAdsorption time/s153360kfz0Initial fracture permeability inz direction/m22.25×10−18pLLangmuir pressure constant/MPa0.93km0Initial matrix permeability/mD0.002VLLangmuir volume constant/m3/kg0.0371pcExtraction pressure/kPa20εLLangmuir strain constant0.012ξ1Increased coefficient of porosity10hKlinkenberg factor/MPa0.76The initial boundary conditions are as follows: the top of the model is affected by the stress of the overlying strata, the bottom is a fixed boundary, and the four sides are sliding boundaries. And the outer boundary of the physical model is the impermeable boundary of coalbed methane and water. The maximum permeability is inx direction. The AB line is the air pressure and water pressure monitoring line, which is in the same plane as the central axis of the horizontal well, and point C is a point on the AB line, which is used to monitor the evolution of permeability. The simulation scheme is divided into two stages. In stage I, the boundary of the horizontal well is set as the water injection pressure boundary, the injection pressure is pw0=10MPa, and the water injection lasts for 480 minutes. In stage II, there are two schemes of coalbed methane extraction before and after hydraulic fracturing. The boundaries of the horizontal well and the vertical well are set as the drainage boundary, and the extraction time is set as 1000 days.


## 4.2. Simulation Results of Hydraulic Fracturing

To intuitively describe the effect of coal seam hydraulic fracturing, we made two cross sections along the central axis of the horizontal well on thexy plane and yz plane of the model and plotted the damage area distribution and water pressure distribution cloud diagram as shown in Figure 14. It shows the distribution damaged area, and water pressure at fracturing time is 5 min, 50 min, and 480 min. In the initial stage of fracturing, the damaged area and water pressure are mainly distributed near the horizontal well. As time increases, the damaged area expands slowly, and the water pressure distribution has expanded from the horizontal well to the boundary of the model.Figure 14
Distribution of the damaged area and water pressure at different fracturing times: (a) distribution of the damaged area and (b) distribution of water pressure.In combination with Figure15(a), it can be seen that the expansion of the damaged area occurs mainly in the initial stage of fracturing. After about 120 minutes, the damaged area size tends to be generally stable, although there may be a small increase. The water pressure at point C increases with time, but its increasing amplitude gradually decreases. Figure 15(b) shows the water pressure distribution curve of line AB at different times. The 25 m distance is the junction of the horizontal well and the vertical well, where the water pressure is 10 MPa. Affected by hydraulic fracturing, the coal seam near here has been damaged, so that the water pressure at different times is higher near here than far away from here, and the water pressure decreases with the distance increasing. As the fracturing time increases, the water pressure on line AB increases significantly.Figure 15
Damaged area size and water pressure: (a) damaged area size and water pressure at point C and (b) water pressure of line AB.
(a)(b)


## 4.3. Simulation Results of Coalbed Methane Extraction

Based on the cloud diagram drawing method in Figure14, we added a section perpendicular to the horizontal well along the xz plane to analyze the gas pressure distribution of coalbed methane extraction before and after hydraulic fracturing, as shown in Figure 16. It can be seen from Figure 16 that after hydraulic fracturing, the gas pressure near the gas well is significantly lower than before fracturing, and as the extraction time increases, the low-pressure gas pressure area gradually extends outward, and the extension range is significantly larger than that before fracturing. In the early stage of coalbed methane extraction, the low-pressure area of the gas pressure is mainly distributed around the gas well. Affected by the anisotropy of permeability, the gas pressure distribution before and after fracturing extends along the direction of maximum permeability with increasing extraction time.Figure 16
Distribution of gas pressure before fracturing and after fracturing at different times: (a) distribution of gas pressure before fracturing and (b) distribution of gas pressure after fracturing.It can be seen from Figure17(a) that the gas pressure after hydraulic fracturing at different extraction times is significantly lower than before fracturing, and the closer the gas well is, the lower the gas pressure. When the extraction time is 1000 days, the gas pressure before fracturing increases sharply and then gradually decreases with the increase of the distance from the gas well. There is an inflection point at about 13 m and 37 m; the influence of coalbed methane extraction is gradually reduced. After fracturing, the gas pressure increases with increasing distance from the gas well, forming an obvious pressure drop funnel, and the range of extraction influence is significantly greater than before fracturing. Figure 17(b) shows that the cumulative production of coalbed methane after fracturing is significantly greater than that before fracturing. When the extraction time is about 700 days, the cumulative production after fracturing has a small inflection point. It can be seen that the impact of hydraulic fracturing is gradually reduced, and the gas production rate is gradually reduced.Figure 17
Gas pressure and cumulative production before and after fracturing: (a) gas pressure of line AB and (b) cumulative production of coalbed methane.
(a)(b)After hydraulic fracturing, the porosity and permeability of the damaged area increase significantly, which promotes the increase of coalbed methane gas production. For point C of the undamaged area before and after fracturing, it can be seen from Equations (15)–(16) and Equation (19) that the permeability of the matrix and fractures is mainly controlled by the initial stress and the adsorption/desorption stress of methane. Due to the increase in gas production after fracturing, the desorption of adsorbed methane increases, which indirectly leads to an increase in permeability. As can be seen from Figure 18, the overall trend of the permeability ratio variation curve is as follows: with the increase of extraction time, both matrix and fracture permeability show a polynomial increase, and the increase of matrix permeability is greater than that of fracture permeability. The permeability of matrix and fracture increases greatly after fracturing, which is the result of hydraulic fracturing and coalbed methane extraction. The variation trend of fracture permeability ratio in xand y directions is basically the same before hydraulic fracturing, indicating that permeability anisotropy has a similar influence on them, while there is a difference between them after hydraulic fracturing. It can be seen that hydraulic fracturing can not only increase coal seam permeability through fracturing damage but also indirectly increase coal seam permeability by increasing the production of coalbed methane.Figure 18
Matrix permeability ratio and fracture permeability ratio before fracturing and after fracturing: (a) matrix permeability ratio and (b) fracture permeability ratio.
(a)(b)


## 5. Conclusion

(1)
An improved coupled modeling method for coalbed methane extraction after hydraulic fracturing is proposed. It mainly includes two parts: hydraulic fracturing modeling and coupling modeling for coalbed methane extraction. In the first part, we introduce a state damage variable based on damage mechanics theory to record the damage changes at different times and combine it with the classical porous elasticity theory to construct the coupling governing equations of the seepage field and the solid mechanics field in the hydraulic fracturing process. In the second part, based on the dual porosity and dual permeability coupling model of the microporous matrix and fractures in coal seams, we consider the influence of the degree of fracture opening in different directions on the evolution of permeability and construct the coupling model of coalbed methane extraction. The relationship between the two is as follows: the numerical simulation results of the first part, such as elastic modulus, damage variables, porosity, and permeability, are used as the initial values of the second part for simulation analysis(2)
The coupling modeling method is realized in COMSOL Multiphysics by adopting the solid mechanics module, Darcy’s law module, and coefficient form PDE module to solve the state damage variable, displacement, and pressure sequentially. And a 2D hydraulic fracturing simulation and a 2D coalbed methane extraction simulation after hydraulic fracturing verify the feasibility of the method. In addition, a 3D coalbed methane extraction simulation is used to analyze the effect of coalbed methane extraction before and after hydraulic fracturing and the evolution characteristics of permeability(3)
All numerical examples indicate that the proposed coupling modeling method can not only simplify the calculation steps of the CwM method and reduce the numerical errors caused, but also characterize the damage, porosity, permeability, and effect on coalbed methane extraction caused by hydraulic fracturing in coal seams with permeability anisotropy



---

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