A numerical simulation method for the coupled water migration in matrix, fractures and roadways

Abstract

The present invention relates to the technical field of intelligent mine exploitation, and particularly to a numerical simulation method for the coupled migration of matrix-fracture-tunnel water. Specifically, it includes the following steps: extracting matrix and tunnel information from the GIS system to generate three-dimensional matrix, two-dimensional fracture, and simplified one-dimensional tunnel grid information; setting physical property parameters, pressure, and velocity boundary conditions according to the actual measured mine conditions; obtaining the permeability of the three-dimensional matrix and two-dimensional fracture; discretely solving the seepage equations of the three-dimensional matrix and two-dimensional fracture and the Bernoulli equation of the tunnel, and setting coupling conditions, namely the same pressure and the principle that the flow rates satisfy superposition for coupling calculation to obtain the flow velocity information in different dimensions; obtaining the water inrush flow rate and water inrush water level after rapid calculation of the one-dimensional tunnel, and on this basis, using the SPH method for local visualization display and early warning. Through this method, the migration situation of groundwater can be completely simulated, thereby achieving the purpose of optimizing engineering design.

Description

Technical Field

[0001] The present invention relates to the technical field of intelligent mine exploitation, and particularly to a numerical simulation method for the coupled water migration in matrix-fracture-roadway. Background Art

[0002] The water migration during the mine exploitation process is a complex process, involving multiple links such as the recharge, flow, influx into the mine, drainage, and discharge of groundwater. It is mainly affected by factors such as geological structures, rock permeability, groundwater level, and mining activities. Among them, geological structures such as faults, joints, and fractures in the geological structure affect the flow path and velocity of groundwater; the permeability of rocks determines the flow velocity and path of groundwater; rocks with higher permeability are prone to cause rapid influx of groundwater into the mine, and the level of the groundwater table directly affects the amount and location of mine water inrush. A high groundwater table will increase the risk of mine water inrush; mining activities change the natural flow path of groundwater, which may lead to a decline in the groundwater level and mine water inrush. Therefore, the analysis and simulation methods of water migration have always been one of the key points in the construction of intelligent mines.

[0003] Due to the complexity of the coupling effect and the high heterogeneity of fractured rock masses, using a mathematical model for numerical simulation is an important means to analyze the seepage and mass transfer process of fractured rock masses. Currently, since the water migration process involves various different regional working conditions, forming a multi-region coupling system, it becomes more complex and changeable, and it is difficult to accurately describe it with a single model. And the simulation for this kind of situation is often just a simple numerical simulation scheme in the seepage field, ignoring the complete water migration situation in the strata mining, resulting in the inability to accurately and comprehensively judge the water disaster risk during the mining process. Summary of the Invention

[0004] The purpose of the present invention is to overcome the deficiencies in the prior art and provide a numerical simulation method for the coupled water migration in matrix-fracture-roadway. This method can analyze the water migration situation based on the pressure and flow rate information in each region according to the actual mining situation, can completely simulate the water migration situation of groundwater, accurately predict the seepage field, reveal the fluid-solid interaction mechanism, and thus achieve the purpose of evaluating environmental risks and optimizing engineering designs.

[0005] To solve the problems in the prior art, the present invention discloses a numerical simulation method for the coupled water migration in matrix-fracture-roadway, which includes the following steps:

[0006] Step 1: Establish the geometric model of the measured mine: According to the actual strata and roadway geometric model of the measured mine on the GIS system (Geographic Information System), determine the size of the calculation domain, as well as the positions of fractures and roadways, the widths of fractures, and the cross-sectional dimensions at various points along the axis of the roadway. Then, determine its permeability based on the physical properties and distribution of the rock matrix of the measured mine, and obtain the boundary conditions and initial conditions according to the survey and on-site information to determine the geometric model of the measured mine.

[0007] Step 2: Generate the grid model of the entire calculation domain: Import the geometric model determined in the first step into the grid algorithm to generate the grid model of the entire calculation domain. Insert two-dimensional elements into the stratum information of this grid model to simulate the fracture surface. As a weak medium in the rock mass, the fracture network has the material properties of high permeability and low strength and must be fully considered in the simulation. That is, a two-dimensional element without thickness is used for calculation. This method has significant advantages in simulating thin-layer structures such as fractures, cracks, and interface transition zones.

[0008] Step 3: Calculate the seepage field and permeability of the rock matrix: Use different three-dimensional seepage models according to the different matrices of the measured mine to calculate its seepage field and permeability.

