A finite-discrete element simulation method for strata behavior of fully-mechanized caving face in super-thick coal seam

Abstract

The present application relates to a thick coal seam fully mechanized working face strata behavior finite discrete element simulation method, according to the research area to establish numerical calculation model;Determine the entity unit size of each rock layer, mesh division is carried out to numerical model;Zero thickness fracture element is generated at the boundary of entity unit;Set analysis step, and set displacement boundary set of node force release under each analysis step;Apply the initial condition of numerical model whole, boundary condition, input material compression, fracture mechanics parameter, set output variable;Through finite discrete element simulation program, numerical solution node force under different analysis step, reveal the advance support pressure, support working resistance and other strata behavior law in the process of thick coal seam fully mechanized working face mining. Effective simulation of thick coal seam fully mechanized working face advance support pressure, support working resistance, and overburden deformation / rupture / movement and the quantitative relationship among the three, for making roof (coal) grouting reinforcement, safe and efficient mining thick coal seam provides the basis.

Description

Technical Field

[0001] This invention belongs to the field of roof disaster mechanism and prevention technology, specifically involving a finite discrete element simulation method for mine pressure manifestation in fully mechanized longwall mining faces of extra-thick coal seams. Background Technology

[0002] In my country's underground coal mining, over 85% of operations utilize longwall mining and caving methods for roof management. Longwall mining necessitates effective roof control. After decades of development, a relatively complete theoretical system has been established regarding the movement of mining-induced overburden strata and the manifestation of mine pressure. Employing reasonable numerical methods to predict mine pressure patterns is crucial for optimizing parameters in fully mechanized longwall mining faces of extra-thick coal seams and formulating roof disaster prevention and control measures.

[0003] Domestic and international scholars mainly employ numerical methods such as the finite element method (FEM), Lagrange multiplier method (LMM), and discrete element method (DEM) to conduct numerical simulation studies of mine pressure manifestation in longwall mining faces. Based on the assumption of a continuous medium, FEM and LMM methods are established, combined with theories of elasticity, plasticity, and damage mechanics to simulate the evolution of stress, displacement, plastic zone, and damage zone in the roof during mining. Based on the assumption of a discrete block, the DEM method is established, often simplifying the discrete block into a rigid body or elastic body, and using Newton's second law to simulate the roof collapse process during mining. However, the roof strata are neither ideal continuous media nor ideal discrete blocks, but rather a quasi-continuous medium. Therefore, a reasonable and effective simulation of the entire process of continuous deformation, fracturing, and movement of the overlying strata during mining is fundamental to accurately simulating mine pressure manifestation in fully mechanized longwall mining faces of extra-thick coal seams. Summary of the Invention

[0004] The purpose of this invention is to provide a finite discrete element simulation method for the mining pressure manifestation in fully mechanized longwall mining faces of extra-thick coal seams.

[0005] To achieve the above objectives, the present invention adopts the following technical solution:

[0006] A finite discrete element method for simulating mine pressure manifestation in fully mechanized longwall mining faces of extra-thick coal seams includes the following steps:

[0007] S1: Establish a numerical calculation model based on the size of the study area;

[0008] S2: Determine the element size and mesh the numerical model;

[0009] S3: Generate zero-thickness crack elements at the boundary of solid elements;

[0010] S4: Set up the analysis step, and set the set of displacement boundaries for nodal force release in each analysis step;

[0011] S5: Apply boundary conditions to the entire numerical model, input material compression and fracture mechanics parameters, and set output variables;

[0012] S6: Numerical calculations are performed using a finite discrete element simulation program, and the nodal forces at the leading working face and the support positions behind the working face are output under different analysis steps.

[0013] S7: Based on nodal force data, reveal the manifestation patterns of mine pressures such as advance support pressure and support working resistance during the mining of extra-thick coal seams in fully mechanized longwall faces. Based on this, and combined with the engineering geological conditions of the working face, formulate roof disaster prevention and control measures to achieve safe and efficient coal mining.

[0014] According to the present invention, in step S1, geometric data such as the dimensions of the working face (length × width = a × b), mining depth (h), and mining thickness (c) are determined. To avoid boundary effects, the model dimensions in the horizontal direction should be appropriately increased, with the increased dimensions being d and e respectively. Thus, a numerical model of length × width × height = (a + d) × (a + e) ​​× c is established.

