Simulation method for predicting water inrush from mining-induced fault

Abstract

The application relates to a simulation method for predicting water inrush from a mining fault, which can quantitatively describe the permeation anisotropy, the permeation nonlinearity and the permeation heterogeneity of a mining fault crack, effectively avoids the subjectivity of the process of adjusting the permeation coefficient, and further accurately describes the distribution evolution of a mining floor fault crack, so as to provide an effective means for formulating optimized fault water prevention measures and parameters.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of mine water prevention and treatment, in particular to a simulation method for predicting water inrush from a mining-induced fault. BACKGROUND

[0002] The Carbon-Permian coal seam is one of the main coal seams of the North China type coalfield, but it is often directly covered by the Ordovician strong karst aquifer. According to statistics, 57 billion tons of coal resources in this region are threatened by water damage and cannot be mined, and more than 80% of the floor water inrush is related to faults. Production practice shows that the fault zone is often filled with sandy and muddy components, and under the influence of mining, tensile or shear fractures occur in the fault filling, which leads to the entry of the floor confined water into the working face along the fault fracture. The development scale, occurrence, coal seam thickness, and water abundance of the aquifer of the fault zone all affect the fault water inrush. It is crucial to predict the water inrush from the mining-induced fault using a reasonable numerical simulation method for working face layout and formulating mine water prevention and control measures.

[0003] In recent years, domestic and foreign scholars mainly use the finite element (FEM) and finite difference (FDM) methods to simulate and study the water inrush from mining-induced faults. Both FEM and FDM methods regard fault mud as a "porous medium", and mining stress will not fundamentally destroy the "pore structure", but only cause changes in "pore size or porosity", thereby changing the "permeability coefficient" and ultimately affecting the fault water inrush. In fact, mining stress not only leads to changes in the size of the fault mud pores, but also causes the formation of fractures. Randomly generated fractures with different openings lead to significant permeability heterogeneity, permeability anisotropy, and permeability nonlinearity in mining-induced faults. Based on the continuous medium assumption, the pore medium model framework, and the concept of "permeability coefficient", it is difficult to completely describe the above three characteristics of the permeability of mining-induced faults using FEM and FDM methods, which leads to a significant deviation of the numerical prediction of water inrush from the actual results. SUMMARY

[0004] In order to overcome at least one of the deficiencies in the prior art, the present application provides a simulation method for predicting water inrush from a mining-induced fault.

[0005] In a first aspect, a simulation method for predicting water inrush from a mining-induced fault is provided, comprising:

[0006] A mining-induced fault numerical model is established, the mining-induced fault numerical model is divided into a plurality of solid elements, and fracture elements are generated between adjacent solid elements; each solid element has elastic and plastic mechanical properties, and each fracture element has ductile fracture mechanical properties;

[0007] The working face recovery process is simulated in the mining-induced fault numerical model according to the working face recovery speed;

[0008] For each solid element, determine the fracture aperture of each fracture element adjacent to the solid element;

[0009] An Euler mesh component is established in the area around the mining fault in the numerical model of the mining fault, and a boundary condition is set; the Euler mesh component includes a plurality of Euler elements, and the Euler elements have fluid dynamics properties;

[0010] Based on the fluid volume VOF method and the fracture aperture of each fracture element, the volume fraction of fluid in each Euler element with gas-water two-way flow in the Euler mesh component is calculated, and the water inrush amount of the mining fault is determined according to the volume fraction of fluid in each Euler element.

[0011] In one embodiment, for each solid element, the fracture aperture of each fracture element adjacent to the solid element is determined, including:

[0012] Determine the solid elements adjacent to the solid element in various directions, and determine all fracture elements between the adjacent solid elements; the solid elements and the fracture elements each include a plurality of nodes;

[0013] Determine the fracture type of the fracture element according to the ductile fracture mechanics properties of the fracture element, and the fracture type includes a pure tensile fracture, a pure shear fracture or a tensile-shear mixed fracture;

[0014] If the fracture type of the fracture element is a pure tensile fracture, the fracture aperture does not need to be corrected by shear dilation, and the fracture aperture of the fracture element is obtained in the following manner:

[0015] Calculate the distance between each pair of corresponding nodes in the fracture element, and obtain the fracture aperture of the fracture element by using a linear interpolation method;

[0016] If the fracture type of the fracture element is a pure shear fracture or a tensile-shear mixed fracture, the fracture aperture needs to be corrected by shear dilation, and the fracture aperture of the fracture element is obtained in the following manner:

[0017] The fracture aperture Δn of the fracture element is calculated by using the following formula when the initial fracture aperture is set as Δn0=0:

[0018]

[0019] Wherein, S is the shear displacement, k1 is a material parameter, σ n is the normal stress, σ U is the uniaxial compressive strength of the adjacent rock mass, and ψ0 is the initial shear dilation angle.

[0020] In one embodiment, the fracture type of the fracture element is determined according to the ductile fracture mechanics properties of the fracture element, including:

[0021] If F e,n > 0 and Fe,s = 0, F e,t = 0, then the crack element is a pure tensile type crack; wherein, F en is the tensile fracture energy, F es , F et is the fracture energy of two shear directions;

[0022] If F e,s > 0, or F e,t > 0, then the crack element is a pure shear type crack, or a tensile-shear mixed type crack.