[0009] The control equation of the seepage field is expressed as formula (1): Where S p represents the storage coefficient of the porous medium, which is calculated through the porosity of the porous medium and the fluid density, p represents pressure, u 达西 represents the Darcy velocity of the fluid in the porous medium, t represents time, q m is the volume source term of the matrix.

[0010] When the matrix of the measured mine is rock, the permeability k of the matrix p is expressed as formula (2): Where ε is the porosity of the matrix, C is the KC constant, and s is the specific surface area.

[0011] When the matrix of the measured mine is a regular structure such as a packed bed or granular soil, the permeability k of the matrix p is expressed as formula (3): Where d p is the effective particle size of the particles constituting the matrix, and ε is the porosity of the matrix.

[0012] Step 4: Calculate the seepage field and permeability of the fractures: Use the discrete fracture model to numerically simulate the fractures. The seepage field equation in the fractures is expressed as formula (4): Where d f is the fracture thickness, S f represents the storage coefficient of the fracture, t is time, uf Denote the fluid flow velocity in the fracture as q f which is the volume source term of the fracture.

[0013] The permeability k of the fracture f is expressed by Equation (5): where d f is the fracture thickness.

[0014] Step 5: Conduct numerical simulations on the matrix and fracture models according to two models of fluid flow through porous media:

[0015] u 达西 The calculation formula is expressed by Equation (6): where k is the permeability of the porous medium or fracture, μ represents the dynamic viscosity of the fluid, and p represents the pressure;

[0016] When the fluid with high-speed flow passes through the complex porous matrix structure, it is non-Darcy flow. According to the Forchheimer equation, its non-Darcy velocity u 非达西 is expressed by Equation (7): μ represents the dynamic viscosity of the fluid. Among them, β is the inertial resistance coefficient, c F is the Forchheimer dimensionless parameter, k is the permeability of the porous medium or fracture, p represents the pressure, and ρ is the fluid density;

[0017] According to the actual situation of the measured mine, substitute each parameter into Equation (6) or Equation (7) to calculate its Darcy velocity u 达西 or non-Darcy velocity u 非达西 .

[0018] Step 6: Conduct numerical simulations on the roadway using the one-dimensional open-channel flow model: Simplify the roadway section into a rectangle, then the water transport in the roadway is simplified into a one-dimensional open-channel flow model with a rectangular cross-section, and the Saint-Venant equations are used for simulation.

[0019] Its continuity equation is expressed by Equation (8):

[0020] Its momentum equation is expressed by Equation (9): where A is the cross-sectional area, determined by the roadway cross-sectional shape and the submerged depth, t is the time, Q is the flow rate through the cross-section, x is the position of the roadway, u is the average cross-sectional velocity, g is the acceleration due to gravity, Z is the effective cross-section, that is, the relative height of the centroid of the cross-section through which the flow passes and the reference horizontal plane, and f is the sum of the frictional head losses along the way and local resistances.

[0021] For a rectangular roadway, the continuity equation is simplified to Equation (10):

[0022] The momentum equation is simplified to Equation (11): where w is the roadway width, h is the water depth, t is the time, Q is the flow rate through the cross-section, u is the average cross-sectional velocity, and f is the sum of the frictional head losses along the length and local losses.

[0023] Step 7: Coupled process of matrix-fracture-roadway water migration:

[0024] Boundary conditions and source terms are established based on the principle of equal pressure and flow rate superposition, specifically divided into three parts: the coupling of matrix-fracture, matrix-roadway, and fracture-roadway:

[0025] A: Coupling process of matrix and fracture: The mass source term of the two-dimensional fracture surface is calculated through the flow rate across the fracture surface. For the two seepage fields of matrix and fracture, they are each other's source / sink terms. Therefore, their volume flow rates are opposite to each other, that is, the source term q f =-q m is satisfied, and then the volume source term q f of the fracture and the flow velocity u fp of the fluid flowing from the matrix into the fracture are obtained, as shown in Equation (12): where n is the fracture normal, k f is the permeability of the fracture, u fp is the flow velocity of the matrix flowing into the fracture, d f is the fracture thickness, A f is the fracture surface area, V f is the fracture volume, μ is the dynamic viscosity of the fluid, and p is the pressure.

[0026] B. Coupling process of matrix and roadway:

[0027] Based on the average Darcy velocity of the seepage field on the roadway wall, the boundary condition of the roadway water inrush velocity is established. The water flow u tm entering the roadway is expressed as Equation (13): u tm =u 达西 ·n t (13), where n t is the inner normal unit vector of the roadway wall, pointing from the matrix to the inside of the roadway, and u 达西 represents the Darcy velocity of the fluid in the matrix in the porous medium.

