Tunnel rockburst numerical simulation method based on node continuous-discontinuous deformation analysis

Abstract

The invention discloses a tunnel rockburst numerical simulation method based on node continuous-discontinuous deformation analysis, and belongs to the field of rockburst numerical simulation. The calculation area is divided into a destructible area and a non-destructible area, adaptive discretization processing is carried out, and joint units are formed; constructing a balance equation in a moment residual form, and performing linearization processing on the balance equation to obtain a displacement increment; adopting a bilinear constitutive model to unify normal and tangential mechanical responses of joint units into a full-quantity linear relation, and defining normal rigidity and tangential rigidity updating principles under different working conditions; constructing a rockburst simulation model based on the updated joint constitutive relationship and contact potential; and driving the rockburst simulation model to operate based on the displacement increment. According to the method, mechanical behavior transformation of rock mass from integrity to damage to block contact is truly reflected, stress waves are effectively absorbed, result distortion caused by boundary reflection waves is avoided, and a reliable stress environment is provided for rockburst simulation in cooperation with scientific initial stress generation.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of numerical simulation of rock burst in tunnel, and particularly to a numerical simulation method of rock burst in tunnel based on node continuous-discontinuous deformation analysis. BACKGROUND

[0002] Rock burst is a dynamic disaster caused by the instantaneous release of elastic strain energy and the violent ejection of rock blocks due to excavation unloading under high ground stress in deep tunnel. Its occurrence process is accompanied by stress redistribution, crack propagation and block movement in surrounding rock, and it has the characteristics of suddenness and strong destructiveness. Therefore, it is of great scientific and engineering significance to establish an efficient and accurate rock burst evaluation model to ensure the safety of deep tunnel.

[0003] Currently, there are four main methods for rock burst research: empirical criterion, field monitoring, laboratory test and numerical simulation. Empirical criterion uses strength-stress ratio, elastic energy index and brittleness coefficient as indicators, which is simple to operate. However, under complex deep geological conditions, a single criterion cannot balance universality and prediction accuracy. Field microseismic monitoring can track the evolution of rock failure and energy release in real time, and has been successfully applied to practical engineering. However, the sensor layout is limited, and signal interpretation relies on experience, which limits its effectiveness in early warning of rock burst. Laboratory tests such as true triaxial loading and acoustic emission monitoring tests can study the failure evolution and energy release mechanism of rock burst under controllable conditions, which helps to reveal its physical nature. However, laboratory tests have the following significant shortcomings: sample size is limited, unloading process is difficult to accurately simulate, control is complex, and repeatability is poor, making it difficult to carry out systematic parameter sensitivity analysis.

[0004] With the development of computing technology, numerical simulation has been widely used in rock burst prediction due to its low cost, repeatability and ease of implementation of complex conditions. The current mainstream methods include finite element method (FEM), finite difference method (FDM) and discrete element method (DEM), etc. However, each method has certain limitations in rock burst simulation: (1) Finite element method (FEM): based on the assumption of continuous medium, suitable for stress field analysis of surrounding rock, but has essential difficulties in simulating the non-continuous behaviors such as crack evolution, block separation and ejection during rock burst; (2) Finite difference method (FDM): has advantages in explicit time integration, suitable for large deformation and dynamic response analysis, but also based on the framework of continuous medium, making it difficult to naturally describe the non-continuous movement of rock mass after cracking; (3) Discrete element method (DEM): uses discrete block model, which can directly simulate the process of rock mass cracking, block movement and energy release, and has strong applicability in rock burst simulation. However, traditional DEM mostly uses explicit integration, and the time step is strictly limited by the critical value, resulting in low computational efficiency; the contact model between blocks is mostly based on simplified assumptions, making it difficult to accurately reflect the real nonlinear mechanical behavior of rock mass. SUMMARY

[0005] The application provides a tunnel rock burst numerical simulation method based on node continuous-discontinuous deformation analysis.