[0023] In an embodiment, based on the fluid volume VOF method and the crack opening of each crack element, the volume fraction of the fluid in each Euler element with gas-water two-way flow in the Euler grid component is calculated, and the water inrush amount of the mining fault is determined according to the volume fraction of the fluid in each Euler element, comprising:

[0024] Based on the crack opening of each crack element, combined with the enhanced immersed boundary algorithm and Hugoniot condition, the flow-solid contact boundary is identified, and the normal velocity, tangential velocity, pressure and density of the fluid on both sides of the flow-solid contact boundary are determined;

[0025] According to the normal velocity, tangential velocity, pressure and density of the fluid, the volume fraction of the fluid in each Euler element with gas-water two-way flow in the Euler grid component is calculated based on the fluid volume VOF method;

[0026] The volume of each Euler element with gas-water two-way flow is multiplied by the volume fraction of the fluid to obtain the water inrush amount in each Euler element;

[0027] The water inrush amounts of all Euler elements are accumulated to obtain the water inrush amount of the mining fault.

[0028] In an embodiment, the mining fault numerical model is a two-dimensional model or a three-dimensional model.

[0029] In a second aspect, a simulation device for predicting the water inrush amount of a mining fault is provided, comprising:

[0030] A model establishing module is configured to divide the mining fault numerical model into a plurality of entity elements, and generate crack elements between adjacent entity elements; each entity element has elastic mechanical properties and plastic mechanical properties, and each crack element has ductile fracture mechanical properties;

[0031] A working face mining process simulation module is configured to simulate the working face mining process in the mining fault numerical model according to the working face mining speed;

[0032] A crack opening determination module is configured to determine the crack opening of each crack element adjacent to each entity element for each entity element;

[0033] an Eulerian mesh part establishing module, configured to establish an Eulerian mesh part around a region of the mining fault numerical model and set a boundary condition; the Eulerian mesh part comprises a plurality of Eulerian units, and each of the Eulerian units has a fluid dynamics attribute;

[0034] a water inrush amount determining module, configured to calculate a volume fraction of fluid in each of the Eulerian units having gas-water two-way flow in the Eulerian mesh part based on a volume of fluid (VOF) method and a fracture opening of each of the crack units, and determine a water inrush amount of the mining fault according to the volume fraction of fluid in each of the Eulerian units.

[0035] In one embodiment, the fracture opening determining module is further configured to:

[0036] determine solid units adjacent to the solid unit in each direction and all crack units between the adjacent solid units; each of the solid units and the crack units comprises a plurality of nodes;

[0037] determine a fracture type of the crack unit according to a ductile fracture mechanics attribute of the crack unit, the fracture type comprising a pure tensile fracture, a pure shear fracture or a tensile-shear mixed fracture;

[0038] if the fracture type of the crack unit is the pure tensile fracture, the fracture opening does not need to be corrected by shear dilation, and the fracture opening of the crack unit is obtained in the following manner:

[0039] calculate a distance between each two nodes corresponding to the crack unit, and obtain the fracture opening of the crack unit by using a linear interpolation method;

[0040] if the fracture type of the crack unit is the pure shear fracture or the tensile-shear mixed fracture, the fracture opening needs to be corrected by shear dilation, and the fracture opening of the crack unit is obtained in the following manner:

[0041] if the initial fracture opening is set as Δn0=0, the fracture opening Δn of the crack unit is calculated by using the following formula:

[0042]

[0043] wherein S is a shear displacement, k1 is a material parameter, σn is a normal stress, σt is a shear stress, and ψ0 is an initial shear dilation angle. n U

[0044] In one embodiment, the water inrush amount determining module is further configured to:

[0045] based on the fracture opening of each of the crack units, identify a fluid-solid contact boundary in combination with an enhanced immersed boundary algorithm and Hugoniot conditions, and determine normal velocity, tangential velocity, pressure and density of fluid on both sides of the fluid-solid contact boundary.​​

[0046] According to the normal velocity, tangential velocity, pressure, density of the fluid, the volume fraction of the fluid in each Euler cell with gas-water two-way flow in the Euler grid component is calculated based on the fluid volume VOF method;

[0047] The volume of each Euler cell with gas-water two-way flow is multiplied by the volume fraction of the fluid to obtain the water inrush amount in each Euler cell;

[0048] The water inrush amounts of all Euler cells are accumulated to obtain the water inrush amount of the mining fault.

[0049] In a third aspect, a computer readable storage medium is provided, and the computer readable storage medium stores a computer program, and the computer program is executed by a processor to implement the simulation method for predicting the water inrush amount of the mining fault.

[0050] In a fourth aspect, a computer program product is provided, and the computer program product includes computer programs / instructions, and the computer programs / instructions are executed by a processor to implement the simulation method for predicting the water inrush amount of the mining fault.

[0051] Compared with the prior art, the present application has the following beneficial effects:

[0052] 1. The present application provides an effective simulation tool for floor fault rupture and water inrush in the process of underground coal mining, which can realize the initiation, expansion, penetration, tension or shear dislocation of the mining floor fault fracture, and at the same time realize the migration of confined water in the fault fracture; the method is a major improvement over the methods based on "continuous medium assumption" and "porous medium seepage model" such as finite element FEM and finite difference FDM, which can directly assign the experimental results of fracture, shear friction and fracture flow to the numerical model without adjusting parameters such as "permeability coefficient"; at the same time, it effectively overcomes the numerical simulation difficulties of heterogeneity, anisotropy and high nonlinearity of fault fracture flow.