[0028] Then the volume flow rate Q tm entering the roadway is expressed as Equation (14): where A tm is the wall area per unit length of the roadway, and u tm is the water flow entering the roadway obtained in Equation (13).

[0029] C. Coupling process of fracture and roadway:

[0030] Determine the cross-sectional area of the fissure based on the width and length of the fissures on the roadway wall surface, and use the fissure flow velocity as the velocity u of water inrush from the roadway tc The boundary condition is expressed by Equation (15): u tc = u f ·n t (15), where u f represents the flow velocity of the fluid in the fissure, and n t is the unit normal vector of the inner surface of the roadway wall, pointing from the matrix to the inside of the roadway.

[0031] Then the volume flow rate Q entering the roadway tm is expressed by Equation (16), where u tf is the fluid velocity entering the roadway from the fissure, s tf is the fissure length per unit length of the roadway, f is the sum of the head losses due to friction along the length and local resistances, and s is the specific surface area.

[0032] Eighth step: Visualization display and setting of threshold warning: Use a three-dimensional visualization tool to display the flow velocity, cross-sectional shape, and dimensions in the roadway in step (S7) to show the water submergence depth, and visualize the simulation results of the matrix-fissure-roadway; and conduct local refined display of water inrush from the roadway based on the SPH method. According to the display situation, analyze the water inrush water level and set the water inrush water level threshold, and set disaster warnings in the formation based on the water inrush water level threshold.

[0033] Preferably, in the third step, according to the different materials of the measured mine matrix, different three-dimensional seepage models are adopted to calculate its seepage field and permeability.

[0034] Preferably, for the effective cross-section Z in Equation (9) and Equation (11) in the sixth step, its calculation formula is expressed as Equation (10): Z = Z b + h / 2, where Z b is the elevation of the roadway bottom, and h is the water depth.

[0035] Preferably, in the sixth step, f in Equation (9) and Equation (11) is obtained by calculating the water flow velocity using the Chezy-Manning formula and then substituting the water flow velocity into the Darcy-Weisbach formula.

[0036] Preferably, in the second step, the two-dimensional element grid for simulating the fissure surface is calculated using a two-dimensional element without thickness.

[0037] Preferably, in the seventh step, the coupling surfaces at the same pressure position are selected for coupling of the coupling surfaces from the three-dimensional matrix to the two-dimensional fissure, from the two-dimensional fissure to the one-dimensional roadway, and from the three-dimensional matrix to the one-dimensional roadway. On the coupling surface, the sum of the flow rates in different dimensions is equal to the total flow rate.

[0038] Preferably, the 3D visualization tool in the eighth step is one of ParaView, VisIt, MATLAB, and Blender.

[0039] Preferably, in the coupling process of the matrix and the roadway in step (S7), the fluid flow of the matrix and the roadway is Darcy flow, so the water flow u entering the roadway tm is expressed by formula (13): u tm = u 达西 ·n t ((13), where n t is the unit normal vector of the inner surface of the roadway wall, pointing from the matrix to the inside of the roadway, and u 达西 represents the Darcy velocity of the fluid in the matrix in the porous medium, and u 达西 is calculated by formula (6) in step (S5).

[0040] More preferably, when the fluid flow in the application scenario is non-Darcy flow, the u in formula (13) 达西 is replaced by the non-Darcy velocity of the fluid in the matrix in the porous medium and is calculated by formula (7) in step (S5).

[0041] In the coupling across dimensions, it generally includes steps such as defining the coupling surface, setting boundary conditions, setting source terms, solving the coupling model, and verifying and optimizing. In this application, we define the coupling surface between the three-dimensional matrix, two-dimensional fracture, and one-dimensional roadway, and obtain the boundary conditions and initial conditions according to the survey and on-site information, that is, the content shown in the first and second steps; while the source terms of the matrix and fracture are set in the third, fourth, and fifth steps respectively; the source term of the roadway is set in the sixth step; in the seventh step, pairwise coupling is performed for these three dimensions, where the flow transfer from the three-dimensional matrix to the two-dimensional fracture is through the mass source term, the flow transfer between the two-dimensional fracture and the one-dimensional roadway is through the average Darcy velocity at the fracture-roadway cross-section, and the flow transfer between the three-dimensional matrix and the one-dimensional roadway is through the average Darcy velocity at the matrix-roadway cross-section. Specific coupling calculations can use numerical methods, such as the finite difference method, finite element method, finite volume method, etc., to solve the coupling model until a definite and optimal solution method is obtained. Then, through the deduction of this method, the water transport situation of the matrix-fracture-roadway is displayed in the 3D visualization tool, and then we can analyze the water inrush level and set the disaster warning in the formation according to the water inrush level threshold, ultimately achieving the purpose of intelligent mining and early warning.

[0042] The beneficial effects of the present invention are as follows:

[0043] 1. It reveals the complex interaction mechanism between the matrix, fracture, and roadway, which helps to more accurately predict the behavior and distribution law of water transport under different conditions (such as different hydraulic gradients, fracture widths, matrix porosities, etc.).

[0044] 2. It can effectively avoid the assessment of environmental risks and provide important basis for engineering design in the fields of groundwater engineering, mine drainage, tunnel engineering, etc.

[0045] 3. By simulating the water migration under different design schemes, the engineering layout can be optimized, water resource waste can be reduced, and engineering benefits can be improved.

[0046] 4. Through intuitive graphical display, comprehensive visualization can be achieved, and the spatio-temporal variation law and characteristic analysis of water migration can be clearly observed, so as to better construct an intelligent mine system. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1 is the flowchart of the present invention;

[0048] Figure 2 is the flowchart of the matrix-fracture-roadway water migration coupling scheme in the present invention;

[0049] Figure 3 is the total Darcy velocity streamline of the topographic section in the present invention;

[0050] Figure 4 is the pressure contour line of the topographic section in the present invention;

[0051] Figure 5 is the local refined display of roadway water inrush in the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0052] The present invention will be further described below with reference to the drawings. The following embodiments are only used to more clearly illustrate the structure of the present invention.