[0015] Specifically, step S2 includes the following:

[0016] S2.1: Determine the solid unit dimensions of each rock stratum;

[0017] S2.2: Mesh the numerical model.

[0018] In step S2.1, determining the dimensions of the solid units for each rock stratum includes:

[0019] Core samples were obtained from each rock stratum using borehole coring. The strata were grouped according to stratigraphy and lithology, and the average length of each core segment was calculated. Simultaneously, the joint and fracture spacing in the roadway and coal face was observed underground. Based on this, the dimensions of each rock stratum solid element in the numerical model were determined as f1, f2, f3, ...

[0020] In step S2.2, the mesh generation for the numerical model includes:

[0021] In the Mesh module, the Seed Part mode is used to seed the numerical model established in S1. Seeding is based on size, using the number f with the highest frequency in the size set {f1, f2, f3, ...} of each rock layer entity unit in step S2.1. iBased on this, the seed spacing for the entire numerical model is set. Then, using the Seed Edges mode, the seed spacing for each rock layer is set to f1, f2, f3, ... The shape of the control unit is selected via Assign Mesh Controls, prioritizing hexahedrons, followed by hexahedrons as the primary shape, and finally tetrahedrons. The entire model is then meshed using Mesh Part. Under the Assign Element Type button, the meshed solid elements are changed to C3D8 solid element type.

[0022] Furthermore, in step S3, the process of generating zero-thickness crack elements at the boundary of the solid element is as follows:

[0023] (1) Locate the maximum node number and cell number in step S2. Set the working directory, create a new job under the Job module, and then locate the .inp script file of the newly created job in that working directory. In the .inp file, find the keyword node set *nset and cell set *elset. Copy the numbers below the keywords into Excel, and use the MAX function to find the maximum node number E of the entity cell. max and the largest unit number N max .

[0024] (2) Re-encode the node numbers of the solid elements. For any hexahedral solid element or tetrahedral solid element (e element number...) j ), with 8 nodes and 4 nodes respectively. Unit e j There are a total of g elements in the surrounding area, generally g = 4 (tetrahedral elements) or g = 6 (hexahedral elements). If element e j The initial number of the i-th node is n. i (i.e., 2≤i≤6 for hexahedral elements and 2≤i≤4 for tetrahedral elements), and copy the node g times. Then modify the node numbers of these copied nodes, changing the node numbers of all elements except the first one to (N) sequentially. max +n i After the above steps, the node number of the entity element can be updated, denoted as n. I The entity unit number remains unchanged and is still denoted as e. j .

[0025] In the script file, add the keyword "*Elset=solid" after "*Part" and before "*Assembly", then sequentially write the node numbers according to the modified cell numbers to obtain the updated entity cell node numbers. Save it as a new file named new.inp, containing the above modifications.

[0026] (3) Generate crack element numbers and node numbers. In the initial script file, copy arrays with different element numbers but three (tetrahedral elements) or four (hexahedral elements) identical node numbers, including both node numbers and element numbers. For two adjacent elements e with shared nodes... j1 and e j2 The shared node is denoted as n. i1 n i2 n i3 n i4 Therefore, the node number n after the root step S3 update is... I1 n I2 n I3 n I4 This refers to the node number of the crack element.

[0027] The element numbering of the crack element is determined as follows: the maximum node number of the solid element is E. max Therefore, the smallest element of the crack element is numbered E. LXmin The value is an integer that is rounded up to the nearest whole number, with the number of digits being the same or greater than 1. Then, the maximum number of the crack element is determined based on the number of elements in the "same node number" array. For example, if the maximum node number of a solid element is 5730, then the initial value of the crack element number is E. LXmin The value is 6000; if the number of arrays with the same node number in the initial .inp file is 1680, it means that there are 1680 fracture elements, and the maximum value of the element number of the fracture element is 6000 + 1680 = 7680.

[0028] Copy the modified element numbers and node numbers to the corresponding locations in the new new.inp file. This allows you to embed crack elements at the boundaries of solid elements.

[0029] Specifically, step S4 includes the following:

[0030] S4.1: Set the analysis step;

[0031] S4.2: Set the displacement boundary for nodal force release in each analysis step;