[0053] 2. The present application proposes a shear dilatancy correction method for the opening of shear-type fractures to solve the problem that the two sides of the compression shear-type fracture are easy to embed each other, which is contrary to the results of shear experiments, and provides an effective method for shear-tension mixed fault water inrush of shear-type faults.

[0054] 3. The present application can realize the arrangement of Euler grid area of any size in any mining fracture development area, thereby greatly speeding up the numerical calculation speed.

[0055] 4. The present application combines the enhanced immersed boundary algorithm, Hugoniot condition and fluid volume VOF method, which not only limits the fluid in the mining fracture, but also realizes the tracking and reconstruction of the free surface of the fluid in the fracture.

[0056] 5、The application can be used for analyzing water flow velocity in different fractures under different engineering and geological parameters, obtaining evolution law of fault water inrush volume with mining parameters, thereby optimizing working face layout and formulating fault water prevention and control measures. BRIEF DESCRIPTION OF DRAWINGS

[0057] The application can be better understood by referring to the following description in conjunction with the accompanying drawings, which are incorporated in and form a part of the specification, wherein:

[0058] Figure 1 A flow chart of a simulation method for predicting water inrush volume of a mining fault according to an embodiment of the application is shown;

[0059] Figure 2 A schematic diagram of grid division of a mining fault numerical model is shown;

[0060] Figure 3 A schematic diagram of simulation of a mining working face mining process is shown;

[0061] Figure 4 A shear dilation correction schematic diagram is shown;

[0062] Figure 5 A mining floor fault water inrush process schematic diagram is shown;

[0063] Figure 6 A structural block diagram of a simulation device for predicting water inrush volume of a mining fault according to an embodiment of the application is shown. DETAILED DESCRIPTION

[0064] In the following, exemplary embodiments of the application will be described with reference to the accompanying drawings. In the description of the embodiments, not all features of the actual embodiments are described. It should be appreciated, however, that many embodiment-specific decisions can be made in the development of any such actual embodiment, in order to achieve the specific goals of the developer, and these decisions can vary from embodiment to embodiment.

[0065] It should also be noted here that, in order to avoid obscuring the application due to unnecessary details, only the device structure closely related to the scheme according to the application is shown in the drawings, and other details not closely related to the application are omitted.

[0066] It should be understood that the application is not limited to the described embodiments by virtue of the following description with reference to the drawings. In this context, embodiments can be combined with each other, features can be replaced or borrowed between different embodiments, and one or more features can be omitted in an embodiment.

[0067] The application provides a simulation method for predicting water inrush from a mining fault, which can quantitatively describe the fracture permeation heterogeneity, permeation anisotropy and permeation nonlinearity of the mining fault, effectively avoid the subjectivity of the process of adjusting the "permeation coefficient", and further accurately describe the distribution evolution of the mining floor fault fracture, thereby providing an effective means for formulating optimized fault water prevention measures and parameters.

[0068] Figure 1 A flow chart of a simulation method for predicting water inrush from a mining fault according to an embodiment of the application is shown, and the method comprises the following steps:

[0069] Step S1, a mining fault numerical model is established, the mining fault numerical model is divided into a plurality of entity units, and a crack unit is generated between adjacent entity units; each entity unit has elastic mechanics properties and plastic mechanics properties, and each crack unit has ductile fracture mechanics properties.

[0070] Here, first, a mining fault numerical model is established according to the occurrence of the mining fault, the position of the aquifer and the characteristics of the working face to be arranged.

[0071] Then, the mining fault numerical model is globally meshed to obtain a plurality of entity units. Here, according to different accuracy requirements of the numerical model, the mining fault numerical model can be a two-dimensional model or a three-dimensional model. When the mining fault numerical model is a two-dimensional model, each entity unit includes four nodes, and when the mining fault numerical model is a three-dimensional model, each entity unit includes eight nodes. The size D (i.e. the side length) of each entity unit satisfies the following condition:

[0072] D≤0.28πEF e / [(1-μ 2 )σ n ]

[0073] Wherein, π=3.14; E is the elastic modulus; F e is the fracture energy; μ is the Poisson's ratio; σ n is the tensile strength.

[0074] Then, a crack unit is generated between adjacent entity units, and the specific method comprises:

[0075] Each entity unit and each node in the entity unit are sequentially numbered to obtain the unit number and the node number of the entity unit;

[0076] The node is copied n times while keeping the coordinates of the nodes of the entity unit unchanged, n represents that the node is shared by n entity units, and the node is renumbered:

[0077] NODE new =NODE+10i+1

[0078] Among them, NODE new NODE is the new node number, where NODE is the original node number and i is the maximum value of the original node number. max The number of decimal places can be determined using the formula w = NODE. max / 10 i Calculate i, w are the set parameters used for decimal number calculation, 0.1≤w<10.

[0079] The initial element number for the crack element is determined to be ELEMENT. max +1, ELEMENT max The largest cell number of the entity;