[0006] The application adopts the following technical scheme: a tunnel rock burst numerical simulation method based on node continuous-discontinuous deformation analysis, comprising the following steps: A calculation region is determined, the calculation region is divided into a breakable region and an unbreakable region according to engineering experience and rock mass damage possibility; the breakable region and the unbreakable region are discretized, and four-node cohesive elements are inserted on the common edges of all blocks in the breakable region to form joint elements; The internal cohesion and internal force are considered based on unbalanced force residual, and a balance equation in the form of residual is constructed; The balance equation is linearized to obtain displacement increment ; The normal mechanical response and the tangential mechanical response of the joint element are unified into a full linear relationship by adopting a bilinear constitutive model, and the updating principle of the normal stiffness and the tangential stiffness under different working conditions is determined; the joint elements in a damaged state are identified and deleted based on the updating principle, and contact modeling is performed between the adjacent blocks of the deleted joint elements to update the joint constitutive relationship and the contact potential; Based on the updated joint constitutive relationship and the contact potential, a static-dynamic unified boundary is adopted to convert the displacement boundary condition of the initial stress state into a viscous boundary condition, and a rock burst simulation model is constructed; The rock burst simulation model is driven by the displacement increment to realize accurate simulation of the non-continuous processes such as large deformation, damage and structure surface cracking of the rock mass involved in the tunnel rock burst, and output simulation result data.

[0007] In further embodiments, the discretization processing steps are as follows: A three-node grid element is selected as a discretization basic element, and grid adaptation design is adopted for the breakable region and the unbreakable region: the breakable region adopts dense three-node element grid elements, and the unbreakable region adopts loose three-node element grid elements; The displacement field of a single discretization basic element is expressed in the following form : ; in the formula, is a horizontal displacement component at a spatial position , a time , is a vertical displacement component at a spatial position , a time ; is a shape function matrix of the i th three-node triangular element, is the displacement unknown vector of the third triangular element.

[0008] In further embodiments, the non-equilibrium force residual is caused by inertial force, contact force and joint element and triangular element deformation; wherein the equilibrium equation in residual form at time instant is represented as follows: wherein, is the residual vector, is the time increment, is the next time instant to be solved, is the external force vector at time instant, is the inertial force vector at time instant, is the contact force vector at time instant, is the cohesive force vector at time instant, is the internal force vector at time instant.

[0009] In further embodiments, the construction procedure of the bilinear constitutive model is as follows: between adjacent discretized base elements, springs are inserted along the normal and tangential directions simultaneously, the bonding effect between adjacent discretized base elements is simulated by the mechanical response of the springs, and the mechanical characterization basis of the joint element is formed; Correspondingly, the working conditions at least include: loading working condition, unloading working condition and failure working condition; the unloading working condition includes: unloading sub-working condition and loading sub-working condition.

[0010] In further embodiments, based on the bilinear constitutive model, the full quantity linear relationship of the joint element in the normal direction is described as follows in the pre-peak loading segment , post-peak loading segment , post-peak unloading segment and failure segment : ; wherein, is the normal stress of the joint element, is the normal penalty parameter, is the normal opening of the joint element; wherein the calculation formula of the normal penalty parameter is as follows: ; wherein, and are the pre-peak and post-peak normal penalty parameters, and ​​​​​respectively represent the opening value when the normal peak tensile strength is reached and the maximum normal opening value before unloading.

[0011] In further embodiments, the full linear relationship of the joint element in the tangential direction based on the bilinear constitutive model is described as follows in the pre-peak loading segment , post-peak loading segment , post-peak unloading segment and failure segment . ; wherein, represents the shear stress of the joint element in the tangential direction, is the tangential penalty parameter, is the opening of the joint element in the tangential direction; wherein, the calculation formula of the tangential penalty parameter is as follows: ; wherein, and respectively represent the pre-peak and post-peak tangential penalty parameters, and respectively represent the slip value when the tangential peak shear strength is reached and the maximum tangential slip value before unloading, is the tangential residual shear stress.

