Coal mining three-dimensional simulation method based on numerical manifold and structured grid

Abstract

The invention discloses a coal mining three-dimensional simulation method based on a numerical manifold and a structured grid, and belongs to the field of coal mining three-dimensional simulation. The method comprises the steps that the initial position of a working face is determined, and a field database is collected; performing grid division on the calculation area, and fusing the field database to obtain a grid model fused with the field data; constructing and coupling a grid model of numerical manifold coverage and fusion field data; coupling of fracture mechanical behaviors and seepage is achieved; based on the CPU-GPU heterogeneous parallel structure, solving a three-dimensional stress field, a displacement field and a pore gas pressure field evolved along with time of the current calculation area; and advancing the coal mine working face in the next step until the whole advancing of the working face is completed, and obtaining the stress distribution change, the fluid migration change and the fracture extension change in the whole advancing process of the working face. The method solves the problems that the existing method is difficult to describe the stress accumulation of the initial microfracture with high precision and is difficult to maintain the topological update of contact / separation after fracture.

Description

Technical Field

[0001] The present invention belongs to the field of three-dimensional simulation of coal mining, and in particular relates to a three-dimensional simulation method of coal mining based on numerical manifolds and structured grids. Background Art

[0002] During the advancement of a coal mine longwall working face, the surrounding rock undergoes a continuous-discontinuous multi-field coupled evolution, which can easily induce dynamic hazards such as roof instability and rock bursts. To reveal the three-dimensional stress wave propagation, fracture propagation, and dynamic load dissipation during the mining process, the engineering community generally uses finite element method (FEM), discrete element method (DEM), or explicit / implicit coupling methods for numerical simulation. However, in highly stressed, multi-fractured coal and rock masses, existing methods generally suffer from the following bottlenecks: (1) Traditional structured grids are difficult to depict complex fracture networks. Multi-scale bedding, faults, and induced fractures are often distributed around coal mine working faces. Traditional FEM relies on single-domain grids with continuous node sharing. If fractures are directly embedded in a three-dimensional structured grid, the grid needs to be locally subdivided or reconstructed, which will dramatically increase the unit distortion and computational scale. This results in "difficulty in three-dimensional structural modeling under the condition of multiple fracture surfaces and difficulty in controlling unit geometric resolution; (2) The efficiency and stability of the dual-medium fluid-solid coupling solution are low. The unloading and reloading caused by the advancement of the working face are accompanied by seepage and thermo-chemical effects. Pure FEM or DEM requires an extremely small time step to ensure convergence when dealing with strong nonlinear problems such as fracture penetration and block slip, resulting in "complex numerical solution of fluid-solid coupling and large amount of calculation; (3) Insufficient capture of continuous-discontinuous transitions. During the top cutting process, the rock mass gradually transitions from a continuous body to a block motion phase. Although explicit DEM can handle discrete blocks, it is difficult to accurately describe the stress accumulation in the initial microcracks. Traditional FEM also struggles to maintain contact / separation topological updates after fracture. Summary of the Invention

[0003] In response to the above-mentioned deficiencies in the prior art, the present invention provides a three-dimensional simulation method for coal mining based on numerical manifolds and structured grids, which solves the problems that the existing methods are difficult to describe the stress accumulation of initial microcracks with high precision and difficult to maintain the topological update of contact / separation after rupture.

[0004] In order to achieve the above-mentioned object, the present invention adopts a technical solution: a three-dimensional simulation method for coal mining based on numerical manifolds and structured grids, comprising: Determine the initial position of the coal mine working face, collect the field data of the coal mine working face, and integrate the field data of the coal mine working face into a unified coordinate system and the same HDF5 file to obtain the field database; Set a fine mesh range, and use the goaf, coal mine working face, and fine mesh range as calculation areas; the fine mesh range is set in the coal mine working face advancing direction and is adjacent to the coal mine working face; Grid the calculation area and integrate it with the field database to obtain a grid model that integrates field data; Construct a numerical manifold overlay and couple it with a grid model that integrates field data; The coverage data is obtained based on the coupling results of the numerical manifold coverage and the grid model fused with the field data; Based on the coverage data, the CPU-GPU heterogeneous parallel structure multi-field coupling is used to solve the time-evolving three-dimensional stress field, displacement field, and pore gas pressure field in the current calculation area; The cube grid where the mined coal body is located is marked as the mined-out state, and the weak material parameters of the collapse area are used instead; The next step of advancing the coal mine working face is to be carried out until the entire advancement of the working face is completed, and the stress distribution changes, fluid migration changes and crack expansion changes of the entire process of the working face advancement are obtained.

