A method and system for calculating the damage depth of surrounding rock in high-geothermal-tunnel blasting excavation

Abstract

The application provides a high-geothermal-tunnel blasting excavation surrounding rock damage depth calculation method and system, relates to the technical field of tunnel engineering, and obtains tunnel section characteristic parameters and blast hole layout parameters, establishes a three-dimensional grid model, and applies a pre-explosion heat exchange boundary condition to calculate an initial temperature field; extracts a pre-explosion temperature dynamic correction detonation speed of a blast hole unit, and solves a blasting pressure time history curve, applies a detonator delay sequence to a grid boundary integral solution to obtain a blasting damage field; applies a post-explosion ventilation cooling heat exchange boundary to calculate a real-time temperature difference and a thermal strain, calculates a thermal damage variable when exceeding a limit, and nonlinearly superimposes the blasting damage variable to obtain a total damage; divides multiple cross sections to construct a local polar coordinate system, searches a farthest grid unit of the total damage exceeding a limit along an azimuth angle and a radial direction, and calculates a distance from the grid unit to a tunnel profile to output as a surrounding rock damage depth. The application realizes accurate engineering quantification mapping of a composite damage field.

Description

Technical Field

[0001] This invention relates to the field of tunnel engineering technology, specifically to a method and system for calculating the damage depth of surrounding rock during high-temperature tunnel blasting excavation. Background Technology

[0002] With the advancement of deep-buried and long tunnel projects in my country, high ground temperatures have become a core bottleneck restricting safe and efficient construction. In traditional tunnel excavation, drill-and-blast methods remain the mainstream means of surrounding rock fracturing and tunneling. To improve the level of refined control of construction plans, using finite element numerical simulation technology to guide the design of traditional blasting schemes has become an inevitable path for the industry's digital and intelligent upgrade. However, current simulation analysis technology is not sufficiently integrated with traditional high-temperature blasting engineering. Conventional numerical calculation models lack a profound deconstruction of the multi-physics field evolution laws under high-temperature environments, resulting in a significant disconnect between software calculation results and actual on-site mechanical responses, failing to achieve high-quality multi-field coupling analysis.

[0003] Existing technology establishes a constitutive model of rock blasting damage based on the physical and mechanical properties of the surrounding rock, compiles it into a user-defined program, and embeds it into dynamic finite element software. In the model establishment stage, by determining the equivalent elastic ring of the blasting of each section of the tunnel face, the explosive pressure is simplified into a triangular pulse load and applied to the equivalent elastic ring. This simulates the micro-delay blasting of the blasting holes at the tunnel face in the software and outputs a command stream K-file displaying the dynamic analysis. Finally, the solver performs transient dynamic solutions by reading the surrounding rock dynamic parameters in the K-file, defines the spatial distance between the boundary of the surrounding rock damage variable being greater than zero and the excavation outline, and thus determines the surrounding rock damage depth at different locations on the tunnel excavation section.

[0004] Existing numerical simulation methods lead to severe distortions in the final calculation results when applied to the extremely complex conditions of high-temperature tunnels. First, these methods fail to reflect the influence of the geothermal environment on the detonation performance of explosives. In actual high-temperature tunneling environments, the high temperature of the surrounding rock in the blast hole causes thermal degradation of the loaded explosives, resulting in a decrease in their actual detonation velocity and detonation pressure. However, this method still uses constant standard parameters when constructing the triangular dynamic load time function, causing a disconnect between the source input of the dynamic load and the actual engineering situation. Second, this method attributes surrounding rock damage solely to mechanical tensile failure caused by the blasting dynamic load, failing to encompass secondary thermal damage caused by post-blast ventilation and cooling. In real high-temperature excavation faces, forced convection heat transfer by cold air will induce transient temperature differences and thermal strain components in the shallow surface layer of the high-temperature exposed rock wall, thereby triggering thermal tensile cracking failure of the rock. Because it is impossible to comprehensively consider the weakening feedback of the high-temperature environment on the detonation load and the superimposed cracking effect of ventilation and rapid cooling on the surrounding rock, the final defined damage variables and output damage depth results cannot truly and completely reflect the actual spatial distribution of the composite damage field under the high-temperature tunnel excavation face.

[0005] The information disclosed in the background section is only intended to enhance the understanding of the background of this disclosure, and therefore may include information that does not constitute prior art known to those skilled in the art. Summary of the Invention

[0006] The purpose of this invention is to provide a method and system for calculating the depth of damage to surrounding rock during high-temperature tunnel blasting excavation, so as to solve the problems mentioned in the background art.

[0007] To achieve the above objectives, the present invention provides the following technical solution: A method and system for calculating the damage depth of surrounding rock during high-temperature tunnel blasting excavation, comprising the following steps: Step 1: Obtain tunnel cross-sectional characteristic parameters and borehole layout parameters, determine the tunnel inner profile based on the tunnel cross-sectional characteristic parameters, establish a three-dimensional spatial model, perform meshing on the three-dimensional spatial model to obtain a three-dimensional mesh model, apply pre-blast heat transfer boundary conditions to the three-dimensional mesh model, and calculate the temperature field distribution of the three-dimensional mesh model. Step 2: Obtain the explosive density and detonator delay sequence. Based on the borehole layout parameters, determine each borehole unit, extract the pre-detonation temperature of each borehole unit, calculate the detonation velocity, calculate the blasting pressure based on the explosive density and convert it into a blasting pressure time history curve. According to the detonator delay sequence, apply the blasting pressure time history curve in stages to the grid boundary of each corresponding borehole layout position for integration and solution to obtain the blasting damage field. Step 3: Apply a post-blast ventilation cooling heat transfer boundary to the three-dimensional mesh model, calculate the mesh element temperature at each cooling time step, calculate the thermal strain component according to the thermoelastic constitutive equation, calculate the thermal damage variable when the thermal strain component exceeds the limit, and nonlinearly superimpose the thermal damage variable and the blast damage variable of the blast damage field to obtain the total damage of each mesh element. Step 4: Divide the tunnel inner profile into multiple detection sections along its geometric axis. Using the intersection of the tunnel inner profile's geometric axis and the detection section as the origin of the polar coordinate system, perform a radial search along the azimuth angle into the depth of the rock mass. Retrieve the farthest grid cell that satisfies the condition that the total damage is greater than a preset critical damage threshold. Calculate the geometric distance of the farthest grid cell from the tunnel inner profile as the surrounding rock damage depth at that azimuth angle.

[0008] Furthermore, the tunnel cross-sectional characteristic parameters include: designed excavation span, excavation height, arch radius, sidewall height, tunnel axis extension direction, and length; The parameters for the borehole layout include: borehole diameter, borehole location, and borehole outward slope angle. Taking the intersection of the tunnel floor centerline and the tunnel face as the origin, the direction of the designed excavation span is defined as... The axis, defining the direction of the excavation height as The direction of the tunnel axis extension is defined as... Establish a global coordinate system based on the axis, and determine the inner contour of the tunnel based on the excavation span, excavation height, sidewall height, and arch radius; Using the center of the tunnel's inner profile as a reference, along the... The positive and negative directions of the axis and the aforementioned The axis extends symmetrically outward in both positive and negative directions by five times the radius of the arch, serving as the boundary dimension, and extends along the axis... Stretch along the axial direction to the stated length to obtain a cubic solid, while simultaneously stretching along the stated length... The tunnel inner contour is stretched along the axial direction according to the length to obtain the tunnel excavation entity. Boolean subtraction is performed between the cube entity and the tunnel excavation entity to obtain a three-dimensional space model. In the area extending outward from each of the aforementioned borehole locations to a predetermined multiple of the borehole diameter, and in the area where the tunnel inner contour extends a predetermined distance into the rock mass, the grid is refined, and the size of the grid unit in this area is controlled within a first predetermined size range. A gradient mesh is performed from the outer boundary of the mesh refinement region to the outermost boundary of the cube entity. The size of the mesh cells in this region is adaptively enlarged as the distance from the outer boundary of the mesh refinement region increases, and is enlarged to a second preset size range at the outermost boundary, wherein the lower limit of the second preset size range is greater than the upper limit of the first preset size range. A common-node topological connection is made at the interface between the refined mesh region and the gradually subdivided mesh region, and the elastic modulus, Poisson's ratio and density of the rock are assigned as mechanical properties to the whole field mesh elements to obtain a three-dimensional mesh model.

[0009] Furthermore, data on ground temperature, ventilation duration, ventilation speed, and airflow temperature are obtained. The ground temperature data is applied to the outermost boundary to bring all grid cells to an initial thermal equilibrium state. The equivalent hydraulic diameter of the tunnel cross-section is calculated based on the tunnel cross-sectional characteristic parameters. The kinematic viscosity, thermal diffusivity, and thermal conductivity of the air are determined based on the airflow temperature. The Reynolds number of the airflow inside the tunnel is obtained by dividing the product of the ventilation velocity and the equivalent hydraulic diameter by the kinematic viscosity. The Prandtl number is obtained by dividing the kinematic viscosity by the thermal diffusivity. The Nusselt number of the tunnel wall is obtained by substituting the Reynolds number and the Prandtl number into the turbulent convection heat transfer correlation. The convective heat transfer coefficient of the tunnel wall is calculated by dividing the product of the Nusselt number and the thermal conductivity by the equivalent hydraulic diameter. The convective heat transfer coefficient and the airflow temperature are applied as convective heat transfer boundary conditions to the grid surface at the location of the tunnel inner profile. Using the ventilation duration as the total solution time, the finite element software is called to calculate the nodal temperature values ​​of each grid element under the convective heat transfer boundary conditions, thus obtaining the temperature field distribution of the three-dimensional grid model.

[0010] Furthermore, based on the location of the borehole and the diameter of the detonator borehole, the grid cells falling within the geometric boundary range of each borehole are divided into corresponding borehole cells. The average temperature of all grid nodes in each borehole unit is extracted as the pre-blast temperature of each borehole. The initial detonation velocity and critical deterioration temperature threshold of the explosive are obtained, and it is determined whether the pre-detonation temperature is greater than the critical deterioration temperature threshold. When the pre-detonation temperature is less than the critical deterioration temperature threshold, the initial detonation velocity is directly used as the detonation velocity. When the pre-detonation temperature is greater than or equal to the critical deterioration temperature threshold, the temperature difference between the pre-detonation temperature and the critical deterioration temperature threshold is calculated, and the difference between the initial detonation velocity and the product of the temperature difference and the preset temperature sensitivity coefficient is taken as the reduced detonation velocity. The corrected blast pressure is obtained by multiplying the explosive density by the square of the detonation velocity.

