A three-dimensional finite element simulation method for excavation unloading of fault fracture zone surrounding rock roadway

Abstract

This invention provides a three-dimensional finite element simulation method for unloading tunnel excavation in fault fracture zone surrounding rock, belonging to the technical field of three-dimensional finite element simulation methods. The method includes constructing a three-dimensional finite element model of the tunnel entirely within the fault fracture zone; applying gravity load to the finite element model; performing geostress balance analysis to obtain the target stable load value; using the average reaction force reaching the target stable load value as the trigger condition for unloading simulation; calculating the radial unloading surface force on the mesh nodes; calculating the unloading damage factor; identifying the plastic zone of the finite element model; and further constructing an exponential softening model. By introducing the unloading damage factor and the exponential softening model, the evolution law of the cohesion and internal friction angle of the fracture zone surrounding rock gradually decreasing with the accumulation of plastic deformation during unloading can be quantitatively described. The unloading simulation is triggered when the average reaction force reaches the target stable load value, making the unloading initiation condition more explicit and specific.

Description

Technical Field

[0001] This invention relates to the field of three-dimensional finite element simulation methods, specifically a three-dimensional finite element simulation method for excavating and unloading roadways in fault fracture zones. Background Technology

[0002] During tunnel excavation, the stress release of the surrounding rock is a typical unloading process, fundamentally different from conventional loading failure. Traditional methods for analyzing the stability of surrounding rock have the following technical shortcomings: First, most methods treat the surrounding rock as an ideal elastoplastic material, ignoring the damage evolution law of the rock mass strength gradually deteriorating with the accumulation of plastic deformation during the unloading process; second, finite element numerical simulation and limit analysis theory are often used separately, making it difficult to directly convert numerical simulation results into engineering-usable safety evaluation indicators; third, for complex geological conditions containing faults, existing methods cannot quantitatively describe the impact of faults on the geometry of the failure surface and the stability of the surrounding rock; fourth, traditional methods lack systematic calibration methods for determining damage softening parameters and cannot perform self-calibration based on field monitoring data.

[0003] In the prior art, patent document CN119442780A discloses a method for simulating stress changes and unloading effects of slopes during excavation by establishing a high-precision three-dimensional finite element model and combining various material parameters and actual boundary conditions. However, this method does not establish a quantitative relationship between the unloading path and the dynamic decay of rock mass cohesion and internal friction angle by introducing unloading damage factors and exponential softening models. It also does not organically combine finite element simulation with the upper limit theory of limit analysis to obtain the distribution of plastic zones and geometric failure parameters through numerical simulation, and then calculate the theoretical ultimate support pressure of the surrounding rock. Furthermore, it does not combine in-situ stress balance analysis to obtain the initial radial stress, and does not conduct a graded assessment of the surrounding rock state through dynamic safety factors and provide corresponding measures. Therefore, there is an urgent need for a three-dimensional finite element simulation method for unloading in roadways excavated in fault fracture zones.

[0004] 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

[0005] The purpose of this invention is to provide a three-dimensional finite element simulation method for excavation and unloading of roadways in fault fracture zones, so as to solve the problems mentioned in the background art.

[0006] To achieve the above objectives, the present invention provides the following technical solution: A three-dimensional finite element simulation method for unloading excavation in a fault-fractured rock tunnel includes the following steps: S1: Construct a three-dimensional finite element model of the tunnel located entirely within the fault fracture zone, divide the surrounding rock area of ​​the fault fracture zone into grid cells, and independently refine the main structural surface within the fault fracture zone. Label each grid node in the surrounding rock area of ​​the fault fracture zone to be excavated, and select the grid nodes at the top of the tunnel as vertex nodes. S2: Apply gravity load to the finite element model, perform ground stress balance analysis, extract the initial radial stress of each grid node on the tunnel excavation boundary, obtain the uniaxial compressive strength of the surrounding rock of the fault fracture zone through indoor uniaxial compression test, take the uniaxial compressive strength of the surrounding rock of the fault fracture zone as the peak strength, further process the peak strength to obtain the target stable load value; S3: After the stress balance analysis is completed, the average reaction force of the top node of the roadway is calculated using finite element software. The average reaction force reaching the target stable load value is used as the trigger condition for unloading simulation. At the same time, the time length of the unloading process is set, and the radial unloading surface force on the grid node is calculated based on the initial radial stress. S4: Apply the calculated radial unloading surface force of the grid node as a dynamic boundary condition to the corresponding grid node on the roadway excavation boundary, and calculate the equivalent plastic strain generated by each grid unit of the surrounding rock of the fault fracture zone to be excavated. Calculate the unloading damage factor based on the time length of the unloading process. S5: Based on the unloading damage factor of the grid element of the surrounding rock of each fault fracture zone, the plastic zone of the finite element model is identified, and an exponential softening model is further constructed to calculate the cohesion and internal friction coefficient of the surrounding rock of the fault fracture zone generated during the unloading process; based on the cohesion and internal friction coefficient of the surrounding rock of the fault fracture zone, the theoretical ultimate support pressure of the surrounding rock is calculated through the upper limit theory of limit analysis, the safety factor is further calculated, and the safety factor is classified and evaluated.

[0007] Furthermore, a three-dimensional finite element model of the tunnel, entirely located within the fault fracture zone, is constructed. The specific steps are as follows: Based on geological exploration data, the spatial distribution range of fault fracture zones is determined. The attitude, spacing, trace length, and angle between the dip angle and the tunnel axis of the structural surfaces within the fault fracture zones are collected through on-site structural surface mapping. The attitude data of the structural surfaces are projected onto stereographic projection maps and grouped. The fractal dimension of each group of structural surfaces is calculated using the box counting method. The fractal dimension of each group of structural surfaces is divided by the average fractal dimension of all structural surface groups, multiplied by the average trace length of the structural surface group divided by the average trace length of all structural surface groups, and multiplied by the sine of the angle between the dip angle and the tunnel axis of the structural surface group to obtain the weighted influence coefficient of the structural surface group. After calculating the weighted influence coefficients of all structural surface groups, they are sorted from largest to smallest, and the two largest groups are selected as the main structural surfaces. A three-dimensional finite element model of the tunnel, entirely located within the fault fracture zone, was constructed. Tunnel parameters, fault fracture zone geometric parameters, and constitutive parameters of the surrounding rock material within the fault fracture zone were imported into the finite element model. The tunnel parameters include the tunnel diameter; the fault fracture zone geometric parameters include the width of the fault fracture zone, the dip angle of the principal structural planes, and the spacing; the constitutive parameters of the surrounding rock material within the fault fracture zone include the initial cohesion, initial internal friction angle, elastic modulus, Poisson's ratio, and dilatation angle of the infill material within the fault fracture zone. The surrounding rock area of ​​the fault fracture zone to be excavated was meshed, and the principal structural planes within the fault fracture zone were independently meshed. Each mesh node in the surrounding rock area to be excavated was calibrated, and the mesh nodes at the top of the tunnel were selected as vertex nodes.

[0008] Furthermore, the initial radial stress is selected and the unloading triggering condition is determined. The specific steps are as follows: Gravity loads corresponding to the mesh elements are applied to the finite element model, and ground stress balance analysis is performed to obtain the initial stress field of the finite element model under its own weight. The initial radial stress of each mesh node in the initial stress field on the boundary of the roadway to be excavated is extracted. The uniaxial compressive strength of the surrounding rock in the fault fracture zone was obtained through indoor uniaxial compression tests and calibrated as the peak strength. The peak strength value was retained as a preset multiple and marked as the target stable load. After the in-situ stress balance analysis was completed, the finite element software automatically calculated the vertical reaction force value of each grid node at the roadway roof position. The average value of these vertical reaction force values ​​was taken to obtain the average reaction force at the roadway roof position. When the average reaction force at the roadway roof position reached the value of the target stable load, the unloading simulation was triggered. The unloading simulation used the initial radial stress obtained from the in-situ stress balance as the initial condition.

[0009] Further, the radial unloading surface force on the grid nodes in the fault fracture zone surrounding rock area is calculated. The specific steps are as follows: After determining the unloading time based on the initial radial force, the radial surface force of each grid node is calculated according to the preset unloading rate based on the initial radial stress of each grid node in the surrounding rock area of ​​the fault fracture zone to be excavated. This allows the radial unloading surface force to gradually decrease over time during the unloading process. When the calculation result is less than 0, it is taken as 0 to indicate that the radial constraint of the node is completely released.

[0010] Further, the unloading damage factor of the mesh element is calculated, and the specific steps are as follows: The calculated radial unloading surface force of the grid node is applied as a time-varying dynamic boundary condition to the corresponding grid node on the roadway excavation boundary. The equivalent plastic strain of each grid element of the fault fracture zone surrounding rock to be excavated at the start of unloading and at the end of unloading is recorded by finite element software. The difference between the two is used as the unloading damage factor of the grid element of the fault fracture zone surrounding rock.

