Numerical simulation method for matrix-fracture-roadway water migration coupling

Abstract

The invention relates to the technical field of intelligent mine mining, in particular to a numerical simulation method for matrix-fracture-roadway water transportation coupling. The method specifically comprises the following steps: extracting matrix and roadway information from a GIS system, and generating simplified one-dimensional roadway grid information including a three-dimensional matrix and a two-dimensional crack; physical property parameters, pressure and speed boundary conditions are set according to actually measured mine conditions; solving the permeability of the three-dimensional matrix and the two-dimensional fracture; carrying out discrete solution on the seepage equation of the three-dimensional matrix and the two-dimensional fracture and the Bernoulli equation of the roadway, and setting coupling conditions, namely, carrying out coupling calculation according to the principle that the pressure intensity is the same and the flow meets superposition, so as to obtain flow velocity information under different dimensions; after one-dimensional roadway rapid calculation, the water inrush flow and the water inrush level are obtained, and on this basis, local visual display and early warning are conducted through an SPH method. By means of the method, the migration condition of underground water can be completely simulated, and then the purpose of optimizing engineering design is achieved.

Description

Technical Field

[0001] The present invention relates to the field of smart mining technology, and in particular to a numerical simulation method for matrix-crack-tunnel water migration coupling. Background Art

[0002] Water migration in the mining process is a complex process, involving multiple links such as groundwater recharge, flow, influx into mines, drainage and excretion. It is mainly affected by factors such as geological structure, rock permeability, groundwater level and mining activities. Among them, geological structures such as faults, joints, and fissures in the geological structure affect the flow path and speed of groundwater; and the permeability of rocks determines the flow speed and path of groundwater; rocks with high permeability are prone to cause groundwater to quickly influx into mines, and the height of the groundwater level directly affects the amount and location of mine water influx. High groundwater levels increase the risk of mine water influx; mining activities change the natural flow path of groundwater, which may lead to a drop in groundwater levels and mine water influx. Therefore, the analysis and simulation methods of water migration have always been one of the key points in the construction of smart mines.

[0003] Due to the complexity of coupling and the high heterogeneity of fractured rock mass, numerical simulation using mathematical models is an important means to analyze the seepage mass transfer process of fractured rock mass. At present, the water migration process involves a variety of different regional conditions, forming a multi-region coupling system, which is more complex and changeable, and difficult to accurately describe with a single model. Moreover, the simulation for such situations is often just a simple numerical simulation scheme within the seepage field, ignoring the complete water migration in the formation mining, resulting in the inability to accurately and comprehensively judge the water hazard risk in the mining process. Summary of the invention

[0004] The purpose of the present invention is to overcome the deficiencies in the prior art and to provide a numerical simulation method for matrix-fracture-tunnel water migration coupling. This method analyzes the water migration situation through the pressure and flow information in each area according to the actual mining situation, and can completely simulate the migration of groundwater, accurately predict the seepage field, and reveal the fluid-solid interaction mechanism, thereby achieving the purpose of assessing environmental risks and optimizing engineering design.

[0005] In order to solve the problems in the prior art, the present invention discloses a numerical simulation method for matrix-crack-tunnel water migration coupling, which comprises the following steps:

