CN113378410A - Method for simulating evolution of mining overburden water guide channel - Google Patents

Method for simulating evolution of mining overburden water guide channel Download PDF

Info

Publication number
CN113378410A
CN113378410A CN202110781104.0A CN202110781104A CN113378410A CN 113378410 A CN113378410 A CN 113378410A CN 202110781104 A CN202110781104 A CN 202110781104A CN 113378410 A CN113378410 A CN 113378410A
Authority
CN
China
Prior art keywords
unit
fracture
rock
node
crack
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110781104.0A
Other languages
Chinese (zh)
Other versions
CN113378410B (en
Inventor
李�浩
武艳霞
武鹏飞
李文达
贺伟
韦婕
姚宏波
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Taiyuan University of Technology
Original Assignee
Taiyuan University of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Taiyuan University of Technology filed Critical Taiyuan University of Technology
Priority to CN202110781104.0A priority Critical patent/CN113378410B/en
Publication of CN113378410A publication Critical patent/CN113378410A/en
Application granted granted Critical
Publication of CN113378410B publication Critical patent/CN113378410B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30—Assessment of water resources

Abstract

The invention belongs to the technical field of water retention coal mining, and discloses a method for simulating the evolution of a water guide channel of a mining overburden, which comprises the following steps: s1, determining geometric parameters of a structural surface in overlying strata, and establishing a three-dimensional network geometric model of the structural surface in the overlying strata; the structural surface comprises a main structural surface and a secondary structural surface; s2, establishing a numerical calculation model of the evolution of the stope overburden water guide channel containing the structural surface; s3, determining the size of a zero-thickness fracture unit, and embedding the zero-thickness fracture unit at the overlying rock structure surface and the entity unit boundary in the matrix of the numerical calculation model respectively so as to simulate the random distribution of a main structure surface and a secondary structure surface in overlying rock; and S4, respectively obtaining the elastic-plastic characteristics of the rock stratum matrix, the brittle/tough fracture mechanical characteristics of the structural surface and the shearing friction characteristics of the completely fractured rock mass, endowing the characteristics to the structural surface and the fracture unit of the numerical calculation model, and adopting FEM and DEM coupling calculation.

Description

Method for simulating evolution of mining overburden water guide channel
Technical Field
The invention belongs to the technical field of water retention coal mining, and particularly relates to a method for simulating the evolution of a water guide channel of a mining overburden.
Background
The Jinshanmeng area contains nearly 60% of coal resources in China, but only has 8% of water resources in China. More seriously, the high-strength underground coal mining can cause the movement and the damage of the overlying rock stratum, thereby affecting the aquifer within the influence range of the water-flowing fractured zone, causing the leakage of underground water and the damage of the ecological environment. This brings great challenges to the goals of ecological protection and high-quality development of the yellow river basin. The engineering practice of Yanzishan mine in great homogeneous mining area shows that: coal-based formations comprise brittle sandstone, limestone, and tough mudstone and kaolinite, and the formations comprise a large number of main structural planes and secondary structural planes. The complex mechanical property and abundant structural surface of the stratum can have great influence on the distribution, height and evolution of the water guide channel of the mining overburden rock, and whether the water in the dwarfism goaf can be discharged to the secondary coal face is determined. A reasonable numerical simulation method is adopted to predict the development height of the water guide channel of the mining overburden rock, which is important for optimizing working surface parameters and preserving water and mining coal.
In recent years, scholars at home and abroad mainly adopt finite element FEM and discrete element DEM methods to carry out mining overburden rock motion simulation research. Wherein:
the FEM method is based on continuous medium hypothesis, and researches the development height and evolution rules of a force field, a plastic region, a damage zone and the like of overlying strata in the mining process by adopting elasticity, plasticity and damage mechanics theories.
The DEM method is based on the discrete medium hypothesis, most structural surfaces in overlying strata are set to be in a regular masonry structural form, elastic brittleness and elastic mechanical constitutive models are adopted for the overlying strata structural surfaces and matrixes respectively, and the development heights of horizontal cracks and vertical cracks in the overlying strata in the coal seam mining process are obtained through numerical simulation. However, the structural surfaces in the overlying strata are randomly distributed and have different scales, and a multi-level structural surface system and a random generation program thereof need to be established in the numerical calculation model; meanwhile, mudstone and kaolinite in the overburden rock have remarkable toughness destruction characteristics.
The DEM method can change the propagation path of the mining overburden rock water guide channel and has a decisive effect on the distribution and evolution regulation of the mining overburden rock water guide channel; the FEM method can consume energy generated by secondary stress to a great extent and influence the development height of the overlying strata fracture.
Disclosure of Invention
In view of the above, in order to solve the problems in the background art, the present invention provides a method for simulating the evolution of a water guide channel for mining overburden.
In order to achieve the purpose, the invention provides the following technical scheme: a simulation method for mining overburden water guide channel modeling comprises the following steps:
s1, determining geometric parameters of a structural surface in overlying strata, and establishing a three-dimensional network geometric model of the structural surface in the overlying strata; the structural surface comprises a main structural surface and a secondary structural surface;
s2, establishing a numerical calculation model of the evolution of the stope overburden water guide channel containing the structural surface;
s3, determining the size of a zero-thickness fracture unit, and embedding the zero-thickness fracture unit at the overlying rock structure surface and the entity unit boundary in the matrix of the numerical calculation model respectively so as to simulate the random distribution of a main structure surface and a secondary structure surface in overlying rock;
s4, respectively obtaining the elastic-plastic characteristics of a rock stratum matrix, the brittle/tough fracture mechanical characteristics of a structural plane and the shearing friction characteristics of a completely fractured rock mass, endowing each characteristic to the structural plane and a fracture unit of the numerical calculation model, and adopting FEM and DEM coupling calculation;
s5, acquiring and inputting material parameters, boundary conditions and initial conditions of the numerical calculation model;
s6, compiling a VUSFLD program for deleting units at the coal seam position along with time and space positions, and embedding the VUSFLD program into the numerical calculation model to simulate the stoping process of the coal face;
s7, solving a simulation result of the numerical calculation model and outputting the simulation result;
and S8, analyzing the simulation result to obtain an evolution rule of the space distribution and the development height of the mining overburden water flowing fractured zone along with the mining parameters, and combining the evolution rule with the distribution state of the surface water-bearing stratum and the underground water-bearing stratum to optimize a coal mining scheme.
Preferably, the determining the geometric parameters of the structural plane in the overlying strata in the step S1 includes:
determining the spacing of the main structural surfaces in each rock stratum according to a resistivity logging method;
determining the spacing of the secondary structural surfaces according to the average length of the rock cores of each rock layer in the drilling and coring process;
and calculating the average spacing, the mean square error of the spacing and the occurrence frequency of the structural surface within a preset inclination angle range.
Preferably, the step S1 of establishing a three-dimensional network geometric model of the structural plane in the overburden rock includes:
determining the number of Voronoi polygons in three orthogonal directions of x, y and z according to the size of a numerical calculation model required to be established;
generating a group of space dispersion points meeting the statistical distribution in python according to the space mean square error of the structural surface;
connecting all the spatial discrete points into a Delaunay triangle;
finding out a fourth point outside three points of the triangle according to a Delaunay space circumsphere criterion, and forming a Delaunay tetrahedron;
a Voronoi polyhedron is enclosed by the perpendicular planes of adjacent tetrahedron faces, and therefore a three-dimensional network geometric model of the structural face in the overlying rock is generated.
Preferably, the determining the size of the zero-thickness crack unit in the step S3 includes:
length of fracture process zone L:
Figure BDA0003157143100000031
wherein, pi is a circumferential rate constant, mu is a Poisson's ratio of the material, D0，cIs modulus of elasticity, GnIn order to be able to break the energy,
Figure BDA0003157143100000032
is the tensile strength of the material;
the average length T of each rock core in the drilling and coring process;
in the large-scale structural surface, the size of the crack unit is less than or equal to L;
in the small-scale structural plane, the size of the crack unit is the smaller of the average length T and the length L.
Preferably, the embedding the zero-thickness crack unit in the step S3 includes:
(1) updating entity unit node number:
performing global grid division on the numerical calculation model established in the step S2, generating an inp script file in a WRITE INPUT under the Job module, respectively finding the maximum node number and the maximum unit number of the entity unit in the keywords NSET and ELSET in the script file, and respectively recording the maximum node number and the maximum unit number as nmax，emax；
Numbering the physical units as ejN of (2)i(niNot equal to 1) a nodes are duplicated, and a is the number of shared physical units;
on the premise that the node coordinates are the same, the node number of the minimum entity unit number is kept unchanged, and the numbers of the rest a-1 nodes are modified as follows: n is a radical ofi＝(a-1)×10Ten(nmax)+ni(ii) a Wherein, Ten (n)max) Is nmaxDecimal digits of (d); specifically, after the above transformation, the n-thiEach node is duplicated a times with the number n1,n2,…,na(ii) a While sharing the n-thiThe unit number of each node is a, and the unit numbers are e in sequence1,e2,…,eaCopying the new node number at the position of the corresponding unit number to form two one-to-one ordered arrays Elem ═ (e)1,e2,…,eb) And Node ═ n1,n2,…,nb) Copying the updated ordered array into the inp script file to update the node number of the entity unit；
(2) Creating crack Unit node number, crack Unit number
In the initial inp script file, the nodes with the same entity unit number are copied and form a new array according to the sequence determined by the following formula:
Group＝(ej,ek,nShared01,nShared02,nShared03,nShared04,Node01,Node02,Node03,Node04) (ii) a Wherein the content of the first and second substances,
e in array GroupjAnd ekNumbering adjacent entity units; n isShared01-nShared04The same numbered nodes between adjacent physical units in the initial inp file; at Node01-Node04In, only remain nShared01-nShared04The updated entity unit node numbers are arranged in reverse order to obtain a unit ejAnd ekThe numbers of the crack unit nodes in between;
the number of the crack unit is determined according to the maximum number of the entity unit: e.g. of the typecoh,min＝k+ei(ii) a Wherein k is an integer power multiple of 10, and k is greater than the maximum number of the entity unit; for example, the maximum node number of the entity unit is 756, the k value can be 1000, and can also be 10000;
combining the crack unit node number with the crack unit number to form a crack unit number set; specifically, a key word ". times.element type ═ COH2D 4" (taking a quadrangular unit as an example, and if the hexahedral unit is adopted, the Element type is modified to "times.element type ═ COH3D 8") is added to the inp file, and ". times.element ═ cohlell" is added between the key word ". times.part" and the "assemby", and then corresponding numbers are written in the order of the first unit number and the second node number, so that all the fracture singlets are combined into a set named "cohlell";
(3) embedding crack cells at a principal structural plane of the numerical calculation model
Duplicating the updated target crack cell number ejAnd corresponding node coordinates nj(xj，yj，zj) And calculating the direction vector d of the target crack unit by inner product methodcoh＝(d1,d2,d3)；
Judgment vector dcohWhether the direction vector d of the main structural plane in the numerical calculation model is related toipCollinear, judging vector dcohIf any node in the target crack unit is located on the main structural plane, determining that the target crack unit is located on the main structural plane.
Specifically, the target crack cell number and the node number satisfying the above conditions are copied in the inp file, and "Elset ═ cohElemDIS" is added between the keywords "· Part" and ". assembiy", so that the cells on all the structural planes are grouped into one set called "cohElemDIS".
Preferably, the step S4 of obtaining the elasto-plastic characteristics of the formation matrix includes:
determining elasto-plastic deformation characteristics of formation matrix during mining
Figure BDA0003157143100000054
Figure BDA0003157143100000051
Wherein:
σsis a stress vector,
Figure BDA0003157143100000052
Is plastic strain, E0Is a matrix of elastic stiffness, epsilonsIs the total strain;
A＝(σb0/σc0-1)/(2σb0/σc0-1)；B＝σs/σt·(1-A)-(1+A)；C＝3(1-Kc)/(2Kc-1)；
peff＝trace(σs,eff)/3,
Figure BDA0003157143100000053
Seff＝σ-peffI；
Figure BDA0003157143100000061
(i ═ 1,2,3) is the maximum principal stress in the three orthogonal directions x, y, z;
σb0/σc0is the ratio of the yield stress under triaxial and uniaxial compression conditions;
σtand σsIs the equivalent tensile and compressive stress; kcThe value is 0.667.
Preferably, the obtaining of the brittle/ductile fracture mechanical characteristics of the structural surface and the shear friction characteristics of the completely fractured rock mass in the step S4 includes:
(1) determining elastic deformation characteristics of non-through structural surface appearance in overburden
σc＝D0,cεc；εc＝S/T0；
Wherein D0，cIs a matrix of elastic stiffness, epsiloncIs a strain vector, T0Is the intrinsic thickness of the slit cell, and T0The value is 0.01 times of the characteristic size of the entity unit, and S is the displacement vector of the crack unit;
(2) determining brittle fracture and ductile fracture characteristics of non-through structural planes in overburden rock
Brittle fracture characteristics:
Figure BDA0003157143100000062
εc＝S/T0(ii) a Wherein σc,lAnd
Figure BDA0003157143100000063
traction and traction peak, respectively, with l being expressed in the normal and two tangential directions;
ductile fracture characteristics:
Figure BDA0003157143100000064
wherein S isn，SsAnd StAre all variables and represent the displacement in the normal and two tangential directions, respectively; gamma-shapedn,ΓSIs the fracture energy constant and is expressed as:
Figure BDA0003157143100000065
GS＝Gs+Gt；
Figure BDA0003157143100000066
χn＝sn,p/sn,χS＝sS,p/sS；
in the formula, GnAnd GSNormal and total tangential fracture energies, respectively; gsAnd GtThe fracture energy in two tangential directions respectively; beta and gamma are stress-strain curve shape parameters of pure I type and pure II type fracture experiments respectively; sn,pAnd sS,pRespectively the peak displacement of the stress-strain curve of the pure I type fracture experiment and the pure II type fracture experiment; snAnd sSThe fracture displacement of the stress-strain curves of the pure I-type fracture experiment and the pure II-type fracture experiment are respectively obtained;
(3) determining shear friction characteristics of a fully fractured rock mass
σs＝μ1σnn(ii) a Wherein σsIs a shear stress; mu.s1Is the coefficient of friction; sigmannIs the normal stress.
Preferably, the calculation of coupling the FEM and the DEM in step S4 includes:
determining boundary conditions, and solving the node force of entity units in the overburden rock matrix in the numerical calculation model through FEM;
transmitting force field data through a shared node of the entity unit and the crack unit, and calculating a displacement component of the crack unit by adopting an explicit solution;
at s > sFJudging that the crack unit is broken and responding through the DEM; wherein
Figure BDA0003157143100000071
sFIs the fracture displacement in the mixed fracture mode;
in DEM, solving by a Newton-Raphson iterative method, wherein the expression of iterative solution is as follows: jacobi (u)j)duj+1＝-Rj(ii) a Wherein, Jacobi (u)j) Is the Jacobian matrix in the jth iteration step, RjResidual errors in the j iteration step are obtained; the iteration convergence condition is as follows:
||R||2< tolerance and | | u | | ventilation2< tolerance; wherein | | | purple hair2Is a norm of two types of vectors, and the tolerance is 10-8Then, the iteration is converged after 2-5 times.
Preferably, the process of simulating coal face extraction in step S6 includes:
defining a state variable parameter Depvar as 1 in a Materials module;
setting a state parameter STATUS containing a completely broken crack unit in an output variable in a Step module;
inputting total time totalTime, analysis step time stepTime, unit number nblock and unit integration point coordinates coordMp in a VUSDFLD program, and when a crack unit in a coal seam meets a set condition, designating stateNew (nblock,1) as 0 to delete the completely broken crack unit;
the setting conditions are as follows:
f(totalTime,coordMp)＝stepTime×v+coordMp(nblock,1)
wherein v represents the coal face advancing speed; "1" in coordMp (nblock,1) indicates the coal seam excavation direction in the numerical calculation model.
Compared with the prior art, the invention has the following beneficial effects:
(1) the invention provides an effective simulation tool for overburden rock breaking, distribution and evolution of the water flowing fractured zone caused by underground coal seam mining. Specifically, by adopting the FEM + DEM method provided by the invention, crack units can be embedded in all the main and secondary structural surfaces of the overlying strata, and the cracks can be expanded in the overlying strata at will in the mining process. In addition, the FEM + DEM method is a great improvement on the traditional finite element, discrete element and Lagrange element method, combines the elastoplasticity characteristics of the overlying strata matrix and reflects the brittleness and toughness fracture mechanical characteristics of the overlying strata structural surface and microcracks, and realizes the numerical simulation of the whole process of complete → fracture → shear slip along the fracture surface of the overlying strata in the mining process from the calculation method and the mechanical theory; in conclusion, the distribution and evolution of the overlying strata collapse zone and the fracture conduction zone in the mining process are effectively obtained, and a simulation method is provided for optimization of mining parameters.
(2) The method adopts an implicit time integral solving algorithm aiming at the entity unit and an explicit time integral solving algorithm aiming at the crack unit, so that the convergence of the numerical solution can be greatly improved, and the calculation efficiency is accelerated.
(3) The method of the invention is based on the logging information and the length of the core, and respectively establishes a main structural surface model and a secondary structural surface model, and crack units with different material parameters are embedded at two positions, so that the shearing slippage of a large structural surface and the random fragmentation process of a complete rock body between the structural surfaces are effectively simulated, and the random expansion and the random distribution of macro and micro cracks in the coal seam overlying strata are realized in the mining process.
(4) The method can be used for analyzing the distribution and evolution law of the water flowing fractures of the mining overburden rock under the influence of different mining parameters and different geological parameters, so that the mining parameters are optimized, and water-retaining coal mining is realized.
Drawings
FIG. 1 is a flow chart of a simulation method of the present invention;
FIG. 2 is a schematic diagram of a three-dimensional network geometry model of a structural plane in overburden rock according to the present invention;
FIG. 3 is a schematic diagram of a method for globally embedding a crack unit according to the present invention;
FIG. 4a is a schematic view of a crack unit embedded at a major structural plane of the present invention;
FIG. 4b is a schematic view of the embedded crack units at the major structural plane of the present invention;
FIG. 4c is a schematic representation of a solid element within a matrix according to the present invention;
FIG. 5 is a diagram showing a numerical calculation result of mining overburden rock fracture distribution;
FIG. 6 is a schematic diagram of distribution positions of main water flowing fractured zones;
FIGS. 7 and 8 are stress cloud charts of overburden Mises and corresponding mining fractures;
Detailed Description
The present invention will be further described with reference to the accompanying drawings, but the scope of the present invention is not limited to the following.
S1, using Yanzishan mine in a great homogeneous mining area as a background;
determining the spacing of the main structural surfaces in each rock stratum according to a resistivity logging method;
determining the spacing of the secondary structural surfaces according to the average length of the rock cores of each rock layer in the drilling and coring process;
calculating the average spacing and the mean square error of the spacing of the structural surfaces and the occurrence frequency of the structural surfaces in a preset inclination angle range;
determining the number of Voronoi polygons in three orthogonal directions of x, y and z according to the size of a numerical calculation model required to be established; specifically, combining the geological conditions of Yanzishan mining engineering, wherein x is 870m, high y is 451m, and z is 1 m;
generating a group of space dispersion points meeting the statistical distribution in python according to the space mean square error of the structural surface;
connecting all the spatial discrete points into a Delaunay triangle;
finding out a fourth point outside three points of the triangle according to a Delaunay space circumsphere criterion, and forming a Delaunay tetrahedron;
the Voronoi polyhedrons are surrounded by the midperpendicular planes of the adjacent tetrahedrons, so that the three-dimensional network geometrical model of the structural surface in the overlying strata is generated as shown in fig. 2.
And S2, establishing a numerical calculation model of the evolution of the stope overburden water guide channel containing the structural plane.
S3, determining the size of the zero-thickness crack unit, and embedding the zero-thickness crack unit at the position of the overlying rock structural surface of the numerical calculation model and the boundary of the solid unit in the matrix respectively so as to simulate the random distribution of the main structural surface and the secondary structural surface in the overlying rock; wherein
When determining the size of the zero thickness fracture cell:
length of fracture process zone L:
Figure BDA0003157143100000101
wherein, pi is a circumferential rate constant, mu is a Poisson's ratio of the material, D0，cIs modulus of elasticity, GnIn order to be able to break the energy,
Figure BDA0003157143100000102
is the tensile strength of the material; wherein the elastic modulus, the fracture energy and the tensile strength are all obtained by a three-point bending experiment;
the average length T of each rock core in the drilling and coring process;
in the large-scale structural plane, the size of the crack unit is less than or equal to L; in the small-scale structural plane, the size of the crack unit is the smaller of the average length T and the length L.
When embedding zero thickness crack units:
adopting a hexahedron grid logarithmic value model to integrally divide grids under a Mesh module; then, by the method of globally embedding crack units as shown in fig. 3, three-dimensional 8-node crack units (as shown in fig. 4a and 4b) are respectively embedded at the primary structural surface (named cohElemDIS) and the secondary structural surface (named cohElemALL) of the numerical model, and the rest is 8-node hexahedral solid units (named baseElem as shown in fig. 4 c).
S4, decomposing the unit sets baseElem, cohemAll and cohemDIS into a baseElem-SY (sandstone entity unit set), a cohemAll-SY (sandstone secondary structure surface fracture unit set), a cohemDIS-SY (sandstone main structure surface fracture unit set) and corresponding fracture units and entity unit sets of mudstone (NY), coal (M), kaolinite (GL) and the like through a subtraction and selection method according to the actual distribution condition of the rock stratum; the material parameters of the matrix of the rock stratum (table 1), the main structural surface (table 2) and the secondary structural surface (table 3) are respectively input, and corresponding solid units and fracture units are endowed.
S5, constraining normal displacements of 5 boundaries of the numerical calculation model except the upper surface under an Assembly module, wherein the upper surface is a free boundary; then, according to the advancing speed and the model size of the coal face, 70 analysis steps (representing 70 excavation steps) are set under a Step module, each analysis Step is set to be geometric non-linear, the time length is 1 (representing 1d), and the time increment is set to be 0.001; setting output variables as stress, strain, frictional displacement, frictional stress, state variables STATUS and the like; the initial stress of the model is determined by means of self-balancing.
S6, defining the state Variable parameter Depvar 1 in Variable number controlling element deletion option of Materials module. The total time totalTime, the analysis step time stepTime, the unit number nblock and the unit integration point coordinate coordMp which are obtained by the calculation of the main program are transmitted into the VUSDFLD subprogram.
Setting key sentences in the subprogram as
DO k＝1,nblock
If(f(totalTime,coordMp).EQ.(stepTime×10m+coordMp(nblock,1)))
stateNew(nblock,1)＝0
END if
END DO
Therefore, the unit deletion at the appointed position under any analysis step can be realized, and the aim of advancing the coal face at the speed of 10m/d is fulfilled.
S7, a task is newly established under a Job module, and a programmed VUSDFLD subprogram is input into a User subcounting file option card; and setting the number of parallel computing cores to be 48 cores in an Edit Job tab, submitting the computation, solving the simulation result of the numerical computing model, and outputting the simulation result.
And S8, analyzing the simulation result to obtain an evolution rule of the space distribution and the development height of the mining overburden water flowing fractured zone along with the extraction parameters, and combining the evolution rule with the distribution state of the surface water-bearing stratum and the underground water-bearing stratum to optimize the coal mining scheme.
Specifically, the step of calculating the output simulation result in step S7 includes:
the mining overburden fracture distribution and evolution law as shown in figure 5; as can be seen from the calculation results and the graphs, the mining of the upper four layers of coal (the dwarfism coal seam) in the model can cause the complete fragmentation of the overlying rock layer, and the water diversion fracture zone passes through the aquifer and directly reaches the earth surface; after the lowest coal seam (the carboniferous coal seam) of the model is mined, corresponding overlying strata cracks are generated, the height of a caving zone is about 26m, and a water-guiding crack zone extends to a Jurassic goaf, so that accumulated water in the goaf can be drained. But the mining fracture of the kaolinite rock stratum area is not developed very much due to the obvious plastic deformation of the kaolinite; meanwhile, the permeability of the kaolinite is low, and the fracture opening is reduced after water absorption and expansion, so that accumulated water in the goaf is judged to be discharged, but the water inflow is limited.
The location of the distribution of the primary water fractured zones as shown in FIG. 6; according to the figure, overlying strata which collapse in the goaf range are compacted tightly by the upper strata, so that the permeability is low; and a water guiding fractured zone with larger width is often formed at the edge of the goaf, so that the position is a water control key area.
Overburden Mises stress cloud charts and corresponding mining fracture distributions as shown in fig. 7 and 8; specifically, fig. 7 is a schematic view of a stress cloud map of overburden Mises and a corresponding mining-induced fracture distribution formed by the 52 th iterative analysis, and fig. 8 is a schematic view of a stress cloud map of overburden Mises and a corresponding mining-induced fracture distribution formed by the 63 th iterative analysis. It can be seen that the proximity of the coal face to the main structural face of the immediate roof results in a rapid increase in Mises stress in the roof; when the coal face, the main structural surface of the direct roof and the main structural surface of the basic roof are close to each other, the stress value of Mises reaches the maximum and the range is the widest, so that the height of the overburden water flowing fracture is increased.
In conclusion, compared with the existing overburden rock water flowing fracture model, on one hand, the influence of the distribution of the primary and secondary structural surfaces on the distribution evolution of the fracture guiding zone is considered, and the invention better accords with the characteristic of discontinuous rock geometry; on the other hand, the elastoplasticity mechanical characteristics of the overburden matrix, the brittle/ductile fracture of the overburden structural surface and the friction effect after complete fracture are considered, and the process of converting the overburden from a continuous medium to a discrete medium is successfully simulated by combining the FEM-DEM numerical simulation method, so that a powerful tool is provided for simulating the evolution of mining overburden water flowing fractures.
Although embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that changes, modifications, substitutions and alterations can be made in these embodiments without departing from the principles and spirit of the invention, the scope of which is defined in the appended claims and their equivalents.
Table 1: mechanical parameters of overlying strata matrix
E/GPa μ σb0/σc0 as /MPa σt /MPa Kc δ ψ/°
Soil for soil 0.02 0.31 2.2 2.5 0.12 0.667 0.1 30
Siltstone 10.7 0.21 1.5 62 3.0 0.667 0.1 40
Fine sandstone 13.8 0.22 1.5 85 4.8 0.667 0.1 42
Middle sandstone 20.6 0.22 1.4 107 5.7 0.667 0.1 42
Kaolin rock 7.4 0.30 1.8 35 2.6 0.667 0.1 30
Coal (coal) 1.6 0.35 1.9 14 1.1 0.667 0.1 34
Mudstone 8.0 0.24 1.6 43 4.6 0.667 0.1 35
Table 2: mechanical parameters of main structural surface
Figure BDA0003157143100000141
Table 3: mechanical parameters of secondary structural surface
Figure BDA0003157143100000151

