Three-dimensional simulation method for coal mining based on numerical flow and structured grid

Abstract

This invention discloses a three-dimensional simulation method for coal mine mining based on numerical manifolds and structured meshes, belonging to the field of three-dimensional simulation of coal mine mining. The method includes determining the initial position of the working face and collecting a field database; meshing the computational domain and fusing the field database to obtain a mesh model incorporating field data; constructing and coupling the numerical manifold coverage with the mesh model incorporating field data; coupling the mechanical behavior of fractures with seepage; solving the three-dimensional stress field, displacement field, and pore gas pressure field of the current computational domain over time based on a CPU-GPU heterogeneous parallel structure; and advancing the coal mine working face to the next stage until the entire working face is advanced, obtaining the stress distribution changes, fluid transport changes, and fracture propagation changes throughout the working face advancement process. This invention solves the problems of existing methods in accurately describing the stress accumulation of initial micro-fractures and maintaining contact / separation topology updates after fracture.

Description

Technical Field

[0001] This invention belongs to the field of three-dimensional simulation of coal mining, and particularly relates to a three-dimensional simulation method for coal mining based on numerical manifolds and structured grids. Background Technology

[0002] During the advancement of longwall coal face, the surrounding rock undergoes continuous-discontinuous multi-field coupled evolution, which can easily induce dynamic disasters such as roof instability and rockburst. To reveal the three-dimensional stress wave propagation, fracture propagation, and dynamic load dissipation mechanisms during mining, the engineering community commonly employs finite element method (FEM), discrete element method (DEM), or explicit / implicit coupling methods for numerical simulation. However, in high-stress, fractured coal and rock masses, existing methods generally suffer from the following bottlenecks: (1) Traditional structured meshes are difficult to characterize complex fracture networks. Multi-scale bedding, faults and induced fractures are often distributed around the coal mine working face. Traditional FEM relies on single-domain meshes shared by continuous nodes. If fractures are directly embedded in a three-dimensional structured mesh, the mesh needs to be subdivided or reconstructed locally. The element distortion and computational scale increase sharply, resulting in "difficulty in three-dimensional structured modeling under the condition of multiple fracture surfaces and difficulty in controlling element geometric resolution; (2) Dual-medium fluid-structure interaction has low solution efficiency and poor stability. The unloading-reloading caused by the advancement of the working face is accompanied by seepage and thermo-chemical effects. Pure FEM or DEM requires a very small time step to ensure convergence when dealing with strongly nonlinear problems such as fracture seepage and block slip, resulting in "complex numerical solution of fluid-structure interaction with large computational load; (3) Insufficient capture of continuous-discontinuous transformation. During the cutting of the working face, the rock mass gradually enters the block movement stage from a continuous body. Although explicit DEM can handle discrete blocks, it is difficult to describe the stress accumulation of the initial microcracks with high accuracy; while traditional FEM is difficult to maintain the topological update of contact / separation after fracturing. Summary of the Invention

[0003] To address the aforementioned shortcomings in existing technologies, this invention provides a three-dimensional simulation method for coal mining based on numerical manifolds and structured meshes. This method solves the problems of existing methods being unable to accurately describe the stress accumulation in the initial microcracks and the difficulty in maintaining contact / separation during topological updates after fracture.

[0004] To achieve the aforementioned objectives, the technical solution adopted by this invention is: a three-dimensional simulation method for coal mining based on numerical manifolds and structured grids, comprising: Determine the initial location of the coal mine working face, collect on-site data of the coal mine working face, unify the coordinate system of the on-site data of the coal mine working face into a single HDF5 file, and obtain the on-site database; The fine mesh range is defined, and the goaf, the coal mine working face, and the fine mesh range are used as the calculation area; the fine mesh range is set in the direction of coal mine working face advancement and is adjacent to the coal mine working face. The computational domain is divided into grids and integrated with the field database to obtain a grid model that incorporates field data. Construct a numerical manifold cover and couple the numerical manifold cover with a mesh model that integrates field data; Coverage data is obtained based on the coupling results of numerical manifold cover and the mesh model fused with field data; Based on the coverage data, the three-dimensional stress field, displacement field and pore gas pressure field of the current computing region are solved by multi-field coupling of CPU-GPU heterogeneous parallel structure; The cubic grid containing the mined coal body is marked as the mined-out state, and the weak material parameters of the collapse zone are used to represent it. The next step is to advance the coal mine working face until the entire working face is advanced, and then obtain the stress distribution changes, fluid migration changes, and crack propagation changes throughout the entire working face advancement process.

