#### Abstract

Risk prediction of dynamic disasters such as rock burst, gas outburst, and water inrush is closely related to the permeability evolution of coal seam. According to the characteristics of the lower protective layer mining, the basic assumption of gas-solid coupling model of the coal was proposed in this paper. The permeability enhancement coefficient of equivalent layer spacing was first put forward. Based on the three-zone-shaped dynamic evolution of the permeability of the overlying protective layer during the lower protective layer mining, the theories of seepage mechanics and damage mechanics were applied to introduce the permeability enhancement coefficient of equivalent layer spacing. A mathematical model of permeability evolution of the protected coal seam in the lower protective layer mining was established. Based on the engineering background of the lower protective layer mining in Changping Coal Mine, the numerical simulation using the proposed mathematical model was performed. The results showed that the stress and permeability of the protected layer in #3 coal seam evolved dynamically with the advancement of the working face of the protective layer in #8 coal seam. When the working face of the protective layer in #8 coal seam advanced to 80 m, the stress reduction rate in the relief area tended to be stable, and the stress in the stress reduction area was about 50% of the original rock stress. When advanced to 80 m, the permeability of the protected layer of the #3 coal seam increased sharply, and the permeability increased by 873 times. With the continuous advancement, the permeability of the protected layer in #3 coal seam tended to increase steadily, and the permeability increased by 1100–1200 times. The calculated magnitude of permeability increment is consistent with that in the engineering practice, indicating that the permeability evolution model is basically reasonable. The research provides a theoretical guidance for the gas drainage field application in the lower protective layer mining for prevention of coal and gas outburst.

#### 1. Introduction

Generally, there is a low permeability of coal seams in most of China’s mining areas. In the process of multicoal seam mining, the pressure relief and permeability increase in adjacent layers are widely used for gas extraction, which is an important direction of simultaneous exploitation of coal and gas. Coal reservoir is a dual porous rock composed of matrix porosity and fracture with its unique cleat system, including matrix pore system, fracture system, and skeleton deformation system. Matrix pore system is the main place of gas storage; fracture system and secondary fracture system, which have high permeability and little gas content, are the main pathways of gas movement.

Gray [1] first established the model which included the comprehensive effect of reservoir pressure and matrix contraction/expansion caused by adsorption in 1987. Palmer and Mansoori [2] proposed a permeability model based on strain effect, namely, P&M model. Pekot [3] proposed the ARI permeability model of coal-rock matrix strain corresponding to Langmuir reservoir pressure strain. Connell et al. [4] established the cubic and exponential expression for the permeability model of coal reservoir under the triaxial stress condition.

Li et al. [5] believed that the influence of the gas slippage effect should be considered in the study of permeability property under the change in effective stress. Li et al. [6] studied the effect of effective stress and coal matrix shrinkage effect on coal seam gas production and established a theoretical model for porosity and permeability of coal seam considering the coal matrix shrinkage effect. Based on the fractal theory of porous media, Fan et al. [7] established the mathematical model of oil-water relative permeability in the low permeability reservoir. Fan et al. [8] constructed a permeability model considering the comprehensive effect of water and gas pressure. Song et al. [3] deduced the expression of coal fracture strain under the coupling effect of mining disturbance and temperature and then constructed the permeability model of coal body in front of the working face under the coupling effect of mining disturbance and temperature. Liu et al. [9] developed the permeability model of hyperbolic function for deep coal body under the mining-induced stress. Li et al. [10] constructed the coal permeability model under the mining-induced stress in the present coal seam. Based on fractal theory and technology, Li et al. [11] put forward the permeability model of Newtonian fluid in porous media and analyzed the influence of microstructure parameters of porous media on the permeability. Based on the cube model of coal, Wang et al. [12] considered the influence of matrix porosity and slippage effect on permeability and established a new model for predicting the dynamic permeability of low rank coal. Besides, sensitivity analysis on the factors affecting absolute permeability and slippage coefficient was conducted, and the influence of methane and nitrogen gas on matrix shrinkage and slippage effect was analyzed.