[0011] Furthermore, the cohesion and internal friction coefficient of the surrounding rock in the fault fracture zone generated by unloading are calculated. The specific steps are as follows: The initial cohesion, initial internal friction angle, cohesion damage softening coefficient, and internal friction angle damage softening coefficient of the surrounding rock of the fault fracture zone are obtained, and the initial cohesion, initial internal friction angle, cohesion damage softening coefficient, and internal friction angle damage softening coefficient are used as input parameters of the finite element model. The cohesive damage softening coefficient and the internal friction angle damage softening coefficient were obtained by conducting indoor triaxial compression tests on the surrounding rock of the fault fracture zone and fitting and calibrating the changes in the rock stress-strain curves obtained from the tests. Obtain the unloading damage factor corresponding to each fault fracture zone surrounding rock grid unit to be excavated; Based on the unloading damage factor and the cohesive damage softening coefficient corresponding to each grid cell, the initial cohesive force is subjected to exponential softening reduction to obtain the unloading cohesive force of each grid cell. Based on the unloading damage factor corresponding to each grid cell and the internal friction angle damage softening coefficient, the internal friction coefficient corresponding to the initial internal friction angle is subjected to exponential softening reduction processing, and the internal friction angle after unloading of each grid cell is determined according to the reduced internal friction coefficient.

[0012] Further, the equivalent cohesion and equivalent internal friction coefficient of the plastic zone are calculated, and the specific steps are as follows: The plastic zone in the finite element model is identified based on the unloading damage factor of the surrounding rock grid unit of each fault fracture zone. The grid unit with the unloading damage factor greater than zero is identified as the plastic zone grid unit, and the various plastic zone grid units are combined into a plastic zone grid unit set. The volume, cohesion, and internal friction angle of each mesh element in the plastic region are obtained respectively. The cohesion of each plastic zone grid element is multiplied by its corresponding volume to obtain multiple cohesion-volume multiplication values. The cohesive volume multiplication values ​​are summed to obtain a weighted sum of cohesive volumes. The volumes of each of the plastic region mesh cells are summed to obtain the total volume of the plastic region; The equivalent cohesive force of the plastic region is obtained by comparing the weighted sum of the cohesive force volume with the total volume of the plastic region. The internal friction angle of each of the plastic zone mesh elements is processed by tangent to obtain the internal friction coefficient of each of the plastic zone mesh elements; The internal friction coefficient and the corresponding volume of each plastic zone mesh element are multiplied to obtain multiple internal friction coefficient-volume multiplication values. The volume multiplication values ​​of each internal friction coefficient are summed to obtain a weighted sum of the internal friction coefficient volumes. The equivalent internal friction coefficient of the plastic region is obtained by comparing the volume weighted sum of the internal friction coefficients with the total volume of the plastic region.

[0013] Further, the angle between the modified potential rupture surface and the horizontal plane is calculated, specifically through the following steps: Multiple vertical monitoring lines were set above the tunnel roof, including a vertical line at the midpoint of the roof, a vertical line at the left 1 / 4 span of the roof, a vertical line at the right 1 / 4 span of the roof, and vertical lines at the left and right edges of the roof. On each vertical monitoring line, starting from the tunnel roof, the line extends vertically upwards to the ground surface. A grid node was selected at fixed intervals as a monitoring node, and the change of the unloading damage factor of each monitoring node over time was recorded. At the same time, horizontal monitoring lines were set at different heights above the tunnel roof to monitor the horizontal extension range of the plastic zone at each height. When all nodes on any vertical monitoring line from the top plate to the ground surface satisfy the unloading damage factor being greater than 0, and all nodes on a horizontal monitoring line at a certain height from left to right in a continuous area satisfy the unloading damage factor being greater than 0, it is determined that the plastic zone is completely penetrated. When the plastic zone is fully penetrated, a distribution cloud map of the plastic zone is generated based on the finite element model using finite element post-processing software. The interface between the plastic zone and the non-plastic zone is determined as the potential fracture surface. Geometric observation parameters are extracted, including the angle between the potential failure surface of the roof and the horizontal plane, the vertical distance from the midpoint of the roadway roof to the intersection line between the potential fracture surface and the ground surface, and the average width of the plastic zone on both sides of the roadway where the unloading damage factor is greater than the preset value. Based on the plastic distribution cloud map, multiple height positions are selected on both sides of the roadway for measurement. At each height, starting from the roadway wall, the measurement is taken horizontally into the surrounding rock until the unloading damage factor drops below the preset value. This horizontal distance is the width of the plastic strain concentration zone at that height. The average width of the width measured at all heights is obtained by averaging the widths at all heights. For the left and right sides of the roadway, the average widths of each side are calculated separately, and the larger value is taken as the final average width. Based on the plastic distribution cloud map, multiple height positions are selected on both sides of the roadway for measurement. At each height, starting from the roadway wall, the measurement is taken horizontally into the surrounding rock until the unloading damage factor drops below the preset value. This horizontal distance is the width of the plastic strain concentration zone at that height. The average width of the width measured at all heights is obtained by averaging the widths at all heights. For the left and right sides of the roadway, the average widths of each side are calculated separately, and the larger value is taken as the final average width. The main structural surfaces within the fault fracture zone in the finite element model are marked, and the spatial relationship between each of the main structural surfaces and the plastic zone is identified. Among the main structural surfaces that intersect with the plastic zone, the structural surface with the longest intersection length with the plastic zone is determined as the control structural surface. The length and inclination angle of the control structural surface are recorded to determine whether the control structural surface intersects with the plastic zone. When the control structure surface intersects with the plastic zone, the angle between the potential fracture surface of the roof and the horizontal plane is corrected based on the angle between the potential fracture surface of the roof and the horizontal plane, the influence coefficient of the fault fracture zone structure surface, the length of the control structure surface, the vertical distance from the midpoint of the roadway roof to the line of intersection between the potential fracture surface and the ground surface, and the dip angle of the control structure surface, to obtain the corrected angle between the potential fracture surface and the horizontal plane. When the control structural surface does not intersect with the plastic zone, the influence coefficient of the fault fracture zone structural surface is set to zero, and the angle between the potential fracture surface of the top plate and the horizontal plane is determined as the corrected angle between the potential fracture surface and the horizontal plane.

[0014] Furthermore, the corrected angle of the failure surface, the equivalent cohesion of the plastic zone, the equivalent internal friction angle of the plastic zone, the average width of the plastic zone with unloading damage factors on both sides of the roadway exceeding the preset value, and the vertical falling velocity of the roof failure block are obtained. The internal friction term is calculated based on the equivalent cohesion and the equivalent internal friction angle to obtain the corresponding internal friction correction value; Based on the vertical falling velocity of the roof failure block, the corrected angle of the failure surface, the average width of the plastic zone on both sides of the roadway where the unloading damage factor is greater than the preset value, and the vertical distance from the midpoint of the roadway roof to the intersection of the potential fracture surface and the ground surface, the geometric energy consumption term is calculated to obtain the corresponding geometric energy consumption value. The equivalent cohesive force, the internal friction correction value, the vertical falling velocity of the top plate failure block, and the geometric energy dissipation value are multiplied to obtain the internal energy dissipation rate of the plastic zone. Obtain the tunnel width and the density of rock per unit volume; The first pressure term is obtained by processing the ratio of the internal energy dissipation rate to the roadway width and the vertical falling velocity of the roof failure block; The second pressure term is obtained by calculating the self-weight influence term based on the unit volume rock mass weight, the tunnel width, the corrected failure surface angle, the vertical distance from the midpoint of the tunnel roof to the intersection of the potential fracture surface and the ground surface, and the average width of the plastic zone on both sides of the tunnel where the unloading damage factor is greater than the preset value. The difference between the first pressure term and the second pressure term is processed to obtain the ultimate support pressure of the surrounding rock.

[0015] Further, the safety factor is calculated as follows: Obtain the ultimate support pressure of the surrounding rock, the surrounding rock grade coefficient of the fault fracture zone, the rock weight per unit volume, and the tunnel burial depth; The surrounding rock load characterization value is obtained by multiplying the surrounding rock grade coefficient of the fault fracture zone, the rock weight per unit volume, and the tunnel burial depth. The difference between the ultimate support pressure of the surrounding rock and the characteristic value of the surrounding rock load is processed to obtain the support pressure margin value. The dynamic safety factor of the roadway is obtained by comparing the support pressure margin value with the ultimate support pressure of the surrounding rock. The roadway safety status is obtained by performing a graded assessment based on the dynamic safety factor. When the dynamic safety factor is greater than the preset first safety threshold, the roadway safety status is determined to be a safe status. When the dynamic safety factor is greater than the preset second safety threshold and less than or equal to the first safety threshold, the roadway safety state is determined to be a critical state, and temporary support is added. When the dynamic safety factor is less than or equal to the second safety threshold, the roadway safety status is determined to be dangerous, and operations in the roadway are immediately stopped.

