Rock spiral tunnel group numerical simulation method, interaction system and equipment

By introducing the smooth Hawke-Brown criterion into the ABAQUS finite element software, a three-dimensional rock plastic damage model was constructed, which solved the limitations of ABAQUS in numerical simulation of rock spiral tunnel groups and achieved an accurate description of the complex stress state and plastic damage characteristics of rock tunnel groups.

CN121960012APending Publication Date: 2026-05-01HEBEI ROAD & BRIDGE GROUP +1
View PDF 0 Cites 0 Cited by

Patent Information

Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
HEBEI ROAD & BRIDGE GROUP
Filing Date
2025-12-24
Publication Date
2026-05-01

AI Technical Summary

Technical Problem

Existing finite element analysis software such as ABAQUS does not perform well in numerical simulation of rock spiral tunnel groups, and it is difficult to effectively describe the complex three-dimensional stress state and the plastic damage characteristics of the surrounding rock.

Method used

A three-dimensional rock plastic damage model based on the smooth Hawke-Brown criterion was adopted. The model was combined with linear, exponential and rational fractional softening functions. The plastic state was handled by a set of nonlinear stress integral equations. The UMAT subroutine in the ABAQUS finite element software was used to perform numerical simulation of the rock spiral tunnel group.

Benefits of technology

A three-dimensional numerical simulation of plastic damage in a rock spiral tunnel group was achieved, which made up for the limitations of ABAQUS software in the analysis of rock tunnel groups and improved the accuracy and precision of the simulation.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN121960012A_ABST
    Figure CN121960012A_ABST
Patent Text Reader

Abstract

The invention discloses a rock spiral tunnel group numerical simulation method, interaction system and equipment, and relates to the technical field of numerical simulation, the method comprises the following steps: obtaining an elastic stiffness matrix, an elastic-plasticity consistent tangent stiffness matrix and a first simulation parameter of the nth step of a target rock spiral tunnel group, based on the parameters and a three-dimensional rock plastic damage model constructed based on a smooth Hook-Brown criterion, determining a rock mass characteristic state of the target rock spiral tunnel group, for an elastic state, calculating a second simulation parameter in the (n + 1) th step based on a linear equation, and for a plastic state, calculating a second simulation parameter in the (n + 1) th step based on a linear equation; second simulation parameters of the (n + 1) th step are calculated on the basis of the nonlinear stress integral equation set, the simulation structure of the target rock spiral tunnel group is obtained on the basis of the second simulation parameters of all the steps, and the limitation of finite element analysis software on numerical simulation of the rock spiral tunnel group can be made up.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This application relates to the field of numerical simulation technology, and in particular to a numerical simulation method, interactive system and equipment for rock spiral tunnel groups. Background Technology

[0002] A spiral tunnel complex is a group of tunnels with small radii and large elevation differences stacked together. This unique design significantly affects the stress state of the surrounding soil and rock mass. Because the spiral tunnels traverse the mountain in a curved manner, the interaction between multiple spiral tunnels causes the soil and rock mass to no longer be in the simple two-dimensional or quasi-three-dimensional stress field of a traditional straight tunnel, but rather to form a complex three-dimensional stress state. This irregular stress distribution significantly increases the difficulty of excavation and support, posing higher requirements for engineering design and numerical simulation.

[0003] Existing analytical methods often fall short when analyzing complex spiral tunnel structures, especially spiral tunnel groups. Simulation analysis of spiral tunnel groups typically requires a combination of experimental and numerical simulation methods. ABAQUS software is a commonly used finite element method for rock mass analysis, but it is primarily used for numerical simulation of soil properties, and its performance on numerical simulations of spiral tunnel groups with specific rock properties is not ideal. Summary of the Invention

[0004] The purpose of this application is to provide a numerical simulation method, interactive system, and equipment for rock spiral tunnel groups, which can overcome the limitations of finite element analysis software in numerical simulation of rock spiral tunnel groups.

[0005] To achieve the above objectives, this application provides the following solution: Firstly, this application provides a numerical simulation method for rock spiral tunnel groups, including: S1. Obtain the elastic stiffness matrix, the elastoplastic uniform tangential stiffness matrix, and the first simulation parameters of the nth step for the target rock spiral tunnel group; the first simulation parameters of the nth step include the first nominal stress, effective stress, strain, plastic strain, equivalent plastic shear strain, and strain increment of the (n+1)th step; n is an integer greater than 0.

[0006] S2. When the rock mass is in an elastic state, the second simulation parameters for step (n+1) are calculated based on the elastic stiffness matrix, the first simulation parameters from step n, and the linear equations. When the rock mass is in a plastic state, the second simulation parameters for step (n+1) are calculated based on the elastic stiffness matrix, the elastoplastic uniform tangential stiffness matrix, the first simulation parameters from step n, and the nonlinear stress integral equations of the second model. The second simulation parameters for step (n+1) include the second nominal stress from step n+1 and the uniform linear stiffness matrix. The second model is a three-dimensional rock plastic damage model constructed based on the smoothed Hawke-Brown criterion. Here, the linear equations refer to linear stress-strain relationships.

[0007] S3. Determine whether the value of n+1 is equal to the preset first threshold. If yes, obtain the simulated structure of the target rock spiral tunnel group based on the second simulation parameters of all steps. Otherwise, update the value of n to n+1, obtain the updated value of n, and return to S1.

[0008] Secondly, this application provides an interactive system comprising a first module and a second module.

[0009] The first module is used to: determine the elastic stiffness matrix, the elastoplastic uniform tangential stiffness matrix, and the first simulation parameters of the target rock spiral tunnel group, and send the elastic stiffness matrix, the elastoplastic uniform tangential stiffness matrix, and the first simulation parameters of the nth step to the second module; the first simulation parameters of the nth step include the first nominal stress, effective stress, strain, plastic strain, equivalent plastic shear strain, and strain increment of the (n+1)th step.

[0010] The process for determining the elastic stiffness matrix, the elastoplastic uniform tangent stiffness matrix, and the first simulation parameter in step n of the target rock spiral tunnel group is as follows: When n=1, the engineering parameters of the target rock spiral tunnel group are obtained, and based on the engineering parameters and the first model, the elastic stiffness matrix, the elastoplastic uniform tangent stiffness matrix, and the first simulation parameters of the nth step are determined; the first model is a three-dimensional solid tunnel model obtained by modeling the target rock spiral tunnel group based on the finite element analysis algorithm.