Claims (10)

1. A method for simulating the evolution of a water guide channel of a mining overburden rock is characterized by comprising the following steps:
s1, determining geometric parameters of a structural surface in overlying strata, and establishing a three-dimensional network geometric model of the structural surface in the overlying strata; the structural surface comprises a main structural surface and a secondary structural surface;
s2, establishing a numerical calculation model of the evolution of the stope overburden water guide channel containing the structural surface;
s3, determining the size of a zero-thickness fracture unit, and embedding the zero-thickness fracture unit at the overlying rock structure surface and the entity unit boundary in the matrix of the numerical calculation model respectively so as to simulate the random distribution of a main structure surface and a secondary structure surface in overlying rock;
s4, respectively obtaining the elastic-plastic characteristics of a rock stratum matrix, the brittle/tough fracture mechanical characteristics of a structural plane and the shearing friction characteristics of a completely fractured rock mass, endowing each characteristic to the structural plane and a fracture unit of the numerical calculation model, and adopting FEM and DEM coupling calculation;
s5, acquiring and inputting material parameters, boundary conditions and initial conditions of the numerical calculation model;
s6, compiling a VUSFLD program for deleting units at the coal seam position along with time and space positions, and embedding the VUSFLD program into the numerical calculation model to simulate a coal face stoping process;
and S7, solving the simulation result of the numerical calculation model and outputting the simulation result.
2. The method for simulating the evolution of the mining overburden water conducting channel according to claim 1, further comprising:
and S8, analyzing the simulation result to obtain an evolution rule of the space distribution and the development height of the mining overburden water flowing fractured zone along with the mining parameters, and combining the evolution rule with the distribution state of the surface water-bearing stratum and the underground water-bearing stratum to optimize a coal mining scheme.
3. The method for simulating the evolution of the water channel for mining overburden as claimed in claim 1 or 2, wherein the step S1 of determining the geometrical parameters of the structural plane in overburden comprises:
determining the spacing of the main structural surfaces in each rock stratum according to a resistivity logging method;
determining the spacing of the secondary structural surfaces according to the average length of the rock cores of each rock layer in the drilling and coring process;
and calculating the average spacing, the mean square error of the spacing and the occurrence frequency of the structural surface within a preset inclination angle range.
4. The method for simulating the evolution of the water guide channel of the mining overburden as claimed in claim 3, wherein the step S1 of establishing the three-dimensional network geometric model of the structural plane in the overburden comprises the following steps:
determining the number of Voronoi polygons in three orthogonal directions of x, y and z according to the size of a numerical calculation model required to be established;
generating a group of space discrete points meeting the statistical distribution in python according to the space mean square error of the structural surface;
connecting all the spatial discrete points into a Delaunay triangle;
finding out a fourth point outside three points of the triangle according to a Delaunay space circumsphere criterion, and forming a Delaunay tetrahedron;
a Voronoi polyhedron is enclosed by the perpendicular planes of adjacent tetrahedron faces, and therefore a three-dimensional network geometric model of the structural face in the overlying rock is generated.
5. The method for simulating the evolution of the mining overburden water conducting channel according to claim 4, wherein the step S3 of determining the size of the zero-thickness fracture unit comprises:
length of fracture process zone L:
Figure FDA0003157143090000021
wherein, pi is a circumferential rate constant, mu is a Poisson's ratio of the material, D0，cIs modulus of elasticity, GnIn order to be able to break the energy,
Figure FDA0003157143090000022
is the tensile strength of the material;
the average length T of each rock core in the drilling and coring process;
in the large-scale structural surface, the size of the crack unit is less than or equal to L;
in the small-scale structural plane, the size of the crack unit is the smaller of the average length T and the length L.
6. The method for simulating the evolution of the mining overburden water conducting channel according to claim 5, wherein the step S3 of embedding the zero-thickness fracture unit comprises:
(1) updating entity unit node number:
performing global grid division on the numerical calculation model established in the step S2, generating an inp script file in a WRITE INPUT under the Job module, respectively finding the maximum node number and the maximum unit number of the entity unit in the keywords NSET and ELSET in the script file, and respectively recording the maximum node number and the maximum unit number as nmax，emax；
Numbering the physical units as ejN of (2)i(niNot equal to 1) a nodes are copied, and a is the number of the shared entity units;
on the premise that the node coordinates are the same, the node number of the minimum entity unit number is kept unchanged, and the numbers of the rest a-1 nodes are modified as follows:
Figure FDA0003157143090000031
wherein, Ten (n)max) Is nmaxDecimal digits of (d);
(2) creating crack Unit node number, crack Unit number
In the initial inp script file, the nodes with the same entity unit number are copied and form a new array according to the sequence determined by the following formula:
Group＝(ej,ek,nShared01,nShared02,nShared03,nShared04,Node01,Node02,Node03,Node04) (ii) a Wherein e in the array GroupjAnd ekNumbering adjacent entity units; n isShared01-nShared04The nodes with the same number between adjacent entity units in the initial inp file; at Node01-Node04In, only remain nShared01-nShared04The updated entity unit node numbers are arranged in reverse order to obtain a unit ejAnd ekThe numbers of the crack unit nodes in between;
the number of the crack unit is determined according to the maximum number of the entity unit: e.g. of the typecoh,min＝k+ei(ii) a Wherein k is an integer power multiple of 10, and k is greater than the maximum number of the entity unit;
combining the crack unit node number with the crack unit number to form a crack unit number set;
(3) embedding crack cells at a principal structural plane of the numerical calculation model
Duplicating the updated target crack cell number ejAnd corresponding node coordinates nj(xj，yj，zj) And calculating the direction vector d of the target crack unit by inner product methodcoh＝(d1,d2,d3)；
Judgment vector dcohWhether the direction vector d of the main structural plane in the numerical calculation model is related toipCollinear, judging vector dcohIf any node in the target crack unit is located on the main structural plane, determining that the target crack unit is located on the main structural plane.
7. The method for simulating the evolution of the mining overburden water conduit according to claim 6, wherein the step S4 of obtaining the elastoplasticity characteristics of the rock stratum matrix comprises the following steps:
determining elasto-plastic deformation characteristics of formation matrix during mining
Figure FDA0003157143090000041
Figure FDA0003157143090000042
Wherein:
σsis a stress vector,
Figure FDA0003157143090000043
Is plastic strain, E0Is a matrix of elastic stiffness, epsilonsIs the total strain;
A＝(σb0/σc0-1)/(2σb0/σc0-1)；B＝σs/σt·(1-A)-(1+A)；C＝3(1-Kc)/(2Kc-1)；
peff＝trace(σs)/3,
Figure FDA0003157143090000044
S＝σs-pI；
Figure FDA0003157143090000045
is the maximum principal stress in the three orthogonal directions x, y, z;
σb0/σc0is the ratio of the yield stress under triaxial and uniaxial compression conditions;
σtand σsIs the equivalent tensile and compressive stress; kcThe value is 0.667.
8. The method for simulating the evolution of the mining overburden water guide channel as claimed in claim 7, wherein the step S4 of obtaining the brittle/ductile fracture mechanical characteristics of the structural surface and the shear friction characteristics of the completely fractured rock mass comprises:
(1) determining elastic deformation characteristics of non-through structural surface appearance in overburden
σc＝D0,cεc；εc＝S/T0；
Wherein D0，cIs a matrix of elastic stiffness, epsiloncIs a strain vector, T0Is the intrinsic thickness of the slit cell, and T0The value is 0.01 times of the characteristic size of the entity unit; and S is the displacement vector of the crack unit.
(2) Determining brittle fracture and ductile fracture characteristics of non-through structural planes in overburden rock
Brittle fracture characteristics:
Figure FDA0003157143090000051
wherein σc,lAnd
Figure FDA0003157143090000052
traction and traction peak, respectively, with l being expressed in the normal and two tangential directions;
ductile fracture characteristics:
Figure FDA0003157143090000053
wherein S isn，SsAnd StAre all variables and represent the displacement in the normal and two tangential directions, respectively; gamma-shapedn,ΓSIs the fracture energy constant and is expressed as:
Figure FDA0003157143090000054
GS＝Gs+Gt；
Figure FDA0003157143090000055
χn＝sn,p/sn,χS＝sS,p/sS；
in the formula, GnAnd GSNormal and total tangential fracture energies, respectively; gsAnd GtThe fracture energy in two tangential directions respectively; beta and gamma are respectively the stress-strain curve shape parameters of the pure I type fracture experiment and the pure II type fracture experiment; sn,pAnd sS,pRespectively the peak displacement of the stress-strain curve of the pure I type fracture experiment and the pure II type fracture experiment; snAnd sSThe fracture displacement of the stress-strain curves of the pure I-type fracture experiment and the pure II-type fracture experiment are respectively obtained;
(3) determining shear friction characteristics of a fully fractured rock mass
σs＝μ1σnn(ii) a Wherein σsIs a shear stress; mu.s1Is the coefficient of friction; sigmannIs the normal stress.
9. The method for simulating the evolution of the mining overburden water channel as claimed in claim 8, wherein the calculation of coupling the FEM and the DEM in step S4 comprises:
determining boundary conditions, and solving the node force of the entity unit in the overburden rock matrix in the numerical calculation model through FEM;
transmitting force field data through a shared node of the entity unit and the crack unit, and calculating a displacement component of the crack unit by adopting an explicit solution;
at s > sFJudging that the crack unit is broken and responding through the DEM; wherein
Figure FDA0003157143090000061
sFFracture displacement in a hybrid fracture mode;
in DEM, solving by a Newton-Raphson iterative method, wherein the expression of iterative solution is as follows: jacobi (u)j)duj+1＝-Rj(ii) a Wherein, Jacobi (u)j) Is the Jacobian matrix in the jth iteration step, RjResidual errors in the j iteration step are obtained; iterative harvestingThe converging conditions are as follows:
||R||2< tolerance and | | u | | ventilation2< tolerance; wherein | | | purple hair2Is a norm of two types of vectors, and the tolerance is 10-8Then, the iteration is converged after 2-5 times.
10. The method for simulating the evolution of the mining overburden water conducting channel according to claim 9, wherein the method comprises the following steps: the process of the step S6 for the recovery of the simulated coal face includes:
defining a state variable parameter Depvar as 1 in a Materials module;
setting a state parameter STATUS containing a completely broken crack unit in an output variable in a Step module;
inputting total time totalTime, analysis step time stepTime, unit number nblock and unit integration point coordinates coordMp in a VUSDFLD program, and when a crack unit in a coal seam meets a set condition, designating stateNew (nblock,1) as 0 to delete the completely broken crack unit;
the setting conditions are as follows:
f(totalTime,coordMp)＝stepTime×v+coordMp(nblock,1)
wherein v represents the coal face advancing speed; "1" in coordMp (nblock,1) represents the coal seam excavation direction in the numerical calculation model.
CN202110781104.0A 2021-07-10 2021-07-10 Simulation method for evolution of mining overburden water guide channel Active CN113378410B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110781104.0A CN113378410B (en) 2021-07-10 2021-07-10 Simulation method for evolution of mining overburden water guide channel

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110781104.0A CN113378410B (en) 2021-07-10 2021-07-10 Simulation method for evolution of mining overburden water guide channel