[0016] Compared with the prior art, the beneficial effects of the present invention are: This scheme establishes a quantitative relationship between the unloading path and the dynamic decay of rock mass cohesion and internal friction angle by introducing unloading damage factors and an exponential softening model. This enables coupled analysis of unloading damage and strength degradation of the surrounding rock, overcoming the shortcomings of traditional methods that neglect unloading damage. By organically combining finite element simulation with the upper limit theory of limit analysis, the scheme obtains the distribution of the plastic zone and geometric failure parameters through numerical simulation, and then calculates the theoretical ultimate support pressure of the surrounding rock, thus improving the accuracy of the scheme. By combining in-situ stress balance analysis to obtain the initial radial stress, the scheme classifies and evaluates the state of the surrounding rock through dynamic safety factors and provides corresponding measures, providing a reliable theoretical basis for tunnel support design. 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 graph showing the relationship between the ultimate support pressure of the surrounding rock and the safety factor of the corresponding roadway. 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 The present invention provides a technical solution: A three-dimensional finite element simulation method for unloading excavation in a fault-fractured rock tunnel includes the following steps: S1: Construct a three-dimensional finite element model of the tunnel located entirely within the fault fracture zone, divide the surrounding rock area of ​​the fault fracture zone into grid cells, and independently refine the main structural surface within the fault fracture zone. Label each grid node in the surrounding rock area of ​​the fault fracture zone to be excavated, and select the grid nodes at the top of the tunnel as vertex nodes. Based on geological exploration data, the spatial distribution range of fault fracture zones is determined. The attitude, spacing, trace length, and angle between the dip angle and the tunnel axis of the structural surfaces within the fault fracture zones are collected through on-site structural surface mapping. The attitude data of the structural surfaces are projected onto stereographic projection maps and grouped. The fractal dimension of each group of structural surfaces is calculated using the box counting method. The fractal dimension of each group of structural surfaces is divided by the average fractal dimension of all structural surface groups, multiplied by the average trace length of the structural surface group divided by the average trace length of all structural surface groups, and multiplied by the sine of the angle between the dip angle and the tunnel axis of the structural surface group to obtain the weighted influence coefficient of the structural surface group. After calculating the weighted influence coefficients of all structural surface groups, they are sorted from largest to smallest, and the two largest groups are selected as the main structural surfaces. In the above process, attitude refers to the orientation and dip state of the structural surface in space, specifically including two parameters: dip direction and dip angle. The dip direction represents the horizontal angle between the dipping direction of the structural surface and the due north direction, and the dip angle represents the maximum acute angle between the structural surface and the horizontal plane. Spacing refers to the vertical distance between two adjacent sets of structural surfaces, and trace length refers to the length of the visible trace that the structural surface can be observed on the exposed rock surface or exploration profile. Attitude determines the orientation and dip degree of the structural surface, spacing characterizes the density of the structural surface, and trace length reflects the extent of the structural surface's extension. First, the attitude data of all structural surfaces collected on-site, including dip and dip angle data, are projected onto a stereographic projection map. Each structural surface corresponds to a pole or a great circle arc on the map. Based on the clustering distribution characteristics of these projection points on the stereographic projection map, structural surfaces with similar attitudes are grouped together by visual identification. For each group of structural surfaces, the fractal dimension of the group is calculated using box counting: specifically, the area where the structural surfaces of the group are distributed is covered with a series of square boxes of different sizes. Starting with the larger box size, the number of boxes containing trace segments of the structural surfaces of the group is counted, and then the box size is gradually reduced, and the number of boxes containing trace segments is counted again. As the box size gradually decreases, the number of boxes containing traces will increase accordingly according to a certain power law relationship. The box size and the corresponding number of boxes are linearly fitted in a double logarithmic coordinate system. The absolute value of the slope of the fitted line is the fractal dimension of the trace distribution of the set of structural surfaces. This fractal dimension quantitatively characterizes the degree of spatial filling and self-similarity complexity of the set of structural surfaces. The larger the dimension value, the denser and more irregular the distribution of structural surfaces. In this context, a trace is a visible line of intersection formed by the intersection of a structural surface and the planar region covered by box counting. Since a structural surface is a planar structure with a certain extended area in three-dimensional space, when this surface is cut by a plane, its intersection line appears as a line, which is called the trace of the structural surface. In box counting analysis, the complete shape of the structural surface in three-dimensional space is not directly measured, but the length, direction and distribution of this trace are measured. The length of the trace is usually approximated as the trace length of the structural surface. By projecting the measured structural surface attitude data onto a stereographic projection map and scientifically grouping them, and using the box counting method to calculate the fractal dimension of each group of structural surfaces, the complexity and self-similarity of the spatial distribution of structural surfaces within the fault fracture zone were quantitatively characterized. Based on this, a weighted influence coefficient was constructed by comprehensively considering the fractal dimension, average trace length, and sine value of the angle between the structural surface dip angle and the roadway axis of each group of structural surfaces. From numerous structural surfaces, the two most important groups that play a controlling role in roadway stability were objectively selected as the main structural surfaces.

[0021] A three-dimensional finite element model of the tunnel, entirely located within the fault fracture zone, was constructed. Tunnel parameters, fault fracture zone geometric parameters, and constitutive parameters of the surrounding rock material within the fault fracture zone were imported into the finite element model. The tunnel parameters include the tunnel diameter; the fault fracture zone geometric parameters include the width of the fault fracture zone, the dip angle of the principal structural planes, and the spacing; the constitutive parameters of the surrounding rock material within the fault fracture zone include the initial cohesion, initial internal friction angle, elastic modulus, Poisson's ratio, and dilatation angle of the infill material within the fault fracture zone. The surrounding rock area of ​​the fault fracture zone to be excavated was meshed, and the principal structural planes within the fault fracture zone were independently meshed. Each mesh node in the surrounding rock area to be excavated was calibrated, and the mesh nodes at the top of the tunnel were selected as vertex nodes.

[0022] In the aforementioned process, when constructing a three-dimensional finite element model of the tunnel entirely located within the fault fracture zone, the two sets of main structural surfaces were independently meshed. Geometric parameters such as the width of the fault fracture zone, the dip angle and spacing of the main structural surfaces, as well as material constitutive parameters such as the initial cohesion, initial internal friction angle, elastic modulus, Poisson's ratio, and dilatation angle of the infill material within the fault fracture zone were accurately imported. Simultaneously, each mesh node in the surrounding rock area to be excavated was calibrated, and mesh nodes at the top of the tunnel were selected as vertex nodes. This processing significantly improved the accuracy of the finite element model in representing the dominant role of the structural surfaces within the fault fracture zone. This allows subsequent in-situ stress balance, unloading simulation, plastic zone identification, and surrounding rock stability assessment to more realistically reflect the control effect of the main structural surfaces on the tunnel excavation response, avoiding simulation distortion caused by arbitrary selection of structural surfaces or coarse mesh generation.

[0023] S2: Apply gravity load to the finite element model, perform ground stress balance analysis, extract the initial radial stress of each grid node on the tunnel excavation boundary, obtain the uniaxial compressive strength of the surrounding rock of the fault fracture zone through indoor uniaxial compression test, take the uniaxial compressive strength of the surrounding rock of the fault fracture zone as the peak strength, further process the peak strength to obtain the target stable load value; The specific steps for selecting the initial radial stress and determining the unloading trigger condition are as follows: Gravity loads corresponding to the mesh elements are applied to the finite element model, and ground stress balance analysis is performed to obtain the initial stress field of the finite element model under its own weight. The initial radial stress of each mesh node in the initial stress field on the boundary of the roadway to be excavated is extracted. The uniaxial compressive strength of the surrounding rock in the fault fracture zone was obtained through indoor uniaxial compression tests and calibrated as the peak strength. The peak strength value was retained as a preset multiple and marked as the target stable load. After the in-situ stress balance analysis was completed, the finite element software automatically calculated the vertical reaction force value of each grid node at the roadway roof position. The average value of these vertical reaction force values ​​was taken to obtain the average reaction force at the roadway roof position. When the average reaction force at the roadway roof position reached the value of the target stable load, the unloading simulation was triggered. The unloading simulation used the initial radial stress obtained from the in-situ stress balance as the initial condition.

[0024] In the above process, the preset multiplier refers to multiplying the uniaxial compressive strength of the surrounding rock obtained from the indoor uniaxial compression test by a coefficient less than 1, based on the engineering safety level and surrounding rock conditions, to obtain the target stable load value used to trigger the unloading simulation. In this scheme, the preset multiplier is usually taken as three levels: 0.6, 0.75, and 0.9, corresponding to low, medium, and high load levels, respectively. Among them, 0.9 is the highest level, used for a safer unloading trigger condition, that is, when the average reaction force at the top of the roadway obtained from the local stress balance analysis reaches 90% of the uniaxial compressive strength of the surrounding rock, the unloading simulation is triggered. In-situ stress equilibrium analysis refers to the static calculation of the initial stress field of the model when it reaches equilibrium under gravity by applying only gravity load, i.e., the self-weight of the rock mass, without applying any excavation or external load in finite element simulation. For the mesh nodes at the top of the tunnel, the finite element software can obtain their vertical reaction force in the following way: During the solution process, the finite element software performs an equivalent transformation of the internal forces of the element containing the mesh node at the top of the tunnel into nodal forces, and obtains the nodal force required by each mesh node to maintain the equilibrium of the element. The vertical component of this nodal force is the vertical reaction force value of the node. Users can directly output the vertical reaction force of these vertex nodes by setting a node set, i.e., pre-selected vertex nodes, in the post-processing module, and then take the average value as the average reaction force at the top of the tunnel.