In the lower protective layer mining, the protected layer is in different degrees of pressure relief deformation, which promotes the development of the pore-fracture system of the mining coal. As a result, the desorption of adsorbed gas is caused, and the permeability of the low permeable coal seam is increased, which is conducive to gas drainage. Recently, research on gas-solid coupling is mostly performed on the basis of Fick’s law and linear Darcy’s law, with the assumption that the unloaded coal is in the elastic stage and has small deformation. However, the difference in gas migration characteristics caused by the permeability change in different areas of the protected layer is not considered, and the objective nature of pressure relief gas flow is not fully studied. In this study, the influence of stress field on gas flow is considered, and the fluid-solid coupling model is deduced based on the constitutive equation and continuity equation of coal deformation to investigate permeability evolution. The numerical simulation using the proposed mathematical model was performed.

#### 2. Influence of Equivalent Layer Spacing on the Permeability of Coal

Equivalent layer spacing is defined as the ratio of the distance between the protective layer and the protected layer to the mining height of the protective layer. According to previous literature and field investigation, statistical analysis of the relationship between the mining parameters of the protective layer and the protection effect is performed. Mining heights of the protective layers are concluded as 1.5 m, 2.0 m, and 3.0 m. In the lower protective layer mining, the protective effect of the protected layer changes with the equivalent layer spacing, and the relationship is systematically analyzed under different mining heights, as shown in Table 1.

Figure 1 shows the protection effect for different mining heights of the protective layer. For the mining height of 1.5 m, the equivalent layer spacing between the protective layer and the protected layer is between 25 and 60 times, and the increased multiple of the permeability coefficient of the protected layer and the reduction rate of the original rock stress decrease linearly with the increase in the equivalent layer spacing. As shown in Figure 1, under a certain mining height, with the decrease in the equivalent layer spacing, the influence of mining on the protected layer increases gradually, and the fracture develops sufficiently so that a pathway is formed for gas flow of pressure-released coal. For the mining height of 2.0 m, with the increase in the equivalent layer spacing, the increased multiple of the permeability coefficient of the protected layer and the decrease rate of the original rock stress monotonously decrease. For the mining height of 3.0 m, with the increase in the equivalent layer spacing, the increased multiple of the permeability coefficient of the protected layer and the decrease rate of the original rock stress monotonously decrease.

The increased multiple of the permeability coefficient of the protected layer tends to stable with the increase in the equivalent layer spacing. This can be explained as follows: when the equivalent layer spacing increases to a certain extent, the influence of stope stress caused by protective layer mining on the protective layer is gradually reduced. Therefore, when the protective layer mining height is 3.0 m, the increased multiple of permeability coefficient has a smaller decrease than that at the mining height of 1.5 m and 2.0 m. When mining height of the protective layer is 1.5 m, 2.0 m, and 3.0 m, with the change in equivalent layer spacing, the permeability evolution relationship is shown in Figure 1.

In Figure 1, when the mining height of the protective layer is about 1.5 m and 2.0 m, the increased multiple of the permeability of the protected layer decreases monotonously with the increase in the equivalent layer spacing. The evolution function of permeability does not decrease monotonously when the protective layer mining height is 3.0 m. After the increase in equivalent layer spacing by about 35 times, the protected permeability shows a trend of increase, which is not consistent with the engineering practice. It is referred that the large equivalent layer spacing value of the protective layer may result in the data rationality distortion.

The above analysis shows that the equivalent layer spacing has an important influence on the permeability of the protected layer in the lower protective layer mining. To establish a reasonable mathematical model of permeability evolution of the protected layer in the lower protective layer mining, the influence of equivalent layer spacing on permeability should be considered in the modeling.

