Continuous-discrete coupling simulation method for floor failure depth of fully mechanized caving mining working face in extra-thick coal seams

Abstract

The present invention discloses a continuous-discrete coupling simulation method for the failure depth of the floor of a fully mechanized caving working face in an extra-thick coal seam. The method comprises the following steps: determining the spacing between natural cracks in the floor and generating a geometric model containing the natural cracks; dividing the geometric model into solid units and embedding crack units at the locations of the natural cracks; experimentally determining the mechanical parameters of rock compression, fracture, and shear friction; compiling a VUSDFLD subroutine to simulate the working face advancement and coal mining process; importing the subroutine into a continuous-discrete coupling simulation main program, setting the numerical model boundary conditions, incremental steps, and output variables, and solving for the spatiotemporal evolution of the mining-induced floor cracks. The continuous-discrete coupling simulation method of the present invention can effectively simulate the static-dynamic loading / unloading conditions experienced by the floor of a fully mechanized caving working face in an extra-thick coal seam, as well as the entire process of the floor transforming from a continuous medium to a discrete medium, providing an effective simulation method for optimizing green coal mining parameters.

Description

Technical Field

[0001] The present invention relates to the field of floor water inrush mechanism and prevention, and green coal mining technology, and in particular to a continuous-discrete coupling simulation method for floor damage depth of a fully mechanized caving mining working face in an extra-thick coal seam. Background Art

[0002] Coal measures are composed of intact blocks of varying lithology and numerous fractures. The abundance of natural fractures within the floor, as well as those caused by mining, is a key factor contributing to floor water inrush and water resource damage. Coal mining practice demonstrates that the depth of floor damage in fully-mechanized caving (Jiangsu) mining faces in extra-thick coal seams is often greater than in medium-thick coal seams. Using appropriate numerical simulation methods to predict mining floor damage is crucial for optimizing parameters for fully-mechanized caving (Jiangsu) mining faces in extra-thick coal seams and promoting environmentally friendly mining.

[0003] In recent years, domestic and foreign scholars have mainly used numerical programs such as finite element and Lagrangian methods to study the depth of mining floor failure. Both finite element and Lagrangian methods are based on the continuum medium assumption and use elasticity, plasticity, and damage mechanics theories to study the stress, deformation, plastic zone depth, and evolution of the floor during mining. However, the rock of the mining floor often cracks rather than plastically deforms; the large number of natural cracks developed in the floor are randomly distributed and of varying scales, and it is necessary to generate randomly distributed natural cracks in numerical modeling; more importantly, under the conditions of fully mechanized caving mining in extremely thick coal seams, the floor rock undergoes quasi-static compression, quasi-static unloading, and dynamic compression in sequence. Only by effectively simulating this special load can the depth of crack development in the mining floor be numerically predicted. Summary of the Invention

[0004] The object of the present invention is to provide a continuous-discrete coupling simulation method for the floor failure depth of a fully mechanized caving mining working face in an extra-thick coal seam.

[0005] In order to achieve the above tasks, the present invention adopts the following technical solutions:

[0006] A continuous-discrete coupling simulation method for floor failure depth of a fully-mechanized caving mining working face in an extra-thick coal seam is characterized by comprising the following steps:

[0007] S1: Determine the spacing of natural cracks in the basement and generate a geometric model including natural cracks;

[0008] S2: Divide the geometric model into units and embed fracture units at the locations of natural fractures;

[0009] S3: Experimentally determine rock compression parameters, type I and type II fracture parameters, and rock shear friction parameters;

[0010] S4: Compile the VUSDFLD subroutine to automatically delete coal aggregate units and simulate the working face advancement process.

[0011] S5: Apply numerical model boundary conditions, input material mechanical parameters, set incremental step size, and set output parameters, including stress, strain, displacement, unit node force, and state variables of crack units.

[0012] S6: The model was solved by continuous-discrete coupling simulation program to obtain the spatiotemporal evolution law of the floor cracks in the fully mechanized caving working face of the extra-thick coal seam.

[0013] S7: Analysis of the evolution of water channels in the mining floor. Based on numerical simulation results, the development patterns of water channels in the mining floor are analyzed under cyclic loading and unloading conditions such as advanced support pressure, floor pressure relief in goaf areas, and dynamic pressure from fully mechanized caving in extremely thick coal seams. This, combined with the distribution of underground aquifers, optimizes coal mining plans and achieves green mining.

[0014] According to the present invention, in step S1, the orientation of the main natural fractures within the coal seam and surrounding rock, as well as the fracture spacing along the borehole axis, is determined using a borehole peek method and the average length of the cores obtained from coring. Based on this measured data, Matlab is used to generate a set of spatially discrete points that satisfy the fracture spacing requirements in the three orthogonal directions of x, y, and z. The corresponding Voronoi polyhedron is then generated using Delaunay triangles. This polyhedron is then imported into the Part module to generate a natural fracture geometric model that approximates the measured spacing and orientation results.

[0015] In step S2, the following contents are included:

[0016] S2.1. Mesh the S1 geometric model and generate solid elements.

[0017] S2.2. Embed crack elements at the boundaries of solid elements.

[0018] In step S2.1, the S1 geometric model is meshed and the content of generating solid elements is:

[0019] In the Mesh module, use the Seed Part mode to seed the geometry created in S1. The spacing between the seeds should not exceed half the edge length of the Voronoi polyhedron. Use the Assign Mesh Controls button to control the shape of the elements to be hexahedrons, tetrahedrons, or primarily hexahedrons. Then, use the Mesh Part button to mesh the entire model. Furthermore, use the Assign Element Type button to change the meshed solid elements to the C3D8 solid element type.

[0020] In step S2.2, the contents of embedding crack elements at the boundaries of solid elements are:

[0021] (1) Update the node numbers of the C3D8 solid elements. After completing step S2.1, create a new job in the Job module and find the .inp file for the corresponding job in the corresponding working directory. Open the .inp file and find the maximum node and maximum element numbers in the numbers below the keywords *nset and *elset, and record them as nset respectively. m , else m ;