[0032] In step S4.1, the content of the analysis step is set as follows:

[0033] Import the new.inp file into the finite discrete element main program. In the Step module, set the analysis step type to Dynamic Explicit. The number of analysis steps is determined as follows: if the total advancing distance of the working face is L and the advancing speed is l, then the maximum number of analysis steps is Step. max =L / l.

[0034] In step S4.2, the content for setting the displacement boundary for nodal force release in each analysis step is as follows:

[0035] In the Assembly module, define the set of nodes that need to release node forces for each analysis step, and the set of nodes that do not need to release node forces for that analysis step, denoted as NBoundary-i (1≤i≤Step). max (This represents the set number where node forces do not need to be released).

[0036] In the Load module, set the corresponding node set NBoundary-i in different analysis steps Step, fix the displacement in the y direction, and deactivate the displacement constraint in the (Step+1)th analysis step.

[0037] Furthermore, in step S6, the contents of the finite discrete element simulation program are:

[0038] The finite discrete element (FME) simulation program works as follows: Under initial boundary conditions and geostress conditions, it first uses a contact pair training algorithm, also known as a simulated through-crack training algorithm, to traverse and search for contact pairs throughout the entire numerical model. For the regions enclosed by the contact pairs, the FME method is used to calculate the nodal forces between the enclosed regions using Newton's second law. Then, the nodal forces at the boundary nodes are transferred to the interior of the regions, and the finite element method is used to calculate the elastoplastic and damage mechanical responses within the regions. Within the regions, the shared nodes of solid elements and crack elements transfer stress, and the finite element method is used to calculate the mechanical response of the crack elements, including stress and displacement. If the maximum displacement W of the numerically obtained crack element is greater than W... 实验 W 实验 To determine the maximum displacement of the crack obtained from the experiment, the element is deleted, thus creating a crack in the numerical model. Simultaneously, the contact pair is activated, and discrete element method (DEM) calculations are used to further perform numerical calculations of shear friction in the through-crack.

[0039] After the numerical calculation is completed, the nodal forces at the leading working face, the nodal forces at the support locations, and the nodal forces in the goaf behind the working face are output for each analysis step. This reveals the evolution law of the leading support pressure and the working resistance of the support during the advancement of the fully mechanized longwall face in an extra-thick coal seam.

[0040] The technological innovation of the finite discrete element simulation method for mine pressure manifestation in ultra-thick coal seam fully mechanized longwall mining faces of the present invention lies in:

[0041] 1. In fully mechanized longwall mining of extra-thick coal seams, the overburden rotation space is larger than in medium-thick coal seams as the top coal is released, and the range of involved strata is wider, making it easier to generate greater mining pressure. The finite discrete element simulation method of this invention can effectively simulate the overburden fracturing and rotation in fully mechanized longwall mining of extra-thick coal seams, output the nodal forces in this process, and obtain the manifestation law of mining pressure such as advance support pressure and support working resistance during mining.

[0042] 2. Releasing node forces step by step according to different analysis steps can effectively avoid writing subroutines and avoid data transfer between subroutines and the main program, thus effectively improving the calculation speed.

[0043] 3. It can be used to analyze the evolution of the advance support pressure and the nodal force at the rear support position of the working face in the fully mechanized longwall mining of extra-thick coal seams under the influence of different mining parameters and geological parameters. In this way, combined with the engineering geological conditions of the working face, the coal mining scheme can be optimized, and a basis can be provided for further formulating the roof (coal) grouting reinforcement scheme. Attached Figure Description

[0044] Figure 1 This is a schematic diagram of the finite discrete element numerical simulation process;

[0045] Figure 2 It is a numerical calculation model;

[0046] Figure 3 This is the 17th analysis step: Mises stress and mining-induced overburden fractures;

[0047] Figure 4 This is the nodal force in the 17th analysis step (vertical axis unit: ×10). 6 N);

[0048] Figure 5 This is step 32 of the analysis, which examines Mises stress and mining-induced overburden fractures.

[0049] Figure 6 This is the nodal force in analysis step 32 (vertical axis unit: ×10). 6 N).

[0050] The present invention will now be described in further detail with reference to the accompanying drawings and embodiments. Detailed Implementation