[0011] Furthermore, the modified blasting pressure is combined with the corresponding explosive detonation time to construct the blasting pressure time history curve; According to the detonator delay sequence, the detonation time node is set for each of the borehole units; At the corresponding detonation time node, the blast pressure time history curve is used as a transient dynamic load and applied to the grid boundary of the corresponding borehole unit; Solve the propagation process of the dynamic stress wave caused by the transient dynamic load in the grid element to obtain the principal stress of each grid element; When the principal stress exceeds the preset rock strength threshold, damage evolution is triggered, the blast damage variables of each grid cell are calculated, and the blast damage field is obtained.

[0012] Furthermore, the convective heat transfer coefficient and the airflow temperature are applied to the grid surface at the location of the tunnel profile to calculate the real-time temperature of each grid cell at the current cooling step. Subtracting the real-time temperature from the pre-explosion temperature of each grid cell yields the transient cooling temperature difference of each grid cell at the current cooling time step. Based on the thermoelastic constitutive equation, the thermal strain components of each grid cell in the corresponding spatial axis direction are calculated by multiplying the transient cooling temperature difference by the preset linear expansion coefficient of the rock. The thermal strain component is compared with a preset rock ultimate tensile strain threshold for determination; When the thermal strain component is greater than the ultimate tensile strain threshold, the thermal damage variable of each grid element is calculated by nonlinear accumulation based on the excess strain of the thermal strain component exceeding the ultimate tensile strain threshold.

[0013] Furthermore, for each of the grid cells in the entire field, the corresponding blast damage variable and thermal damage variable are extracted respectively; Subtract the blast damage variable and the thermal damage variable from the constant one respectively to obtain the first integrity coefficient characterizing the proportion of grid cells that are not damaged after blasting, and the second integrity coefficient characterizing the proportion of grid cells that are not damaged after rapid cooling. Multiply the first complete coefficient with the second complete coefficient to obtain the total complete coefficient; The total damage of each mesh cell is obtained by subtracting the total integrity coefficient from the constant.

[0014] Furthermore, a geometric axis perpendicular to the inner contour of the tunnel is used to divide the tunnel into multiple detection sections. On each detection section, a polar coordinate system is established with the intersection of the geometric axis and the corresponding detection section as the origin of the polar coordinate system. The azimuth angle is cyclically changed with a preset angle step size, and a search ray is emitted from the origin of the polar coordinate system to the deep part of the rock mass along the current azimuth angle; Along the search ray, all grid cells traversed by the search ray are traversed, and the total damage corresponding to each grid cell is extracted; The total damage is compared with a preset critical damage threshold to filter out damaged mesh cells whose total damage is greater than the critical damage threshold. Among the screened damaged mesh cells, the target mesh cell that is farthest from the origin of the polar coordinate system is identified; Obtain the center coordinates of the target grid cell, and obtain the boundary coordinates of the intersection point of the tunnel inner profile and the search ray at the current azimuth angle; Calculate the spatial geometric distance between the center coordinates and the boundary coordinates, and output the spatial geometric distance as the surrounding rock damage depth at the current azimuth angle on the current detection section.

[0015] A system for calculating the depth of damage to surrounding rock during high-temperature tunnel blasting excavation, comprising: The model building module is used to obtain tunnel cross-sectional characteristic parameters and borehole layout parameters, determine the tunnel inner profile based on the tunnel cross-sectional characteristic parameters, build a three-dimensional spatial model, perform meshing on the three-dimensional spatial model to obtain a three-dimensional mesh model, apply pre-blast heat transfer boundary conditions to the three-dimensional mesh model, and calculate the temperature field distribution of the three-dimensional mesh model. The blasting damage module is used to acquire explosive density and detonator delay sequence, determine each borehole unit based on the borehole layout parameters, extract the pre-blast temperature of each borehole unit, calculate the detonation velocity, calculate the blasting pressure based on the explosive density and convert it into a blasting pressure time history curve, and apply the blasting pressure time history curve in stages to the grid boundary of each corresponding borehole layout position according to the detonator delay sequence for integration solution to obtain the blasting damage field; The heat loss superposition module is used to apply the post-explosion ventilation cooling heat transfer boundary to the three-dimensional mesh model, calculate the mesh unit temperature at each cooling time step, calculate the thermal strain component according to the thermoelastic constitutive equation, calculate the thermal damage variable when the thermal strain component exceeds the limit, and nonlinearly superimpose the thermal damage variable and the blast damage variable of the blast damage field to obtain the total damage of each mesh unit. The depth retrieval module is used to divide the tunnel inner profile into multiple detection sections along the geometric axis. Taking the intersection of the geometric axis of the tunnel inner profile and the detection section as the origin of the polar coordinate system, it performs a radial search along the azimuth angle into the deep rock mass, retrieves the farthest grid cell that satisfies the condition that the total damage is greater than a preset critical damage threshold, and calculates the geometric distance of the farthest grid cell from the tunnel inner profile as the surrounding rock damage depth at that azimuth angle.

[0016] Compared with the prior art, the beneficial effects of the present invention are: This invention addresses the impact of ambient temperature on explosive detonation performance in high geothermal environments. It establishes a three-dimensional mesh model using acquired tunnel cross-sectional characteristic parameters and borehole layout parameters, and calculates the temperature field distribution by applying pre-detonation heat transfer boundary conditions. Based on this, the pre-detonation temperature of each borehole element is extracted to dynamically calculate the corresponding detonation velocity and blasting pressure, thereby constructing a blasting pressure time-history curve reflecting the thermal degradation characteristics of the explosive. This curve is then applied to the mesh boundary at each borehole location for integration to obtain the true blasting damage field. Compared to existing technologies that use constant standard detonation parameters, leading to a disconnect between the dynamic load source input and the actual conditions of high geothermal engineering, and resulting in distorted calculation results, this improvement achieves deep coupling between heat conduction and dynamic calculation, significantly enhancing the realism of initial dynamic damage calculations. This invention also incorporates the secondary thermal cracking mechanism induced by post-blast ventilation and rapid cooling. A post-blast ventilation cooling heat transfer boundary is applied to the three-dimensional mesh model. The temperature of the mesh unit at each cooling step is calculated, and the thermal strain component is calculated based on the thermoelastic constitutive equation. When the thermal strain component exceeds the limit, the thermal damage variable is calculated. The thermal damage variable and the blast damage variable of the blast damage field are nonlinearly superimposed to obtain the total damage of each mesh unit. Subsequently, multiple detection sections are divided along the geometric axis of the tunnel's inner profile. The intersection of the geometric axis of the tunnel's inner profile and the detection section is taken as the origin of the polar coordinate system. A radial search is performed along the azimuth angle towards the depth of the rock mass. The farthest mesh unit that satisfies the condition that the total damage is greater than a preset critical damage threshold is retrieved. The geometric distance of the farthest mesh unit from the tunnel's inner profile is calculated as the surrounding rock damage depth at that azimuth angle. Compared to existing technologies that attribute damage solely to mechanical tensile damage caused by explosive dynamic loads, ignore thermal cracking caused by sudden cooling, and suffer from fuzzy depth extraction mapping, this invention, through the nonlinear superposition of dual damage and multi-section spatial coordinate amplitude retrieval, not only ensures the integrity of the evolution of composite damage but also achieves accurate mapping of three-dimensional spatial damage data to the engineering quantitative level. Attached Figure Description

[0017] Figure 1 This is a schematic diagram of the overall method flow of the present invention; Figure 2 This is a block diagram of the overall system modules of the present invention; Figure 3 This is a comparison diagram of the multi-field coupled composite damage between existing technologies and the present invention. Detailed Implementation

[0018] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments.

[0019] It should be noted that, unless otherwise defined, the technical or scientific terms used in this invention should have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in this invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed following the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0020] Example: Please see Figure 1 and Figure 3 The present invention provides a technical solution: A method for calculating the damage depth of surrounding rock during high-temperature tunnel blasting excavation, comprising the following steps: Step 1: Obtain tunnel cross-sectional characteristic parameters and borehole layout parameters, determine the tunnel inner profile based on the tunnel cross-sectional characteristic parameters, establish a three-dimensional spatial model, perform meshing on the three-dimensional spatial model to obtain a three-dimensional mesh model, apply pre-blast heat transfer boundary conditions to the three-dimensional mesh model, and calculate the temperature field distribution of the three-dimensional mesh model.

[0021] In this embodiment, the tunnel cross-sectional characteristic parameters include: designed excavation span, excavation height, arch radius, sidewall height, tunnel axis extension direction, and length; The parameters for the borehole layout include: borehole diameter, borehole location, and borehole outward slope angle. In the specific engineering application scenario of this embodiment, the specific values ​​of the tunnel cross-section characteristic parameters are: the designed excavation span is 14.2. The excavation height is 11.5 meters. The radius of the arch is 7.1. The side wall height is 4.4 meters. The tunnel axis extends in a unidirectional excavation direction, and the tunnel length is 40.0 meters. The specific values ​​for the borehole layout parameters are as follows: the diameter of the detonator borehole is 45 mm. The borehole depth is 4.2 meters. The spatial relative coordinates of each blast hole in the two-dimensional cross-section of the working face are input according to the layout parameter matrix; the outward inclination angle of the blast holes for the surrounding smooth blasting holes is... .

[0022] Taking the intersection of the tunnel floor centerline and the tunnel face as the origin, the direction of the designed excavation span is defined as... The axis, defining the direction of the excavation height as The direction of the tunnel axis extension is defined as... Establish a global coordinate system based on the axis, and determine the inner contour of the tunnel based on the excavation span, excavation height, sidewall height, and arch radius; In the underlying layer of the 3D geometric modeling software, the origin O(0,0,0) is set at the intersection of the transverse centerline of the tunnel floor and the tunnel face. The transverse axis of symmetry of the designed excavation span is defined as... The axis, with the left side wall of the tunnel facing outwards as negative and the right side wall facing outwards as positive, defines the vertical upward direction of the excavation height as... The direction of the tunnel axis extension, i.e., the longitudinal excavation advance direction, is defined as... The axes are used to establish a global Cartesian coordinate system. Based on the input control segments such as excavation span, excavation height, sidewall height, and arch radius, a global Cartesian coordinate system is established. Curve fitting is performed on a two-dimensional cross-section with a radius of 0 to uniquely determine the tunnel's inner profile. This provides a unique and direction-free spatial mathematical reference for the spatial positioning of all subsequent finite element mesh nodes, the application of surface force directions to transient dynamic impact loads, and the multidimensional geometric retrieval of multi-section polar coordinate rays.