The permeability enhancement coefficient of equivalent layer spacing is defined and introduced in this paper. The physical definition is as follows: the permeability enhancement coefficient of equivalent layer spacing is constantly greater than 0, so as to describe the comprehensive influence of the interval (between the protective layer and the protected layer) and the mining thickness of the protective layer on the permeability of the protected layer during the lower protective layer mining.

As shown in Figure 1, the equivalent layer spacing and the permeability of the protected layer present reduction relation of quadratic function. Under the conditions of three different mining heights, the same terms in this function are combined, and then, coefficients of each factor are averaged to obtain quadratic function of the equivalent layer spacing and permeability evolution in the lower protective layer mining. Due to the distortion of the quadratic function relationship at the height of 3.0 m, the function between the equivalent layer spacing and permeability at 3.0 m is eliminated. Therefore, in the lower protective layer mining, the general function between the equivalent layer spacing and the permeability of the protected layer is as follows:where is the permeability enhancement coefficient of equivalent layer spacing and is the equivalent layer spacing between the protective layer and the protected layer in the lower protective layer mining.

#### 3. Permeability Evolution Model for Mining Disturbed Coal and Rock

##### 3.1. Basic Assumption

With the exploitation of the protective layer, the original balance of rock stress in the protected layer is broken by the mining. As a result, the coal body of the protected layer has different pressure relief deformations or even fractures; the pore pressure is reduced, and part of the adsorbed gas in the matrix pore system is converted into free gas, which flows into the fracture system in the form of seepage and diffusion. Free gas in the fracture system flows into the mining fracture, resulting in a decrease in the gas content in the protected layer. At the same time, the gas content in the protective layer is reduced, which has an impact on the stress balance of the coal in the protected layer. In other words, there is a strong dynamic coupling between the coal seam deformation and gas migration. In this paper, the following assumptions are made:(1)Coal is an isotropic dual pore structure.(2)There is only gas in the gas-bearing coal body, and the influence of water is not considered.(3)Gas is an ideal gas, and the temperature in the seepage field is constant.(4)Gas-bearing coal is a linear elastic material, and the mining-induced stress-strain relationship conforms to the elastic constitutive model of mining.(5)Gas-solid coupling system of coal and gas is composed of coal-rock solid skeleton and unidirectional gas.(6)The gas flow in fracture system accords with Darcy law.(7)The adsorbed gas content follows the Langmuir equation.(8)The gas desorbed in the pressure-released coal is extracted timely during the lower protective layer mining, and the fracture space of the protected layer is large enough for the lower protective layer. Therefore, the gas pressure of the fracture system is not considered in the stress field and deformation field analysis of the coal.(9)The influence of adsorption strain on fracture volume is not considered.

##### 3.2. Governing Equations of Stress-Strain Relation

The equations of stress field and deformation field of gas-bearing coal affected by mining are composed of four parts: stress balance equation, geometric equation, constitutive equation, and continuity equation. According to the characteristics of stress and fracture development in the lower protected layer mining [31], the permeability of the protected layer is divided into “three areas,” as shown in Figure 2. The application scope of gas-solid coupling equation is limited in the dotted box in Figure 2.

In Figure 3, differential elements are taken in the pressure-released coal of the protected layer in #3 coal seam. The average normal stress acting on the back is . We can get the force balance equation in the *x* direction:where is the normal stress, Pa; *τ* is the shear stress, Pa; and *F*_{x} is the force of gravity in *x* direction, Pa.

In the same way, the differential equations in three directions of the differential element are finally obtained:

In the protective layer mining project, pressure relief gas drainage project is usually equipped. Therefore, it is assumed that the gas desorbed in the pressure-released coal is extracted in time, and for the lower protective layer, the fracture space of the protected layer is large enough; the time of matrix gas desorption diffusion to the fracture system is slow, and the free gas pressure in the fracture system in the pressure relief coal body is slightly released. Therefore, the total stress acting on the pressure relief coal is mainly the skeleton stress of coal and the gas pressure in the matrix pore.