[0005] Furthermore, the field database includes a list of mechanical parameters for each rock stratum, a list of seepage parameters for each rock stratum, a three-dimensional fracture surface file of the working face, a pore water pressure field, a pore permeability coefficient field, and a gas seepage parameter table.

[0006] Furthermore, the process of obtaining the mesh model that integrates field data specifically involves: dividing the computational region into several cubic grids based on the rock stratum thickness, further subdividing the cubic grids within the fine mesh range, and removing cubic grids with deformation greater than the deformation threshold to obtain a multi-scale cubic grid file; mapping the data from the field database to the grid nodes in the multi-scale cubic grid file using a unified coordinate transformation matrix to obtain the mesh model that integrates field data.

[0007] Furthermore, the construction of the numerical manifold cover specifically involves: extracting the fracture centerlines of each fracture surface in the mesh model, and arranging several layers of cover elements on both sides of the fracture centerlines to describe the opening and closing of fractures and the movement of the block; the coverage influence radius of the cover element is:

[0008]

[0009] in, The coverage influence radius of the coverage unit; This is an empirical coefficient; The elastic modulus of the rock; This is the critical value for the fracture energy release rate; Tensile strength of rock; The maximum stress difference between the elements before and after subdivision; The reference maximum stress.

[0010] Furthermore, the covering unit specifically comprises several virtual nodes arranged along the normal direction on both sides of the crack centerline, with a crack constitutive model set between each virtual node; The covering unit is sparsely mapped and coupled with the grid nodes of the grid model that integrates field data through a numerical manifold method to simulate the crack opening and closing process and block movement. The displacement degrees of freedom of the covering unit are independent of the mesh model that integrates field data.

[0011] Furthermore, the numerical manifold cover is coupled with the mesh model fused with field data, specifically by mapping the displacement information of the cover cells onto the mesh nodes of the mesh model fused with field data:

[0012] in, It is a sparse mapping matrix; The displacement vector of the grid nodes; This is the displacement vector of the covering element node.

[0013] Furthermore, the overlay data includes a normal opening / closing model to describe the stiffness change during crack opening and closing, a shear slip model to describe the relationship between shear force and slip displacement on the crack surface, and crack permeability:

[0014]

[0015]

[0016] in, The normal stress in the crack; The initial normal stiffness; The crack opening; This is the initial closing opening. This refers to the residual normal stress; Shear stress; It is the friction angle; It is cohesive force; This represents the current crack permeability. Initial penetration rate; The penetration-aperture coupling coefficient; This refers to the crack opening.

[0017] Furthermore, the method of solving the three-dimensional stress field, displacement field, and pore gas pressure field of the current computational region based on CPU-GPU heterogeneous parallel structure multi-field coupling specifically involves: On the GPU, the central difference method is used to explicitly integrate the structured mesh nodes and update the displacement field:

[0018]

[0019]

[0020]

[0021] in, For the displacement increment of the cube mesh; Density; For time step; For gradient operators; For stress tensor; This is the stress acceleration vector; Minimum unit size; It is the elastic wave velocity; The elastic modulus of the rock; Here is the tangent stiffness matrix; The displacement of the covering element; It is a nonlinear residue; On the GPU, the stress field is calculated using an elastic constitutive model based on the strain tensor. At the same time, the local stress within the cover element is updated by combining the normal opening and closing model and the shear slip model of the crack region. The gas pressure field is implicitly solved on the CPU based on a dual-medium seepage model, which incorporates permeability controlled by fracture aperture and gas viscosity corrected by the Klinkenberg slip effect.

[0022]

[0023] in, The volumetric flow velocity density in the fracture; This represents the current crack permeability. This represents the pressure gradient within the fracture. For gradient operators; This refers to the fracture pressure; The corrected gas viscosity; Reference viscosity; The Klinkenberg coefficient; The three-dimensional stress field, displacement field, and pore gas pressure field are iteratively solved on both GPU and CPU based on the coupling of fracture aperture, stress state, and seepage parameters. After each iteration at time step, the overall residual and energy conservation error are calculated. If the overall residual is lower than the set threshold and the energy conservation error is less than the error threshold, a three-dimensional stress field, displacement field and pore gas pressure field that evolve over time are formed; otherwise, the adaptive evolution process is entered.

[0024] Furthermore, the overall residual is:

