Numerical simulation method for rock mass matrix-fracture-roadway gas diffusion

Abstract

The invention relates to the field of underground resource mining, in particular to a numerical simulation method for rock mass matrix-fracture-roadway gas diffusion. A numerical simulation method for rock mass matrix-fracture-roadway gas diffusion comprises the following steps: (1) defining a boundary range of a calculation area according to an actual stratum structure and a roadway geometrical shape provided by a GIS (Geographic Information System), and meanwhile, determining accurate positions of a fracture and a roadway, width parameters of the fracture and cross section sizes of all points of the roadway in the axial direction; meanwhile, the permeability attribute of the rock matrix is determined by combining the physical characteristics and the spatial distribution rule of the rock matrix; the method disclosed by the invention has the beneficial effects that a complex interaction mechanism among the matrix, the crack and the roadway is disclosed, and the behavior and the distribution rule of gas diffusion under different conditions can be more accurately predicted; the gas diffusion process in the complex geologic structure is simulated in detail through a high-precision mathematical model and algorithm, and meanwhile, numerical simulation has higher efficiency and lower cost.

Description

Technical Field

[0001] The invention relates to the field of underground resource mining, and in particular to a numerical simulation method for rock matrix-crack-tunnel gas diffusion. Background Art

[0002] With the development and utilization of underground resources, safety and environmental issues such as gas leakage in tunnels and rock matrix during mining are becoming increasingly prominent. With the improvement of people's scientific and technological level, coalbed methane has become an important energy source during mining. Due to the complexity of the mining process and the region, the use of mathematical models for numerical simulation is currently an important means to analyze the mass transfer process of fractures, rock mass and tunnels. At present, the numerical simulation method of gas diffusion process is often only the analysis under the rock matrix-fracture or fracture-tunnel model, ignoring the complete gas diffusion of rock matrix-fracture-tunnel in stratum mining, resulting in the inability to accurately and comprehensively judge the gas leakage risk in the mining process and the inability to efficiently and reasonably utilize gas resources. The main reason is that the process has: cross-dimensional coupling of three-dimensional strata, two-dimensional fracture surface and one-dimensional simplified tunnel; there is seepage and convection coupling between rock matrix and tunnel, so it is difficult to form a unified cross-dimensional solution. Comprehensively improve the accuracy of gas diffusion process simulation, optimize resource utilization and effectively evaluate gas leakage risks. Summary of the invention

[0003] In order to comprehensively simulate the underground gas diffusion phenomenon, accurately predict the seepage field distribution, deeply reveal the interaction mechanism between fluid and solid medium, effectively evaluate environmental risks and optimize related engineering designs, this paper proposes a coupled numerical simulation method for matrix, fracture and tunnel gas diffusion. This method can accurately grasp the specific situation of gas diffusion by analyzing the pressure distribution and flow dynamic information inside each area according to the actual mineral mining conditions. The specific scheme is as follows: A numerical simulation method for rock matrix-crack-tunnel gas diffusion comprises the following steps: (1) Based on the geological model and tunnel geometry provided by the Geographic Information System (GIS) platform, the boundary range of the calculation area is defined, and the precise location of the fractures and tunnels, the width parameters of the fractures, and the cross-sectional dimensions of the tunnels at each point along their axis are determined; at the same time, the permeability properties of the rock matrix are determined in combination with its physical properties and spatial distribution. In addition, the boundary conditions and initial conditions required for the simulation are accurately obtained with the help of geological survey data and field monitoring data; (2) The geological model in (1) is converted into a high-precision grid model covering the entire calculation area; and low-dimensional cells are embedded inside the stratum to accurately simulate the geometric characteristics of the fracture surface; (3) Determine the gas control equations for the three systems of rock matrix, fractures, and tunnels; (4) Use the discrete fracture model to numerically simulate the fracture system, and use Darcy or non-Darcy seepage and skin coefficient to numerically simulate the fracture and matrix model according to the needs; (5) The gas diffusion in the tunnel is simplified into a one-dimensional convection-diffusion model with a rectangular cross section for numerical simulation; the state equations of the gas in the matrix-fracture-tunnel system are coupled to obtain the gas concentration in the tunnel; (6) Compare the gas concentration at different times in the tunnel obtained by numerical simulation with the actual gas threshold to conduct disaster warning in the formation, and visualize the simulation results of the matrix-fracture-tunnel.

