An underground chamber surrounding rock thermal-mechanical coupling fatigue phase field simulation method and system

Abstract

This invention belongs to the field of geological energy storage technology, specifically providing a method and system for simulating the thermal-mechanical coupling fatigue phase field of surrounding rock in underground chambers. The method includes: constructing phase field evolution equations using a phase field fracture model and fatigue degradation functions; calculating heat exchange in the tunnel walls based on the calculation model and Newton's law of heat exchange according to the air temperature and pressure inside the tunnel, thus constructing a tunnel wall heat conduction equation; constructing a rock mass deformation equation using rock mass displacement and surrounding rock strain based on the calculation model; numerically discretizing the phase field evolution equations, tunnel wall heat conduction equations, and rock mass deformation equations on a finite element mesh, and solving all discretized equations using a thermal-mechanical coupling method to obtain the evolution results of the temperature field, stress field, displacement field, and crack phase field of the target area of ​​the surrounding rock in the underground chamber. This invention applies the phase field method to the thermal-mechanical coupling simulation of crack propagation in the surrounding rock of underground gas storage tunnels, maximizing the realism of the tunnel crack simulation.

Description

Technical Field

[0001] This invention belongs to the field of geological energy storage technology, specifically relating to a method and system for simulating the thermal-mechanical coupling fatigue phase field of surrounding rock in underground caverns. Background Technology

[0002] Large-scale compressed air energy storage (CASS) offers advantages such as large scale, long lifespan, and low cost, making it a preferred choice for addressing the challenges of large-scale utilization of wind / solar energy fluctuations in new power systems. CASS storage devices typically utilize artificially excavated chambers to avoid limitations imposed by geological resources. In large-scale CASS systems, the storage unit is the core component, and its performance directly determines the system's operational efficiency and safety stability. To overcome the dependence of natural geological gas storage facilities on geological resources, artificially excavated chambers have become the mainstream construction method for large-scale CASS storage facilities. These storage facilities require long-term operation under high-frequency circulation (typically once a day for both inflation and deflation) and high-pressure conditions (working pressure range 5-16 MPa). During this process, the surrounding rock of the storage facility continuously bears the coupled effects of cyclic internal pressure loads and thermal expansion loads caused by temperature field changes. Long-term thermo-mechanical coupled cyclic loading can easily lead to fatigue damage evolution behaviors such as microcrack initiation, propagation, and macro-crack penetration within the surrounding rock of gas storage facilities. This significantly reduces the mechanical properties and integrity of the surrounding rock, directly threatening the long-term operational safety of the gas storage facility and severely affecting its design service life. Among existing technologies, the extended finite element method (EFEM) is widely used in numerical simulation of crack propagation in the surrounding rock of underground compressed air storage chambers. While this method can effectively simulate the geometric propagation path of cracks, it does not consider the influence of fatigue damage accumulation on the crack propagation rate and mode, making it difficult to accurately reflect the fatigue failure mechanism of the surrounding rock under long-term cyclic loading. Physical model tests, while able to directly reflect the damage evolution process of the surrounding rock, have limitations such as long testing cycles, high testing costs, and difficulty in accurately simulating boundary conditions, failing to meet the practical needs of rapid prediction of gas storage facility lifespan during the engineering design phase.

[0003] In summary, a fatigue simulation method for gas storage surrounding rock that can take into account the thermo-mechanical coupling effect and fatigue damage accumulation is established, enabling accurate prediction of the evolution law of fatigue damage in surrounding rock and the service life of gas storage. Summary of the Invention

[0004] To overcome the shortcomings of existing technologies, this invention provides a method and system for simulating the thermal-mechanical coupling fatigue phase field of surrounding rock in underground chambers. The phase field method is applied to the thermal-mechanical coupling simulation of crack propagation in the surrounding rock of underground gas storage chambers, so as to reflect the authenticity of the chamber cracks to the greatest extent.

[0005] To achieve the above objectives, the present invention provides the following solution: A method for simulating the thermo-mechanically coupled fatigue phase field of surrounding rock in underground caverns includes: A computational model of the target area of ​​the surrounding rock of the underground chamber is constructed; the computational model includes the geometric model of the target area, attribute parameters, boundary conditions, and initial conditions; The geometric model is meshed to obtain finite element mesh elements; Based on the aforementioned calculation model, a phase field evolution equation for describing the changes in surrounding rock fractures is constructed using the phase field fracture model and fatigue degradation function. Based on the aforementioned calculation model, and according to the air temperature and pressure inside the tunnel, Newton's law of heat exchange is used to calculate the heat exchange of the tunnel wall, and a tunnel wall heat conduction equation is constructed to describe the temperature field distribution of the surrounding rock. Based on the aforementioned calculation model, a rock mass deformation equation is constructed to describe the mechanical response of the surrounding rock using rock mass displacement and surrounding rock strain. The phase field evolution equation, the cave wall heat conduction equation, and the rock mass deformation equation are numerically discretized on the finite element mesh. The discretized equations are then solved by thermo-mechanical coupling using an iterative method to obtain the temperature field, stress field, displacement field, and crack phase field evolution results of the target area of ​​the underground cavern surrounding rock. Based on the crack phase field evolution results, the scale of crack formation in the cavern surrounding rock and the service life of the gas storage tank are determined.

[0006] Preferably, the method for constructing the phase field evolution equation includes: Based on the variational principle of fractured rock mass, a total energy functional is established that includes the elastic strain energy of the rock mass and the dissipated energy of the fracture. A degradation function describing the stiffness of the rock mass during fracture propagation and a fracture energy density function describing the characteristics of the fracture surface are introduced into the total energy functional to construct a phase-field fracture model. Based on the bulk modulus, Poisson's ratio, and deviatoric strain tensor of the rock mass, the positive and negative parts of the strain energy density are calculated. Combined with the degradation function, the rock mass strain energy density during fracture propagation is obtained. A fatigue degradation function describing the fatigue loss of the rock mass is introduced into the phase-field fracture model. Combining the phase-field fracture model, the rock mass strain energy density, and the fracture phase field, the phase-field evolution equation is constructed.

[0007] Preferably, the method for numerically discretizing the phase field evolution equation includes: The equivalent weak integrals of the phase field evolution equation and phase field boundary conditions are calculated. Based on the phase field stiffness matrix, phase parameters, and phase source vector of the finite element mesh element, the phase parameter stiffness equation of the equivalent weak integral in each finite element mesh element is constructed. By combining the phase parameter stiffness equations of all finite element mesh elements, the overall phase parameter stiffness equation is obtained, and the numerical discretization of the phase field evolution equation is completed.