[0025] in, For the overall residual; Here is the tangent stiffness matrix; The displacement of the covering element; This is the external load vector; The energy conservation error for:

[0026] in, The internal energy of the system; The kinetic energy of the system; Work is done by external forces acting on the system.

[0027] Furthermore, the adaptive evolution process specifically includes: Calculate the current energy release rate of the fracture:

[0028] in, Energy release rate; For strain energy; For external skills; The area of ​​the crack; When the current crack energy release rate is greater than the release rate threshold At that time, crack propagation is triggered; the crack extension amount of the crack propagation is:

[0029] in, Extending for a single-step forward; For control coefficients; The side length of the fine mesh unit; Based on the amount of crack extension, a covering unit is automatically generated at the new crack location, and the mapping relationship is updated simultaneously. Update crack permeability:

[0030] in, The updated crack permeability; The permeability of the cracks before the update; The penetration-aperture coupling coefficient; The crack opening; The multi-field coupling calculation is then re-performed in the updated mesh and crack field.

[0031] The beneficial effects of this invention are as follows: This invention seamlessly couples regular cubic meshes with numerical manifold overlay elements. The continuous medium still employs an efficient structured mesh, where cracks and block movements only appear locally in overlay form, and the two maintain consistent displacement and stress through a mapping matrix. This avoids element distortion caused by forcibly inserting cracks into the mesh while preserving space for free crack growth. Practical examples show that, within the same cross-sectional area, the total mesh size can be reduced by more than 40% compared to pure finite element methods, and the computation time is compressed to one-third of the original. Simultaneously, the error between crack initiation and extension locations and the on-site microseismic positioning is controlled within a few tens of meters. For the dynamic process of "advancement-unloading" at the working face, this method proposes a modeling strategy of "local fine mesh at the leading edge + real-time deletion / activation of the goaf area." The leading edge region is automatically subdivided to ensure spatial resolution of the unloading wave process; while the mined-out area is immediately deleted and replaced with collapsed material, avoiding memory explosion caused by global fine mesh deployment. Extensive testing showed that, at the same level of precision, memory usage was only 30-40% of that of traditional global fine mesh solutions, and the wall clock time for a single advance step remained stable within one to two hours, meeting shift-level decision-making requirements. This invention employs a heterogeneous parallel framework of "GPU explicit integration + CPU implicit iteration." Linear elastic-plastic calculations on the structured mesh are entirely processed by the parallel cores of the graphics card, while the strongly nonlinear equations within the coverage domain are solved using Newton-Krylov iteration on a multi-core CPU. The two types of results exchange residuals in sub-iterations until convergence. Compared to a pure CPU solver, this division of labor can further shorten the solution time by approximately three-quarters and ensure numerical stability during the violent crack opening phase. Crack propagation uses an energy release rate criterion, and incremental coverage is generated instantaneously at the propagation location. Whether a crack continues to grow no longer depends on the local instability of mesh elements, but rather on whether the energy pulled apart by the rock exceeds a critical value. Therefore, the crack extension direction is not limited by the mesh axis, avoiding the problem of "cracks running at right angles" in traditional finite element methods. Compared with the true triaxial loading test in the laboratory, the simulated fracture dip angle error is less than five degrees, and the separation path of the roof slab is smoother and more reliable. Regarding fluid coupling, the model simultaneously solves for two pressure fields: the matrix pores and the fracture channels, and applies Klinkenberg correction to the viscosity of low-pressure gas flow. This improvement significantly enhances the prediction accuracy of gas-rich areas: in typical mine tests, the spatial positioning accuracy of high-concentration gas areas has increased from 70% to over 90%, and can be directly used to guide the layout of drainage boreholes. From the outset, this invention uses field-measured data from four major categories—drilling, geophysical exploration, hydrology, and gas—to drive mesh and overlay generation, reducing a significant amount of secondary calibration work relying on empirical parameters. The initial stress and pore pressure deviations from the measured values ​​are controlled within 10%, laying a reliable foundation for subsequent dynamic evolution calculations. Attached Figure Description

[0032] Figure 1 This is a flowchart of the method of the present invention. Detailed Implementation

[0033] The specific embodiments of the present invention are described below to enable those skilled in the art to understand the present invention. However, it should be understood that the present invention is not limited to the scope of the specific embodiments. For those skilled in the art, various changes are obvious as long as they are within the spirit and scope of the present invention as defined and determined by the appended claims. All inventions utilizing the concept of the present invention are protected.

