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.

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 [36]. 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 [915]. 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 [1821] 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 [2931]. 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 Section 2. 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 where and are the initial elastic modulus of REVs before and after damage, MPa, and is the state damage variable, which is a historical state variable related to time and strain tensor. The relationship between and the damage-driven evolution function satisfies the Kuhn Tucker condition, which can be described as, where is the time, , and is the strain tensor at time .

It can be seen from Equation (2) that when , (0, ), 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, (0, ) = 0. When , is the larger value between the damage-driven evolution function at time and at time And the damage-driven evolution function is defined as where and are the maximum tensile and compressive principal strain, , ; and are the ultimate tensile strength and compressive strength, Pa; and are the first and third principal strains; and and 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: in which , , where , , and are the first, second, and third principal stress, Pa; is the internal friction angle, °; and 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: where is the material density, kg/m3; is the variable of the displacement field, ; is the time, ; is the Cauchy stress tensor, Pa; is the water Biot coefficient; is the water pressure, Pa; is the acceleration of gravity, m/s2; is the fracture porosity; is the compressibility factor, ; is the water flow velocity, m/s; and 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 follows Darcy’s law, and it is defined by the following: where is the water dynamic viscosity, Pa.s, and 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: where and are the initial fracture porosity and permeability; is the increased coefficient of porosity after damage; and is the increased coefficient of permeability after damage.

Substituting Equation (7) and Equation (8) into Equation (6), Equation (6) can be written as

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 () in different directions.

As the cuboid fracture model in Figure 1, the initial fracture porosity and permeability can be obtained by the following (Equation (19)) where is the matrix size, m; is the initial fracture permeability in direction, m2; i, j = x, y, z, represent the three directions of the Cartesian coordinate; and and are the opening degree of the fracture perpendicular to the directions and , m.

Taking direction as an example, the relationship between stress and strain per unit volume is where is the volumetric strain per unit volume in direction; and are the Langmuir strain constant and pressure constant; is the matrix methane pressure increment, Pa; is the improved fracture stiffness, Pa; is the bulk modulus, Pa; and and 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 ; that is, [33], and then, Equation (11) can be modified to

In coal seam, the fracture opening degree increment is closely related to the effective stress increment and the initial fracture opening degree, and its value can be obtained by

Then, the fracture porosity in the coal seam is

According to the permeability cubic law, the fracture permeability in direction is

Similarly, the permeability in and directions is where represents the fracture permeability in direction, m2, and and are the fracture opening degree and the opening degree increments perpendicular to direction, respectively, m.

The porosity of the coal matrix is defined as where where is the initial matrix porosity; is the matrix methane Biot coefficient; is the matrix bulk modulus, Pa; is the Poisson’s ratio; is the volumetric strain; and is the matrix methane pressure, Pa.

The permeability of the coal matrix can be described as where 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 where is the Langmuir volume constant, m3/kg; is the dynamic viscosity of methane, Pa.s; and are the temperature and the methane pressure under standard condition; is the temperature of the coal seam; is the methane desorption time, s; is the methane molar mass, kg/mol; is the molar constant; is the pressure methane in fracture, Pa; and is the Klinkenberg factor, Pa.

Similarly, based on the Navier equation, the stress field governing equation for coalbed methane migration can be modified to where is the shear modulus, Pa; is the fracture methane Biot coefficient; and 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 Section 2 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:

3.1.1. Stage I: Hydraulic Fracturing Stage

In this stage, the water pressure of the seepage field is applied to the mechanical field in the form of body load. Under the combined action of water pressure 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 . Meanwhile, the state damage variable records the damage at this moment and transmits it to the elastic modulus of the REV. As the elastic modulus of some element nodes weakened, the stress field will be adjusted. And the porosity and the permeability 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 variable , elastic modulus , porosity , and permeability 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.

Figure 3 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 ; (3) create 2D or 3D geometry; (4) set up a solid mechanics module, a Darcy’s law module, a coefficient form PDE (), and a coefficient form PDE (), 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 , state damage variable , porosity , and permeability ; (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).

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 Figure 4. 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.

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: , , , , , , , and . The water injection lasts for 120 minutes, and we select the water injection pressure , , , and for simulation analysis. Grid generation of the domain, including 8408 elements, is done using the free triangular mesh method.

The state damage variable is used to record the damaged area of the coal seam, so the surface integral of 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.

Based on the geometry in Figure 4, the proposed method is compared with the CwM method [1821] 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 and the proposed method when the injection pressure . 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 of the CwM method is similar to the proposed method. And the damaged area size is shown as CwM, , , , . The damage area size of the CwM method increases with the time step . 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 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.

3.2.2. Coalbed Methane Extraction Simulation after Hydraulic Fracturing

According to Section 3.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 . 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: , , , , , , and . 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, , λ∈(0,1). The anisotropy coefficients , , , and are adopted for simulation analysis. The coalbed methane extraction simulation lasts for 180 days.

It can be seen from Figure 9 that damaged area size mainly develops in direction of coal seam under different anisotropy coefficients. With increasing anisotropy coefficient , the damaged area size in direction decreases gradually, while the damaged area size in 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 direction anddirection 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 , the development of the damaged area is slower, and the development time of the damage lasts about 100 minutes. And when , the development of the damaged area is faster and lasts about 100 minutes. This is because the permeability in direction in the model is small, and 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 12 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 direction is greater than that in direction. As the anisotropy coefficient increases, the ratio of matrix and fracture permeability both increase significantly.

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 Figure 13 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.

The 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 in 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 , 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 the xy 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.

In combination with Figure 15(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.

4.3. Simulation Results of Coalbed Methane Extraction

Based on the cloud diagram drawing method in Figure 14, 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.

It can be seen from Figure 17(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.

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 and 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.

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

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This research was financially supported by the National Natural Science Foundation of China (Grant Nos. 52174117, 52004117), the Basic Research Project of Key Laboratory of Liaoning Provincial Education Department (Grant No. LJ2020JCL005), and the project supported by the Postdoctoral Science Foundation of China (Grant Nos. 2021T140290, 2020M680975). This article benefited from valuable comments and suggestions by co-editors and other anonymous reviewers.