[0004] In step (2), in the process of converting the preprocessed geometric model into a high-precision grid model covering the entire calculation area, two-dimensional units are embedded in the formation to accurately simulate the geometric characteristics of the fracture surface. At the same time, unique material property numbers are assigned to the fracture area and different types of porous media formations to facilitate the subsequent implementation of differentiated seepage calculation strategies to ensure the high accuracy and applicability of the simulation process; In step (4), when using the discrete fracture model to numerically simulate the fracture system, the fracture network needs to be regarded as the weak part inside the rock mass, which is characterized by high permeability and low mechanical strength. These characteristics must be fully considered in the simulation analysis. For this reason, a two-dimensional unit without thickness is used as a calculation tool to accurately capture the geometric and physical characteristics of the fracture. This method shows significant advantages in simulating thin-layer structures such as fractures, fine cracks, and interface transition areas. In step (5), for cross-dimensional couplings such as three-dimensional matrix to two-dimensional fracture, two-dimensional fracture to one-dimensional tunnel, the principle of equal pressure and flow superposition is used to establish boundary conditions and source terms. The specific schemes include matrix-fracture coupling, matrix-tunnel coupling, and fracture-tunnel coupling.

[0005] Preferably, the three-dimensional seepage model of the rock matrix in step (3) is performed under the default gas saturation state, and the seepage field control equation is as follows: (1) in ε p is the porosity of the matrix, dimensionless, u g Represents the Darcy velocity of the fluid in the porous medium, unit: m / s, ρ g Indicates the density of gas, unit: kg / m 3 , t Indicates time, unit: s, q m is the volume source term of the matrix, unit: 1 / s, ▽ is the gradient; The governing equation of the fracture two-dimensional seepage model described in step (3) is as follows: (2) in ε f is the porosity of the fracture, dimensionless, u f Represents the Darcy velocity of the fluid in the fracture, unit: m / s, q f is the volume source term of the crack, unit: 1 / s; The governing equation of the one-dimensional seepage model of the roadway described in step (3) is as follows: (3) in C is the gas concentration, unit: mol / m³, K is the diffusion coefficient, dimensionless, u t It is the Darcy velocity / average velocity of gas in the tunnel, unit: m / s.

[0006] When using discrete fracture models to numerically simulate fracture systems, the fracture network must be considered as a weak part of the rock mass, characterized by high permeability and low mechanical strength. These characteristics must be fully considered in the simulation analysis. To this end, a two-dimensional unit with no thickness is used as a computational tool to accurately capture the geometric and physical properties of the fractures. This method shows significant advantages in simulating thin-layer structures such as fractures, fine cracks, and interface transition areas.

[0007] Preferably, the coupling of the three systems in step (5) includes: coupling between the matrix and the cracks, coupling between the matrix and the tunnels, and coupling between the cracks and the tunnels. The specific coupling process is as follows: 1) Coupling process between matrix and cracks: The volume source term of the two-dimensional fracture surface is calculated by the flow through the fracture surface, that is, by q m and q f To couple the matrix with the cracks, q f The calculation formula is as follows: (4) in n 1 is the crack normal direction, k f is the permeability of the fracture, unit: m 2 , u fp is the flow rate of gas from the matrix into the cracks, unit: m / s,A f is the crack surface area, unit: m 2 , V f is the volume of the crack, unit: m 3 , P Represents the pressure of the crack gas, unit: Pa, μ Represents the fluid dynamic viscosity coefficient, unit: Pa·s. For two seepage fields, the two are source / sink terms for each other, so the volume flow rates of the two are inverse to each other; therefore, the source term should be satisfied q f=- q m , to satisfy the coupling between matrix and cracks; 2) Coupling process between matrix and lane: According to the average Darcy velocity of the tunnel wall seepage field, the boundary conditions of the tunnel gas diffusion velocity are established, and the flow velocity entering the tunnel is u tm The calculation formula is: (5) in n t is the inner normal unit vector of the tunnel wall, u p is the Darcy velocity of the matrix at the tunnel surface, in m / s, from the matrix to the inside of the tunnel; therefore, the average flow velocity entering the tunnel u t1 for: (6) In the formula is the wall area per unit length of the tunnel; 3) Coupling process between cracks and tunnels: If the crack happens to be connected to the tunnel: determine the cross-sectional area of ​​the crack by the crack width and length of the tunnel wall, and calculate the cross-sectional area of ​​the crack according to the flow rate of gas in the crack. u tf As the velocity boundary condition for tunnel gas diffusion, the calculation formula is: (7) In the formula u f represents the Darcy velocity of the fluid in the fracture, and the average flow velocity entering the tunnel is u t2 for: (8) In the formula It is the crack length per unit length of tunnel, unit: m.