[0080] Number the new nodes of multiple nodes with the same coordinates in adjacent solid elements. new By reversing the order and using it as the node number in the crack element, crack elements can be generated between adjacent solid elements.

[0081] The elastic mechanical properties of solid elements satisfy the following relationship:

[0082] The expression for elastic deformation is:

[0083] σ ij =λε kk δ ij +2Gε ij

[0084] Where, σ ij For stress components, ε ij Let λ be the strain components (i,j = 1, 2, 3), i be the direction of the surface, j be the direction of the force, λ be the Lamé constant, G be the shear modulus, and δ be the strain components (i,j = 1, 2, 3). ij For a unit tensor, ε kk (k = 1, 2, 3) is ε ij The strain components at i = j are the elements on the main diagonal of the strain matrix.

[0085] The plastic mechanical properties of solid elements satisfy the following relationship:

[0086] ① Yield criterion:

[0087]

[0088] Where c is the cohesive force. Let θ be the internal friction angle, and τ and σ be the shear stress and principal stress on the slip surface, respectively.

[0089] ② The Law of Flow:

[0090]

[0091] wherein, is the plastic strain component, Φ is the plastic potential function, and dθ is the plastic factor;

[0092] Expression of Φ:

[0093] Φ = [(δσ t tanψ) 2 +(Bq) 2 ] 1 / 2 -ptanψ

[0094] wherein q is the deviatoric stress, p is the spherical stress, δ is the sharp-point curvature of the plastic potential function meridian in the tension zone of the q-p plane, and takes the value 0.1; ψ is the dilatancy angle; B is a parameter for controlling the shape of the plastic potential function G in the π plane, and σ t represents the tensile strength of the rock.

[0095]

[0096] wherein α is the azimuth angle on the π plane, and e is the eccentricity, and e = (3-sinα) / (3+sinα).

[0097] The crack element has the property of ductile fracture mechanics, and satisfies the following relationship:

[0098] In the elastic deformation stage of fracture, the relationship between the traction force and the separation displacement is:

[0099]

[0100] wherein t is the load of the crack element, t n , t s , t t are the loads in the normal direction and two tangential directions of the crack element, E nn , E ss , E tt are the stiffnesses in the normal direction and two tangential directions of the crack element, δ n , δ s , δ t are the displacements in the normal direction and two tangential directions of the crack element, and Eδ is the traction force.

[0101] When the stress satisfies the following condition, the crack element begins to fracture:

[0102]

[0103] wherein t n 0 , t s 0 , tt 0 are the peak loads in the normal direction and two tangential directions of the crack element, respectively, n > represents that only tensile load is considered in the normal direction of the crack element, when t n <0, it is a compressive load, and take t n > = 0, that is, the compressive stress has no effect on the fracture of the crack element.

[0104] With the continuous increase of the separation displacement, the contact surface area between the rock blocks on both sides of the fracture surface gradually decreases, which leads to the decrease of the tensile and shear strength of the material. Therefore, damage is introduced to describe the above phenomenon. The relationship between the mixed type fracture displacement and the damage evolution equation is established:

[0105]

[0106] where d c is the damage variable, is the mixed type total displacement of the crack tip at a certain time, and the value is the vector sum of the pure tensile and pure shear displacement at that time. is the total displacement of the crack element when it is initially fractured. is the total fracture displacement of the crack element, which can be obtained by experiment.

[0107] The tensile and shear strengths of the crack element after damage are:

[0108]

[0109]

[0110]

[0111] where t′ n is the tensile strength of the crack element after damage, is the normal force on the crack element, t′ s is the shear strength of the first shear direction of the crack element after damage, is the shear force on the first shear direction of the crack element, t′ t is the shear strength of the second shear direction of the crack element after damage, is the shear force on the second shear direction of the crack element.

[0112] Whether the crack element is completely fractured is judged by the fracture critical displacement δ F under tensile and shear load:

[0113]

[0114] where δ n represents the fracture displacement in the normal direction; δ s , δt The fracture displacement of two shear directions (the first shear direction and the second shear direction).

[0115] F e,n The fracture displacement of two shear directions (the first shear direction and the second shear direction). e,s The fracture displacement of two shear directions (the first shear direction and the second shear direction). e,t The fracture displacement of two shear directions (the first shear direction and the second shear direction). The fracture displacement of two shear directions (the first shear direction and the second shear direction).

[0116] Figure 2 A schematic diagram of meshing of a numerical model of mining fault is shown. The model size is long x = 760 m, and high y = 388 m. The entire model is globally meshed by using 8-node hexahedral solid elements, and the size D of the solid element at the fault is controlled to be below 1 m.

[0117] The maximum node number of the solid element is NODE max = 63430, w = 5. The element nodes in the corner points in the model are copied once, the nodes at the boundary (non-corner points) are copied twice, and the nodes of the elements in the middle part are copied four times. Then, the nodes of all the elements are renumbered. For example, the original node number of the solid element is 3270, and there are four surrounding shared elements, then the new node numbers are 132701, 132702, 132703, and 132704 in turn.