[0023] Using the center of the tunnel's inner profile as a reference, along the... The positive and negative directions of the axis and the aforementioned The axis extends symmetrically outward in both positive and negative directions by five times the radius of the arch, serving as the boundary dimension, and extends along the axis... Stretch along the axial direction to the stated length to obtain a cubic solid, while simultaneously stretching along the stated length... The tunnel inner contour is stretched along the axial direction according to the length to obtain the tunnel excavation entity. Boolean subtraction is performed between the cube entity and the tunnel excavation entity to obtain a three-dimensional space model. Using the center of the cross-section of the tunnel's inner profile as the geometric reference point, along the... The positive and negative directions of the axis and the aforementioned The axis extends symmetrically outward in both positive and negative directions by five times the radius of the arch, resulting in a total span of 71.0 mm for the large surrounding rock sections in the cross-section. The total height extends outwards by 35.5 meters both below the base plate and above the vault. This forms the outer control boundary scale of the plane; subsequently along the... The topology is translated and stretched along the tunnel length along the axial direction to construct a cubic entity that encloses the entire tunnel excavation disturbance area. Simultaneously, the two-dimensional closed inner contour of the tunnel is separately transferred along the... The tunnel is stretched along the axial direction by the stated length to obtain a tunnel excavation entity representing the space removed during excavation. Finally, a Boolean subtraction operation is performed on the cube entity and the tunnel excavation entity in the software background to remove the cavity volume inside the tunnel, generating a final three-dimensional spatial model containing the surrounding rock retaining structure.

[0024] By strictly limiting the boundary dimensions of the cubic entity to more than five times the radius of the arch, in finite element dynamics simulations, when the intense stress wave generated by the explosive explosion propagates to the artificially cut geometric boundary, if the boundary is too close, it will cause a severe artificial boundary reflection effect, resulting in the false wavefront rebounding inward and superimposing with the real surrounding rock strain field, thus compromising the accuracy of damage assessment. Simultaneously, in heat conduction simulations, the convective boundary temperature gradient caused by post-blast ventilation will penetrate deep into the rock mass. By extending the far-field outward by five times the radius, it can be ensured that the dynamic shock wave naturally attenuates into a weak elastic sound wave before reaching the outer boundary, and the original geothermal field at the far-field boundary maintains a constant initial thermal equilibrium state throughout the ventilation duration, thereby eliminating the computational boundary distortion caused by geometric truncation. Performing Boolean subtraction is to accurately deduct the volume of rock that has been excavated and removed before meshing, providing the cavity inner wall as the true mechanical boundary for the loading of the blasting pressure time history curve.

[0025] In the area extending outward from each of the aforementioned borehole locations to a predetermined multiple of the borehole diameter, and in the area where the tunnel inner contour extends a predetermined distance into the rock mass, the grid is refined, and the size of the grid unit in this area is controlled within a first predetermined size range. Furthermore, non-uniform adaptive mesh discretization is performed on the three-dimensional spatial model. Within each of the stated borehole locations, extending outwards to a preset multiple of the borehole diameter—in this embodiment, the preset multiple is specifically set to 3 times the borehole diameter, i.e., extending outwards by 0.135 degrees from the geometric center line of each borehole. Within the cylindrical envelope area, and within the area where the tunnel outline extends a predetermined distance into the deeper rock mass, in this embodiment, the predetermined distance is specifically set to 2.5. That is, extending 2.5 kilometers from the cave wall into the depths of the mountain. High-density mesh refinement is performed within the annular thickness region, and the size of the finite element mesh elements within this refined region is strictly controlled within a first preset size range, specifically 0.03. Up to 0.06 .

[0026] In blasting mechanics, at the moment of explosive detonation, the rock mass inside the borehole is first bombarded by a shock wave of extremely high pressure, forming a pulverized ring tightly attached to the borehole wall. Based on engineering analogy and rock blasting and fragmentation theory, for the weakly weathered hard granite medium selected in this embodiment, the radius of the pulverized ring typically extends outward to 2 to 3 times the initial radius of the borehole. Therefore, setting the preset multiple to 3 times the borehole diameter can precisely encompass the extremely high strain rate core response zone where the explosive detonation transitions from the detonation chemical reaction to irreversible near-field mechanical damage in the solid. If this multiple is set too small, such as less than 2 times, the detonation shock wave will prematurely jump from the fine grid into the outer coarse grid before attenuating into an elastic-plastic stress wave, causing severe numerical dissipation and waveform pseudo-reflection at the grid interface, leading to a large-scale distortion in the determination of the source damage variables.

[0027] In this embodiment, the depth of the blast hole is 4.2 mm. Following micro-delay blasting at the tunnel face, due to the longitudinal and lateral propagation of stress waves from multiple pores, a loosening zone extending deep into the rock mass is formed around the tunnel's inner contour. Based on experience excavating hard rock, the macroscopic development depth of this loosening zone is generally around 1.5 meters. Up to 2.2 Between. During continuous ventilation of 14400.0 During its lifecycle, 22.0 The cold air is effective against 65.0. Forced convection heat transfer occurs through the high-temperature cave walls, with heat being conducted self-generated from the shallow surface to the deeper interior of the bedrock. Based on the unsteady-state heat transfer temperature gradient derivation, the physical penetration depth of the zone caused by convective contraction and thermal stress disturbance and abrupt temperature changes within 4.0 hours is exactly 2.0. Up to 2.5 Within the annular thickness range. Therefore, the preset distance is strictly locked at 2.5. It can perfectly intercept the complete multi-field coupled shear failure zone where the mechanical cracks from blasting and the secondary cooling, contraction, and thermal stretching cracks from ventilation overlap and interact.

[0028] The area around the blast hole is a high-strain-rate, extremely high-pressure response zone at the moment of blasting, while the tunnel wall surface along the tunnel profile is the region with the highest frequency of temperature gradients and tensile thermal stress gradients caused by forced convection heat transfer from cold air ventilation after blasting. If the mesh elements here are too coarse, due to the constraints of the stable time step between the characteristic wave velocity and the characteristic length of the mesh in explicit dynamic calculations, a coarse mesh will not only fail to capture the abrupt jumps in the stress wave front and the large plastic strain gradients, but will also cause serious numerical dissipation and temperature field smoothing errors, directly resulting in the inability to accurately activate and determine blast damage variables and thermal damage variables at the micro-mesh scale.

[0029] A gradient mesh is performed from the outer boundary of the mesh refinement region to the outermost boundary of the cube entity. The size of the mesh cells in this region is adaptively enlarged as the distance from the outer boundary of the mesh refinement region increases, and is enlarged to a second preset size range at the outermost boundary, wherein the lower limit of the second preset size range is greater than the upper limit of the first preset size range. A gradient mesh is applied from the outermost boundary of the refined mesh region towards the outermost isothermal boundary of the cube. The mesh cell size in this region is adaptively scaled up with increasing spatial geometric distance from the outermost boundary of the refined mesh region. At the outermost boundary of the cube, the mesh cell size is scaled up to a second preset size range, specifically 0.6. Up to 1.2 The lower limit of the second preset size range is significantly greater than the upper limit of the first preset size range, satisfying the gradient topology constraint of monotonically increasing mesh size.

[0030] As the blasting dynamic wave propagates deep into the mountain, it is accompanied by geometric diffusion and material damping absorption, at a distance of 2.5 kilometers from the tunnel center. In the far-field surrounding rock, the stress wave energy has dissipated to below the rock strength threshold, and the surrounding rock only undergoes weak linear elastic deformation without any irreversible mechanical damage. Simultaneously, the strong cold airflow from ventilation and cooling cannot cause any substantial convective temperature degradation to the deep rock mass in the far field. Using gradually enlarging coarse mesh elements here can geometrically reduce the node and element matrix size of the entire finite element model while ensuring computational accuracy, shortening the iterative time of the solver during time-step progression, and achieving an optimal balance between computational accuracy and efficiency.

[0031] A common-node topological connection is made at the interface between the refined mesh region and the gradually subdivided mesh region, and the elastic modulus, Poisson's ratio and density of the rock are assigned as mechanical properties to the whole field mesh elements to obtain a three-dimensional mesh model.

[0032] Finally, strict common-node topology connections are performed at the interface between the refined mesh region and the gradient mesh region to eliminate any artificial physical discontinuities caused by misaligned meshes. Next, the basic physical and mechanical properties of real rock obtained from indoor rock mechanics tests are uniformly assigned as fundamental mechanical properties to all meshed cells throughout the field. In this embodiment, the specific parameter value assigned to the mechanical property is: the elastic modulus of the rock is 38.5. Its Poisson's ratio is 0.23 and its density is 2680. By using a shared-node topology connection, the continuous transmission of momentum and energy flux of stress waves at the interface of meshes with different densities is ensured, avoiding the occurrence of artificial reflection faults; while the accurate assignment of elastic modulus, Poisson's ratio and density provides an essential constitutive carrier for the stiffness and mass matrix, resulting in a standard three-dimensional mesh model.

[0033] In this embodiment, ground temperature data, ventilation duration, ventilation speed, and airflow temperature are acquired. In the specific high geothermal engineering application scenario of this embodiment, the geothermal data, i.e., the original thermal reservoir temperature of the deep, undisturbed rock mass, is 65. The ventilation duration is 14400. This refers to 4.0 hours of ventilation before drilling and loading the explosive; the ventilation velocity is 4.5. The airflow temperature, i.e., the temperature of the cold airflow transported from outside the tunnel to the excavation face through the ventilation duct, is 22.0°C. .

[0034] The ground temperature data is applied to the outermost boundary to bring all grid cells to an initial thermal equilibrium state. The measured ground temperature data was used as a first-type isothermal boundary condition and directly applied to the outermost boundary of the three-dimensional mesh model. Subsequently, under the condition that the ventilation inside the tunnel was not turned on, the finite element heat conduction solver was called to perform a steady-state thermal equilibrium pre-solution, so that all mesh elements in the field reached a uniform initial thermal equilibrium state under the spontaneous conduction of the heat conduction control equation.