[0034] like Figure 1 As shown, in one embodiment of the present invention, a three-dimensional simulation method for coal mining based on numerical manifolds and structured grids includes: Determine the initial location of the coal mine working face, collect on-site data of the coal mine working face, unify the coordinate system of the on-site data of the coal mine working face into a single HDF5 file, and obtain the on-site database; The fine mesh range is defined, and the goaf, the coal mine working face, and the fine mesh range are used as the calculation area; the fine mesh range is set in the direction of coal mine working face advancement and is adjacent to the coal mine working face. The computational domain is divided into grids and integrated with the field database to obtain a grid model that incorporates field data. Construct a numerical manifold cover and couple the numerical manifold cover with a mesh model that integrates field data; Coverage data is obtained based on the coupling results of numerical manifold cover and the mesh model fused with field data; Based on the coverage data, the three-dimensional stress field, displacement field and pore gas pressure field of the current computing region are solved by multi-field coupling of CPU-GPU heterogeneous parallel structure; The cubic grid containing the mined coal body is marked as the mined-out state, and the weak material parameters of the collapse zone are used to represent it. The next step is to advance the coal mine working face until the entire working face is advanced, and then obtain the stress distribution changes, fluid migration changes, and crack propagation changes throughout the entire working face advancement process.

[0035] The field database includes a list of mechanical parameters for each rock stratum, a list of seepage parameters for each rock stratum, a three-dimensional fracture surface file of the working face, a pore water pressure field, a pore permeability coefficient field, and a gas seepage parameter table.

[0036] In this embodiment, drilling data is collected. Experimental data such as formation thickness, core photographs, and rock strength and permeability are obtained through advance drilling or tunnel boreholes. An initial list of mechanical and seepage parameters for each rock stratum is established. Geophysical data collection. Collection of geophysical data including transient electromagnetic, seismic wave, and ground-penetrating radar data. Generation of 3D fracture / fault surface files for the working face; Hydrogeological and gas geological data collection. Pumping or injection tests to determine aquifer pressure and permeability coefficient, obtaining the original pore water pressure field and permeability coefficient field. Measurement of gas content, gas pressure, and adsorption constant, establishing a gas seepage parameter table; Data integration. A unified coordinate system is established, and multi-source data is integrated into a single HDF5 file. A standardized field database is generated for use in subsequent steps.

[0037] The process of obtaining the mesh model that integrates field data is as follows: the computational region is divided into several cubic grids according to the rock layer thickness, and the cubic grids within the fine mesh range are further subdivided. Cube grids with deformation greater than the deformation threshold are removed to obtain a multi-scale cubic grid file. The data in the field database is mapped to the grid nodes in the multi-scale cubic grid file through a unified coordinate transformation matrix to obtain the mesh model that integrates field data.

[0038] In this embodiment, the generation of the three-dimensional structured mesh on the working surface is specifically as follows: Determine the calculation scope: cover the goaf, the fine mesh area in front, and the affected surrounding rock within the stress wave range.

[0039] Generate the parent grid: Divide the entire area into cubic (hexahedral) grids according to the rock layer thickness; the grid side length is taken as the smaller layer thickness or about two meters to ensure accuracy.

[0040] Front-end subdivision: Within 30 to 40 meters in front of the working face, the cube is further subdivided into smaller grids of one-quarter to one-eighth the size of the parent grid, making the unloading wave simulation more refined.

[0041] Quality inspection: Remove highly deformed grids, ensure that each grid has a good shape and is not severely stretched or twisted, and obtain a qualified multi-scale cubic mesh file.

[0042] Multi-source data is automatically mapped to grid / cover nodes through a unified coordinate transformation matrix, with a spatial error of ≤0.2m, avoiding errors from manual mapping.

[0043] The construction of the numerical manifold cover specifically involves: extracting the fracture centerlines of each fracture surface in the mesh model, and arranging several layers of cover elements on both sides of the fracture centerlines to describe the opening and closing of fractures and the movement of the block; the coverage influence radius of the cover element is:

[0044]

[0045] in, The coverage influence radius of the coverage unit; This is an empirical coefficient; The elastic modulus of the rock; This is the critical value for the fracture energy release rate; Tensile strength of rock; The maximum stress difference between the elements before and after subdivision; The reference maximum stress.

[0046] In this embodiment, 4–8 layers of "covering units" are arranged on both sides of each crack. These act like meshless patches, freely describing crack opening and closing and block movement on the structural mesh. The influence radius of the cover is adaptively set according to the rock toughness and tensile strength. By associating the nodes of the covering units and the structural mesh, it is ensured that displacement and force can be mutually transferred between the two discrete systems.