[0012] In further embodiments, the identification method of the joint element in the failure state is as follows: real-time monitoring of the normal opening and the tangential opening of the integral points of the joint element, combining the bilinear constitutive model to give the failure threshold and establish the determination condition, the determination condition includes the normal determination condition and the tangential determination condition; if the joint element satisfies at least one of the normal determination condition or the tangential determination condition, it is determined that the joint element as a whole is in the failure state; wherein, the normal determination condition is: if both integral points have , it indicates that the normal stress exceeds the tensile strength of the rock mass , and it is determined that the corresponding joint element is normal failure; represents the normal opening value when the normal stress reaches the tensile strength of the rock mass ; the tangential determination condition is: if both integral points have , it indicates that the tangential shear stress has exceeded the shear strength of the rock mass , and it is determined that the corresponding joint element is tangential failure; represents the tangential residual slip when the tangential stress reaches the shear strength of the rock mass .

[0013] In a further embodiment, the contact potential is updated as follows: After deleting the joint units in a damaged state, identify adjacent blocks. and blocks Overlapping areas In overlapping areas The total normal contact force is obtained using the following formula. : ;in, Indicates the normal contact penalty parameter, and They are blocks and blocks Contact potential function, For gradient operators, and These are the areas of the infinitesimal elements. Relative to block and blocks Local coordinates, It is an overlapping region The area; Based on the total normal contact force The tangential contact force is updated using the following formula. : In the formula, express Friction at all times For tangential contact penalty parameters, The tangential relative velocity, For time increments, The length of the contact area; The tangential contact force is calculated. Make judgments and corrections if The tangential contact force is then corrected using the following formula. : In the formula, Let be the coefficient of friction of the contact surface. It is a norm.

[0014] In a further embodiment, the construction process of the rockburst simulation model is as follows: Based on the bilinear constitutive model, artificial damping is introduced to carry out dynamic analysis in order to generate the initial stress equilibrium state; By applying displacement boundary conditions using a spring, the initial stress equilibrium equation is constructed using the spring force equation, which is as follows: In the formula, and respectively the stiffness coefficient of the spring component and component, denotes the directional spring reaction force, denotes the directional spring reaction force; The spring that removes the edge of the calculation area for applying displacement conditions calculates the corresponding reaction force at the above-mentioned edge by using the initial stress equilibrium state equation; The pair of dampers are introduced to construct the viscous boundary by removing the calculation area edge of the spring, and the viscous boundary is expressed by the normal traction force and the tangential traction force Equation, specifically: , wherein, is the density of the rock mass, and is the absorption coefficient of the viscous boundary and is uniformly set, and are respectively the wave and wave speed in the rock mass, is the normal relative velocity, is the tangential relative velocity.

[0015] In further embodiments, the linearization expression of the balance equation is as follows: ; wherein, is the equivalent stiffness matrix, , is the elastic stiffness matrix, is the integral coefficient, is the time increment, is the mass matrix, is the joint element stiffness matrix; is the net force after deducting all resistances, which is the actual driving force of the rock mass to generate a new displacement increment ; ; and are respectively the velocity vector and the acceleration vector of all nodes in the calculation area at time; and are respectively the contact force vector and the joint element cohesion force vector at time.

[0016] The beneficial effects of the present application are: ​The application proposes a node continuous-discontinuous deformation analysis (NCDDA), which divides the calculation domain into an unbreakable area (continuous medium) and a breakable area (joint element connected breakable block), and simulates the progressive failure of surrounding rock through a bilinear constitutive joint element; a general contact model based on contact potential avoids complex searching; a viscous boundary and element failure technology realizes high ground stress instantaneous unloading and wave absorption. NCDDA can efficiently reproduce the whole process of rock burst of a circular tunnel, and can quantitatively evaluate the influence of strength parameters and structural planes on the timing and area of rock burst.

[0017] The application divides the calculation domain into an unbreakable area and a breakable area, and uses a dense three-node grid and a four-node cohesive joint element for the area which is significantly disturbed by tunnel excavation and easy to break, and combines a bilinear constitutive model to accurately describe the stress-strain relationship of the joint element in the normal and tangential directions, so that the whole process of rock burst from crack initiation (joint element unbroken) to crack propagation (joint element gradually broken) to block ejection (contact between adjacent blocks of broken joint elements) can be accurately restored. At the same time, by monitoring the normal opening and tangential opening in real time and combining with the damage threshold to determine the joint failure, the simulation distortion caused by the excessive simplification of the traditional method can be avoided, and the gradual process of the tunnel wall from local small damage area to overall rock burst can be accurately captured, which is consistent with the actual engineering rock burst evolution characteristics.