[0035] Because the surrounding rock strata beyond a radius of 5 are located deep within the Earth's crust, inside a large mountain, external airflow cannot disturb the temperature of such deep rock masses during the tunnel construction ventilation lifespan of several hours. Therefore, the outermost boundary is physically a typical isothermal infinity boundary. (The last sentence appears to be incomplete and possibly refers to a different topic.) Applying heat to the outermost boundary and solving for thermal equilibrium allows the three-dimensional mesh model to naturally form a continuous and undisturbed initial uniform temperature field that highly replicates the actual mountain environment. This eliminates the technical defects caused by artificially setting non-uniform initial values ​​for the internal mesh, which leads to numerical divergence or false transient heat flow in the early stages of time progression of the partial differential equation of heat conduction.

[0036] The equivalent hydraulic diameter of the tunnel cross-section is calculated based on the tunnel cross-sectional characteristic parameters. The kinematic viscosity, thermal diffusivity, and thermal conductivity of the air are determined based on the airflow temperature. The Reynolds number of the airflow inside the tunnel is obtained by dividing the product of the ventilation velocity and the equivalent hydraulic diameter by the kinematic viscosity. The Prandtl number is obtained by dividing the kinematic viscosity by the thermal diffusivity. The Nusselt number of the tunnel wall is obtained by substituting the Reynolds number and the Prandtl number into the turbulent convection heat transfer correlation. The convective heat transfer coefficient of the tunnel wall is calculated by dividing the product of the Nusselt number and the thermal conductivity by the equivalent hydraulic diameter. Subsequently, dimensionless hydrodynamic analytical calculations were performed on the forced convection heat transfer coefficient inside the tunnel. Based on the designed excavation span, excavation height, arch radius, and sidewall height, the pre-processing backend automatically calculated the excavation cross-sectional area of ​​the current straight-wall semi-circular arch horseshoe-shaped section to be 141.65. The wetted perimeter of the cross-section, i.e., the total closed perimeter of the tunnel's inner outline on the cross-section, is 45.31. Dividing four times the excavation cross-sectional area by the wetted perimeter of the cross-section yields the equivalent hydraulic diameter of the tunnel cross-section, which is calculated to be 12.50 mm. Simultaneously, based on the input airflow temperature, the system queries a standard air thermal property table to automatically determine the various thermophysical parameters of the cooling air fluid at the current temperature: the kinematic viscosity of the air is 1.55. The thermal diffusivity is 2.14. The thermal conductivity is 0.026. .

[0037] Furthermore, the Reynolds number of the airflow inside the tunnel is calculated by dividing the product of the ventilation velocity and the calculated equivalent hydraulic diameter by the kinematic viscosity; the specific value is 3.63. Meanwhile, the Prandtl number is calculated by dividing the kinematic viscosity of the air by the thermal diffusivity, and its specific value is 0.72.

[0038] Subsequently, the calculated Reynolds number and Prandtl number are substituted into the preset turbulent convection heat transfer correlation to perform power function calculation, and the Nusselt number, which characterizes the dimensionless temperature gradient of the tunnel wall, is calculated, with a specific value of 551.3.

[0039] The preset correlation for turbulent convection heat transfer is: In the formula, The Nusselt number of the tunnel wall. The Reynolds number is the airflow number inside the tunnel. For Prandtl numbers, It is a characteristic index of convective heat transfer state.

[0040] Finally, the convective heat transfer coefficient of the tunnel wall is calculated by multiplying the Nusselt number and the thermal conductivity by the equivalent hydraulic diameter, and the final calculated value is 24.5. .

[0041] The convective heat transfer coefficient and the airflow temperature are applied as convective heat transfer boundary conditions to the grid surface at the location of the tunnel inner profile. The calculated convective heat transfer coefficient and the airflow temperature are combined to form a convective heat transfer boundary condition, which is applied to the inner wall grid surface of the cavity at the location of the tunnel inner profile, thereby constructing a dynamic heat dissipation channel between the cooling airflow and the high-temperature bedrock on the grid physical surface.

[0042] Using the ventilation duration as the total solution time, the finite element software is called to calculate the nodal temperature values ​​of each grid element under the convective heat transfer boundary conditions, thus obtaining the temperature field distribution of the three-dimensional grid model.

[0043] The input ventilation duration is used as the total solution time, controlling the solver's sequential loop upon startup. At each time step, the solver calculates the heat flux carried away by the airflow based on the convective heat transfer coefficient on the tunnel wall grid surface and the real-time wall temperature difference, and dynamically updates the heat conduction difference matrix of each grid cell within the surrounding rock. This iterative loop continues until the total ventilation duration is reached. The solver then calculates and outputs the nodal temperature values ​​of each grid cell under the convective heat transfer boundary conditions, ultimately outputting the temperature field distribution of the three-dimensional mesh model.

[0044] By applying convective heat transfer boundary conditions and performing transient calculations, the physical process of continuous cold air blowing across the tunnel walls before drilling and charging high-temperature tunnels was reconstructed. This process caused a deep cooling of the surrounding rock from the shallow surface to the bedrock interior, resulting in a steep transient temperature gradient. The final output temperature field distribution not only marks the end of the thermal evolution during the pre-detonation preparation stage but also provides an aligned spatial three-dimensional thermal state data foundation for subsequent data processing, including capturing the pre-detonation temperature of each borehole assembly just before detonation and calculating the transient cooling temperature difference caused by ventilation.

[0045] Step 2: Obtain the explosive density and detonator delay sequence. Determine each borehole unit based on the borehole layout parameters. Extract the pre-detonation temperature of each borehole unit and calculate the detonation velocity. Calculate the blasting pressure based on the explosive density and convert it into a blasting pressure time history curve. According to the detonator delay sequence, apply the blasting pressure time history curve in stages to the grid boundary of each corresponding borehole layout position for integration and solution to obtain the blasting damage field.

[0046] In this embodiment, based on the location of the borehole and the diameter of the detonator borehole, the grid cells falling within the geometric boundary range of each borehole are divided into corresponding borehole cells. The positions of each borehole and the diameter of the detonator borehole are read from the control boundary input. In the three-dimensional mesh model, each borehole is represented in three-dimensional space as having a diameter of 0.045. Depth is 4.2 The spatial geometric boundary of the cylinder is determined by a spatial geometric topology inclusion relation algorithm. The algorithm retrieves and filters out the grid cells whose spatial geometric centroid or grid node coordinates fall completely within the spatial geometric boundary of each cylinder. These selected grid cells are then uniformly marked and divided into corresponding borehole cells in the internal data structure.

[0047] The average temperature of all grid nodes in each borehole unit is extracted as the pre-blast temperature of each borehole. The process iterates through each set of marked borehole units, querying and extracting the real-time temperature values ​​from the temperature field distribution calculated by all grid cells falling within that borehole unit. Subsequently, the average temperature of all grid nodes in that borehole unit is calculated, and the calculated statistical average is used as the pre-blast temperature of the corresponding borehole. In this embodiment, due to the significant spatial differences in the cooling depth of pre-blast convection heat transfer affecting boreholes at different locations, such as slotted holes near the center of the tunnel cross-section and peripheral smooth blasting holes near the outer tunnel wall, the extracted pre-blast temperatures are represented as a characteristic temperature array, with a specific preset value range of 53.5. Up to 61.2 .

[0048] The initial detonation velocity and critical deterioration temperature threshold of the explosive are obtained, and it is determined whether the pre-detonation temperature is greater than the critical deterioration temperature threshold. The initial detonation velocity of the explosive under standard room temperature static conditions and the critical degradation temperature threshold used to characterize the deterioration of the explosive's chemical stability are obtained. The extracted pre-detonation temperatures are then numerically compared with the critical degradation temperature thresholds to determine whether the explosive triggers a high-temperature performance degradation response.

[0049] When the pre-detonation temperature is less than the critical deterioration temperature threshold, the initial detonation velocity is directly used as the detonation velocity. When the pre-detonation temperature of a specific borehole is determined to be less than the critical degradation temperature threshold, it indicates that the surrounding rock of the borehole is subjected to excellent convective cooling, and the local temperature has dropped to within the safe thermal stability range of the explosive. The high temperature has not caused any physical weakening of the molecular structure of the explosive or the detonation reaction kinetic chain. At this time, the control does not trigger reduction, and the initial detonation velocity is directly input into the dynamic solver as the actual detonation velocity of the borehole.

[0050] When the pre-detonation temperature is greater than or equal to the critical deterioration temperature threshold, the temperature difference between the pre-detonation temperature and the critical deterioration temperature threshold is calculated, and the difference between the initial detonation velocity and the product of the temperature difference and the preset temperature sensitivity coefficient is taken as the reduced detonation velocity. When the pre-detonation temperature of a specific borehole is determined to be greater than or equal to the critical degradation temperature threshold, in this embodiment, taking a deep, incompletely cooled cut hole as an example, its calculated pre-detonation temperature is 58.5°C. Greater than 50.0 This indicates that the explosives loaded at that location underwent thermal degradation due to prolonged exposure to high-temperature bedrock. At this point, the temperature difference between the pre-detonation temperature and the critical degradation temperature threshold was calculated, yielding 8.5. Subsequently, a preset temperature sensitivity coefficient was invoked. Based on experimental data from high-temperature thermorheological reaction kinetics tests of industrial explosives, the preset temperature sensitivity coefficient was specifically set to 12.5. The product of the temperature difference and the temperature sensitivity coefficient is calculated, yielding an absolute decrease in the detonation velocity of 106.25. Finally, by subtracting the attenuation from the initial detonation velocity, the detonation velocity after borehole reduction is calculated, specifically 4093.75. .

[0051] The corrected blast pressure is obtained by multiplying the explosive density by the square of the detonation velocity.

[0052] The density of the input explosive is obtained, and its preset value is 1200.0. In the background algorithm, the explosive density is multiplied by the square of the detonation velocity of the corresponding borehole calculated in the previous steps to calculate the blast pressure after high-temperature thermal degradation correction. The final calculated pressure value is 20.11. The modified blasting pressure is then assigned as the characteristic peak value to the corresponding borehole entity boundary.

[0053] At the leading edge of a discontinuous detonation wave in hydrodynamics, the peak transient detonation pressure generated by explosive detonation is directly proportional to the initial charge density of the explosive and strictly proportional to the square of the detonation wave propagation velocity. Existing technologies, by blindly adopting a constant standard detonation velocity, result in significantly inflated detonation pressure values ​​under extreme high-temperature environments, causing severe distortion of the mechanical energy input at the simulation source. This invention, through the product operation of this square term, accurately pinpoints the peak microscopic dynamic pressure actually acting on the borehole wall surrounding rock under thermo-mechanical coupling, ensuring absolute determinism in subsequent dynamic calculations and detonation damage field output at the data link source.