According to the modified Terzaghi effective stress principle, the equilibrium equation of elastic deformation of coal skeleton is obtained as follows [17]:where is the total stress of gas coal body, MPa; is the Kronecker symbol; *p*_{m} is the pressure of coalbed methane in the matrix, MPa; is the Biot effective stress coefficient corresponding to the pore; is the volume force, Pa; and *i*, *j* = *x*, *y*, *z*.

In the lower protective layer, the relationship between the strain and displacement of the pressure-released coal is described by the geometric equation as follows:where is the strain tensor of coal seam and and are displacement components of coal body, *m*.

The linear elastic constitutive relation of the solid skeleton of the coal follows the generalized Hooke’s law. If coal is regarded as the isotropic material, the elastic constitutive relation is expressed as follows:where is the effective stress, Pa; *G* is the shear modulus, Pa; is the coal elastic modulus, Pa; K is the coal bulk modulus, Pa; is Poisson’s ratio; and

Coal is regarded as a nonuniform elastic medium with dual structure of matrix and fracture, as well as dual permeability. Coal is affected by matrix pore pressure, gas adsorption/desorption, and overburden stress. Then, the stress field control equation is as follows:where is the Kronecker factor; is the adsorption strain, ; and is the corresponding gas pressure at 0.5.

The continuity equation of coal skeleton, namely, the mass conservation equation of skeleton is as follows:where is the skeleton density of coal, kg/m^{3}; is the density of coal body; is the porosity of coal matrix; is the porosity of coal fracture; and is the velocity vector of coal and rock skeleton, m/s.

The pore-fracture equation of the dual-porosity medium is as follows:

The continuity equation of coal skeleton can be obtained:where is the volume strain of coal skeleton.

##### 3.3. Analyses on Porosity and Permeability of Matrix System

Under mining conditions, mining stress has an important influence on the deformation of coal body and gas adsorption-desorption. Mining stress can cause the expansion and contraction of pores in the coal matrix system so that the porosity changes. Mining fracture provides sufficient space for gas desorption and makes the gas pressure in the original coal seam sharply reduce, so as to promote gas desorption. According to the assumption, the model is simplified without considering the influence of temperature on the adsorption stress and strain; the matrix permeability is used to deduce and optimize the matrix permeability model [32]. According to previous research [33], the matrix porosity can be expressed aswhere and .

Permeability in the pore medium is related to the size of porosity, and the relationship between permeability and porosity in the pore medium is as follows:where is the permeability of matrix system, ; is the effective diameter of coal matrix particles, *m*; and is the permeability of matrix pore system.

The change in matrix permeability is

##### 3.4. Analyses on Fracture Degree and Permeability of Fracture System

The porosity and permeability of fracture system are mainly controlled by the pore size of the model. To analyze the change in fracture aperture, representative elements are selected from the two porosity models, as shown in Figure 4. The square filled with solid at the center is a matrix block, and the fracture surrounds the matrix body. The solid line is the opening contour before deformation, and the dotted line is the contour after deformation.

The deformation of matrix and fracture is defined as follows:where *ε*_{eff} is the total strain in a certain direction; *ε*_{m} and *ε*_{f} are the volume strain of matrix system and fracture system, respectively; _{b} is the width of fracture, *m*; and *a* is the width of matrix block, *m*; *d* = +*b*. Based on the assumption of the double porosity medium model, *b* << *a* is obtained; therefore, *d* ≈ *a*.