[0047] In this embodiment, the "affected range (the radius of influence of the covering unit)" is determined based on the balance between the attenuation of disturbance energy or stress propagation in the rock mass and the rock mass's resistance to damage (strength, toughness). This is determined by rock sample tests, and the error in the crack tip strength factor can be controlled within 5%.

[0048] The covering unit is specifically a number of virtual nodes arranged along the normal direction on both sides of the crack centerline, and a crack constitutive model is set between each virtual node. The covering unit is sparsely mapped and coupled with the grid nodes of the grid model that integrates field data through a numerical manifold method to simulate the crack opening and closing process and block movement. The displacement degrees of freedom of the covering unit are independent of the mesh model that integrates field data.

[0049] In this embodiment, the covering element describes the normal opening and closing of the crack and shear slip by arranging multiple virtual nodes along the normal direction on both sides of the crack centerline and introducing a dedicated crack constitutive model between these nodes. Its displacement degree of freedom is independent of the background mesh, which can capture the relative displacement and opening changes between the blocks on both sides of the crack. It is coupled with the main mesh nodes through sparse mapping using the numerical manifold method, thereby achieving high-precision simulation of the crack opening and closing process and block motion without modifying the original structural mesh.

[0050] The coupled numerical manifold cover and the fused field data mesh model specifically involve mapping the displacement information of the cover cells onto the mesh nodes of the fused field data mesh model:

[0051] in, It is a sparse mapping matrix; The displacement vector of the grid nodes; This is the displacement vector of the covering element node.

[0052] In this embodiment, the displacement information on the main region (cover) is mapped to the computational grid nodes through interpolation, projection, or constraint rules.

[0053] The overlay data includes a normal opening / closing model to describe the stiffness change during crack opening and closing, a shear slip model to describe the relationship between shear force and slip displacement on the crack surface, and crack permeability:

[0054]

[0055]

[0056] in, The normal stress in the crack; The initial normal stiffness; The crack opening; This is the initial closing opening. This refers to the residual normal stress; Shear stress; It is the friction angle; It is cohesive force; This represents the current crack permeability. Initial penetration rate; The penetration-aperture coupling coefficient; This refers to the crack opening.

[0057] In this embodiment, the nonlinear mechanical response of a crack during compressive closure or tensile opening is simulated by setting the normal stiffness of the crack to vary with its opening and closing states. The relationship between shear force and slip displacement on the crack surface is described using parameters such as the friction coefficient, bond strength, and residual strength, reflecting the bond-slip mechanism before and after crack failure. The initial crack permeability is set to increase exponentially with crack opening. The permeability-stress coupling formula expresses the nonlinear response relationship where permeability decreases exponentially with increasing normal stress and increases rapidly with crack opening, through the stress-controlled change in crack opening. Fitted by borehole transient test.

[0058] The method for solving the three-dimensional stress field, displacement field, and pore gas pressure field of the current computational region based on CPU-GPU heterogeneous parallel structure multi-field coupling is as follows: On the GPU, the central difference method is used to explicitly integrate the structured mesh nodes and update the displacement field:

[0059]

[0060]

[0061]

[0062] in, For the displacement increment of the cube mesh; Density; For time step; For gradient operators; For stress tensor; This is the stress acceleration vector; Minimum unit size; It is the elastic wave velocity; The elastic modulus of the rock; Here is the tangent stiffness matrix; The displacement of the covering element; It is a nonlinear residue; On the GPU, the stress field is calculated using an elastic constitutive model based on the strain tensor. At the same time, the local stress within the cover element is updated by combining the normal opening and closing model and the shear slip model of the crack region. The gas pressure field is implicitly solved on the CPU based on a dual-medium seepage model, which incorporates permeability controlled by fracture aperture and gas viscosity corrected by the Klinkenberg slip effect.

[0063]

[0064] in, The volumetric flow velocity density in the fracture; This represents the current crack permeability. This represents the pressure gradient within the fracture. For gradient operators; This refers to the fracture pressure; The corrected gas viscosity; Reference viscosity; The Klinkenberg coefficient; The three-dimensional stress field, displacement field, and pore gas pressure field are iteratively solved on both GPU and CPU based on the coupling of fracture aperture, stress state, and seepage parameters. After each iteration at time step, the overall residual and energy conservation error are calculated. If the overall residual is lower than the set threshold and the energy conservation error is less than the error threshold, a three-dimensional stress field, displacement field and pore gas pressure field that evolve over time are formed; otherwise, the adaptive evolution process is entered.