[0011] When n > 1, the second simulation parameters of step n are obtained, and the first simulation parameters of step n are determined based on the second simulation parameters of step n and the first model; the second simulation parameters of step n include the second nominal stress and the uniform linear stiffness matrix of step n.

[0012] The second module is used to: execute the numerical simulation method for rock spiral tunnel groups described in any of the above-mentioned methods.

[0013] Thirdly, this application provides a computer device, including: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the steps of the numerical simulation method for rock spiral tunnel groups described above.

[0014] According to the specific embodiments provided in this application, this application has the following technical effects: This application provides a numerical simulation method, interactive system, and device for rock spiral tunnel groups. The method involves: S1. Obtaining the elastic stiffness matrix, elastoplastic uniform tangent stiffness matrix, and first simulation parameters for step n of the target rock spiral tunnel group; the first simulation parameters for step n include the first nominal stress, effective stress, strain, plastic strain, equivalent plastic shear strain, and strain increment for step n+1; where n is an integer greater than 0; S2. When the rock mass is in an elastic state, calculating the second simulation parameters for step n+1 based on the elastic stiffness matrix, the first simulation parameters for step n, and the linear equations; when the rock mass is in a plastic state, calculating the second simulation parameters for step n+1 based on the elastic stiffness matrix, the elastoplastic uniform tangent stiffness matrix, and the nonlinear stress integral equations of the second model. The second simulation parameters; the second simulation parameters of step n+1 include the second nominal stress and the uniform linear stiffness matrix of step n+1; the second model is a three-dimensional rock plastic damage model constructed based on the smooth Hawke-Brown criterion; S3. Determine whether the value of n+1 is equal to the preset first threshold. If yes, obtain the simulated structure of the target rock spiral tunnel group based on the second simulation parameters of all steps. Otherwise, update the value of n to n+1, obtain the updated value of n, and return to S1. This allows the three-dimensional rock plastic damage model constructed based on the smooth Hawke-Brown criterion to describe the three-dimensional strength and deformation characteristics of the surrounding rock in the spiral tunnel, and to process the rock mass in the plastic state through a set of nonlinear stress integral equations, effectively making up for the limitations of finite element analysis software in numerical simulation of rock spiral tunnel groups. Attached Figure Description

[0015] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0016] Figure 1 This is an application environment diagram of a numerical simulation method for a rock spiral tunnel group according to an embodiment of this application; Figure 2 A flowchart illustrating a numerical simulation method for a group of rock spiral tunnels provided in an embodiment of this application; Figure 3 Shape factor Schematic diagram of the effect on the intensity curve; Figure 4 This is a schematic diagram illustrating the evolution of a three-dimensional rock plastic damage model based on the smooth Hawke-Brown criterion as a function of the softening function. Figure 5 A schematic diagram of the stress update process of a numerical simulation method for a rock spiral tunnel group provided in an embodiment of this application; Figure 6 This is a schematic diagram of the structure of an interactive system provided in one embodiment of this application; Figure 7 This is a schematic diagram of the structure of a computer device provided in an embodiment of this application.

[0017] Figure reference numerals: 102 terminal, 104 server, 1 first module, 2 second module. Detailed Implementation

[0018] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0019] To make the above-mentioned objectives, features and advantages of this application more apparent and understandable, the application will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0020] The numerical simulation method for rock spiral tunnel groups provided in this application can be applied to, for example... Figure 1In the application environment shown, terminal 102 communicates with server 104 via a network. A data storage system can store the data that server 104 needs to process. The data storage system can be set up independently, integrated into server 104, or placed in the cloud or on other servers. Terminal 102 can send the elastic stiffness matrix, elastoplastic uniform tangent stiffness matrix, and the first simulation parameters of step n of the target rock spiral tunnel group to server 104. After receiving the first simulation parameters of step n of the target rock spiral tunnel group, server 104, based on the elastic stiffness matrix, the first simulation parameters of step n, and the second model of the target rock spiral tunnel group, determines the rock mass characteristic state of the target rock spiral tunnel group. When the rock mass characteristic state is elastic, server 104 calculates the second model of step n+1 based on the elastic stiffness matrix, the first simulation parameters of step n, and the linear equation. The simulated parameters are calculated as follows: when the rock mass is in a plastic state, the second simulation parameters for step n+1 are calculated based on the elastic stiffness matrix, the elastoplastic uniform tangent stiffness matrix, the first simulation parameters of step n, and the nonlinear stress integral equations of the second model. The second simulation parameters for step n+1 include the second nominal stress and the uniform linear stiffness matrix of step n+1. The second model is a three-dimensional rock plastic damage model constructed based on the smooth Hawke-Brown criterion. It is determined whether the value of n+1 is equal to a preset first threshold. If so, the simulated structure of the target rock spiral tunnel group is obtained based on the second simulation parameters of all steps. Otherwise, the value of n is updated to n+1, the updated value of n is obtained, and the step of obtaining the first simulation parameters of step n of the target rock spiral tunnel group is returned. The server 104 can feed back the obtained simulated structure or the second simulation parameters of step n+1 to the terminal 102. Furthermore, in some embodiments, the numerical simulation method for rock spiral tunnel groups can also be implemented separately by the server 104 or the terminal 102. For example, the terminal 102 can obtain the first simulation parameters of the nth step of the target rock spiral tunnel group to be processed and process the first simulation parameters of the nth step of the target rock spiral tunnel group to be processed. Alternatively, the server 104 can obtain the first simulation parameters of the nth step of the target rock spiral tunnel group to be processed from the data storage system and process the first simulation parameters of the nth step of the target rock spiral tunnel group to be processed.

[0021] The terminal 102 can be, but is not limited to, various desktop computers, laptops, smartphones, tablets, IoT devices, and portable wearable devices. IoT devices can include smart speakers, smart TVs, smart air conditioners, and smart in-vehicle devices. Portable wearable devices can include smartwatches, smart bracelets, and head-mounted devices. The server 104 can be implemented using a standalone server or a server cluster composed of multiple servers, or it can be a cloud server.

[0022] In one exemplary embodiment, such as Figure 2 As shown, a numerical simulation method for rock spiral tunnel groups is provided. This method is executed by computer equipment, specifically by a terminal or server alone, or by both a terminal and a server. In this embodiment, the method is applied to... Figure 1 Taking server 104 as an example, the following steps are included: S1. Obtain the elastic stiffness matrix, the elastoplastic uniform tangential stiffness matrix, and the first simulation parameters of the nth step for the target rock spiral tunnel group; the first simulation parameters of the nth step include the first nominal stress, effective stress, strain, plastic strain, equivalent plastic shear strain, and strain increment of the (n+1)th step; n is an integer greater than 0.

