Underground space excavation simulation method, device, equipment and medium

Abstract

According to the underground space excavation simulation method, device and equipment and the medium, a smooth fluid dynamics method is adopted to disperse a soil body of a to-be-simulated area into a particle system, and the spatial position and physical parameters of boundary virtual particles are determined; according to the spatial position and the physical parameter of the boundary virtual particle, determining the interaction between the particles, and determining the spin rate tensor and the relative speed of the particles; according to the spin rate tensor and the relative speed of the particles, the stress rate of the soil body particles is calculated in combination with the constitutive relation of the particles; performing time integration on the stress rate, and updating the stress state of each particle; judging whether the updated stress state exceeds the damage surface or not so as to determine whether the particles are damaged or not; marking the damaged particles, and outputting a simulation result. According to the scheme, under complex geological conditions and discontinuous deformation scenes, the simulation precision and stability are high, and the safety of the construction process can be guaranteed.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of rock mechanics numerical simulation, and in particular to a method, device, equipment and medium for simulating underground space excavation. BACKGROUND

[0002] Pumped storage power stations have become an important part of modern power systems due to their ability to respond quickly, adjust peak load flexibly, and operate stably for a long period, providing key support for efficient utilization of renewable energy. Abandoned coal mines are considered ideal sites for building pumped storage underground reservoirs due to their vast underground space resources and natural elevation conditions. These mines not only significantly reduce reservoir construction costs, but also have the dual benefits of environmental remediation and resource reuse. In this context, the "mine-energy synergy" model based on mine resource utilization has emerged, providing a new direction for the development of abandoned mines and opening up new avenues for the storage and deployment of renewable energy. However, the transformation of abandoned coal mine underground reservoirs involves complex engineering technical problems, including the damage behavior and mechanical properties of surrounding rock, which have a significant impact on the safety and stability of the reservoir.

[0003] In the prior art, numerical simulation methods have been widely used in the analysis of the mechanical response of underground space excavation. However, existing numerical simulation methods have poor precision and stability under complex geological conditions and non-continuous deformation scenarios, making it difficult to ensure the safety of the construction process. SUMMARY

[0004] To address the above-mentioned defects, the present application provides a method, device, equipment and medium for simulating underground space excavation, which can achieve high precision and stability in simulating complex geological conditions and non-continuous deformation scenarios, ensuring the safety of the construction process.

[0005] The present application provides a method for simulating underground space excavation, which comprises:

[0006] The soil in the region to be simulated is discretized into a particle system using a smooth fluid dynamics method, and the spatial position and physical parameters of the boundary virtual particles are determined;

[0007] The interaction between particles is determined based on the spatial position and physical parameters of the boundary virtual particles, and the spin rate tensor and relative velocity of the particles are determined;

[0008] The stress rate of the soil particles is calculated based on the spin rate tensor and relative velocity of the particles, combined with the particle constitutive relationship;

[0009] The stress state of each particle is updated by time integration of the stress rate;

[0010] It is determined whether the updated stress state exceeds the failure surface, thereby determining whether the particle is damaged;

[0011] marking the damaged particles and outputting the simulation results.

[0012] Preferably, the determination of the interaction between particles, the spin rate tensor and the relative velocity of the particles according to the spatial position and physical parameters of the boundary virtual particles comprises:

[0013] The neighborhood search according to the spatial position and physical parameters of the boundary virtual particles determines the interacting particle pairs and determines the kernel function.

[0014] The boundary application according to the determined particle pairs and the interpolation determine the velocity and stress of the boundary particles.

[0015] The interaction between particles according to the velocity and stress of the boundary particles determines the spin rate tensor and the relative velocity of the particles.

[0016] Preferably, the determination of whether the updated stress state exceeds the failure surface, thereby determining whether the particle is damaged, comprises:

[0017] When the maximum tensile stress of a particle on the failure surface is greater than the preset tensile strength of the rock, the particle is determined to be a tensile damage particle.

[0018] When the maximum tensile stress of a particle on the failure surface is not greater than the tensile strength of the rock, the yield criterion of the particle is calculated.

[0019] When the yield criterion of the particle is zero, the particle is determined to be a shear damage particle.

[0020] Preferably, the yield criterion is

[0021] Where σ is the total stress tensor, c is the cohesion, p(σ) is the hydrostatic stress, J2(s) is the deviatoric stress, and η and ζ are shape parameters controlling the yield surface under the DP criterion.

[0022] Preferably, the marking of the damaged particles and the output of the simulation results comprise:

[0023] The particles marked as damaged are corrected by a preset damage parameter to cut off their connection with surrounding particles and are removed from the calculation domain.

[0024] The damaged particles are filtered using a threshold function in open source software to simulate crack propagation of the damaged particles.

[0025] Preferably, the control equations of the particle system comprise continuity equations and momentum equations.

[0026] The continuity equation is

[0027] The momentum equation is