[0008] Preferably, the method for numerically discretizing the rock mass deformation equation includes: Based on the elastic and strain matrices of the finite element mesh elements, an element stiffness matrix is ​​constructed. Based on the elastic and strain matrices of the finite element mesh elements, the initial strain generated by thermal deformation, and the number of element nodes and shape functions, an equivalent nodal load array is constructed. Based on the element stiffness matrix, the equivalent nodal load array, and the rock mass displacement vector, an element discretization scheme for the rock mass deformation equation is constructed, and based on the element discretization scheme, the mechanical stiffness equations of all finite element mesh elements are obtained. By combining the mechanical stiffness equations of all finite element mesh elements, the overall mechanical stiffness equation is obtained, thus completing the numerical discretization of the rock mass deformation equation.

[0009] Preferably, the method for numerically discretizing the heat conduction equation of the cave wall includes: The equivalent weak integral of the heat conduction equation of the cave wall is calculated; based on the temperature stiffness matrix and the temperature equivalent node matrix of the finite element mesh element, the heat conduction stiffness equation of the equivalent weak integral of the heat conduction equation of the cave wall is obtained; the heat conduction stiffness equations of all finite element mesh elements are combined to obtain the overall heat conduction stiffness equation, thus completing the numerical discretization of the heat conduction equation of the cave wall.

[0010] The present invention also provides a thermal-mechanical coupling fatigue phase field simulation system for the surrounding rock of underground chambers, used to implement the method, comprising: The computational model construction module is used to construct a computational model of the target area of ​​the surrounding rock of the underground chamber; the computational model includes the geometric model of the target area, attribute parameters, boundary conditions, and initial conditions; The mesh generation module is used to perform mesh generation on the geometric model to obtain finite element mesh elements; The phase field evolution module is used to construct a phase field evolution equation to describe the changes in surrounding rock fractures based on the calculation model, using the phase field fracture model and fatigue degradation function. The heat conduction module is used to calculate the heat exchange of the tunnel wall based on the calculation model and the air temperature and pressure inside the tunnel, using Newton's law of heat exchange, and to construct a tunnel wall heat conduction equation to describe the temperature field distribution of the surrounding rock. The rock mass deformation module is used to construct a rock mass deformation equation to describe the mechanical response of the surrounding rock based on the calculation model, using rock mass displacement and surrounding rock strain. The coupled solution module is used to numerically discretize the phase field evolution equation, the cave wall heat conduction equation, and the rock mass deformation equation on the finite element mesh, and to perform thermo-mechanical coupled solution on all the discretized equations through an iterative method to obtain the temperature field, stress field, displacement field, and crack phase field evolution results of the target area of ​​the underground cavern surrounding rock; based on the crack phase field evolution results, the scale of crack generation in the cavern surrounding rock and the service life of the gas storage tank are determined.

[0011] Preferably, the phase field evolution module includes: a phase field fracture model construction unit, used to establish a total energy functional including the elastic strain energy of the rock mass and the dissipated energy of the fracture based on the variational principle of fractured rock mass, and to introduce a degradation function for describing the stiffness of the rock mass during fracture propagation and a fracture energy density function for describing the characteristics of the fracture surface into the total energy functional to construct a phase field fracture model; a strain energy density calculation unit, used to calculate the positive and negative parts of the strain energy density based on the bulk modulus, Poisson's ratio, and deviatoric strain tensor of the rock mass, and to obtain the rock mass strain energy density during fracture propagation by combining the degradation function; and a phase field evolution equation construction unit, used to introduce a fatigue degradation function for describing the fatigue loss of the rock mass into the phase field fracture model, and to construct the phase field evolution equation by combining the phase field fracture model, the rock mass strain energy density, and the fracture phase field.

[0012] Preferably, the coupled solution module includes a phase field discretization unit for numerically discretizing the phase field evolution equation; the phase field discretization unit includes: an equivalent weak integral calculation sub-unit for calculating the equivalent weak integral of the phase field evolution equation and the phase field boundary conditions; a stiffness equation construction sub-unit for constructing the phase parameter stiffness equation of the equivalent weak integral in each finite element mesh element based on the phase field stiffness matrix, phase parameters, and phase source vector of the finite element mesh element; and a stiffness equation combination sub-unit for obtaining the overall phase parameter stiffness equation by combining the phase parameter stiffness equations of all finite element mesh elements, thereby completing the numerical discretization of the phase field evolution equation.

[0013] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. Applying the phase-field method to the numerical simulation of surrounding rock in underground chambers improves the realism of the simulation. Traditional methods such as the finite element method cannot automatically generate initial microcracks at the start of the simulation, requiring pre-setting the potential locations of initial microcracks. This invention, with the phase-field method as its core, systematically describes the application of the phase-field method combined with heat conduction theory in the simulation of surrounding rock in underground compressed air storage chambers. This method can automatically generate initial microcracks through phase-field variables during the simulation process, and the continuous crack propagation process is visualized, allowing for complete observation of cracks and heat diffusion during the working process of the chamber's surrounding rock. 2. To make the simulation results more closely resemble the actual inflation and deflation process of underground chambers, the phase-field method is combined with heat transfer theory. The heat transfer field and mechanical field of this invention are not simply superimposed, but rather represent the interaction of two physical fields. The thermo-mechanical coupled physical field realistically reflects the actual environment during the inflation and deflation process of underground chambers, reconstructing the interaction process of "temperature change - stress action - internal damage evolution" in the chamber's surrounding rock. This coupling is not a simple superposition of two physical fields, but rather an interaction between them, making the simulation results closer to the actual working conditions of underground chambers. During the inflation and deflation of the chamber, temperature changes occur, affecting heat transfer through the chamber walls and consequently the mechanical properties of the surrounding rock. The interaction of the two physical fields jointly influences crack propagation in the surrounding rock and the service life of the chamber. 3. Applying the phase-field method to thermo-mechanical coupling ensures that the interactions of each field are based on explicit physical laws, avoiding artificial assumptions about individual physical fields. This makes the simulation results more consistent with the physical process and improves their reliability. Attached Figure Description