[0053] As Figure 1-2 shown, a numerical simulation method for matrix-fracture-roadway water migration coupling includes the following steps:

[0054] The first step: Establish a geometric model of the measured mine: According to the actual strata and roadway geometric model of the measured mine on the GIS system (Geographic Information System), determine the size of the calculation domain, as well as the positions of fractures and roadways, the width of fractures, and the cross-sectional dimensions of the roadways at various positions along the axis direction. Then, determine its permeability according to the physical properties and distribution of the rock matrix of the measured mine, and obtain the boundary conditions and initial conditions based on the survey and on-site information to determine the geometric model of the measured mine.

[0055] Step 2: Generate the grid model of the entire computational domain: Import the geometric model determined in Step 1 into the grid algorithm to generate the grid model of the entire computational domain. Insert two-dimensional elements into the formation information of this grid model to simulate the fracture surface. As a weak medium in the rock mass, the fracture network has the material properties of high permeability and low strength and must be fully considered in the simulation. That is, a two-dimensional element without thickness is used for calculation. This method has significant advantages in simulating thin-layer structures such as fractures, cracks, and interface transition zones.

[0056] Step 3: Calculate the seepage field and permeability of the rock matrix: Use different three-dimensional seepage models to calculate the seepage field and permeability according to the different measured mine matrices.

[0057] The control equation of the seepage field is expressed as Equation (1): where S p represents the storage coefficient of the porous medium, which is calculated through the porosity of the porous medium and the fluid density, p represents the pressure, u 达西 represents the Darcy velocity of the fluid in the porous medium, t represents time, q m is the volume source term of the matrix.

[0058] When the matrix of the measured mine is rock, the permeability k of the matrix p is expressed as Equation (2): where ε is the porosity of the matrix, C is the KC constant, and s is the specific surface area.

[0059] When the matrix of the measured mine is a regular structure such as a packed bed or granular soil, the permeability k of the matrix p is expressed as Equation (3): where d p is the effective particle size of the particles constituting the matrix, and ε is the porosity of the matrix.

[0060] In actual operation, different models can be established according to the different matrix materials and different fluids, and then different models can be represented by different labels. When used later, they can be directly selected correspondingly.

[0061] Step 4: Calculate the seepage field and permeability of the fractures: Use the discrete fracture model to numerically simulate the fractures. The seepage field equation in the fractures is expressed as Equation (4): where d f is the fracture thickness, S f represents the storage coefficient of the fracture, t is time, u f represents the flow velocity of the fluid in the fracture, q f is the volume source term of the fracture.

[0062] The permeability k of the fracture f is expressed as Equation (5): where d f is the fracture thickness.

[0063] Step 5: Numerically simulate the matrix and fracture models according to two models of fluid flow through porous media:

[0064] u 达西 The calculation formula is expressed as Formula (6): where k is the permeability of the porous medium or fracture, μ represents the dynamic viscosity coefficient of the fluid, and p represents the pressure;