[0054] In this embodiment, the modified blasting pressure is combined with the corresponding explosive detonation time to construct the blasting pressure time history curve; Input the required characteristic time control parameters and rock damage judgment threshold in the kinetic preprocessing control panel. In the specific numerical simulation scenario of this embodiment, the specific preset values ​​of each parameter are as follows: the corrected blasting pressure is 20.11. The specific time sequence for the explosive detonation is set as follows: the load rise phase duration is 20 seconds. The load descent phase lasts for 100 seconds. The total duration of the detonation gas effect was 120. The detonator delay sequence is specifically set as follows: the detonation delay of the first segment slotted hole is 0. The detonation delay for the second auxiliary hole is 50 seconds. The detonation delay of the third-stage peripheral blast holes is 110 seconds. The preset rock strength threshold specifically includes: a rock dynamic tensile strength threshold of 12.5. The dynamic compressive strength threshold of the rock is 165.0. .

[0055] The modified burst pressure, which serves as the peak pressure, is spatiotemporally coupled with the explosive detonation time series. On the time axis... =0 to =20 Within the range, the pressure value linearly increased from 0 to a maximum of 20.11. ;exist =20 to =120 Within the range, the pressure value increased from 20.11. The pressure decays exponentially or linearly to zero, thus constructing a classic triangular pulse-shaped time-history curve of the blast pressure at the data level. In finite element explicit dynamic analysis, a single static pressure peak cannot drive the transient fluctuation solver; the solver must rely on a continuous, time-varying force function to accurately calculate the momentum impulse of the work done by the blast shock wave on the surrounding rock of the borehole. This triangular waveform perfectly fits the real detonation thermodynamic physical scenario where the explosive undergoes an instantaneous chemical detonation reaction within a microscopic cavity, the pressure surges to a peak value, and then rapidly decays due to the expansion and leakage of detonation product gases.

[0056] According to the detonator delay sequence, the detonation time node is set for each of the borehole units; Based on the input detonator delay sequence, a unique detonation time node is assigned and set for each borehole unit set corresponding to each segment. Specifically, the detonation time node of the borehole unit set corresponding to the slotted hole is set to 0 of the global time history. Set the auxiliary hole set to 50. Set the surrounding light burst hole set to 110. Modern tunnel excavation strictly prohibits simultaneous blasting across the entire tunnel face; millisecond-level detonators and micro-delay blasting must be used. The first blasted slot hole creates a cavity in the center of the tunnel face, thus creating a crucial, unrestricted first free face for subsequent blasting of auxiliary and surrounding holes. This reduces the shear resistance of the fractured rock and suppresses seismic wave peaks propagating deep into the mountain.

[0057] At the corresponding detonation time node, the blast pressure time history curve is used as a transient dynamic load and applied to the grid boundary of the corresponding borehole unit; When the explicit time history cycle progresses to the corresponding detonation time node, the detonation pressure time history curve corresponding to the auxiliary hole is dynamically activated and treated as a transient dynamic load, applied strictly as a surface force to the inner wall mesh boundary of the corresponding auxiliary hole unit set. This reconstructs the mechanical transmission boundary process after detonation, where the high-pressure detonation fluid products directly and violently impact and push against the solid rock surface of the hole wall, realizing the conversion of chemical energy into solid mechanical kinetic energy and strain energy. This provides an indisputable source impact boundary condition for solving the partial differential equations of full-field elastic dynamics.

[0058] Solve the propagation process of the dynamic stress wave caused by the transient dynamic load in the grid element to obtain the principal stress of each grid element; A central difference explicit time integration algorithm is employed, with microsecond-level time steps, to calculate the outward anisotropic diffusion, reflection, and attenuation process of the impact stress wave excited by the transient dynamic load within a continuous grid cell. At each calculation time step, the three-dimensional strain tensor of each grid cell is extracted, and the corresponding three-dimensional real-time principal stress tensor is calculated through constitutive equation transformation. The first and third principal stresses are extracted and used as the real-time principal stresses of each grid cell at the current time step. The explosion energy travels through the rock mass grid in a wave-like manner, accompanied by geometric diffusion, with the high stress at the wavefront spreading outward. To dynamically monitor the crack initiation state of each spatial grid location at a microscopic level, its complex stress tensor must be transformed into three mutually perpendicular principal stress vectors, thus providing a clean data source for subsequent strength criterion comparisons.

[0059] When the principal stress exceeds the preset rock strength threshold, damage evolution is triggered, the blast damage variables of each grid cell are calculated, and the blast damage field is obtained.

[0060] At the end of each iteration at each time step, the real-time principal stress calculated for each grid cell is compared in real time with the preset rock strength threshold. When the real-time tensile stress of any grid cell is determined to be greater than 12.5... Or, the real-time compressive stress is greater than 165.0. This indicates that the rock at that location was violently torn or compressed by the explosion wave, and the strain exceeded the limit of brittle fracture of the material, thus automatically triggering the damage evolution in the preset material constitutive model. In this embodiment, the damage evolution follows the following dynamic crack propagation accumulation equation: In the formula, For the aforementioned blast damage variable, To accumulate equivalent plastic tensile strain, The critical plastic strain threshold. and For explosive damage material constants, For the analysis time, It is an exponential function. It is a singular function.

[0061] The legal physical range of blast damage variables in continuum mechanics is strictly limited to Within the equation, 0 represents intact rock mass, and 1 represents complete fragmentation and loss of stiffness. This is due to the material constant... And time exponent term The output value is naturally locked to Within the interval, thus making The derived blast damage variable spontaneously converges to [a specific value] during linear time history progression. The legal interval completely eliminates the underlying defect of numerical iteration causing damage to exceed the limit and leading to the collapse of the finite element solver due to the singularity of the stiffness matrix from the mathematical framework.

[0062] When the rock mesh elements are subjected to the cumulative equivalent plastic tensile strain generated by the explosive stress wave... If the material's inherent elastic limit has not yet been exceeded, the output inside the parentheses is negative, and the program automatically forces it to be truncated to zero. This perfectly recreates the real physical scenario where rocks are in a purely elastic or slightly plastic stable state in the far field of an explosion, without producing irreversible mechanical brittle damage.

[0063] In this embodiment, the specific preset control value of the critical plastic strain threshold is set as follows: This parameter characterizes the plastic inflection point in hard surrounding rock materials under dynamic high strain rate tensile loading, where internal microcracks begin to unstablely propagate, converge, and transition to irreversible macroscopic damage. Due to the extremely brittle characteristics of hard granite in high geothermal environments, its elastic limit range is extremely narrow, and its ultimate tensile deformation capacity is extremely low. Based on indoor Hopkinson bar dynamic rock splitting test data, The dynamic plastic tensile initiation strain of hard rock of grade II hardness is concentrated at... to Within the range. Therefore, this invention locks in its core target value. Once the plastic deformation of the mesh exceeds this watershed, it means that substantial tearing has occurred at the rock grain boundaries, and the damage process is legally activated.

[0064] In this embodiment, the explosive damage material constant is used as a scaling factor. The specific preset control value is set as follows This constant term primarily controls the rapid evolution rate of the damage variable towards the ultimate failure state after the plastic strain exceeds the limit. Since explosive explosions are typical transient ultra-high strain rate events, the interaction time of the explosion shock wave and dynamic stress wave within the borehole wall surrounding rock is extremely short, typically on the order of microseconds to milliseconds. To facilitate the explicit dynamic solver's accurate capture of stiffness degradation within extremely short physical time steps, the rate constant... The magnitude of the explosion damage variable must be extremely high to allow the near-field mesh elements subjected to severe dynamic mechanical tension to accumulate rapidly within the total dynamic analysis time. The destruction state is thus accurately matched to the transient time scenario of explosive detonation in the time dimension.

[0065] In this embodiment, the explosive damage material constant is used as an evolution index. Specific preset control value settings This parameter determines the nonlinear curvature of the damage accumulation curve. If This represents a simple linear damage evolution and fails to reflect the physical nature of nonlinear brittle fracture in rocks. Under real-world conditions, when rocks are torn apart by explosive dynamic loads, the propagation of microcracks exhibits a nonlinear S-shaped or exponentially accelerated characteristic of "initial incubation, abrupt acceleration in the middle stage, and saturation in the later stage." When the coefficient is set to 1.8, the equation mathematically forms a power function with an accelerating slope, meaning that the stiffness reduction effect caused by the convergence and penetration of microcracks is nonlinearly amplified as the excess plastic strain increases. This perfectly aligns with the microscopic damage characteristic in dynamic damage mechanics that "mechanical damage to rocks deteriorates at an accelerated rate with plastic expansion," eliminating numerical smoothing errors caused by linear dissipation from within the structure and significantly improving the spatial prediction accuracy of the initial blast damage field.

[0066] Based on the excess plastic strain increment of the grid cell exceeding the threshold, a nonlinear damage accumulation iteration is performed to calculate the blast damage variable of the grid cell at the current time step. When the explicit total analysis time expires, the blast damage variables finally accumulated by all grid cells in the entire field are combined into a matrix sequence according to the spatial three-dimensional coordinate topological relationship, and the final output is the standard blast damage field.

[0067] Under strong dynamic impact, hard rock undergoes irreversible microcrack convergence and stiffness softening, and the values ​​of damage variables directly correspond to the spatial density of rock microcracks. By outputting the spatial combination of damage variables across the entire grid, the originally abstract dynamic stress wave propagation process is successfully and intuitively transformed into a three-dimensional, visualized blast damage field containing "fracture zone, crack propagation zone, and elastically undamaged zone," providing an isomorphic and aligned data carrier for subsequent introduction of ventilation-induced rapid cooling thermal damage and execution of nonlinear matrix superposition.

[0068] Step 3: Apply a post-blast ventilation cooling heat transfer boundary to the three-dimensional mesh model, calculate the mesh element temperature at each cooling time step, calculate the thermal strain component according to the thermoelastic constitutive equation, calculate the thermal damage variable when the thermal strain component exceeds the limit, and nonlinearly superimpose the thermal damage variable and the blast damage variable of the blast damage field to obtain the total damage of each mesh element.

[0069] In this embodiment, the convective heat transfer coefficient and the airflow temperature are applied to the grid surface at the location of the tunnel inner profile to calculate the real-time temperature of each grid cell at the current cooling step. The convective heat transfer coefficient of the tunnel wall and the input airflow temperature are combined to form a third type of transient fluid heat transfer boundary load. After the blasting dynamics calculation is completed and the debris is removed, this heat transfer boundary load is precisely and in real-time applied to the exposed cavity wall grid surface at the location of the tunnel inner profile. Subsequently, a thermal unsteady-state time-step progressive solution is initiated to calculate the real-time temperature of each grid cell in the entire field at the current cooling time step. In the simulation scenario of this embodiment, taking the current cooling time step to half an hour of ventilation as an example, the real-time temperature of a certain circumferential grid cell near the tunnel wall surface is solved and output as a specific value of 32.4. .

[0070] The blasting and excavation disrupted the original high-temperature thermal equilibrium deep within the mountain's surrounding rock. (The last sentence appears to be incomplete and unrelated to the preceding text.) When a large volume of cooling airflow directly blows onto the newly exposed high-temperature wall surface after the explosion, the bedrock surface releases energy into the air through strong convective heat flux. By continuously solving for the real-time temperature of each grid cell at each specific cooling time step, the transient evolution process of the high-temperature surrounding rock gradually decreasing in temperature from the outside to the inside and the drastic temperature fluctuations inside the heat transfer boundary layer can be captured with extremely high temporal resolution, thus pinpointing the precise temperature state source for subsequent thermal stress deformation calculations.

[0071] Subtracting the real-time temperature from the pre-explosion temperature of each grid cell yields the transient cooling temperature difference of each grid cell at the current cooling time step. At the current cooling time step, retrieve the pre-detonation temperature of each grid cell at the instant before detonation. Using a pure numerical subtraction algorithm, subtract the real-time temperature calculated at the current time step from the pre-detonation temperature of each grid cell to calculate the transient cooling temperature difference of each grid cell at the current cooling time step.

[0072] Based on the thermoelastic constitutive equation, the thermal strain components of each grid cell in the corresponding spatial axis direction are calculated by multiplying the transient cooling temperature difference by the preset linear expansion coefficient of the rock. Obtain the preset linear expansion coefficient of the rock, which is based on the indoor thermal expansion test results of hard granite, and the specific preset physical constant value is set to 7.5. In the background solution matrix, the transient cooling temperature difference calculated for each grid cell is multiplied by the preset linear expansion coefficient to calculate the thermal strain component of each grid cell along the corresponding spatial axis. For this specific grid cell in this embodiment, the calculated thermal strain component values ​​are as follows: .

[0073] In a macroscopically unconstrained state, solids spontaneously contract upon cooling. However, in tunnel engineering, the surface mesh cells are tightly topologically connected to the deep, unaffected bedrock mesh. When the surface mesh, driven by rapid cooling, attempts to contract inward, it is inevitably constrained and stretched by the stiffness of the large, uncooled rock mass in the far field. The thermal strain component calculated through this product operation is, in mechanical essence, a typical example of constrained contraction-tension strain. It successfully and seamlessly transforms the microscopic gas-fluid heat dissipation effect into the macroscopic tensile deformation caused by deformation constraints within the surrounding rock solid structure, providing the most direct stiffness matrix independent variable for subsequent crack damage assessment.

[0074] The thermal strain component is compared with a preset rock ultimate tensile strain threshold for determination; The preset rock ultimate tensile strain threshold is retrieved, which characterizes the elastic critical point at which the current hard rock material will undergo tensile fracture. Based on tensile fracture tests, its specific preset control value is set to 120. The thermal strain components calculated for each mesh element along its corresponding spatial axis are compared numerically item by item with the ultimate tensile strain threshold. In this embodiment, if the thermal strain component of a specific mesh element is determined to be significantly greater than the ultimate tensile strain threshold, then the mesh element is determined to have triggered cold shrinkage and thermal stretching cracking failure along that spatial axis.

[0075] When the thermal strain component is greater than the ultimate tensile strain threshold, the thermal damage variable of each grid element is calculated by nonlinear accumulation based on the excess strain of the thermal strain component exceeding the ultimate tensile strain threshold.

[0076] When it is determined that the thermal strain component of a certain mesh element is greater than the ultimate tensile strain threshold, a difference operation is performed to calculate the excess strain of the thermal strain component exceeding the ultimate tensile strain threshold. In this embodiment, the excess strain of the specific element is specifically... Subsequently, a nonlinear cumulative calculation was performed on the excess strain: In the formula, For the aforementioned thermal damage variable, The thermal damage sensitivity coefficient, The thermal brittleness index, To exceed the required response, It is an exponential function.

[0077] In continuum mechanics and fracture mechanics, the range of values ​​for thermal damage variables must be strictly limited. Internally. Due to the thermal damage sensitivity coefficient. and excess strain exponent term The output value is always negative or zero, resulting in The evaluation range of the function converges to... Within the interval. Passing through. The mathematical framework transformation allows the calculated thermal damage variable to spontaneously and absolutely converge to a certain value during the time progression of ventilation and heat dissipation. Within the legal physical boundaries, the underlying technical defect of numerical solvers, which leads to stiffness matrix degradation and computational divergence collapse due to overflow of damage variables, is solved from the source of the formula structure.

[0078] This study reconstructs the transient dynamics of shallow brittle crystalline fracture initiation in high-temperature rock masses subjected to forced convection heat transfer by cold wind. Once the strain response induced by sudden cooling crosses a critical point, the thermal damage variable exhibits a nonlinear evolution trajectory of "rapid initial increase, nonlinear transition in the middle stage, and tendency to saturate in the later stage" with increasing excess strain. This mathematical form closely matches the macroscopic physical manifestations of internal microscopic thermal stress concentration and instantaneous unstable expansion of the fracture network in rock materials under thermal shock.

[0079] In this embodiment, the excess strain is the difference between the thermal strain component of the grid element in the corresponding spatial axis direction and the preset rock ultimate tensile strain threshold. Specifically, the calculated net excess cold shrinkage tensile deformation is: According to the thermoelastic boundary conditions, when high-temperature bedrock encounters strong cooling airflow, the surface rock spontaneously undergoes contraction and deformation. However, this is limited by the stiffness of the uncooled rock mass deep within the mountain, thus transforming into intense restrictive contraction-tension thermal stress. As long as the thermal strain component does not exceed the inherent tensile-brittle fracture threshold of the rock, only reversible elastic contraction occurs within the rock, and the thermal damage variable remains constant at 0. This invention extracts the excess strain obtained by subtracting the two as the core independent variable driving damage evolution, precisely identifying the effective energy release source that truly induces irreversible thermal tensile cracking and failure in the shallow surface layer of the rock mass within the data chain.

[0080] In this embodiment, the preset control value of the thermal damage sensitivity coefficient, which serves as a proportional adjustment factor, is set to 4500.0. This constant term primarily controls the sensitivity and severity of the thermal damage curve to fluctuations in excess strain. Due to the extremely low thermal conductivity and strong "tense" thermal brittleness of rock materials, the tensile failure response caused by the temperature gradient induced by ventilation cooling in the shallow surface layer of the cave wall is extremely rapid and sensitive, representing a typical high-gradient transient stress release event. To accurately capture material degradation caused by thermal cracking in each cooling step using the fluid-thermal-solid coupled solver, this sensitivity coefficient must be maintained in the high order of thousands. Setting it to 4500.0 ensures that when the excess cooling tensile strain only experiences a tiny jump at the micrometer level, the program can sensitively drive the thermal damage variable upwards to accumulate to a certain value. The reasonable thermal degradation state was reproduced, which is a real physical engineering scenario of instantaneous thermal spalling of the tunnel wall after the tunnel is exposed by excavation at high ground temperature.

[0081] In this embodiment, the preset control value of the thermal brittleness index, which serves as the nonlinear curvature control term, is set to 1.05. This index term determines the nonlinear curvature when excess strain is converted into thermal damage. If This means that the damage dissipates linearly or at a decelerating rate with cooling and contraction, which severely contradicts the physical reality that brittle rocks experience accelerated convergence and connection of cracks under thermal shock. Based on the stress-strain response curves of indoor unsteady-state temperature impact tensile fracture tests on rocks, hard granite exhibits a weak nonlinear accelerated propagation characteristic of internal microcracks when subjected to drastic temperature drops. By... The precise value of 1.05 ensures that the evolution slope of the equation exhibits a mild nonlinear accelerated softening as the excess strain increases. This not only precisely matches the microscopic evolution law of thermal brittle fracture in rocks under high geothermal excavation faces, but also ensures the continuity of the derivative of the curve at zero, significantly improving the spatial accuracy and numerical stability of the secondary thermal damage variable calculation output.

[0082] Finally, the thermal damage variable corresponding to the grid cell is calculated and output, and its value is strictly limited to between 0 and 1. In this embodiment, the specific value of the thermal damage variable finally calculated and output by the cell is 0.28.

[0083] The crack propagation induced by thermal shock in rocks is not a simple linear and uniform process of increasing stress, but rather a nonlinear abrupt change accompanied by the penetration of microcracks and accelerated softening of stiffness. By extracting the core independent variable of "excess strain" and performing nonlinear accumulation, the damage evolution of rocks caused by excessive cooling contraction is perfectly and quantitatively characterized, leading to a sharp reduction in the effective load-bearing area and a decrease in micro-fracture toughness within the material. Ultimately, the thermal damage variables output by each grid cell across the entire field achieve spatial and typological isomorphism with the explosive mechanical damage variables, providing a fully ready input source for the nonlinear matrix superposition of multiple fields of dual damage.

[0084] In this embodiment, for each grid cell in the entire field, the corresponding blast damage variable and thermal damage variable are extracted respectively; Subtract the blast damage variable and the thermal damage variable from the constant one respectively to obtain the first integrity coefficient characterizing the proportion of grid cells that are not damaged after blasting, and the second integrity coefficient characterizing the proportion of grid cells that are not damaged after rapid cooling. Subtract the blast damage variable and thermal damage variable corresponding to each grid cell from the obtained constant to calculate the first integrity coefficient representing the proportion of the grid cell that is not damaged after blasting, and the second integrity coefficient representing the proportion of the grid cell that is not damaged after rapid cooling.

[0085] Multiply the first complete coefficient with the second complete coefficient to obtain the total complete coefficient; The calculated first integrity coefficient and the second integrity coefficient are multiplied to obtain the total integrity coefficient, which characterizes the proportion of the effective bearing area remaining after the double multiphysics field overlap degradation of the mesh cell. For the specific mesh cell in this embodiment, the first integrity coefficient and the second integrity coefficient are multiplied to determine and output the specific value of the total integrity coefficient of the specific mesh cell as 0.396. The total damage of each mesh cell is obtained by subtracting the total integrity coefficient from the constant.

[0086] Subtracting the total integrity coefficient of the corresponding mesh element calculated in the previous step from the constant, the final total damage of the mesh element is calculated by reverse calculation. For this specific mesh element in this embodiment, the specific parameter value of the total damage of the mesh element is output as 0.604. By traversing all mesh elements in the entire field and repeating the above calculation, a uniform composite total damage matrix for the entire field is finally generated.

[0087] Step 4: Divide the tunnel inner profile into multiple detection sections along its geometric axis. Using the intersection of the tunnel inner profile's geometric axis and the detection section as the origin of the polar coordinate system, perform a radial search along the azimuth angle into the depth of the rock mass. Retrieve the farthest grid cell that satisfies the condition that the total damage is greater than a preset critical damage threshold. Calculate the geometric distance of the farthest grid cell from the tunnel inner profile as the surrounding rock damage depth at that azimuth angle.

[0088] In this embodiment, a geometric axis perpendicular to the inner contour of the tunnel is used to divide multiple detection sections. On each detection section, a polar coordinate system is established with the intersection of the geometric axis and the corresponding detection section as the origin of the polar coordinate system. Read the global Cartesian coordinate system, along the... The axis, i.e., the geometric axis extending from the inner contour of the tunnel, is set according to a preset spatial axial step size. In this embodiment, the preset value of the spatial axial step size is set to 2.0. Throughout the entire 40.0 Within the long model range, a series of perpendicular to The plane of the axis geometrically cuts the three-dimensional mesh model containing the total damage matrix of all mesh elements, thereby dividing it into multiple parallel transverse inspection sections along the entire tunneling advance direction.

[0089] The azimuth angle is cyclically changed with a preset angle step size, and a search ray is emitted from the origin of the polar coordinate system to the deep part of the rock mass along the current azimuth angle; On each of the segmented transverse detection sections, the geometric axis of the tunnel's inner profile and the spatial geometric three-dimensional intersection point with the current detection section plane are calculated. This intersection point is then used as the origin of the polar coordinate system in the current plane, and an independent local two-dimensional polar coordinate system is established within each corresponding transverse detection section, defining the horizontal direction to the right as the initial zero direction of the local polar coordinate polar axis.

[0090] With a preset angle step size, in this embodiment, the specific preset value of the angle step size is set as follows: In the local polar coordinate system Control the azimuth angle of the cyclic variable during travel within the entire circumferential angle range. ,Right now At the current azimuth angle, starting from the origin of the polar coordinate system, a one-dimensional search ray is emitted from the current transverse detection section plane towards the deeper part of the surrounding rock mass.

[0091] Along the search ray, all grid cells traversed by the search ray are traversed, and the total damage corresponding to each grid cell is extracted; Using a ray-polyhedron intersection penetration topology algorithm, the algorithm extends outward along the search ray at the current azimuth angle, sequentially traversing and identifying all discrete mesh elements traversed by the search ray. Then, the data reading thread extracts the corresponding total damage values ​​of these traversed mesh elements from the final calculated composite total damage matrix. The total damage is compared with a preset critical damage threshold to filter out damaged mesh cells whose total damage is greater than the critical damage threshold. The preset critical damage threshold is retrieved. Based on the indoor damage rheology and residual strength tests of high-temperature hard granite, in this embodiment, the preset control value of the critical damage threshold is set to 0.25, representing the characteristic limit scalar value of the rock losing its self-bearing capacity. In the judgment algorithm, the total damage value extracted from each grid cell in the above-mentioned penetration chain is compared with the critical damage threshold. When it is determined that the total damage of some cells is greater than 0.25, the grid is determined to have been filtered and classified into a set of damaged grid cells.

[0092] Among the screened damaged mesh cells, the target mesh cell that is farthest from the origin of the polar coordinate system is identified; Among all the damaged mesh cells selected at the current azimuth angle, the spatial distance between each cell and the origin of the polar coordinate system is calculated by comparing the three-dimensional spatial coordinates of its center. The cells are then sorted from largest to smallest to uniquely identify the mesh cell with the greatest geometric distance from the origin of the polar coordinate system at the spatial topology level. This mesh cell is defined as the target mesh cell. The center geometric coordinates of this target mesh cell are then extracted.

[0093] Obtain the center coordinates of the target grid cell, and obtain the boundary coordinates of the intersection point of the tunnel inner profile and the search ray at the current azimuth angle; Calculate the spatial geometric distance between the center coordinates and the boundary coordinates, and output the spatial geometric distance as the surrounding rock damage depth at the current azimuth angle on the current detection section.

[0094] The intersection point of the search ray at the current azimuth angle and the two-dimensional closed curve of the tunnel inner profile, which serves as the excavation reference tunnel wall boundary, is determined in the plane. The spatial coordinates of this intersection point are extracted and defined as boundary coordinates. Subsequently, the standard three-dimensional spatial geometric distance formula is used to calculate the spatial geometric distance between the geometric coordinates of the locked target grid cell center and the boundary coordinates. This calculated spatial geometric distance is then used as the final surrounding rock damage depth at the current azimuth angle on the current transverse detection section, and the data is output as a matrix format stream of a text file.

[0095] Table 1: Comparison of Multi-Field Coupled Composite Damage between Existing Technology and This Application This technical solution, through deep coupling and rigorous solution of multiphysics fields, intuitively demonstrates the three core technical advantages of this invention under extreme high-temperature tunnel conditions: First, it has the ability to adapt to the environment and reconstruct the dynamic load source with high fidelity. Compared with the energy input distortion caused by blindly using constant detonation parameters in existing technologies, this method can keenly capture the thermal degradation effect of high temperature on the detonation performance of explosives. By extracting the pre-detonation temperature and dynamically correcting the time history curves of detonation velocity and blasting pressure, it fundamentally eliminates the source distortion of the initial mechanical boundary and realizes the deep and fundamental coupling of heat conduction and dynamic evolution. Second, it has the ability to completely capture the multiphysics field composite fracture mechanism and nonlinearly extrapolate it. Traditional assessment methods often attribute damage solely to mechanical tensile failure caused by blasting dynamic load, thus leaving dangerous secondary disaster blind spots. This method introduces secondary cooling-contraction and thermal-tension cracks induced by post-blast cold air convection and nonlinearly superimposes mechanical damage and sudden cooling thermal damage to intercept and reproduce the composite failure zone under multi-field interactive evolution, ensuring the absolute integrity of the surrounding rock deterioration simulation in terms of physical mechanisms. Finally, it possesses a high-precision topological mapping capability from three-dimensional abstract damage to engineering quantitative space. Addressing the fuzziness and coarseness of traditional methods in extracting failure boundaries, this method relies on multi-section segmentation and polar coordinate amplitude radial retrieval technology, using the critical damage threshold as a guideline to target and lock the farthest target grid cell that meets the instability conditions. It successfully bridges the data gap between microscopic grid damage variables and macroscopic safety control distances, achieving a substantial leap from abstract three-dimensional failure maps to precise on-site support parameter guidance.

[0096] Please see Figure 2The present invention also provides a system for calculating the depth of damage to surrounding rock in high-temperature tunnel blasting excavation. This system is used to implement the aforementioned method for calculating the depth of damage to surrounding rock in high-temperature tunnel blasting excavation, and includes: The model building module is used to obtain tunnel cross-sectional characteristic parameters and borehole layout parameters, determine the tunnel inner profile based on the tunnel cross-sectional characteristic parameters, build a three-dimensional spatial model, perform meshing on the three-dimensional spatial model to obtain a three-dimensional mesh model, apply pre-blast heat transfer boundary conditions to the three-dimensional mesh model, and calculate the temperature field distribution of the three-dimensional mesh model. The blasting damage module is used to acquire explosive density and detonator delay sequence, determine each borehole unit based on the borehole layout parameters, extract the pre-blast temperature of each borehole unit, calculate the detonation velocity, calculate the blasting pressure based on the explosive density and convert it into a blasting pressure time history curve, and apply the blasting pressure time history curve in stages to the grid boundary of each corresponding borehole layout position according to the detonator delay sequence for integration solution to obtain the blasting damage field; The heat loss superposition module is used to apply the post-explosion ventilation cooling heat transfer boundary to the three-dimensional mesh model, calculate the mesh unit temperature at each cooling time step, calculate the thermal strain component according to the thermoelastic constitutive equation, calculate the thermal damage variable when the thermal strain component exceeds the limit, and nonlinearly superimpose the thermal damage variable and the blast damage variable of the blast damage field to obtain the total damage of each mesh unit. The depth retrieval module is used to divide the tunnel inner profile into multiple detection sections along the geometric axis. Taking the intersection of the geometric axis of the tunnel inner profile and the detection section as the origin of the polar coordinate system, it performs a radial search along the azimuth angle into the deep rock mass, retrieves the farthest grid cell that satisfies the condition that the total damage is greater than a preset critical damage threshold, and calculates the geometric distance of the farthest grid cell from the tunnel inner profile as the surrounding rock damage depth at that azimuth angle.

[0097] The above formulas are all dimensionless calculations. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0098] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented in software, the above embodiments can be implemented, in whole or in part, as a computer program product. Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented by electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution.

[0099] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment, depending on actual needs.

[0100] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A method for calculating the damage depth of surrounding rock during high-temperature tunnel blasting excavation, characterized in that, The specific steps include: Step 1: Obtain tunnel cross-sectional characteristic parameters and borehole layout parameters, determine the tunnel inner profile based on the tunnel cross-sectional characteristic parameters, establish a three-dimensional spatial model, perform meshing on the three-dimensional spatial model to obtain a three-dimensional mesh model, apply pre-blast heat transfer boundary conditions to the three-dimensional mesh model, and calculate the temperature field distribution of the three-dimensional mesh model. Step 2: Obtain the explosive density and detonator delay sequence. Based on the borehole layout parameters, determine each borehole unit, extract the pre-detonation temperature of each borehole unit, calculate the detonation velocity, calculate the blasting pressure based on the explosive density and convert it into a blasting pressure time history curve. According to the detonator delay sequence, apply the blasting pressure time history curve in stages to the grid boundary of each corresponding borehole layout position for integration and solution to obtain the blasting damage field. Step 3: Apply a post-blast ventilation cooling heat transfer boundary to the three-dimensional mesh model, calculate the mesh element temperature at each cooling time step, calculate the thermal strain component according to the thermoelastic constitutive equation, calculate the thermal damage variable when the thermal strain component exceeds the limit, and nonlinearly superimpose the thermal damage variable and the blast damage variable of the blast damage field to obtain the total damage of each mesh element. Step 4: Divide the tunnel inner profile into multiple detection sections along its geometric axis. Using the intersection of the tunnel inner profile's geometric axis and the detection section as the origin of the polar coordinate system, perform a radial search along the azimuth angle into the depth of the rock mass. Retrieve the farthest grid cell that satisfies the condition that the total damage is greater than a preset critical damage threshold. Calculate the geometric distance of the farthest grid cell from the tunnel inner profile as the surrounding rock damage depth at that azimuth angle.