[0023] S2. Based on the elastic stiffness matrix, the first simulation parameters of step n, and the second model, determine the rock mass characteristic state of the target rock spiral tunnel group. When the rock mass characteristic state is elastic, calculate the second simulation parameters of step n+1 based on the elastic stiffness matrix, the first simulation parameters of step n, and the linear equation. When the rock mass characteristic state is plastic, calculate the second simulation parameters of step n+1 based on the elastic stiffness matrix, the elastoplastic uniform tangential stiffness matrix, the first simulation parameters of step n, and the nonlinear stress integral equations of the second model. The second simulation parameters of step n+1 include the second nominal stress and the uniform linear stiffness matrix of step n+1. The second model is a three-dimensional rock plastic damage model constructed based on the smooth Hawke-Brown criterion.

[0024] S3. Determine whether the value of n+1 is equal to the preset first threshold. If yes, obtain the simulated structure of the target rock spiral tunnel group based on the second simulation parameters of all steps. Otherwise, update the value of n to n+1, obtain the updated value of n, and return to the step of obtaining the first simulation parameters of the nth step of the target rock spiral tunnel group.

[0025] The effective stress and / or nominal stress in the first or second simulation parameters refer to the average stress acting on the undamaged configuration in actual engineering. The difference in names is only to distinguish different process stages and to serve as an intermediate variable for easy calculation.

[0026] The elastic stiffness matrix is ​​the inherent elastic stiffness matrix of the material. It is a fourth-order tensor and depends on the elastic modulus and Poisson's ratio of the target rock spiral tunnel group.

[0027] The elastoplastic uniform tangential stiffness matrix depends on the direction of plastic flow, hardening law, etc., and is only defined after plasticity is confirmed to have occurred.

[0028] Implementing the above steps enables a three-dimensional rock plastic damage model based on the smooth Hawke-Brown criterion to describe the three-dimensional strength and deformation characteristics of the surrounding rock in the spiral tunnel. Furthermore, the nonlinear stress integral equation system is used to process the rock mass in the plastic state, effectively compensating for the limitations of finite element analysis software in numerical simulation of rock spiral tunnel groups.

[0029] In another embodiment, the total number of incremental steps n (incremental step number) is determined by the analysis step setting in the finite element software ABAQUS. A computational model can contain multiple analysis steps, the number of which is directly determined by the user, and the total time for each analysis step can also be set by the user.

[0030] Each analysis step contains multiple increment steps (i.e., calculations from n to n+1), and the total number of increment steps can be automatically determined by the program based on the convergence of the nonlinear solution. ABAQUS allows users to set the initial increment step size, maximum increment step size, and minimum increment step size to control computational stability. Therefore, the total time corresponding to the entire analysis step will contain several increment steps, each increment step corresponding to an n+1 update process in a subroutine.

[0031] The origin of the strain increment from n to n+1: In ABAQUS, during the finite element numerical simulation of a rock spiral tunnel group, the user controls the overall deformation or stress level of the structure (the structure in the first model) within an analysis step by applying displacement boundary conditions or load boundary conditions. This overall deformation level corresponds to the total strain change within that analysis step. Under displacement control, the local strain increment is calculated from the element strain caused by the nodal displacement of the current increment step; under force control, the local strain increment is calculated from the nodal displacement increment caused by external forces. The entire analysis step is divided into multiple increment steps. For each increment step, the finite element solver calculates the local strain increment (strain increment from n to n+1) at each integration point based on the nodal displacement of the current increment step, and uses this strain increment as input to the second model for stress prediction, plasticity correction, and updating of material state variables.

[0032] In another exemplary embodiment of this application, the process for confirming the rock mass characteristic state involved in S2 includes: Substituting the first simulation parameters of step n into the yield function of the second model, the yield degree of step n is obtained. The expression of the yield function is: ; ; in, f This represents the yield function, where the value of the yield function is the degree of yield. , , These are, respectively, effective deviatoric stress, effective mean stress, and Lod angle of stress; It is a shape function; The uniaxial compressive strength of intact rock; These are empirical parameters used to control the shape of the curve in the failure criterion, which is related to GSI (Geological Strength Index of Rock Mass). These are dimensionless empirical constants; For dimensionless parameters, To measure the cohesion or integrity of rock masses; For shape factor, The values ​​of effective deviatoric stress, effective average stress, and stress Lod angle are obtained based on elastic trial stress and equivalent plastic shear strain (this is existing technology and will not be elaborated here).

[0033] like Figures 3-4 As shown, Figures 3-4 Yield surface of a three-dimensional rock plastic damage model based on the smoothed Hawke-Brown criterion, which is provided as an embodiment of the numerical simulation method for a group of rock spiral tunnels in this application; Figure 3 Shape factor Schematic diagram of the effect on the intensity curve. Figure 4 This is a schematic diagram illustrating the evolution of a three-dimensional rock plastic damage model based on the smooth Hawke-Brown criterion as a function of the softening function. Figure 3 In the middle, when or At that time, the intensity curves on the deviated plane degenerate into Drucker-Prager circles and curvilinear triangles, respectively, the latter circumscribed by the six corner points of the original Hawke-Brown criterion.

[0034] when At that time, the rock mass characteristic state is determined to be elastic.

[0035] when At that time, the rock mass characteristic state was determined to be plastic.

[0036] in, These are preset parameters, representing the tolerance used to determine whether the rock mass has entered a plastic state; This represents the equivalent plastic shear strain at step n; This represents the elastic trial stress at step n+1. ; This represents the elastic stiffness matrix. This represents the first nominal stress at step n. This represents the strain increment at step n+1.

[0037] In another exemplary embodiment of this application, S2, when the rock mass characteristic state is elastic, calculates the second simulation parameters for the (n+1)th step based on the elastic stiffness matrix, the first simulation parameters of the nth step, and the linear equation, specifically including: When the rock mass is in an elastic state, calculate the first intermediate parameter in step n+1. The first intermediate parameter in step n+1 includes the effective stress, plastic strain, equivalent plastic shear strain, and strain in step n+1. The calculation formula for the first intermediate parameter in step n+1 is as follows: ; ; ; ; ; in, This represents the effective stress at step n+1. This represents the elastic test stress at the (n+1)th step; This represents the elastic stiffness matrix. This represents the first nominal stress at step n. This represents the strain increment at step n+1; This represents the plastic strain at step n. This represents the plastic strain at step n+1; This represents the equivalent plastic shear strain at step n; This represents the equivalent plastic shear strain at step n+1; This represents the strain at step n+1. This represents the strain at step n.