[0065] When the fluid with high-speed flow passes through the complex porous matrix structure, it is non-Darcy flow. According to the Forchheimer equation, its non-Darcy velocity u 非达西 is expressed as Formula (7): μ represents the dynamic viscosity coefficient of the fluid. Among them, β is the inertial resistance coefficient, and c F is the Forchheimer dimensionless parameter, k is the permeability of the porous medium or fracture, p represents the pressure, and ρ is the fluid density;

[0066] According to the actual situation of the measured mine, substitute each parameter into Formula (6) or Formula (7) to calculate its Darcy velocity u 达西 or non-Darcy velocity u 非达西 .

[0067] Step 6: Numerically simulate the roadway using the one-dimensional open-channel flow model: Simplify the roadway section into a rectangle, then simplify the water transport at the roadway into a one-dimensional open-channel flow model with a rectangular cross-section, and use the Saint-Venant equations for simulation.

[0068] Its continuity equation is expressed as Formula (8):

[0069] Its momentum equation is expressed as Formula (9): where A is the cross-sectional area, determined by the shape of the roadway cross-section and the submergence depth, t is the time, Q is the flow rate through the cross-section, x is the position of the roadway, u is the average cross-sectional velocity, g is the acceleration due to gravity, Z is the effective cross-section, that is, the relative height of the centroid of the cross-section through which the flow passes and the reference horizontal plane, and f is the sum of the frictional head losses along the way and local resistances.

[0070] For a rectangular roadway, the continuity equation is simplified to Formula (10):

[0071] And the momentum equation is simplified to Formula (11): where w is the roadway width, h is the water depth, t is the time, Q is the flow rate through the cross-section, u is the average cross-sectional velocity, and f is the sum of the frictional head losses along the way and local resistances.

[0072] Step 7: Matrix-fracture-tunnel water migration coupling process:

[0073] Boundary conditions and source terms are established using the principle of equal pressure and superposed flow rates, specifically divided into three parts: matrix-fracture coupling, matrix-tunnel coupling, and fracture-tunnel coupling:

[0074] A: Matrix-fracture coupling process: The mass source term of the two-dimensional fracture surface is calculated through the flow rate passing through the fracture surface. For the two seepage fields of the matrix and the fracture, they are each other's source / sink terms. Therefore, their volume flow rates are opposite to each other, that is, the source term q f =-q m is satisfied, and then the volume source term q f of the fracture and the flow velocity u fp of the fluid flowing from the matrix into the fracture are obtained, as shown in formula (12): where n is the fracture normal direction, k f is the permeability coefficient of the fracture, u fp is the flow velocity of the matrix flowing into the fracture, d f is the fracture thickness, A f is the fracture surface area, V f is the fracture volume, μ is the fluid dynamic viscosity coefficient, and p is the pressure.

[0075] B. Matrix-tunnel coupling process:

[0076] Based on the average Darcy velocity of the seepage field on the tunnel wall surface, the boundary condition of the tunnel water inrush velocity is established. The water flow into the tunnel is expressed as formula (13): u tm =u 达西 ·n t (13), where n t is the inner normal unit vector of the tunnel wall surface, pointing from the matrix to the inside of the tunnel, and u 达西 represents the Darcy velocity of the fluid in the matrix in the porous medium.

[0077] Then the volume flow rate Q tm flowing into the tunnel is expressed as formula (14): where A tm is the wall surface area per unit length of the tunnel.

[0078] C. Fracture-tunnel coupling process:

[0079] The cross-sectional area of the fracture is determined by the slit width and length of the tunnel wall surface, and the fracture flow velocity is used as the velocity boundary condition of the tunnel water inrush, expressed as formula (15): u tc =u f ·n t (15), where u f represents the flow velocity of the fluid in the fracture, and n t$\vec{n}$ is the unit normal vector of the roadway wall, pointing from the matrix to the inside of the roadway.

[0080] Then the volumetric flow rate $Q$ entering the roadway tm is expressed as Equation (16), where $u$ tf is the fluid velocity entering the roadway from the fracture, and $s$ tf is the fracture length per unit length of the roadway.

[0081] Eighth step: Visualization display and setting of threshold warning: Use a three-dimensional visualization tool to display the flow velocity, cross-sectional shape, and size in the roadway in step (S7) to show the water inundation depth, and visualize the simulation results of the matrix-fracture-roadway; and perform local refinement display of roadway water inrush based on the SPH method. According to the display situation, analyze the water inrush water level and set the water inrush water level threshold, and set the disaster warning in the formation according to the water inrush water level threshold.