[0014] To more clearly illustrate the technical solution of the present invention, the drawings used in the embodiments are briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0015] Figure 1 This is a schematic diagram illustrating the model calculation method in an embodiment of the present invention; Figure 2 This is a graph showing the relationship between air pressure and the number of cycles in an embodiment of the present invention; Figure 3 This is a flowchart illustrating the fatigue mechanics solution process in an embodiment of the present invention. Figure 4 This is a flowchart illustrating the heat conduction solution process according to an embodiment of the present invention. Figure 5 This is a flowchart of the thermo-mechanical coupling solution in an embodiment of the present invention; Figure 6 This is a schematic diagram of the model dimensions in an embodiment of the present invention; Figure 7 This is a schematic diagram of a grid according to an embodiment of the present invention; Figure 8 This is a crack propagation diagram at 7.5s in an embodiment of the present invention; Figure 9 This is a crack propagation diagram at 102.5 s in Embodiment 1 of the present invention; Figure 10 This is a crack propagation diagram at 250s in an embodiment of the present invention; Figure 11 This is a stress-strain curve diagram of an embodiment of the present invention; Figure 12 This is a diagram showing the stress-strain curves at the valley and peak values ​​according to an embodiment of the present invention; wherein, Figure 12 (a) in the figure is the stress-strain curve at the trough of the sinusoidal cycle; Figure 12 (b) in the figure is the stress-strain curve at the peak of the sinusoidal cycle; Figure 13 This is a heat diffusion diagram at 5 seconds in an embodiment of the present invention; Figure 14 This is a heat diffusion diagram at 7.5s in an embodiment of the present invention; Figure 15 This is a heat diffusion diagram for 250 seconds in an embodiment of the present invention. Detailed Implementation

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

[0017] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0018] Example 1: like Figure 1 As shown, a method for simulating the thermal-mechanical coupling fatigue phase field of surrounding rock in an underground chamber includes: S1: Construct a computational model of the target area of ​​the surrounding rock of the underground chamber; the computational model includes the geometric model of the target area, attribute parameters, boundary conditions, and initial conditions; in this embodiment, before starting the simulation calculation, necessary basic assumptions are made to the model. The basic assumptions in this embodiment are as follows: 1. Rock deformation satisfies the small deformation assumption. 2. The rock is assumed to be an isotropic linear elastic material. 3. The energy dissipation process is irreversible. 4. The loading method is uniform loading over time. 5. It is assumed that there is heat exchange between the wall temperature of the air filling and de-filling chamber and the air temperature, satisfying Newton's law of heat exchange. 6. It is assumed that the temperature at the end of air filling and de-filling forms a stable propagation within the surrounding rock.

[0019] The method of this invention is mainly used to simulate the fracturing and deformation characteristics of surrounding rock under the gas filling and releasing cycle of underground hard rock gas storage, such as... Figure 1 As shown. Figure 2 The diagram illustrates the changes in pressure and temperature with the number of cycles. During the filling and high-pressure storage of the chamber, the increased pressure leads to greater deformation of the surrounding rock and internal damage. Additionally, the compressed air increases its temperature, which propagates through the tunnel walls to the surrounding rock, causing further deformation and damage due to thermal expansion. Further venting and low-pressure storage reduce the pressure and decrease the deformation of the surrounding rock, but the damage intensifies. The air temperature decreases, and the surrounding rock temperature propagates into the chamber. Under long-term, full-cycle operation, the damage to the surrounding rock further worsens, leading to fracturing and crack propagation. This numerical calculation method simulates the fracturing and deformation patterns of the surrounding rock in the entire gas storage facility, providing support for the safety and service life evaluation of the gas storage facility.

[0020] S2: Mesh the geometric model to obtain finite element mesh elements.

[0021] S3: Based on the computational model, using the phase-field fracture model and fatigue degradation function, a phase-field evolution equation is constructed to describe the changes in surrounding rock fractures; a further implementation method includes the following steps for constructing the phase-field evolution equation: Based on the variational principle of fractured rock mass, a total energy functional is established, incorporating the elastic strain energy of the rock mass and the dissipated energy of the fractures. A degradation function describing the rock mass stiffness during fracture propagation and a fracture energy density function describing the characteristics of the fracture surface are introduced into this total energy functional to construct a phase-field fracture model. Specifically, the phase-field fracture model is introduced from three aspects: the internal energy during fracture propagation, the rock mass strain energy density during fracture propagation, and stress. The calculation model is a rock mass region Ω containing discrete fracture surfaces Γ. The fracture surface Γ is expandable under external conditions and is defined as the inner boundary. According to the variational principle of fractured rock mass, the internal energy during fracture propagation is divided into the elastic deformation energy of the rock mass and the dissipated energy of the fractures. The total energy functional L of these two components is calculated as follows: , (1) In the formula, For elastic strain energy density, Let p represent the strain of the rock, p represent the direction of the normal to the surface acting on it, and q represent the direction of the strain. Let V represent the critical energy release rate of the rock material, V represent the surrounding rock volume, and A represent the crack area. A phase field d containing the fracture and the crack width is used. The fracture energy density function γ(d,▽d) describes the characteristics of the fracture surface Γ.

[0022] Introducing the degenerate function g(d) to describe the decrease in rock mass stiffness during fracture propagation, equation (1) becomes: , (2) In the formula, the degenerate function g(d) = (1-d) 2 The fracture energy density function is constructed using the AT2 model. .

[0023] Based on the bulk modulus, Poisson's ratio, and deviatoric strain tensor of the rock mass, the positive and negative components of the strain energy density are calculated. Combined with a degenerate function, the strain energy density of the rock mass during fracture propagation is obtained. Specifically, tensile and shear stresses are the main factors driving fracture propagation. Therefore, for linearly elastic isotropic materials, the strain energy density can be decomposed into a positive component... and negative part , , (3) , (4) In the formula, K is the bulk modulus of the rock, and μ is the Poisson's ratio of the rock. The brackets are for Macaulay. For the partial strain tensor, The trace represents the strain matrix, which is the sum of the diagonals.