[0025] S3: After the stress balance analysis is completed, the average reaction force of the top node of the roadway is calculated using finite element software. The average reaction force reaching the target stable load value is used as the trigger condition for unloading simulation. At the same time, the time length of the unloading process is set, and the radial unloading surface force on the grid node is calculated based on the initial radial stress. The specific steps for calculating the radial unloading surface force on the grid nodes in the fault fracture zone surrounding rock area are as follows: After determining the unloading time based on the initial radial force, the radial surface force of each grid node is calculated according to the preset unloading rate based on the initial radial stress of each grid node in the surrounding rock area of ​​the fault fracture zone to be excavated. This allows the radial unloading surface force to gradually decrease over time during the unloading process. When the calculation result is less than 0, it is taken as 0 to indicate that the radial constraint of the node is completely released.

[0026] The formula used in the above process is: in, express At the moment Radial unloading surface force of each grid node; Indicates the unloading rate; Indicates the first Initial radial stress of each grid node; This indicates a positive operator, meaning that the value inside the parentheses is 0 when it is negative. Indicates the uninstallation completion time.

[0027] In the above process, radial unloading surface force From initial radial stress Initially, the unloading rate decreases to zero linearly over time, where the unloading rate is... Controlling the rate of stress release, the technical effect of this linear unloading method is that it can truly reflect the continuous evolution process of the surrounding rock stress gradually decreasing from the original ground stress state to zero after the tunnel excavation. This makes the unloading path on each grid node match its initial stress state, that is, the grid node with the larger initial stress has a larger total unloading amount. At the same time, through a unified unloading rate parameter, it can flexibly simulate the impact of different excavation speeds, such as rapid blasting excavation or slow mechanical excavation, on the stability of the surrounding rock, providing dynamic boundary conditions that conform to physical reality for the accurate calculation of subsequent unloading damage factors. Radial unloading surface force It is the dependent variable, representing the process of unloading. Time applied to the first The normal pressure on each grid node, the independent variable in the formula includes the initial radial stress. Unloading rate Current time Uninstallation completion time The initial radial stress determines the initial amplitude of the unloading surface force, the unloading rate controls how quickly the surface force decreases, and the current time and the unloading end time together determine the unloading process. The relationships between the variables and the dependent variable are as follows: the radial unloading surface force is positively correlated with the initial radial stress; the greater the initial radial stress, the greater the initial unloading surface force. The radial unloading surface force is negatively correlated with the unloading rate and unloading time; the surface force decreases linearly as the unloading time or unloading rate increases. (Positive operator is used.) Ensure that the calculated surface force value is automatically reset to zero when it is negative, to avoid unnatural tensile stress states.

[0028] S4: Apply the calculated radial unloading surface force of the grid node as a dynamic boundary condition to the corresponding grid node on the roadway excavation boundary, and calculate the equivalent plastic strain generated by each grid unit of the surrounding rock of the fault fracture zone to be excavated. Calculate the unloading damage factor based on the time length of the unloading process. The specific steps for calculating the unloading damage factor of the mesh element are as follows: The calculated radial unloading surface force of the grid node is applied as a time-varying dynamic boundary condition to the corresponding grid node on the roadway excavation boundary. The equivalent plastic strain of each grid element of the fault fracture zone surrounding rock to be excavated at the start of unloading and at the end of unloading is recorded by finite element software. The difference between the two is used as the unloading damage factor of the grid element of the fault fracture zone surrounding rock.

[0029] In the above process, by applying the time-varying radial unloading surface force as a dynamic boundary condition to the roadway excavation boundary, the actual physical process of the gradual release of surrounding rock stress during unloading can be realistically simulated. By recording the equivalent plastic strain of each grid element at the start and end of unloading and calculating the difference, the degree of plastic deformation accumulated by each element throughout the unloading process can be quantitatively obtained, i.e., the unloading damage factor. This damage factor directly reflects the degree of damage to the surrounding rock under unloading.

[0030] S5: Based on the unloading damage factor of the grid element of the surrounding rock of each fault fracture zone, the plastic zone of the finite element model is identified, and an exponential softening model is further constructed to calculate the cohesion and internal friction coefficient of the surrounding rock of the fault fracture zone generated during the unloading process; based on the cohesion and internal friction coefficient of the surrounding rock of the fault fracture zone, the theoretical ultimate support pressure of the surrounding rock is calculated through the upper limit theory of limit analysis, the safety factor is further calculated, and the safety factor is classified and evaluated.

[0031] The specific steps for calculating the cohesion and internal friction coefficient of the surrounding rock in the fault fracture zone generated by unloading are as follows: The initial cohesion, initial internal friction angle, cohesion damage softening coefficient, and internal friction angle damage softening coefficient of the surrounding rock of the fault fracture zone are obtained, and the initial cohesion, initial internal friction angle, cohesion damage softening coefficient, and internal friction angle damage softening coefficient are used as input parameters of the finite element model. The cohesive damage softening coefficient and the internal friction angle damage softening coefficient were obtained by conducting indoor triaxial compression tests on the surrounding rock of the fault fracture zone and fitting and calibrating the changes in the rock stress-strain curves obtained from the tests. Obtain the unloading damage factor corresponding to each fault fracture zone surrounding rock grid unit to be excavated; Based on the unloading damage factor and the cohesive damage softening coefficient corresponding to each grid cell, the initial cohesive force is subjected to exponential softening reduction to obtain the unloading cohesive force of each grid cell. Based on the unloading damage factor corresponding to each grid cell and the internal friction angle damage softening coefficient, the internal friction coefficient corresponding to the initial internal friction angle is subjected to exponential softening reduction processing, and the internal friction angle after unloading of each grid cell is determined according to the reduced internal friction coefficient.

[0032] The formula upon which the above process is based is: The initial cohesion of the surrounding rock Initial internal friction angle The cohesion damage softening coefficient of the surrounding rock was calibrated by indoor triaxial compression tests. and internal friction angle damage softening coefficient As input parameters to the finite element model, based on the unloading damage factor of the mesh elements in each surrounding rock area to be excavated, an exponential softening model is constructed to calculate the cohesion of the unloaded rock: in, The first [section] represents the surrounding rock area to be excavated. The cohesion of the rock in each grid cell; The cohesive damage softening coefficient is obtained by fitting the stress-strain curve of the rock through triaxial compression tests on the tunnel rock. The first [section] represents the surrounding rock area to be excavated. The unloading damage factor of each grid cell; This represents the initial cohesive force; Calculate the internal friction coefficient: in, The first [section] represents the surrounding rock area to be excavated. The internal friction angle of each grid cell; The internal friction angle damage softening coefficient is obtained by fitting the stress-strain curve of the rock through a triaxial compression test on the tunnel rock. This represents the initial internal friction angle.

[0033] In the above process, the dependent variable is the rock cohesion under the current unloading state. and internal friction angle These parameters reflect the current shear strength and cohesion of the surrounding rock after unloading damage. Using an exponential softening model The calculations are based on the post-peak strength decay observed in triaxial compression tests of rocks: after reaching peak strength, the rock's strength decreases with plastic deformation, i.e., the unloading damage factor. The cohesion of the rock decays exponentially rather than nonlinearly due to the accumulation of microcracks within the rock, because the initiation, propagation, and connection of microcracks is an accelerated process. The technical advantage of this formula lies in: calibrating the damage softening coefficient through indoor triaxial compression tests. It can accurately describe the different decay rates of cohesion during the damage process of different rock types, such as hard and brittle rocks and soft and plastic rocks, so that the current cohesion of each mesh element in the finite element model matches its accumulated plastic deformation. internal friction angle Using formula The calculation is based on the physical nature of rock friction characteristics: the internal friction angle reflects the interlocking and friction between rock particles. With the accumulation of plastic deformation and the propagation of microcracks, the contact area between particles decreases, and the degree of interlocking weakens, leading to a gradual decrease in the internal friction angle, and its tangent value... It also follows the law of exponential decay. The technical advantage of this formula lies in: by converting the change in the internal friction angle into an exponential decay of its tangent, it ensures... The physical plausibility of a monotonically decreasing trend with increasing damage factor is that it is always greater than 0, while also being supported by an independent damage softening coefficient. The decay rates of cohesion and internal friction angle can be controlled separately, enabling the model to more precisely describe the strength evolution characteristics of different types of rocks during the unloading damage process. The independent variable in the formula includes the initial cohesion. Initial internal friction angle Cohesion damage softening coefficient Internal friction angle damage softening coefficient and unloading damage factors The initial cohesion and initial internal friction angle serve as the baseline strength values ​​under undamaged conditions. The damage softening coefficient reflects the rate of rock strength decay with accumulated damage; it is calibrated through indoor triaxial compression tests. The unloading damage factor characterizes the degree of accumulated plastic deformation in the mesh element. The relationship between each variable and the dependent variable is as follows: current cohesion... With initial cohesion There is a positive correlation between the initial cohesion and the current cohesion; the greater the initial cohesion, the greater the current cohesion. The current cohesion is also positively correlated with the damage softening coefficient. and unloading damage factors Negative correlation or The larger the exponent, the greater the value of the exponent. The smaller the value, the more severe the current decrease in cohesion. Similarly, the tangent of the current internal friction angle... tangent of the initial internal friction angle It is positively correlated with the damage softening coefficient. and unloading damage factors It shows a negative correlation.