[0022] For C3D8 solid elements, any element (assuming the number is else j ) has 6 nodes; there are usually 2 to 6 solid elements around the unit. Assume that the nth i (n i ≠1, and n i ≤6) nodes are originally numbered b i , copy the node of the unit a times (2≤a≤6), a is the number of units with the same node. Then modify the node numbers of these copies and modify them in counterclockwise order to ((nset m ) 10 +b i ). Among them, (nset m ) 10 nset m The value obtained by rounding up to one digit. For example, nset m =126, then (nset m ) 10 After adding one digit, (nset m ) 10 = 200. After this node number modification method, the node number of the entity element can be updated, which is recorded as N i ; The entity element number remains unchanged and is still recorded as e j .

[0023] Add the keyword *Elset=solid after the keyword "*Part" and before "*Assembly", and then write the corresponding numbers in the modified element number sequence to obtain a solid element set.

[0024] (2) Create the crack element nodes and element numbers. In the .inp file, copy the array with different element numbers but 2 or 4 nodes with the same number. For this array, the element numbers of the two adjacent elements are recorded as e. j1 and e j2 , the node number shared by two adjacent units is recorded as n i1 , n i2 , n i3 , n i4The node number is the node number of the fracture element.

[0025] The number of the crack unit is based on the maximum number of solid units. m To be determined, and the initial value of the crack unit number is (elset m ) 10 , and then determine the maximum number of crack elements according to the number of "same node number" arrays. For example, if the maximum node number of the solid element is 115, the initial value of the crack element number is (elset m ) 10 =200; in the initial .inp file, if the number of the "same node number" array is 208, it means that there are 208 crack elements and the maximum crack element number is 200+208=408.

[0026] Save the original .inp file as .inp and add the keyword *Element type=coh3d8 to the new .inp file. Also, add the keyword *Elset=coh after the keyword "*Part" and before the keyword "*Assembly." Then, enter the corresponding numbers in the modified element and node numbering sequence. This will embed fracture elements at the boundaries of the solid elements and create a fracture element set.

[0027] In step S3, the following contents are included:

[0028] S3.1. Experimental determination of rock compression parameters;

[0029] S3.2. Experimentally determine the fracture mechanics parameters of Type I and Type II;

[0030] S3.3. Experimentally determine the rock shear friction parameters.

[0031] In step S3.1, the contents of experimentally determining the rock compression mechanical parameters are:

[0032] The laboratory temperature was controlled at (21±0.5)°C to minimize temperature disturbances on the test data acquisition equipment. A rock specimen with a diameter of 50 mm and a height of 100 mm was placed in a triaxial compression testing machine. Using displacement control, the loading rate was 0.02 mm / min, and three confining pressure conditions of 0 MPa, 5 MPa, and 10 MPa were considered. The loading process was divided into two stages. First, the confining pressure was applied to the predetermined confining pressure value through the hydraulic oil in the pressure chamber. Then, the deviatoric stress was generated by axial displacement control until the rock sample failed. The stress-strain curve of the rock during the compression failure process was recorded, and the elastic modulus and strength were calculated.

[0033] In step S3.2, the contents of experimentally determining the fracture mechanics parameters of type I and type II are:

[0034] Compact tensile tests were carried out on coal rock specimens to determine the type I fracture parameters. The minimum geometric dimension in the coal rock specimen is at least ten times the size of the crystal in the coal rock. It is recommended that the coal rock specimen used in the experiment be an oblate cylinder with a diameter of 100 mm and a height of 35 mm. A 0.5 mm thick diamond cutting blade was used to prefabricate a crack along the symmetry axis of the specimen. The crack was 15 mm long and the two symmetrical loading holes had a diameter of 13 mm. The distance from the center of the loading hole to the straight edge of the specimen was 13.0 mm, and the distance to the symmetry center of the specimen was 8 mm. The specimen was connected to the testing machine by a loading roller. During the experiment, the displacement loading rate was 0.02 mm / min, and the load-crack tip opening displacement curve was recorded during the deformation and destruction of the specimen to obtain the type I fracture strength, fracture energy, and fracture displacement w max .

[0035] A shear penetration test was carried out on the coal-rock specimen to determine the type II fracture parameters. The shear penetration specimen is a cylinder with a diameter of 50 mm and a height of 50 mm. Two concentric cracks of different depths were prefabricated on two parallel sections of the specimen along the axis of the cylinder, and a load was applied to the complete cylindrical area at the center. The prefabricated crack radius was 25 mm, the width was 1.5 mm, the upper prefabricated crack depth was 10 mm, and the lower prefabricated crack depth was 30 mm. When conducting the shear penetration test, a cylindrical loading steel head with a diameter of 25 mm was used to apply a load to the complete small cylindrical part in the center of the coal-rock specimen. The displacement loading method was adopted, and the loading speed was 0.02 mm / min. The load-displacement curve was recorded during the experiment. The type II fracture strength, fracture energy, and fracture displacement were calculated.

[0036] In step S3.3, the contents of the experimental determination of rock shear friction parameters are:

[0037] Prepare a cylindrical specimen with a diameter of 50mm and a height of 100mm. Then, use Brazilian splitting along the axial direction of the cylinder to prepare a specimen with a single rough crack. Place the specimen in a triaxial pressure chamber and increase the confining pressure to 0.2MP to ensure that all pressure and displacement sensors are normal. According to the measured ground stress data, the confining pressure and axial pressure are simultaneously increased to the test set value. Then, use the displacement loading method to perform shear displacement loading at a constant rate of 0.002mm / s, and set the shear displacement to 1 / 2 of the maximum loading displacement of the experimental equipment. Record the shear load-shear displacement curve during the experiment. The shear strength and shear modulus are calculated from this.

[0038] In step S4, the following contents are included:

[0039] S4.1: Set the set of cells to be deleted in the pre-processor;

[0040] S4.2: Prepare the VUSDFLD subroutine;

[0041] S4.3: Import the subroutine into the main program to simulate the working face advancement process;

[0042] Among them, in step S4.1, the content of the unit set to be deleted is set in the pre-processing program:

[0043] In the numerical model of step S2, the unit set to be deleted is set as follows:

[0044] (1) First, determine the range of coal seam units that need to be deleted in the numerical model;

[0045] (2) Copy the unit sets "*Elset=solid" and "*Elset=coh" created in step S2 respectively, and then select the coal seam units to be deleted by subtraction;

[0046] (3) According to the average size of the solid unit and the actual advancement speed of the working surface, the number of units to be deleted at each advancement distance is determined. Then, through the reduction method, the solid units and fracture units to be deleted at each advancement distance are further determined and their respective unit numbers are obtained.

[0047] In step S4.2, the contents of the VUSDFLD subroutine are:

[0048] First, in the compiler, define the parameters for importing the VUSDFLD subroutine. Use the "subroutinevusdfld" subroutine to import the parameters from the main program, including the element number, integration point coordinates, state variable parameters, stress, displacement, total time, analysis step time, time increment, direction, etc., and define each parameter as a scalar, vector, or tensor.

[0049] Second, define the set of cells that need to be deleted in each excavation step. After defining the total time totalTime, the time deleteTime required for each advance of a certain distance, the cell number nblock, the cell integration point coordinates coordMp, and the direction direct, fix the starting cell integration point coordinates coordMp(0) to be deleted. On this basis, when the ratio of deleteTime to totalTime is c (c is a constant), define all cells that meet this condition as "set 1"; when the ratio of deleteTime to totalTime is 2c, define all cells that meet this condition as "set 2"; and so on.

[0050] Third, set the initial state variable stateOld = 1 for all cells at the initial moment. Then, use the command stateNew = 0 to delete the cell sets "Set 1" and "Set 2" in sequence. This achieves the goal of cell deletion and work surface advancement.

[0051] In step S4.3, the content of importing the subroutine into the main program is:

[0052] In the Job module, create a new job task, select the subroutine file under the "General" button, and import the VUSDFLD.for subroutine in step S4.2 into the main program.

[0053] In step S6, the contents of the continuous-discrete coupling simulation program are:

[0054] The idea behind the continuous-discrete coupling simulation program is to calculate the nodal forces at all nodes of the numerical model using the finite element method under initial boundary conditions and ground stress conditions. After deleting the elements through the VUSDFLD.for subroutine in step S4.2, the finite element method is used to calculate the nodal forces of the solid elements. The stress is then transferred through the shared nodes of the solid elements and the fracture elements, and the displacement w of the fracture elements is calculated using the finite element method. The shear friction effect on the fracture surface is calculated using the discrete element method. The mechanical parameters of the solid elements (simulating complete coal blocks), fracture elements (simulating non-through fractures), and contact pairs (simulating through fractures) are determined according to step S3. Through this continuous-discrete coupling simulation program, the stress, strain, and displacement patterns of the floor rock mass during the advancement of the working face are solved.

[0055] The mathematical criterion for the generation of mining-induced floor cracks is as follows: Based on the above continuous-discrete coupling simulation program, the node forces and node displacements of the crack unit are calculated. The fracture energy of the crack unit at any moment is calculated by the following formula:

[0056]

[0057] Among them, G n , G s , G t are the fracture energies in the normal direction, the first tangential direction, and the direction perpendicular to the first tangential direction within the unit plane; σ n , σ s , σ t is the nodal force in three directions; S n , S s , S t are the node displacements in three directions.

[0058] When there is G n >G n实验 , G s >G s实验 , G t >G t实验 If any of the inequalities holds true, the crack unit is deleted, and the spatiotemporal evolution process of the initiation, expansion, and penetration of mining-induced floor cracks is numerically realized.

[0059] Compared with traditional numerical simulation methods based on the continuum medium assumption, such as finite element and Lagrangian methods, the continuous-discrete coupled simulation method of the floor failure depth of fully-mechanized caving mining working faces in extra-thick coal seams can monitor the development patterns of mining floor cracks in real time and reproduce the dynamic loads caused by the mining floor being subjected to advanced support pressure loading, unloading in the goaf, and overburden movement. This provides a simulation method suitable for the evolution of water-conducting cracks in the floor of fully-mechanized caving mining working faces in extra-thick coal seams. The technical innovations brought about by this method are:

[0060] 1. In fully mechanized top-coal caving working faces in extremely thick coal seams, the floor is subjected to a combination of static and dynamic loading and unloading conditions. Specifically, it is subjected to dynamic loads caused by the loading of advance support pressure, unloading in the goaf, and overburden movement. The continuous-discrete coupled simulation method of this invention can effectively simulate this special load, representing a significant improvement over traditional finite element, discrete element, and Lagrangian methods.

[0061] 2. It provides a powerful numerical simulation tool for studying the temporal and spatial evolution of floor damage and cracks caused by underground coal mining. It can numerically simulate the entire process of the floor rock transforming from a continuous medium to a discrete medium under mining stress conditions, enabling real-time monitoring of the development of floor cracks.

[0062] 3. It can be used to analyze the spatiotemporal evolution of water-conducting cracks in the mining floor under the influence of different mining parameters and different geological parameters, thereby optimizing mining parameters, realizing green mining of underground coal, and providing a basis for further floor grouting transformation. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 It is a schematic diagram of the continuous-discrete coupling numerical simulation process;

[0064] Figure 2 It is a set of numerical calculation models and corresponding fracture units;

[0065] Figure 3 It is a type I fracture specimen;

[0066] Figure 4 It is a type II fracture specimen;

[0067] Figure 5 is the distribution of cracks in the mining floor;

[0068] Figure 6 It is the surrounding rock displacement field under the mining conditions of fully mechanized caving mining working face in extra-thick coal seams;

[0069] Figure 7 It is the surrounding rock stress field under the mining conditions of fully mechanized caving mining working face in extra-thick coal seams.

[0070] The present invention will be further described in detail below with reference to the accompanying drawings and embodiments. DETAILED DESCRIPTION