[0118] ELEMENT max = 48321, and the initial number of the crack element is 48322. For the eight nodes of two adjacent solid elements (the element numbers of which are 22355 and 22357) and the same coordinates, the copied node numbers are 172421, 172417, 172404, 172423, 172419, 172406, 172402 in turn, and the crack element 49124 is assigned after being arranged counterclockwise. The keyword “*Element type = COH3D8” is added before the element number in EditKeyword, so that the crack element is generated between the adjacent solid elements.

[0119] The material parameters in the constitutive theory are assigned to the numerical calculation model, and the details are as follows:

[0120] Table 1 Mechanical parameters of complete rock mass in coal measures strata

[0121]

[0122]

[0123] Table 2 Fracture mechanics parameters

[0124]

[0125] Step S2, according to the working face mining speed, simulating the mining process of the coal mining face in the numerical model of the mining faulted strata.

[0126] Here, if the numerical model of the mining faulted strata is a two-dimensional model, the entire working face range is determined according to the two corner points of the lower left and the upper right, i.e. a total of 2 (x, y) coordinates, if the numerical model of the mining faulted strata is a three-dimensional model, the entire working face range is determined according to the two corner points of the front left and the rear right, i.e. a total of 2 (x, y, z) coordinates, according to the working face mining speed, part of the solid elements and the crack elements are deleted to simulate the mining process of the coal mining face.

[0127] Specifically, the range of the working face (i.e. the coal seam being mined) is determined by the diagonal corner points (x1, y1) = (647.76, -220) and (x2, y2) = (1256, -249), and in combination with the working face advancing speed of 4 m / d, the working face is advanced by 4 m every 1 d. Figure 3 A schematic diagram of simulating the mining process of the coal mining face is shown.

[0128] Step S3, for each solid element, the crack opening of each crack element adjacent to the solid element is determined.

[0129] Step S4, an Euler mesh component is established around the faulted area of the numerical model of the mining faulted strata, and a boundary condition is set, such as Figure 2 The white dashed box shown in the figure is the Euler network component; the Euler mesh component includes a plurality of Euler elements, and the Euler element has fluid dynamics properties. Here, the area around the fault can be a range of 100 m in front of the fault and 80 m behind the fault. The boundary condition can be to set the confined water pressure of the bottom boundary of the Euler mesh component to 1.4 MPa, and to set the Euler boundary with a speed of 0 on the top, front, back, left and right.

[0130] Specifically, the Euler mesh component has fluid dynamics properties, which satisfy the following relationship:

[0131] ① The gas-water two-phase flow mass conservation equation is:

[0132]

[0133] Where t is time, u is velocity vector, and p is the density of the Euler element.

[0134] The expression of p is:

[0135] p = p w V p + p a (1 - V p )

[0136] Where pw ρ is the density of water. a V is the density of air. p Let V be the volume fraction of water in the Euler unit, and 0 ≤ V. p ≤1, is a variable.

[0137] ②The equation for the conservation of gas-water two-phase flow is:

[0138]

[0139] Where, p wa ρ is the fluid pressure; g is the acceleration due to gravity; p is the fluid pressure; F st Let be the surface tension of water; τ be the stress tensor caused by viscous forces, and have...

[0140]

[0141] Where μ is the dynamic viscosity, and has

[0142] μ = μ w V p +μ a (1-V p )

[0143] Where, μ w μ represents the dynamic viscosity of water. a This indicates the dynamic viscosity of air.

[0144] ③ Equations of state:

[0145]

[0146] Where, p wa Where c is the fluid pressure, c0 is the material constant, η is the nominal volumetric compressive strain, η = 1 - ρ0 / ρ, and ρ0 is the reference density, taken as 1000 kg / m³. 3 .

[0147] Step S5: Based on the fluid volume VOF method and the fracture aperture of each fracture element, calculate the fluid volume fraction in each Euler element with gas-water bidirectional flow in the Euler grid component, and determine the mining-induced fault water inrush volume according to the fluid volume fraction in each Euler element.

[0148] In one embodiment, step S3, for each solid element, determines the crack aperture of each crack element adjacent to the solid element, including:

[0149] Step S31: Determine the solid elements that are adjacent to the solid elements in each direction, and determine all crack elements between adjacent solid elements; both solid elements and crack elements include multiple nodes.

[0150] Specifically, taking the integral point (x1, y1) of any one entity unit (denoted as ELEMENT0) in the working face mining process as the center, the radius r gradually expands from 0 until it first completely contains a certain entity unit (denoted as ELEMENT1, the integral point of which is (x2, y2)) around it, the integral points of the two units are connected to obtain the direction vector a1 = (x2-x1, y2-y1). Thus, ELEMENT1 is the entity unit adjacent to ELEMENT0 in the direction of a1. Similarly, a2, a3, a4 and other entity units adjacent to ELEMENT0 in multiple directions can be obtained; and all crack units between the adjacent entity units are determined, for example, four crack units between the adjacent entity units are determined, and the nodes contained in each crack unit are represented by a node set, which are {NODE0-1}, {NODE0-2}, {NODE0-3}, and {NODE0-4}, respectively.

[0151] Step S32, determining the crack type of the crack unit according to the toughness fracture mechanics properties of the crack unit, the crack type including a pure tensile crack, a pure shear crack, or a tensile-shear mixed crack.

[0152] Specifically, if F e,n > 0 and F e,s = 0, F e,t = 0, the crack unit is a pure tensile crack; wherein F e,n is a tensile fracture energy, F e,s , F e,t are two shear direction fracture energies.