[0028] Where i represents the i-th particle being calculated; j represents the j-th neighboring particle of the i-th particle within the particle support domain; N is the total number of neighboring particles of the i-th particle; α, β = 1, 2, 3 represent integer indices of two spatial directions, ρ i represents the density of the i-th particle, t is the calculation time, m j represents the mass of the j-th neighboring particle of the i-th particle, v β ij represents the velocity of the j-th neighboring particle of the i-th particle in the β spatial direction, x β i represents the spatial position coordinate of the i-th particle in the β spatial direction, Π ij represents the artificial viscosity term of the j-th neighboring particle of the i-th particle; W ij，β represents the smoothing kernel function of the j-th neighboring particle of the i-th particle in the β spatial direction, v α i represents the velocity of the i-th particle in the α spatial direction, ρ j represents the density of the j-th neighboring particle of the i-th particle, is the stress tensor of the i-th particle in two spatial directions, represents the stress tensor of the j-th neighboring particle of the i-th particle in two spatial directions, I is the total number of particles.

[0029] Preferably, the spin rate tensor is

[0030] The relative velocity is

[0031] Where α, β = 1, 2, 3 represent integer indices of two spatial directions, v α represents the velocity in the α spatial direction, v β represents the velocity in the β spatial direction, x β represents the spatial position coordinate in the β spatial direction, x α represents the spatial position coordinate in the α spatial direction.

[0032] Another embodiment of the present application provides a device for simulating underground space excavation, the device comprising:

[0033] A modeling module is configured to discretize the soil in the region to be simulated into a particle system using a smoothed fluid dynamics method, determine the spatial position and physical parameters of a boundary virtual particle.

[0034] An interaction module is configured to determine the interaction between particles according to the spatial position and physical parameters of the boundary virtual particles, determine the spin rate tensor and relative velocity of the particles;

[0035] A stress module is configured to calculate the stress rate of the soil particles according to the spin rate tensor and relative velocity of the particles, and in combination with the particle constitutive relation;

[0036] An integration module is configured to time-integrate the stress rate to update the stress state of each particle;

[0037] A judgment module is configured to determine whether the updated stress state exceeds the failure surface, so as to determine whether the particle is damaged;

[0038] A result module is configured to mark the damaged particle and output the simulation result.

[0039] The embodiment of the present application further provides a terminal device, including a processor, a memory and a computer program stored in the memory and configured to be executed by the processor, when the processor executes the computer program, the underground space excavation simulation method in any one of the above embodiments is realized.

[0040] The embodiment of the present application further provides a computer readable storage medium, the computer readable storage medium includes a stored computer program, wherein, when the computer program runs, the device where the computer readable storage medium is located executes the underground space excavation simulation method in any one of the above embodiments.

[0041] The underground space excavation simulation method, device, equipment and medium provided by the present application adopt the smooth fluid dynamics method to disperse the soil in the region to be simulated into a particle system, determine the spatial position and physical parameters of the boundary virtual particles, determine the interaction between particles according to the spatial position and physical parameters of the boundary virtual particles, determine the spin rate tensor and relative velocity of the particles, calculate the stress rate of the soil particles according to the spin rate tensor and relative velocity of the particles, in combination with the particle constitutive relation, time-integrate the stress rate to update the stress state of each particle, determine whether the updated stress state exceeds the failure surface, so as to determine whether the particle is damaged, mark the damaged particle and output the simulation result. The present application can simulate under complex geological conditions and non-continuous deformation scenes, and has high precision and stability, which can guarantee the safety of the construction process. BRIEF DESCRIPTION OF DRAWINGS

[0042] Figure 1 is a flowchart of the underground space excavation simulation method provided by the embodiment of the present application;

[0043] Figure 2 is a modeling schematic diagram of the particle system provided by the embodiment of the present application;

[0044] Figure 3 is another flowchart of a method for simulating underground space excavation provided by an embodiment of the present application;

[0045] Figure 4 is a structural diagram of a waste coal mine expansion excavation model provided by an embodiment of the present application;

[0046] Figure 5 is a diagram of roof evolution and stress distribution in a rock mass excavation process provided by an embodiment of the present application;

[0047] Figure 6 is a diagram of fracture development in different rock mass damage degrees provided by an embodiment of the present application;

[0048] Figure 7 is a structural diagram of a device for simulating underground space excavation provided by an embodiment of the present application;

[0049] Figure 8 is a structural diagram of a terminal device provided by an embodiment of the present application. DETAILED DESCRIPTION

[0050] The technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only some of the embodiments of the present application, but not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative work fall within the protection scope of the present application.

[0051] Referring to Figure 1 is a flowchart of a method for simulating underground space excavation provided by an embodiment of the present application, and the method comprises steps S1-S6:

[0052] In step S1, a soil body in a region to be simulated is discretized into a particle system by using a smooth fluid dynamics method, and the spatial position and physical parameters of a boundary virtual particle are determined.

[0053] In step S2, the interaction between particles is determined according to the spatial position and physical parameters of the boundary virtual particle, and the spin rate tensor and relative velocity of the particles are determined.

[0054] In step S3, the stress rate of the soil body particles is calculated according to the spin rate tensor and relative velocity of the particles and in combination with the particle constitutive relation.

[0055] In step S4, the stress state of each particle is updated by time integration of the stress rate.