[0038] Based on the first intermediate parameters of step (n+1), the second simulation parameters of step (n+1) are obtained; the calculation formula for the second simulation parameters of step (n+1) is as follows: ; ; in, This represents the second nominal stress at step n+1. This represents the effective stress at step n+1. Represents the damage variable at step n; Let represent the uniform linear stiffness matrix at step n+1; This represents the strain at step n+1; This represents the elastic stiffness matrix.

[0039] The damage variable is obtained based on a preset softening function, and the calculation formula for the damage variable is: ; in, Indicates damage variables; This represents the value of a preset softening function, which is a correlation function between the softening modulus and the equivalent plastic shear strain. The softening modulus is a parameter determined based on the target rock spiral tunnel.

[0040] If the rock mass enters a plastic state, plastic correction is required. This application proposes to perform plastic correction by iteratively solving the nonlinear stress integral equation system.

[0041] In another exemplary embodiment of this application, when the rock mass characteristic state is plastic, the calculation of the second simulation parameters in step n+1 based on the elastic stiffness matrix, the elastic-plastic consistent tangent stiffness matrix, the first simulation parameters in step n, and the nonlinear stress integral equations of the second model in step S2 specifically includes: When the rock mass is in a plastic state, the nonlinear stress integral equations of the second model are obtained, and based on the nonlinear stress integral equations and the first simulation parameters of the nth step, the second intermediate parameters of the (n+1)th step are calculated; the second intermediate parameters of the (n+1)th step include: the effective stress, equivalent plastic shear strain, and plastic multiplier of the (n+1)th step; the expression of the nonlinear stress integral equations is: ; in, denoted as a nonlinear stress integral equation system; x represents an iteration point, which includes the effective stress, equivalent plastic shear strain, and plastic multiplier at step n+1; This represents the effective stress at step n+1, with its initial value being the elastic test stress at step n+1. This represents the effective stress at step n; This represents the elastic-plastic uniform tangent stiffness matrix; This represents the strain increment at step n+1; This represents the plastic multiplier at step n+1, with an initial value of 0; The plastic potential function is represented by g, the value of which is obtained based on the elastic trial stress and the equivalent plastic shear strain (this is existing technology and will not be elaborated here). This represents the value of the plastic potential function at step n+1; This represents the effective deviatoric stress at step n+1; This represents the equivalent plastic shear strain at step n+1, with its initial value being the equivalent plastic shear strain at step n. This represents the equivalent plastic shear strain at step n; f Represents the yield function. This represents the second nominal stress at step n+1.

[0042] Based on the second intermediate parameters obtained in step (n+1), the second simulation parameters for step (n+1) are obtained, and the calculation formula for the second simulation parameters in step (n+1) is as follows: ; ; ; ; in, This represents the strain at step n+1. Indicates the strain at step n. This represents the strain increment at step n+1; This represents the plastic strain at step n. This represents the plastic strain at step n+1; This represents the second nominal stress at step n+1. This represents the effective stress at step n+1. This represents the damage variable at step n+1; Let represent the uniform linear stiffness matrix at step n+1.

[0043] The damage variable is obtained based on a preset softening function, and the calculation formula is: ; in, Indicates damage variables; This represents the value of a preset softening function, which is a correlation function between the softening modulus and the equivalent plastic shear strain. The softening modulus is a parameter determined based on the target rock spiral tunnel.

[0044] In another exemplary embodiment of this application, the softening function involved in the above embodiments is a rational number, exponential, or linear function, and the expression of the softening function is: ; From top to bottom, these correspond to rational, exponential, or linear function expressions, respectively. Represents the softening function. Indicates the softening modulus. It represents the equivalent plastic shear strain.

[0045] For example, when calculating the damage variable in step n, any softening function can be selected to represent the equivalent plastic shear strain in step n. Substitute the softening function above (in the alternative formula) The value of the softening function in step n is obtained, and then the damage variable in step n is calculated. When calculating the damage variable in step n+1, any softening function is selected, and the equivalent plastic shear strain in step n+1 is... Substitute the softening function above (in the alternative formula) The value of the softening function at step n+1 is obtained, and then the damage variable at step n+1 is calculated.

[0046] In another exemplary embodiment of this application, the calculation of the second intermediate parameters for the (n+1)th step based on the nonlinear stress integral equations and the first simulation parameters of the nth step, as described in the above embodiments, specifically includes: S100. Construct an initial iteration point and determine it as the first iteration point, wherein the initial iteration point is represented as: ; in, This represents the initial iteration point corresponding to solving the nonlinear stress integral equations in step n+1. This represents the effective stress at the initial iteration point, with an initial value of [value missing]. , , This represents the elastic trial stress at step n+1. This represents the elastic stiffness matrix. This represents the first nominal stress at step n. This represents the strain increment at step n+1; Let represent the equivalent plastic shear strain at the initial iteration point, whose initial value is the equivalent plastic shear strain at the nth step. The plastic multiplier at the initial iteration point is denoted by 0, and its initial value is 0.

[0047] S200. Based on the line search method and the first iteration point, obtain the search step size, update the first iteration point based on the search step size, and then determine the updated first iteration point as the second iteration point; the search step size includes the search size and the search direction; the update formula for the first iteration point is: ; in, This represents the second iteration point obtained by solving the nonlinear stress integral equations in the (n+1)th step, where k+1 represents the (k+1)th iteration. This represents the first iteration point. k Indicates the first k The next iteration; Indicates the search direction. Indicates the search size; jThe number of optimizations represents the number of times the line search method has been executed.

[0048] S300. Update the first iteration point to the second iteration point to obtain the updated first iteration point.

[0049] S400. When the first iteration point does not satisfy... When that happens, return to S200; This is the first iteration point.

[0050] When the first iteration point satisfies At that time, the second intermediate parameter for the (n+1)th step is obtained based on the first iteration point.

[0051] In another exemplary embodiment of this application, the search step size obtained in S200 of the above embodiment based on the line search method and the first iteration point specifically includes: S201. Obtain the initial value of the search direction and determine it as the first search direction. Then, based on the first iteration point and the first search direction, obtain the first set of iteration parameters. The first set of iteration parameters includes: a first search size, a merit function corresponding to the first iteration point, and a merit function corresponding to the first search direction. The merit function corresponding to the first search direction is the merit function obtained by updating the merit function corresponding to the first iteration point according to the first search direction. The calculation formula for the first set of iteration parameters is as follows: ; ; ; in, This indicates the first search direction. The initial value is 1, and the superscript is... k Indicates the current iteration number, subscript j Indicates the number of times the search step size has been optimized. j The initial value is 0 in each iteration; This represents the first search size; The Jacobian matrix represents the nonlinear stress integral equation system corresponding to the first iteration point x in the (n+1)th step. This represents the set of nonlinear stress integral equations corresponding to the first iteration point x in the (n+1)th step; This represents the set of nonlinear stress integral equations corresponding to the first iteration point x in the nth step; This represents the merit function corresponding to the first search direction; This represents the merit function corresponding to the first iteration point; This represents the first iteration point x at step n+1, following the first search direction. The nonlinear stress integral equations updated with the first search size d.