[0024] , (5) In the formula; For Kronecker notation, represents the principal strain; n is the dimension. The degenerate function only affects the positive part of the strain energy density; therefore, the strain energy density of the rock mass during fracture propagation... and stress They are respectively: , (6) , (7) A fatigue degradation function describing rock mass fatigue loss is introduced into the phase-field fracture model. Combining the phase-field fracture model, rock mass strain energy density, and fracture phase field, a phase-field evolution equation is constructed. Specifically, to incorporate the fatigue behavior of rock materials into the phase-field solution framework, the fatigue degradation function α(t) is used to characterize the critical energy release rate of rock materials under cyclic loading. The decrease. The fatigue degradation function α(t) is defined as, (8) In the formula, He() is the Heaviside function. This represents the threshold that triggers fatigue failure. represents the derivative of the positive part of the strain energy density, and t represents the cycle time. When When >0, He()=1; otherwise, He()=0. Its physical meaning is that when... When the load is greater than 0, the rock mass is under loading and exhibits irreversible fatigue damage. The degree of fatigue damage gradually increases with each loading cycle. When <0, the rock mass is in an unloaded state and there is no fatigue damage. After introducing the fatigue degradation function, equation (2) becomes: (9) The fracture phase field d can be described by the Allen-Cahn equation. , (10) In the formula, η is the shift parameter (η≥0). Represents partial differential equations. Representing the variation. Integrating equations (2), (6), and (9), we obtain: (11) When η=0, it is considered quasi-static. Since compressive stress does not lead to crack propagation, a long-term variable H is used: , (12) Then equation (11) becomes: (13) Meanwhile, the boundary conditions for the phase field are: , (14) In the formula, n is the normal vector of the boundary. It represents the divergence of the phase field variables.

[0025] S4: Based on the computational model, and according to the air temperature and pressure inside the chamber, Newton's law of heat exchange is used to calculate the heat exchange in the tunnel walls, constructing a heat conduction equation for the tunnel walls to describe the temperature field distribution of the surrounding rock. Specifically, the rapid dynamic changes in the air inside the gas storage tank affect the temperature inside the tank. The temperature inside the gas storage tank is analyzed using Kushnir. and air pressure for: Inflation and high-pressure storage stage: , (15) Venting and low-pressure storage stage: , (16) The density of air inside the gas storage chamber after it begins operation for: , (17) The air pressure inside the chamber is as follows: , (18) In the formula, The air pressure inside the chamber. Let V be the volume of the chamber. The temperature inside the tunnel. Let R be the specific heat capacity of air, R be the gas constant of air, and Z be the compressibility factor of air. For the air inflation rate, c is the air release rate. p Where N is the isobaric specific heat capacity, and N is the number of charge-discharge cycles. This represents the time taken for the number of iterations N.

[0026] There is heat exchange between the cave walls and the hot air, which satisfies Newton's law of heat exchange. , (19) In the formula, Q represents the amount of heat exchanged between the hot air and the cave wall. The heat exchange coefficient between air and cave walls. The contact area between air and the cave wall is assumed to be the ground temperature. .

[0027] The temperature conduction from the cave wall into the surrounding rock satisfies the following equation: , (20) In the formula, β is the porosity. is the thermal conductivity. Divergence representing the dimension of the surrounding rock.

[0028] S5: Based on the computational model, a rock mass deformation equation is constructed using rock mass displacement and surrounding rock strain to describe the mechanical response of the surrounding rock; the deformation of the rock mass satisfies the linear elastic equation. , (twenty one) , (twenty two) In the formula, u represents the rock mass displacement, D is the elastic matrix, and ε is the surrounding rock strain. This is the strain caused by thermal expansion. Represents stress. Boundary conditions are displacement boundaries and stress boundaries.

[0029] S6: The phase field evolution equation, the cave wall heat conduction equation, and the rock mass deformation equation are numerically discretized on a finite element mesh. The discretized equations are then solved by thermo-mechanical coupling using an iterative method to obtain the temperature field, stress field, displacement field, and crack phase field evolution results of the target area of ​​the underground cavern surrounding rock. Based on the crack phase field evolution results, the scale of crack generation in the cavern surrounding rock and the service life of the gas storage tank are determined.

[0030] A further implementation method involves numerically discretizing the phase-field evolution equations, including: Calculate the equivalent weak integrals of the phase field evolution equation and the phase field boundary conditions; specifically, for the phase field evolution equation (13) and the boundary condition equation (14), their equivalent weak integral forms can be expressed as follows: (twenty three) In the formula, This is a trial function for the phase-field evolution equation. It employs... As a trial function For shape functions, triangular units are used here.

[0031] Based on the phase field stiffness matrix, phase parameters, and phase source vector of the finite element mesh element, the equivalent weak integral phase parameter stiffness equation for each finite element mesh element is constructed; specifically, equation (23) is the phase parameter stiffness equation for each finite element mesh element as follows: ,(twenty four) , , , These are the phase field stiffness matrices of the element. R is the phase parameter of the element. e is the phase source vector of the element.

[0032] By combining the phase parameter stiffness equations of all finite element mesh elements, the overall phase parameter stiffness equation is obtained, thus completing the numerical discretization of the phase field evolution equation. Further derivation yields the expressions for the phase field stiffness matrix of the finite element mesh elements as follows: , (25) , (26) , (27) , (28) , (29) In the formula, Represents the direction of coordinate i. The j-direction component of the shape function. The m-direction component of the shape function. The normal vector represents the direction of the shape function i.

[0033] By combining the phase parameter stiffness equations of all elements, the overall phase parameter stiffness equation is obtained. , (30) The differential represents the phase field variable.

[0034] Since the overall stiffness equation contains a time term, equation (30) is rewritten in an iterative format. , (31) The calculation can be performed using the post-Euler explicit solution method.

[0035] A further implementation method involves numerically discretizing the rock mass deformation equation, including: Based on the elastic and strain matrices of the finite element mesh elements, the element stiffness matrix is ​​constructed. Based on the elastic and strain matrices of the finite element mesh elements, the initial strain generated by thermal deformation, and the number of element nodes and shape functions, the element equivalent nodal load array is constructed. Based on the element stiffness matrix, the element equivalent nodal load array, and the rock mass displacement vector, the element discretization scheme of the rock mass deformation equation is constructed, and based on the element discretization scheme, the mechanical stiffness equations of all finite element mesh elements are obtained. By combining the mechanical stiffness equations of all finite element mesh elements, the overall mechanical stiffness equation is obtained, thus completing the numerical discretization of the rock mass deformation equation.

[0036] Specifically, considering the heat conduction effect, a discrete solution format (element stiffness matrix) for the rock mass deformation equation is established according to the finite element analysis process. The element discretization format is as follows: , (32) In the formula, The element stiffness matrix, The element is the equivalent nodal load array. This is the displacement vector.

[0037] use As a trial function For shape functions, triangular units are used here. This yields... , They are respectively: , (33) , (34) In the formula, D is the elasticity matrix and B is the strain matrix. Represents physical strength. Represents surface force.

[0038] , (35) Let be the initial strain caused by thermal expansion. For two-dimensional deformation, its expression is: , (36) Represents the temperature of the surrounding rock. This represents the coefficient of thermal expansion.