[0071] The design idea of ​​the present invention is to establish a geometric model of random distribution of cracks in the coal seam floor by using the average length data of the cores of each rock layer obtained based on the borehole peek method and coring; divide the geometric model into solid units, and embed crack units at the boundary positions between natural cracks and solid units; obtain rock compression mechanics parameters, type I and type II fracture mechanics parameters, and rock shear friction parameters through experiments; compile a Vusdfld subroutine that reflects the advancement process of the coal mining working face, import it into the continuous-discrete coupling simulation main program, and numerically realize the entire process of shear sliding of the floor of the fully mechanized caving mining working face in ultra-thick coal seams from intact, fragmented, and structural surfaces.

[0072] This embodiment provides a finite discrete element simulation method for the mine pressure manifestation of a fully-mechanized caving working face in an extremely thick coal seam, which specifically includes the following steps:

[0073] S1: Determine the spacing of natural cracks in the basement and generate a geometric model including natural cracks;

[0074] In this embodiment, the borehole peek method and the average length of the core of each rock layer obtained by coring are used to determine the orientation of the main natural fractures in the coal seam and surrounding rock, and the fracture spacing along the axial direction of the borehole; based on this measured data, a set of spatial discrete points that meet the fracture spacing in the three orthogonal directions of x, y, and z are generated through Matlab; then the Delaunay triangle is used to generate the corresponding Voronoi polyhedron; the polyhedron is imported into the Part module to generate a natural fracture geometric model that is close to the measured spacing and orientation results.

[0075] By drilling and coring, the spacing of the main structural planes in each type of rock layer is determined. Figure 2 As shown. Through the Voronoi polyhedron generation program, a three-dimensional geometric model with length x = 200m, height y = 174m, and thickness z = 1m is established, as shown Figure 2 shown.

[0076] S2: Divide the geometric model into units and embed fracture units at the locations of natural fractures;

[0077] Specifically include the following:

[0078] S2.1: Mesh the S1 geometric model and generate solid elements. The specific contents are:

[0079] In the Mesh module, use the Seed Part mode to seed the geometry created in S1. The spacing between the seeds should not exceed half the edge length of the Voronoi polyhedron. Use the Assign Mesh Controls button to control the shape of the elements to hexahedron, tetrahedron, or primarily hexahedron. Then, use the Mesh Part function to mesh the entire model. Furthermore, use the Assign Element Type button to change the solid element type to C3D8.

[0080] S2.2: Embed crack elements at the boundaries of solid elements. Specifically:

[0081] (1) Update the node number of C3D8 solid element

[0082] After completing step S2.1, create a new job in the Job module and find the .inp file for the corresponding job in the corresponding working directory. Open the .inp file and find the maximum node and maximum element numbers in the numbers below the keywords *nset and *elset, and record them as nset respectively. m , else m ;

[0083] For C3D8 solid elements, assume the number is elset j Each unit has 6 nodes; there are usually 2 to 6 solid units around the unit; assuming that the nth unit i The original number of the node is b i , n i ≠1, and n i ≤6, copy the node of the unit a times, a is the number of units with shared nodes, and 2≤a≤6; then modify the numbers of the copied nodes and modify them in counterclockwise order to ((nset m ) 10 +b i ); where (nset m ) 10 nset m The value obtained by rounding up to one digit; after this modification of the node number, the node number of the entity element can be updated, which is recorded as N i ; The entity element number remains unchanged and is still recorded as e j ;

[0084] Add the keyword *Elset=solid after the keyword "*Part" and before "*Assembly", and then write the corresponding numbers in the modified unit number sequence to obtain a solid unit set.

[0085] (2) Create crack element nodes and element numbers

[0086] In the .inp file, copy an array with different element numbers but 2 or 4 identical node numbers; for such an array, the element numbers of the two adjacent elements are recorded as e j1 and e j2 , the node number shared by two adjacent units is recorded as n i1 , n i2 , n i3 , n i4 , the node number is the node number of the crack element;

[0087] The number of the crack unit is based on the maximum number of solid units. m To be determined, and the initial value of the crack unit number is (elset m ) 10 , and then determine the maximum number of crack elements based on the number of "same node number" arrays; in the initial .inp file, if the number of "same node number" arrays is 208, it means there are 208 crack elements and the maximum number of crack elements is 200 + 208 = 408;

[0088] Save the initial .inp file and add the keyword *Element type=coh3d8 in the new .inp file. At the same time, add the keyword *Elset=coh after the keyword "*Part" and before "*Assembly", and then write the corresponding numbers in the modified unit and node numbering order; thereby embedding the fracture unit at the boundary of the solid unit and obtaining a fracture unit set.

[0089] In the Mesh module, use the Seed Part mode to seed the geometry. Set the seed spacing to 3 meters. Use the Assign Mesh Controls command to set the element shape to a hexahedron, and use the Mesh Part command to mesh the entire model. Then, use the Assign Element Type command to change the solid element type to a C3D8 solid element. Group all solid elements into a unit collection named "solid."