[0052] S202. When the first search direction and the corresponding merit function satisfy the preset iteration condition, or the number of times the search step size is optimized. j When the preset second threshold is reached, the search step size is obtained based on the first search direction and the first search size, and the optimization ends.

[0053] S203. When the first search direction and the merit function corresponding to the first search direction do not satisfy the preset iteration conditions, and the number of times the search step size is optimized... j If the preset second threshold is not reached, the first search direction is optimized and updated to obtain the second search direction and the corresponding good value function.

[0054] Based on the second search direction and the corresponding merit function, the first search direction and the corresponding merit function are updated respectively to obtain the updated first search direction and the updated merit function.

[0055] Update the number of times the search step size optimization was performed. j The value is j +1.

[0056] Return to S202.

[0057] The method in this embodiment operates in each iteration of the calculation process, using a line search method to optimize the search step size and improve the convergence of the solution.

[0058] In another exemplary embodiment of this application, the preset iteration condition in the above embodiments is: ; in, It is a preset experience value.

[0059] The calculation formula for the second search direction is as follows: ; in, Indicates the second search direction; It is a preset experience value; The first search direction is indicated by the formula for calculating the second search direction, which is used to avoid slow convergence caused by an excessively small search step size.

[0060] The merit function for the second search direction is as follows: ; in, The merit function represents the second search direction; This represents the first iteration point x at step n+1, following the second search direction. The nonlinear stress integral equations updated with the first search size d.

[0061] In this embodiment, and For the parameters of the inaccurate line search algorithm, to avoid slow convergence due to excessively small search step size, the values ​​of the second search direction must meet the following requirements: ; In this embodiment, for the case where the rock mass enters a plastic state, plastic correction is performed. Newton's iterative method is used to solve the nonlinear stress integral equations, and an inaccurate line search method is used to optimize the iteration step size. During the iteration process, the stress point is gradually pulled back to the yield surface at the (n+1)th step. Figure 5 As shown.

[0062] In another exemplary embodiment of this application, the first simulation parameters for the nth step of obtaining the target rock spiral tunnel group involved in S1 include: When n=1, obtain the first simulation parameters of the nth step from the output of the first model.

[0063] When n > 1, the second simulation parameters of step n are passed to the first model, and the first simulation parameters of step n are obtained from the output of the first model.

[0064] The first model is a three-dimensional solid tunnel model obtained by modeling the target rock spiral tunnel group based on the finite element analysis algorithm.

[0065] To verify the effectiveness of the method in this embodiment, a rock mass modeling analysis was performed on the Hankou Tunnel project case. The specific process and data results are as follows: Technical terms: User Material, Mohr Coulomb Plasticity, job, and Depvar all refer to functional modules in finite element software.

[0066] Based on the project (Hankou Tunnel), the entrance end is a steep slope with an incline of approximately 70°, and the valley is generally V-shaped. The tunnel elevation ranges from 781.46 to 1340.22 meters, with a maximum relative elevation difference of approximately 560 meters. The left tunnel is 4457 meters long, and the right tunnel is 4366 meters long.

[0067] A double-helix tunnel model (a three-dimensional solid tunnel model, corresponding to the first model in the embodiment) was established in ABAQUS. The radius of curvature was 700m. The positive step ring excavation method was adopted, and concrete lining was arranged. The anchor bolt type was quincunx.

[0068] In the material input interface of the ABAQUS software, User Material is used to customize user material properties and assign interface properties to the rock mass. At the same time, the properties of the lining are defined using ABAQUS's built-in Mohr Coulomb Plasticity, while the anchor bolts are in an elastic state.

[0069] Set the model boundary conditions, gravity field, and ground loads. Mesh the rock mass model, concrete lining, and anchor bolts using C3D8 element type. Submit the user subfile (the subroutine corresponding to the method in this embodiment) in the job interface and submit the job for computation. Verify the effectiveness of the method in this embodiment through numerical simulation results, as follows: Table 1 Comparison of settlement results for the arch of the model with a curvature of 700m

[0070] The method in this embodiment successfully integrates a three-dimensional rock plastic damage model based on the smooth Hawke-Brown criterion into the ABAQUS finite element software, realizing the three-dimensional plastic damage numerical simulation and analysis of rock spiral tunnel groups under complex stress states, thus making up for the shortcomings of the built-in material model in rock mass analysis in ABAQUS.

[0071] This invention verifies the correctness and effectiveness of the methods described in the above embodiments through rock mass modeling and analysis of the Hankou Tunnel project, providing a reliable tool for practical engineering analysis and possessing significant scientific value and application significance.

[0072] This application also provides an interactive system (such as...). Figure 6 As shown in the figure, the interactive system includes a first module and a second module.

[0073] The first module is used to: determine the elastic stiffness matrix, the elastoplastic uniform tangential stiffness matrix, and the first simulation parameters of the target rock spiral tunnel group, and send the elastic stiffness matrix, the elastoplastic uniform tangential stiffness matrix, and the first simulation parameters of the nth step to the second module; the first simulation parameters of the nth step include the first nominal stress, effective stress, strain, plastic strain, equivalent plastic shear strain, and strain increment of the (n+1)th step.

[0074] The process for determining the elastic stiffness matrix, the elastoplastic uniform tangent stiffness matrix, and the first simulation parameter in step n of the target rock spiral tunnel group is as follows: When n=1, the engineering parameters of the target rock spiral tunnel group are obtained, and based on the engineering parameters and the first model, the elastic stiffness matrix, the elastoplastic uniform tangent stiffness matrix, and the first simulation parameters of the nth step are determined; the first model is a three-dimensional solid tunnel model obtained by modeling the target rock spiral tunnel group based on the finite element analysis algorithm.

[0075] When n > 1, the second simulation parameters of step n are obtained, and the first simulation parameters of step n are determined based on the second simulation parameters of step n and the first model; the second simulation parameters of step n include the second nominal stress and the uniform linear stiffness matrix of step n.