[0039] By combining the mechanical stiffness equations of all elements, the overall mechanical stiffness equation is obtained. (37) The above equations are a large-scale linear system of equations in a band, which can be solved using the preconditional conjugate gradient method.

[0040] A further implementation method involves numerically discretizing the heat conduction equation of the cave wall, including: The equivalent weak integral of the cave wall heat conduction equation is calculated; based on the temperature stiffness matrix and temperature equivalent node matrix of the finite element mesh element, the heat conduction stiffness equation of the equivalent weak integral of the cave wall heat conduction equation is obtained; the heat conduction stiffness equations of all finite element mesh elements are combined to obtain the overall heat conduction stiffness equation, thus completing the numerical discretization of the cave wall heat conduction equation. Specifically, the cave wall heat conduction equation (20) is also discretized using the finite element method, and the equivalent weak integral form of the cave wall heat conduction equation (20) is first obtained: , (38) In the formula, This is a trial function. It uses... As a trial function The shape function can be selected according to the principle of finite element shape function construction; here, triangular elements are used. Represents porosity. Let x represent thermal conductivity, and y represent the coordinate axes. Then, the discrete form of equation (32) on each element is: , (39) In the formula Here is the temperature stiffness matrix. Let be the equivalent node matrix for temperature.

[0041] , (40) , (41) The j-direction component of the shape function. This represents the m-direction component of the shape function.

[0042] By combining the heat conduction equations of all the units in the cavity wall, the overall heat conduction stiffness equation is obtained. , (42) The above equations are a large-scale linear system of equations in a band, which can be solved using the preconditional conjugate gradient method.

[0043] The specific methods for obtaining the evolution results of the temperature field, stress field, displacement field, and fracture phase field of the target area of ​​the surrounding rock of the underground chamber by solving all discretized equations through thermo-mechanical coupling using iterative methods include: like Figure 3 The diagram shows the solution process for fatigue mechanics. First, the rock displacement field is defined. Phase field parameters Fatigue degradation function and long-term evolution variables Then, the elastic deformation equation is solved, and the stress σ and displacement u are initialized (as defined). =0, =0). Calculate the element stiffness matrix and equivalent nodal load array according to equations (26), (27), and (28), and then integrate all element stiffness matrices according to equation (30) to obtain the overall stiffness matrix. Check whether the displacement converges by calculating using the preconditional conjugate gradient method. If it does not converge, re-initialize and recalculate the stress and displacement. If it converges, solve the phase field evolution equation.

[0044] Using the updated stress field and displacement field as the initial values ​​for fatigue iteration, the long-term variable H, fatigue degradation function α(t), element phase field matrix, and global phase field matrix are calculated according to equations (12), (8), (18), and (25), respectively. The convergence of the above variables is checked by the preconditional conjugate gradient method. If they do not converge, iterative calculation cannot be performed. At this time, the stress field needs to be calculated and updated according to equation (7). If they converge, the phase field variables are updated according to the above calculation results and the stress field is updated according to equation (7). Iterative calculation is performed for the next time step until the last time step is reached, the fatigue cycle ends, and the calculation results are output.

[0045] like Figure 4The solution process for the heat conduction module is shown below. First, the temperature solution is initialized according to the actual working conditions, such as determining the initial temperature value by combining the surface temperature and geothermal gradient. Then, appropriate boundary conditions are applied to each boundary: for example, the temperature of the chamber wall is set to a temperature that alternates between hot and cold with the sinusoidal change of pressure or a boundary with a continuously fixed temperature; the outer boundary of the model is set to an insulation layer or a temperature boundary with a fixed initial value. The stiffness matrix of the unit and the equivalent node array are calculated according to equations (40) and (41), and combined into the overall stiffness matrix by equation (30). The precondition conjugate gradient method is used for calculation to check whether the calculated temperature value converges: if it does not converge, the surrounding rock temperature field is updated and the overall stiffness matrix is ​​recalculated; if it converges, the calculation is ended and the final temperature field distribution is obtained.

[0046] like Figure 5 The diagram shows the flow chart for solving the thermo-mechanical coupling problem. First, the basic structure of the model is constructed, such as geometric dimensions, parameter system, materials, functions, and periods. Based on the accuracy requirements, the region is meshed. Then, physical and mechanical parameters such as elastic modulus, Poisson's ratio, and critical energy release rate, as well as thermodynamic parameters such as porosity, thermal conductivity, and heat exchange coefficient between air and the chamber wall, are input. After the basic construction is completed, the two physical fields are coupled. The calculation step size is defined according to the period to meet the simulation requirements of the preset number of cycles N. The internal pressure and temperature of the chamber thermodynamic parameters are calculated using equations (17) to (19). The temperature field, stress field, and displacement field are calculated using the above fatigue mechanics and heat conduction methods. If the above three fields do not converge, the temperature field, displacement field, and fatigue phase field are updated and iteratively calculated until the results of the temperature field, stress field, and displacement field all converge. Then, the calculation is stopped, and the final solution of the multi-field coupling is output.

[0047] Example 2: This embodiment uses the thermo-mechanical coupling solution method of Embodiment 1 to simulate the propagation of surrounding rock fissures during the gas filling and releasing process of an underground gas storage facility.

[0048] Example Introduction: This example simulates the gas filling and releasing process of a chamber located in the middle of the model. The surrounding rock mass is Late Hercynian intrusive granodiorite, with local carbonate interlayers of poor integrity. This model is simplified to simulate ordinary granite with good integrity.

[0049] This model is a two-dimensional numerical calculation model, and the calculation and analysis are carried out using the working condition of a burial depth of 135m and a tunnel diameter of 10m. First, as... Figure 6As shown, a rectangular computational domain with a length of 215m and a width of 140m was established, filled with granite material. A circular region with a radius of 10m was constructed, centered at a point 145m from the top surface of the computational domain and 70m from the left, right, and bottom surfaces, to simulate the chamber; this region was filled with air. The cyclic loading model adopted a load and temperature cycling mode with a period of T=5s, and the total calculation time was set to 250s, with the first 5s being a linear increase phase of pressure and temperature. The model completed a total of 49 cyclic loading processes.

[0050] In the mechanics module, the circular chamber wall is divided into 12 equal element segments, and boundary loads are applied to each element segment. The left, right, and lower surfaces of the rectangular computational domain are constrained and fixed with axial supports, while the upper surface is set as a free constraint. Sinusoidal cyclic loads with an initial load of 0 MPa and a period of 5 s, ranging from 2 MPa to 10 MPa, are applied around the circular chamber.