[0005] Furthermore, the field database includes a list of mechanical parameters of each rock layer, a list of seepage parameters of each rock layer, 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 grid model obtained by integrating the field data is specifically as follows: the calculation area is divided into a number of cubic grids according to the thickness of the rock layer, and the cubic grids in the fine mesh range are subdivided, and the cubic grids with deformation greater than the deformation threshold are eliminated to obtain a multi-scale cubic grid file; the data in the field database are mapped to the grid nodes in the multi-scale cubic grid file through a unified coordinate transformation matrix to obtain the grid model integrating the field data.

[0007] Furthermore, the construction of the numerical manifold cover is specifically as follows: extracting the fracture centerline of each fracture surface of the grid model, arranging several layers of cover units on both sides of the fracture centerline for describing the opening and closing of the fracture and the movement of the block; the coverage influence radius of the cover unit is:

[0008]

[0009] in, is the coverage influence radius of the coverage unit; is the empirical coefficient; is the elastic modulus of rock; is the critical value of fracture energy release rate; is the tensile strength of rock; is the maximum stress difference of the unit before and after subdivision; is the baseline maximum stress.

[0010] Furthermore, the covering unit is specifically a number of virtual nodes arranged along the normal direction on both sides of the crack center line, 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 integrated with field data through the numerical manifold method to simulate the crack opening and closing process and block movement; The displacement degrees of freedom of the covering element are independent of the grid model integrated with the field data.

[0011] Furthermore, the numerical manifold overlay is coupled with the grid model of the fused field data, specifically by mapping the displacement information of the overlay unit to the grid nodes of the grid model of the fused field data:

[0012] in, is a sparse mapping matrix; is the mesh node displacement vector; is the displacement vector of the covering element node.

[0013] Furthermore, the coverage data includes a normal opening and closing model for describing the stiffness change when the crack opens and closes, a shear slip model for describing the relationship between shear force and slip displacement on the crack surface, and crack permeability:

[0014]

[0015]

[0016] in, is the crack normal stress; is the initial normal stiffness; is the crack opening; is the initial closed opening; is the residual normal stress; is the shear stress; is the friction angle; for cohesion; is the current fracture permeability; is the initial permeability; is the permeability-openness coupling coefficient; is the crack opening.

[0017] Furthermore, the CPU-GPU heterogeneous parallel structure multi-field coupling is used to solve the time-evolving three-dimensional stress field, displacement field, and pore gas pressure field in the current calculation area, specifically: The displacement field is updated by explicitly integrating the structured grid nodes using the central difference method on the GPU:

[0018]

[0019]

[0020]

[0021] in, is the displacement increment of the cube grid; is the density; is the time step; is the gradient operator; is the stress tensor; is the stress acceleration vector; is the minimum unit size; is the elastic wave velocity; is the elastic modulus of rock; is the tangent stiffness matrix; is the displacement of the covering unit; is the nonlinear residual; The stress field is calculated on the GPU based on the strain tensor through the elastic constitutive model, and the local stress in the cover element is updated by combining the normal opening and closing model and the shear slip model in the crack area. The gas pressure field is implicitly solved on the CPU based on the dual-medium seepage model, which introduces the permeability controlled by the fracture aperture and the gas viscosity corrected by the Klinkenberg slippage effect:

[0022]

[0023] in, is the volume flow density in the fracture; is the current fracture permeability; is the pressure gradient in the fracture; is the gradient operator; is the fracture pressure; is the corrected gas viscosity; is the base viscosity; is the Klinkenberg coefficient; On the GPU and CPU side, the three-dimensional stress field, displacement field, and pore gas pressure field are iteratively solved based on the coupling of fracture aperture, stress state, and seepage parameters: After the iteration of each time step is completed, 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 with time are formed. Otherwise, the adaptive evolution process begins.

[0024] Furthermore, the overall residual is:

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

[0026] in, is the internal energy of the system; is the kinetic energy of the system; The work done by the external force on the system.

[0027] Furthermore, the adaptive evolution process is specifically as follows: Calculate the current crack energy release rate:

[0028] in, is the energy release rate; is the strain energy; For external gong; is the crack area; When the current crack energy release rate is greater than the release rate threshold When , the crack expansion is triggered; the crack extension amount of the crack expansion is:

[0029] in, Extended for a single-step forward; is the control coefficient; is the side length of the fine mesh unit; Based on the crack extension amount of the crack expansion, the cover unit is automatically generated at the new crack position and the mapping relationship is updated at the same time; Update fracture permeability:

[0030] in, is the updated fracture permeability; is the fracture permeability before updating; is the permeability-openness coupling coefficient; is the crack opening; Then re-execute the multi-field coupling calculation in the updated mesh and fracture field.

[0031] The beneficial effects of the present invention are as follows: the present invention seamlessly couples the regular cubic grid with the numerical manifold covering unit. The continuous medium still uses an efficient structured grid, and the cracks and block motions only appear in the form of covering locally, and the two maintain the consistency of displacement and force through the mapping matrix. This not only avoids the unit distortion caused by forcibly cutting into the cracks in the grid, but also reserves space for the free growth of cracks. Actual examples show that within the same cross-sectional range, the total amount of grid can be reduced by more than 40% compared with pure finite elements, and the calculation time is compressed to one-third of the original. At the same time, the errors of the crack initiation and extension positions and the on-site microseismic positioning are controlled within more than ten meters. In response to the dynamic process of "advancement-unloading" of the working face, this method proposes a modeling strategy of "local fine mesh at the front edge + real-time deletion / activation of the goaf". The front area is automatically subdivided to ensure the spatial resolution of the unloading wave process; and the mined area is immediately deleted and replaced with collapsed materials, avoiding the memory explosion caused by the global laying of fine mesh. Large-scale testing has shown that, at the same accuracy, graphics memory usage is only 30% to 40% of that of traditional global fine-mesh solutions. The wall-clock time for a single advance step remains stable at within one to two hours, meeting shift-level decision-making requirements. This invention utilizes a heterogeneous parallel framework of "GPU explicit integration + CPU implicit iteration." Linear elastic-plastic calculations on the structured mesh are entirely handled by the graphics card's parallel cores, while the strongly nonlinear equations within the coverage domain are solved using Newton-Krylov iterations on the multi-core CPU. The two types of results are exchanged within sub-iterations until convergence. Compared to pure CPU solvers, this division of labor reduces solution time by approximately three-quarters and ensures numerical stability during the intense crack opening phase. Crack propagation is criterioned by the energy release rate, with incremental coverage generated instantly at the propagation location. Whether a crack continues to grow no longer depends on local grid cell instability, but rather on whether the energy exerted on the rock exceeds a critical value. Therefore, the crack extension direction is not constrained by the grid axis, avoiding the "perpendicular crack" problem of traditional finite element methods. Compared with the indoor true triaxial loading test, the simulated fracture inclination error is less than five degrees, and the roof plate separation path is smoother and more reliable. In terms of fluid coupling, the model simultaneously solves two sets of pressure fields for the matrix pores and fracture channels, and performs Klinkenberg correction on the viscosity of low-pressure gas flow. This improvement significantly improves the prediction accuracy of gas-rich areas: in typical mine tests, the spatial positioning accuracy of high-gas concentration areas has increased from 70% to more than 90%, which can be directly used to guide the layout of extraction drilling holes. From the very beginning, the present invention uses four major types of field measured data from drilling, geophysical exploration, hydrology, and gas to drive grid and coverage generation, reducing a large amount of secondary adjustment work that relies on empirical parameters. The deviation of the model's initial stress and pore pressure from the actual measurement is controlled within ten percent, laying a reliable foundation for subsequent dynamic evolution calculations. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] Figure 1 Flow chart of the method of the present invention. DETAILED DESCRIPTION