Publications (2)

Publication Number Publication Date
CN113378410A true CN113378410A (en) 2021-09-10
CN113378410B CN113378410B (en) 2023-06-16

Family

ID=77581676

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110781104.0A Active CN113378410B (en) 2021-07-10 2021-07-10 Simulation method for evolution of mining overburden water guide channel

Country Status (1)

Country Link
CN (1) CN113378410B (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116244794A (en) * 2022-12-31 2023-06-09 中铁隧道勘察设计研究院有限公司 Calculation method for minimum safe overlying strata thickness of underwater non-blasting undercut tunnel

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2006323569A (en) * 2005-05-18 2006-11-30 Foundation Of River & Basin Integrated Communications Japan Disaster prevention/evacuation action simulation system
CN204731234U (en) * 2015-07-10 2015-10-28 贵州理工学院 Close-in seams water protection mining solid-liquid coupling analog simulation device
CN107044289A (en) * 2017-06-22 2017-08-15 中国矿业大学 A kind of water damage prevention and controls of bored grouting closure overlying strata water producing fractures main channel
CN108763650A (en) * 2018-04-28 2018-11-06 湘潭大学 A kind of overlying strata mining induced fissure network model construction method

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2006323569A (en) * 2005-05-18 2006-11-30 Foundation Of River & Basin Integrated Communications Japan Disaster prevention/evacuation action simulation system
CN204731234U (en) * 2015-07-10 2015-10-28 贵州理工学院 Close-in seams water protection mining solid-liquid coupling analog simulation device
CN107044289A (en) * 2017-06-22 2017-08-15 中国矿业大学 A kind of water damage prevention and controls of bored grouting closure overlying strata water producing fractures main channel
CN108763650A (en) * 2018-04-28 2018-11-06 湘潭大学 A kind of overlying strata mining induced fissure network model construction method

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116244794A (en) * 2022-12-31 2023-06-09 中铁隧道勘察设计研究院有限公司 Calculation method for minimum safe overlying strata thickness of underwater non-blasting undercut tunnel
CN116244794B (en) * 2022-12-31 2023-12-29 中铁隧道勘察设计研究院有限公司 Calculation method for minimum safe overlying strata thickness of underwater non-blasting undercut tunnel

Also Published As

Publication number Publication date
CN113378410B (en) 2023-06-16

Similar Documents

Publication Publication Date Title
Lu et al. Numerical simulation of mining-induced fracture evolution and water flow in coal seam floor above a confined aquifer
Li et al. A numerical investigation of the hydraulic fracturing behaviour of conglomerate in Glutenite formation
Ye et al. Numerical simulation on tendency mining fracture evolution characteristics of overlying strata and coal seams above working face with large inclination angle and mining depth
Li et al. A possible prediction method to determine the top concealed karst cave based on displacement monitoring during tunnel construction
Deng et al. 3D finite element modeling of directional hydraulic fracturing based on deformation reinforcement theory
CN113378410B (en) Simulation method for evolution of mining overburden water guide channel
Deng et al. Numerical simulation on the evolution of mining-induced fracture network in a coal seam and its overburden under the top coal caving method
Xia et al. Mechanism of excavation-induced cracking of the protective layer of a rock bench in a large underground powerhouse under high tectonic stress
Li et al. Influence of preexisting discontinuities on hydraulic fracture complexity in a naturally fractured reservoir
Yang et al. Prediction of residual gas content during coal roadway tunneling based on drilling cuttings indices and BA-ELM algorithm
Wang et al. Dynamic propagation behaviors of hydraulic fracture networks considering hydro-mechanical coupling effects in tight oil and gas reservoirs: A multi-thread parallel computation method
Yuan et al. Numerical modeling on hydraulic fracturing in coal-rock mass for enhancing gas drainage
Liu et al. Mining method selection and optimization for hanging-wall ore-body at Yanqianshan Iron Mine, China
Ye et al. A multi-objective optimization approach for a fault geothermal system based on response surface method
Niu et al. An intelligent prediction method of the karst curtain grouting volume based on support vector machine
CN112883661B (en) Fracturing simulation method of crushed soft low-permeability hydrocarbon reservoir
Lu et al. Prediction of fracture evolution and groundwater inrush from karst collapse pillars in coal seam floors: a micromechanics-based stress-seepage-damage coupled modeling approach
Mohajerani et al. An efficient computational model for simulating stress-dependent flow in three-dimensional discrete fracture networks
Zhou et al. Feasibility study on fully mechanized large mining height long wall top-coal caving mining in ultra-thick (20–30 m), parting-rich coal seams: A case study of the Laosangou mining field in China
Liu et al. Numerical modeling on anisotropy of seepage and stress fields of stratified rock slope
Nie et al. Coupled hydro-mechanical modelling on hydraulic fracturing using DDA
Zhang et al. Optimum Layout of Multiple Tree-type Boreholes in Low-Permeability Coal Seams to Improve Methane Drainage Performance
Sakhno et al. Numerical simulation of the surface subsidence evolution caused by the flooding of the longwall goaf during excavation of thin coal seams
Nan et al. A predictive model of mining collapse extent and its application
Xia et al. Network Design Mode of In-Seam Gas Extraction Parameters Using Mathematical Modelling—Take Tangan Colliery as an Example

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant