Near-field dynamic simulation method and system for dynamic fracture of underground cavern surrounding rock under thermo-hydraulic coupling

Abstract

The present application relates to the technical field of underground engineering surrounding rock mechanics simulation, and specifically discloses a near-field dynamics simulation method and system for underground cavern surrounding rock thermal-hydraulic coupling dynamic fracture, which comprises the following steps: using RSCFS algorithm to control randomization of the obtained surrounding rock parameter information to obtain a heterogeneous parameter field; building a dynamic thermal-hydraulic coupling near-field dynamics theoretical model, inputting the heterogeneous parameter field into the dynamic thermal-hydraulic coupling near-field dynamics theoretical model, and obtaining an initial material point cloud through numerical discretization; solving the initial material point cloud through an explicit numerical solver to obtain a discrete damage point cloud; analyzing the field variable dynamic data and crack propagation characteristics of multiple physical fields to obtain the simulation result of the underground cavern surrounding rock thermal-hydraulic coupling dynamic fracture. The near-field dynamics simulation method of the present application realizes high-precision and high-efficiency simulation of the underground cavern surrounding rock thermal-hydraulic coupling dynamic fracture, and can accurately reproduce the crack morphology and multi-field evolution characteristics.

Description

Technical Field

[0001] This invention relates to the field of simulation technology of surrounding rock mechanics in underground engineering, and in particular to a near-field dynamic simulation method and system for dynamic fracture of surrounding rock in underground caverns caused by hydrothermal coupling. Background Technology

[0002] In major underground engineering projects such as deep resource extraction, deep-earth laboratory construction, and geothermal energy development, the surrounding rock is often subjected to a multi-field coupled environment of high ground stress, high pore water pressure, and high geothermal temperature (THM). It is also susceptible to dynamic impact loads such as drilling and blasting, and earthquakes, which can induce rapid crack initiation and propagation, or even catastrophic failure, seriously threatening engineering safety. Therefore, accurate simulation of the THM-coupled dynamic fracture behavior of the surrounding rock in underground caverns under impact loads is a core key to the stability assessment of deep underground engineering projects.

[0003] Currently, numerical simulation methods for dynamic rock fracture mainly include the discrete element method (DEM), finite element method (FEM), boundary element method (BEM), and traditional near-field dynamics methods. However, all of these methods have significant technical limitations: the DEM can describe the separation and movement of granular media, but it struggles to handle continuous medium deformation and complex fluid-thermal coupling; the FEM and its extensions can describe multi-field coupling in continuous media, but the governing equations based on the continuity assumption exhibit singularities at the crack tip, and the algorithms for dynamic multi-crack propagation are complex and computationally expensive; while traditional near-field dynamics methods can naturally describe the fracture process from continuous to discontinuous materials and avoid stress singularities at the crack tip, Existing research mainly focuses on single mechanical fields or simple coupling scenarios, lacking a complete theoretical framework that can accurately capture the fully coupled THM effect under high strain rate impact. Furthermore, it suffers from three major problems: First, the parameter field randomization method easily disrupts the cylindrical / radial symmetry of the surrounding rock and the spatial correlation of the physical field, leading to distorted parameter characterization. Second, the characteristic timescales of elastic wave propagation, thermal diffusion, and fluid seepage in the THM coupling system differ greatly, and a single time step results in extremely low computational efficiency. Third, the discrete damage data output from simulations is difficult to directly convert into a continuous crack network, making it impossible to achieve quantitative analysis of crack propagation driven by energy.

[0004] In summary, existing simulation methods cannot meet the high-precision, high-efficiency, and quantitative requirements of near-field dynamic simulation of THM coupled dynamic fracture of surrounding rock in deep underground caverns. There is an urgent need to develop a near-field dynamic simulation method and system that adapts to the geometric characteristics and multi-field coupling properties of surrounding rock in underground caverns. Summary of the Invention

[0005] The purpose of this invention is to overcome the shortcomings of the prior art and provide a near-field dynamic simulation method and system for the dynamic fracture of the surrounding rock of underground caverns by hydrothermal coupling. This invention solves the technical problems of existing methods being unable to accurately capture the full coupling effect of THM under high strain rate impact, parameter characterization distortion, low computational efficiency, and difficulty in quantitative analysis of crack networks using discrete data. This invention achieves high-precision, high-efficiency, and quantitative simulation of the dynamic fracture of the surrounding rock of underground caverns by THM coupling.

[0006] To achieve the above objectives, the present invention adopts the following technical solution: On the one hand, this invention provides a near-field dynamic simulation method for the dynamic fracture of surrounding rock in underground caverns caused by hydrothermal coupling, comprising the following steps: S1. The obtained surrounding rock parameter information is processed by the RSCFS radial symmetric correlation field shuffling algorithm in a controllable randomization manner to obtain a heterogeneous parameter field subject to geometric and physical constraints. S2. Construct a dynamic hydrothermal coupling near-field dynamics theoretical model. Input the heterogeneous parameter field into the dynamic hydrothermal coupling near-field dynamics theoretical model and obtain the initial material point cloud characterizing the surrounding rock through numerical discretization calculation. The initial material point cloud includes the initial state of material points and bonds. S3. An explicit numerical solver is constructed using a multi-time-step coupling strategy and a domain decomposition parallel algorithm. The initial material point cloud is solved by the explicit numerical solver in the time domain to obtain a discrete damage point cloud. The discrete damage point cloud includes dynamic data of field variables of multi-physics fields. The dynamic data of field variables of multi-physics fields include the mechanical field, thermal field, and flow field data of the surrounding rock, as well as discrete damage data determined by the damage function. S4. Based on the triangulation method, the discrete damage point cloud is reconstructed into a continuous crack mesh, and the parameters obtained by the dynamic hydrodynamic coupling near-field dynamic theory model are mapped to the continuous crack mesh through weight quantization and multi-field parameter mapping to obtain crack propagation characteristics. S5. Based on the dynamic data of the field variables of the multi-physics field and the crack propagation characteristics, generate the spatiotemporal evolution curves of the multi-physics field, crack morphology diagrams, crack mesh topology parameters and energy evolution reports, respectively, and obtain the simulation results of the dynamic fracture of the surrounding rock of the underground cavern by hydrothermal coupling.

[0007] Furthermore, in step S1, the obtained surrounding rock parameter information is subjected to controllable randomization processing using the RSCFS radial symmetric correlation field shuffling algorithm, specifically including: S11. Obtain surrounding rock parameter information, perform spatial parameterization on the surrounding rock parameter information to obtain a physical parameter field, and perform exponential covariance function and Karhunen-Loève expansion processing on the physical parameter field to generate a random field of initial parameters of the surrounding rock. S12. Perform a three-layer constraint shuffling operation on the random field of the initial parameters of the surrounding rock. The three-layer constraint shuffling operation adopts a layered randomization strategy, and applies geometric and physical constraints in the angular dimension, radial dimension and axial dimension in sequence to obtain the processed random field of the initial parameters of the surrounding rock. Wherein, the angular dimension is to perform sector rotation to maintain cylindrical symmetry, the radial dimension is to employ a distance-constrained Fisher-Yates variant to maintain parametric gradient characteristics, and the axial dimension allows for medium-distance exchanges to simulate interlayer variability; S13. Apply physical limit constraints and gradient smoothing constraints to the processed initial parameter random field of the surrounding rock, and obtain a controllable randomized heterogeneous parameter field through iterative correction. The heterogeneous parameter field includes mechanical field, thermal field and flow field.

[0008] Furthermore, the expression for the physical parameter field is:

[0009] in, For physical parameter fields; For the material domain of the surrounding rock, ,in The radial coordinates of the material point. Let be the angular coordinates of the material point. Let be the axial coordinates of the material point. The radius of the surrounding rock sample is denoted as . The height of the surrounding rock sample; The dimension of the physical parameter field represents the total number of physical parameters carried by each material point; The physical parameter field satisfies , For the position of the matter point The physical parameter vector at that location; The expression for the exponential covariance function is:

[0010] in, It is an exponential covariance function. , Let them be two material points within the physical parameter field; The radial coordinate difference between two material points. This represents the difference in the axial coordinates of two material points; , These are the radial length and axial length of the physical parameter field, respectively. The variance of the physical parameter field; The Karhunen-Loève expanded expression is:

[0011] in, The physical parameters of the random field are reconstructed by Karhunen-Loève expansion; is the mean function of the physical parameter field; M is the number of truncation terms selected according to the calculation accuracy requirements; The eigenvalues ​​are the covariance operators of the physical parameter field; The characteristic function of the covariance operator; Let be a set of independent and identically distributed standard normal random variables.

[0012] Furthermore, the dynamic thermodynamic coupling near-field dynamics theoretical model includes a nonlocal integral THM coupling control equation set, as well as a rate-thermal coupling dynamic damage criterion and a field parameter evolution criterion. A dynamic hydrodynamic coupling peri-field dynamics theoretical model is constructed, specifically including: A nonlocal integral type THM coupled control equation set based on conventional near-field dynamics is established, which includes mechanical field control equations, flow field control equations and thermal field control equations. A rate-thermal coupling dynamic damage criterion based on critical bond elongation is established, wherein the rate-thermal coupling dynamic damage criterion includes the critical fracture elongation function. The critical fracture elongation function Introducing damage function Determine the bond breakage state; where, Let t be the temperature variable and t be the time variable; Based on damage function Define the local damage field of the initial material point cloud With equivalent crack aperture This establishes the permeability of the initial material point cloud. With damage, it evolves exponentially, and the thermal conductivity... Evolution criteria of field parameters that decay linearly with damage.

[0013] Furthermore, the nonlocal integral THM coupled control equation set includes mechanical field control equations, flow field control equations, and thermal field control equations, specifically including: The expression for the governing equations of the mechanical field is as follows:

[0014] in, A unit volume of substance point x Inertial force; x For matter points, for x The ethnic point; For matter points The near field; To be relative to position The force state vector is a variable; For matter points x and The relative positions between them; For matter points x right The applied near-field force state vector; For matter points x right The applied reaction force; To be related to matter points The associated representative volume element; t is the integral variable representing the infinitesimal volume element; t is the time variable. The density of external forces in the volume; The expression for the flow field control equation is:

[0015] For matter points x The rate of change of pore water pressure at a given location over time; The fluid diffusion coefficient, The mechanical coupling coefficient is... Thermodynamic coupling coefficient; The weight function; Pore ​​water pressure; For matter points x right The pore water pressure difference between them; For matter points Mechanical volumetric strain rate at time t For matter points The rate of temperature change at time t; The expression for the thermal field control equation is:

[0016] in, The rate of change of internal energy per unit volume; Thermal conductivity; For matter points The temperature value at time t; Matter Point At any moment t Temperature value; This is a mechanically coupled heat source term; This is a heat source term for seepage and convection.

[0017] Furthermore, step S3 specifically includes: S31. An explicit numerical solver is constructed using a multi-time-step coupling strategy and a domain decomposition parallel algorithm. S32. Using an explicit numerical solver, the parameters of the mechanical field of the initial material point cloud are determined with a mechanical field time step constrained by CFL conditions. The parameters of the thermal field of the initial material point cloud are processed by using a temperature field time step. The parameters of the seepage field of the initial material point cloud are processed by using a seepage field time step. Processing, among which , ,and The computational domain of the initial material point cloud is obtained by using an integer as the integer and ensuring coupling consistency through nested sub-loops and lightweight Picard iterations. S33. Divide the computational domain into multiple non-overlapping subdomains using an explicit numerical solver; and construct a thickness of [missing information] for each subdomain. The ghost layer is used to handle nonlocal interactions across subdomains; before the field variables are updated at each time step, all state information of matter points in the ghost layer is synchronized through the asynchronous MPI communication message passing interface to complete the state evaluation of all bonds. S34. Introduce a dynamic load rebalancing mechanism based on recursive coordinate bisection to redivide the point set of each subdomain computational load in real time to obtain a discrete damage point cloud, wherein the discrete damage point cloud includes dynamic data of field variables of multi-physics fields.

[0018] Further, step S4 specifically includes: S41. The discrete damage point cloud is reconstructed into a continuous crack mesh using the triangulation method; the expression for the continuous crack mesh is: ,in It is a triangular unit. k The triangular elements are indexed; and the geometric parameters of the continuous crack mesh are calculated; the continuous crack network includes a two-dimensional continuous crack network and a three-dimensional continuous crack network, the three-dimensional continuous crack mesh is constructed based on the triangular elements of the two-dimensional continuous crack mesh; wherein, the geometric parameters of the continuous crack mesh include at least: the area of ​​the triangular elements of the two-dimensional continuous crack network and the crack volume of the closed surface elements of the three-dimensional continuous crack network; The expression for the area of ​​the triangular unit is:

[0019] in, Let be the area of ​​the triangular unit. , , These are the three vertices of a triangular unit, and , , All are material points; The expression for the crack volume of the closed surface element is:

[0020] in, The crack volume is the volume of the closed surface element. A set of triangular elements for closed surface elements; This is the normal vector of the triangular element; This represents the centroid displacement of the triangular element; S42. The continuous crack mesh is processed by weight quantization, and the parameters obtained by the simulation of the dynamic hydrodynamic coupling near-field dynamic theory model are mapped to the triangular elements of the continuous crack mesh by linear interpolation to obtain the mechanical strain energy, thermal strain energy and fluid potential energy of the triangular elements. Based on the mechanical strain energy, thermal strain energy and fluid potential energy of the triangular elements, the element energy weight and scale transition factor are obtained. S43. Define an energy release rate criterion based on the unit energy weight and the scale transition factor, and obtain crack propagation characteristics based on the energy release rate criterion.

[0021] Furthermore, the expression for the unit energy weight is:

[0022] in, For unit energy weights, , , These are the coupling coefficients for the mechanical field, thermal field, and flow field, respectively. , , These represent the mechanical strain energy, thermal strain energy, and fluid potential energy of the triangular element, respectively. For unit area, The characteristic length of the triangular unit; The energy release rate criterion is as follows: ,in, Energy release rate, The preset energy release rate threshold; Wherein, the energy release rate The expression is:

[0023] It is a measure of unit energy. , The scale transition factor; For time steps, This represents the change in the crack area of ​​the triangular unit.

[0024] Furthermore, the simulation results include at least the spatiotemporal evolution curves of multiple physical quantities, crack morphology diagrams, crack mesh topology parameters, and energy evolution reports; The spatiotemporal evolution curves of the multiple physical quantities are curves formed by iteratively updating the time step and fitting the dynamic evolution data of the physical field variables. The crack morphology map is generated by extracting material points that have been damaged at a specified time step and reconstructing them using a triangulation method. The crack mesh topology parameters are obtained by updating the time step successively to acquire dynamic statistical data on crack propagation characteristics. The energy evolution report is generated by summarizing the mechanical strain energy, thermal strain energy, fluid potential energy, and energy release rate of the triangular unit at each time step.

[0025] On the other hand, the present invention provides a near-field dynamic simulation system for the hydrothermal-hydraulic coupled dynamic fracture of surrounding rock in underground caverns, used to execute the above-mentioned near-field dynamic simulation method for the hydrothermal-hydraulic coupled dynamic fracture of surrounding rock in underground caverns, the system comprising: The data input module is used to generate a physical parameter field based on the geometric parameters, material physical parameters, initial temperature and pressure stress conditions, impact load parameters, and boundary conditions of the surrounding rock of the underground cavern. The parameter field generation module has a built-in RSCFS algorithm and is used to generate a non-homogeneous parameter field based on the physical parameter field. The simulation calculation module has a built-in dynamic hydrodynamic coupling near-field dynamics theoretical model and an explicit numerical solver, which is used to solve the parameters of the heterogeneous parameter field of the parameter field generation module to obtain discrete damage point cloud and dynamic data of field variables of multiphysics field of discrete damage point cloud. The crack analysis module has a built-in triangulation framework. It uses the triangulation method to reconstruct the discrete damage point cloud into a continuous crack mesh and realizes quantitative analysis of crack propagation characteristics and energy evolution. The results output module is used to receive and integrate the data from the simulation calculation module and the crack analysis module, and visualize the simulation results of the dynamic fracture of the surrounding rock coupled with hydrothermal and hydraulic forces. The simulation results include one or more of the following: spatiotemporal evolution curves of multiple physical quantities, crack morphology diagrams, crack mesh topology parameters, and energy evolution reports. Both the simulation calculation module and the crack analysis module achieve dynamic data updates through iterative updates of the time step.

[0026] Compared with the prior art, the beneficial effects of the present invention are as follows: The dynamic hydrothermal coupling near-field dynamics theoretical model constructed in this invention is the first to combine conventional state near-field dynamics with the full coupling depth of hydrothermal forces. Through the three-field decomposition of force state vector, rate-thermal coupling damage criterion and closed-loop evolution law of field parameters, it accurately captures the THM full coupling effect of surrounding rock under high strain rate impact, solves the technical defects of traditional PD method that can only handle single mechanical fields or simple couplings, and has high prediction accuracy. The RSCFS algorithm proposed in this invention targets the cylindrical / radial symmetry characteristics of the surrounding rock of underground caverns. Through three-layer constraint shuffling and iterative correction of physical and gradient constraints, it ensures the heterogeneity of the parameter field, maintains geometric symmetry and physical spatial correlation, solves the problem of parameter representation distortion in traditional randomization methods, and improves the physical realism of the simulation. The explicit numerical solver developed in this invention, which combines multiple time steps and domain decomposition parallelism, matches the characteristic time scales of different physical fields through a multi-scale time integration strategy. Combined with domain decomposition parallelism, dynamic load rebalancing, and computation-communication overlap technology, it significantly improves computational efficiency and solves the problem of high computational cost in THM coupled simulation. It is suitable for three-dimensional simulation of surrounding rock in large-scale deep underground caverns. The triangulation framework constructed in this invention realizes the geometric reconstruction from discrete damage data to a continuous crack network. Through the energy weight quantization module, multi-field energy data is mapped to crack elements. Based on the energy release rate criterion, quantitative analysis of crack propagation direction, velocity, stability and topological evolution is realized. This solves the technical problem that traditional near-field dynamic simulation cannot quantitatively study crack networks and provides a quantitative tool for the analysis of surrounding rock fracture mechanism. The simulation system of this invention forms an integrated simulation system from parameter characterization, fully coupled dynamic simulation to quantitative analysis of crack network energy. It can accurately output the multi-field evolution law, crack characteristics and energy evolution data of THM coupled dynamic fracture of surrounding rock, providing scientific simulation tools and data support for stability assessment and disaster prevention of surrounding rock in deep underground caverns. It has broad engineering application prospects in the fields of energy extraction, underground storage and geological disaster prevention. Attached Figure Description

[0027] Figure 1 This is a flowchart of a near-field dynamic simulation method for dynamic fracture of surrounding rock in underground caverns using hydrothermal coupling, according to an embodiment of the present invention. Detailed Implementation

[0028] The present invention will be further described in detail below with reference to specific embodiments. The scope of protection of the present invention is not limited to the following embodiments. All equivalent transformations made based on the content of the present invention are within the scope of protection of the present invention.

[0029] This embodiment provides a near-field dynamic simulation method and system for dynamic fracturing of surrounding rock in underground caverns caused by hydrothermal coupling. A deep green sandstone underground cavern is used as the surrounding rock sample. The surrounding rock has a cylindrical symmetrical structure with a maximum diameter of 50 mm and a height of 25 mm, and is subjected to a 350s test. - ¹ High strain rate impact load, initial temperature 25-80℃, confining pressure / permeability 0-15MPa. In the following embodiments, the maximum diameter and height of the surrounding rock sample are constants, and all material points generated discretely have their own stored spatial location, thermal field, pressure field, damage state, and other multi-field physical parameters, which are dynamic variables updated in real time with the time step.

[0030] This invention provides a near-field dynamic simulation method for the dynamic fracture of surrounding rock in underground caverns caused by hydrothermal coupling, such as... Figure 1 As shown, it includes the following steps: S1. The RSCFS algorithm (Radially Symmetric Correlated Field Shuffling) is used to generate controllable randomization of the acquired surrounding rock parameter information to obtain a heterogeneous parameter field that satisfies geometric and physical constraints. The geometric and physical constraints include cylindrical geometric symmetry, exponential spatial correlation, and physical limits and gradient smoothing constraints. In S1, the RSCFS radial symmetric correlation field shuffling algorithm is used to perform controllable randomization processing on the acquired surrounding rock parameter information, specifically including the following steps: S11. Obtain surrounding rock parameter information, perform spatial parameterization on the surrounding rock parameter information to obtain physical parameter field, perform exponential covariance function and Karhunen-Loève expansion (KL expansion) on the physical parameter field to generate a random field of initial parameters of surrounding rock with cylindrical geometric symmetry. In this embodiment, the surrounding rock parameter information first considers the material domain in cylindrical coordinates. The physical parameter field is defined according to the material domain, and the surrounding rock parameter information is updated as the crack evolves. The expression is:

[0031] in, The radial coordinates of the material point. Let be the angular coordinates of the material point. Let be the axial coordinates of the material point. The radius of the surrounding rock sample is denoted as . The height of the surrounding rock sample; The set of positions of matter points; The physical parameter field is defined as a mapping function, and the physical parameter field is expressed as: ,in The dimension of the physical parameter field represents the total number of physical parameters carried by each material point, and represents each position in the material domain. Each corresponds to an n-dimensional physical parameter vector, including physical parameters such as spatial location, thermal field, pressure field, and damage state. Indicates the position of a substance point The physical parameter vector at a given location, and to ensure physical consistency along the circumference, the physical parameter field should satisfy:

[0032] To characterize the spatial heterogeneity of the surrounding rock physical parameters, the physical parameter field is constructed into a random field of initial parameters of the surrounding rock based on the exponential covariance function and the Karhunen-Loève expansion.

[0033] The expression for the exponential covariance function is:

[0034] in, , Let be the position vector of two material points within the physical parameter field; The radial coordinate difference between two material points. This represents the difference in the axial coordinates of two material points; , These are the radial length and axial length of the physical parameter field, respectively. Let V be the variance of the physical parameter field.

[0035] The Karhunen-Loève expanded expression is:

[0036] in, The physical parameters of the random field are reconstructed by Karhunen-Loève expansion; is the mean function of the physical parameter field; M is the number of truncation terms selected according to the calculation accuracy requirements; The eigenvalues ​​are the covariance operators of the physical parameter field; The characteristic function of the covariance operator; Let be a set of independent and identically distributed standard normal random variables.

[0037] S12. A three-layer constrained shuffling operation is performed on the random field of initial parameters of the surrounding rock. The three-layer constrained shuffling operation adopts a hierarchical randomization strategy, applying geometric and physical constraints sequentially in the angular, radial, and axial dimensions to obtain the processed random field of initial parameters of the surrounding rock. The shuffling operation further increases the exponential spatial correlation of the random field of initial parameters of the surrounding rock: the angular dimension undergoes sector rotation to maintain cylindrical symmetry, the radial dimension uses a distance-constrained Fisher-Yates variant to maintain parameter gradient characteristics, and the axial dimension allows for medium-distance exchanges to simulate inter-layer variability. Specifically, this includes: material domain Discretize into a regular grid:

[0038] Where the radial step size Angle step size Axial step size , , , These represent the number of discrete points in the radial, angular, and axial directions, respectively.

[0039] In the angular dimension, to maintain cylindrical symmetry, sector rotation is used instead of element-level exchange; a sector rotation operator is defined. for Then we have:

[0040] in For randomly selected rotation offsets, from the set {0,1,…, Evenly extract from}; ,j, These are indices for the radial, circumferential, and axial layers, respectively. This operation maintains continuity in the circumferential direction while introducing randomness.

[0041] In the radial dimension, a Fisher-Yates variant algorithm with distance constraints is adopted to maintain the radially gradual trend of the parameters. Its radial distance constraint condition is as follows:

[0042] in, The index distance between the two radial layers; These are indices for two different radial layers, corresponding to mesh layers at different radius positions; This represents the maximum permissible radial exchange distance. This ensures that parameter exchange can only occur between adjacent or nearby radial layers.

[0043] To reflect the distance attenuation effect, the formula for calculating the radial parameter exchange probability is:

[0044] in radial layer and The probability of exchange between them; Radial attenuation coefficient, For radial basis exchange probability.

[0045] A wide range of parameter exchanges are allowed in the axial dimension to simulate interlayer variability in real materials:

[0046] in For indexes of two different axial layers The index distance between two axial layers; This represents the maximum permissible axial exchange distance. =4 or 5, allowing for axial interlayer exchange at medium distances.

[0047] The formula for calculating the probability of axial parameter exchange is:

[0048] in, Axial layer and The probability of exchange between them; For the axial base exchange probability, The axial attenuation length, This is an axial interlayer weighting term, which can adjust the exchange tendency between different layers according to the actual physical properties.

[0049] S13. Apply physical limit constraints and gradient smoothing constraints to the random field of initial parameters of the surrounding rock after step 12. Obtain a controllable randomized heterogeneous parameter field through iterative correction. The heterogeneous parameter field includes mechanical field, thermal field and flow field. After completing the shuffling operation in step 12, it is necessary to verify and ensure that the physical parameter field satisfies the physical limit constraints:

[0050] in, This represents the minimum physically feasible range of the physical parameter field. This represents the maximum physically feasible range of the physical parameter field.

[0051] To maintain the smoothness of the physical parameter field, a gradient smoothing constraint is applied to the physical parameter field:

[0052] , These are the absolute values ​​of the gradients of the physical parameter field in the radial and axial directions, respectively. , These are the maximum allowable gradient values ​​for the radial and axial directions, respectively. This constraint ensures that the parameter changes of the enhanced parametric field are continuous and smooth, preserving the heterogeneity while avoiding numerical abrupt changes.

[0053] Further discretization of the physical parameter field after gradient smoothing constraint yields the discrete form of gradient smoothing constraint, which can be expressed as:

[0054]

[0055] The physical parameter field is iteratively modified multiple times to simultaneously satisfy the physical limit constraints of the surrounding rock material and the gradient smoothing constraints of the discrete form, thus obtaining a controllable randomized heterogeneous parameter field.

[0056] In this embodiment, the processing of the initial parameters of the surrounding rock using a random field involves setting the initial and boundary conditions of the surrounding rock in the underground cavern. These boundary conditions include mechanical boundary conditions, thermal boundary conditions, and seepage boundary conditions, in order to construct a hydrodynamic coupling simulation environment that conforms to physical reality. (1) Mechanical boundary conditions: Simulate the in-situ stress field and dynamic disturbances of the surrounding rock. Apply a depth-dependent gradient stress field, including vertical principal stresses, to the outer surface of the model. and horizontal principal stress The stress state after excavation and unloading is applied to the inner wall of the tunnel; viscous or viscoelastic artificial boundaries are set at the bottom and around the model to absorb stress waves generated under dynamic impact loads and reduce false reflections; on the impact load surface of the surrounding rock of the tunnel, a time-varying pressure boundary is applied according to the measured shock wave waveform. Such as incident wave, reflected wave and transmitted wave signals measured by the Split Hopkinson Bar (SHPB) test.