[0056] In step S5, it is determined whether the updated stress state exceeds a failure surface, so as to determine whether the particle is damaged.

[0057] Step S6, mark the damaged particles, and output the simulation result.

[0058] In the implementation of the embodiment, the soil in the region to be simulated is discretized into a particle system by using the SPH (Smoothed Particle Hydrodynamics) method. In the SPH, each particle represents a micro unit of the soil and has physical attributes such as mass, position, and velocity. For the setting of the boundary virtual particles, the virtual particles are generated outside the physical boundary within a certain range according to the geometric boundary conditions of the simulation region. The spatial positions of the virtual particles are strictly determined according to the boundary shape and the extension rules, and the physical parameters (such as density and elastic modulus) of the virtual particles are usually assigned according to the attributes of the adjacent real particles and the boundary conditions. For example, at the fixed boundary, the velocity of the virtual particle is set to zero, and the density and elastic modulus are consistent with the immediately adjacent real particle, so as to simulate the constraint effect of the boundary on the soil.

[0059] Referring to Figure 2 is a modeling schematic diagram of the particle system provided by the embodiment of the application. The motion and properties of the research object are approximately described by dividing the research object in the problem domain into many micro particles and establishing a mathematical model on the particles. When the SPH method is used for modeling, the soil in the region to be simulated is discretized into a particle system according to the physical and geometric parameters of the region to be simulated. Each particle represents a part of the research object and has its position, velocity, density, and other attributes. The boundary conditions are established to introduce the boundary virtual particles. The spatial positions and physical parameters of the boundary virtual particles are determined according to the types of the simulation boundary (such as the fixed wall boundary, friction boundary, or free boundary of the soil), and the part of the soil interacting with the boundary is processed.

[0060] According to the spatial positions and physical parameters of the boundary virtual particles and the real particles, the interaction between the particles is determined by using the kernel function. The kernel function defines the range and weight of the influence domain of the particles, and the interaction exists only between the particles in the influence domain. By calculating the distance between the particles and the relationship with the kernel function, the interaction force between the particles is obtained, and the spin rate tensor and the relative velocity of the particles are determined. The spin rate tensor reflects the rotation characteristics of the particle itself, and the relative velocity embodies the motion difference between the particles. In the calculation process, based on the Newtonian motion law and the principles of fluid mechanics, the mass, position, and velocity information of the particles are combined to gradually and accurately solve these parameters by iterative calculation.

[0061] According to the spatial positions and physical parameters of the boundary virtual particles and real particles, the interaction between particles is determined by using a kernel function. The kernel function defines the range and weight of the particle influence domain, and interaction only exists between particles within the influence domain. By calculating the distance between particles and the relationship with the kernel function, the interaction force between particles is obtained, and the spin rate tensor and relative velocity of the particles are determined. The spin rate tensor reflects the rotation characteristics of the particles themselves, while the relative velocity embodies the motion difference between particles. In the calculation process, based on Newton's law of motion and the principles of fluid mechanics, combined with the mass, position and velocity information of the particles, these parameters are gradually accurately solved through iterative calculation.

[0062] According to the determined spin rate tensor and relative velocity of the particles, the stress rate of the soil particles is calculated by combining the constitutive relationship of the particles (such as elastic-plastic constitutive model, viscoelastic constitutive model, etc.). The constitutive relationship describes the relationship between stress and strain of the soil material under force, and the motion state parameters of the particles are substituted into the constitutive equation to obtain the stress rate through mathematical operation. For example, in the elastic-plastic constitutive model, considering the elastic deformation stage and plastic deformation stage of the soil, the corresponding stress rate change is calculated according to the strain increment and yield condition of the particles, to reflect the mechanical response of the soil in the process of force.

[0063] The stress rate is time-integrated to update the stress state of each particle. In the process of time discretization, the finite difference method or other numerical integration methods are usually used to divide the time into multiple small time steps. In each time step, the stress state at the current time is calculated by integrating the stress state at the previous time and the stress rate calculated at the current time.

[0064] The updated stress state is judged to determine whether it exceeds the failure surface, which is usually determined by the strength criterion of the soil such as the Mohr-Coulomb strength criterion, the Drucker-Prager strength criterion, etc., to determine whether the particle is damaged. If the stress state of the particle meets the failure condition, i.e. exceeds the ultimate stress range that the soil material can withstand, it is determined that the particle is damaged. For example, under the Mohr-Coulomb strength criterion, by comparing the relationship between the shear stress and normal stress of the particle and the internal friction angle and cohesion of the material, it is judged whether the particle reaches the failure state.

[0065] The damaged particles are marked for subsequent analysis and visual display. The marking method can be to assign a specific identifier or change the color, shape and other attributes of the particles. Finally, the particle state information (including position, velocity, stress, damage condition, etc.) in the entire simulation process is sorted and output, and the simulation results are presented through data files or visualization software, helping researchers to intuitively understand the mechanical behavior and failure process of the soil under different working conditions.