[0051] In the heat transfer module, the initial temperature of all parts outside the circular chamber in the rectangular computational domain is set to 293.15K. The cavity wall temperature is calculated based on the inflation and deflation parameters and equations (15) to (18), and the heat conduction boundary is set as the output condition.

[0052] Grid division: such as Figure 7 As shown, the meshing of the surrounding rock of the chamber is carried out using the "very fine" meshing method in COMSOL Multiphysics 6.3.

[0053] Simulation results: (1) Fatigue crack propagation analysis: such as Figure 8 As shown, at the first peak moment after the start of cyclic loading (7.5s), the fatigue crack propagation is limited, and its affected area is relatively narrow; as Figure 9 As shown, when the loading time reaches 102.5 s, the damage area affected by fatigue shows a significant expansion trend compared to 7.5 s; as Figure 10 As shown, at the end of the cyclic loading (250s), the rock damage area was still significantly larger than at 102.5s, but its propagation rate slowed down. (2) Stress-strain curve: as shown Figure 11 As shown, the stress-strain curve clearly exhibits two stages. In the first stage, the load and temperature increase linearly, and the strain increases linearly accordingly; this stage is dominated by elastic deformation. Subsequently, the fatigue cyclic stage begins, where the strain changes are subtle, and the changes are difficult to show on a conventional scale. Therefore, the valleys and peak values ​​during the cyclic loading process are analyzed in a magnified manner. Figure 12 As shown, in Figure 12 (a) and Figure 12In (b), the stress-strain relationship shows a good linear trend. At the trough and peak of the sinusoidal cycle, the rock mass around the chamber is still in the elastic deformation stage that conforms to Huke's law. In addition, as the number of cycles continues to increase, the strain of the rock mass around the chamber has a cumulative effect. When the cumulative strain reaches a certain critical value, the gas storage chamber will rupture, that is, the chamber will fail. (3) Temperature propagation characteristics: such as Figure 13 As shown, the temperature field linearly increases from 293.15K to 298.15K at 5s, followed by a cyclic loading phase; as... Figure 14 As shown, at 7.5s, the temperature is at its peak within the sinusoidal cycle, exhibiting significant temperature changes. The temperature around the circular chamber rises to a maximum of 348.15K, showing a gradual decrease from the periphery outwards. Figure 15 As shown, the minimum point of the cycle is at the end of the cycle (250s). Therefore, the temperature field distribution at this time is very similar to the temperature field distribution at 5s. This is because the thermal conductivity of granite is only 2.9W / (m·K). After only 49 cycles, the range of heat propagation is still relatively limited.

[0054] The above describes a phase-field method for simulating the thermo-mechanical coupling fatigue of surrounding rock in underground tunnels. Compared with the discrete element method and the finite element method (FEM), the phase-field method has significant advantages in simulating fatigue in underground tunnels: 1. Underground tunnels are generally located in hard rock strata more than 100 meters underground. The rock properties in these areas are complex, making it difficult to determine potential cracks and failures. This invention provides a detailed simulation method for simulating the fatigue of surrounding rock in underground tunnels based on the phase-field method. Compared with the traditional discrete element method and finite element method, this method does not require pre-setting cracks; it uses damage field variables to describe crack propagation. Cracks, as a result of damage field evolution, automatically initiate when damage accumulates to a certain threshold. This characteristic better reflects the randomness of initial fatigue crack formation in underground tunnels. 2. Based on the thermo-mechanical coupling phenomenon that occurs during fatigue loading in underground tunnels, this invention introduces the heat transfer field into the phase-field method, which acts together with the mechanical field on the surrounding rock of the tunnel. The thermo-mechanical coupling physical field realistically reflects the actual environment during the gas filling and releasing process of underground chambers, reconstructing the interaction process of "temperature change-stress action-internal damage evolution" in the surrounding rock of the chamber. This coupling is not a simple superposition of two physical fields, but rather an interaction between the two fields, making the simulation results closer to the actual working conditions of underground chambers. 3. By applying the phase-field method to thermo-mechanical coupling, the interaction of each field is based on explicit physical laws, avoiding artificial assumptions about each physical field, making the simulation results more consistent with the physical process, and improving the reliability of the simulation results.

[0055] Example 3: This invention also provides a thermal-mechanical coupling fatigue phase-field simulation system for surrounding rock in underground chambers, used to implement the method of Embodiment 1, comprising: a computational model construction module for constructing a computational model of the target area of ​​the surrounding rock in the underground chamber; the computational model includes a geometric model of the target area, attribute parameters, boundary conditions, and initial conditions; a mesh generation module for meshing the geometric model to obtain finite element mesh elements; a phase-field evolution module for constructing a phase-field evolution equation describing the changes in fractures in the surrounding rock based on the computational model, using a phase-field fracture model and fatigue degradation function; and a heat conduction module for simulating the thermal-mechanical coupling fatigue phase-field ... The system employs a calculation module to determine the heat exchange in the tunnel walls and construct a heat conduction equation describing the temperature field distribution of the surrounding rock. A rock mass deformation module, based on the calculation model and utilizing rock mass displacement and surrounding rock strain, constructs a rock mass deformation equation describing the mechanical response of the surrounding rock. A coupled solution module numerically discretizes the phase field evolution equation, tunnel wall heat conduction equation, and rock mass deformation equation on a finite element mesh, and then uses an iterative method to perform a thermo-mechanical coupled solution on all discretized equations to obtain the temperature field, stress field, displacement field, and crack phase field evolution results of the target area of ​​the underground tunnel's surrounding rock. Based on the crack phase field evolution results, the scale of crack formation in the tunnel's surrounding rock and the service life of the gas storage tank are determined.

[0056] A further implementation method is that the phase field evolution module includes: a phase field fracture model construction unit, used to establish a total energy functional containing the elastic strain energy of the rock mass and the dissipated energy of the fracture based on the variational principle of fractured rock mass, and to introduce a degradation function to describe the stiffness of the rock mass during fracture propagation and a fracture energy density function to describe the characteristics of the fracture surface into the total energy functional to construct a phase field fracture model; a strain energy density calculation unit, used to calculate the positive and negative parts of the strain energy density based on the bulk modulus, Poisson's ratio, and deviatoric strain tensor of the rock mass, and to obtain the rock mass strain energy density during fracture propagation by combining the degradation function; and a phase field evolution equation construction unit, used to introduce a fatigue degradation function to describe the fatigue loss of the rock mass into the phase field fracture model, and to construct a phase field evolution equation by combining the phase field fracture model, the rock mass strain energy density, and the fracture phase field.