According to the above assumptions and definitions, the total strain of the protected layer can be obtained:where is the permeability enhancement coefficient of the distance between the protective layer and the protected layer in the lower protective layer mining, and = 0.1456*H*_{R}^{2} − 19.7525*H*_{R} + 725.7; *H*_{R} is the distance between the protective layer and the protected layer under the lower protective layer mining; △*p*_{m} is the matrix pressure increment, MPa; is the Langmuir adsorption volume strain coefficient; *K*_{f} is the improved crack rigidity and *K*_{f} = *aK*_{n}; *K*_{n} is the fracture rigidity, MPa; *K*_{m} is the matrix rigidity, MPa; and is the effective stress increment of matrix and fracture, respectively; and *P*_{L} is the Langmuir stress coefficient, MPa.

Under the condition of equilibrium mechanics, the total stress is the same. Assuming that the effective stress increment of matrix and fracture is the same as the total effective stress [34], i.e., , equation (15) can be modified as follows:

The effective stress increment can be deduced from equation (16):

The change in effective stress is mainly caused by the adsorption effect. According to the change in effective stress, the increment of fracture width can be obtained as follows [35]:

According to the definition of fracture porosity, the change in fracture porosity can be written as follows:

According to the relationship between porosity and permeability, the change in permeability of fracture system is as follows:

The permeability of the overlying protective layer is in a dynamic process in the lower protective layer mining. Therefore, based on research results of mining stress and permeability, the evolution model of the mining permeability of the lower protective layer is established in line with the high permeability area, transition area, and low permeability area. In the high permeability area, coal is fully decompressed, the permeability is mainly determined by the fracture system, and then, the matrix permeability of the decompressed coal can be ignored. In the transition area, coal is in the area of stress concentration and rich longitudinal fracture development, which belongs to the scope of full pressure relief. Thus, the permeability of the decompressed coal in this area is mainly determined by the permeability of the matrix system and the permeability of the fracture system. In the low permeability area, coal is affected by the peak value of mining stress, and the coal is compressed and the fracture is closed. The permeability of the coal in this area is mainly determined by the permeability of the matrix system, and the permeability of the fracture system can be ignored. Therefore, the mathematical model of permeability of the protected layer in the lower protective layer mining is established as follows [36]:where *k*_{h} is the permeability within the high permeability area, m^{2}; *k*_{t} is the permeability within the transition permeability area, m^{2}; and *k*_{l} is the permeability within the low permeability area, m^{2}.

#### 4. Numerical Simulation of Permeability Evolution of the Protected Layer

##### 4.1. Model Parameters

By using COMSOL multiphysics numerical simulation software, the mathematical model of permeability evolution and gas flow continuity equation of mining coal in the lower protective layer mining is developed. The prevention and control of gas disaster of #3 coal seam in Changping Mine is taken as the engineering background [31], and the numerical calculation model is established, with the geometric dimension of 600 m × 400 m × 70 m. The mesh is encrypted at the interface and boundary between layers. The model and mesh are shown in Figure 5.

Through the results of numerical model, the adaptability and rationality of permeability evolution model of a lower protective layer are verified. The mining thickness of #8 protective layer is 2 m and that of #3 protected layer is 5.6 m. The model is mined from left to right. Considering the boundary effect, 50 m coal pillars are reserved on the left and right sides, the length of working face is 500 m, and the excavation step distance is 20 m. In order to ensure the convergence of the calculation, the coal strata in the simulation area are simplified. Table 2 shows the mechanical and other physical parameters of the simulated protective layer, the protected layer, and the strata.

##### 4.2. Initial and Boundary Conditions

The initial gas pressure in #8 protective layer is 0.42 MPa, the initial gas pressure in #3 protected layer is 0.55 MPa, the initial time is *t* = 0 s, and the simulation time is 30 d.

The left and right ends of the model are roller supports, the lower end is fixed boundary, and the unit area of the upper boundary of the model is applied with a stable load of 10.0 MPa. The gas pressure at both ends of the #8 protected layer is 0.42 MPa, and the gas pressure at the right end of the #3 protective layer is 0.55 MPa. The working face side gas pressure is 0.1 MPa standard atmospheric pressure, no gas is considered in each rock strata, and the gas pressure is 0.1 MPa.