[0066] By discretizing soil into a particle system and precisely setting boundary virtual particles and calculating the interaction between particles, the mechanical response of soil under complex stress and deformation conditions can be accurately captured. Compared with the traditional continuum mechanics method, this scheme can better handle the nonlinear behavior of soil such as large deformation, fracture and flow, and provide more reliable basis for the design and analysis of soil engineering. For example, when simulating the liquefaction process of soil under earthquake action or the slope instability problem, the stress distribution, deformation development and failure path inside the soil can be clearly shown. By relying on the SPH method for numerical simulation, this method effectively overcomes the technical bottleneck of grid distortion and calculation divergence often encountered in traditional grid methods when dealing with large deformation problems, and is particularly suitable for simulating complex engineering problems in the fields of tunnel excavation, mineral exploitation and pile construction.

[0067] In another embodiment provided by the application, the determination of the interaction between particles, the determination of the spin rate tensor and the relative velocity of the particles according to the spatial position and the physical parameters of the boundary virtual particles comprises:

[0068] According to the spatial position and the physical parameters of the boundary virtual particles, the neighborhood search is performed to determine the particle pairs for interaction, and the kernel function is determined;

[0069] According to the determined particle pairs, the boundary is applied, and the velocity and stress of the boundary particles are determined by interpolation;

[0070] According to the velocity and stress of the boundary particles, the interaction between particles is determined, and the spin rate tensor and the relative velocity of the particles are calculated.

[0071] In the implementation of the present embodiment, refer to Figure 3 Fig. 2 is another flowchart of the simulation method of underground space excavation provided by the embodiment of the application, and specifically:

[0072] First, the soil is discretized, that is, the particle discretization of the problem domain and the initialization of physical information are performed.

[0073] Then, the SPH calculation cycle is executed to determine the spin rate tensor and the relative velocity of the particles, and specifically:

[0074] Based on the neighborhood search, the interaction particle pairs are determined, the kernel function and the kernel function gradient field are calculated; based on the spatial position information of the boundary virtual particles, efficient data structures such as spatial hash table and octree are constructed, and the simulation area is spatially divided. Taking each boundary virtual particle as the center, according to the set search radius (usually related to the influence range of the kernel function), the real particles and other virtual particles in its neighborhood range are quickly located, so as to determine the interaction particle pairs. For example, in the spatial hash table, by calculating the hash value corresponding to the particle coordinates, the particles in the same grid or adjacent grid are quickly found, greatly reducing the search time. According to the material properties of the soil and the simulation requirements, a suitable kernel function is selected, such as Gaussian kernel function, cubic spline kernel function, etc. By analyzing the distance distribution and physical parameter difference between particle pairs, combined with empirical formula or optimization algorithm, the key parameters of the kernel function are determined, such as the smoothing length. The smoothing length determines the influence range of the kernel function, and affects the strength and range of the interaction between particles. Based on the selected kernel function and its parameters, the derivative operation of the kernel function is carried out, and the kernel function gradient field is calculated. The kernel function gradient field describes the change trend of the kernel function in space, and reflects the change of the interaction strength between particles with distance. In actual calculation, numerical differentiation methods such as finite difference method can be used to discretize the position of each particle pair, and the gradient value of the kernel function at each point is solved, so as to construct the kernel function gradient field of the entire simulation area.

[0075] Boundary application, interpolation to determine boundary particle velocity, stress and other information.

[0076] According to the actual boundary conditions of the simulation scene (such as fixed boundary, free boundary, stress boundary, etc.), the boundary conditions are applied to the determined particle pairs. At the fixed boundary, the displacement of the boundary virtual particles and the related real particles is set to zero, limiting their movement in a specific direction; at the stress boundary, the boundary particles are given corresponding stress values according to the design requirements. Using weighted average interpolation, radial basis function interpolation and other methods, the velocity and stress of the boundary particles are determined according to the velocity and stress of the interaction particle pairs. Taking weighted average interpolation as an example, the velocity of the boundary particle is calculated as the distance between the boundary particle and the adjacent particles according to the inverse of the distance between the particles as the weight. In this way, the parameters of the boundary particles can reasonably reflect the state of the surrounding particles, ensuring the accurate transmission of the boundary conditions.

[0077] The inter-particle interaction is determined, the spin rate tensor and the relative velocity are calculated according to the particle velocity gradient. According to the velocity and stress determined by the boundary particles, the inter-particle interaction force is calculated by combining Newton's law of motion and the basic equations of fluid mechanics. Considering the influence of particle mass, acceleration and nuclear function gradient field on the interaction force, the size and direction of the interaction force between particles are solved through the force balance equation and the momentum conservation equation. Based on the calculated inter-particle interaction force, the spin rate tensor and the relative velocity of the particle are further derived. The spin rate tensor reflects the rotation characteristics of the particle itself, which is calculated by analyzing the torque and moment of inertia acting on the particle; the relative velocity is obtained by calculating the difference between the velocity vectors of the particles. In the calculation process, an iterative algorithm is used to gradually update the position, velocity and other parameters of the particles, ensuring the accuracy and convergence of the calculation results, so as to accurately describe the motion and mechanical response of the particles under the influence of boundary conditions.