[0057] A further implementation method is that the coupled solution module includes a phase field discretization unit for numerically discretizing the phase field evolution equation; the phase field discretization unit includes: an equivalent weak integral calculation sub-unit for calculating the equivalent weak integral of the phase field evolution equation and the phase field boundary conditions; a stiffness equation construction sub-unit for constructing the equivalent weak integral phase parameter stiffness equation of each finite element mesh element based on the phase field stiffness matrix, phase parameters, and phase source vector of the finite element mesh element; and a stiffness equation combination sub-unit for obtaining the overall phase parameter stiffness equation by combining the phase parameter stiffness equations of all finite element mesh elements, thus completing the numerical discretization of the phase field evolution equation.

[0058] The above embodiments are merely descriptions of preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Various modifications and improvements made by those skilled in the art to the technical solutions of the present invention without departing from the spirit of the present invention should fall within the protection scope defined by the claims of the present invention.

Claims

1. A method for simulating the thermo-mechanical coupling fatigue phase field of surrounding rock in underground chambers, characterized in that, include: Construct a computational model of the target area of ​​the surrounding rock of the underground chamber; The computational model includes the geometric model of the target region, attribute parameters, boundary conditions, and initial conditions; The geometric model is meshed to obtain finite element mesh elements; Based on the aforementioned calculation model, a phase field evolution equation for describing the changes in surrounding rock fractures is constructed using the phase field fracture model and fatigue degradation function. Based on the aforementioned calculation model, and according to the air temperature and pressure inside the tunnel, Newton's law of heat exchange is used to calculate the heat exchange of the tunnel wall, and a tunnel wall heat conduction equation is constructed to describe the temperature field distribution of the surrounding rock. Based on the aforementioned calculation model, a rock mass deformation equation is constructed to describe the mechanical response of the surrounding rock using rock mass displacement and surrounding rock strain. The phase field evolution equation, the cave wall heat conduction equation, and the rock mass deformation equation are numerically discretized on the finite element mesh. The discretized equations are then solved by thermo-mechanical coupling using an iterative method to obtain the temperature field, stress field, displacement field, and crack phase field evolution results of the target area of ​​the underground cavern surrounding rock. Based on the crack phase field evolution results, the scale of crack formation in the cavern surrounding rock and the service life of the gas storage tank are determined. The method for constructing the phase field evolution equation includes: Based on the variational principle of fractured rock mass, a total energy functional including the elastic strain energy of rock mass and the dissipated energy of fracture is established. A degradation function for describing the stiffness of rock mass during fracture propagation and a fracture energy density function for describing the characteristics of fracture surface are introduced into the total energy functional to construct a phase field fracture model. Based on the bulk modulus, Poisson's ratio, and deviatoric strain tensor of the rock mass, the positive and negative parts of the strain energy density are calculated. Combined with the degradation function, the strain energy density of the rock mass during the fracture propagation process is obtained. A fatigue degradation function for describing rock mass fatigue loss is introduced into the phase field fracture model. The phase field evolution equation is constructed by combining the phase field fracture model, the rock mass strain energy density, and the fracture phase field. The fatigue degradation function α(t) is defined as: In the formula, He() is the Heaviside function. This represents the threshold that triggers fatigue failure. The derivative of the positive part of the strain energy density is represented by t, where t represents the cycle time; when When >0, He()=1, otherwise He()=0; the physical meaning is that when When the load is greater than 0, the rock mass is under loading and exhibits fatigue damage, which is irreversible; with increasing loading cycles, the degree of fatigue damage gradually increases; when... When <0, the rock mass is in an unloaded state and there is no fatigue damage; after introducing the fatigue degradation function, we get: ， g(d) represents the degradation function, Γ represents the discrete fracture surface, and Ω represents the rock mass region. For elastic strain energy density, Let p represent the strain of the rock, p represent the direction of the normal to the surface acting on it, and q represent the direction of the strain. The critical energy release rate of the rock material is given by V, where V represents the volume of the surrounding rock and A represents the area of ​​the crack. The method incorporates the fracture phase field d and the crack width. The fracture energy density function γ(d,▽d) describes the characteristics of the fracture surface Γ. The fracture phase field d is described using the Allen-Cahn equation. ， where η is a shift parameter, denotes partial differentiation, denotes variation; yields When η=0, it is considered quasi-static; since compressive stress does not lead to crack propagation, a long-term variable H is adopted: ， get: Temperature inside the gas storage facility analyzed using Kushnir and air pressure for: Inflation and high pressure storage phase: Deflation and low pressure storage phase: , The density of the air in the chamber at the start of operation of the gas storage is: ， The air pressure inside the chamber is as follows: ， In the formula, The air pressure inside the chamber. Let V be the volume of the chamber. The temperature inside the tunnel. Let R be the specific heat capacity of air, R be the gas constant of air, and Z be the compressibility factor of air. For the air inflation rate, c is the air release rate. p Where N is the isobaric specific heat capacity, and N is the number of charge-discharge cycles. This represents the time taken for the number of iterations N. The methods for numerically discretizing the phase field evolution equations include: Calculate the equivalent weak integrals of the phase field evolution equation and the phase field boundary conditions; Based on the phase field stiffness matrix, phase parameters, and phase source vector of the finite element mesh element, the equivalent weak integral phase parameter stiffness equation of each finite element mesh element is constructed. By combining the phase parameter stiffness equations of all finite element mesh elements, the overall phase parameter stiffness equation is obtained, thus completing the numerical discretization of the phase field evolution equation. The methods for numerically discretizing the rock mass deformation equation include: The element stiffness matrix is ​​constructed based on the elasticity matrix and strain matrix of the finite element mesh element. Based on the elastic matrix, strain matrix, initial strain generated by thermal deformation, and the number of element nodes and shape function of the finite element mesh, an equivalent nodal load array is constructed for the element. Based on the element stiffness matrix, the element equivalent nodal load array, and the rock mass displacement vector, an element discretization scheme for the rock mass deformation equation is constructed, and based on the element discretization scheme, the mechanical stiffness equations of all finite element mesh elements are obtained. By combining the mechanical stiffness equations of all finite element mesh elements, the overall mechanical stiffness equation is obtained, thus completing the numerical discretization of the rock mass deformation equation. The methods for numerically discretizing the heat conduction equation of the cave wall include: Calculate the equivalent weak integral of the heat conduction equation of the cave wall; Based on the temperature stiffness matrix and temperature equivalent node matrix of the finite element mesh, the equivalent weak integral thermal conduction stiffness equation of the cavity wall heat conduction equation is obtained. By combining the thermal conduction stiffness equations of all finite element mesh elements, the overall thermal conduction stiffness equation is obtained, thus completing the numerical discretization of the thermal conduction equation of the cavity wall. The specific methods for obtaining the evolution results of the temperature field, stress field, displacement field, and fracture phase field of the target area of ​​the surrounding rock of the underground chamber by solving all discretized equations through thermo-mechanical coupling using iterative methods include: First, the rock displacement field u0, phase field parameter d0, fatigue degradation function α(t)0, and long-term evolution variable H0 are defined. Then, the elastic deformation equation is solved, and the stress σ and displacement u are initialized. The element stiffness matrix and equivalent nodal load array are calculated, and then all element stiffness matrices are integrated to obtain the global stiffness matrix. The convergence of displacement is checked by calculating using the preconditional conjugate gradient method. If it does not converge, the stress and displacement are re-initialized and recalculated. If it converges, the phase field evolution equation is solved. Using the updated stress field and displacement field as the initial values ​​for fatigue iteration, the long-term variable H, fatigue degradation function α(t), element phase field matrix, and global phase field matrix are calculated respectively. The convergence of the above variables is checked by the preconditional conjugate gradient method. If they do not converge, iterative calculation cannot be performed, and the stress field needs to be calculated and updated. If they converge, the phase field variables and stress field are updated according to the above calculation results, and iterative calculation is performed for the next time step until the last time step is reached, at which point the fatigue cycle ends and the calculation results are output. The temperature solution is initialized based on actual working conditions, and the initial temperature value is determined by combining the surface temperature and geothermal gradient. Subsequently, appropriate boundary conditions are applied to each boundary: the temperature of the chamber wall is set as an alternating hot and cold temperature that varies sinusoidally with pressure or as a continuously fixed temperature boundary; the outer boundary of the model is set as an insulation layer or a temperature boundary with a fixed initial value; the stiffness matrix of the calculated element and the equivalent node array are combined into the global stiffness matrix, and the preconditioned conjugate gradient method is used for calculation to check whether the calculated temperature value converges: if it does not converge, the surrounding rock temperature field is updated and the global stiffness matrix is ​​recalculated; if it converges, the calculation ends and the final temperature field distribution is obtained. Complete the basic structure of the model, including geometry, materials, and period. Mesh the region based on accuracy requirements, and then input physical and mechanical parameters such as elastic modulus, Poisson's ratio, critical energy release rate, porosity, thermal conductivity, and heat exchange coefficient between air and chamber walls. After completing the basic structure, couple the temperature field and stress field, defining the calculation step size according to the period to meet the simulation requirements of the preset number of cycles N. Calculate the internal pressure and temperature of the chamber's thermodynamic parameters, and calculate the temperature field, stress field, and displacement field by combining fatigue mechanics and heat conduction equations. If the three fields do not converge, update the temperature field, displacement field, and fatigue phase field for iterative calculation until the results of the temperature field, stress field, and displacement field all converge. Then stop the calculation and output the final solution of the multi-field coupling.