[0033] The specific embodiments of the present invention are described below to facilitate understanding of the present invention by those skilled in the art. However, it should be clear that the present invention is not limited to the scope of the specific embodiments. For those skilled in the art, as long as various changes are within the spirit and scope of the present invention as defined and determined by the appended claims, these changes are obvious, and all inventions and creations utilizing the concepts 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 position of the coal mine working face, collect the field data of the coal mine working face, and integrate the field data of the coal mine working face into a unified coordinate system and the same HDF5 file to obtain the field database; Set a fine mesh range, and use the goaf, coal mine working face, and fine mesh range as calculation areas; the fine mesh range is set in the coal mine working face advancing direction and is adjacent to the coal mine working face; Grid the calculation area and integrate it with the field database to obtain a grid model that integrates field data; Construct a numerical manifold overlay and couple it with a grid model that integrates field data; The coverage data is obtained based on the coupling results of the numerical manifold coverage and the grid model fused with the field data; Based on the coverage data, the CPU-GPU heterogeneous parallel structure multi-field coupling is used to solve the time-evolving three-dimensional stress field, displacement field, and pore gas pressure field in the current calculation area; The cube grid where the mined coal body is located is marked as the mined-out state, and the weak material parameters of the collapse area are used instead; The next step of advancing the coal mine working face is to be carried out until the entire advancement of the working face is completed, and the stress distribution changes, fluid migration changes and crack expansion changes of the entire process of the working face advancement are obtained.

[0035] The field database includes a list of mechanical parameters of each rock layer, a list of seepage parameters of each rock layer, 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 collection is performed. Experimental data such as stratum thickness, core photographs, and rock strength and permeability are obtained through advance drilling or tunnel drilling. A list of initial mechanical and seepage parameters for each rock layer is established. Geophysical data collection. Transient electromagnetic, seismic wave, geological radar and other geophysical data collection. Generate three-dimensional fracture / fault surface files of the working face; Collect hydrogeological and gas geological data. Conduct pumping or injection tests to determine the aquifer water pressure and permeability, and obtain the original pore water pressure and permeability fields. Determine gas content, gas pressure, and adsorption constants, and establish a gas seepage parameter table. Data consolidation. Unify the coordinate system and integrate multiple data sources into a single HDF5 file. Generate a standardized field database for subsequent steps.

[0037] The grid model for fusing field data is obtained by specifically dividing the calculation area into a plurality of cubic grids according to the thickness of the rock layer, subdividing the cubic grids in the fine mesh range, and eliminating the cubic grids with deformation greater than the deformation threshold to obtain a multi-scale cubic grid file; and mapping the data in the field database to the grid nodes in the multi-scale cubic grid file through a unified coordinate transformation matrix to obtain the grid model for fusing field data.

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

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

[0040] Front subdivision: Within thirty to forty meters in front of the working face, the cube is subdivided into small grids that are one-quarter to one-eighth the size of the main grid to make the unloading wave simulation more detailed.

[0041] Quality check: Remove highly deformed grids to ensure that each grid has a good shape without severe stretching or distortion, and obtain qualified multi-scale cube mesh files.

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

[0043] The construction of the numerical manifold cover is specifically as follows: extracting the fracture centerline of each fracture surface of the grid model, arranging several layers of cover units on both sides of the fracture centerline for describing the opening and closing of the fracture and the movement of the block; the coverage influence radius of the cover unit is:

[0044]

[0045] in, is the coverage influence radius of the coverage unit; is the empirical coefficient; is the elastic modulus of rock; is the critical value of fracture energy release rate; is the tensile strength of rock; is the maximum stress difference of the unit before and after subdivision; is the baseline maximum stress.

[0046] In this implementation, four to eight layers of "covering elements" are placed on either side of each fracture. These layers act like meshless patches, freely describing the crack opening and closing and block motion on the structural grid. The cover's influence radius is adaptively set based on the rock's toughness and tensile strength. The cover elements are linked to the nodes of the structural grid to ensure that displacements and forces are transmitted between the two discrete systems.

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

[0048] The covering units are 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 integrated with field data through the numerical manifold method to simulate the crack opening and closing process and block movement; The displacement degrees of freedom of the covering element are independent of the grid model integrated with the field data.

[0049] In this embodiment, the covering unit describes the normal opening and closing of the crack and shear slip by arranging multiple layers of 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 degrees of freedom are independent of the background grid, and can capture the relative displacement and opening changes between the blocks on both sides of the crack. It is sparsely mapped and coupled with the main grid nodes through the numerical manifold method, thereby achieving high-precision simulation of the crack opening and closing process and block motion without modifying the original structural grid.

[0050] The coupling of numerical manifold overlay and grid model of fused field data is to map the displacement information of overlay cells to the grid nodes of the grid model of fused field data:

[0051] in, is a sparse mapping matrix; is the mesh node displacement vector; is the displacement vector of the covering element node.

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

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