[0153] If F e,s > 0, or F e,t > 0, the crack unit is a pure shear crack or a tensile-shear mixed crack.

[0154] Step S33, if the crack type of the crack unit is a pure tensile crack, the crack opening degree does not need to be corrected by shear dilation, and the crack opening degree of the crack unit is obtained in the following manner:

[0155] The distance between the corresponding two nodes in the crack unit is calculated, and the crack opening degree of the crack unit is obtained by linear interpolation; here, the corresponding two nodes belong to two adjacent entity units.

[0156] If the crack type of the crack unit is a pure shear crack or a tensile-shear mixed crack, the crack opening degree needs to be corrected by shear dilation, Figure 4 a shear dilation correction schematic diagram is shown. The crack opening degree of the crack unit is obtained in the following manner:

[0157] The initial crack opening degree is set as Δn0=0, and the crack opening degree Δn of the crack unit is calculated by the following formula:

[0158]

[0159] wherein S is the shear displacement, k1 is a material parameter, reflecting that the influence of shear displacement S on Δn gradually weakens with the increase of loading and unloading times, and has k1 = 0.89-0.23e -(Num-1.12) / 2.87 , Num represents the number of shear cycles, σ n is the normal stress, σ U is the uniaxial compressive strength of the adjacent rock mass, and ψ0 is the initial dilatancy angle.

[0160] In one embodiment, in step S5, based on the fluid volume VOF method and the fracture opening of each fracture element, the volume fraction of fluid in each Euler element with gas-water two-phase flow in the Euler grid component is calculated, and the water inrush amount of the mining fault is determined according to the volume fraction of fluid in each Euler element, including:

[0161] Step S51, based on the fracture opening of each fracture element, combined with the enhanced immersed boundary algorithm, Hugoniot condition, identify the flow-solid contact boundary, and determine the normal velocity, tangential velocity, pressure and density of the fluid on both sides of the flow-solid contact boundary.

[0162] Here, in the foregoing embodiment, after calculating the distance between the corresponding two nodes in the fracture element, the minimum value of the distance can be selected for use in the enhanced immersed boundary algorithm.

[0163] Step S52, according to the normal velocity, tangential velocity, pressure and density of the fluid, based on the fluid volume VOF method, the volume fraction of fluid in each Euler element with gas-water two-phase flow in the Euler grid component is calculated;

[0164] Step S53, the volume of each Euler element with gas-water two-phase flow is multiplied by the volume fraction of fluid to obtain the water inrush amount in each Euler element;

[0165] Step S54, the water inrush amounts of all Euler elements are accumulated to obtain the water inrush amount of the mining fault.

[0166] Figure 5 A schematic diagram of the water inrush process of the mining floor fault is shown, and the water inrush process of the mining floor fault is analyzed Figure 5The results show that with the forward advance of the working face, the mining damage envelope of the coal seam floor fault group is generally in the form of w, and the fault zone and its hanging wall are severely damaged. The deepest floor damage depth is located at the F2 fault and its hanging wall, with a damage depth of 48.6 m. The damage depth is the shallowest at the F2 fault footwall, with a damage depth of 23 m. The F1 fault damage depth is the second, with a damage depth of 42.3 m at the F1 fault and its hanging wall. The damage depth is 24.6 m at the F1 fault footwall. The damage depth of the non-structural floor is between 24-30.7 m. In addition, the maximum water flow velocity in the fissures of the fault zone and its hanging wall is 44.25 m / h and 87.49 m / h, respectively. Combined with the number of Euler elements and the fluid volume fraction of each Euler element, it is determined that the water inrush amount of Ordovician limestone water into the working face goaf is 118.74 m 3 / h. The above damage depth has exceeded the distance between the Ordovician limestone aquifer and the coal seam (30 m), and the above water inrush amount (2849 m 3 of water inrush in 20 h) has exceeded the upper limit of the working face drainage amount (2760 m 3 of drainage in 24 h), indicating that the mining floor fault has a water inrush risk, and there is a risk of flooding the working face. According to the numerical simulation results, a waterproof coal pillar with a size of 35 m needs to be arranged; or floor grouting is carried out at a depth of 40 m below the coal seam floor (i.e. 10 m or more below the top boundary of the Ordovician limestone), so as to prevent and control the Ordovician limestone water inrush disaster.

[0167] Based on the same inventive concept as the simulation method for predicting the water inrush amount of the mining fault, the embodiment also provides a simulation device for predicting the water inrush amount of the mining fault, Figure 6 a structural block diagram of the simulation device for predicting the water inrush amount of the mining fault according to the embodiment of the application is shown, which comprises:

[0168] The model establishing module 61 is used to divide the mining fault numerical model into a plurality of entity elements, and generate crack elements between adjacent entity elements. Each entity element has elastic mechanical properties and plastic mechanical properties, and each crack element has ductile fracture mechanical properties.

[0169] The working face mining process simulation module 62 is used to simulate the working face mining process in the mining fault numerical model according to the working face mining speed.

[0170] The fissure opening determination module 63 is used to determine the fissure opening of each crack element adjacent to the entity element for each entity element.