[0018] The linearization of the balance equation is , the equivalent stiffness matrix is formed by integrating the stiffness matrix K, the mass matrix M and the joint element stiffness matrix , and the displacement increment is quickly solved by combining the Gaussian elimination method and other mature numerical methods, avoiding the time-consuming problem of nonlinear equation iteration. In addition, the contact potential method calculates the normal contact force by overlapping area integration, without point-by-point edge searching of the contact relationship, further improving the calculation efficiency when the number of blocks increases, and meeting the demand for comparison analysis of multiple parameters such as internal friction angle and structural plane inclination angle in engineering. BRIEF DESCRIPTION OF DRAWINGS

[0019] Figure 1 is a flowchart of the numerical simulation method of tunnel rock burst.

[0020] Figure 2 is an identification diagram of the unbreakable area and the breakable area.

[0021] Figure 3 is a normal constitutive relationship diagram.

[0022] Figure 4 is a tangential constitutive relationship diagram.

[0023] Figure 5This is the definition diagram of the contact potential at point P in the triangular unit.

[0024] Figure 6 This is a contact model diagram based on the contact potential method. Detailed Implementation

[0025] The present invention will now be further described with reference to the accompanying drawings and embodiments.

[0026] Example 1 like Figure 1 As shown, this embodiment discloses a numerical simulation method for tunnel rockburst based on nodal continuous-discontinuous deformation analysis, including the following steps: The computational domain is determined, and based on engineering experience and the possibility of rock mass failure, the computational domain is divided into destructible and non-destructible regions. The destructible and non-destructible regions are discretized accordingly, and four-node cohesive elements are inserted on the common edges of all blocks in the destructible region to form joint elements. Based on the residual of non-equilibrium force, construct The equilibrium equations in the form of time residuals are obtained; the equilibrium equations are then linearized to obtain the displacement increments. ; A bilinear constitutive model is employed to unify the normal and tangential mechanical responses of joint elements into a fully linear relationship, clarifying the update principles for normal and tangential stiffness under different working conditions. Based on these update principles, joint elements in a failed state are identified and deleted. Contact modeling is then performed between adjacent blocks of the deleted joint elements to update the joint constitutive relationship and contact potential. It should be noted that the joint constitutive relationship mentioned here is a stress-strain relationship. Joint elements are a type of element specifically used to simulate weak surfaces such as joints, fissures, and faults in discontinuous media such as rock masses and concrete. Their core function is to describe the mechanical response (such as deformation and failure) of these weak surfaces under load, and the mathematical expression of this mechanical response is the joint element constitutive relationship. It determines how joint elements transmit forces and generate deformation, and is key to simulating the mechanical behavior of discontinuous media.

[0027] Based on the updated joint constitutive relation and contact potential, a rockburst simulation model is constructed by transforming the displacement boundary condition that generates the initial stress state into a viscous boundary condition using a unified static and dynamic boundary. Combined with the displacement increment Drive the rockburst simulation model to accurately simulate the discontinuous processes involved in tunnel rockburst, such as large deformation, failure, and structural surface rupture of the rock mass, and output simulation results data.

[0028] It should be noted that the purpose of dividing the calculation region into the breakable region and the unbreakable region in the present embodiment is as follows: rock mass damage is usually limited to a very small part of the entire calculation domain. For example, the excavation disturbed zone (EDZ) of a circular tunnel is usually estimated to be 1-3 times the radius of the tunnel. Therefore, in order to improve the calculation efficiency and ensure the calculation accuracy in the numerical simulation of rock mass damage in the tunneling process, only the concerned region where cracking and breaking may occur needs to be finely processed, and the other regions which remain intact during the entire simulation process need to be appropriately simplified. As shown in Figure 2 , the unbreakable region is assumed to remain intact, and a standard finite element discretization method can be used. However, the breakable region needs to be finely simulated to consider the initiation and propagation of cracks.

[0029] In order to improve the efficiency in processing numerical integration and contact modeling, the processing steps of the discretization processing in the present embodiment are as follows: As shown in Figure 2 , since the contact force is to be calculated using the contact potential, a three-node triangular element must be used in the two-dimensional case. Therefore, a three-node mesh element is selected as the discretization basic element, and a mesh adaptation design is adopted for the breakable region and the unbreakable region: the breakable region selects a dense three-node element mesh, and in the calculation, a non-continuous deformation analysis method is used, which is based on the rigid body-spring model, and the dense three-node element is regarded as a rigid body, and the elements are connected through joint elements (which can be regarded as springs with normal and tangential stiffness). By establishing the balance equation and the motion equation of each rigid body, considering the non-continuous behaviors such as contact, sliding and even separation of the joint elements, the interaction force between the elements can be calculated, which can accurately capture the non-continuous deformation phenomena such as cracking and sliding that may occur in the breakable region. The unbreakable region selects a loose three-node element mesh, and in the calculation, a finite element method is used, which is based on the variational principle or the weighted residual method, and the stiffness matrix of the element is established, and the node displacement is obtained by solving the linear equation set, and the deformation and stress of the material are described by using the constitutive relationship of the continuous medium.

[0030] The present embodiment uses the following form to express the displacement field of a single discretization basic element : ; in the formula, is the horizontal displacement component at the spatial position , time , is the vertical displacement component at the spatial position , time ; is the shape function matrix of the th three-node triangular element, is the displacement unknown vector of the th three-node triangular element.

[0031] In further embodiments, the shape function matrix is expressed as follows: wherein, is a 2x2 unit matrix, is the shape function corresponding to the th node, .

[0032] The displacement unknown vector is expressed as follows: , and is the horizontal displacement component and the vertical displacement component of the th node, .

[0033] Further, the unbalanced force residual is caused by the inertial force, the contact force and the deformation of the joint element and the triangular element; wherein the residual form of the balance equation at the moment is expressed as follows: wherein, is the residual vector, is the time increment, is the next moment to be solved, is the external force vector at the moment, represents the inertial force vector at the moment, is the contact force vector at the moment, is the cohesive force vector at the moment, represents the internal force vector at the moment.

[0034] Before the rock mass starts to be damaged, the independent DDA blocks are connected together through the joint elements embedded on the common block edges to simulate the continuous mechanical deformation. When the applied stress exceeds its tensile or shear strength, the joint element is regarded as being damaged and is deleted to simulate the rock mass rupture. Therefore, the mechanical response and the damage model of the joint element determine the accuracy of the continuous-discontinuous deformation analysis.

[0035] Therefore, the construction process of the bilinear constitutive model according to the embodiment is as follows: Between adjacent discretized basic elements, springs are inserted along the normal and tangential directions at the same time, the bonding effect between adjacent discretized basic elements is simulated by the mechanical response of the springs, and the mechanical characterization basis of the joint element is formed; Correspondingly, the working conditions at least include: a loading working condition, an unloading working condition and a failure working condition. The unloading working condition includes: an unloading sub-working condition and a loading sub-working condition.

[0036] Further, based on a bilinear constitutive model combined with Figure 3 , the full linear relationship of the joint element in the normal direction in the pre-peak loading segment , the post-peak loading segment , the post-peak unloading segment and the failure segment is described as follows: ; in the formula, denotes the normal stress of the joint element, is a normal penalty parameter, is the normal opening of the joint element; wherein the calculation formula of the normal penalty parameter is as follows: ; in the formula, and are the pre-peak and post-peak normal penalty parameters, and respectively denote the opening value when the normal peak tensile strength is reached and the maximum normal opening value reached before unloading.

[0037] Take Figure 4 for example, the full linear relationship of the joint element in the tangential direction in the pre-peak loading segment , the post-peak loading segment , the post-peak unloading segment and the failure segment is described as follows: ; in the formula, denotes the shear stress of the joint element in the tangential direction, is a tangential penalty parameter, is the tangential opening of the joint element; wherein the calculation formula of the tangential penalty parameter is as follows: ; in the formula, and respectively denote the pre-peak and post-peak tangential penalty parameters, and respectively denote the slip value when the tangential peak shear strength is reached and the maximum tangential slip value reached before unloading, is the tangential residual shear stress.

[0038] The identification method of the joint element in the failure state is as follows: real-time monitoring the normal opening and the tangential opening of the joint element , the bilinear constitutive model is combined to give the failure threshold and establish the criterion, which includes the normal criterion and the tangential criterion; If the joint element satisfies at least one of the normal criterion or the tangential criterion, it is determined that the joint element as a whole is in a failure state; Wherein, the normal criterion is: if both integral points have , it indicates that the normal stress exceeds the tensile strength of the rock mass , and it is determined that the corresponding joint element is normal failure; represents the normal opening value when the normal stress reaches the tensile strength of the rock mass ; The tangential criterion is: if both integral points have , it indicates that the tangential shear stress has exceeded the shear strength of the rock mass , and it is determined that the corresponding joint element is tangential failure; represents the tangential residual slip when the tangential stress reaches the shear strength of the rock mass .

[0039] Further, Figure 3 and Figure 4 represent the normal constitutive relationship and the tangential constitutive relationship, respectively. Among them, Figure 3 The calculation formula of the normal residual opening in is as follows: , wherein, is the peak normal stress related opening, is the energy release rate, and its calculation formula is: .

[0040] Correspondingly, Figure 4 the calculation formula of the tangential residual slip in is as follows: , wherein, is the peak shear stress related opening, is the shear energy release rate, and its calculation formula is: .

[0041] After the joint element is destroyed, it is deleted and contact modeling is performed between adjacent blocks, as shown in Figure 5 . The specific contact potential updating process is as follows: After deleting the joint element in a failure state, the overlapping area between the adjacent blocks and is identified; in the overlapping area ,The corresponding total normal contact force is obtained by using the following formula : ; wherein, represents the normal contact penalty parameter, and are the contact potential functions of the block and the block , is the gradient operator, and are the local coordinates of the micro-element area relative to the block and the block , is the area of the overlapping area .

[0042] Based on the total normal contact force , the tangential contact force is updated by using the following formula: ; wherein, represents the friction force at the moment, is the tangential contact penalty parameter, is the tangential relative velocity, is the time increment, is the length of the contact area; The calculated tangential contact force is judged and corrected. If , the tangential contact force is corrected by using the following formula: ; wherein, is the friction coefficient of the contact surface, is the norm.

[0043] Wherein, in order to prevent the influence of irregular unit shape on the smoothness of the contact force, the calculation of the contact potential function is based on the normal embedded displacement, instead of the traditional area-based calculation, as shown in Figure 6 . The calculation formula of the contact potential function is as follows: , wherein P is a point in the triangular element, , and are the areas of the triangles , and , is the area of the triangle , ​All are three-node triangular units.

[0044] Rockburst is a result of both high initial stress fields and sudden unloading. The key to rockburst modeling lies in correctly applying internal unloading conditions and non-reflection conditions at the outer boundaries of the computational domain. Typically, sudden unloading generates waves that propagate within the rock mass. If not handled properly, these waves can reflect off the boundaries of the computational domain, causing tensile stress waves and leading to unrealistic rock failure near those boundaries.

[0045] Therefore, the construction process of the rockburst simulation model in this embodiment is as follows: Based on the bilinear constitutive model, artificial damping is introduced to carry out dynamic analysis in order to generate the initial stress equilibrium state; By applying displacement boundary conditions using a spring, the initial stress equilibrium equation is constructed using the spring force equation, which is as follows: In the formula, and These are the stiffness coefficients of the springs. Components and Quantity, express Directional spring reaction force, express Directional spring reaction force.

[0046] Remove the springs at the edges of the computational region used to apply displacement conditions, and calculate the corresponding reaction forces at the aforementioned edges using the initial stress equilibrium equation. Paired dampers are introduced at the edge of the computational region where the spring is removed to construct a viscous boundary, which is subjected to a normal traction force. and tangential traction force The equation is expressed as follows: In the formula, For rock mass density, and The viscous boundary absorption coefficient is set uniformly. and They are respectively in the rock mass wave and wave speed, The normal relative velocity, The tangential relative velocity.