[0034] The specific steps for calculating the equivalent cohesion and equivalent internal friction coefficient in the plastic zone are as follows: The plastic zone in the finite element model is identified based on the unloading damage factor of the surrounding rock grid unit of each fault fracture zone. The grid unit with the unloading damage factor greater than zero is identified as the plastic zone grid unit, and the various plastic zone grid units are combined into a plastic zone grid unit set. The volume, cohesion, and internal friction angle of each mesh element in the plastic region are obtained respectively. The cohesion of each plastic zone grid element is multiplied by its corresponding volume to obtain multiple cohesion-volume multiplication values. The cohesive volume multiplication values ​​are summed to obtain a weighted sum of cohesive volumes. The volumes of each of the plastic region mesh cells are summed to obtain the total volume of the plastic region; The equivalent cohesive force of the plastic region is obtained by comparing the weighted sum of the cohesive force volume with the total volume of the plastic region. The internal friction angle of each of the plastic zone mesh elements is processed by tangent to obtain the internal friction coefficient of each of the plastic zone mesh elements; The internal friction coefficient and the corresponding volume of each plastic zone mesh element are multiplied to obtain multiple internal friction coefficient-volume multiplication values. The volume multiplication values ​​of each internal friction coefficient are summed to obtain a weighted sum of the internal friction coefficient volumes. The equivalent internal friction coefficient of the plastic region is obtained by comparing the volume weighted sum of the internal friction coefficients with the total volume of the plastic region.

[0035] The above process is based on the following formula: The equivalent cohesion and equivalent internal friction coefficient of all plastic zones are calculated based on the cohesion and internal friction coefficient of the rock in each grid cell of the surrounding rock. in, This represents the equivalent cohesive force across all plastic regions; This represents the equivalent internal friction angle for all plastic zones; Indicates the first The volume of a grid cell in a rock mass area to be excavated; This represents the set of mesh elements with an unloading damage factor greater than 0.

[0036] In the above process, the range of the plastic zone in the finite element model is accurately identified by determining the region where the unloading damage factor is greater than 0. The volume, current cohesion and current internal friction angle of each grid element in the plastic zone are then volume-weighted averaged. The non-uniformly damaged plastic zone is equivalent to a homogeneous material with uniform equivalent cohesion and equivalent internal friction angle. This allows for the use of reasonable equivalent strength parameters as input when calculating the ultimate support pressure of the surrounding rock based on the upper limit theory of limit analysis. This establishes a quantitative bridge between the damage of the micro-element and the analysis of macroscopic failure, significantly improving the theoretical rigor and engineering reliability of the calculation of the ultimate support pressure of the surrounding rock.

[0037] The specific steps for calculating the angle between the modified potential rupture surface and the horizontal plane are as follows: Multiple vertical monitoring lines were set above the tunnel roof, including a vertical line at the midpoint of the roof, a vertical line at the left 1 / 4 span of the roof, a vertical line at the right 1 / 4 span of the roof, and vertical lines at the left and right edges of the roof. On each vertical monitoring line, starting from the tunnel roof, the line extends vertically upwards to the ground surface. A grid node was selected at fixed intervals as a monitoring node, and the change of the unloading damage factor of each monitoring node over time was recorded. At the same time, horizontal monitoring lines were set at different heights above the tunnel roof to monitor the horizontal extension range of the plastic zone at each height. When all nodes on any vertical monitoring line from the top plate to the ground surface satisfy the unloading damage factor being greater than 0, and all nodes on a horizontal monitoring line at a certain height from left to right in a continuous area satisfy the unloading damage factor being greater than 0, it is determined that the plastic zone is completely penetrated. When the plastic zone is fully penetrated, a distribution cloud map of the plastic zone is generated based on the finite element model using finite element post-processing software. The interface between the plastic zone and the non-plastic zone is determined as the potential fracture surface. Geometric observation parameters are extracted, including the angle between the potential failure surface of the roof and the horizontal plane, the vertical distance from the midpoint of the roadway roof to the intersection line between the potential fracture surface and the ground surface, and the average width of the plastic zone on both sides of the roadway where the unloading damage factor is greater than the preset value. Based on the plastic distribution cloud map, multiple height positions are selected on both sides of the roadway for measurement. At each height, starting from the roadway wall, the measurement is taken horizontally into the surrounding rock until the unloading damage factor drops below the preset value. This horizontal distance is the width of the plastic strain concentration zone at that height. The average width of the width measured at all heights is obtained by averaging the widths at all heights. For the left and right sides of the roadway, the average widths of each side are calculated separately, and the larger value is taken as the final average width. In the above process, by setting up multiple vertical monitoring lines, including the midpoint of the top plate, the left and right 1 / 4 span points and the left and right edge points, and horizontal monitoring lines at different heights, the expansion pattern of the plastic zone in both vertical and horizontal directions can be comprehensively captured. This overcomes the technical defect that using only a single vertical line may miss asymmetric failure or local penetration. When all nodes from the top plate to the ground surface on any vertical monitoring line meet the unloading damage factor greater than 0, and all nodes in the continuous area from left to right on the horizontal monitoring line at a certain height meet the unloading damage factor greater than 0, the plastic zone is judged to be completely penetrated. This comprehensive judgment criterion ensures that the plastic zone of the top plate not only penetrates to the ground surface in the vertical direction, but also connects with the plastic zones of the two side slopes in the horizontal direction to form a continuous failure surface. After the plastic zone is penetrated, the interface between the plastic zone and the non-plastic zone is determined as the potential fracture surface by the distribution cloud map of the plastic zone. The main structural surfaces within the fault fracture zone in the finite element model are marked, and the spatial relationship between each of the main structural surfaces and the plastic zone is identified. Among the main structural surfaces that intersect with the plastic zone, the structural surface with the longest intersection length with the plastic zone is determined as the control structural surface. The length and inclination angle of the control structural surface are recorded to determine whether the control structural surface intersects with the plastic zone. When the control structure surface intersects with the plastic zone, the angle between the potential fracture surface of the roof and the horizontal plane is corrected based on the angle between the potential fracture surface of the roof and the horizontal plane, the influence coefficient of the fault fracture zone structure surface, the length of the control structure surface, the vertical distance from the midpoint of the roadway roof to the line of intersection between the potential fracture surface and the ground surface, and the dip angle of the control structure surface, to obtain the corrected angle between the potential fracture surface and the horizontal plane. When the control structural surface does not intersect with the plastic zone, the influence coefficient of the fault fracture zone structural surface is set to zero, and the angle between the potential fracture surface of the top plate and the horizontal plane is determined as the corrected angle between the potential fracture surface and the horizontal plane.

[0038] The formula upon which the above process is based is: in, Indicates the angle between the modified potential fracture surface and the horizontal plane; Indicates the angle between the potential failure surface of the top plate and the horizontal plane; This represents the fault influence coefficient, which is 0 when the fault does not intersect with the plastic zone. Indicates the length of the controlling fault; It represents the vertical distance from the midpoint of the tunnel roof to the line of intersection between the potential fracture surface and the ground surface; Indicates the dip angle of the controlling fault.

[0039] In the above process, by marking faults in the finite element model and identifying their spatial relationship with the plastic zone, the controlling fault with the longest intersection length with the plastic zone is identified, thereby quantitatively describing the influence of the fault on the dip angle of the potential fracture surface. When the controlling fault intersects with the plastic zone, it acts as a "traction" or "deflection" force, causing the original fracture surface to tend to extend along the fault direction, resulting in an increase in the angle between the fracture surface and the horizontal plane. The corrected angle is positively correlated with the fault length, the sine of the dip angle of the controlling fault, and the fault influence coefficient, and negatively correlated with the height of the failure range. That is, the longer the controlling fault, the steeper the fault, and the shallower the failure range, the more significant the deflection effect on the fracture surface. When the controlling fault does not intersect with the plastic zone, the correction term disappears, and the fracture surface angle remains unchanged. Through this correction method, the quantitative adjustment of the dip angle of the potential fracture surface under fault conditions is achieved. The dependent variable is the angle between the modified potential rupture surface and the horizontal plane. It reflects the actual dip angle of the potential failure surface of the roof after considering the influence of the fault. The independent variables in the formula include the uncorrected fracture surface angle. Fault influence coefficient Controlling fault length Height of the damaged top slab and fault dip angle The uncorrected fracture surface angle This reflects the fault surface tendency and fault length when there is no fault influence. and fault dip angle It reflects the geometric characteristics of the controlling fault, and the height of the top plate failure range. It reflects the degree of upward extension of the plastic zone, and the fault influence coefficient. The intensity of the combined influence of the control fault on the fracture surface angle is determined when the control fault does not intersect with the plastic zone. The correction item disappeared.