2. A thermal-mechanical coupled fatigue phase-field simulation system for surrounding rock in underground chambers, used to implement the method described in claim 1, characterized in that, include: The computational model construction module is used to construct a computational model of the target area of ​​the surrounding rock of the underground chamber; the computational model includes the geometric model of the target area, attribute parameters, boundary conditions, and initial conditions; The mesh generation module is used to perform mesh generation on the geometric model to obtain finite element mesh elements; The phase field evolution module is used to construct a phase field evolution equation to describe the changes in surrounding rock fractures based on the calculation model, using the phase field fracture model and fatigue degradation function. The heat conduction module is used to calculate the heat exchange of the tunnel wall based on the calculation model and the air temperature and pressure inside the tunnel, using Newton's law of heat exchange, and to construct a tunnel wall heat conduction equation to describe the temperature field distribution of the surrounding rock. The rock mass deformation module is used to construct a rock mass deformation equation to describe the mechanical response of the surrounding rock based on the calculation model, using rock mass displacement and surrounding rock strain. The coupled solution module is used to numerically discretize the phase field evolution equation, the cave wall heat conduction equation, and the rock mass deformation equation on the finite element mesh, and to perform thermo-mechanical coupled solution on all the discretized equations through an iterative method to obtain the temperature field, stress field, displacement field, and crack phase field evolution results of the target area of ​​the underground cavern surrounding rock; based on the crack phase field evolution results, the scale of crack generation in the cavern surrounding rock and the service life of the gas storage tank are determined.

3. The system of claim 2, wherein, The phase field evolution module includes: The phase-field fracture model construction unit is used to establish a total energy functional that includes the elastic strain energy of the rock mass and the dissipated energy of the fracture based on the variational principle of fractured rock mass. The unit introduces a degradation function to describe the stiffness of the rock mass during fracture propagation and a fracture energy density function to describe the characteristics of the fracture surface into the total energy functional to construct the phase-field fracture model. The strain energy density calculation unit is used to calculate the positive and negative parts of the strain energy density based on the bulk modulus, Poisson's ratio, and deviatoric strain tensor of the rock mass, and to obtain the strain energy density of the rock mass during the fracture propagation process by combining the degradation function. The phase field evolution equation construction unit is used to introduce a fatigue degradation function to describe the fatigue loss of rock mass into the phase field fracture model, and to construct the phase field evolution equation by combining the phase field fracture model, the strain energy density of the rock mass and the fracture phase field.

4. The system of claim 3, wherein, The coupled solution module includes a phase field discretization unit, which is used to numerically discretize the phase field evolution equation; The phase field discrete unit includes: The equivalent weak integral computation sub-unit is used to calculate the equivalent weak integral of the phase field evolution equation and the phase field boundary conditions. The stiffness equation is used to construct sub-elements, which are used to construct the equivalent weak integral phase parameter stiffness equation of each finite element mesh element based on the phase field stiffness matrix, phase parameters and phase source vector of the finite element mesh element. The stiffness equation combination sub-element is used to obtain the overall phase parameter stiffness equation by combining the phase parameter stiffness equations of all finite element mesh elements, thereby completing the numerical discretization of the phase field evolution equation.

Images