[0171] The Euler grid component establishing module 64 is used to establish an Euler grid component around the fault zone of the mining fault numerical model and set a boundary condition. The Euler grid component comprises a plurality of Euler elements, and the Euler element has fluid dynamics properties.

[0172] The water inrush amount determination module 65 is configured to calculate the volume fraction of the fluid in each Euler element with gas-water two-way flow in the Euler grid component based on the fluid volume VOF method and the fracture aperture of each fracture element, and determine the water inrush amount of the mining fault according to the volume fraction of the fluid in each Euler element.

[0173] The simulation device for predicting the water inrush amount of the mining fault in the embodiment has the same inventive concept as the simulation method for predicting the water inrush amount of the mining fault described above, and therefore the specific implementation of the device can be seen from the embodiment part of the simulation method for predicting the water inrush amount of the mining fault described above, and the technical effects thereof correspond to those of the method described above, which will not be repeated here.

[0174] The embodiment of the present application provides a computer readable storage medium, which stores a computer program, and the computer program is executed by a processor to implement the simulation method for predicting the water inrush amount of the mining fault.

[0175] The embodiment of the present application provides a computer program product, which includes computer programs / instructions, and the computer programs / instructions are executed by a processor to implement the simulation method for predicting the water inrush amount of the mining fault.

[0176] In summary, the present application has the following technical effects:

[0177] 1. The present application provides an effective simulation tool for floor fault rupture and water inrush in the process of underground coal mining, which can realize the initiation, expansion, penetration, tension or shear dislocation of the mining floor fault fracture, and the migration of confined water in the fault fracture at the same time. The method is a major improvement over the methods based on the "continuum hypothesis" and "porous medium seepage model" such as finite element FEM and finite difference FDM. The results of the experiment, such as fracture, shear friction and fracture flow, can be directly assigned to the numerical model without adjusting parameters such as "permeability coefficient". At the same time, the numerical simulation difficulties of the heterogeneity, anisotropy and high nonlinearity of the fault fracture flow are effectively overcome.

[0178] 2. The present application proposes a shear dilation correction method for the opening degree of the shear type fracture to solve the problem that the two rock blocks on both sides of the compression shear type fracture are easy to embed each other, which is contrary to the results of the shear experiment, and provides an effective method for the tension-shear mixed type fault water inrush of the shear type fault.

[0179] 3. The present application can realize the arrangement of the Euler grid area of any size in any mining fracture development area, thereby greatly speeding up the numerical calculation speed.

[0180] 4. The present application combines the enhanced immersed boundary algorithm, Hugoniot condition and fluid volume VOF method, which not only limits the fluid in the mining fracture, but also realizes the tracking and reconstruction of the free surface of the fluid in the fracture.

[0181] 5、The application can be used for analyzing water flow velocities in different fractures under different engineering and geological parameters, obtaining evolution rules of fault water inrush volume with mining parameters, and thus optimizing working face arrangement and formulating fault water prevention and control measures.

[0182] The above merely describes various embodiments of the application, but the protection scope of the application is not limited thereto, and any person skilled in the art can easily think of changes or replacements within the technical scope disclosed by the application, which should be encompassed in the protection scope of the application. Therefore, the protection scope of the application should be subject to the protection scope of the claims.

Claims

1. A simulation method for predicting water inrush amount of mining fault, characterized in that, The method comprises the following steps: establishing a numerical model of a mining fault, and dividing the numerical model of the mining fault into a plurality of solid elements; generating crack elements between adjacent solid elements; each of the solid elements has elastic mechanical properties and plastic mechanical properties, and each of the crack elements has ductile fracture mechanical properties; simulating a mining process of a coal mining face in the numerical model of the mining fault according to a mining speed of the coal mining face; for each of the solid elements, determining a crack opening of each of the crack elements adjacent to the solid element; establishing an Euler mesh component in a surrounding area of the fault in the numerical model of the mining fault, and setting a boundary condition; the Euler mesh component comprises a plurality of Euler elements, and the Euler elements have fluid dynamics properties; based on a volume of fluid (VOF) method and the crack opening of each of the crack elements, calculating a volume fraction of fluid in each of the Euler elements with gas-water two-way flow in the Euler mesh component, and determining a water inrush amount of the mining fault according to the volume fraction of fluid in each of the Euler elements; for each of the solid elements, determining a crack opening of each of the crack elements adjacent to the solid element, comprising: determining solid elements adjacent to the solid element in each direction, and determining all crack elements between the adjacent solid elements; the solid elements and the crack elements each comprise a plurality of nodes; determining a crack type of the crack element according to the ductile fracture mechanical properties of the crack element, wherein the crack type comprises a pure tensile crack, a pure shear crack or a tensile-shear mixed crack; if the crack type of the crack element is the pure tensile crack, the crack opening does not need to be corrected by shear dilation, and the crack opening of the crack element is obtained in the following manner: calculating a distance between two nodes corresponding to each other in the crack element, and obtaining the crack opening of the crack element by using a linear interpolation method; if the crack type of the crack element is the pure shear crack or the tensile-shear mixed crack, the crack opening needs to be corrected by shear dilation, and the crack opening of the crack element is obtained in the following manner: The initial opening of the fracture is set as Δ n 0 = 0, the fracture opening of the fracture element Δ n is calculated by the following formula: where, S is the shear displacement, k 1 is a material parameter, σ n is the normal stress, σ U is the uniaxial compressive strength of the adjacent rock mass, ψ 0 is the initial dilatancy angle.