[0082] Obviously, in the coupling process between the matrix and the roadway in step (S7), the fluid flow in the matrix and the roadway is Darcy flow. Therefore, the water flow $u$ entering the roadway tm is expressed as Equation (13): $u$ tm $= u$ 达西 $\cdot\vec{n}$ t ((13), where $\vec{n}$ t is the unit normal vector of the roadway wall, pointing from the matrix to the inside of the roadway, and $u$ 达西 represents the Darcy velocity of the fluid in the matrix in the porous medium, and $u$ 达西 is calculated by Equation (6) in step (S5).

[0083] When the fluid flow in the application scenario is non-Darcy flow, $u$ in Equation (13) 达西 is replaced by the non-Darcy velocity of the fluid in the matrix in the porous medium and is calculated by Equation (7) in step (S5).

[0084] Briefly speaking, we extract formation and roadway information from the GIS system, and thus generate information including three-dimensional porous matrix, two-dimensional planar fractures, and simplified one-dimensional roadway grids; then we set the physical properties of the matrix, fluid, and fractures according to the specific mining area conditions, specifically including the position and size information of the formation, fractures, and roadways, the density of the fluid, the dynamic viscosity coefficient, the porosity of the rock, etc., and set the pressure or velocity boundary conditions according to the influence of the overlying aquifer; then according to the porosity and rock pore size, we set the permeability of the porous matrix according to Kozeny-Carman, obtain the fracture thickness and set its permeability according to the cubic law; through the constructed three-dimensional formation model, two-dimensional planar fractures, and simplified one-dimensional roadway, as well as the corresponding physical property parameters and boundary conditions, we discretely solve the seepage equations of the matrix and fractures and the Bernoulli equation of the roadway respectively, and set the coupling condition that the pressure is the same and the flow rate satisfies the superposition principle for coupling calculation to obtain the flow velocity information in different dimensions; finally, we analyze the water inrush level and set the disaster warning in the formation and conduct visual display according to the water inrush level threshold.

[0085] Specifically, we extract the full-height mining coal seam and roadway information from the GIS system. After scaling and simplifying, the calculation domain size of the formation is about 1m * 6m * 0.5m, the roadway length is about 0.3m, and the cross-sectional size of the fracture is 3.5E -3 m * 4E -3 m.

[0086] The cross-sectional area of the fracture is about 0.01m 2 , and the fracture thickness is 0.02m.

[0087] The density of the fluid is 1000kg / m 3 , and the dynamic viscosity coefficient is 1E -3 (Pa·s).

[0088] The porosity of the formation matrix is set to 0.15 according to the physical properties of fine-grained rocks.

[0089] According to the porosity and rock pore size, the permeability k of the porous matrix is calculated according to formula (2) p to be 1.875E - 9 m 2 .

[0090] According to formula (5), the permeability k of the fracture f is 3.33E -5 m 2 .

[0091] Set the velocity of the upper boundary of the formation to 1E -7 m / s and the seepage direction is downward.

[0092] The pressure at the roadway position is 1 atm, which is one standard atmosphere.

[0093] After model deduction, at 1E 3 s, select the formation cross-section in the yz direction of the formation and display its total Darcy velocity streamline Figure 3 and pressure isobars Figure 4 .

[0094] The calculated normal total flow rate through the roadway surface and pointing into the roadway per unit time is 2.68518E - 9 m 3 / s.

[0095] After one-dimensionalizing the roadway, the calculated water depth of the roadway is 2.2E -3 m. Through visualization, the water inrush situation of the roadway can be displayed, as Figure 5 shown.

[0096] Then we can set the water inrush water level threshold to (2.2E -3 ~0.3) m. During actual mining, when the water inrush water level threshold is exceeded, a warning is issued to evacuate relevant personnel.

[0097] The above are only the preferred embodiments of the present invention. It should be noted that for those of ordinary skill in the art, without departing from the technical principle of the present invention, several improvements and deformations can be made, and these improvements and deformations should also be regarded as the protection scope of the present invention.

Claims

1. A numerical simulation method for matrix-fracture-tunnel water transport coupling, characterized by: The following steps are involved: (S1) Establishing a geometric model of the measured mine: According to the actual strata and tunnel geometric model of the measured mine on the GIS system, determine the size of the calculation domain, the location of the cracks and tunnels, the width of the cracks and the cross-sectional dimensions of the tunnels at various locations along the axis direction, and then determine its permeability according to the physical properties and distribution of the rock matrix of the measured mine, and obtain boundary conditions and initial conditions according to the survey and field information to determine the geometric model of the measured mine; (S2) generating a grid model of the entire computational domain: importing the geometric model determined in step (S1) into a grid algorithm to generate a grid model of the entire computational domain, and inserting two-dimensional cells into the stratigraphic information of the grid model to simulate the fracture surface; (S3) Calculating the seepage field and permeability of the mine matrix: Selecting a three-dimensional seepage model based on the measured mine matrix material to calculate its seepage field and permeability, The governing equation of the seepage field is expressed as formula (1): Where S p represents the water storage coefficient of porous media, which is calculated by the porosity of porous media and fluid density, p represents pressure, u 达西 represents the Darcy velocity of the fluid in the porous medium, t represents time, and q m The volume source term of the matrix, When the matrix of the measured mine is rock, the permeability k of the matrix is p It is expressed as formula (2): Where ε is the porosity of the matrix, C is the KC constant, and s is the specific surface area. When the matrix of the measured mine is a regular structure such as a packed bed or granular soil, the permeability rate k of the matrix is p It is expressed as formula (3): where d p is the effective particle size of the particles constituting the matrix, ε is the porosity of the matrix; (S4) Calculate the seepage field and permeability of the fracture: The fracture is numerically simulated using a discrete fracture model, and the seepage field equation in the fracture is expressed as formula (4): where d f is the crack thickness, S f represents the water storage coefficient of the crack, t is the time, u f represents the flow velocity of the fluid in the fracture, q f is the volume source term of the fracture, The permeability k of the fracture f It is expressed as formula (5): where d f is the crack thickness; (S5) Numerical simulations of matrix and fracture models are performed based on two models of fluid flow through porous media: When a low-speed fluid flows through a uniform porous matrix structure, it is a Darcy flow, and its Darcy velocity u 达西 The calculation formula is expressed as formula (6): Where k is the permeability of the porous medium or fracture, μ represents the fluid dynamic viscosity coefficient, and p represents the pressure; When a high-speed fluid flows through a complex porous matrix structure, it is a non-Darcy flow. According to the Forchheimer equation, its non-Darcy velocity u 非达西 It is expressed as formula (7): μ represents the fluid dynamic viscosity coefficient, β is the inertial resistance coefficient, c F is the Forchheimer dimensionless parameter, k is the permeability of the porous medium or fracture, p represents the pressure, and ρ is the fluid density; According to the actual situation of the measured mine, the parameters are substituted into formula (6) or formula (7) to calculate the Darcy speed u 达西 or non-Darcy speed u 非达西 ; (S6) Numerical simulation of the tunnel using a one-dimensional open channel flow model: The tunnel cross section is simplified to a rectangle, and the water movement in the tunnel is simplified to a one-dimensional open channel flow model with a rectangular cross section, and the Saint-Venant equations are used for simulation. Its continuity equation is expressed as formula (8): Its momentum equation is expressed as formula (9): Where A is the cross-sectional area, which is determined by the cross-sectional shape and flooding depth of the tunnel, t is the time, Q is the flow rate through the cross-section, x is the position of the tunnel, u is the average flow velocity of the cross-section, g is the gravitational acceleration, Z is the effective cross-sectional area, that is, the relative height between the centroid of the flow cross-sectional area and the reference horizontal plane, and f is the sum of the longitudinal and local resistances caused by friction; For rectangular tunnels, the continuity equation is simplified to formula (10): The momentum equation is simplified to formula (11): Where w is the width of the tunnel, h is the water depth, t is the time, Q is the flow rate through the section, u is the average flow velocity of the section, and f is the sum of the longitudinal and local resistance caused by friction; (S7) Matrix-fracture-tunnel water migration coupling process: The boundary conditions and source terms are established based on the principle of equal pressure and flow superposition. They are divided into three parts: matrix-crack coupling, matrix-roadway coupling, and crack-roadway coupling: A: Coupling process between matrix and fracture: The mass source term of the two-dimensional fracture surface is calculated by the flow through the fracture surface. For the two seepage fields of matrix and fracture, the two are source / sink terms for each other, so the volume flow rates of the two are opposite to each other, that is, the source term q is satisfied. f =-q m , then we get the volume source term q of the crack f , the flow rate u of the fluid flowing from the matrix into the fracture fp， As shown in formula (12), Where n is the crack normal, k f is the permeability coefficient of the fracture, u fp is the flow rate of matrix into the fracture, d f is the crack thickness, A f is the crack surface area, V f is the volume of the crack, μ is the fluid dynamic viscosity coefficient, and p is the pressure; B. Coupling process between matrix and lane: According to the average Darcy velocity of the tunnel wall seepage field, the tunnel water inrush velocity boundary condition is established, and the water flow u entering the tunnel is tm Expressed as formula (13)u tm =u 达西 ·n t (13), where n t is the inner normal unit vector of the tunnel wall, pointing from the matrix to the inside of the tunnel, u 达西 represents the Darcy velocity of the fluid in the matrix in the porous medium, The volume flow rate entering the tunnel is Q tm It is expressed as formula (14): Among them A tm is the wall area per unit length of the tunnel, A is the cross-sectional area, u tm is the water flow entering the tunnel obtained in formula (13); C. The coupling process between cracks and tunnels: The crack cross-sectional area is determined by the crack width and crack length of the tunnel wall, and the crack flow velocity is used as the speed of water inrush in the tunnel u. tc The boundary condition is expressed as formula (15): tc =u f ·n t (15), where u f represents the flow velocity of the fluid in the fracture, n t is the inner normal unit vector of the tunnel wall, pointing from the matrix to the inside of the tunnel, The volume flow rate entering the tunnel is Q tm It is expressed as formula (16), Among them, u tf is the velocity of the fluid entering the tunnel from the fracture, s tf is the crack length per unit length of the tunnel, f is the sum of the longitudinal and local resistance caused by friction, and s is the specific surface; (S8) Visualization and setting of threshold warning: Using a three-dimensional visualization tool, the flow velocity, cross-sectional shape, and size in the tunnel of step (S7) are used to display the water flooding depth, and the simulation results of the matrix-fracture-tunnel are visualized. Based on the SPH method, a local refined display of tunnel water inrush is performed. According to the display situation, the water inrush level is analyzed and a water inrush level threshold is set. According to the water inrush level threshold, a disaster warning in the formation is set.

2. The method for numerical simulation of matrix-fracture-tunnel water transport coupling according to claim 1, characterized in that: In step (S3), different three-dimensional seepage models are adopted to calculate the seepage field and permeability according to the different materials of the measured mine matrix.

3. The method for numerical simulation of matrix-fracture-tunnel water transport coupling according to claim 1, characterized in that: The effective cross section Z in formula (9) and formula (11) in step (S6) is calculated as formula (10): Z = Z b +h / 2, where Z b is the bottom elevation of the tunnel and h is the water depth.

4. The method for numerical simulation of matrix-fracture-tunnel water transport coupling according to claim 1, characterized in that: In the step (S6), the water flow velocity is calculated by the Xie Cai-Manning formula from the formula (9) and the formula (11), and then the water flow velocity is substituted into the Darcy-Weisbach formula to obtain the value.

5. The method for numerical simulation of matrix-fracture-tunnel water transport coupling according to claim 1, characterized in that: In step (S2), the two-dimensional unit grid simulating the fracture surface is calculated using two-dimensional units without thickness.

6. The method for numerical simulation of matrix-fracture-tunnel water transport coupling according to claim 1, characterized in that: In step (S7), coupling surfaces at the same pressure position are selected to couple coupling surfaces from the three-dimensional matrix to the two-dimensional crack, from the two-dimensional crack to the one-dimensional lane, and from the three-dimensional matrix to the one-dimensional lane; On the coupling surface, the sum of the flows in different dimensions equals the total flow.

7. The method for numerical simulation of matrix-fracture-tunnel water transport coupling according to claim 1, characterized in that: The three-dimensional visualization tool in the step (S8) is one of ParaView, VisIt, MATLAB, and Blender.

8. The method for numerical simulation of matrix-fracture-tunnel water transport coupling according to claim 1, characterized in that: In the coupling process between the matrix and the lane in step (S7), the fluid flow in the matrix and the lane is Darcy flow, so the water flow u entering the lane is tm It is expressed as formula (13): tm =u 达西 ·n t (13), where n t is the inner normal unit vector of the tunnel wall, pointing from the matrix to the inside of the tunnel, u 达西 represents the Darcy velocity of the fluid in the matrix in the porous medium, u 达西 Calculated by formula (6) in step (S5).

9. A method for numerical simulation of matrix-fracture-tunnel water transport coupling according to claim 8, characterized in that: When the fluid flow in the application scenario is non-Darcy flow, u in formula (13) 达西 is replaced by the non-Darcy velocity of the fluid in the matrix in the porous medium and calculated by formula (7) in step (S5).

Images