[0047] In another embodiment, the linearized expression of the equilibrium equation is as follows: In the formula, The equivalent stiffness matrix, , Here is the stiffness matrix. is the integral coefficient, is the time increment, is the mass matrix, is the joint element stiffness matrix; is the net force acting on the rock mass after all the resistances are subtracted from the external force; ; and are the velocity vector and acceleration vector of all the nodes in the calculation region at time and are the contact force vector and cohesion force vector at time .​​

Claims

1. A numerical simulation method for rock burst in tunnel based on node continuous-discontinuous deformation analysis, characterized in that, The method comprises the following steps: determining a calculation region, dividing the calculation region into a destructible region and a non-destructible region according to engineering experience and rock mass destruction possibility, discretizing the destructible region and the non-destructible region, and inserting four-node cohesive elements on common edges of all blocks in the destructible region to form joint elements; Based on the unbalanced force residual, construct the balance equation in the form of time residual; linearize the balance equation to obtain displacement increment ; unifying normal mechanical response and tangential mechanical response of the joint element method into a total linear relationship by using a bilinear constitutive model, and determining the updating principle of normal stiffness and tangential stiffness under different working conditions; identifying and deleting joint elements in a destruction state based on the updating principle, and updating the joint constitutive relationship and contact potential between adjacent blocks of the deleted joint elements through contact modeling; based on the updated joint constitutive relationship and contact potential, converting the displacement boundary condition generating an initial stress state into a viscous boundary condition by using a static-dynamic unified boundary to construct a rockburst simulation model; based on the displacement increment The driving rock burst simulation model runs, realizes the accurate simulation of the non-continuous processes such as large deformation, damage and structure surface fracture of rock mass involved in the rock burst of the tunnel, and outputs simulation result data.

2. The tunnel rockburst numerical simulation method based on the node continuous-discontinuous deformation analysis according to claim 1, characterized in that, the discretization processing comprises the following steps: selecting a three-node grid element as a discretization basic element, and designing grid adaptation for the destructible region and the non-destructible region: selecting a dense three-node element grid element for the destructible region, and selecting a loose three-node element grid element for the non-destructible region; The displacement field of a single discretized base element is expressed in the following form : wherein is the horizontal displacement component at spatial location , time is the vertical displacement component at spatial location , time is the shape function matrix of the th three-node triangular element, is the displacement unknown vector of the th three-node triangular element.​​ 3. The tunnel rockburst numerical simulation method based on the node continuous-discontinuous deformation analysis according to claim 1, characterized in that, The non-equilibrium force residuals are caused by inertial forces, contact forces, and joint element and triangular element deformations, which are The equilibrium equation in the form of the time residual is represented as follows: wherein, is the residual vector, is the time increment, is the next instant time to be solved, is the external force vector at time instant, denotes inertia force vector at time instant, is the contact force vector at time instant, is the cohesion force vector at time instant, denotes internal force vector at time instant.

4. The tunnel rockburst numerical simulation method based on the continuum-discontinuum deformation analysis of nodes according to claim 1, wherein, the construction process of the bilinear constitutive model comprises the following steps: inserting springs along the normal direction and the tangential direction between adjacent discretization basic elements, equivalent simulating the bonding effect between adjacent discretization basic elements through the mechanical response of the springs to form a mechanical characterization basis of the joint element; correspondingly, the working conditions at least include a loading working condition, an unloading working condition and a destruction working condition; the unloading working condition includes an unloading sub-working condition and a loading sub-working condition.

5. The tunnel rockburst numerical simulation method based on the continuum-discontinuum deformation analysis of nodes according to claim 1, wherein, Based on the bilinear constitutive model, the full linear relationship of the joint element in the normal direction is described as follows: in the pre-peak loading segment , post-peak loading segment , post-peak unloading segment and failure segment . ; wherein, represents the normal stress of the joint element, is a normal penalty parameter, is the normal opening of the joint element; wherein, the normal penalty parameter is calculated as follows: ； wherein and are the pre-peak and post-peak normal penalty parameters, and respectively represent the opening value at the normal peak tensile strength and the maximum normal opening value reached before unloading.