[0054]

[0055]

[0056] in, is the crack normal stress; is the initial normal stiffness; is the crack opening; is the initial closed opening; is the residual normal stress; is the shear stress; is the friction angle; for cohesion; is the current fracture permeability; is the initial permeability; is the permeability-openness coupling coefficient; is the crack opening.

[0057] In this embodiment, the nonlinear mechanical response of the crack during compression closure or tension opening is simulated by setting the normal stiffness of the crack to change with the crack opening or closing state. The relationship between the shear force and the slip displacement on the crack surface is described by parameters such as the friction coefficient, bonding strength, and residual strength, reflecting the bonding-slip mechanism before and after the crack surface is destroyed. The initial permeability of the crack is set to increase exponentially with the opening. The permeability-stress coupling formula expresses the nonlinear response relationship of permeability exponentially decaying with the increase of normal stress and rapidly increasing with the opening of the crack through the change of the crack opening controlled by stress. Fitted by borehole transient test.

[0058] The CPU-GPU heterogeneous parallel structure multi-field coupling solution is used to solve the time-evolving three-dimensional stress field, displacement field and pore gas pressure field in the current calculation area. Specifically, The displacement field is updated by explicitly integrating the structured grid nodes using the central difference method on the GPU:

[0059]

[0060]

[0061]

[0062] in, is the displacement increment of the cube grid; is the density; is the time step; is the gradient operator; is the stress tensor; is the stress acceleration vector; is the minimum unit size; is the elastic wave velocity; is the elastic modulus of rock; is the tangent stiffness matrix; is the displacement of the covering unit; is the nonlinear residual; The stress field is calculated on the GPU based on the strain tensor through the elastic constitutive model, and the local stress in the cover element is updated by combining the normal opening and closing model and the shear slip model in the crack area. The gas pressure field is implicitly solved on the CPU based on the dual-medium seepage model, which introduces the permeability controlled by the fracture aperture and the gas viscosity corrected by the Klinkenberg slippage effect:

[0063]

[0064] in, is the volume flow density in the fracture; is the current fracture permeability; is the pressure gradient in the fracture; is the gradient operator; is the fracture pressure; is the corrected gas viscosity; is the base viscosity; is the Klinkenberg coefficient; On the GPU and CPU side, the three-dimensional stress field, displacement field, and pore gas pressure field are iteratively solved based on the coupling of fracture aperture, stress state, and seepage parameters: After the iteration of each time step is completed, 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 with time are formed. Otherwise, the adaptive evolution process begins.

[0065] In this implementation, an explicit central difference method updates the displacements and stresses of a cube mesh on the graphics card. By explicitly integrating the mass acceleration obtained by dividing the nodal force by the mass, the time step is explicitly integrated to update the displacements and velocities point by point, achieving efficient and parallel explicit dynamics calculations. The time step must satisfy the CFL condition to ensure that the wave velocity propagation does not exceed the mesh size within a single time step, thereby ensuring numerical stability and the correct capture of physical causal relationships. The implicit iterative formula simultaneously solves for unknown quantities within the current time step, making the numerical solution of the system stable and capable of handling large time steps.

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

[0067] In this embodiment, a multi-field coupling solution is performed based on CPU-GPU heterogeneous parallelism: 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 evolving with time are obtained. In this method, the time evolution process of the three-dimensional stress field, displacement field, and pore gas pressure field is realized through multi-field coupling iteration: the central difference method is used on the GPU to explicitly integrate the structured grid nodes and update the displacement field; the stress field is calculated according to the strain tensor through the elastic or elastoplastic constitutive model, and the local stress in the covering unit is updated in combination with the normal opening and closing model and the shear slip model of the fracture area; the gas pressure field is implicitly solved on the CPU based on the dual-medium seepage model, and the permeability controlled by the fracture aperture and the viscosity corrected by Klinkenberg are introduced; the three are iteratively coupled with each other through the fracture 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 with time. Otherwise, the adaptive evolution process begins: first, the current fracture energy release rate is calculated. If it exceeds a critical threshold, fracture extension is triggered, new cover cells are generated along the fracture front, and the fracture permeability is updated. At the same time, the numerical manifold cover structure and node mapping relationship are adjusted according to the distribution of the latest cover cells. Then, the multi-field coupling calculation is re-executed in the updated mesh and fracture field, and iterations are continued 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 fracture extension and fluid migration.