[0008] Preferably, the calculation formula for the Darcy velocity u of the fluid in the porous medium is: (9) in k is the permeability, unit: m 2 ,use k p Expressed as the permeability of the rock matrix, unit: m 2 , k f Expressed as the permeability of the fracture, unit: m 2 , μ Represents the dynamic viscosity coefficient of the fluid, unit: Pa·s, p is the gas pressure, unit: Pa, ▽p is the pressure gradient, unit: Pa / m; in, u g Represents the Darcy velocity of gas in porous media, unit: m / s, u f Represents the Darcy velocity at the crack surface, unit: m / s, u t is the Darcy velocity of gas in the tunnel, unit: m / s, all are applicable to this calculation; Permeability of rock matrix k p It can be derived from the Kozeny-Carman relationship: (10) in, C 1 is the Kozeny-Carman constant, s For specific surface, in a packed bed or regular structure of granular soil, the rock matrix permeability can be derived from the Kozeny-Carman formula: (11) in, d p is the effective particle size of the particles that make up the matrix, unit: m; Fracture permeability k f It is expressed as: (12) in d f is the crack thickness, unit: m; Will u t2 as well as u t1， Substitute them into formula (3) to obtain the gas concentration in the tunnel under the two conditions.

[0009] Combining formula (2), formula (4), formula (9), and formula (12), we can obtain the gas pressure P in the cracks. Substituting the obtained pressure into formula (4), we can obtain the flow rate of gas flowing from the matrix into the cracks: u fp , further substituted into formula (2) to obtain the Darcy flow velocity in the fracture u f , further based on q f=— q m The Darcy velocity in the matrix can be obtained by formula (1): u g ;get u fp Then, substituting into formula (8) we can obtain u t2， By formula (9) u p, Further gain u t1 ， Will u t2 as well as u t1 ， Substitute them into formula (3) to obtain the gas concentration in the tunnel under the two conditions.

[0010] Preferably, the skin factor is also used S k Revised permeability, skin factor S k The relationship with permeability is: (13) in k , k d are the formation permeabilities before and after damage, unit: m 2 ; r e Indicates the radius of the contaminated area of ​​the nearby strata, unit: m, relative to the center of the fracture or tunnel; r w Indicates the area radius of the tunnel or crack, unit: m; in k is the formation permeability before damage, and the following formula is obtained: (14) in β is the inertial drag coefficient, unit: 1 / m, c Fis the Forchheimer dimensionless parameter, where Forchheimer is Forchheimer.

[0011] When calculating gas leakage in this scheme, the dispersion of clay, the existence of mud cakes and the blocking of pores by fine particles lead to the deterioration of the permeability of nearby formations, which makes the gas encounter additional resistance during the seepage process, thereby generating additional pressure drop, thus causing the leakage to the cracks or tunnels to weaken. The reason for this change is called the skin effect.

[0012] The dimensionless parameter that indicates the size of the skin effect is called the skin coefficient S k The skin coefficient is an important parameter, which reflects the influence of the change of permeability of nearby formations on seepage. The permeability can be revised through the skin coefficient.