[0040] The relationship between each variable and the dependent variable is as follows: Corrected fracture surface angle Angle with the uncorrected Positively correlated The larger the value, the larger the corrected angle. With fault influence coefficient Controlling fault length The sine value of the dip angle of the controlling fault All are positively correlated; the larger these parameters are, the greater the correction magnitude. Height of the top plate damage range Negative correlation The larger the fault length and The smaller the ratio, the smaller the correction magnitude. This correction formula enables quantitative adjustment of the dip angle of the potential rupture surface under the influence of faults.

[0041] The specific steps for calculating the ultimate support pressure of the surrounding rock are as follows: The corrected angle of the failure surface, the equivalent cohesion of the plastic zone, the equivalent internal friction angle of the plastic zone, the average width of the plastic zone with unloading damage factors greater than the preset value on both sides of the roadway, and the vertical falling velocity of the roof failure block are obtained. The internal friction term is calculated based on the equivalent cohesion and the equivalent internal friction angle to obtain the corresponding internal friction correction value; Based on the vertical falling velocity of the roof failure block, the corrected angle of the failure surface, the average width of the plastic zone on both sides of the roadway where the unloading damage factor is greater than the preset value, and the vertical distance from the midpoint of the roadway roof to the intersection of the potential fracture surface and the ground surface, the geometric energy consumption term is calculated to obtain the corresponding geometric energy consumption value. The equivalent cohesive force, the internal friction correction value, the vertical falling velocity of the top plate failure block, and the geometric energy dissipation value are multiplied to obtain the internal energy dissipation rate of the plastic zone. Obtain the tunnel width and the density of rock per unit volume; The first pressure term is obtained by processing the ratio of the internal energy dissipation rate to the roadway width and the vertical falling velocity of the roof failure block; The second pressure term is obtained by calculating the self-weight influence term based on the unit volume rock mass weight, the tunnel width, the corrected failure surface angle, the vertical distance from the midpoint of the tunnel roof to the intersection of the potential fracture surface and the ground surface, and the average width of the plastic zone on both sides of the tunnel where the unloading damage factor is greater than the preset value. The formula upon which the above process is based is: in, This represents the internal energy dissipation rate of all plastic regions; This indicates the average width of the plastic zone on both sides of the roadway where the unloading damage factor is greater than the preset value; This indicates the vertical falling velocity of the top plate failure block; The dependent variable is the internal energy dissipation rate in the plastic region. It reflects the energy consumed per unit time by the surrounding rock of the roadway under ultimate failure state due to plastic deformation and friction. Its technical effect is to quantify the complex failure process in the plastic zone into a scalar energy rate index, providing a dissipation term in the energy balance equation for solving the ultimate support pressure of the surrounding rock based on the upper bound theorem. The independent variables in the formula include equivalent cohesion. , equivalent internal friction angle Vertical falling speed of the top plate failure block Height of the damaged top slab Corrected fracture surface angle and the average width of the plastic zone on both sides of the tunnel The equivalent cohesion and equivalent internal friction angle reflect the overall shear strength of the plastic zone. The falling velocity of the roof failure block is a kinematically permissible velocity field parameter. The height of the roof failure range and the corrected fracture surface angle together determine the geometry of the roof failure surface. The average width of the plastic zone on both sides of the roadway reflects the horizontal extension range of the side failure zone.

[0042] The relationship between each variable and the dependent variable is as follows: internal energy dissipation rate With equivalent cohesion The cosine value of the equivalent internal friction angle Top plate falling speed All of these parameters are positively correlated; the larger these parameters are, the more energy is dissipated. Height of the top plate damage range They are positively correlated; the larger the area of ​​destruction, the more energy is dissipated. Angle with the corrected fracture surface The relationship is determined by two factors, the first of which is... Follow Increase and decrease, the second term Follow As it increases, the overall relationship becomes non-monotonic. Average width of the plastic zone on both sides of the tunnel They are positively correlated; the larger the area of ​​destruction, the more energy is dissipated.

[0043] Based on the upper limit theory of limit analysis and the internal energy dissipation rate of the plastic zone, the ultimate support pressure of the surrounding rock can be obtained. in, Indicates the ultimate support pressure of the surrounding rock; The weight of rock per unit volume; Indicates the width of the alleyway.

[0044] In the above process, the dependent variable is the ultimate support pressure of the surrounding rock. It represents the minimum support pressure required to maintain the stability of the surrounding rock in a roadway. Its technical effect is to transform the internal energy dissipation rate into a support pressure index that can be directly used in engineering design through the energy balance equation of the upper bound theorem, providing a benchmark value for calculating the dynamic safety factor. The independent variable in the formula includes the internal energy dissipation rate. Lane width The falling speed of the top plate failure block Rock mass Height of the damaged top slab Uncorrected fracture surface angle and the average width of the plastic zone on both sides of the tunnel Among them, the internal energy dissipation rate represents the ability of the surrounding rock to resist damage, the rock weight reflects the self-weight load of the damaged block, and the height of the top plate failure range and the angle of the fracture surface determine the geometry and volume of the plastic zone. The relationship between each variable and the dependent variable is as follows: ultimate support pressure of surrounding rock With internal energy dissipation rate There is a positive correlation; the greater the internal energy dissipation rate, the stronger the resistance of the surrounding rock itself, and the greater the required support pressure. With the width of the alley The correlation is negative; the wider the tunnel, the less energy is dissipated per unit width, and the larger the area affected by the self-weight load. With the speed of the top plate falling Negative correlation The larger the first item The smaller; With rock weight There is a negative correlation: the heavier the rock mass, the greater its own weight load and the smaller the required support pressure, because its own weight itself provides part of the resistance. Height of the top plate damage range , fracture surface angle and the average width of the plastic zone This is related to the expression for the volume of the plastic zone in the second term.

[0045] The safety factor is calculated as follows: Obtain the ultimate support pressure of the surrounding rock, the surrounding rock grade coefficient of the fault fracture zone, the rock weight per unit volume, and the tunnel burial depth; The surrounding rock load characterization value is obtained by multiplying the surrounding rock grade coefficient of the fault fracture zone, the rock weight per unit volume, and the tunnel burial depth. The difference between the ultimate support pressure of the surrounding rock and the characteristic value of the surrounding rock load is processed to obtain the support pressure margin value. The dynamic safety factor of the roadway is obtained by comparing the support pressure margin value with the ultimate support pressure of the surrounding rock. The roadway safety status is obtained by performing a graded assessment based on the dynamic safety factor. When the dynamic safety factor is greater than the preset first safety threshold, the roadway safety status is determined to be a safe status. When the dynamic safety factor is greater than the preset second safety threshold and less than or equal to the first safety threshold, the roadway safety state is determined to be a critical state, and temporary support is added. When the dynamic safety factor is less than or equal to the second safety threshold, the roadway safety status is determined to be dangerous, and operations in the roadway are immediately stopped.

[0046] The formula used in the above process is: in, Indicates the safety factor of the tunnel; Indicates the surrounding rock grade coefficient; Indicates the depth of the tunnel; In the above process, the dependent variable is the dynamic safety factor of the roadway. It reflects the safety margin of the current support resistance relative to the ultimate support pressure of the surrounding rock. Its technical effect is to transform the complex results of surrounding rock mechanics analysis into an intuitive, dimensionless safety evaluation index, facilitating engineers to quickly determine the stability state of the surrounding rock and take appropriate measures. The independent variable in the formula includes the ultimate support pressure of the surrounding rock. Rock grade coefficient Rock mass and the depth of the tunnel The ultimate support pressure of the surrounding rock represents the theoretical minimum support pressure required to maintain the stability of the surrounding rock. This represents the support resistance that the current support structure can provide. The relationship between each variable and the dependent variable is as follows: safety factor Ultimate support pressure of surrounding rock Positively correlated The larger the molecule The larger the size, the higher the safety factor; With current support resistance There is a negative correlation; the greater the support resistance, the smaller the molecule and the lower the safety factor. Among them, the surrounding rock grade coefficient is a quantitative parameter determined by comprehensively evaluating the quality grade of the surrounding rock according to the national "Engineering Rock Mass Classification Standard" and relevant industry specifications. Specifically, the surrounding rock grade coefficient... The value of this value directly corresponds to the classification standard of surrounding rock from level I to level V. The higher the level of surrounding rock, the more complete and harder the rock mass is, and the stronger its self-stabilizing ability. The larger the value, the lower the surrounding rock grade, indicating that the rock mass is more fragmented and weaker, and has a poorer self-stabilizing ability. The smaller the value, the better. Among them, Class I surrounding rock is extremely hard rock with an intact rock mass. The highest value is typically between 0.9 and 1.0; Class II surrounding rock is hard rock with a relatively intact rock mass. The value ranges from 0.7 to 0.9; Class III surrounding rock is relatively soft rock with generally poor rock mass integrity. The value is between 0.5 and 0.7; Class IV surrounding rock is soft rock and the rock mass is fractured. The value ranges from 0.3 to 0.5; Class V surrounding rock is extremely soft rock and the rock mass is extremely fractured. The value is the lowest, usually between 0.1 and 0.3; In the above embodiments, 20 sets of data on the ultimate support pressure of the surrounding rock and the corresponding safety factor of the roadway are given to reflect the change of the roadway safety factor with the change of the ultimate support pressure of the surrounding rock, as shown in Table 1: Table 1: Relationship between ultimate support pressure of surrounding rock and safety factor of corresponding roadway In Table 1 above, the surrounding rock grade coefficient is set. Corresponding to Class III surrounding rock, rock mass weight The tunnel is buried deep You can see the safety factor. Ultimate support pressure of surrounding rock Positively correlated The larger the molecule The larger the size, the higher the safety factor.