[0057] (2) Thermal boundary conditions: Simulating the heat exchange process under deep high geothermal environment and dynamic loads. To simulate the initial geothermal field, a uniform or linearly distributed thermal field with depth is applied throughout the entire model domain. At the outer boundary of the model, based on temperature monitoring data from the engineering site... Apply an isothermal boundary (Dirichlet boundary). To simulate short-term, transient loading characteristics during dynamic impact, adiabatic boundary conditions (Neumann boundaries) are applied to the sides and loading end faces of the model. ( (This is the outer normal vector) to prevent non-physical heat dissipation.

[0058] (3) Seepage boundary conditions: Simulating the occurrence and migration of groundwater. To simulate the initial pore water pressure field, a hydrostatic pressure that increases linearly with depth is applied throughout the entire model domain. On the inner walls of the cavern and the outer boundary of the model, based on the groundwater recharge and discharge conditions... Apply constant pore pressure boundary (Dirichlet boundary). During dynamic impact, to simulate the rapid loading effect on the loading end face, it is set as a zero-flow boundary (Neumann boundary). .

[0059] S2. Construct a dynamic hydrothermal coupled peri-field dynamics theoretical model and input the heterogeneous parameter field into the dynamic hydrothermal coupled peri-field dynamics theoretical model (THM-PD). The dynamic THM-PD model obtains the initial material point cloud characterizing the surrounding rock through numerical discretization. The initial material point cloud includes the initial state of material points and bonds, specifically including: S21. The dynamic THM-PD model includes a nonlocal integral THM (thermal hydraulic) coupled control equation set, as well as a rate-thermal coupled dynamic damage criterion and field parameter evolution criterion. The specific content of the dynamic THM-PD model construction includes: S211. Establish a nonlocal integral type THM coupled control equation set based on conventional near-field dynamics. The nonlocal integral type THM coupled control equation set includes mechanical field control equation, flow field control equation and thermal field control equation to characterize the coupling relationship between the three fields of mechanical field, thermal field and flow field. The specific steps include: The force state vector state of OSB-PD Decomposed into mechanical elastic force state Force state caused by pore fluid pressure and the force state caused by thermal stress The nonlocal integral expressions for each force state are derived, and the mechanical field control equations for the coupling of hot water and force are established. Under THM coupling, the expression for the force state vector state is:

[0060] Wherein, the mechanical elastic force state Force state caused by pore fluid pressure Force state caused by thermal stress The expressions are as follows:

[0061]

[0062]

[0063] in, Bulk modulus For volumetric strain, For weighted volume, For the influence function, For parameters related to shear modulus, For strain tensor, For key vectors, It is a relative displacement vector. For Biot coefficient, Pore ​​pressure, The fluid-structure interaction coefficient is... The coefficient of thermal expansion is... Relative to reference temperature Temperature changes, This is the thermo-coupling coefficient.

[0064] mechanical elastic force state Force state caused by pore fluid pressure Force state caused by thermal stress Substituting the expression for the force state vector state into the expression for the force state vector state, and then substituting the force state vector state into the basic equations of motion, we obtain the control equations of the mechanical field, where the expression for the control equations of the mechanical field is:

[0065] in, A unit volume of substance point x Inertial force; x For matter points, for x The ethnic point; For matter points The near field; To be relative to position The force state vector is a variable; For matter points x and The relative positions between them, that is x and The key vector between them; For matter points x right The applied near-field force state vector; For matter points x right The applied reaction force; To be related to matter points The associated representative volume element; t is the integral variable representing the infinitesimal volume element; t is the time variable. The density of external forces in the volume; The flow field governing equations are established based on the law of conservation of fluid mass. The expression of the flow field governing equations is as follows:

[0066] For matter points x The rate of change of pore water pressure at a given location over time; For matter points Mechanical volumetric strain rate at time t For matter points The rate of temperature change at time t; It is the osmotic pressure; For matter points x right The pore water pressure difference between them; The weight function; The fluid diffusion coefficient, The mechanical coupling coefficient is... Let be the thermo-coupling coefficient, where Dynamic viscosity; Here, M is the permeability tensor; and Biot modulus. For Biot coefficients; It is porosity; The coefficient of thermal expansion of the fluid volume; The coefficient of thermal expansion of the surrounding rock skeleton is given. This refers to the water pressure storage capacity of porous media.

[0067] The flow field control equations are established based on the classical thermal energy balance equations. The expression for the flow field control equations is as follows:

[0068] in, The rate of change of internal energy per unit volume. For matter points The temperature value at time t; Matter Point At any moment t Temperature value; For mechanically coupled heat source terms, It is the Taylor-Quinni coefficient, which represents the proportion of inelastic work converted into heat; For convective heat transfer energy, It is the fluid density. For specific heat capacity, It is the fluid velocity divergence; is the thermal conductivity.

[0069] S212. Establish a rate-thermal coupling dynamic damage criterion based on critical bond elongation, wherein the rate-thermal coupling dynamic damage criterion includes a critical fracture elongation function. Introducing damage function Determine the bond breakage state, where, Let t be the temperature variable and t be the time variable; Based on damage function Define the local damage field of the initial material point cloud With equivalent crack aperture This establishes the permeability of the initial material point cloud. With damage, it evolves exponentially, and the thermal conductivity... Evolution criteria of field parameters that decay linearly with damage.

[0070] Wherein the critical fracture elongation function The expression is:

[0071] in, This is the quasi-static critical elongation. For rate sensitivity coefficient, For reference strain rate, The coefficient of thermal softening, For reference temperature, The variable is temperature.

[0072] Damage function State variables representing the key, damage function Defined as:

[0073] when When the time is 0, it indicates that the bond is broken and no longer transmits interaction forces.

[0074] When numerically solving the control equations of the mechanical field, a damage function is introduced. The force state vector state is corrected to obtain the effective force state. ,when When the interaction force of the bond disappears, the crack surface naturally initiates and expands.

[0075] based on Define the matter point Local damage field at the location Its near field Weighting of internal fractured bonds, local damage field The formula for calculation is:

[0076] Define equivalent crack aperture Elongation of the broken bond The average value, equivalent crack aperture The formula for calculation is:

[0077] in, For matter points The number of broken bonds in the near-field region. Among them, The elongation of the bond; Crack initiation and propagation cause dynamic changes in permeability and thermal conductivity, forming a THM-coupled closed-loop feedback. The evolution follows the cube law:

[0078] in, Initial penetration rate, The discrete distance between the material points.

[0079] thermal conductivity The evolution follows a linear decay:

[0080] in, The initial thermal conductivity, is the thermal conductivity of the fluid within the fissure.

[0081] S22, the heterogeneous parameter field is input into the dynamic THM-PD model. The dynamic THM-PD model obtains the initial material point cloud characterizing the surrounding rock through numerical discretization. The initial material point cloud includes the initial state of material points and bonds.

[0082] During the crack evolution process of the surrounding rock, the permeability and thermal conductivity are updated in real time through the statistical results of the damage function. The physical parameters are updated by updating the damage function, and the development of damage is further analyzed. Moreover, updating the permeability and thermal conductivity will affect the flow field and the thermal field, realizing the bidirectional coupling between damage and multi-physics field.

[0083] S3. An explicit numerical solver is constructed using a multi-time-step coupling strategy and a domain decomposition parallel algorithm. This explicit numerical solver performs time-domain progressive solving on the initial material point cloud to obtain a discrete damage point cloud. The discrete damage point cloud includes dynamic data of field variables from multiple physics fields. This dynamic data includes the mechanical, thermal, and flow field data of the surrounding rock, as well as discrete damage data determined by the damage function. S31. An explicit numerical solver is constructed using a multi-time-step coupling strategy and a domain decomposition parallel algorithm. S32. Using an explicit numerical solver, the parameters of the mechanical field of the initial material point cloud are determined with a mechanical field time step constrained by CFL conditions. The parameters of the thermal field of the initial material point cloud are processed using a thermal field time step. The parameters of the flow field of the initial material point cloud are processed by using a flow field time step. Processing; among which , ,and The integer represents the heat field time step. Time step, flow field time step All are mechanical field time steps The amplification is an integer multiple; the coupling consistency is ensured by nested sub-loops and lightweight Picard iterations to obtain the computational domain of the initial material point cloud, so as to complete the time-domain progressive solution of the discrete computational domain; time step Inside, from arrive The mechanical field update sequence is as follows:

[0084]

[0085]

[0086]

[0087] in, Indicates displacement. Indicates speed, Indicates acceleration; The density of a substance's point mass. For matter points The near field region, , For matter points volume, It represents the density of body force.

[0088] The mechanical field time step in this embodiment For the minimum time step, =1×10 8 s, warm-up time step =10 Flow field time step =5 By using nested sub-loops and lightweight Picard iterations, convergence is achieved in 2-3 iterations, ensuring coupling consistency.

[0089] S33. Divide the computational domain into multiple non-overlapping subdomains using an explicit numerical solver, and construct a thickness of [missing information] for each subdomain. A ghost layer is used to handle nonlocal interactions across subdomains; before the field variables are updated at each time step, all state information of matter points within the ghost layer is synchronized through an asynchronous MPI communication message passing interface to complete the evaluation of all bonds.

[0090] This embodiment divides the computational domain into four non-overlapping subdomains, constructing a domain with a thickness of... =1.5mm ghost layer.

[0091] S34. A dynamic load rebalancing mechanism based on recursive coordinate bisection is introduced. The point set is re-divided in real time according to the calculated load of each subdomain. In this embodiment, when the load imbalance is >0.2, the point set is re-divided based on the number of bond calculations of material points. At the same time, non-blocking MPI communication is used to achieve overlap between calculation and communication, and the discrete damage point cloud is obtained by solving. The discrete damage point cloud contains dynamic data of field variables of multiphysics, including displacement, temperature, pore pressure of surrounding rock, and discrete damage data determined by the damage function. The discrete damage point cloud is expressed as follows: ; S4. Based on the triangulation method, the discrete damage point cloud is reconstructed into a continuous crack network. Then, through weight quantization and multi-field parameter mapping, the parameters obtained from the dynamic THM-PD model simulation are mapped to the continuous crack network to obtain crack propagation features. S41. The damage point cloud of the surrounding rock is reconstructed into a continuous crack mesh using the triangulation method, whereby... The geometric entity of the crack network is realized by using triangular elements; and the geometric parameters of the continuous crack mesh are calculated; the continuous crack network includes a two-dimensional continuous crack network and a three-dimensional continuous crack network, the three-dimensional continuous crack mesh is constructed based on the triangular elements of the two-dimensional continuous crack mesh; wherein, the geometric parameters of the continuous crack mesh include at least: the area of ​​the triangular elements of the two-dimensional continuous crack network and the crack volume of the closed surface elements of the three-dimensional continuous crack network.

[0092] Continuous crack mesh representation is Each triangular unit Consists of three vertices , , constitute, , , All of them are material points.

[0093] Area of ​​each triangular unit Derived from the cross product of position vectors:

[0094] The expression for the crack volume of a closed surface element is:

[0095] in, For the crack volume of a closed surface element, A set of triangular elements for closed surface elements. Let be the normal vector of the closed surface element. This represents the centroid displacement of the triangular element. For matter points The displacement vector.

[0096] S42. By weight quantization, the parameters calculated by the dynamic THM-PD model are mapped to the triangular elements of the continuous crack mesh using linear interpolation. The mechanical strain energy, thermal strain energy, and fluid potential energy of the triangular elements are obtained. Based on the mechanical strain energy, thermal strain energy, and fluid potential energy of the triangular elements, the element energy weight and scale transition factor are obtained. The scale transition factor bridges the micro and macro scales. Specifically, the dynamic THM-PD model is simulated using a nonlocal integral THM coupled control equation system to obtain the vertices of the triangular elements. mechanical strain energy density Thermal strain energy density Fluid potential energy density Mechanical strain energy density Thermal strain energy density Fluid potential energy density The calculation formulas are as follows:

[0097]

[0098]

[0099] in, for The near field; It is the force state vector; The strain state vector; For the vertex The associated representative volume element; t is the time variable; for The race point, of which the near field region is excluded All points outside the family are family points; The relative positions of vertices and family points. = - ; Represents vertices Influenced by its ethnic points The force; Represents vertices Influenced by its ethnic points The stress; The coefficient of thermal expansion; These are thermoelastic parameters; For material points At any moment Temperature increment; For Biot coefficient, Pore ​​pressure, For volumetric strain;

[0100] And based on the vertex mechanical strain energy density Thermal strain energy density Fluid potential energy density Calculate the mechanical strain energy of the triangular element Thermal strain energy Fluid potential energy The calculation formulas are as follows:

[0101]

[0102]

[0103] in As vertices The associated volume is determined using the Voronoi diagram partitioning method, with vertices as the vertices. The volume of the spatial region at the initial moment. .

[0104] Based on the mechanical strain energy of the triangular element Thermal strain energy Fluid potential energy Computational unit energy weight Unit energy weight The calculation formula is:

[0105] in, , , These are the coupling coefficients of the mechanical field, thermal field, and flow field, respectively. For unit area, The feature length of the unit;

[0106] in, , , It represents the length between any two vertices of a triangular unit.

[0107] And introduce a scale transition factor Used to bridge micro and macro scales, defining the characteristic scale of cracks. Energy-weighted average crack size:

[0108] in, The characteristic scale of the crack, The characteristic length of the material can be taken as the near-field radius. The expression for the crack characteristic scale is:

[0109] S43. An energy release rate criterion is defined based on unit energy weights and scale transition factors, and crack propagation characteristics are obtained based on this criterion, achieving a quantitative characterization of crack network topology evolution. The crack propagation characteristics include at least the crack propagation direction, crack propagation velocity, and crack propagation stability. Specifically, these include: The energy release rate criterion is defined based on the unit energy weight and scale transition factor. The energy release rate criterion is determined by the energy release rate and the energy release rate threshold.

[0110] The energy release rate is defined as the partial derivative of the energy weight with respect to the crack area for a triangular element. Considering crack propagation increment Energy release rate The expression is:

[0111] It is a measure of unit energy. ;in For time steps, This represents the change in the crack area of ​​the triangular unit.

[0112] Energy release rate threshold Determined by the material's fracture toughness, and set as a preset value, the energy release rate criterion is:

[0113] Using energy release rate With the preset energy release rate threshold To determine whether a triangular element participates in the crack propagation process, the crack propagation condition is: If the energy release rate Greater than or equal to the preset energy release rate threshold If the condition is met, the triangular element is determined to be an active crack element and participates in the crack propagation process. The calculation method for crack propagation characteristics is as follows: The crack propagation direction is determined by the spatial gradient of the energy weight. Follow the direction of maximum energy gradient to maximize energy release:

[0114] Actual expansion direction Defined as the displacement direction of the centroid of the crack front in adjacent time steps:

[0115] The two deviation angles The calculation is as follows:

[0116] Crack propagation rate It has a power-law relationship with the energy release rate:

[0117] in, , Rayleigh wave speed, The index is yet to be determined.

[0118] Crack propagation stability is determined by the second derivative of the energy weight, and the energy curvature is defined. for:

[0119] Based on this, a crack propagation stability index is defined. :

[0120] when When the energy release rate decreases with crack propagation, crack propagation becomes stable; when At this time, crack propagation may become unstable.

[0121] Multiscale analysis using the fractal dimension of crack networks Achieving topological representation and fractal dimension The expression is:

[0122] fractal dimension The discrete form can be calculated and dynamically updated through linear regression. and It allows for tracing crack propagation paths and reveals the evolution of crack networks under THM coupling.

[0123] S5. Based on the dynamic data of the field variables of the multiphysics field and the crack propagation characteristics, the spatiotemporal evolution curves of the multiphysics field quantities, crack morphology diagrams, crack mesh topology parameters, and energy evolution reports are generated respectively, obtaining the simulation results of the dynamic fracturing of the surrounding rock of the underground cavern under hydrothermal coupling. Among them: The spatiotemporal evolution curves of multiple physical quantities are curves formed by iteratively updating the time step and fitting the dynamic evolution data of the physical field variables.

[0124] The crack morphology image is generated by extracting material points that have been damaged at a specified time step and reconstructing them using Delaunay triangulation.

[0125] The crack mesh topology parameters are obtained by dynamically calculating the crack area, crack volume, and fractal dimension by updating the time step successively. The Energy Evolution Report is a report generated by statistically analyzing the mechanical strain energy, thermal strain energy, fluid potential energy, and energy release rate of triangular units at each time step.

[0126] Secondly, the present invention provides a near-field dynamic simulation system for the dynamic fracture of surrounding rock in underground caverns caused by hydrothermal coupling, for executing the simulation method described above, including a data input module, a parameter field generation module, a simulation calculation module, a crack analysis module, and a result output module: The data input module is used to input the geometric parameters, material physical parameters, initial temperature and pressure stress conditions, impact load parameters, and boundary conditions of the surrounding rock of the underground cavern. The parameter field generation module has a built-in RSCFS algorithm, which is used to generate a heterogeneous parameter field that conforms to the cylindrical symmetry of the surrounding rock based on the parameters of the data input module. It supports setting radial / axial correlation length, maximum exchange distance, physical limit and gradient constraint parameters. The simulation calculation module incorporates a dynamic THM-PD model and an explicit numerical solver that operates in parallel with multiple time steps and domain decomposition. The explicit numerical solver includes a time integration submodule, a parallel computing submodule, and a load balancing submodule. The time integration submodule executes a multi-scale time integration strategy, the parallel computing submodule executes a parallel domain decomposition algorithm and constructs a ghost layer, and the load balancing submodule monitors and adjusts the computational load in real time. The solution obtains the discrete damage point cloud and the dynamic data of the field variables of the multiphysics field of the discrete damage point cloud, which are used for the spatiotemporal evolution of the damage field. The crack analysis module has a built-in triangulation framework for reconstructing the discrete damage point cloud into a continuous crack network. It includes a mesh reconstruction submodule and an energy analysis submodule. The mesh reconstruction submodule performs triangulation and verifies mesh quality. The energy analysis submodule maps multi-field energy data to triangular elements, calculates energy weights, and quantitatively analyzes crack propagation characteristics for quantitative analysis of energy evolution. The results output module is used to receive and integrate the data from the simulation calculation module and the crack analysis module, and to visualize and output the simulation results of the dynamic fracture of the surrounding rock coupled with hydrothermal and hydraulic forces, including the spatiotemporal evolution curves of multiple physical quantities, crack morphology diagrams, crack network topology characteristics and energy analysis reports. Both the simulation calculation module and the crack analysis module achieve dynamic data updates through iterative updates of the time step.

[0127] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Any simple modifications, equivalent changes and alterations made to the above embodiments based on the technical essence of the present invention shall still fall within the scope of the technical solution of the present invention.

Claims

1. A near-field dynamic simulation method for hydrothermal coupling dynamic fracture of surrounding rock in underground caverns, characterized in that, Includes the following steps: S1. The obtained surrounding rock parameter information is processed by the RSCFS radial symmetric correlation field shuffling algorithm in a controllable randomization manner to obtain a heterogeneous parameter field subject to geometric and physical constraints. In step S1, the RSCFS radial symmetric correlation field shuffling algorithm is used to perform controllable randomization processing on the acquired surrounding rock parameter information, specifically including: S11. Spatial parameterization is performed on the surrounding rock parameter information to obtain a physical parameter field. The physical parameter field is processed by exponential covariance function and Karhunen-Loève expansion to generate a random field of initial parameters of the surrounding rock. S12. Perform a three-layer constraint shuffling operation on the random field of the initial parameters of the surrounding rock. The three-layer constraint shuffling operation adopts a layered randomization strategy, and applies geometric and physical constraints in the angular dimension, radial dimension and axial dimension in sequence to obtain the processed random field of the initial parameters of the surrounding rock. Wherein, the angular dimension is to perform sector rotation to maintain cylindrical symmetry, the radial dimension is to employ a distance-constrained Fisher-Yates variant to maintain parametric gradient characteristics, and the axial dimension allows for medium-distance exchanges to simulate interlayer variability; S13. Apply physical limit constraints and gradient smoothing constraints to the processed initial parameter random field of the surrounding rock, and obtain a controllable randomized heterogeneous parameter field through iterative correction. The heterogeneous parameter field includes a mechanical field, a thermal field, and a flow field. S2. Construct a dynamic hydrothermal coupling near-field dynamics theoretical model, input the heterogeneous parameter field into the dynamic hydrothermal coupling near-field dynamics theoretical model, and obtain the initial material point cloud characterizing the surrounding rock through numerical discretization calculation. The initial material point cloud includes the initial state of material points and bonds. The dynamic thermodynamic coupling near-field dynamics theoretical model includes a nonlocal integral THM coupling control equation set, as well as a thermodynamic coupling dynamic damage criterion and a field parameter evolution criterion. A dynamic hydrodynamic coupling peri-field dynamics theoretical model is constructed, specifically including: A nonlocal integral type THM coupled control equation set based on conventional near-field dynamics is established, which includes mechanical field control equations, flow field control equations and thermal field control equations. A rate-thermal coupling dynamic damage criterion based on critical bond elongation is established. The rate-thermal coupling dynamic damage criterion includes a function, which is introduced into the damage function to determine the bond fracture state; where is a temperature variable and t is a time variable. Based on the damage function definition, the local damage field and equivalent crack aperture of the initial material point cloud are defined, and a field parameter evolution criterion is established in which the permeability of the initial material point cloud evolves exponentially with damage and the thermal conductivity decreases linearly with damage. S3. An explicit numerical solver is constructed using a multi-time-step coupling strategy and a domain decomposition parallel algorithm. The initial material point cloud is solved by the explicit numerical solver in the time domain to obtain a discrete damage point cloud. The discrete damage point cloud includes dynamic data of field variables of a multi-physics field, which includes the mechanical field, thermal field, and flow field data of the surrounding rock, as well as discrete damage data determined by the damage function. S4. Based on the triangulation method, the discrete damage point cloud is reconstructed into a continuous crack mesh, and the parameters obtained by the dynamic hydrodynamic coupling near-field dynamic theory model are mapped to the continuous crack mesh through weight quantization and multi-field parameter mapping to obtain crack propagation characteristics. S5. Based on the dynamic data of the field variables of the multi-physics field and the crack propagation characteristics, generate the spatiotemporal evolution curves of the multi-physics field, crack morphology diagrams, crack mesh topology parameters and energy evolution reports, respectively, and obtain the simulation results of the dynamic fracture of the surrounding rock of the underground cavern by hydrothermal coupling.

2. The near-field dynamic simulation method for hydrothermal coupling dynamic fracture of surrounding rock in underground caverns according to claim 1, characterized in that, The expression for the physical parameter field is: in, For physical parameter fields; For the material domain of the surrounding rock, ,in The radial coordinates of the material point. Let be the angular coordinates of the material point. Let be the axial coordinates of the material point. The radius of the surrounding rock sample is denoted as . The height of the surrounding rock sample; The dimension of the physical parameter field represents the total number of physical parameters carried by each material point; The physical parameter field satisfies , For the position of the matter point The physical parameter vector at that location; The expression for the exponential covariance function is: in, It is an exponential covariance function. , Let them be two material points within the physical parameter field; The radial coordinate difference between two material points. This represents the difference in the axial coordinates of two material points; , These are the radial length and axial length of the physical parameter field, respectively. The variance of the physical parameter field; The Karhunen-Loève expanded expression is: in, The physical parameters of the random field are reconstructed by Karhunen-Loève expansion; is the mean function of the physical parameter field; M is the number of truncation terms selected according to the calculation accuracy requirements; The eigenvalues ​​are the covariance operators of the physical parameter field; The characteristic function of the covariance operator; Let be a set of independent and identically distributed standard normal random variables.

3. The near-field dynamic simulation method for hydrothermal coupling dynamic fracture of surrounding rock in underground caverns as described in claim 1, characterized in that, The nonlocal integral THM coupled control equation set includes mechanical field control equations, flow field control equations, and thermal field control equations, specifically including: The expression for the governing equations of the mechanical field is as follows: in, A unit volume of substance point x Inertial force; x For matter points, for x The ethnic point; For matter points The near field; To be relative to position The force state vector is a variable; For matter points x and The relative positions between them; For matter points x right The applied near-field force state vector; For matter points x right The applied reaction force; To be related to matter points The associated representative volume element; t is the integral variable representing the infinitesimal volume element; t is the time variable. The density of external forces in the volume; The expression for the flow field control equation is: For matter points x The rate of change of pore water pressure at a given location over time; The fluid diffusion coefficient, The mechanical coupling coefficient is... Thermodynamic coupling coefficient; The weight function; Pore ​​water pressure; For matter points x right The pore water pressure difference between them; For matter points Mechanical volumetric strain rate at time t For matter points The rate of temperature change at time t; The expression for the thermal field control equation is: in, for ; Thermal conductivity; For matter points The temperature value at time t; Matter Point At any moment t Temperature value; This is a mechanically coupled heat source term; This is a heat source term for seepage and convection.

4. The near-field dynamic simulation method for hydrothermal coupling dynamic fracture of surrounding rock in underground caverns as described in claim 1, characterized in that, Step S3 specifically includes: S31. An explicit numerical solver is constructed using a multi-time-step coupling strategy and a domain decomposition parallel algorithm. S32. Using an explicit numerical solver, the parameters of the mechanical field of the initial material point cloud are determined with a mechanical field time step constrained by CFL conditions. The parameters of the thermal field of the initial material point cloud are processed by using a temperature field time step. The parameters of the seepage field of the initial material point cloud are processed by using a seepage field time step. Processing, among which , ,and The computational domain of the initial material point cloud is obtained by using an integer as the integer and ensuring coupling consistency through nested sub-loops and lightweight Picard iterations. S33. Divide the computational domain into multiple non-overlapping subdomains using an explicit numerical solver; and construct a thickness of [missing information] for each subdomain. The ghost layer is used to handle nonlocal interactions across subdomains; before the field variables are updated at each time step, the state information of all matter points in the ghost layer is synchronized through asynchronous MPI message passing interface communication to complete the state evaluation of all bonds. S34. Introduce a dynamic load rebalancing mechanism based on recursive coordinate bisection to redivide the point set of each subdomain computational load in real time to obtain a discrete damage point cloud, wherein the discrete damage point cloud includes dynamic data of field variables of multi-physics fields.

5. The near-field dynamic simulation method for hydrothermal coupling dynamic fracture of surrounding rock in underground caverns as described in claim 1, characterized in that, Step S4 specifically includes: S41. The discrete damage point cloud is reconstructed into a continuous crack mesh using the triangulation method; the expression for the continuous crack mesh is: ,in It is a triangular unit. k The triangular elements are indexed; and the geometric parameters of the continuous crack mesh are calculated; the continuous crack network includes a two-dimensional continuous crack network and a three-dimensional continuous crack network, the three-dimensional continuous crack mesh is constructed based on the triangular elements of the two-dimensional continuous crack mesh; wherein, the geometric parameters of the continuous crack mesh include at least: the area of ​​the triangular elements of the two-dimensional continuous crack network and the crack volume of the closed surface elements of the three-dimensional continuous crack network; The expression for the area of ​​the triangular unit is: in, Let be the area of ​​the triangular unit. , , These are the three vertices of a triangular unit, and , , All are material points; The expression for the crack volume of the closed surface element is: in, The crack volume is the volume of the closed surface element. A set of triangular elements for closed surface elements; This is the normal vector of the triangular element; This represents the centroid displacement of the triangular element; S42. The continuous crack mesh is processed by weight quantization, and the parameters obtained by the simulation of the dynamic hydrodynamic coupling near-field dynamic theory model are mapped to the triangular elements of the continuous crack mesh by linear interpolation to obtain the mechanical strain energy, thermal strain energy and fluid potential energy of the triangular elements. Based on the mechanical strain energy, thermal strain energy and fluid potential energy of the triangular elements, the element energy weight and scale transition factor are obtained. S43. Define an energy release rate criterion based on the unit energy weight and the scale transition factor, and obtain crack propagation characteristics based on the energy release rate criterion.

6. The near-field dynamic simulation method for hydrothermal coupling dynamic fracture of surrounding rock in underground caverns according to claim 5, characterized in that, The expression for the unit energy weight is: in, For unit energy weights, , , These are the coupling coefficients for the mechanical field, thermal field, and flow field, respectively. , , These represent the mechanical strain energy, thermal strain energy, and fluid potential energy of the triangular element, respectively. For unit area, The characteristic length of the triangular unit; The energy release rate criterion is as follows: ,in, Energy release rate, The preset energy release rate threshold; Wherein, the energy release rate The expression is: It is a measure of unit energy. , The scale transition factor; For time steps, This represents the change in the crack area of ​​the triangular unit.

7. The near-field dynamic simulation method for hydrothermal coupling dynamic fracture of surrounding rock in underground caverns according to any one of claims 1-6, characterized in that, The simulation results include at least the spatiotemporal evolution curves of multiple physical quantities, crack morphology diagrams, crack mesh topology parameters, and energy evolution reports; The spatiotemporal evolution curves of the multiple physical quantities are curves formed by iteratively updating the time step and fitting the dynamic evolution data of the physical field variables. The crack morphology map is generated by extracting material points that have been damaged at a specified time step and reconstructing them using a triangulation method. The crack mesh topology parameters are obtained by updating the time step successively to acquire dynamic statistical data on crack propagation characteristics. The energy evolution report is generated by summarizing the mechanical strain energy, thermal strain energy, fluid potential energy, and energy release rate of the triangular unit at each time step.

8. A near-field dynamic simulation system for hydrothermal coupling dynamic fracture of surrounding rock in underground caverns, characterized in that, A near-field dynamic simulation method for performing dynamic fracturing of surrounding rock in underground caverns according to any one of claims 1-7, the system comprising: The data input module is used to generate a physical parameter field based on the geometric parameters, material physical parameters, initial temperature and pressure stress conditions, impact load parameters, and boundary conditions of the surrounding rock of the underground cavern. The parameter field generation module has a built-in RSCFS algorithm and is used to generate a non-homogeneous parameter field based on the physical parameter field. The simulation calculation module has a built-in dynamic hydrodynamic coupling near-field dynamics theoretical model and an explicit numerical solver, which is used to solve the parameters of the heterogeneous parameter field of the parameter field generation module to obtain discrete damage point cloud and dynamic data of field variables of multiphysics field of discrete damage point cloud. The crack analysis module has a built-in triangulation framework, which is used to reconstruct the discrete damage point cloud into a continuous crack mesh and realize quantitative analysis of crack propagation characteristics and energy evolution. The results output module is used to receive and integrate the data from the simulation calculation module and the crack analysis module, and to visualize and output the simulation results of the dynamic fracture of the surrounding rock coupled with hydrothermal and hydraulic forces. The simulation results include the spatiotemporal evolution curves of multiple physical quantities, crack morphology diagrams, crack mesh topology parameters, and energy evolution reports. Both the simulation calculation module and the crack analysis module achieve dynamic data updates through iterative updates of the time step.

Images