[0013] Preferably, in high Reynolds number flow, a nonlinear term, i.e., non-Darcy flow, needs to be introduced. According to the Forchheimer equation, the calculation formula of the non-Darcy velocity is: (15) in β is the inertial drag coefficient, unit: 1 / m, c F is the Forchheimer dimensionless parameter.

[0014] When the gas in the rock matrix, cracks, and tunnels flows at a high Reynolds number, the Darcy velocity mentioned above is replaced by the non-Darcy velocity.

[0015] Preferably, the density of the gas ρ g The calculation is as follows: 1) Under ideal conditions, the ideal state equation is used to calculate the gas density. The equation is as follows: (16) Converted to density form: (17) in n is the number of moles of gas, unit: mol, V is the gas volume, unit: m 3 , P Represents the pressure of the gas, unit: Pa, M Represents the molar mass of the gas, kg / m 3 , that is, the mass of each mole of gas, R is the ideal gas constant of gas, dimensionless, for different gases and unit systems, R The values ​​of are different.T Represents the absolute temperature of the gas, usually in Kelvin K, introducing the compression factor Z The ideal gas state equation is corrected and the corrected equation is as follows: (18) in Z is the compressibility factor, which is used to modify the ideal gas state equation to reflect the difference between real gas and ideal gas; 2) Under non-ideal conditions, considering the adsorption characteristics of gas in coal, the Langmuir equation, a warm adsorption equation, is used for calculation. The equation is as follows: (19) in V g is the adsorption volume, V L , P L They respectively represent the maximum gas quantity and pressure that can be adsorbed by unit mass or unit volume of coal under saturated adsorption state; The van der Waals equation is a semi-empirical equation modified based on the ideal gas state equation: (20) in a , b is a constant determined by experiment.

[0016] There are differences in the calculation of gas density under ideal and non-ideal conditions. Therefore, the gas density in both cases should be considered and selected according to the situation during calculation.

[0017] Preferably, the step (2) further comprises assigning different material property numbers to fracture areas and different types of porous medium formations.

[0018] Unique material property numbers are assigned to fracture areas and different types of porous media formations to facilitate the subsequent implementation of differentiated seepage calculation strategies and ensure the high accuracy and applicability of the simulation process.

[0019] Beneficial effects of the present invention: This method reveals the complex interaction mechanism between matrix, cracks and tunnels, which helps to more accurately predict the behavior and distribution law of gas diffusion under different conditions, such as rock physical properties, crack distribution, gas pressure and temperature; and through high-precision mathematical models and algorithms, the diffusion process of gas in complex geological structures is simulated in detail; at the same time, compared with traditional experimental methods, numerical simulation has higher efficiency and lower cost, and can quickly obtain simulation results in a short time, which provides a strong scientific basis and decision-making support for the prediction and prevention of gas disasters, the optimization of coal mining design, the reduction of production costs and the utilization of coalbed methane; in addition, the simulation method can also provide intuitive three-dimensional visualization results, enabling researchers to clearly observe the temporal and spatial variation laws and characteristic analysis of gas in rock matrix, cracks and tunnels, thereby building a complete intelligent mining system. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] Figure 1 This is the flow chart of numerical simulation of gas diffusion in rock matrix-crack-tunnel; Figure 2 It is the coupled process of gas diffusion in rock matrix, cracks and tunnels; Figure 3 For the stratigraphic model; Figure 4 is the velocity streamline diagram of the yz section of the tunnel; Figure 5 It is the concentration contour line of the yz section of the tunnel. DETAILED DESCRIPTION

[0021] Example 1 S1: Extract stratum and tunnel information from the GIS system, thereby generating a three-dimensional porous matrix, two-dimensional plane fractures, and a simplified one-dimensional tunnel grid information, and determine the seepage field control equations of the matrix, fractures, and tunnels respectively; S2: Set the physical parameters of matrix, gas and fracture according to the specific mining area, including the location and size information of strata, fractures and tunnels, the ideal gas constant R of gas, the molar mass of gas, the absolute temperature of gas, the fluid dynamic viscosity coefficient, and the porosity of matrix; S3: Based on the porosity and rock pore size, the permeability of the porous matrix is ​​determined according to Kozeny-Carman, the fracture thickness is obtained and its permeability is determined according to the cubic law; S4: Based on the three-dimensional formation model, two-dimensional plane fractures, and simplified one-dimensional tunnels, as well as the corresponding physical parameters and boundary conditions, the seepage equations of the matrix and fractures and the convection-diffusion equations of the tunnels are discretely solved, and the coupling conditions are determined, that is, the pressure is the same and the flow rate satisfies the principle of superposition for coupling calculation to obtain the flow velocity information in different dimensions; S5: Analyze the gas concentration in the tunnel and set disaster warnings in the formation based on the concentration threshold and perform visual display.