[0047] Based on dynamic safety factor The values ​​are graded and evaluated: when The time is a safe state; when The situation is critical, requiring the installation of temporary supports; when The situation is dangerous; all work inside the tunnel must be stopped immediately.

[0048] In the above process, the classification threshold is determined based on the relative relationship between the ultimate support pressure of the surrounding rock and the current support resistance, as well as engineering practice experience. When the dynamic safety factor SF is greater than 0.3, it indicates that the current support resistance is greater than 70% of the ultimate support pressure of the surrounding rock, meaning that the actual support force provided is much greater than the theoretical minimum required value, and the surrounding rock has sufficient safety reserves, thus it is judged as a safe state. When SF is between 0 and 0.3, it indicates that the current support resistance is in the range of 70% to 100% of the ultimate support pressure of the surrounding rock, the safety margin is insufficient, and the surrounding rock is close to the ultimate equilibrium state. At this time, the roof displacement and plastic zone may accelerate, and if no measures are taken, it is very easy to turn into instability and failure, thus it is judged as a critical state and temporary support needs to be added. When SF is less than or equal to 0, it indicates that the current support resistance can no longer meet the basic requirements of the ultimate support pressure of the surrounding rock, meaning that the actual support force provided is less than or equal to the theoretical minimum required value, the surrounding rock load has exceeded the support capacity, and instability accidents such as roof collapse and spalling may occur at any time, so work must be stopped immediately and emergency reinforcement measures must be taken. The division of these three thresholds reflects the gradual evolution from sufficient safety reserves to insufficient safety margins and then to complete failure.

[0049] 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.

[0050] 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.

[0051] 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.

[0052] 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 three-dimensional finite element simulation method for excavation and unloading of tunnels in fault-fractured rock zones, characterized by the following steps: include: S1: Construct a three-dimensional finite element model of the tunnel located entirely within the fault fracture zone, divide the surrounding rock area of ​​the fault fracture zone into grid cells, and independently refine the main structural surface within the fault fracture zone. Label each grid node in the surrounding rock area of ​​the fault fracture zone to be excavated, and select the grid nodes at the top of the tunnel as vertex nodes. S2: Apply gravity load to the finite element model, perform ground stress balance analysis, extract the initial radial stress of each grid node on the tunnel excavation boundary, obtain the uniaxial compressive strength of the surrounding rock of the fault fracture zone through indoor uniaxial compression test, take the uniaxial compressive strength of the surrounding rock of the fault fracture zone as the peak strength, further process the peak strength to obtain the target stable load value; S3: After the stress balance analysis is completed, the average reaction force of the top node of the roadway is calculated using finite element software. The average reaction force reaching the target stable load value is used as the trigger condition for unloading simulation. At the same time, the time length of the unloading process is set, and the radial unloading surface force on the grid node is calculated based on the initial radial stress. S4: Apply the calculated radial unloading surface force of the grid node as a dynamic boundary condition to the corresponding grid node on the roadway excavation boundary, and calculate the equivalent plastic strain generated by each grid unit of the surrounding rock of the fault fracture zone to be excavated. Calculate the unloading damage factor based on the time length of the unloading process. S5: Based on the unloading damage factor of the grid element of the surrounding rock of each fault fracture zone, the plastic zone of the finite element model is identified, and an exponential softening model is further constructed to calculate the cohesion and internal friction coefficient of the surrounding rock of the fault fracture zone generated during the unloading process; based on the cohesion and internal friction coefficient of the surrounding rock of the fault fracture zone, the theoretical ultimate support pressure of the surrounding rock is calculated through the upper limit theory of limit analysis, the safety factor is further calculated, and the safety factor is classified and evaluated.

2. The three-dimensional finite element simulation method for excavation and unloading of tunnels in fault-fractured rock zones according to claim 1, characterized in that, The specific steps for constructing a three-dimensional finite element model of the tunnel located entirely within the fault fracture zone are as follows: Based on geological exploration data, the spatial distribution range of fault fracture zones is determined. The attitude, spacing, trace length, and angle between the dip angle and the tunnel axis of the structural surfaces within the fault fracture zones are collected through on-site structural surface mapping. The attitude data of the structural surfaces are projected onto stereographic projection maps and grouped. The fractal dimension of each group of structural surfaces is calculated using the box counting method. The fractal dimension of each group of structural surfaces is divided by the average fractal dimension of all structural surface groups, multiplied by the average trace length of the structural surface group divided by the average trace length of all structural surface groups, and multiplied by the sine of the angle between the dip angle and the tunnel axis of the structural surface group to obtain the weighted influence coefficient of the structural surface group. After calculating the weighted influence coefficients of all structural surface groups, they are sorted from largest to smallest, and the two largest groups are selected as the main structural surfaces. A three-dimensional finite element model of the tunnel, entirely located within the fault fracture zone, was constructed. Tunnel parameters, fault fracture zone geometric parameters, and constitutive parameters of the surrounding rock material within the fault fracture zone were imported into the finite element model. The tunnel parameters include the tunnel diameter; the fault fracture zone geometric parameters include the width of the fault fracture zone, the dip angle of the principal structural planes, and the spacing; the constitutive parameters of the surrounding rock material within the fault fracture zone include the initial cohesion, initial internal friction angle, elastic modulus, Poisson's ratio, and dilatation angle of the infill material within the fault fracture zone. The surrounding rock area of ​​the fault fracture zone to be excavated was meshed, and the principal structural planes within the fault fracture zone were independently meshed. Each mesh node in the surrounding rock area to be excavated was calibrated, and the mesh nodes at the top of the tunnel were selected as vertex nodes.

3. The three-dimensional finite element simulation method for excavation and unloading of tunnels in fault-fractured rock zones according to claim 2, characterized in that, The specific steps for selecting the initial radial stress and determining the unloading trigger condition are as follows: Gravity loads corresponding to the mesh elements are applied to the finite element model, and ground stress balance analysis is performed to obtain the initial stress field of the finite element model under its own weight. The initial radial stress of each mesh node in the initial stress field on the boundary of the roadway to be excavated is extracted. The uniaxial compressive strength of the surrounding rock in the fault fracture zone was obtained through indoor uniaxial compression tests and calibrated as the peak strength. The peak strength value was retained as a preset multiple and marked as the target stable load. After the in-situ stress balance analysis was completed, the finite element software automatically calculated the vertical reaction force value of each grid node at the roadway roof position. The average value of these vertical reaction force values ​​was taken to obtain the average reaction force at the roadway roof position. When the average reaction force at the roadway roof position reached the value of the target stable load, the unloading simulation was triggered. The unloading simulation used the initial radial stress obtained from the in-situ stress balance as the initial condition.

4. The three-dimensional finite element simulation method for excavation and unloading of tunnels in fault-fractured rock zones according to claim 1, characterized in that, The specific steps for calculating the radial unloading surface force on the grid nodes in the surrounding rock region of the fault fracture zone are as follows: After determining the unloading time based on the initial radial force, the radial surface force of each grid node is calculated according to the preset unloading rate based on the initial radial stress of each grid node in the surrounding rock area of ​​the fault fracture zone to be excavated. This allows the radial unloading surface force to gradually decrease over time during the unloading process. When the calculation result is less than 0, it is taken as 0 to indicate that the radial constraint of the node is completely released.

5. The three-dimensional finite element simulation method for excavation and unloading of roadways in fault-fractured zones according to claim 4, characterized in that, The specific steps for calculating the unloading damage factor of the mesh element are as follows: The calculated radial unloading surface force of the grid node is applied as a time-varying dynamic boundary condition to the corresponding grid node on the roadway excavation boundary. The equivalent plastic strain of each grid element of the fault fracture zone surrounding rock to be excavated at the start of unloading and at the end of unloading is recorded by finite element software. The difference between the two is used as the unloading damage factor of the grid element of the fault fracture zone surrounding rock.