[0078] In another embodiment provided by the application, the damaged particles are marked, and the simulation results are output, including:

[0079] The particles marked as damaged are modified by a preset damage parameter to correct the nuclear function in the contact pair, cut off the connection with the surrounding particles, and removed from the calculation domain;

[0080] The damaged particles are filtered using a threshold function in open source software to simulate crack propagation of the damaged particles.

[0081] In the specific implementation of the embodiment, in the calculation of the elastic-plastic fracture model, the stress rate of the soil particles is calculated according to the constitutive relation of the particles, and specifically:

[0082] In classical elastic-plastic mechanics, the strain of a material is usually divided into reversible elastic strain and irreversible plastic strain. In order to calculate the stress-strain relationship of soil, concepts such as yield function and plastic potential function need to be introduced.

[0083] A general elastic-plastic constitutive model is used to calculate the total stress tensor, and the constitutive model equation describes the relationship between stress, strain and internal variables of the material state. The expression of the total stress tensor σ and the plasticity index λ are:

[0084]

[0085] Wherein, wherein, σ αβ =s αβ +1 / 3σ γγ δ αβ ; And are the deviatoric stress rate tensor and the deviatoric strain rate tensor respectively; K is the bulk modulus; G is the shear modulus; δ αβis the Kronecker function symbol, δ αβ = 1 when α ≠ β, δ αβ = 0; g is the plastic potential function defined with respect to the velocity.

[0086] Time integration is performed to update the physical state of each particle.

[0087] Then, it is judged whether the stress state of the particle exceeds the failure surface;

[0088] The brittle fracture process of rock can be characterized by a brittle damage model. When the stress state of a point inside the rock satisfies the failure criterion, the point fails due to damage. In the simulation process, first, the stress state information of each particle on the failure surface is obtained, and the maximum tensile stress value is extracted. The maximum tensile stress is compared with the preset tensile strength of the rock. If the maximum tensile stress of a particle on the failure surface is greater than the preset tensile strength of the rock, according to the judgment rule, the particle is directly judged to be a tensile damage particle, and is marked accordingly for subsequent analysis and visualization processing. This step quickly screens out the particles damaged due to tensile action, providing a basis for subsequent analysis.

[0089] When the maximum tensile stress of a particle on the failure surface is not greater than the tensile strength of the rock, it indicates that the particle has not been damaged by tension, but may still be damaged by shear. At this time, the yield criterion of the particle needs to be calculated. Commonly used yield criteria include the Mohr-Coulomb criterion, the Drucker-Prager criterion, etc. The appropriate criterion is selected according to the characteristics of the rock material in the simulation scene and the actual needs.

[0090] In the calculation, the normal stress and shear stress suffered by the particle are combined with the cohesion and internal friction angle parameters of the rock material to calculate the yield criterion value of the particle by substituting the formula. In actual operation, the stress state data of the particle can be obtained by finite element analysis, numerical simulation calculation, etc. to ensure the accuracy of the calculation.

[0091] The calculated yield criterion value of the particle is compared with zero. When the yield criterion of the particle is zero, it means that the stress suffered by the particle has reached the shear failure limit of the rock material, satisfying the shear failure condition. At this time, the particle is judged to be a shear damage particle and is marked. In this way, further damage judgment of the particles that have not been damaged by tension is completed, and the damage type of the particles is determined comprehensively.

[0092] In another embodiment provided by the application, the yield criterion is

[0093] where σ is the total stress tensor, c is the cohesion, p(σ) is the hydrostatic stress; J2(s) is the deviatoric stress, and η and ζ are shape parameters controlling the yield surface under the DP criterion.

[0094] The D-P yield surface is tangent to the Mohr-Coulomb yield surface in the pi plane, and the relevant parameters are as follows:

[0095]

[0096] Wherein, φ is the internal friction angle; α φ and k c are DP constants.

[0097] In still another embodiment provided by the application, the particle damaged is marked, and the output simulation result comprises:

[0098] The particle damaged is marked, the kernel function in the contact pair is corrected through a preset damage parameter, the connection with surrounding particles is cut off, and the particle is removed from the calculation domain;

[0099] The damaged particle is filtered by using a threshold function in open source software, and crack propagation of the damaged particle is simulated.

[0100] In the embodiment, the particle damaged is marked, the kernel function in the contact pair is corrected through a damage parameter mu, the connection with surrounding particles is cut off, and the particle is removed from the calculation domain.

[0101]

[0102] In each loop process of the program, the fracture state of each particle is judged, when the particle does not satisfy the damage condition, mu is recorded as 1, and when the particle satisfies the damage condition, mu is recorded as 0.

[0103] The algorithm simulation system calculates field variable information such as stress, strain, displacement, velocity and acceleration of each particle by using a smoothed particle hydrodynamics method, the obtained information is digital text information, the digital text information is converted into a binary VTK file to realize visualization. The particles in different partitions are output through particle type control, so that material output classification can be realized, and finally binary files of soil and boundary particles can be output, and the binary files can be displayed in some general open source software such as Paraview. In each calculation process, the damaged particle is marked through a self-defined variable, so that the output result in post-processing can be controlled, in the example, the threshold function of Paraview is used for filtering based on the damage parameter mu, so that the display of the damaged particle is controlled, and then the crack propagation phenomenon of the rock mass is simulated.