[0006] Step 1: Establish the geometric model of the measured mine: According to the actual strata and tunnel geometric model of the measured mine on the GIS system (Geographic Information 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 along the axis. Then, determine the permeability based on 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 field information to determine the geometric model of the measured mine.

[0007] Step 2: Generate a mesh model for the entire computational domain: Import the geometric model determined in the first step into the mesh algorithm to generate a mesh model for the entire computational domain, and insert two-dimensional units into the stratigraphic information of the mesh model to simulate the fracture surface. As a weak medium in the rock mass, the fracture network has high permeability and low strength material properties, which must be fully considered in the simulation, that is, using two-dimensional units without thickness 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: Different three-dimensional seepage models are used to calculate the seepage field and permeability according to the different measured mine matrices.

[0009] 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 p represents the Darcy velocity of the fluid in the porous medium, t represents time, and q m Volume source term for the matrix.

[0010] 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.

[0011] 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, and ε is the porosity of the matrix.

[0012] Step 4: Calculate the seepage field and permeability of the fracture: The discrete fracture model is used to perform numerical simulation on the fracture, 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, uf represents the flow velocity of the fluid in the fracture, q f is the volume source term of the fracture.

[0013] The permeability k of the fracture f It is expressed as formula (5): where d f is the crack thickness.

[0014] Step 5: Numerical simulation of the matrix and fracture models based on two models of fluid flow through porous media:

[0015] 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;

[0016] 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;

[0017] 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 非达西 .

[0018] Step 6: Use the one-dimensional open channel flow model to numerically simulate the tunnel: Simplify the tunnel cross section into 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 simulate it using the Saint-Venant equations.

[0019] Its continuity equation is expressed as formula (8):

[0020] Its momentum equation is expressed as formula (9): Where A is the cross-sectional area, which is determined by the cross-sectional shape of the tunnel and the flooding depth, t is the time, Q is the flow rate through the section, x is the position of the tunnel, u is the average flow velocity of the section, g is the gravitational acceleration, Z is the effective section, that is, the relative height between the centroid of the flow section and the reference horizontal plane, and f is the sum of the along-the-way and local resistance caused by friction.

[0021] For rectangular tunnels, the continuity equation is simplified to formula (10):

[0022] The momentum equation is simplified to formula (11): Where w is the tunnel width, h is the water depth, t is the time, Q is the flow rate through the section, u is the average flow velocity in the section, and f is the sum of the along-the-way and local resistance caused by friction.

[0023] Step 7: Matrix-fracture-tunnel water transport coupling process:

[0024] 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:

[0025] 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.

[0026] B. Coupling process between matrix and lane:

[0027] 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 It is expressed as formula (13): tm =u p ·n t (13), where n t is the inner normal unit vector of the tunnel wall, pointing from the matrix to the interior of the tunnel.

[0028] 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, u tm is the water flow entering the tunnel obtained in formula (13).

[0029] C. The coupling process between cracks and tunnels:

[0030] 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 interior of the tunnel.

[0031] 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 along-the-way and local resistance caused by friction, and s is the specific surface.

[0032] Step 8: Visualization display and setting threshold warning: Use three-dimensional visualization tools to display the water flooding depth by the flow velocity, cross-sectional shape and size in the tunnel of step (S7), and visualize the simulation results of matrix-crack-tunnel; and perform local refined display of tunnel water inrush based on the SPH method. According to the display situation, analyze the water inrush level and set the water inrush level threshold, and set the disaster warning in the formation according to the water inrush level threshold.

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

[0034] Preferably, the effective cross section Z in formula (9) and formula (11) in the sixth step 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.

[0035] Preferably, in the ninth step, the water velocity is calculated by the Xie Cai-Manning formula in formula (9) and formula (11), and then the water velocity is substituted into the Darcy-Weisbach formula to obtain f.

[0036] Preferably, the two-dimensional unit grid simulating the fracture surface in the second step is calculated using two-dimensional units without thickness.

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

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

[0039] Preferably, 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 p ·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 p represents the Darcy velocity of the fluid in the matrix in the porous medium, u p The formula (6) in step (S5) is used for calculation.

[0040] More preferably, when the fluid flow in the application scenario is non-Darcy flow, u in formula (13) p 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).