[0068] The overall residual is:

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

[0070] in, is the internal energy of the system; is the kinetic energy of the system; The work done by the external force on the system.

[0071] In this embodiment, the overall residual formula reflects the errors of the system mechanics and fluid conservation equations by measuring the inconsistency between the current iterative solution and the equilibrium equation, thereby determining the degree of convergence of the numerical solution.

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

[0073] in, is the energy release rate; is the strain energy; For external gong; is the crack area; When the current crack energy release rate is greater than the release rate threshold When , the crack expansion is triggered; the crack extension amount of the crack expansion is:

[0074] in, Extended for a single-step forward; is the control coefficient; is the side length of the fine mesh unit; Based on the crack extension amount of the crack expansion, the cover unit is automatically generated at the new crack position and the mapping relationship is updated at the same time; Update fracture permeability:

[0075] in, is the updated fracture permeability; is the fracture permeability before updating; is the permeability-openness coupling coefficient; is the crack opening; Then re-execute the multi-field coupling calculation in the updated mesh and fracture 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 extend. The crack extension formula proportionally determines the distance the crack front advances based on how much the energy release rate exceeds the critical value, enabling adaptive adjustment of the crack extension step size.

[0077] The fracture front extends proportionally forward by one or more grid lengths at a time, with the specific extension determined by the degree to which the energy release rate exceeds a critical value. A patch cover is automatically generated at the new fracture location, and the mapping relationship is updated simultaneously. New cover nodes inherit the frontier, and the mapping matrix is ​​adaptively updated to ensure that the global bandwidth increase is less than 3%. Once a fracture opens, the permeability increases exponentially; it decreases as the fracture closes. The permeability update formula adjusts the permeability based on the exponential change in the fracture aperture, achieving a dynamic response with a rapid increase in permeability when the fracture opens and a rapid decrease when it closes.

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

[0079] Goaf deletion: After completing each step, the grid where the mined coal body is located is deleted and replaced by the weak material parameters of the collapse area.

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

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

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

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

Claims

1. A three-dimensional simulation method for coal mining based on numerical manifolds and structured grids, characterized in that: include: Determine the initial position of the coal mine working face, collect the field data of the coal mine working face, and integrate the field data of the coal mine working face into a unified coordinate system and the same HDF5 file to obtain the field database; Set a fine mesh range, and use the goaf, coal mine working face, and fine mesh range as calculation areas; the fine mesh range is set in the coal mine working face advancing direction and is adjacent to the coal mine working face; Grid the calculation area and integrate it with the field database to obtain a grid model that integrates field data; Construct a numerical manifold overlay and couple it with a grid model that integrates field data; Obtain coverage data based on the coupling results of numerical manifold coverage and grid models fused with field data; Based on the coverage data, the CPU-GPU heterogeneous parallel structure multi-field coupling is used to solve the time-evolving three-dimensional stress field, displacement field, and pore gas pressure field in the current calculation area; The cube grid where the mined coal body is located is marked as the mined-out state, and the weak material parameters of the collapse area are used instead; The next step of advancing the coal mine working face is to be carried out until the entire advancement of the working face is completed, and the stress distribution changes, fluid migration changes and crack expansion changes of the entire process of the working face advancement are obtained.

2. The three-dimensional simulation method for coal mining based on numerical manifold and structured grid according to claim 1, characterized in that: The field database includes a list of mechanical parameters of each rock layer, a list of seepage parameters of each rock layer, 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 manifold and structured grid according to claim 1, characterized in that: The grid model obtained by integrating the field data is specifically as follows: the calculation area is divided into a plurality of cubic grids according to the thickness of the rock layer, and the cubic grids in the fine mesh range are subdivided, and the cubic grids with deformation greater than the deformation threshold are eliminated to obtain a multi-scale cubic grid file; The data in the field database are mapped to the grid nodes in the multi-scale cube grid file through a unified coordinate transformation matrix to obtain a grid model that integrates the field data.

4. The three-dimensional simulation method for coal mining based on numerical manifold and structured grid according to claim 1, characterized in that: The construction of the numerical manifold cover is specifically as follows: extracting the fracture centerline of each fracture surface of the grid model, arranging several layers of cover units on both sides of the fracture centerline for describing the opening and closing of the fracture and the movement of the block; the coverage influence radius of the cover unit is: in, is the coverage influence radius of the coverage unit; is the empirical coefficient; is the elastic modulus of rock; is the critical value of fracture energy release rate; is the tensile strength of rock; is the maximum stress difference of the unit before and after subdivision; is the baseline maximum stress.