6. The three-dimensional finite element simulation method for excavation and unloading of tunnels in fault-fractured rock zones according to claim 5, characterized in that, The specific steps for calculating the cohesion and internal friction coefficient of the surrounding rock in the fault fracture zone generated by unloading are as follows: The initial cohesion, initial internal friction angle, cohesion damage softening coefficient, and internal friction angle damage softening coefficient of the surrounding rock of the fault fracture zone are obtained, and the initial cohesion, initial internal friction angle, cohesion damage softening coefficient, and internal friction angle damage softening coefficient are used as input parameters of the finite element model. The cohesive damage softening coefficient and the internal friction angle damage softening coefficient were obtained by conducting indoor triaxial compression tests on the surrounding rock of the fault fracture zone and fitting and calibrating the changes in the rock stress-strain curves obtained from the tests. Obtain the unloading damage factor corresponding to each fault fracture zone surrounding rock grid unit to be excavated; Based on the unloading damage factor and the cohesive damage softening coefficient corresponding to each grid cell, the initial cohesive force is subjected to exponential softening reduction to obtain the unloading cohesive force of each grid cell. Based on the unloading damage factor corresponding to each grid cell and the internal friction angle damage softening coefficient, the internal friction coefficient corresponding to the initial internal friction angle is subjected to exponential softening reduction processing, and the internal friction angle after unloading of each grid cell is determined according to the reduced internal friction coefficient.

7. The three-dimensional finite element simulation method for excavation and unloading of tunnels in fault-fractured rock zones according to claim 6, characterized in that, The specific steps for calculating the equivalent cohesion and equivalent internal friction coefficient in the plastic zone are as follows: The plastic zone in the finite element model is identified based on the unloading damage factor of the surrounding rock grid unit of each fault fracture zone. The grid unit with the unloading damage factor greater than zero is identified as the plastic zone grid unit, and the various plastic zone grid units are combined into a plastic zone grid unit set. The volume, cohesion, and internal friction angle of each mesh element in the plastic region are obtained respectively. The cohesion of each plastic zone grid element is multiplied by its corresponding volume to obtain multiple cohesion-volume multiplication values. The cohesive volume multiplication values ​​are summed to obtain a weighted sum of cohesive volumes. The volumes of each of the plastic region mesh cells are summed to obtain the total volume of the plastic region; The equivalent cohesive force of the plastic region is obtained by comparing the weighted sum of the cohesive force volume with the total volume of the plastic region. The internal friction angle of each of the plastic zone mesh elements is processed by tangent to obtain the internal friction coefficient of each of the plastic zone mesh elements; The internal friction coefficient and the corresponding volume of each plastic zone mesh element are multiplied to obtain multiple internal friction coefficient-volume multiplication values. The volume multiplication values ​​of each internal friction coefficient are summed to obtain a weighted sum of the internal friction coefficient volumes. The equivalent internal friction coefficient of the plastic region is obtained by comparing the volume weighted sum of the internal friction coefficients with the total volume of the plastic region.

8. The three-dimensional finite element simulation method for excavation and unloading of tunnels in fault-fractured rock zones according to claim 1, characterized in that, The specific steps for calculating the angle between the modified potential rupture surface and the horizontal plane are as follows: Multiple vertical monitoring lines were set above the tunnel roof, including a vertical line at the midpoint of the roof, a vertical line at the left 1 / 4 span of the roof, a vertical line at the right 1 / 4 span of the roof, and vertical lines at the left and right edges of the roof. On each vertical monitoring line, starting from the tunnel roof, the line went vertically upwards to the ground surface. A grid node was selected at fixed intervals as a monitoring node, and the change of the unloading damage factor of each monitoring node over time was recorded. At the same time, horizontal monitoring lines were set at different heights above the tunnel roof to monitor the horizontal extension range of the plastic zone at each height. When all nodes on any vertical monitoring line from the top plate to the ground surface satisfy the unloading damage factor being greater than 0, and all nodes on a horizontal monitoring line at a certain height from left to right in a continuous area satisfy the unloading damage factor being greater than 0, it is determined that the plastic zone is completely penetrated. When the plastic zone is fully penetrated, a distribution cloud map of the plastic zone is generated based on the finite element model using finite element post-processing software. The interface between the plastic zone and the non-plastic zone is determined as the potential fracture surface. Geometric observation parameters are extracted, including the angle between the potential failure surface of the roof and the horizontal plane, the vertical distance from the midpoint of the roadway roof to the intersection line between the potential fracture surface and the ground surface, and the average width of the plastic zone on both sides of the roadway where the unloading damage factor is greater than the preset value. Based on the plastic distribution cloud map, multiple height positions are selected on both sides of the roadway for measurement. At each height, starting from the roadway wall, the measurement is taken horizontally into the surrounding rock until the unloading damage factor drops below the preset value. This horizontal distance is the width of the plastic strain concentration zone at that height. The average width of the width measured at all heights is obtained by averaging the widths at all heights. For the left and right sides of the roadway, the average widths of each side are calculated separately, and the larger value is taken as the final average width. The main structural surfaces within the fault fracture zone in the finite element model are marked, and the spatial relationship between each of the main structural surfaces and the plastic zone is identified. Among the main structural surfaces that intersect with the plastic zone, the structural surface with the longest intersection length with the plastic zone is determined as the control structural surface. The length and inclination angle of the control structural surface are recorded to determine whether the control structural surface intersects with the plastic zone. When the control structure surface intersects with the plastic zone, the angle between the potential fracture surface of the roof and the horizontal plane is corrected based on the angle between the potential fracture surface of the roof and the horizontal plane, the influence coefficient of the fault fracture zone structure surface, the length of the control structure surface, the vertical distance from the midpoint of the roadway roof to the line of intersection between the potential fracture surface and the ground surface, and the dip angle of the control structure surface, to obtain the corrected angle between the potential fracture surface and the horizontal plane. When the control structural surface does not intersect with the plastic zone, the influence coefficient of the fault fracture zone structural surface is set to zero, and the angle between the potential fracture surface of the top plate and the horizontal plane is determined as the corrected angle between the potential fracture surface and the horizontal plane.

9. A three-dimensional finite element simulation method for excavation and unloading of tunnels in fault-fractured rock zones according to claim 8, characterized in that, The specific steps for calculating the ultimate support pressure of the surrounding rock are as follows: The corrected angle of the failure surface, the equivalent cohesion of the plastic zone, the equivalent internal friction angle of the plastic zone, the average width of the plastic zone with unloading damage factors greater than the preset value on both sides of the roadway, and the vertical falling velocity of the roof failure block are obtained. The internal friction term is calculated based on the equivalent cohesion and the equivalent internal friction angle to obtain the corresponding internal friction correction value; Based on the vertical falling velocity of the roof failure block, the corrected angle of the failure surface, the average width of the plastic zone on both sides of the roadway where the unloading damage factor is greater than the preset value, and the vertical distance from the midpoint of the roadway roof to the intersection of the potential fracture surface and the ground surface, the geometric energy consumption term is calculated to obtain the corresponding geometric energy consumption value. The equivalent cohesive force, the internal friction correction value, the vertical falling velocity of the top plate failure block, and the geometric energy dissipation value are multiplied to obtain the internal energy dissipation rate of the plastic zone. Obtain the tunnel width and the density of rock per unit volume; The first pressure term is obtained by processing the ratio of the internal energy dissipation rate to the roadway width and the vertical falling velocity of the roof failure block; The second pressure term is obtained by calculating the self-weight influence term based on the unit volume rock mass weight, the tunnel width, the corrected failure surface angle, the vertical distance from the midpoint of the tunnel roof to the intersection of the potential fracture surface and the ground surface, and the average width of the plastic zone on both sides of the tunnel where the unloading damage factor is greater than the preset value. The difference between the first pressure term and the second pressure term is processed to obtain the ultimate support pressure of the surrounding rock.

10. A three-dimensional finite element simulation method for excavation and unloading of tunnels in fault-fractured rock zones according to claim 3, characterized in that, The safety factor is calculated as follows: Obtain the ultimate support pressure of the surrounding rock, the surrounding rock grade coefficient of the fault fracture zone, the rock weight per unit volume, and the tunnel burial depth; The surrounding rock load characterization value is obtained by multiplying the surrounding rock grade coefficient of the fault fracture zone, the rock weight per unit volume, and the tunnel burial depth. The difference between the ultimate support pressure of the surrounding rock and the characteristic value of the surrounding rock load is processed to obtain the support pressure margin value. The dynamic safety factor of the roadway is obtained by comparing the support pressure margin value with the ultimate support pressure of the surrounding rock. The roadway safety status is obtained by performing a graded assessment based on the dynamic safety factor. When the dynamic safety factor is greater than the preset first safety threshold, the roadway safety status is determined to be a safe status. When the dynamic safety factor is greater than the preset second safety threshold and less than or equal to the first safety threshold, the roadway safety state is determined to be a critical state, and temporary support is added. When the dynamic safety factor is less than or equal to the second safety threshold, the roadway safety status is determined to be dangerous, and operations in the roadway are immediately stopped.

Images