[0065] In this embodiment, the explicit central difference method updates the displacement and stress of the cubic mesh on the graphics card. The particle acceleration is obtained by dividing the nodal force by the mass, and explicit integration is performed using a time step to update displacement and velocity point by point, achieving efficient and parallel explicit dynamic calculations. The time step must satisfy the CFL condition to ensure that wave velocity propagation does not exceed the mesh size within a single time step, thereby ensuring numerical stability and accurate capture of physical causality. The implicit iterative formula solves for unknowns simultaneously within the current time step, making the numerical solution of the system stable and capable of handling larger time steps.

[0066] In this embodiment, the dual-medium flow rate formula describes the bidirectional exchange process of fluid between microscopic pores and macroscopic fractures by separately solving the fluid pressure fields in the pores and fracture channels of the coal and rock mass, and coupling the two with a leakage term. The gas viscosity correction formula is based on the Klinkenberg effect, and adjusts the gas viscosity to more accurately reflect the actual seepage characteristics by considering the slip flow phenomenon of gas under low pressure.

[0067] In this embodiment, a multi-field coupling solution based on CPU-GPU heterogeneous parallelism is used: if the overall residual is lower than a set threshold and the energy conservation error is less than the error threshold, the three-dimensional stress field, displacement field, and pore gas pressure field evolving over time are obtained. In this method, the evolution of the three-dimensional stress field, displacement field, and pore gas pressure field over time is realized through multi-field coupling iteration: on the GPU, the central difference method is used to explicitly integrate the structured mesh nodes to update the displacement field; the stress field is calculated based on the strain tensor using an elastic or elastoplastic constitutive model, while the local stress within the covering unit is updated by combining the normal opening and closing model and the shear slip model of the crack region; on the CPU, the gas pressure field is implicitly solved based on a dual-medium seepage model, and the permeability controlled by the crack aperture and the viscosity corrected by Klinkenberg are introduced; the three are solved iteratively by coupling the crack aperture, stress state, and seepage parameters until the residual and energy conservation error meet the convergence conditions, forming a three-dimensional stress, displacement, and gas pressure field evolving over time. Otherwise, an adaptive evolution process is initiated: First, the current crack energy release rate is calculated. If it exceeds the critical threshold, crack propagation is triggered, new covering cells are generated along the crack leading edge, and the crack permeability is updated. At the same time, the numerical manifold covering structure and node mapping relationship are adjusted according to the distribution of the latest covering cells. Subsequently, multi-field coupling calculations are re-executed in the updated mesh and crack field, and iteration continues until the residual and energy conservation conditions are met, ensuring that the evolution of the three-dimensional stress field, displacement field, and pore gas pressure field reflects the physical response of crack propagation and fluid transport.

[0068] The overall residual is:

[0069] in, For the overall residual; Here is the tangent stiffness matrix; The displacement of the covering element; This is the external load vector; The energy conservation error for:

[0070] in, The internal energy of the system; The kinetic energy of the system; Work is done by external forces acting on the system.

[0071] In this embodiment, the overall residual formula measures the magnitude of the inconsistency between the current iterative solution and the equilibrium equation, reflects the error of the system's mechanical and fluid conservation equations, and thus determines the degree of convergence of the numerical solution.

[0072] The adaptive evolution process is specifically as follows: Calculate the current energy release rate of the fracture:

[0073] in, Energy release rate; For strain energy; For external skills; The area of ​​the crack; When the current crack energy release rate is greater than the release rate threshold At that time, crack propagation is triggered; the crack extension amount of the crack propagation is:

[0074] in, Extending for a single-step forward; For control coefficients; The side length of the fine mesh unit; Based on the amount of crack extension, a covering unit is automatically generated at the new crack location, and the mapping relationship is updated simultaneously. Update crack permeability:

[0075] in, The updated crack permeability; The permeability of the cracks before the update; The penetration-aperture coupling coefficient; The crack opening; The multi-field coupling calculation is then re-performed in the updated mesh and crack field.

[0076] In this embodiment, the energy release rate criterion calculates the elastic energy released per unit crack extension at the crack tip and compares it with the material's critical fracture energy to determine whether the crack will continue to propagate. The crack extension formula determines the advance distance of the crack front proportionally based on the magnitude of the energy release rate exceeding the critical value, achieving adaptive adjustment of the crack extension step size.