[0051] This embodiment presents a finite discrete element (FEMI) simulation method for mine pressure manifestation in fully mechanized longwall mining faces of extra-thick coal seams. The method establishes a numerical calculation model based on the dimensions of the working face under study; determines the solid element size of each rock stratum based on core drilling and measured data of coal and rock fractures underground, and then meshes the numerical model; generates zero-thickness fracture elements at the boundaries of solid elements by searching for the maximum node / element number and re-encoding the solid element number; sets a set of displacement boundaries for nodal force release and implements them sequentially in each analysis step; and numerically solves the nodal forces under different analysis steps using a finite discrete element simulation program, revealing the mine pressure manifestation laws such as advance support pressure and support working resistance during the mining of extra-thick coal seams in fully mechanized longwall mining faces. This method can not only numerically reproduce the rock plastic deformation and micro-fracture evolution and failure mechanism, but also effectively simulate the advance support pressure, support working resistance, and the quantitative relationship between overburden deformation / fracture / movement in fully mechanized longwall mining faces of extra-thick coal seams, providing a basis for further developing roof (coal) grouting reinforcement schemes and for the safe and efficient mining of extra-thick coal seams.

[0052] Specifically, the following steps are included:

[0053] S1: Establish a numerical calculation model based on the size of the study area;

[0054] The geometric data of the working face are determined as follows: length × width of the working face, represented by a × b; mining depth, represented by h; and mining thickness, represented by c. To avoid boundary effects, the model size in the horizontal direction is appropriately increased, with the added dimensions being d and e, thereby establishing a numerical model of length × width × height = (a + d) × (a + e) ​​× c.

[0055] Taking the Yushupo Coal Mine as an example, the 5107 working face has a mining depth of 380m, a mining thickness of 16m, and a working face advance distance of 1350m. Considering the large size of the calculation model and the periodic collapse of the overburden, the working face advance distance is shortened to 560m. Considering the boundary effect of the numerical model, a 100m boundary is designated as a retained coal pillar. A numerical calculation model with length × height = 760m × 380m is established.

[0056] S2: Determine the element size and mesh the numerical model;

[0057] S2.1: Determine the solid unit size of each rock layer

[0058] Core samples were obtained from each rock stratum using the borehole coring method. The rock strata were grouped according to their stratigraphy and lithology, and the average length of each core sample was calculated. At the same time, the joint and fracture spacing of the coal face in the underground tunnel was observed. Based on this, the dimensions of each rock stratum solid unit in the numerical model were determined to be f1, f2, f3, ...

[0059] S2.2: Mesh generation of the numerical model

[0060] In the Mesh module, the Seed Part mode is used to perform seeding operations on the numerical model established in S1; the seeding method is based on size, using the number f with the highest frequency in the size set {f1, f2, f3, ...} of each rock layer entity unit in step S2.1. i Based on this, the seed spacing of the overall numerical model is set; then, the Seed Edges mode is used to set the seed spacing of each rock layer, namely f1, f2, f3, ...; the shape of the control unit is controlled by Assign Mesh Controls, with hexahedrons being the preferred choice, followed by hexahedrons as the main shape, and finally tetrahedrons; then, the entire model is meshed using Mesh Part; under the Assign Element Type button, the meshed solid elements are modified to C3D8 solid element type.

[0061] In this embodiment, core samples are obtained from each rock stratum using a borehole coring method. The strata are grouped according to stratigraphy and lithology, and the average length of each core segment is calculated. Simultaneously, the joint and fracture spacing of the coal face and tunnels are observed underground. The fracture spacing for various rock mass types is determined considering computational cost. In the Mesh module, the Seed Part mode is used to seed the geometric model. The Seed Edges mode is used to set the spacing of the seeds for each rock stratum. Then, the numerical model is meshed, such as... Figure 2 As shown.

[0062] S3: Generate zero-thickness crack elements at the boundary of solid elements;

[0063] The process of generating zero-thickness crack elements at the boundary of solid elements is as follows:

[0064] (1) Find the maximum node number and cell number in step S2.

[0065] Set the working directory, create a new job under the Job module, and then locate the .inp script file of the newly created job in that working directory. In the .inp file, find the keyword node set *nset and the cell set *elset. Copy the numbers below the keywords into Excel, and use the MAX function to find the largest node number E of each entity cell. max and the largest unit number N max ;