[0022] Example 2 S1: Extract stratum and tunnel information from GIS system, scale and simplify to three stratum layers and connect tunnels. Figure 3 As shown in the figure, the calculation domain size is about 100m*40m*60m, the size of each stratum is about 100m*40m*20m, and the ideal gas state equation is used to simulate the physical properties of gas, that is, the molar mass of gas is 0.016kg / m 3 , the dynamic viscosity coefficient is 1e -3 Pa·s; S2: According to the actual working conditions, the boundary velocity above the formation is determined to be 1e -4 The m / s direction diffuses downward, and the tunnel is connected to the outside world, that is, the boundary pressure on the left, right, and sides is set to 1 [atm]; S3: According to the GIS system, different strata information is obtained and numbered, that is, 1, 2, and 3 from top to bottom, which are coarse-grained rock, coal seam, and fine-grained rock, respectively. The porosity and permeability parameters are set according to the physical properties of the matrix, as shown in Table 1. The data in Table 1 are substituted into the three-dimensional seepage equation of the matrix, and the volume source term q of the matrix is m The inflow velocity u through the fracture tf and the crack thickness d f Calculate, that is , to carry out numerical simulation of gas migration in the matrix; Table 1: Formation permeability parameters Strata 1 Strata 2 Strata 3 Porosity 0.32 0.4 0.34 Particle size [mm] 0.08 0.05 0.06 <![CDATA[Permeability [m 2 > <![CDATA[2.5196e -6 ]]> <![CDATA[1.2346e -7 ]]> <![CDATA[1.0828e -7 ]]> S4: The cross-sectional area of ​​the crack is about 0.01m 2 The fracture thickness is 0.02m, the fracture porosity is determined to be 0.8, and the fracture permeability is calculated according to the cubic law to be 3.33e -5 m 2 ; The physical parameters of the fracture are brought into the two-dimensional fracture seepage equation to perform numerical simulation of gas migration in the fracture, that is, substitute:

[0023] S5: Determine that the tunnel length is about 100m, the cross-sectional dimensions are 6m*6m, and the diffusion coefficient in the tunnel is 1.0e - 9 m 2 / s, the default gas concentration of the matrix diffusing into the tunnel is 1 mol / L, and the tunnel parameters are brought into the one-dimensional gas convection diffusion equation to perform numerical simulation of gas diffusion in the tunnel, that is, substitute:

[0024] S6: Calculate the volume source term of the matrix and the crack at different times, as well as the cross-sectional flow velocity between the crack and the roadway, and between the matrix and the roadway. 3 After s, the average flow velocity u that diffuses to the upper, lower, left, right, and rear surfaces of the tunnel and points to the interior of the tunnel is calculated. t They are: 238553e -4 m / s, 7.24681e -5 m / s, 1.53664e -4 m / s, 157997e - 5 m / s.

[0025] S7: Select the stratum section in the yz direction of the tunnel and display its concentration contours and total Darcy velocity streamlines. The results are as follows: Figure 4 and Figure 5 As shown, the gas diffusion inside the tunnel can be displayed through visualization.

Claims

1. A numerical simulation method for rock matrix-crack-tunnel gas diffusion, characterized by: The following steps are involved: (1) Based on the geological model and tunnel geometry provided by the GIS platform, the boundary range of the calculation area is defined, and the precise location of the fractures and tunnels, the width parameters of the fractures, and the cross-sectional dimensions of the tunnels at each point along their axis are determined; at the same time, the permeability properties of the rock matrix are determined in combination with its physical properties and spatial distribution. In addition, the boundary conditions and initial conditions required for the simulation are accurately obtained with the help of geological survey data and field monitoring data; (2) The geological model in (1) is converted into a high-precision grid model covering the entire calculation area; and low-dimensional cells are embedded inside the stratum to accurately simulate the geometric characteristics of the fracture surface; (3) Determine the gas control equations for the three systems of rock matrix, fractures, and tunnels; (4) Use the discrete fracture model to numerically simulate the fracture system, and use Darcy or non-Darcy seepage and skin coefficient to numerically simulate the fracture and matrix model according to the needs; (5) The gas diffusion in the tunnel is simplified into a one-dimensional convection-diffusion model with a rectangular cross section for numerical simulation; the state equations of the gas in the matrix-fracture-tunnel system are coupled to obtain the gas concentration in the tunnel; (6) Compare the gas concentration at different times in the tunnel obtained by numerical simulation with the actual gas threshold to conduct disaster warning in the formation, and visualize the simulation results of the matrix-fracture-tunnel.

2. The numerical simulation method for rock matrix-crack-tunnel gas diffusion according to claim 1, characterized in that: The three-dimensional seepage model of the rock matrix described in step (3) is carried out under the default gas saturation state, and its seepage field control equation is as follows: （1） in ε p is the porosity of the matrix, dimensionless, u g Represents the Darcy velocity of the fluid in the porous medium, unit: m / s, ρ g Indicates the density of gas, unit: kg / m 3 , t represents time, unit is s, q m is the volume source term of the matrix, unit: 1 / s, ▽ is the gradient; The governing equation of the fracture two-dimensional seepage model described in step (3) is as follows: （2） in ε f is the porosity of the fracture, dimensionless, u f Represents the Darcy velocity of the fluid in the fracture, unit: m / s, q f is the volume source term of the crack, unit: 1 / s; The governing equation of the one-dimensional seepage model of the roadway described in step (3) is as follows: （3） Where C is the gas concentration, unit: mol / m³, K is the diffusion coefficient, dimensionless, u t It is the Darcy velocity / average velocity of gas in the tunnel, unit: m / s.

3. The numerical simulation method for rock matrix-crack-tunnel gas diffusion according to claim 2 is characterized in that: The coupling of the three systems in step (5) includes: coupling between matrix and cracks, coupling between matrix and tunnels, and coupling between cracks and tunnels. The specific coupling process is as follows: 1) Coupling process between matrix and cracks: The volume source term of the two-dimensional fracture surface is calculated by the flow through the fracture surface, that is, by q m and q f To couple the matrix with the cracks, q f The calculation formula is as follows: （4） in n 1 is the crack normal direction, k f is the permeability of the fracture, unit: m 2 , u fp is the flow rate of gas from the matrix into the cracks, unit: m / s, A f is the crack surface area, unit: m 2 , V f is the volume of the crack, unit: m 3 , P Represents the pressure of the crack gas, unit: Pa, μ Represents the fluid dynamic viscosity coefficient, unit: Pa·s. For two seepage fields, the two are source / sink terms for each other, so the volume flow rates of the two are inverse to each other; therefore, the source term should be satisfied q f=- q m , to satisfy the coupling between matrix and cracks; 2) Coupling process between matrix and lane: According to the average Darcy velocity of the tunnel wall seepage field, the boundary conditions of the tunnel gas diffusion velocity are established, and the flow velocity entering the tunnel is u tm The calculation formula is: （5） in n t is the inner normal unit vector of the tunnel wall, u p is the Darcy velocity of the matrix at the tunnel surface, from the matrix to the inside of the tunnel; therefore, the average flow velocity entering the tunnel is u t1 for: （6） In the formula is the wall area per unit length of the tunnel; 3) Coupling process between cracks and tunnels: If the crack happens to be connected to the tunnel: determine the cross-sectional area of ​​the crack by the crack width and length of the tunnel wall, and calculate the cross-sectional area of ​​the crack according to the flow rate of gas in the crack. u tf As the velocity boundary condition for tunnel gas diffusion, the calculation formula is: （7） In the formula u f represents the Darcy velocity of the fluid in the fracture, and the average flow velocity entering the tunnel is u t2 for: （8） in is the crack length per unit length of tunnel, unit: .