[0090] Crack elements are embedded at the boundaries of solid elements. The element type is coh3d8, and all crack elements are combined into an element set named "coh" ( Figure 2 (indicated by the broken line).

[0091] S3: Experimentally determine rock compression parameters, type I and type II fracture parameters, and rock shear friction parameters;

[0092] Specifically include the following:

[0093] S3.1: Experimental determination of rock compression parameters

[0094] The test chamber temperature was controlled at 21°C ± 0.5°C to minimize temperature disturbances on the test data acquisition equipment. A rock specimen with a diameter of 50 mm and a height of 100 mm was placed in a triaxial compression testing machine. Displacement control was used with a loading rate of 0.02 mm / min. Three confining pressures, 0 MPa, 5 MPa, and 10 MPa, were considered. The loading process was divided into two stages. First, confining pressure was applied to a predetermined value using hydraulic oil in the pressure chamber. Then, deviatoric stress was generated through axial displacement control until the rock sample failed. The stress-strain curve of the rock during the compression failure process was recorded, and the elastic modulus and strength were calculated.

[0095] S3.2: Experimental determination of type I and type II fracture mechanics parameters

[0096] Compact tensile tests were carried out on coal and rock specimens to determine the mode I fracture parameters:

[0097] The minimum geometric dimension of the coal rock specimen is at least ten times the size of the crystals within the coal rock. The coal rock specimen used in the experiment is an oblate cylinder with a diameter of 100 mm and a height of 35 mm. A 15 mm long crack is prefabricated along the symmetry axis of the specimen using a 0.5 mm thick diamond cutting blade. The diameter of the two symmetrical loading holes is 13 mm. The distance from the center of the loading hole to the straight edge of the specimen is 13.0 mm, and the distance to the symmetry center of the specimen is 8 mm. The specimen is connected to the testing machine via a loading roller.

[0098] During the experiment, the displacement loading rate was 0.02 mm / min, and the load-crack tip opening displacement curve was recorded during the deformation and destruction of the specimen to obtain the mode I fracture strength, fracture energy, and fracture displacement w max ;

[0099] Shear penetration tests were carried out on coal-rock specimens to determine the type II fracture parameters:

[0100] The shear penetration specimen is a cylinder with a diameter of 50 mm and a height of 50 mm. Two concentric cracks of different depths are prefabricated on two parallel sections of the specimen along the axis of the cylinder. Load is applied to the complete cylindrical area at the center. The prefabricated crack radius is 25 mm, the width is 1.5 mm, the upper prefabricated crack depth is 10 mm, and the lower prefabricated crack depth is 30 mm. When conducting the through-shear test, a cylindrical loading steel head with a diameter of 25 mm is used to apply load to the complete small cylindrical part in the center of the coal and rock specimen. Displacement loading is adopted, and the loading speed is 0.02 mm / min.

[0101] During the experiment, the load-displacement curve was recorded, from which the mode II fracture strength, fracture energy, and fracture displacement were calculated;

[0102] S3.3: Experimental determination of rock shear friction parameters

[0103] Prepare a cylindrical specimen with a diameter of 50mm and a height of 100mm. Then, use Brazilian splitting along the axial direction of the cylinder to prepare a specimen with a single rough crack. Place the specimen in a triaxial pressure chamber and increase the confining pressure to 0.2MP. Ensure that all pressure and displacement sensors are functioning properly. Based on the measured in-situ stress data, increase the confining pressure and axial pressure to the test set values ​​simultaneously. Then, use displacement loading, apply shear displacement loading at a constant rate of 0.002mm / s, and set the shear displacement to 1 / 2 of the maximum loading displacement of the experimental equipment.

[0104] During the experiment, the shear load-shear displacement curve was recorded, from which the shear strength and shear modulus were calculated.

[0105] In this embodiment, cores were taken from typical rock formations on the bottom plate and prepared into standard cylindrical specimens (50 mm in diameter and 100 mm in height), cylindrical specimens (50 mm in diameter and 200 mm in height with prefabricated cracks), and oblate cylindrical specimens (50 mm in diameter and 50 mm in height with prefabricated annular cracks). Figure 3 and Figure 4 As shown. Compression tests, mode I fracture, mode II fracture, and shear friction tests were performed to obtain load-displacement curves, and parameters such as elastic modulus, peak strength, fracture energy, and fracture displacement were calculated. The material parameters are shown in Tables 1 to 3 below.

[0106] Table 1: Mechanical parameters of solid elements

[0107] Elastic modulus / GPa Poisson's ratio Compression strength / MPa siltstone 10.7 0.21 32 medium sandstone 20.6 0.22 47 coal 1.6 0.35 6 mudstone 8.0 0.24 19

[0108] Table 2: Mechanical parameters of fracture elements

[0109]

[0110] Table 3: Contact pair mechanical parameters

[0111] Shear modulus / GPa Shear peak strength / MPa Shear residual strength / MPa 7.3 12.1 6.7

[0112] S4: Compile the VUSDFLD subroutine to automatically delete coal aggregate units and simulate the working face advancement process;

[0113] Specifically include the following:

[0114] S4.1: Set the set of cells to be deleted in the pre-processor

[0115] In the numerical model of step S2, the unit set to be deleted is set as follows:

[0116] (1) First, determine the range of coal seam units that need to be deleted in the numerical model;

[0117] (2) Copy the unit sets "*Elset=solid" and "*Elset=coh" created in step S2 respectively, and then select the coal seam units to be deleted by subtraction;

[0118] (3) According to the average size of the solid unit and the actual advancement speed of the working surface, the number of units that need to be deleted at each advancement distance is determined. Then, through the reduction method, the solid units and fracture units that need to be deleted at each advancement distance are further determined, and their respective unit numbers are obtained;

[0119] S4.2: Compile VUSDFLD subroutine

[0120] First, in the compiler, define the parameters for importing the VUSDFLD subroutine. Use the "subroutinevusdfld" subroutine to import the parameters in the main program, including element number, integration point coordinates, state variable parameters, stress, displacement, total time, analysis step time, time increment, direction, etc., and define each parameter as a scalar, vector, and tensor.

[0121] Second, define the set of cells that need to be deleted in each excavation step; based on the definition of total time totalTime, the time deleteTime required for each advance of a certain distance, the cell number nblock, the cell integration point coordinates coordMp, and the direction direct, fix the starting cell integration point coordinates coordMp(0) to be deleted;

[0122] On this basis, when the ratio of deleteTime to totalTime is c, c is a constant, and all cells that meet this condition are defined as "set 1";

[0123] When the ratio of deleteTime to totalTime is 2c, define all cells that meet this condition as "set 2"; and so on;

[0124] Third, set the initial state variable stateOld=1 of all units at the initial moment, and then delete the unit sets such as "Set 1" and "Set 2" in sequence through the stateNew=0 command; thereby achieving the goal of unit deletion and working surface advancement.

[0125] S4.3: Import the subroutine into the main program to simulate the working face advancement process

[0126] In the Job module, create a new job task, select the subroutine file under the "General" button, and import the VUSDFLD.for subroutine in step S4.2 into the main program.

[0127] In this embodiment, in the numerical model, the number of units that need to be deleted and their corresponding boundary positions are determined for each working face advancement distance based on the average size of the solid units and the actual advancement speed of the working face. The unit sets "*Elset=solid" and "*Elset=coh" set in step S2 are copied respectively. Then, the coal seam units that need to be deleted are framed by subtraction, and the solid units that need to be deleted are named solid-KW01 and solid-KW02 respectively; the fracture units that need to be deleted are named coh-KW01 and coh-KW02 respectively.

[0128] In the Materials module of the preprocessing process, define the material state variable parameter Depvar = 1. In the VUSDFLD subroutine, define total time totalTime = 10 and the time required for each advance of a certain distance deleteTime = 5. The first time, the solid elements to be deleted are numbered nblock = 287 to 309, and the fracture elements to be deleted are numbered nblock = 10103 to 10135. The second time, the solid elements to be deleted are numbered nblock = 462 to 490, and the fracture elements to be deleted are numbered nblock = 10294 to 10326. When the ratio of deleteTime to totalTime is 0.5, stateNew("1", 1) = 0, deleting the first "unit set". When the ratio of deleteTime to totalTime is 1, stateNew("2", 1) = 0, deleting the second "unit set". This achieves the goal of advancing the coal mining face.

[0129] In the Job module, create a new job task Job15, select the subroutine file under the "General" button, and import the VUSDFLD.for subroutine into the main program.

[0130] S5: Apply numerical model boundary conditions, input material mechanical parameters, set incremental step size, and set output parameters, including stress, strain, displacement, unit node force, and state variables of crack units;

[0131] In this embodiment, the material parameters (Tables 1 to 3) are input in the Materials module. In the Assembly module, the normal displacements of the five boundaries of the numerical calculation model, except for the upper surface, are constrained. The upper surface is a free boundary. Then, according to the advancement speed of the coal mining face and the model size, one analysis step is set in the Step module, geometric nonlinearity is turned on, the time length is 10, and the time increment is set to 0.001; the output variables are set to stress, strain, friction displacement, friction stress, state variable STATUS, etc. The initial stress of the model is determined by self-balancing.

[0132] S6: The model was solved by continuous-discrete coupling simulation program to obtain the spatiotemporal evolution of floor cracks in fully mechanized caving working faces in extra-thick coal seams;

[0133] The contents of the continuous-discrete coupled simulation program are:

[0134] Under the initial boundary conditions and ground stress conditions, the nodal forces of all nodes of the numerical model are calculated using the finite element method. After deleting the unit through the VUSDFLD.for subroutine in step S4.2, the nodal forces of the solid unit are calculated using the finite element method. Then, the stress is transferred through the shared nodes of the solid unit and the fracture unit, and the displacement w of the fracture unit is calculated using the finite element method.

[0135] The shear friction effect of the fracture surface is calculated using the discrete element method. The mechanical parameters of the solid unit, i.e., simulating the complete coal rock block, the fracture unit, i.e., simulating the non-through fracture, and the contact pair, i.e., simulating the through fracture, are determined according to step S3. The stress, strain, displacement, and other laws of the floor rock mass during the advancement of the working face are solved through a continuous-discrete coupling simulation program.

[0136] The mathematical criterion for the generation of mining-induced floor cracks is as follows: Based on the above continuous-discrete coupling simulation program, the node forces and node displacements of the crack unit are calculated. The fracture energy of the crack unit at any moment is calculated by the following formula:

[0137]

[0138] Among them, G n , G s , G t are the fracture energies in the normal direction, the first tangential direction, and the direction perpendicular to the first tangential direction within the unit plane; σ n , σ s , σ t is the nodal force in three directions; S n , S s , S t is the node displacement in three directions;

[0139] When there is G n >Gn实验 , G s >G s实验 , G t >G t实验 If any of the inequalities holds true, the crack unit is deleted, and the spatiotemporal evolution process of the initiation, expansion, and penetration of mining-induced floor cracks is numerically realized.

[0140] In this embodiment, under the Job module, for the newly created job task Job15, the number of parallel computing cores is set to 48 in the Edit Job tab, and the calculation is submitted.

[0141] S7: Analysis of the evolution of water channels in mining-induced floor

[0142] Based on the results of numerical simulation, the development law of the floor water guide channel under cyclic loading and unloading conditions such as advanced support pressure, goaf floor unloading, and dynamic pressure of ultra-thick coal seam fully mechanized caving is analyzed. Then, combined with the distribution of underground aquifers, the coal mining plan is optimized to achieve green coal mining.

[0143] In this example, numerical calculations were performed to determine the spatiotemporal evolution of floor cracks in fully mechanized top-coal caving working faces in extra-thick coal seams. Using a continuous-discrete coupling program, the nodal forces and displacements of fracture units in each rock layer of the floor under the influence of mining were numerically calculated, and the fracture energies of each unit were integrated to determine the appropriate fracture energy. Using experimental results for Type I and Type II fracture energies of siltstone, medium sandstone, coal, and mudstone as mathematical criteria, the numerically calculated fracture energies were compared with the experimental results to determine whether to delete the fracture units. This numerically captures the spatiotemporal evolution of the initiation, expansion, and penetration of floor cracks under mining.

[0144] The analysis results show that:

[0145] The distribution and evolution of mining floor cracks, such as Figure 5 As shown by the black line. From the calculation results, it can be seen that the mining of the extra-thick coal seam in the model will cause its floor damage. The floor within the goaf and close to the working face will rupture as a whole according to the plastic slip line theory; the floor slightly away from the working face within the goaf is in a pressure-relief state; and the extra-thick coal seam roof will collapse ( Figure 6 and Figure 7 Under this cyclic loading and unloading condition (advance support pressure loading-goaf unloading-dynamic load caused by overburden movement), the damage and rupture of the floor will be further aggravated.

[0146] Outside the goaf and within 15 meters of the coal face, the coal seam experienced tensile failure, while the floor exhibited a combined tension / shear failure mode. Water-conducting fissures in the floor advanced 25 meters ahead of the coal wall at the working face, resulting in floor failure depths of up to 30 meters. Compared to fully mechanized caving mining in medium-thick seams, fully mechanized caving mining in extra-thick seams exhibited greater floor failure depths.

[0147] Compared with existing floor failure numerical methods, the finite discrete element simulation method for rock pressure manifestation in fully-mechanized caving working faces in extra-thick coal seams developed in this embodiment not only demonstrates the process by which floor rock transforms from a continuous medium to a discrete medium under mining-induced stress conditions, but also deeply reflects the static-dynamic loading and unloading conditions under which the floor is subjected to dynamic loads caused by advance support pressure loading, unloading in the goaf, and overburden movement, as well as the resulting floor damage and fracture. This provides a powerful numerical simulation tool for the evolution of water-conducting fractures in the floor of fully-mechanized caving working faces in extra-thick coal seams.

[0148] It should be noted that the above embodiments are preferred examples of the present invention and the present invention is not limited to the above embodiments. Any addition, improvement, modification, or equivalent changes made by any person skilled in the art without departing from the scope of the technical solution of the present invention shall fall within the scope of protection defined by the claims of the present invention.

Claims

1. A continuous-discrete coupling simulation method for floor failure depth of fully mechanized caving mining working face in extra-thick coal seams, characterized by: The following steps are involved: S1: Determine the spacing of natural cracks in the basement and generate a geometric model including natural cracks; S2: Divide the geometric model into units and embed fracture units at the locations of natural fractures; S3: Experimentally determine rock compression parameters, type I and type II fracture parameters, and rock shear friction parameters; S4: Compile the VUSDFLD subroutine to automatically delete coal aggregate units and simulate the working face advancement process; S5. Apply numerical model boundary conditions, input material mechanical parameters, set the incremental step size, and set output parameters, including stress, strain, displacement, element node force, and state variables of crack elements; S6: Solve the model through the continuous-discrete coupling simulation program to obtain the spatiotemporal evolution law of the floor cracks in the fully mechanized caving working face of the extra-thick coal seam. In step S6, the content of the continuous-discrete coupling simulation program is: Under the initial boundary conditions and ground stress conditions, the nodal forces of all nodes of the numerical model are calculated by the finite element method. After deleting the unit through the VUSDFLD.for subroutine, the nodal forces of the solid unit are calculated by the finite element method. Then, the stress is transferred through the shared nodes of the solid unit and the fracture unit, and the displacement w of the fracture unit is calculated by the finite element method. The shear friction effect of the fracture surface is calculated using the discrete element method. The mechanical parameters of the solid unit, i.e., simulating the complete coal rock block, the fracture unit, i.e., simulating the non-through fracture, and the contact pair, i.e., simulating the through fracture, are determined according to step S3. The stress, strain, and displacement laws of the floor rock mass during the advancement of the working face are solved through a continuous-discrete coupling simulation program. The mathematical criterion for determining the generation of mining-induced floor cracks is as follows: the nodal forces and nodal displacements of the crack units are calculated according to the above-mentioned continuous-discrete coupling simulation program; the fracture energy of the crack unit at any moment is calculated using the following formula: Among them, G n , G s , G t are the fracture energies in the normal direction, the first tangential direction, and the direction perpendicular to the first tangential direction within the unit plane; σ n , σ s , σ t is the nodal force in three directions; S n , S s , S t is the node displacement in three directions; When there is G n >G n实验 , G s > G s实验 , G t > G t实验 If any of the inequalities holds true, the crack unit is deleted, and the spatiotemporal evolution process of the initiation, expansion, and penetration of mining-induced floor cracks is numerically realized. S7: Analysis of the evolution of water channels in mining-induced floor Based on the results of numerical simulation, the development law of the floor water guide channel under cyclic loading and unloading conditions such as advanced support pressure, goaf floor unloading, and dynamic pressure of ultra-thick coal seam fully mechanized caving is analyzed. Then, combined with the distribution of underground aquifers, the coal mining plan is optimized to achieve green coal mining.

2. The method according to claim 1, wherein In step S1, the borehole peek method and the average length of the core of each rock layer obtained by drilling coring are used to determine the orientation of the main natural cracks in the coal seam and surrounding rock, and the crack spacing along the axial direction of the borehole; based on this measured data, a set of equations that meet the requirements are generated by Matlab. x , y , z The spatial discrete points that meet the crack spacing in three orthogonal directions are then generated using Delaunay triangles to generate the corresponding Voronoi polyhedrons. The polyhedrons are imported into the Part module to generate a natural crack geometry model that is close to the measured spacing and orientation results.

3. The method according to claim 1, wherein In step S2, the following contents are included: S2.1: Mesh the S1 geometric model and generate solid elements. The specific contents are: In the Mesh module, use the Seed Part mode to seed the geometry created in S1. The spacing between the seeds should not exceed 1 / 2 of the edge length of the Voronoi polyhedron. Use Assign Mesh Controls to control the shape of the elements to hexahedron, tetrahedron, or mainly hexahedron. Then, use Mesh Part to mesh the entire model. Furthermore, under the Assign Element Type button, change the meshed solid elements to C3D8 solid element type. S2.2: Embed crack elements at the boundaries of solid elements. Specifically: (1) Update the node number of C3D8 solid element After completing step S2.1, create a new job in the Job module and find the .inp file for the corresponding job in the corresponding working directory; Open the .inp file and find the maximum node and maximum element numbers in the numbers below the keywords *nset and *elset, and record them as nset m , else m ; For C3D8 solid element, assuming the number is else j Each unit has 6 nodes; there are usually 2 to 6 solid units around the unit; n i The original number of nodes is b i , n i ≠1, and n i ≤6, copy the node of the unit a Second-rate, a is the number of elements with common nodes, and 2≤ a ≤6; then modify the node numbers of these copied nodes and modify them in counterclockwise order to (( nset m ) 10 + b i ); in, ( nset m ) 10 for nset m The value obtained by rounding up one digit; after this modification of the node number, the node number of the entity element can be updated, which is recorded as N i ; The entity element number remains unchanged and is still recorded as e j ; Add the keyword *Elset=solid after the keyword "*Part" and before "*Assembly", and then write the corresponding numbers in the modified element number sequence to obtain the solid element set; (2) Create fracture element nodes and element numbers In the .inp file, copy an array with different element numbers but 2 or 4 identical node numbers; for such an array, record the element numbers of the two adjacent elements as e j1 and e j2 , the node numbers shared by two adjacent units are recorded as n i1 , n i2 , n i3 , n i4 , the node number is the node number of the crack element; The crack unit number is based on the maximum number of solid units else m The initial value of the crack unit number is determined as ( else m ) 10 , and then determine the maximum number of crack elements based on the number of "same node number" arrays; in the initial .inp file, if the number of "same node number" arrays is 208, it means there are 208 crack elements, and the maximum number of crack elements is 200 + 208 = 408; Save the original .inp file and add the keyword *Element type= coh3d8 to the new .inp file. Also, add the keyword *Elset=coh after the keyword "*Part" and before "*Assembly." Then, write the corresponding numbers in the modified element and node numbering sequence. This will embed fracture elements at the boundaries of the solid elements and create a fracture element set.

4. The method according to claim 1, wherein In step S3, the following contents are included: S3.1: Experimental determination of rock compression parameters The test chamber temperature was controlled at 21°C ± 0.5°C to minimize temperature disturbances on the test data acquisition equipment. A rock specimen with a diameter of 50 mm and a height of 100 mm was placed in a triaxial compression testing machine. Displacement control was used with a loading rate of 0.02 mm / min. Three confining pressures (0 MPa, 5 MPa, and 10 MPa) were considered, and the loading process was divided into two stages. First, confining pressure was applied to a predetermined value using hydraulic oil in the pressure chamber. Then, deviatoric stress was generated through axial displacement control until the rock specimen failed. The rock stress-strain curve during the compression failure process was recorded, and the elastic modulus and strength were calculated. S3.2: Experimental determination of type I and type II fracture mechanics parameters Compact tensile tests were carried out on coal and rock specimens to determine the mode I fracture parameters: The minimum geometric dimension of the coal rock specimen is at least ten times the size of the crystals within the coal rock. The coal rock specimen used in the experiment is an oblate cylinder with a diameter of 100 mm and a height of 35 mm. A 15 mm long crack is prefabricated along the symmetry axis of the specimen using a 0.5 mm thick diamond cutting blade. Two symmetrical loading holes have a diameter of 13 mm. The distance from the center of the loading hole to the straight edge of the specimen is 13.0 mm, and the distance to the center of symmetry is 8 mm. The specimen is connected to the testing machine via a loading roller. During the experiment, the displacement loading rate was 0.02 mm / min, and the load-crack tip opening displacement curve was recorded during the deformation and destruction of the specimen to obtain the mode I fracture strength, fracture energy, and fracture displacement w max ; Shear penetration tests were carried out on coal-rock specimens to determine the type II fracture parameters: The shear penetration specimen is a cylinder with a diameter of 50 mm and a height of 50 mm. Two concentric cracks of different depths are prefabricated on two parallel sections of the specimen along the axis of the cylinder. Load is applied to the complete cylindrical area at the center. The prefabricated crack radius is 25 mm, the width is 1.5 mm, the upper prefabricated crack depth is 10 mm, and the lower prefabricated crack depth is 30 mm. During the shear penetration test, a cylindrical loading steel head with a diameter of 25 mm is used to apply load to the complete small cylindrical part in the center of the coal and rock specimen. Displacement loading is adopted at a loading speed of 0.02 mm / min. During the experiment, the load-displacement curve was recorded, from which the mode II fracture strength, fracture energy, and fracture displacement were calculated; S3.3: Experimental determination of rock shear friction parameters Prepare a cylindrical specimen with a diameter of 50 mm and a height of 100 mm. Then, use Brazilian splitting along the axial direction of the cylinder to prepare a specimen with a single rough crack. Place the specimen in a triaxial pressure chamber and increase the confining pressure to 0.2 MPa to ensure that all pressure and displacement sensors are functioning properly. Based on the measured ground stress data, the confining pressure and axial pressure were simultaneously increased to the test set values. Then, displacement loading was adopted, and shear displacement loading was performed at a constant rate of 0.002 mm / s. The shear displacement was set to 1 / 2 of the maximum loading displacement of the experimental equipment. During the experiment, the shear load-shear displacement curve was recorded; the shear strength and shear modulus were calculated from it.

5. The method according to claim 1, wherein In step S4, the following contents are included: S4.1: Set the set of cells to be deleted in the pre-processor In the numerical model of step S2, the unit set to be deleted is set as follows: (1) First, determine the range of coal seam units that need to be deleted in the numerical model; (2) Copy the unit sets "*Elset=solid" and "*Elset=coh" created in step S2 respectively, and then select the coal seam units to be deleted by subtracting them; (3) According to the average size of the solid unit and the actual advancement speed of the working face, the number of units that need to be deleted at each advancement distance is determined; then, through the reduction method, the solid units and fracture units that need to be deleted at each advancement distance are further determined, and their respective unit numbers are obtained; S4.2: Compile VUSDFLD subroutine First, in the compiler, define the parameters for importing the VUSDFLD subroutine. Use the "subroutine vusdfld" subroutine to import the parameters in the main program, including the element number, integration point coordinates, state variable parameters, stress, displacement, total time, analysis step time, time increment, and direction. Define each parameter as a scalar, vector, or tensor. Second, define the set of cells that need to be deleted in each excavation step; based on the definition of total time totalTime, the time deleteTime required for each advance of a certain distance, the cell number nblock, the cell integration point coordinates coordMp, and the direction direct, fix the starting cell integration point coordinates coordMp(0) to be deleted; On this basis, when the ratio of deleteTime to totalTime is c, c is a constant, and all cells that meet this condition are defined as "set 1"; When the ratio of deleteTime to totalTime is 2c, define all cells that meet this condition as "set 2"; and so on; Third, set the initial state variable stateOld=1 for all cells at the initial moment, and then delete the cell sets "set 1" and "set 2" in sequence through the stateNew=0 command; thereby achieving the goal of cell deletion and working face advancement; S4.3: Import the subroutine into the main program to simulate the working face advancement process In the Job module, create a new job task, select the subroutine file under the "General" button, and import the VUSDFLD.for subroutine in step S4.2 into the main program.

Images