##### 4.3. Evolution Law of Stress and Permeability of #3 Coal Seam in the Protected Layer

As shown in Figure 6, with the continuous advancement of the #8 protective layer, the stress of the #3 coal seam decreases and increases in varying degrees, and the permeability also evolves dynamically. When the working face of #8 protective layer is advanced to 40 m, the maximum stress of #3 coal seam is 9.7 MPa, the minimum stress is 7.0 MPa, and the maximum stress of local area is reduced by 22%; the maximum permeability is about 30 × 10^{−15} m^{2}, and the maximum permeability of #3 coal seam affected by mining is 6 times of the initial permeability, as shown in Figure 6(a).

**(a)**

**(b)**

**(c)**

**(d)**

When the working face of the protective layer in #8 coal seam is advanced to 80 m, the maximum stress of #3 coal seam is 10.6 MPa, the minimum stress is 5.0 MPa, and the maximum stress of local area is reduced by 44%; the maximum permeability is about 400 × 10^{−15} m^{2}, and the maximum permeability of #3 coal seam affected by mining is 368 times of the initial permeability, as shown in Figure 6(b).

When the working face of the protective layer in #8 coal seam is advanced to 120 m, the maximum stress of #3 coal seam is 12.0 MPa, the minimum stress is about 4.7 MPa, and the maximum stress of local area is reduced by 48%; the maximum permeability is about 1100 × 10^{−15} m^{2}, and the maximum permeability of #3 coal seam affected by mining is 733 times of the initial permeability, as shown in Figure 6(c).

When the working face of the protective layer in #8 coal seam is advanced to 480 m, the maximum mining stress is almost unchanged and the distribution law is stable; the maximum permeability is about 1840 × 10^{−15} m^{2}, and the maximum permeability of #3 coal seam affected by mining is 1232 times of the initial permeability, as shown in Figure 6(d).

When the working face of #8 coal seam is advancing, the permeability evolution rule of #3 coal seam is shown in Figure 7. It can be seen that when the working face is advancing to 80 m, the permeability of #3 coal seam increases sharply, which is caused by the initial pressure of the working face top of #8 coal seam and the development of mining fissures of #3 coal seam. When the working face advances within 120–160 m, the permeability increase in #3 coal seam slows down, and the permeability increases multiple is about 873 times. When the working face advances 200 m, the permeability increase in #3 coal seam slows down, and the permeability increase multiple is about 1000 times. This is caused by the square working face and the violent overburden movement. After that, with the continuous advance of #8 coal seam, the permeability of #3 coal seam tends to be stable, and the permeability increases about 1100–1200 times.

In the advancing of the working face of the protective layer in #8 coal seam, the evolution rule of mining stress of the protected layer in #3 coal seam is shown in Figure 8. It can be seen that, with the advancing of #8 coal seam working face, the mining stress of protected layer in #3 coal seam evolves dynamically; when the working face of protective layer in #8 coal seam advances to 80 m, the stress reduction rate in the pressure relief zone tends to be stable, and the stress in the stress reduction zone is about 50% of the original rock stress, and the relief effect is obvious.

Using the numerical simulation method, the distribution characteristics of stress field and permeability characteristics of the protected layer are shown in Figures 9 and 10 . Numerical simulation shows that the permeability evolution model of the protected layer is basically reasonable, and the calculated permeability characteristics and increase multiple are basically consistent with that in the engineering practice. The increased multiple of permeability obtained by numerical calculation is slightly larger than those obtained by engineering practice, but they are consistent in the magnitude. The reason for the larger increased multiple of permeability in numerical simulation is explained as follows: the mined-out area will be recompacted, while this cannot be reflected in the numerical simulation. In fact, with the increase in time, coal seam #3 is compacted by the overburden subsidence. As a result, the permeability of coal seam in the high permeability area gradually decreases and finally tends to the initial permeability.

#### 5. Conclusions