[0077] The fracture front extends forward proportionally by one or several mesh lengths each time, the specific extension depending on the degree to which the energy release rate exceeds the critical value. A patch is automatically generated at the new fracture location, simultaneously updating the mapping relationship. New cover nodes inherit from the front edge, adaptively updating the mapping matrix to ensure that the overall bandwidth increase is less than 3%. Once a fracture opens, permeability increases exponentially; it decreases when the fracture closes. The permeability update formula adjusts permeability through the exponential change in fracture opening, achieving a dynamic response where permeability increases rapidly when the fracture opens and decreases rapidly when it closes.

[0078] In this embodiment, the working face advancement segmentation and mesh activation are specifically as follows: Loading and advancing plan: Calculate the total number of advancing steps based on the total length of the working face and the step distance ΔL.

[0079] Goaf Reduction: After each step is completed, the grid containing the mined coal body is deleted and replaced with the weak material parameters of the collapsed area.

[0080] Fine mesh activation: Simultaneously activate the fine mesh area at the leading edge of the next step, and inherit the residual stress, water pressure, and other states from the end of the previous step to the new mesh.

[0081] Adaptive time step setting: If the energy release rate approaches the rock critical value, the calculation time step is shortened to ensure that the crack propagation is captured sufficiently accurately.

[0082] The output generates a "restart file" for each step, recording deleted grids, activated grids, and remaining field variables.

[0083] In this embodiment, the cubic grid containing the mined coal body is marked as a goaf state in the calculation model and no longer participates in bearing rigidity, but is still retained in the new calculation area. The weak material parameters of the collapse area are assigned to accurately simulate the stress disturbance of the goaf on the overlying rock mass, the induction effect of crack propagation path and the change of gas seepage channels. This ensures that the goaf, working face and fine mesh area are dynamically updated and reasonably considered in the multi-field coupled solution that progresses over time, and are represented by the weak material parameters of the collapse area.

Claims

1. A three-dimensional simulation method for coal mining based on numerical manifolds and structured grids, characterized in that, include: Determine the initial location of the coal mine working face, collect on-site data of the coal mine working face, unify the coordinate system of the on-site data of the coal mine working face into a single HDF5 file, and obtain the on-site database; The fine mesh range is defined, and the goaf, the coal mine working face, and the fine mesh range are used as the calculation area; the fine mesh range is set in the direction of coal mine working face advancement and is adjacent to the coal mine working face. The computational domain is divided into grids and integrated with the field database to obtain a grid model that incorporates field data. A numerical manifold cover is constructed and coupled with a mesh model fused with field data. Specifically, the construction of the numerical manifold cover involves extracting the fracture centerlines of each fracture surface in the mesh model, and arranging several layers of cover elements on both sides of the fracture centerlines to describe fracture opening and closing and block movement. The coverage influence radius of each cover element is: in, The coverage influence radius of the coverage unit; This is an empirical coefficient; The elastic modulus of the rock; This is the critical value for the fracture energy release rate; Tensile strength of rock; The maximum stress difference between the elements before and after subdivision; The reference maximum stress; Coverage data is obtained based on the coupling results of numerical manifold cover and the mesh model fused with field data; Based on the coverage data, the three-dimensional stress field, displacement field, and pore gas pressure field of the current computational region are solved using multi-field coupling based on a CPU-GPU heterogeneous parallel structure; specifically, the solution of the three-dimensional stress field, displacement field, and pore gas pressure field of the current computational region over time using multi-field coupling based on a CPU-GPU heterogeneous parallel structure is as follows: On the GPU, the central difference method is used to explicitly integrate the structured mesh nodes and update the displacement field: in, For the displacement increment of the cube mesh; Density; For time step; For gradient operators; For stress tensor; This is the stress acceleration vector; Minimum unit size; It is the elastic wave velocity; The elastic modulus of the rock; Here is the tangent stiffness matrix; The displacement of the covering element; It is a nonlinear residue; On the GPU, the stress field is calculated using an elastic constitutive model based on the strain tensor. At the same time, the local stress within the cover element is updated by combining the normal opening and closing model and the shear slip model of the crack region. The gas pressure field is implicitly solved on the CPU based on a dual-medium seepage model, which incorporates permeability controlled by fracture aperture and gas viscosity corrected by the Klinkenberg slip effect. in, The volumetric flow velocity density in the fracture; This represents the current crack permeability. This represents the pressure gradient within the fracture. For gradient operators; This refers to the fracture pressure; The corrected gas viscosity; Reference viscosity; The Klinkenberg coefficient; The three-dimensional stress field, displacement field, and pore gas pressure field are iteratively solved on both GPU and CPU based on the coupling of fracture aperture, stress state, and seepage parameters. After each iteration at time step, the overall residual and energy conservation error are calculated. If the overall residual is lower than a set threshold and the energy conservation error is less than the error threshold, a three-dimensional stress field, displacement field, and pore gas pressure field that evolve over time are formed; otherwise, an adaptive evolution process is initiated. The adaptive evolution process is specifically as follows: Calculate the current energy release rate of the fracture: in, Energy release rate; For strain energy; For external skills; The area of ​​the crack; When the current crack energy release rate is greater than the release rate threshold At that time, crack propagation is triggered; the crack extension amount of the crack propagation is: in, Extending for a single-step forward; For control coefficients; The side length of the fine mesh unit; Based on the amount of crack extension, a covering unit is automatically generated at the new crack location, and the mapping relationship is updated simultaneously. Update crack permeability: in, The updated crack permeability; The permeability of the cracks before the update; The penetration-aperture coupling coefficient; The crack opening; If the current energy release rate of the fracture is close to the critical value of the rock, shorten the calculation time step; Based on the updated fracture permeability, shortened computation time step, and updated overlay cells, the multi-field coupling solution is re-executed in the updated mesh and fracture field; The cubic grid containing the mined coal body is marked as the mined-out state, and the weak material parameters of the collapse zone are used to represent it. The next step is to advance the coal mine working face until the entire working face is advanced, and then obtain the stress distribution changes, fluid migration changes, and crack propagation changes throughout the entire working face advancement process.