[0104] In still another embodiment provided by the application, the control equation of the particle system comprises a continuity equation and a momentum equation;

[0105] The soil deformation is calculated in the SPH framework, and the soil adopts a subplasticity constitutive model to update the stress state. The SPH is a Lagrangian meshless method, and the calculation domain is discretized into a series of material particles, which carry field variable information such as density, velocity, stress, strain, etc. By solving the control equation, the field variable information corresponding to different time steps can be obtained. The continuity equation and the momentum equation are as follows:

[0106]

[0107] After discretization, the following can be obtained:

[0108] The continuity equation is

[0109] The momentum equation is

[0110] Where i represents the i th particle being calculated; j represents the j th adjacent particle of the i th particle in the particle support domain; N is the total number of adjacent particles of the i th particle; alpha, beta = 1, 2, 3, represent integer indexes of two spatial directions, rho i represents the density of the i th particle, t is the calculation time, m j represents the mass of the j th adjacent particle of the i th particle, v β ij represents the velocity of the j th adjacent particle of the i th particle in the beta spatial direction, x β i represents the spatial position coordinate of the i th particle in the beta spatial direction, Pi ij represents the artificial viscosity term of the j th adjacent particle of the i th particle; W ij，β represents the smoothing kernel function of the j th adjacent particle of the i th particle in the beta spatial direction, v α i represents the velocity of the i th particle in the alpha spatial direction, rho j represents the density of the j th adjacent particle of the i th particle, is the stress tensor of the i th particle in two spatial directions, represents the stress tensor of the j th adjacent particle of the i th particle in two spatial directions, and I is the total number of particles.

[0111] In another embodiment provided by the application, the spin rate tensor is

[0112] The relative velocity is

[0113] Where alpha, beta = 1, 2, 3, represent integer indexes of two spatial directions, v α represents the velocity in the alpha spatial direction, vβ x represents the velocity in the β space direction. β x represents the spatial coordinates of the location in the β-space direction. α This represents the spatial position coordinates in the α-space direction.

[0114] This application's scheme, by embedding the Drucker-Prager yield criterion and particle fracture model into the SPH numerical framework, can accurately capture the large deformation behavior of rock masses during excavation, precisely predict the stability and fracture evolution of surrounding rock under complex geological conditions, and provide an efficient numerical tool for predicting the mechanical response of surrounding rock in coal mine underground reservoir expansion excavation projects. It can significantly improve the safety and feasibility assessment of underground reservoir expansion excavation projects. Through numerical simulation based on the SPH method, this method effectively overcomes the technical bottlenecks of mesh distortion and computational non-convergence often encountered by traditional mesh methods when dealing with large deformation problems. It achieves accurate simulation of the entire process of excavation, fracture development, and collapse of abandoned coal mine rock strata, and can further combine the spatial variability characteristics of soil materials to use Monte Carlo simulation technology for reliability analysis of mining engineering, tunnel construction, and underground space expansion problems, providing a scientific basis for engineering decisions. This application's scheme can consider the impact of construction excavation disturbance on rock masses under different damage levels, accurately simulate the fracture initiation, propagation, and rock mass instability behavior during excavation, and further optimize surrounding rock stability analysis and construction risk assessment.

[0115] In another embodiment of the present invention, an expansion excavation model for an abandoned coal mine is established based on on-site survey data. The model is divided into seven geological strata, see [link to relevant documentation]. Figure 4 This is a schematic diagram of the structure of the abandoned coal mine expansion excavation model provided in this embodiment of the invention. The soil and rock mass is simulated using an elastoplastic model combined with a particle fracture model. The material parameters are shown in Table 1, and the boundaries are simulated using virtual particles.

[0116] Table 1 Material Parameters

[0117]

[0118] According to the algorithm provided by this invention, the roof evolution and stress distribution during the excavation of intact rock mass in abandoned coal mine areas are predicted. (See also...) Figure 5 This is a schematic diagram of the roof evolution and stress distribution during the rock mass excavation process provided in an embodiment of the present invention.

[0119] Affected by geological structure and mining disturbance, the roof rock mass has initial damage, and there is a certain safety risk in the secondary excavation of large chamber in the rock mass. The patent can also consider the historical damage of the coal mine rock mass by the multi-stratum rock mass strength reduction method, simulate the large deformation problem of the secondary excavation of the underground reservoir roof, and apply 0.5, 0.7, 0.9 and 1.0 times strength reduction to the friction angle, cohesion, and tensile strength parameters of the siltstone, fine-grained sandstone, coarse-grained sandstone, and medium-grained sandstone from bottom to top for the slightly damaged rock mass (simulating weak mining influence). The corresponding reduction coefficients are 0.3, 0.6, 0.9 and 1.0 for the severely damaged rock mass (simulating strong mining influence). The crack development of rock mass with different damage degrees is as shown in Figure 6 , which is a crack development schematic diagram of rock mass with different damage degrees provided by the embodiment of the present application.