4. The numerical simulation method for rock matrix-crack-tunnel gas diffusion according to claim 3 is characterized by: The calculation formula of Darcy velocity u of fluid in porous media is: （9） in k is the permeability, unit: m 2 ,use k p Expressed as the permeability of the rock matrix, unit: m 2 , k f Expressed as the permeability of the fracture, unit: m 2 , μ Represents the fluid dynamic viscosity coefficient, unit: Pa·s, p is the gas pressure, unit: Pa, ▽p is the pressure gradient, unit: Pa / m; in, u g Represents the Darcy velocity of gas in porous media, unit: m / s, u f Represents the Darcy velocity at the crack surface, unit: m / s, u t is the Darcy velocity of gas in the tunnel, unit: m / s, all are applicable to this calculation; Permeability of rock matrix k p Derived from the Kozeny-Carman relationship: （10） in, C 1 is the Kozeny-Carman constant, s For specific surface, in a packed bed or regular structure of granular soil, the rock matrix permeability is derived from the Kozeny-Carman formula: （11） in, d p is the effective particle size of the particles that make up the matrix, unit: m; The permeability k of the fracture f It is expressed as: （12） in d f is the crack thickness, unit: m; Will u t2 as well as u t1 ， Substitute them into formula (3) to obtain the gas concentration in the tunnel under the two conditions.