5. The three-dimensional simulation method for coal mining based on numerical manifold and structured grid according to claim 4, characterized in that: The covering units are 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 integrated with field data through the numerical manifold method to simulate the crack opening and closing process and block movement; The displacement degrees of freedom of the covering element are independent of the grid model integrated with the field data.

6. The three-dimensional simulation method for coal mining based on numerical manifold and structured grid according to claim 5, characterized in that: The coupling of numerical manifold overlay and grid model of fused field data is to map the displacement information of overlay cells to the grid nodes of the grid model of fused field data: in, is a sparse mapping matrix; is the mesh node displacement vector; is the displacement vector of the covering element node.

7. The three-dimensional simulation method for coal mining based on numerical manifold and structured grid according to claim 1, characterized in that: The coverage data includes a normal opening and closing model for describing the stiffness change during crack opening and closing, a shear slip model for describing the relationship between shear force and slip displacement on the crack surface, and crack permeability: in, is the crack normal stress; is the initial normal stiffness; is the crack opening; is the initial closed opening; is the residual normal stress; is the shear stress; is the friction angle; for cohesion; is the current fracture permeability; is the initial permeability; is the permeability-openness coupling coefficient; is the crack opening.

8. The three-dimensional simulation method for coal mining based on numerical manifold and structured grid according to claim 1, characterized in that: The CPU-GPU heterogeneous parallel structure multi-field coupling solution is used to solve the time-evolving three-dimensional stress field, displacement field and pore gas pressure field in the current calculation area. Specifically, The displacement field is updated by explicitly integrating the structured grid nodes using the central difference method on the GPU: in, is the displacement increment of the cube grid; is the density; is the time step; is the gradient operator; is the stress tensor; is the stress acceleration vector; is the minimum unit size; is the elastic wave velocity; is the elastic modulus of rock; is the tangent stiffness matrix; is the displacement of the covering unit; is the nonlinear residual; The stress field is calculated on the GPU based on the strain tensor through the elastic constitutive model, and the local stress in the cover element is updated by combining the normal opening and closing model and the shear slip model in the crack area. The gas pressure field is implicitly solved on the CPU based on the dual-medium seepage model, which introduces the permeability controlled by the fracture aperture and the gas viscosity corrected by the Klinkenberg slippage effect: in, is the volume flow density in the fracture; is the current fracture permeability; is the pressure gradient in the fracture; is the gradient operator; is the fracture pressure; is the corrected gas viscosity; is the base viscosity; is the Klinkenberg coefficient; On the GPU and CPU side, the three-dimensional stress field, displacement field, and pore gas pressure field are iteratively solved based on the coupling of fracture aperture, stress state, and seepage parameters: After the iteration of each time step is completed, 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 with time are formed. Otherwise, the adaptive evolution process begins.

9. The three-dimensional simulation method for coal mining based on numerical manifold and structured grid according to claim 8, characterized in that: The overall residual is: in, is the overall residual; is the tangent stiffness matrix; is the displacement of the covering unit; is the external load vector; The energy conservation error for: in, is the internal energy of the system; is the kinetic energy of the system; The work done by the external force on the system.

10. The three-dimensional simulation method for coal mining based on numerical manifold and structured grid according to claim 1, characterized in that: The adaptive evolution process is specifically as follows: Calculate the current crack energy release rate: in, is the energy release rate; is the strain energy; For external gong; is the crack area; When the current crack energy release rate is greater than the release rate threshold When , the crack expansion is triggered; the crack extension amount of the crack expansion is: in, Extended for a single-step forward; is the control coefficient; is the side length of the fine mesh unit; Based on the crack extension amount of the crack expansion, the cover unit is automatically generated at the new crack position and the mapping relationship is updated at the same time; Update fracture permeability: in, is the updated fracture permeability; is the fracture permeability before updating; is the permeability-openness coupling coefficient; is the crack opening; If the current crack energy release rate approaches the rock critical value, shorten the calculation time step; Based on the updated fracture permeability, shortened computational time step, and updated cover cells, the multi-field coupling solution is re-executed in the updated mesh and fracture field.

Images