[0076] The second module is used to execute the numerical simulation method for rock spiral tunnel groups described in the above embodiments.

[0077] In another embodiment, the first model is configured as a finite element numerical model processing main program built based on the finite element software ABAQUS; the second model is configured as follows: the UMAT subroutine written in Fortran is embedded into the ABAQUS finite element numerical model processing main program. The second model returns the second simulation parameters of step n to the ABAQUS finite element numerical model processing main program through the variables STRESS and DDSDDE, and simultaneously stores the first simulation parameters of step n in the state variable STATEV for subsequent analysis.

[0078] The construction process of the first model includes: Based on the actual rock spiral tunnel group project, the geometric parameters of the tunnel and surrounding rock structure were extracted, and a 1:1 scale three-dimensional solid tunnel model was established in ABAQUS. For different spiral angle sections, corresponding tunnel structure modules were constructed, and the complete three-dimensional solid tunnel model of the rock spiral tunnel group was realized by segment splicing.

[0079] The first model is also used to calibrate model parameters based on experimental data of the mechanical properties of actual rock masses in the field. In the ABAQUS software's material input interface, the User Material parameter is used to customize user material properties; the model's material parameters include the elastic modulus. Poisson's ratio internal friction angle Expansion angle Cohesive strength Softening modulus Material compressive strength Rock disturbance coefficient Geological strength indicators At the same time, use Depvar to set the number of STATEVs.

[0080] Before the first module can function, initialization configuration is required, including: The analysis steps are set up, mainly including the initial analysis step, the geostress analysis step, and the static general analysis step. The analysis time is set according to the simulation needs. Preferably, the total analysis step time is set to 1s, the initial increment step is set to 1E-03s, the minimum increment step size is set to 1E-8s, and the maximum increment step is set to 0.01s. These parameters are suggested values ​​and users can adjust them according to specific circumstances. After the analysis steps are set up, the field variables and historical variables to be output are specified.

[0081] In the LOAD module of ABAQUS, the boundary conditions of the module are set, and the initial ground load and gravity field are applied. In the Interaction module, the step-ring excavation method is adopted, and the excavation sequence is defined in different analysis steps.

[0082] In ABAQUS Module-Mesh, the first model was meshed, and the meshing should balance computational accuracy and efficiency. To improve the accuracy of tunnel stress analysis, eight-node linear three-dimensional solid elements C3D8 were selected for all rock elements.

[0083] In the Job module, before submitting the task, select the subroutine corresponding to the second model, submit the model, and wait for the numerical analysis calculation to complete.

[0084] This invention addresses the numerical simulation problem of plastic damage in rock spiral tunnel groups. Based on the ABAQUS finite element software platform, it delves into a three-dimensional numerical simulation method for rock tunnels. A three-dimensional rock plastic damage model based on the smoothed Hawke-Brown criterion is employed, combined with linear, exponential, and rational fractional softening functions, to accurately describe the three-dimensional strength and deformation characteristics of the surrounding rock during spiral tunnel construction. This is achieved by introducing a line search-implicit return mapping process (corresponding to the substitution and solution of the yield function and based on...). This invention proposes a numerical simulation method for rock helical tunnel groups by efficiently solving the nonlinear stress integral equations (constitutive governing equations) and relying on the UMAT subroutine developed in Fortran. This invention effectively overcomes the limitations of the embedded material model in ABAQUS software for rock tunnel engineering analysis, significantly improving the computational efficiency and accuracy of three-dimensional plastic damage simulation in rock engineering.

[0085] In one exemplary embodiment, a computer device is provided, which may be a server or a terminal, and its internal structure diagram may be as follows. Figure 7As shown, the computer device includes a processor, memory, input / output (I / O) interfaces, and a communication interface. The processor, memory, and I / O interfaces are connected via a system bus, and the communication interface is also connected to the system bus via the I / O interfaces. The processor provides computational and control capabilities. The memory includes non-volatile storage media and internal memory. The non-volatile storage media stores the operating system, computer programs, and a database. The internal memory provides the environment for the operation of the operating system and computer programs stored in the non-volatile storage media. The database stores numerical simulation data of rock spiral tunnel groups. The I / O interfaces are used for information exchange between the processor and external devices. The communication interface is used for communication with external terminals via a network connection. When the computer program is executed by the processor, it implements a numerical simulation method for rock spiral tunnel groups.

[0086] Those skilled in the art will understand that Figure 7 The structures shown are merely block diagrams of some structures related to the present application and do not constitute a limitation on the computer device to which the present application is applied. Specific computer devices may include more or fewer components than shown in the figures, or combine certain components, or have different component arrangements. In an exemplary embodiment, a computer device is provided, including a memory and a processor. The memory stores a computer program, and the processor executes the computer program to implement the steps in the above-described method embodiments.

[0087] In one exemplary embodiment, a computer-readable storage medium is provided storing a computer program that, when executed by a processor, implements the steps in the above-described method embodiments.

[0088] In one exemplary embodiment, a computer program product is provided, including a computer program that, when executed by a processor, implements the steps in the above-described method embodiments.

[0089] It should be noted that the user information (including but not limited to user device information, user personal information, etc.) and data (including but not limited to data used for analysis, data stored, data displayed, etc.) involved in this application are all information and data authorized by the user or fully authorized by all parties, and the collection, use and processing of the relevant data must comply with relevant regulations.

[0090] Those skilled in the art will understand that all or part of the processes in the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium. When executed, the computer program can include the processes of the embodiments described above. Any references to memory, databases, or other media used in the embodiments provided in this application can include at least one of non-volatile and volatile memory. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, optical memory, high-density embedded non-volatile memory, resistive random access memory (ReRAM), magnetic random access memory (MRAM), ferroelectric random access memory (FRAM), phase change memory (PCM), graphene memory, etc. Volatile memory can include random access memory (RAM) or external cache memory, etc. By way of illustration and not limitation, RAM can take many forms, such as Static Random Access Memory (SRAM) or Dynamic Random Access Memory (DRAM).

[0091] The databases involved in the embodiments provided in this application may include at least one type of relational database and non-relational database. Non-relational databases may include, but are not limited to, blockchain-based distributed databases. The processors involved in the embodiments provided in this application may be general-purpose processors, central processing units, graphics processing units, digital signal processors, programmable logic devices, quantum computing-based data processing logic devices, etc., and are not limited to these.

[0092] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0093] This document uses specific examples to illustrate the principles and implementation methods of this application. The descriptions of the above embodiments are only for the purpose of helping to understand the methods and core ideas of this application. Furthermore, those skilled in the art will recognize that, based on the ideas of this application, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of this application.