[0120] The scheme can consider the construction disturbance influence of the rock mass under different damage degrees, and accurately simulate the crack initiation, expansion and rock mass instability behavior in the excavation process.

[0121] The embodiment of the present application also provides an underground space excavation simulation device, as shown in Figure 7 , which is a structural schematic diagram of the underground space excavation simulation device provided by the embodiment of the present application, and the device comprises:

[0122] A modeling module is used to discretize the soil in the region to be simulated into a particle system by using a smooth fluid dynamics method, determine the spatial position and physical parameters of the boundary virtual particles;

[0123] An interaction module is used to determine the interaction between particles according to the spatial position and physical parameters of the boundary virtual particles, determine the spin rate tensor and relative velocity of the particles;

[0124] A stress module is used to calculate the stress rate of the soil particles in combination with the particle constitutive relation according to the spin rate tensor and relative velocity of the particles;

[0125] An integration module is used to time-integrate the stress rate and update the stress state of each particle;

[0126] A judgment module is used to judge whether the updated stress state exceeds the failure surface, so as to determine whether the particle is damaged;

[0127] A result module is used to mark the damaged particles and output the simulation result.

[0128] It should be noted that the underground space excavation simulation device provided by the embodiment of the present application can execute the underground space excavation simulation method described in any of the above embodiments, and the specific functions of the underground space excavation simulation device will not be repeated here.

[0129] , which is a structural schematic diagram of the underground space excavation simulation device provided by the embodiment of the present application, and the device comprises: Figure 8Fig. 1 is a structural schematic diagram of a terminal device according to an embodiment of the present application. The terminal device includes a processor, a memory, and a computer program stored in the memory and executable on the processor, for example, an underground space excavation simulation program. The processor implements the steps of each of the underground space excavation simulation method embodiments described above when executing the computer program, for example, the steps S1-S6 shown in Fig. 1. Alternatively, the processor implements the functions of each module in each of the device embodiments described above when executing the computer program. Figure 1

[0130] For example, the computer program can be divided into one or more modules / units, which are stored in the memory and executed by the processor to complete the present application. The one or more modules / units can be a series of computer program instruction segments capable of completing a function, which are used to describe the execution process of the computer program in the terminal device. For example, the computer program can be divided into various modules, and the specific functions of the modules are not described again.

[0131] The terminal device can be a desktop computer, a notebook computer, a palm computer, a cloud server, and the like. The terminal device can include, but is not limited to, a processor and a memory. Those skilled in the art can understand that the schematic diagram is only an example of the terminal device, and does not limit the terminal device, which can include more or fewer components than those shown, or combine certain components, or different components, for example, the terminal device can also include an input / output device, a network access device, a bus, and the like.

[0132] The processor can be a central processing unit (CPU), and can also be other general-purpose processors, digital signal processors (DSP), application specific integrated circuits (ASIC), field programmable gate arrays (FPGA) or other programmable logic devices, discrete gates or transistor logic devices, discrete hardware components, or the like. The general-purpose processor can be a microprocessor or the processor can also be any conventional processor or the like. The processor is a control center of the terminal device, which connects all parts of the terminal device through various interfaces and lines.

[0133] ​The memory can be used to store the computer program and / or modules, and the processor realizes various functions of the terminal device by running or executing the computer program and / or modules stored in the memory, and calling data stored in the memory. The memory can mainly include a program storage area and a data storage area, wherein the program storage area can store an operating system, at least one application program required by a function (such as a sound playing function, an image playing function, etc.), etc.; and the data storage area can store data created according to the use of the mobile phone (such as audio data, a phone book, etc.), etc. In addition, the memory can include a high-speed random access memory, and can also include a non-volatile memory, for example, a hard disk, a memory, a plug-in hard disk, a smart media card (SMC), a secure digital (SD) card, a flash card, at least one disk storage device, a flash memory device, or other volatile solid-state memory devices.

[0134] The modules / units integrated in the terminal device, if realized in the form of software function units and sold or used as independent products, can be stored in a computer readable storage medium. Based on this understanding, all or part of the processes in the above-mentioned embodiment methods can also be completed by a computer program instructing related hardware, and the computer program can be stored in a computer readable storage medium. When the processor executes the computer program, the steps of the above-mentioned various method embodiments can be realized. The computer program includes computer program code, which can be in the form of code, object code, executable files or some intermediate forms, etc. The computer readable medium can include any entity or device capable of carrying the computer program code, recording medium, U disk, mobile hard disk, magnetic disk, optical disk, computer memory, read-only memory (ROM), random access memory (RAM), electrical carrier signal, telecommunication signal and software distribution medium, etc.

[0135] The above is the preferred embodiment of the present application, and it should be pointed out that for those skilled in the art, without departing from the principles of the present application, a number of improvements and refinements can be made, which are also considered within the scope of protection of the present application.

Claims