[0041] When coupling across dimensions, it generally includes steps such as defining coupling surfaces, setting boundary conditions, setting source terms, solving coupling models, and verifying optimization. In this application, we define the coupling surfaces between the three-dimensional matrix, two-dimensional cracks, and one-dimensional tunnels, and obtain boundary conditions and initial conditions based on survey and field information, that is, the contents shown in the first and second steps; the third, fourth, and fifth steps set the source terms of the matrix and cracks respectively; the sixth step sets the source term of the tunnel; the seventh step couples the three dimensions in pairs, where the three-dimensional matrix to the two microcracks is used for flow transfer through the mass source term, the two microcracks and the one-dimensional tunnel are used for flow transfer through the average Darcy velocity at the cross section of the crack-tunnel, and the three-dimensional matrix and the one-dimensional tunnel are used for flow transfer through the average Darcy velocity at the cross section of the matrix-tunnel. The specific coupling calculation can use numerical methods, such as finite difference method, finite element method, finite volume method, etc. to solve the coupling model until a definite and optimal solution is obtained. Then, through the deduction of this method, the water migration situation of the matrix-fracture-tunnel is displayed in the three-dimensional visualization tool. Then we can analyze the water inrush level and set disaster warning in the formation according to the water inrush level threshold, ultimately achieving the purpose of smart mining and early warning.

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

[0043] 1. The complex interaction mechanism between matrix, fractures and tunnels is revealed, which helps to more accurately predict the behavior and distribution of water migration under different conditions (such as different hydraulic gradients, fracture widths, matrix porosity, etc.).

[0044] 2. It can effectively avoid and assess 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 movement under different design schemes, we can optimize the project layout, reduce water resource waste and improve project efficiency.

[0046] 4. Through intuitive graphic display and comprehensive visualization, the temporal and spatial variation patterns and characteristic analysis of water migration can be clearly observed, thereby better building an intelligent mining system. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1 is a flow chart of the present invention;

[0048] Figure 2 It is a flow chart of the matrix-fracture-tunnel water transport coupling scheme in the present invention;

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

[0050] Figure 4 is the pressure contour of the terrain cross section in the present invention;

[0051] Figure 5 It is a local and detailed display of the tunnel water inrush in the present invention. DETAILED DESCRIPTION

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

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