2. The method for calculating the damage depth of surrounding rock during high-temperature tunnel blasting excavation according to claim 1, characterized in that: The tunnel cross-sectional characteristic parameters include: designed excavation span, excavation height, arch radius, sidewall height, tunnel axis extension direction, and length; The parameters for the borehole layout include: borehole diameter, borehole location, and borehole outward slope angle. Taking the intersection of the tunnel floor centerline and the tunnel face as the origin, the direction of the designed excavation span is defined as... The axis, defining the direction of the excavation height as The direction of the tunnel axis extension is defined as... Establish a global coordinate system based on the axis, and determine the tunnel's inner outline based on the excavation span, excavation height, sidewall height, and arch radius; Using the center of the tunnel's inner profile as a reference, along the... The positive and negative directions of the axis and the aforementioned The axis extends symmetrically outward in both positive and negative directions by five times the radius of the arch, serving as the boundary dimension, and extends along the axis... Stretch along the axial direction to the stated length to obtain a cubic solid, while simultaneously stretching along the stated length... The tunnel inner contour is stretched along the axial direction according to the length to obtain the tunnel excavation entity. Boolean subtraction is performed between the cube entity and the tunnel excavation entity to obtain a three-dimensional space model. In the area extending outward from each of the aforementioned borehole locations to a predetermined multiple of the borehole diameter, and in the area where the tunnel inner contour extends a predetermined distance into the rock mass, the grid is refined, and the size of the grid unit in this area is controlled within a first predetermined size range. A gradient mesh is performed from the outer boundary of the mesh refinement region to the outermost boundary of the cube entity. The size of the mesh cells in this region is adaptively enlarged as the distance from the outer boundary of the mesh refinement region increases, and is enlarged to a second preset size range at the outermost boundary, wherein the lower limit of the second preset size range is greater than the upper limit of the first preset size range. A common-node topological connection is made at the interface between the refined mesh region and the gradually subdivided mesh region, and the elastic modulus, Poisson's ratio and density of the rock are assigned as mechanical properties to the whole field mesh elements to obtain a three-dimensional mesh model.

3. The method for calculating the damage depth of surrounding rock during high-temperature tunnel blasting excavation according to claim 1, characterized in that: Acquire ground temperature data, ventilation duration, ventilation speed, and airflow temperature; The ground temperature data is applied to the outermost boundary to bring all grid cells to an initial thermal equilibrium state. The equivalent hydraulic diameter of the tunnel cross-section is calculated based on the tunnel cross-sectional characteristic parameters. The kinematic viscosity, thermal diffusivity, and thermal conductivity of the air are determined based on the airflow temperature. The Reynolds number of the airflow inside the tunnel is obtained by dividing the product of the ventilation velocity and the equivalent hydraulic diameter by the kinematic viscosity. The Prandtl number is obtained by dividing the kinematic viscosity by the thermal diffusivity. The Nusselt number of the tunnel wall is obtained by substituting the Reynolds number and the Prandtl number into the turbulent convection heat transfer correlation. The convective heat transfer coefficient of the tunnel wall is calculated by dividing the product of the Nusselt number and the thermal conductivity by the equivalent hydraulic diameter. The convective heat transfer coefficient and the airflow temperature are applied as convective heat transfer boundary conditions to the grid surface at the location of the tunnel inner profile. Using the ventilation duration as the total solution time, the finite element software is called to calculate the nodal temperature values ​​of each grid element under the convective heat transfer boundary conditions, thus obtaining the temperature field distribution of the three-dimensional grid model.

4. The method for calculating the damage depth of surrounding rock during high-temperature tunnel blasting excavation according to claim 1, characterized in that: Based on the location of the borehole and the diameter of the detonator borehole, the grid cells falling within the geometric boundary range of each borehole are divided into corresponding borehole cells. The average temperature of all grid nodes in each borehole unit is extracted as the pre-blast temperature of each borehole. The initial detonation velocity and critical deterioration temperature threshold of the explosive are obtained, and it is determined whether the pre-detonation temperature is greater than the critical deterioration temperature threshold. When the pre-detonation temperature is less than the critical deterioration temperature threshold, the initial detonation velocity is directly used as the detonation velocity. When the pre-detonation temperature is greater than or equal to the critical deterioration temperature threshold, the temperature difference between the pre-detonation temperature and the critical deterioration temperature threshold is calculated, and the difference between the initial detonation velocity and the product of the temperature difference and the preset temperature sensitivity coefficient is taken as the reduced detonation velocity. The corrected blast pressure is obtained by multiplying the explosive density by the square of the detonation velocity.

5. The method for calculating the damage depth of surrounding rock during high-temperature tunnel blasting excavation according to claim 4, characterized in that: The time history curve of the blasting pressure is constructed by combining the corrected blasting pressure with the corresponding explosive detonation time. According to the detonator delay sequence, the detonation time node is set for each of the borehole units; At the corresponding detonation time node, the blast pressure time history curve is used as a transient dynamic load and applied to the grid boundary of the corresponding borehole unit; Solve the propagation process of the dynamic stress wave caused by the transient dynamic load in the grid element to obtain the principal stress of each grid element; When the principal stress exceeds the preset rock strength threshold, damage evolution is triggered, the blast damage variables of each grid cell are calculated, and the blast damage field is obtained.

6. The method for calculating the damage depth of surrounding rock during high-temperature tunnel blasting excavation according to claim 3, characterized in that: The convective heat transfer coefficient and the airflow temperature are applied to the grid surface at the location of the tunnel profile to calculate the real-time temperature of each grid cell at the current cooling step. Subtracting the real-time temperature from the pre-explosion temperature of each grid cell yields the transient cooling temperature difference of each grid cell at the current cooling time step. Based on the thermoelastic constitutive equation, the thermal strain components of each grid cell in the corresponding spatial axis direction are calculated by multiplying the transient cooling temperature difference by the preset linear expansion coefficient of the rock. The thermal strain component is compared with a preset rock ultimate tensile strain threshold for determination; When the thermal strain component is greater than the ultimate tensile strain threshold, the thermal damage variable of each grid element is calculated by nonlinear accumulation based on the excess strain of the thermal strain component exceeding the ultimate tensile strain threshold.

7. The method for calculating the damage depth of surrounding rock during high-temperature tunnel blasting excavation according to claim 1, characterized in that: For each grid cell in the entire field, the corresponding blast damage variable and thermal damage variable are extracted respectively. Subtract the blast damage variable and the thermal damage variable from the constant one respectively to obtain the first integrity coefficient characterizing the proportion of grid cells that are not damaged after blasting, and the second integrity coefficient characterizing the proportion of grid cells that are not damaged after rapid cooling. Multiply the first complete coefficient with the second complete coefficient to obtain the total complete coefficient; The total damage of each mesh cell is obtained by subtracting the total integrity coefficient from the constant.

8. The method for calculating the damage depth of surrounding rock during high-temperature tunnel blasting excavation according to claim 1, characterized in that: Multiple detection sections are divided by a geometric axis perpendicular to the inner contour of the tunnel. On each detection section, a polar coordinate system is established with the intersection of the geometric axis and the corresponding detection section as the origin of the polar coordinate system. The azimuth angle is cyclically changed with a preset angle step size, and a search ray is emitted from the origin of the polar coordinate system to the deep part of the rock mass along the current azimuth angle; Along the search ray, all grid cells traversed by the search ray are traversed, and the total damage corresponding to each grid cell is extracted; The total damage is compared with a preset critical damage threshold to filter out damaged mesh cells whose total damage is greater than the critical damage threshold. Among the screened damaged mesh cells, the target mesh cell that is farthest from the origin of the polar coordinate system is identified; Obtain the center coordinates of the target grid cell, and obtain the boundary coordinates of the intersection point of the tunnel inner profile and the search ray at the current azimuth angle; Calculate the spatial geometric distance between the center coordinates and the boundary coordinates, and output the spatial geometric distance as the surrounding rock damage depth at the current azimuth angle on the current detection section.

9. A system for calculating the depth of damage to surrounding rock during high-temperature tunnel blasting excavation, wherein the system is used to implement the method for calculating the depth of damage to surrounding rock during high-temperature tunnel blasting excavation as described in any one of claims 1-8, characterized in that, include: The model building module is used to obtain tunnel cross-sectional characteristic parameters and borehole layout parameters, determine the tunnel inner profile based on the tunnel cross-sectional characteristic parameters, build a three-dimensional spatial model, perform meshing on the three-dimensional spatial model to obtain a three-dimensional mesh model, apply pre-blast heat transfer boundary conditions to the three-dimensional mesh model, and calculate the temperature field distribution of the three-dimensional mesh model. The blasting damage module is used to acquire explosive density and detonator delay sequence, determine each borehole unit based on the borehole layout parameters, extract the pre-blast temperature of each borehole unit, calculate the detonation velocity, calculate the blasting pressure based on the explosive density and convert it into a blasting pressure time history curve, and apply the blasting pressure time history curve in stages to the grid boundary of each corresponding borehole layout position according to the detonator delay sequence for integration solution to obtain the blasting damage field; The heat loss superposition module is used to apply the post-explosion ventilation cooling heat transfer boundary to the three-dimensional mesh model, calculate the mesh unit temperature at each cooling time step, calculate the thermal strain component according to the thermoelastic constitutive equation, calculate the thermal damage variable when the thermal strain component exceeds the limit, and nonlinearly superimpose the thermal damage variable and the blast damage variable of the blast damage field to obtain the total damage of each mesh unit. The depth retrieval module is used to divide the tunnel inner profile into multiple detection sections along the geometric axis. Taking the intersection of the geometric axis of the tunnel inner profile and the detection section as the origin of the polar coordinate system, it performs a radial search along the azimuth angle into the deep rock mass, retrieves the farthest grid cell that satisfies the condition that the total damage is greater than a preset critical damage threshold, and calculates the geometric distance of the farthest grid cell from the tunnel inner profile as the surrounding rock damage depth at that azimuth angle.

Images