6. The tunnel rockburst numerical simulation method based on the continuum-discontinuum deformation analysis of nodes according to claim 1, wherein, Based on the bilinear constitutive model, the full linear relationship of joint elements in the tangential direction is described as follows: in the pre-peak loading segment , post-peak loading segment , post-peak unloading segment and failure segment . ; wherein, represents the shear stress of the joint element in the tangential direction, is a tangential penalty parameter, is the tangential opening of the joint element; wherein the tangential penalty parameter is calculated as follows: ； wherein and represent the pre- and post-peak tangential penalty parameters, respectively, and represent the slip value at the tangential peak shear strength and the maximum tangential slip value reached prior to unloading, respectively, is the tangential residual shear stress.

7. The tunnel rockburst numerical simulation method based on the continuum-discontinuum deformation analysis of nodes according to claim 1, wherein, the identification method of the joint element in the destruction state comprises the following steps: Real-time monitoring of normal opening at integral points of joint elements and tangential opening in combination with a bilinear constitutive model to give a failure threshold and to establish a criterion, which includes a normal criterion and a tangential criterion if the joint element satisfies at least one of the normal judgment condition or the tangential judgment condition, it is determined that the joint element is in a destruction state as a whole; Wherein, the normal determination condition is: if both integral points have , it indicates that the normal stress exceeds the rock mass tensile strength , and it is determined that the corresponding joint element is normal damage. represents the normal opening value when the normal stress reaches the rock mass tensile strength . The tangential judging condition is: if both integral points have , it means that the tangential shear stress has exceeded the rock mass shear strength , and it is judged that the corresponding joint element is tangential failure. , it means that the tangential stress reaches the rock mass shear strength , and the tangential residual slip is represented. 8.The tunnel rock burst numerical simulation method based on the node continuous-discontinuous deformation analysis according to claim 1, wherein, the updating process of the contact potential comprises the following steps: After deleting the joint units in a damaged state, identify adjacent blocks. and blocks Overlapping areas In overlapping areas The total normal contact force is obtained using the following formula. : ; wherein, represents a normal contact penalty parameter, and are contact potential functions of the block and the block respectively, is the gradient operator, and are local coordinates of the infinitesimal area with respect to the block and the block respectively, is the area of the overlapping region ; Based on the total normal contact force The tangential contact force is updated using the following equation : ; where, represents the friction force at the time instant, is the tangential contact penalty parameter, is the tangential relative velocity, is the time increment, is the length of the contact area; The calculated tangential contact force The determination and correction are performed, and if The tangential contact force is corrected using the following equation : ; where is the friction coefficient of the contact surface, is the norm. 9.The tunnel rock burst numerical simulation method based on the continuum-discontinuum deformation analysis of nodes according to claim 1, wherein, the construction process of the rockburst simulation model comprises the following steps: based on the bilinear constitutive model, introducing artificial damping to carry out dynamic analysis to generate an initial stress equilibrium state; applying displacement boundary conditions by using springs, and constructing an initial stress equilibrium state equation by using a spring force equation, the spring force equation being: ; where and are the stiffness coefficients of the springs components and components, denotes the spring reaction force in the direction denotes the spring reaction force in the direction removing the springs at the edges of the calculation region for applying displacement conditions, and calculating the corresponding reaction forces at the edges by using the initial stress equilibrium state equation; Pairs of dampers are introduced at the edges of the computational domain to model the removal of the spring, which creates a viscous boundary by normal drag forces and tangential drag forces equations, in particular: wherein is the density of the rock mass, and is the viscous boundary absorption coefficient and is set uniformly, and are the velocities of the and waves in the rock mass, respectively, is the normal relative velocity, is the tangential relative velocity.

10. The tunnel rockburst numerical simulation method based on the continuum-discontinuum deformation analysis of nodes according to claim 3, wherein, the linear expression of the equilibrium equation is as follows: ; where, is the equivalent stiffness matrix, is the net force of the external force minus all the resistance forces, which is the actual driving force to make the rock mass produce new displacement increment .

Images