Claims

1. A numerical simulation method for a group of spiral tunnels in rock, characterized in that, The numerical simulation method for the rock spiral tunnel group includes: S1. Obtain the elastic stiffness matrix, elastoplastic uniform tangential stiffness matrix, and the first simulation parameters of the target rock spiral tunnel group; the first simulation parameters of the nth step include the first nominal stress, effective stress, strain, plastic strain, equivalent plastic shear strain, and strain increment of the (n+1)th step; n is an integer greater than 0; S2. When the rock mass is in an elastic state, the second simulation parameters for step n+1 are calculated based on the elastic stiffness matrix, the first simulation parameters from step n, and the linear equations. When the rock mass is in a plastic state, the second simulation parameters for step n+1 are calculated based on the elastic stiffness matrix, the elastoplastic uniform tangent stiffness matrix, the first simulation parameters from step n, and the nonlinear stress integral equations of the second model. The second simulation parameters for step n+1 include the second nominal stress from step n+1 and the uniform linear stiffness matrix. The second model is a three-dimensional rock plastic damage model constructed based on the smooth Hawke-Brown criterion. S3. Determine whether the value of n+1 is equal to the preset first threshold. If yes, obtain the simulated structure of the target rock spiral tunnel group based on the second simulation parameters of all steps. Otherwise, update the value of n to n+1, obtain the updated value of n, and return to S1.

2. The numerical simulation method for rock spiral tunnel groups according to claim 1, characterized in that, In S2, the process of confirming the rock mass characteristic state includes: Substituting the first simulation parameters of step n into the yield function of the second model, the yield degree of step n is obtained. The expression of the yield function is: ; ; in, f This represents the yield function, where the value of the yield function is the degree of yield. , , These are, respectively, effective deviatoric stress, effective mean stress, and Lod angle of stress; It is a shape function; The uniaxial compressive strength of intact rock; These are empirical parameters; These are dimensionless empirical constants; For dimensionless parameters, ; For shape factor, The values ​​of effective deviatoric stress, effective average stress, and stress Lode angle are obtained based on elastic trial stress and equivalent plastic shear strain. when At that time, the rock mass property state is determined to be elastic; when At that time, the rock mass property state was determined to be plastic state; in, These are preset parameters, representing the tolerance used to determine whether the rock mass has entered a plastic state; This represents the equivalent plastic shear strain at step n; This represents the elastic trial stress at step n+1. ; This represents the elastic stiffness matrix. This represents the first nominal stress at step n. This represents the strain increment at step n+1.

3. The numerical simulation method for rock spiral tunnel groups according to claim 1, characterized in that, When the rock mass is in an elastic state, the second simulation parameters for the (n+1)th step are calculated based on the elastic stiffness matrix, the first simulation parameters of the nth step, and the linear equation. Specifically, this includes: When the rock mass is in an elastic state, calculate the first intermediate parameter in step n+1. The first intermediate parameter in step n+1 includes the effective stress, plastic strain, equivalent plastic shear strain, and strain in step n+1. The calculation formula for the first intermediate parameter in step n+1 is as follows: ; ; ; ; ; in, This represents the effective stress at step n+1. This represents the elastic test stress at the (n+1)th step; This represents the elastic stiffness matrix. This represents the first nominal stress at step n. This represents the strain increment at step n+1; This represents the plastic strain at step n. This represents the plastic strain at step n+1; This represents the equivalent plastic shear strain at step n; This represents the equivalent plastic shear strain at step n+1; This represents the strain at step n+1. This represents the strain at step n; Based on the first intermediate parameters of step (n+1), the second simulation parameters of step (n+1) are obtained; the calculation formula for the second simulation parameters of step (n+1) is as follows: ; ; in, This represents the second nominal stress at step n+1. This represents the effective stress at step n+1. Represents the damage variable at step n; Let represent the uniform linear stiffness matrix at step n+1; This represents the strain at step n+1; This represents the elastic stiffness matrix; The damage variable is obtained based on a preset softening function, and the calculation formula for the damage variable is: ; in, Indicates damage variables; This represents the value of a preset softening function, which is a correlation function between the softening modulus and the equivalent plastic shear strain. The softening modulus is a parameter determined based on the target rock spiral tunnel.

4. The numerical simulation method for rock spiral tunnel groups according to claim 1, characterized in that, When the rock mass is in a plastic state, the second simulation parameters for the (n+1)th step are calculated based on the elastic stiffness matrix, the elastic-plastic consistent tangent stiffness matrix, the first simulation parameters of the nth step, and the nonlinear stress integral equations of the second model. Specifically, this includes: When the rock mass is in a plastic state, the nonlinear stress integral equations of the second model are obtained, and based on the nonlinear stress integral equations and the first simulation parameters of the nth step, the second intermediate parameters of the (n+1)th step are calculated; the second intermediate parameters of the (n+1)th step include: the effective stress, equivalent plastic shear strain, and plastic multiplier of the (n+1)th step; the expression of the nonlinear stress integral equations is: ; in, denoted as a nonlinear stress integral equation system; x represents an iteration point, which includes the effective stress, equivalent plastic shear strain, and plastic multiplier at step n+1; This represents the effective stress at step n+1, with its initial value being the elastic test stress at step n+1. This represents the effective stress at step n; This represents the elastic-plastic uniform tangent stiffness matrix; This represents the strain increment at step n+1; This represents the plastic multiplier at step n+1, with an initial value of 0; This represents the plastic potential function, where the value of g is obtained based on the elastic trial stress and the equivalent plastic shear strain. This represents the value of the plastic potential function at step n+1; This represents the effective deviatoric stress at step n+1; This represents the equivalent plastic shear strain at step n+1, with its initial value being the equivalent plastic shear strain at step n. This represents the equivalent plastic shear strain at step n; f Represents the yield function. This represents the second nominal stress at step n+1; Based on the second intermediate parameters obtained in step (n+1), the second simulation parameters for step (n+1) are obtained, and the calculation formula for the second simulation parameters in step (n+1) is as follows: ; ; ; ; in, This represents the strain at step n+1. Indicates the strain at step n. This represents the strain increment at step n+1; This represents the plastic strain at step n. This represents the plastic strain at step n+1; This represents the second nominal stress at step n+1. This represents the effective stress at step n+1. This represents the damage variable at step (n+1). Let represent the uniform linear stiffness matrix at step n+1; The damage variable is obtained based on a preset softening function, and the calculation formula is: ; in, Indicates damage variables; This represents the value of a preset softening function, which is a correlation function between the softening modulus and the equivalent plastic shear strain. The softening modulus is a parameter determined based on the target rock spiral tunnel.