5. The numerical simulation method for rock matrix-crack-tunnel gas diffusion according to claim 3 is characterized by: Skin Factor S k Revised permeability, skin factor S k The relationship with permeability is: （13） where k, k d are the formation permeabilities before and after damage, unit: m 2 ; r e Indicates the radius of the contaminated area of ​​the nearby strata, unit: m, relative to the center of the fracture or tunnel; r w Indicates the area radius of the tunnel or crack, unit: m; Where k is the formation permeability before damage, the following formula is obtained: （14） in β is the inertial drag coefficient, unit: 1 / m, c F is the Forchheimer dimensionless parameter.

6. The numerical simulation method for rock matrix-crack-tunnel gas diffusion according to claim 5, characterized in that: In high Reynolds number flow, nonlinear terms, namely non-Darcy flow, need to be introduced. According to the Forchheimer equation, the calculation formula of its non-Darcy velocity is: （15） in β is the inertial drag coefficient, unit: 1 / m, c F is the Forchheimer dimensionless parameter.

7. The numerical simulation method for rock matrix-crack-tunnel gas diffusion according to claim 2 is characterized by: Density of gas ρ g The calculation is as follows: 1) Under ideal conditions, the ideal state equation is used to calculate the gas density. The equation is as follows: （16） Converted to density form: （17） Where n is the number of moles of gas, unit: mol, V is the volume of gas, unit: m 3 , P represents the gas pressure, unit: Pa, M represents the molar mass of the gas, kg / m 3 , that is, the mass of each mole of gas, R is the ideal gas constant of gas, dimensionless, and the value of R is different for different gases and unit systems. T Represents the absolute temperature of the gas, usually in Kelvin K. The compression factor Z is introduced to correct the ideal gas state equation. The corrected equation is as follows: （18） in Z is the compressibility factor, which is used to modify the ideal gas state equation to reflect the difference between real gas and ideal gas; 2) Under non-ideal conditions, considering the adsorption characteristics of gas in coal, the Langmuir temperature adsorption equation is used for calculation, and the equation is as follows: （19） in V g is the adsorption volume, V L , P L They respectively represent the maximum gas quantity and pressure that can be adsorbed by unit mass or unit volume of coal under saturated adsorption state; The van der Waals equation is a semi-empirical equation modified based on the ideal gas state equation: （20） in a , b is a constant determined by experiment.

8. The numerical simulation method for rock matrix-crack-tunnel gas diffusion according to claim 1 is characterized in that : The step (2) also includes assigning different material attribute numbers to fracture areas and different types of porous medium formations.

Images