2. The three-dimensional simulation method for coal mining based on numerical manifolds and structured grids according to claim 1, characterized in that, The field database includes a list of mechanical parameters for each rock stratum, a list of seepage parameters for each rock stratum, a three-dimensional fracture surface file of the working face, a pore water pressure field, a pore permeability coefficient field, and a gas seepage parameter table.

3. The three-dimensional simulation method for coal mining based on numerical manifolds and structured grids according to claim 1, characterized in that, The obtained mesh model based on the field data is specifically as follows: the computational area is divided into several cubic grids according to the rock layer thickness, and the cubic grids within the fine mesh range are further subdivided, and cubic grids with deformation greater than the deformation threshold are removed to obtain a multi-scale cubic mesh file; The data in the field database is mapped to the grid nodes in the multi-scale cubic grid file through a unified coordinate transformation matrix to obtain a grid model that integrates field data.

4. The three-dimensional simulation method for coal mining based on numerical manifolds and structured grids according to claim 1, characterized in that, The covering unit is specifically a number of virtual nodes arranged along the normal direction on both sides of the crack centerline, and a crack constitutive model is set between each virtual node. The covering unit is sparsely mapped and coupled with the grid nodes of the grid model that integrates field data through a numerical manifold method to simulate the crack opening and closing process and block movement. The displacement degrees of freedom of the covering unit are independent of the mesh model that integrates field data.

5. The three-dimensional simulation method for coal mining based on numerical manifolds and structured grids according to claim 4, characterized in that, The coupled numerical manifold cover and the fused field data mesh model specifically involve mapping the displacement information of the cover cells onto the mesh nodes of the fused field data mesh model: in, It is a sparse mapping matrix; The displacement vector of the grid nodes; This is the displacement vector of the covering element node.

6. The three-dimensional simulation method for coal mining based on numerical manifolds and structured grids according to claim 1, characterized in that, The overlay data includes a normal opening / closing model to describe the stiffness change during crack opening and closing, a shear slip model to describe the relationship between shear force and slip displacement on the crack surface, and crack permeability: in, The normal stress in the crack; The initial normal stiffness; The crack opening; This is the initial closing opening. This refers to the residual normal stress; Shear stress; It is the friction angle; It is cohesive force; This represents the current crack permeability. Initial penetration rate; The penetration-aperture coupling coefficient; This refers to the crack opening.

7. The three-dimensional simulation method for coal mining based on numerical manifolds and structured grids according to claim 1, characterized in that, The overall residual is: in, For the overall residual; Here is the tangent stiffness matrix; The displacement of the covering element; This is the external load vector; The energy conservation error for: in, The internal energy of the system; The kinetic energy of the system; Work is done by the external forces acting on the system.

Images