5. The numerical simulation method for rock spiral tunnel groups according to claim 4, characterized in that, Based on the nonlinear stress integral equations and the first simulation parameters of step n, the second intermediate parameters of step n+1 are calculated, specifically including: S100. Construct an initial iteration point and determine it as the first iteration point, wherein the initial iteration point is represented as: ; in, This represents the initial iteration point corresponding to solving the nonlinear stress integral equations in step n+1. This represents the effective stress at the initial iteration point, with an initial value of [value missing]. , , This represents the elastic trial stress at step n+1. This represents the elastic stiffness matrix. This represents the first nominal stress at step n. This represents the strain increment at step n+1; Let represent the equivalent plastic shear strain at the initial iteration point, whose initial value is the equivalent plastic shear strain at the nth step. The plastic multiplier representing the initial iteration point has an initial value of 0; S200. Based on the line search method and the first iteration point, obtain the search step size, update the first iteration point based on the search step size, and then determine the updated first iteration point as the second iteration point; the search step size includes the search size and the search direction; the update formula for the first iteration point is: ; in, This represents the second iteration point obtained by solving the nonlinear stress integral equations in the (n+1)th step, where k+1 represents the (k+1)th iteration. This represents the first iteration point. k Indicates the first k The next iteration; Indicates the search direction. Indicates the search size; j This indicates the number of optimization attempts, which is the number of times the line search method has been executed. S300. Update the first iteration point to the second iteration point to obtain the updated first iteration point; S400. When the first iteration point does not satisfy... When that happens, return to S200; This is the first iteration point; When the first iteration point satisfies At that time, the second intermediate parameter for the (n+1)th step is obtained based on the first iteration point.

6. The numerical simulation method for rock spiral tunnel groups according to claim 5, characterized in that, Based on the line search method and the first iteration point, the search step size is obtained, specifically including: S201. Obtain the initial value of the search direction and determine it as the first search direction. Then, based on the first iteration point and the first search direction, obtain the first set of iteration parameters. The first set of iteration parameters includes: a first search size, a merit function corresponding to the first iteration point, and a merit function corresponding to the first search direction. The merit function corresponding to the first search direction is the merit function obtained by updating the merit function corresponding to the first iteration point according to the first search direction. The calculation formula for the first set of iteration parameters is as follows: ; ; ; in, This indicates the first search direction. The initial value is 1, and the superscript is... k Indicates the current iteration number, subscript j Indicates the number of times the search step size has been optimized. j The initial value is 0 in each iteration; Both d and d represent the first search size, where d is abbreviation; The Jacobian matrix represents the nonlinear stress integral equation system corresponding to the first iteration point x in the (n+1)th step. The nonlinear stress integral equation system corresponding to the first iteration point x in the (n+1)th step; The nonlinear stress integral equations corresponding to the first iteration point x in the nth step are represented by the following: This represents the merit function corresponding to the first search direction; This represents the merit function corresponding to the first iteration point; This represents the first iteration point x at step n+1, following the first search direction. The nonlinear stress integral equations updated with the first search size d; S202. When the first search direction and the corresponding merit function satisfy the preset iteration condition, or the number of times the search step size is optimized. j When the preset second threshold is reached, the search step size is obtained based on the first search direction and the first search size, and the optimization ends. S203. When the first search direction and the merit function corresponding to the first search direction do not satisfy the preset iteration conditions, and the number of times the search step size is optimized... j If the preset second threshold is not reached, the first search direction is optimized and updated to obtain the second search direction and the corresponding good value function. Based on the second search direction and the merit function corresponding to the second search direction, the first search direction and the merit function corresponding to the first search direction are updated respectively to obtain the updated first search direction and the updated merit function corresponding to the first search direction. Update the number of times the search step size optimization was performed. j The value is j +1; Return to S202.

7. The numerical simulation method for rock spiral tunnel groups according to claim 6, characterized in that, The preset iteration condition is: ; in, It is a preset experience value; The calculation formula for the second search direction is as follows: ; in, Indicates the second search direction; It is a preset experience value; Indicates the first search direction; The merit function for the second search direction is as follows: ; in, The merit function represents the second search direction; This represents the first iteration point x at step n+1, following the second search direction. The nonlinear stress integral equations are updated with the first search size d.

8. The numerical simulation method for rock spiral tunnel groups according to claim 2, characterized in that, Obtain the first simulation parameters for the nth step of the target rock spiral tunnel group, including: When n=1, obtain the first simulation parameters of the nth step output by the first model; When n > 1, the second simulation parameters of step n are passed to the first model, and the first simulation parameters of step n are obtained from the output of the first model; The first model is a three-dimensional solid tunnel model obtained by modeling the target rock spiral tunnel group based on the finite element analysis algorithm.

9. An interactive system, characterized in that, The interactive system includes a first module and a second module; The first module is used to: determine the elastic stiffness matrix, the elastoplastic uniform tangential stiffness matrix, and the first simulation parameters for the target rock spiral tunnel group, and send the elastic stiffness matrix, the elastoplastic uniform tangential stiffness matrix, and the first simulation parameters for the nth step to the second module; the first simulation parameters for the nth step include the first nominal stress, effective stress, strain, plastic strain, equivalent plastic shear strain, and the strain increment for the (n+1)th step; The process for determining the elastic stiffness matrix, the elastoplastic uniform tangent stiffness matrix, and the first simulation parameter in step n of the target rock spiral tunnel group is as follows: When n=1, the engineering parameters of the target rock spiral tunnel group are obtained, and based on the engineering parameters and the first model, the elastic stiffness matrix, the elastoplastic uniform tangent stiffness matrix and the first simulation parameters of the target rock spiral tunnel group are determined. The first model is a three-dimensional solid tunnel model obtained by modeling the target rock spiral tunnel group based on the finite element analysis algorithm; When n > 1, the second simulation parameters of step n are obtained, and the first simulation parameters of step n are determined based on the second simulation parameters of step n and the first model; the second simulation parameters of step n include the second nominal stress and the uniform linear stiffness matrix of step n. The second module is used to: execute the numerical simulation method for rock spiral tunnel groups as described in any one of claims 1-8.

10. A computer device, comprising: A memory, a processor, and a computer program stored in the memory and capable of running on the processor, characterized in that the processor executes the computer program to implement the numerical simulation method for rock spiral tunnel groups according to any one of claims 1-8.