[0066] (2) Re-encode the node number of the entity unit.

[0067] For any hexahedral or tetrahedral solid element, i.e., element number e j There are 8 nodes and 4 nodes respectively; unit e jThere are a total of g elements. For typical tetrahedral elements, g = 4; for hexahedral elements, g = 6. If element e... j The initial number of the i-th node is n. i That is, for hexahedral elements, 2≤i≤6, and for tetrahedral elements, 2≤i≤4, the node is copied g times; then the node numbers of these copied nodes are modified, except for the first element, the node numbers of the remaining elements are modified sequentially to (N max +n i After the above steps, the node number of the entity element can be updated, denoted as n. I The entity unit number remains unchanged and is still denoted as e. j ;

[0068] Add the keyword "*Elset=solid" after "*Part" and before "*Assembly" in the script file, and then write the node numbers of the updated entity cells in sequence; save it as a new file new.inp and save the above modified information.

[0069] (3) Generate fracture element numbers and node numbers

[0070] In the initial script file, arrays with different element numbers but three tetrahedral elements or four hexahedral elements sharing the same node number are copied, including both the node number and the element number; for two adjacent elements e with shared nodes... j1 and e j2 The shared node is denoted as n. i1 n i2 n i3 n i4 Therefore, the node number n after the root step S3 update is... I1 n I2 n I3 n I4 This refers to the node number of the crack element;

[0071] The element numbering of the crack element is determined as follows: the maximum node number of the solid element is E. max Therefore, the smallest element of the crack element is numbered E. LXmin The value is an integer that is rounded up to the nearest whole number and has the same or greater than 1 digit. Then, the maximum number of the crack element is determined based on the number of identical node numbers in the array. That is, if the maximum node number of the solid element is 5730, then the initial value of the crack element number is E. LXmin The value is 6000; if the number of arrays with the same node number in the initial .inp file is 1680, it means that there are 1680 fracture elements, and the maximum value of the element number of the fracture element is 6000 + 1680 = 7680.

[0072] Copy the modified element number and node number to the corresponding location in the new new.inp file, thereby embedding the crack element at the boundary of the solid element.

[0073] In the Mesh module, the predefined solid elements are modified to hexahedral C3D8 solid element type using the Assign Element Type setting. All solid elements are then grouped into a single element set named "solid". A script file named YSP.inp is generated in the Job module.

[0074] In this embodiment, a series of processes are performed, including finding the maximum node number and element number in the script file, copying nodes and re-encoding the node numbers of the entity elements, and generating the gap element number and node number, to embed gap elements at the boundary of the entity elements. The gap element type is modified to COH3D8 using the Assign Element Type function, and all gap elements are combined into an element set named "coh".

[0075] S4: Set up the analysis step, and set the set of displacement boundaries for nodal force release in each analysis step;

[0076] Includes the following:

[0077] S4.1, Set up the analysis step, the content of which is:

[0078] Import the new.inp file into the finite discrete element main program. Under the Step module, set the analysis step type to Dynamic Explicit. The number of analysis steps is determined as follows: If the total advancing distance of the working face is L and the advancing speed is l, then the maximum number of analysis steps is Step. max =L / l;

[0079] S4.2, Set the displacement boundary for nodal force release in each analysis step, the content of which is:

[0080] Within the Assembly module, define the set of nodes that need to release node forces for each analysis step, and the set of nodes that do not need to release node forces for that analysis step, denoted as NBoundary-i, i.e., 1≤i≤Step. max , represents the set number that does not require the release of node forces;

[0081] In the Load module, set the corresponding node set NBoundary-i in different analysis steps Step, fix the displacement in the y direction, and deactivate the displacement constraint in the (Step+1)th analysis step.

[0082] In this embodiment, the YSP.inp file is imported into the finite discrete element main program. First, in the Step module, the analysis step type is set to Dynamic Analytic. Based on the working face advancing speed of 10m / d and the working face advancing distance of 560m, 56 analysis steps are set.

[0083] In the Assembly module, starting from the cut point on the working face, a set of nodes that release nodal forces is defined every 10m, named JieDian01, JieDian02, ..., JieDian56. The set of nodes that do not require nodal force release is named NBoundary-01, NBoundary-02, ..., NBoundary-56.

[0084] In the Load module, in the first analysis step, the y-direction displacement of node set NBoundary-01 is set to 0, and this constraint is invalidated in the second and subsequent analysis steps; in the second analysis step, the y-direction displacement of node set NBoundary-02 is set to 0, and this constraint is invalidated in the third and subsequent analysis steps; and so on.

[0085] S5: Apply boundary conditions to the entire numerical model, input material compression and fracture mechanics parameters, and set output variables;

[0086] In this embodiment, chain constraints are applied to the boundaries of the numerical model to constrain the rigid body displacement of the model. The input parameters are the compressive mechanics parameters of the rock strata, including elastic modulus, Poisson's ratio, internal friction angle, and cohesion; and the input parameters are the material fracture mechanics parameters, including Mode I fracture strength, Mode II fracture strength, Mode I fracture displacement, Mode II fracture displacement, Mode I fracture energy, and Mode II fracture energy, as shown in Tables 1 to 3. Output variables are set, including stress components, strain components, displacement components, fracture energy, frictional stress, and frictional displacement.

[0087] Table 1: Mechanical parameters of solid elements

[0088] Elastic modulus / GPa Poisson's ratio Compressive 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

[0089] Table 2: Mechanical parameters of fracture elements

[0090]

[0091] Table 3: Contact Pair Mechanical Parameters

[0092] Shear modulus / GPa Peak shear strength / MPa Shear residual strength / MPa 7.3 12.1 6.7

[0093] S6: Numerical calculations are performed using a finite discrete element simulation program, and the nodal forces at the leading working face and the support positions behind the working face are output under different analysis steps.

[0094] The contents of the finite discrete element simulation program are:

[0095] The finite discrete element (FME) simulation program works as follows: Under initial boundary conditions and geostress conditions, it first uses a contact pair training algorithm, also known as a simulated through-crack training algorithm, to traverse and search for contact pairs throughout the entire numerical model. For the regions enclosed by the contact pairs, the discrete element method is used to calculate the nodal forces between the enclosed regions using Newton's second law. Then, the nodal forces at the boundary nodes are transferred to the interior of the regions, and the finite element method is used to calculate the elastoplastic and damage mechanical responses within the regions. Inside the regions, stress is transferred through shared nodes between solid elements and crack elements, and the mechanical response of the crack elements, including stress and displacement, is calculated using the finite element method. If the maximum displacement W of the crack element obtained numerically is greater than W0, the simulation is considered complete. 实验 W 实验 To obtain the maximum displacement of the crack from the experiment, the element is deleted, thus forming a crack in the numerical model; at the same time, the contact pair is activated, and the discrete element method is used to carry out the numerical calculation of shear friction of the through crack.

[0096] After the numerical calculation is completed, the nodal forces of the advanced working face, the nodal forces at the support positions, and the nodal forces of the goaf behind the working face are output for each analysis step; revealing the evolution law of the advanced support pressure and the working resistance of the support during the advancement of the fully mechanized longwall face in extra-thick coal seams.

[0097] In this embodiment, under the Job module, for the newly created job task YSP, the number of parallel computing cores is set to 52 in the Edit Job tab, and the computation is submitted.

[0098] S7: Based on nodal force data, reveal the manifestation patterns of mine pressures such as advance support pressure and support working resistance during the mining of extra-thick coal seams in fully mechanized longwall mining faces; on this basis, and in combination with the engineering geological conditions of the working face, formulate roof disaster prevention and control measures to achieve safe and efficient coal mining.

[0099] In this embodiment, numerical calculations are used to output nodal forces, revealing the manifestation patterns of mine pressures such as advance support pressure and support working resistance during the mining of extra-thick coal seams in fully mechanized longwall faces. The analysis results show that the distribution and evolution patterns of mining-induced roof fractures, such as... Figure 3 and Figure 5 The black lines indicate this. Specifically, the advanced support pressure is extracted from the nodes of the coal wall in the advanced working face; while the working resistance on the hydraulic support is extracted from the nodes of the coal wall in the lagging working face, such as... Figure 4 and Figure 6 As shown.

[0100] contrast Figure 4 and Figure 6It can be seen that when the working face advances to the length of the working face (210m), the key layer in the overburden fractures, and the pressure on the advance support increases dramatically, from 40×10 6 N rapidly increased to 80 × 10 6 N causes coal face fracturing at the working face, resulting in roof collapse and other roof disasters. Numerical results of mine pressure indicate that when the working face advances to this location, measures such as roof grouting and basic roof decompression should be taken in advance.

[0101] The finite discrete element method (FEMI) simulation of mine pressure manifestation in fully mechanized longwall mining faces of extra-thick coal seams presented in this embodiment, compared with existing numerical methods for mine pressure manifestation, not only realizes the process of roof rock transforming from a continuous medium to a discrete medium under mining-induced stress, but more importantly, it quantitatively describes the relationship between advance support pressure, support working resistance, and overburden deformation / fracture / movement under fully mechanized longwall mining conditions of extra-thick coal seams. The study of the evolution of advance support pressure and support working resistance with the analysis step (or advance distance) provides a powerful numerical simulation tool for formulating roof disaster prevention and control measures for fully mechanized longwall mining faces of extra-thick coal seams.

[0102] 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 additions, improvements, modifications, or equivalent changes made by those skilled in the art without departing from the scope of the technical solution of the present invention shall fall within the protection scope defined by the claims of the present invention.

Claims

1. A finite discrete element simulation method for mine pressure manifestation in fully mechanized longwall mining faces of extra-thick coal seams, characterized in that, Includes the following steps: S1: Establish a numerical calculation model based on the size of the study area; S2: Determine the element size and mesh the numerical model; this includes the following: S2.1: Determine the solid unit size of each rock layer Core samples were obtained from each rock stratum using the borehole coring method. The rock strata were grouped according to their stratigraphy and lithology, and the average length of each core sample was calculated. At the same time, the joint and fracture spacing of the coal face in the underground tunnel was observed. Based on this, the dimensions of each rock stratum solid unit in the numerical model were determined to be f1, f2, f3, ... S2.2: Mesh generation of the numerical model In the Mesh module, the Seed Part mode is used to perform seeding operations on the numerical model established in S1; the seeding method is based on size, using the number f with the highest frequency in the size set {f1, f2, f3, ...} of each rock layer entity unit in step S2.

1. i Based on this, the seed spacing for the entire numerical model is set; Then, using the Seed Edges mode, the spacing of the seeds for each rock layer is set to f1, f2, f3, ...; the shape of the control unit is controlled by Assign Mesh Controls, prioritizing hexahedrons, followed by hexahedrons as the main shape, and finally tetrahedrons; then, the entire model is meshed using Mesh Part; under the Assign Element Type button, the meshed solid elements are modified to C3D8 solid element type; S3: Generate zero-thickness crack elements at the boundary of solid elements; S4: Set up the analysis step, and set the set of displacement boundaries for nodal force release in each analysis step; S5: Apply boundary conditions to the entire numerical model, input material compression and fracture mechanics parameters, and set output variables; S6: Numerical calculations are performed using a finite discrete element simulation program, and the nodal forces at the leading working face and the support positions behind the working face are output under different analysis steps. The process of finite discrete element simulation is as follows: Under the initial boundary conditions and geostress conditions, the contact pair training algorithm, also known as the simulated through-crack training algorithm, is first used to traverse and search for contact pairs throughout the entire numerical model. The discrete element method is used between the regions enclosed by the contact pairs to calculate the nodal forces between the enclosed regions using Newton's second law. Then, the nodal forces of the boundary nodes are transferred to the interior of the region, and the elastoplastic and damage mechanical responses inside the region are calculated using the finite element method. Within the region, stress is transferred through shared nodes of solid elements and crack elements. The mechanical response of the crack elements, including stress and displacement, is calculated using the finite element method. If the maximum displacement W of the fracture element obtained from the numerical analysis is greater than W 实验 W 实验 To determine the maximum displacement of the crack obtained from the experiment, the element is deleted, thus forming a crack in the numerical model. At the same time, the contact pair is activated, and the discrete element method is used to carry out numerical calculations of shear friction of the through crack. After the numerical calculation is completed, the nodal forces of the working face ahead, the nodal forces at the support positions, and the nodal forces of the goaf behind the working face are output for each analysis step. To reveal the evolution law of advance support pressure and support working resistance during the advancement of fully mechanized longwall mining faces in extra-thick coal seams; S7: Based on the nodal force data, reveal the advance support pressure and support working resistance during the mining of extra-thick coal seams in fully mechanized longwall mining. On this basis, combined with the engineering geological conditions of the working face, formulate roof disaster prevention and control measures to achieve safe and efficient coal mining.

2. The method as described in claim 1, characterized in that, In step S1, the geometric data of the working face are determined as follows: the length × width of the working face, represented by a × b, the mining depth, represented by h, and the mining thickness, represented by c. To avoid boundary effects, the model size in the horizontal direction is appropriately increased, with the added dimensions being d and e, thereby establishing a numerical model of length × width × height = (a + d) × (a + e) ​​× c.

3. The method as described in claim 1, characterized in that, In step S3, the process of generating zero-thickness crack elements at the boundary of solid elements is as follows: (1) Find the maximum node number and cell number in step S2. Set the working directory, create a new job under the Job module, and then find the .inp script file of the newly created job in that working directory; Find the set of keyword nodes in the .inp file. nset and unit set else; Copy the numbers below the above keywords into Excel, and use the MAX function to find the largest node number E of each entity cell. max and the largest unit number N max ; (2) Re-encode the node number of the entity unit For any hexahedral solid element or tetrahedral solid element, i.e., element numbered as e j There are 8 nodes and 4 nodes respectively; unit e j There are a total of g elements. For typical tetrahedral elements, g=4; for hexahedral elements, g=6. If the elements... e j The i The initial number of each node is n i That is, in a hexahedral element, 2≤ i ≤6, 2≤ in tetrahedral elements i ≤4, copy the node g Then, the node numbers of these copied nodes are modified. Except for the first unit, the node numbers of the remaining units are modified sequentially to (N... max + n i After the above steps, the node number of the entity element can be updated, denoted as . n I The entity unit number remains unchanged and is still denoted as... e j ; The keyword "in the script file" After "Part", Add keywords before "Assembly" Then, set Elset=solid, and then sequentially write the modified cell numbers to obtain the updated node numbers of the solid cells; save it as a new file new.inp, and save the above modified information. (3) Generate fracture element numbers and node numbers In the initial script file, arrays with different element numbers but three tetrahedral elements or four hexahedral elements sharing the same node number are copied, including both the node number and the element number; for two adjacent elements sharing a node... e j1 and e j2 The shared node number is denoted as n i1 , n i2 , n i3 , n i4 Therefore, the node number updated in root step S3 is... n I1 , n I2 , n I3 , n I4 This refers to the node number of the crack element; The element numbering of the crack element is determined as follows: the maximum node number of the solid element is E. max Therefore, the smallest element of the crack element is numbered E. LXmin The value is an integer that is rounded up without rounding down, and the number of digits is the same or greater than 1; then the maximum number of the crack unit is determined based on the number of "same node numbers" in the array; That is: if the maximum node number of the solid element is 5730, then the initial value of the crack element number E is... LXmin The value is 6000; if the number of arrays with the same node number in the initial .inp file is 1680, it means that there are 1680 fracture elements, and the maximum value of the element number of the fracture element is 6000 + 1680 = 7680. Copy the modified element number and node number to the corresponding location in the new new.inp file, thereby embedding the crack element at the boundary of the solid element.

4. The method as described in claim 1, characterized in that, Step S4 includes the following: S4.1, Set up the analysis step, the content of which is: Import the new.inp file into the finite discrete element main program. Under the Step module, set the analysis step type to Dynamic Explicit. The number of analysis steps is determined as follows: If the total advancing distance of the working face is L and the advancing speed is l, then the maximum number of analysis steps is Step. max =L / l; S4.2, Set the displacement boundary for nodal force release in each analysis step, the content of which is: Within the Assembly module, define the set of nodes that need to release node forces for each analysis step, and the set of nodes that do not need to release node forces for that analysis step, denoted as NBoundary-i, i.e., 1≤i≤Step. max , represents the set number that does not require the release of node forces; In the Load module, set the corresponding node set NBoundary-i in different analysis steps Step, fix the displacement in the y direction, and deactivate the displacement constraint in the (Step+1)th analysis step.

Images