[0054] Step 1: Establish the geometric model of the measured mine: According to the actual strata and tunnel geometric model of the measured mine on the GIS system (Geographic Information 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 along the axis. Then, determine the permeability based on 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 field information to determine the geometric model of the measured mine.

[0055] Step 2: Generate a mesh model for the entire computational domain: Import the geometric model determined in the first step into the mesh algorithm to generate a mesh model for the entire computational domain, and insert two-dimensional units into the stratigraphic information of the mesh model to simulate the fracture surface. As a weak medium in the rock mass, the fracture network has high permeability and low strength material properties, which must be fully considered in the simulation, that is, using two-dimensional units without thickness 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: Different three-dimensional seepage models are used to calculate the seepage field and permeability according to the different measured mine matrices.

[0057] 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 p represents the Darcy velocity of the fluid in the porous medium, t represents time, and q m Volume source term for the matrix.

[0058] 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.

[0059] 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, and ε is the porosity of the matrix.

[0060] In actual operation, different models can be established according to different matrix materials and different classifications of fluids, and then different numbers are used to represent the models in different situations. When used later, the corresponding models can be directly selected.

[0061] Step 4: Calculate the seepage field and permeability of the fracture: The discrete fracture model is used to perform numerical simulation on the fracture, 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.

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

[0063] Step 5: Numerical simulation of the matrix and fracture models based on 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 fluid dynamic viscosity coefficient, and p represents the pressure;

[0065] 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;

[0066] 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 非达西 .

[0067] Step 6: Use the one-dimensional open channel flow model to numerically simulate the tunnel: Simplify the tunnel cross section into 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 simulate it using the Saint-Venant equations.

[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, which is determined by the cross-sectional shape of the tunnel and the flooding depth, t is the time, Q is the flow rate through the section, x is the position of the tunnel, u is the average flow velocity of the section, g is the gravitational acceleration, Z is the effective section, that is, the relative height between the centroid of the flow section and the reference horizontal plane, and f is the sum of the along-the-way and local resistance caused by friction.

[0070] For rectangular tunnels, the continuity equation is simplified to formula (10):

[0071] The momentum equation is simplified to formula (11): Where w is the tunnel width, h is the water depth, t is the time, Q is the flow rate through the section, u is the average flow velocity in the section, and f is the sum of the along-the-way and local resistance caused by friction.

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

[0073] 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:

[0074] 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.

[0075] B. Coupling process between matrix and lane:

[0076] 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 entering the tunnel is expressed as formula (13): u tm =u p ·n t (13), where n t is the inner normal unit vector of the tunnel wall, pointing from the matrix to the interior of the tunnel.

[0077] 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.

[0078] C. The coupling process between cracks and tunnels:

[0079] 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 velocity boundary condition of the tunnel water inrush as expressed in formula (15): u 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 interior of the tunnel.

[0080] 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 tunnel.

[0081] Step 8: Visualization display and setting threshold warning: Use three-dimensional visualization tools to display the water flooding depth by the flow velocity, cross-sectional shape and size in the tunnel of step (S7), and visualize the simulation results of matrix-crack-tunnel; and perform local refined display of tunnel water inrush based on the SPH method. According to the display situation, analyze the water inrush level and set the water inrush level threshold, and set the disaster warning in the formation according to the water inrush level threshold.

[0082] Obviously, 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 p ·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 p represents the Darcy velocity of the fluid in the matrix in the porous medium, u p The formula (6) in step (S5) is used for calculation.

[0083] When the fluid flow in the application scenario is non-Darcy flow, u in formula (13) p 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).

[0084] In simple terms, we extract stratum and tunnel information from the GIS system, thereby generating a grid of three-dimensional porous matrix, two-dimensional plane fractures and simplified one-dimensional tunnel information; then we set the physical parameters of the matrix, fluid and fracture according to the specific mining area, including the location and size information of the stratum, fracture and tunnel, 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 upper aquifer; then according to the porosity and rock pore size, the permeability of the porous matrix is ​​set according to Kozeny-Carman, the fracture thickness is obtained and its permeability is set according to the cubic law; through the constructed three-dimensional stratum model, two-dimensional plane fracture and simplified one-dimensional tunnel, as well as the corresponding physical parameters and boundary conditions, the seepage equations of the matrix and fracture and the Bernoulli equation of the tunnel are discretely solved, and the coupling conditions are set, that is, the pressure is the same and the flow satisfies the principle of superposition for coupling calculation, and the flow velocity information in different dimensions is obtained; finally, the water level of the water inrush is analyzed and the disaster warning in the stratum is set according to the water inrush level threshold and visualized.

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

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

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

[0088] The porosity of the formation matrix was set to 0.15 based on the physical properties of fine-grained rock.

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

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

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

[0092] The pressure in the tunnel is 1atm, which is one standard atmosphere.

[0093] After model deduction, in 1E 3 After s, select the formation section in the yz direction to display its total Darcy velocity streamline Figure 3 And pressure contours Figure 4 .

[0094] The total normal flow rate from the tunnel surface to the inside of the tunnel per unit time is 2.68518E - 9 m 3 / s.

[0095] The water depth of the tunnel is calculated to be 2.2E after one-dimensionalization. -3 m, the tunnel water inrush situation can be displayed through visualization, such as Figure 5 shown.

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

[0097] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the technical principles of the present invention. These improvements and modifications should also be regarded as the scope of protection 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 p 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 It is expressed as formula (13): tm =u p ·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 p 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 (S9), the water flow velocity is calculated by the Xie Cai-Manning formula and the f in the formula (9) and the formula (11) is then substituted into the Darcy-Weisbach formula to obtain the water flow velocity.

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 p ·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 p represents the Darcy velocity of the fluid in the matrix in the porous medium, u p The formula (6) in step (S5) is used for calculation.

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) p 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