1. A method of simulating excavation of an underground space, characterized by, The method comprises: Discretizing the soil in the region to be simulated into a particle system by using a smooth fluid dynamics method, determining the spatial position and physical parameters of boundary virtual particles; Determining the interaction between particles, the spin rate tensor and the relative velocity of the particles according to the spatial position and physical parameters of the boundary virtual particles; Calculating the stress rate of the soil particles according to the spin rate tensor and the relative velocity of the particles and combining the particle constitutive relation; Time-integrating the stress rate to update the stress state of each particle; Determining whether the updated stress state exceeds the failure surface, thereby determining whether the particle is damaged; Marking the damaged particles and outputting the simulation result.

2. The underground space excavation simulation method according to claim 1, characterized by, The determination of the interaction between particles, the spin rate tensor and the relative velocity of the particles according to the spatial position and physical parameters of the boundary virtual particles comprises: Performing neighborhood search according to the spatial position and physical parameters of the boundary virtual particles, determining the particle pairs for interaction and determining the kernel function; Determining the velocity and stress of the boundary particles by interpolation according to the boundary application of the determined particle pairs; Determining the interaction between particles, calculating the spin rate tensor and the relative velocity of the particles according to the velocity and stress of the boundary particles.

3. The underground space excavation simulation method according to claim 1, characterized by, The determination of whether the updated stress state exceeds the failure surface, thereby determining whether the particle is damaged, comprises: When the maximum tensile stress of a particle on the failure surface is greater than the preset rock tensile strength, the particle is determined to be a tensile damage particle; When the maximum tensile stress of a particle on the failure surface is not greater than the rock tensile strength, the yield criterion of the particle is calculated; When the yield criterion of the particle is zero, the particle is determined to be a shear damage particle.

4. The underground space excavation simulation method according to claim 3, characterized by, The yield criterion is Wherein, σ is the total stress tensor, c is the cohesion, p(σ) is the hydrostatic stress; J2(s) is the deviatoric stress, η and ζ are the shape parameters of the yield surface under the DP criterion.

5. The underground space excavation simulation method according to claim 1, characterized by, The marking of the damaged particles and the output of the simulation result comprise: The particles marked as damaged are corrected by a preset damage parameter to modify the kernel function in the contact pair, cut off the connection with the surrounding particles and removed from the calculation domain; The threshold function is used in the open source software to filter the damaged particles and simulate the crack propagation of the damaged particles.

6. The underground space excavation simulation method according to claim 1, wherein, The control equation of the particle system comprises a continuity equation and a momentum equation. The continuity equation is The momentum equation is where i denotes the i-th particle being calculated; j denotes the j-th neighboring particle of the i-th particle within the particle support domain; N is the total number of neighboring particles of the i-th particle; a, b = 1, 2, 3 denote integer indices of the two spatial directions; p i denotes the density of the i-th particle, t is the calculation time, m j denotes the mass of the j-th neighboring particle of the i-th particle, v β ij denotes the velocity of the j-th neighboring particle of the i-th particle in the b spatial direction, x β i denotes the spatial position coordinate of the i-th particle in the a spatial direction, Π ij denotes the artificial viscosity term of the j-th neighboring particle of the i-th particle; W ij，β denotes the smoothing kernel function of the j-th neighboring particle of the i-th particle in the b spatial direction, v α i denotes the velocity of the i-th particle in the a spatial direction, p j denotes the density of the j-th neighboring particle of the i-th particle, is the stress tensor of the i-th particle in the two spatial directions, denotes the stress tensor of the j-th neighboring particle of the i-th particle in the two spatial directions, I is the total number of particles.

7. The underground space excavation simulation method according to claim 1, wherein, The spin rate tensor is The relative speed is where a, b = 1, 2, 3, represent integer indices of two spatial directions, v α represents velocity in a spatial direction a, v β represents velocity in a spatial direction b, x β represents spatial position coordinate in a spatial direction b, x α represents spatial position coordinate in a spatial direction a.

8. An underground space excavation simulation apparatus characterized by comprising: The device comprises: A modeling module for discretizing the soil in the region to be simulated into a particle system by using a smooth fluid dynamics method, determining the spatial position and physical parameters of boundary virtual particles; An interaction module for determining the interaction between particles, the spin rate tensor and the relative velocity of the particles according to the spatial position and physical parameters of the boundary virtual particles; A stress module for calculating the stress rate of the soil particles according to the spin rate tensor and the relative velocity of the particles and combining the particle constitutive relation; An integration module for time-integrating the stress rate to update the stress state of each particle; A judgment module for determining whether the updated stress state exceeds the failure surface, thereby determining whether the particle is damaged; A result module for marking the damaged particles and outputting the simulation result.

9. A terminal device, comprising: An underground space excavation simulation method as claimed in any one of claims 1 to 7, comprising a processor, a memory, and a computer program stored in the memory and configured to be executed by the processor, the computer program being implemented when executed by the processor.

10. A computer-readable storage medium, characterized in that, A computer readable storage medium comprising a stored computer program, wherein the computer readable storage medium, when the computer program is run, controls a device in which the computer readable storage medium is located to perform an underground space excavation simulation method as claimed in any one of claims 1 to 7.

Images