(1)The basic assumption of fluid-solid coupling modeling of mining coal is put forward for coal seam protected by the lower layer mining. The permeability enhancement coefficient of equivalent layer spacing is first proposed with the physical meaning. Through the data statistical analysis, the expression of the permeability enhancement coefficient of equivalent layer spacing under the lower protective layer mining is established.(2)Based on the dynamic evolution characteristics of the permeability of the overlying protective layer in the lower protective layer mining, a mathematical model of permeability evolution partition of the protected layer is established by introducing the permeability enhancement coefficient of equivalent layer spacing.(3)The results show that, in the process of advancing the working face of the lower protective layer, the permeability of the protected layer evolves dynamically. When the working face of #8 coal seam is advanced to 80 m, the permeability of the protective layer of #3 coal seam increases sharply, and the permeability increases by 873 times. When the working face is advanced to 200 m, the permeability increase in #3 coal seam slows down, and the permeability increases by 1000 times. After that, as the working face of #8 coal seam continues to advance, the permeability of #3 coal seam tends to be stable, and the permeability increases about 1100–1200 times.(4)The numerical simulation shows that the permeability evolution model of the protected layer is basically reasonable, and the calculated permeability characteristics and increase multiple are consistent with that in engineering practice.

The mining of the protective layer is time dependent, which should be considered in the study of permeability evolution in the future.

#### Nomenclature

a: | Width of matrix block |

_{b}: | Width of fracture |

d _{e}: | Effective diameter of coal matrix particles |

E: | Coal elastic modulus |

F _{i}: | Volume force |

F _{x}: | Force of gravity in x direction |

H _{R}: | Equivalent layer spacing between the protective layer and the protected layer in the lower protective layer mining |

K: | Coal bulk modulus |

K _{f}: | Improved crack rigidity |

k _{h}: | Permeability within the high permeability area |

k _{l}: | Permeability within the low permeability area |

K _{m}: | Matrix rigidity |

k _{m}: | Permeability of matrix system |

K _{n}: | Fracture rigidity |

k _{t}: | Permeability within the transition permeability area |

P _{L}: | Langmuir stress coefficient |

p _{l1}: | Corresponding gas pressure at 0.5ε_{l} |

p _{m}: | Pressure of coalbed methane in the matrix |

u _{ij}: | Displacement components of coal body |

α _{m}: | Biot effective stress coefficient |

δ _{ij}: | Kronecker factor |

ε _{eff}: | Total strain |

ε _{ij}: | Strain tensor of coal seam |

ε _{L}: | Langmuir adsorption volume strain coefficient |

ε _{m}: | Volume strain of matrix system |

ε _{f}: | Volume strain of fracture system |

ε _{s}: | Adsorption strain |

: | Volume strain of coal skeleton |

: | Lame constant |

G: | Shear modulus |

: | Permeability enhancement coefficient of equivalent layer spacing |

: | Skeleton density of coal |

ρ _{s}: | Density of coal body |

σ: | Normal stress |

σ _{ij}: | Effective stress |

σ _{ij}: | Total stress of gas coal body |

: | Shear stress |

: | Poisson’s ratio |

: | Velocity vector of coal and rock skeleton |

φ _{f}: | Porosity of coal fracture |

φ _{m}: | Permeability of matrix pore system |

*Subscript*

0: | Initial value of the variable |

m: | Matrix |

f: | Fracture. |

#### Data Availability

The data used in the article are from the simulation results, and we wrote the equations as codes independently to conduct these simulations. More information is available from the corresponding author on reasonable request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest with respect to the research, authorship, and/or publication of this article.

#### Acknowledgments

This research was financially supported by the National Natural Science Foundation of China (Grant no. 51504127), Science and Technology Research Support Project of Liaoning Provincial Department of Education (Grant no. LJ2019FL007), and Shanxi Province Science and technology plan announced bidding project (Grant no. 20191101015).