2. The method of claim 1, wherein, wherein, determining the crack type of the crack element according to the ductile fracture mechanical properties of the crack element, comprising: If F e,n > 0 and F e,s = 0, F e,t = 0, the crack element is a pure tensile crack; wherein, F e,n is the tensile fracture energy, F e,s , F e,t is the fracture energy of two shear directions; If F e,s > 0 ， Or F e,t > 0 ， The crack unit is a pure shear type fracture, or a tensile-shear mixed type fracture.

3. The method of claim 1, wherein, wherein, based on the volume of fluid (VOF) method and the crack opening of each of the crack elements, calculating a volume fraction of fluid in each of the Euler elements with gas-water two-way flow in the Euler mesh component, and determining a water inrush amount of the mining fault according to the volume fraction of fluid in each of the Euler elements, comprising: based on the crack opening of each of the crack elements, combining a reinforced immersed boundary algorithm and Hugoniot conditions to identify a flow-solid contact boundary, and determining normal velocity, tangential velocity, pressure and density of fluid on both sides of the flow-solid contact boundary; based on the volume of fluid (VOF) method, calculating the volume fraction of fluid in each of the Euler elements with gas-water two-way flow in the Euler mesh component according to the normal velocity, the tangential velocity, the pressure and the density of the fluid; multiplying the volume of each of the Euler elements with gas-water two-way flow by the volume fraction of the fluid to obtain a water inrush amount in each of the Euler elements; accumulating the water inrush amounts of all the Euler elements to obtain the water inrush amount of the mining fault.

4. The method of claim 1, wherein, The mining fault numerical model is a two-dimensional model or a three-dimensional model.

5. A simulation device for predicting water inrush from a mining-induced fault, characterized in that, Comprise: The model establishment module is used for dividing the mining fault numerical model into a plurality of entity units; And generate crack units between adjacent entity units; Each of the entity units has elastic mechanical properties and plastic mechanical properties, and each of the crack units has ductile fracture mechanical properties; The working face mining process simulation module is used for simulating the mining process of the coal mining face in the mining fault numerical model according to the working face mining speed; The crack opening degree determination module is used for determining the crack opening degree of each crack unit adjacent to the entity unit for each of the entity units; The Euler grid component establishment module is used for establishing an Euler grid component around the area of the fault in the mining fault numerical model and setting a boundary condition; The Euler grid component comprises a plurality of Euler units, and the Euler units have fluid dynamics properties; The water inrush volume determination module is used for calculating the volume fraction of fluid in each Euler unit with gas-water two-way flow in the Euler grid component based on a fluid volume VOF method and the crack opening degree of each crack unit, and determining the mining fault water inrush volume according to the volume fraction of fluid in each Euler unit. The crack opening degree determination module is also used for: Determining the entity units adjacent to the entity unit in each direction and determining all crack units between the adjacent entity units; The entity units and the crack units each comprise a plurality of nodes; Determining the crack type of the crack unit according to the ductile fracture mechanical properties of the crack unit, wherein the crack type comprises a pure tensile crack, a pure shear crack or a tensile-shear mixed crack; If the crack type of the crack unit is a pure tensile crack, the crack opening degree does not need to be corrected by shear dilation, and the crack opening degree of the crack unit is obtained in the following manner: The distance between each two nodes corresponding to the crack unit is calculated, and the crack opening degree of the crack unit is obtained by using a linear interpolation method; If the crack type of the crack unit is a pure shear crack or a tensile-shear mixed crack, the crack opening degree needs to be corrected by shear dilation, and the crack opening degree of the crack unit is obtained in the following manner: The initial opening of the fracture is set as Δ n 0 = 0, the fracture opening Δ of the fracture element n is calculated by the following formula: where, S is the shear displacement, k 1 is a material parameter, σ n is the normal stress, σ U is the uniaxial compressive strength of the adjacent rock mass, ψ 0 is the initial dilatancy angle.

6. The apparatus of claim 5, wherein, The water inrush volume determination module is also used for: Based on the crack opening degree of each crack unit, in combination with an enhanced immersed boundary algorithm and Hugoniot conditions, identifying a fluid-solid contact boundary and determining the normal velocity, tangential velocity, pressure and density of fluid on both sides of the fluid-solid contact boundary; Based on a fluid volume VOF method, calculating the volume fraction of fluid in each Euler unit with gas-water two-way flow in the Euler grid component according to the normal velocity, tangential velocity, pressure and density of fluid; Multiplying the volume of each Euler unit with gas-water two-way flow by the volume fraction of fluid to obtain the water inrush volume in each Euler unit; Accumulating the water inrush volumes of all Euler units to obtain the mining fault water inrush volume.

7. A computer readable storage medium characterized in that, The computer readable storage medium stores a computer program, and the computer program is executed by the processor to implement the simulation method for predicting the mining fault water inrush volume according to any one of claims 1-4.

8. A computer program product, characterised in that, Computer programs / instructions, when executed by a processor, implement the simulation method for predicting water inrush from a mining-induced fault according to any one of claims 1-4.

Images