Continuous-discontinuous numerical simulation method considering rock mass thermal coupling problem

Abstract

The invention discloses a continuous-discontinuous numerical simulation method considering a rock mass thermal coupling problem, and the method comprises the steps: dividing a simulation region into a potential damage region and a non-damage region, and carrying out the grid division of the simulation region; then generating a thermal calculation node system and a force calculation node system for analog calculation, wherein a mapping relationship exists between the thermal calculation node system and the force calculation node system; jointing units are arranged on the grid boundary of the potential damage area, and the rock mass fracture process is simulated through damage of the jointing units; the simulation process is decomposed into a plurality of calculation steps, the heat conduction process and the mechanical calculation process in each calculation step are solved in sequence, whether the joint unit is damaged or not is judged according to the deformation and stress state of the joint unit, the mapping relation is updated until the last calculation step is finished, and the thermal fracture process simulation of the rock mass is completed. According to the method, artificial heat transfer parameters can be prevented from being introduced in the heat conduction simulation process, the fracture evolution rule of the local area of the rock mass is researched in a targeted mode, the calculation cost of the model is reduced, and the dynamic fracture process in the rock mass thermal coupling problem is effectively simulated.

Description

Technical Field

[0001] The invention belongs to the technical field of geotechnical engineering numerical simulation, and in particular relates to a continuous-discontinuous numerical simulation method considering the thermal-mechanical coupling problem of rock mass. Background Art

[0002] Thermal stress caused by temperature changes can lead to the initiation and expansion of cracks in the rock mass, further causing instability and damage to engineering facilities. Although physical tests can intuitively demonstrate the thermal cracking process of rocks, they are time-consuming, costly, and limited in scale.

[0003] Using numerical methods to simulate the thermal cracking process of rock mass is an effective way to study the destructive behavior of deep rock facilities, which can overcome the difficulties encountered by indoor experimental methods. However, in order to ensure the accuracy of numerical simulation, the coverage area of ​​the calculation model is usually much larger than the area of ​​concern in engineering design. In order to simplify the modeling cost and facilitate program writing, the existing thermo-solid coupling methods, such as the papers "Y.Jiao, X.Zhang, H.Zhang, H.Li, S.Yang, J.Li, 2015, A coupled thermo-mechanical discontinuum model for simulating rock cracking induced bytemperature stresses, Computers and Geotechnics" and "C.Yan, Y.Jiao, 2020, A 2Ddiscrete heat transfer model considering the thermal resistance effect offractures for simulating the thermal cracking of brittle materials, ActaGeotechnica", usually do not perform additional processing on the study area when simulating the thermal cracking process, resulting in unnecessary computational consumption. At the same time, in order to simplify the heat conduction calculation process, additional artificial heat transfer parameters are introduced, which requires additional calibration of the algorithm when simulating rock heat transfer problems. Summary of the invention

[0004] The purpose of the present invention is to provide a continuous-discontinuous numerical simulation method that takes into account the thermal-mechanical coupling problem of rock mass, so as to solve the technical problem that in the prior art, when simulating the thermal cracking process, the study area is usually not additionally processed, resulting in unnecessary computing consumption; at the same time, in order to simplify the heat conduction calculation process, additional artificial heat transfer parameters are additionally introduced, resulting in the need for additional calibration of the algorithm when simulating the rock mass heat transfer problem.

[0005] In order to solve the above technical problems, the present invention adopts the following technical solutions:

[0006] A continuous-discontinuous numerical simulation method considering the thermal-mechanical coupling problem of rock mass comprises the following steps:

[0007] Step S1. According to the stress characteristics of the rock mass, the simulation area of ​​the rock mass model is divided into a potential damage area and an indestructible area, and the potential damage area and the indestructible area are meshed;

[0008] Step S2. setting a joint unit at the boundary of the potential damage zone grid, the joint unit connects the grids on both sides thereof, and simulates the rock mass fracture process through the destruction of the joint unit;

[0009] Step S3. Based on the model mesh division and joint unit insertion results, a thermal calculation node system and a force calculation node system for simulation calculation are generated, and there is a mapping relationship between the two. The mapping relationship is continuously and dynamically updated during the process of simulating rock mass fracture;

[0010] Step S4. Decompose the simulation process of rock mass fracture into several calculation steps, and solve the heat conduction process and the mechanical calculation process in each calculation step in turn;

[0011] Step S5. After each calculation step is completed, determine whether the joint unit is damaged based on its deformation and stress state, and update the mapping relationship between the thermal calculation node system and the force calculation node system;

[0012] Step S6. Repeat the above step S5 until the last calculation step is completed, and the thermal cracking process simulation of the rock mass is completed.

[0013] Further optimization, in step S1, the stress concentration area caused by internal defects or external disturbances of the rock mass is set as a potential damage area, and the remaining areas are set as indestructible areas; for the two-dimensional rock mass model, a triangular computational grid is used to discretize the model.

[0014] Further optimization, in step S2, the computational grids of the potential damage zone are independent of each other, and adjacent computational grids are connected together by inserting joint units, and the failure of the joint units represents the initiation and expansion of real rock cracks.

[0015] Further optimization, in step S3, based on the model mesh division and joint unit insertion results, a thermal calculation node system and a force calculation node system for simulation calculation are generated, specifically:

[0016] When meshing the model in step S1, a large number of nodes are first inserted inside and on the boundary of the model, and three adjacent nodes are connected to each other to form a triangular computational mesh.

[0017] All nodes are divided into two sets of node systems: thermal calculation and force calculation. The nodes in the thermal calculation system are called thermal calculation nodes, and the nodes in the force calculation system are called force calculation nodes. In the indestructible area, the thermal calculation nodes and the force calculation nodes are exactly the same, and there are multiple adjacent triangular calculation meshes sharing the same node. In the potential damage area, due to the insertion of joint units, the number of force calculation nodes is greater than the number of thermal calculation nodes, and one thermal calculation node is associated with at least one or more force calculation nodes.

[0018] Further optimization, in step S4, solving the heat conduction process and the mechanical calculation process in each calculation step specifically includes the following steps;

[0019] Step S4.1. Calculate the heat flowing into each heat calculation node at this moment, specifically:

[0020] 1) When the joint unit is not destroyed, the crack does not exist, and the node corresponding to the thermal calculation node Ti is shared by multiple triangular calculation grids. Equal to the heat that flows into Ti through multiple triangles that share the node. sum;

[0021] 2) When the joint unit is destroyed and there is a crack, the crack divides the multiple triangular computational meshes that share the node into two parts, one on each side of the crack. Calculated by formula (1):

[0022]

[0023] In the above formula, The heat flowing into Ti for the triangles located on the same side of the thermal calculation node Ti is calculated. The heat flowing into node Ti through the crack is calculated for the triangle mesh on the other side of the crack.

[0024] Step S4.2. Calculated according to step S4.1 Substitute into the rock mass heat conduction formula (2) and calculate the temperature of each heat calculation node at this time;

[0025]

[0026] In the above formula, and are the temperatures of the i-th thermal calculation node Ti in the model at time t and t+Δt respectively; is the heat flowing into the i-th node; c p is the specific heat capacity of the rock mass; M i is the node mass corresponding to the i-th thermal calculation node, which is equal to the sum of the masses of all associated force calculation nodes.

[0027] Step S4.3. According to the temperature of each thermal calculation node at this moment calculated in step S4.2, the current temperature field distribution of the rock mass model is obtained, and the load caused by the temperature change is calculated; for each triangle, the grid ΔF is calculated. j F j+1 F j+2 , calculate the load f caused by temperature change according to formula (3) T :

[0028]

[0029] In the above formula, is the total temperature difference of the corresponding node compared to the beginning of the simulation, j is a positive integer, a=0,1,2; K * is the coefficient.

[0030] Step S4.4. According to the latest temperature field distribution of the rock model obtained in step S4.3, the load caused by the temperature change is substituted into formula (4) to obtain the acceleration of each force calculation node, and then the velocity and deformation of each force calculation node of the rock model are updated to perform rock mechanics simulation;

[0031] m k a k =f ext -f int -f T (4)

[0032] In the above formula, a k is the acceleration vector of the kth force calculation node, m k is the mass of the kth force calculation node, which is associated with the ith thermal calculation node; f ext is the external force on the node, including body force, external load, contact force and joint force, etc.; f int is the internal force on the node, f T is the temperature load.

[0033] The velocity and deformation calculation formulas of the kth force calculation node are as follows (5) and (6):

[0034]

[0035] In the above formula, and are the velocity vectors of the kth force calculation node of the rock mass model at time t and t+Δt respectively; and are the displacement vectors of the kth force calculation node of the calculation model at time t and t+Δt respectively.

[0036] Further optimization, in step S4.1, Through the common node F j The sum of the heat flowing into the computing node Ti from several triangular computing grids, where each triangular computing grid ΔF j F j+1 F j+2 The heat flowing into Ti is calculated using formula (4.1):

[0037]

[0038] In the above formula, q x and q y are the heat flow rates in the x and y directions respectively; L x and L y Calculate mesh ΔF for triangles j F j+1 F j+2 Edge F j F j+1 and edge F j F j+2 The x- and y-direction components of the midpoint coordinate interpolation;

[0039] q x and q y The calculation formula (4.2) is as follows:

[0040]

[0041] In the above formula, k is the thermal conductivity coefficient, which is known; and They are the x and y coordinates of the corresponding nodes respectively; is the corresponding node temperature, a=0,1,2.

[0042] Further optimization, in step S4.1, the heat flowing into the heat calculation node Ti through the crack heat exchange Equal to the common node F on the other side of the crack p The sum of the heat flowing into Ti through heat exchange in the triangular computational grid of the side; among them, the common node F p Calculate the mesh ΔF for each triangle p F p+1 F p+2 The heat flowing into Ti through the cracks for heat exchange is calculated using formula (4.3):

[0043]

[0044] In the above formula, k J is the crack heat conductivity coefficient; and Node F p and F jtemperature; L is the length of the crack; the node F j and F p When the joint unit is not inserted, they are overlapped. After the joint unit is destroyed, the two are located on both sides of the crack and correspond to each other. p is a positive integer.

[0045] Further optimization, in step S4.3, the coefficient K * The value of is related to the plane assumption. Under the plane stress assumption,

[0046]

[0047] In the above formula, E is the elastic modulus, α is the thermal expansion coefficient, and μ is the Poisson's ratio, all of which can be measured experimentally.

[0048] Under the plane strain assumption,

[0049]

[0050] Further optimization, in step S5, after each calculation step is completed, the damaged joint unit is obtained and removed according to the deformation and stress state of the joint unit, and the joint unit damage mode is as follows:

[0051] 1) Tensile failure: d n ≥d nc And d s <d s0 ;

[0052] 2) Shear failure: d s ≥d sc And d n <d n0 ;

[0053] 3) Mixed destruction: d s0 <d s <d sc And d n0 <d n <d n0 ,satisfy

[0054] Among them, d n and d s are the normal deformation and tangential deformation of the joint element respectively; the other variables satisfy the following form:

[0055] d n0 =f t / k n (5.1)

[0056] d nc =d n0 +2G I / f t(5.2)

[0057] d s0 =f s / k s (5.3)

[0058] d sc =d s0 +2G II / f s (5.4)

[0059] In the formula, f t and f s are the tensile and shear strength of the rock mass, respectively; k n and k s are the normal and tangential joint stiffnesses respectively; G I and G II are the mode I and II fracture energies, respectively, and are known.

[0060] Compared with the prior art, the present invention has the following beneficial effects:

[0061] The present invention specifically divides the calculation model into a potential damage zone and an indestructible zone, and only considers the evolution process of cracks in the potential damage zone, thereby reducing the degree of freedom of the calculation system and saving calculation time. At the same time, two different sets of thermal and mechanical calculation nodes are introduced to avoid the introduction of artificial heat transfer parameters, and the entire thermal-mechanical coupling process of the complex rock system can be simulated. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 A flow chart of the continuous-discontinuous numerical simulation method considering the thermal-solid coupling problem according to the present invention;

[0063] Figure 2 A schematic diagram of the division of the potential damage area and the indestructible area of ​​the present invention;

[0064] Figure 3 A schematic diagram of setting joint units at the grid boundary of the potential damage zone according to the present invention;

[0065] Figure 4 A schematic diagram of the dynamic update of the mapping relationship between the thermal calculation node system and the force calculation node system of the present invention;

[0066] Figure 5 This is a schematic diagram of heat transfer calculation according to the present invention;

[0067] Figure 6 This is a schematic diagram of a rock column model containing joint units according to the second embodiment of the present invention;

[0068] Figure 7 The temperature evolution curve of the measuring line obtained by simulating the method 1 in the second embodiment of the present invention;

[0069] Figure 8 The temperature evolution curve of the measuring line obtained by simulating the method 2 described in the second embodiment of the present invention;

[0070] Fig. 9 This is a schematic diagram of a rock mass model containing complex fractures according to the third embodiment of the present invention;

[0071] Fig.10 This is a schematic diagram of the evolution of temperature of a rock mass containing complex fractures over time obtained by simulation in Example 3 of the present invention;

[0072] Fig.11 This is a schematic diagram of the conglomerate formation in the fourth embodiment of the present invention;

[0073] Fig.12 This is a schematic diagram of the evolution of conglomerate formation temperature and cracks over time simulated in Example 4 of the present invention. DETAILED DESCRIPTION

[0074] In order to make the purpose, technical solution and advantages of the embodiments of the present invention clearer, the technical solution of the present invention will be clearly and completely described below in conjunction with the accompanying drawings. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0075] like Figure 1 As shown, a continuous-discontinuous numerical simulation method considering the thermal-mechanical coupling problem of rock mass includes the following steps:

[0076] Step S1. According to the stress characteristics of the rock mass, the simulation area of ​​the rock mass model is divided into a potential damage area and an indestructible area, and the potential damage area and the indestructible area are meshed.

[0077] Specifically, the stress concentration area caused by internal defects of the rock mass or external disturbance is considered and set as the potential damage area, and the rest of the area is set as the indestructible area. In the case of two-dimensional calculation, the triangular calculation grid is used to discretize the model. This is because the triangular grid calculation is the simplest and easy to implement. For details, refer to the schematic diagram. Figure 2 .

[0078] in, Figure 2 (a) is a schematic diagram of dividing the potential simulation area of ​​a circular tunnel into a potential damage area and an indestructible area; Figure 2 (c) The result of meshing the potential simulation area of ​​the circular tunnel. Figure 2 (b) is a schematic diagram of dividing the potential simulation area of ​​the arch tunnel into a potential damage area and an indestructible area; Figure 2(d) Result diagram of meshing of the potential simulation area of ​​the arch tunnel.

[0079] Step S2. Joint units are set at the grid boundaries of the potential damage zone. The joint units connect the grids on both sides thereof, and the rock mass fracture process is simulated through the destruction of the joint units.

[0080] The computational grids in the potential damage zone are independent of each other. Adjacent computational grids are connected together by inserting several joint units. The failure of the joint unit represents the initiation and expansion of the real crack. For details, refer to the schematic diagram. Figure 3 , the blue part in the figure is the joint unit, which is distributed inside and on the boundary of the potential damage area. Figure 3 (a) Schematic diagram of inserting joint elements into the potential failure zone of the circular tunnel potential model; Figure 3 (b) Schematic diagram of inserting joint elements into the potential failure zone of the arch tunnel potential model.

[0081] Step S3. Based on the model meshing and joint unit insertion results, a thermal calculation node system and a force calculation node system for simulation calculation are generated, and there is a mapping relationship between the two. The mapping relationship is continuously and dynamically updated during the simulation of rock mass fracture. Specifically, when the model is meshed in step S1, a large number of nodes are first inserted inside and on the boundary of the model, and three adjacent nodes are connected to each other to form a triangular calculation grid; all nodes are divided into two sets of node systems for thermal calculation and force calculation. The nodes in the thermal calculation system are called thermal calculation nodes, and the nodes in the force calculation system are called force calculation nodes; in the indestructible area, the thermal calculation nodes and the force calculation nodes are completely consistent, and there are multiple adjacent triangular calculation grids that share the same node; inside the potential damage area, due to the insertion of joint units, the number of force calculation nodes is greater than the number of thermal calculation nodes, and one thermal calculation node is at least associated with one or more force calculation nodes.

[0082] Before the joint is broken, the adjacent triangular calculation meshes share a node, and one thermal calculation node is associated with 4, 3 or 2 force calculation nodes; after the joint is broken, the adjacent triangular calculation meshes no longer share a node, and one thermal calculation node is associated with a corresponding calculation node.

[0083] The above node system settings can refer to the schematic diagram Figure 4(a), there are 4 triangular computational grids in the potential damage zone, and joint units are inserted at the grid interface. For the convenience of display, the joint units are enlarged in the figure. Before the simulation starts, the joint units actually have no thickness. All vertices of the 4 triangular units are numbered in sequence, and a total of 12 mechanical calculation nodes from F1 to F12 are obtained. Among them, the coordinates of F1, F4, F7 and F10 are the same. A thermal calculation node T1 for heat conduction calculation is set at this coordinate. T1 is associated with the four mechanical calculation nodes F1, F4, F7 and F10. The mapping relationship can be expressed as a counterclockwise linked list, T1: {F1→F4→F7→F10→F1}, which can be used for subsequent heat conduction calculations, and the temperature of each point on the linked list is the same, satisfying the following form:

[0084] T T1 = T F1 = T F4 = T F7 = T F10 (1)

[0085] Step S4. Decompose the simulation process of rock mass fracture into several calculation steps, and solve the heat conduction process and mechanical calculation process in each calculation step. Specifically, it includes the following steps:

[0086] Step S4.1. Calculate the heat flowing into each heat calculation node at this moment, specifically:

[0087] 1) When the joint unit is not destroyed, the crack does not exist, such as Figure 4 (a), the node corresponding to the hot calculation node Ti is shared by four triangular calculation grids, then Equal to the heat flowing into the Ti node through the 4 triangles that share this node. sum.

[0088] 2) When the joint unit is destroyed and cracks exist, the cracks divide the triangular computational grid of the four common nodes into two parts, one on each side of the crack, such as Figure 4 (c) shows that Calculated by formula (2):

[0089]

[0090] In the above formula, The heat flowing into Ti for the triangles located on the same side of the thermal calculation node Ti is calculated. The heat flowing into node Ti through the crack is calculated for the triangle mesh on the other side of the crack.

[0091] Reference Figure 5 , further explains the calculation process of heat conduction, where Figure 5(a) is an example diagram of heat flow into F1 in the triangular computational grid ΔF1F2F3; Figure 5 (b) Schematic diagram of heat inflow through cracks for heat exchange. Assume that the joint units connecting ΔF1F2F3 and ΔF10F11F12 and connecting ΔF4F5F6 and ΔF7F8F9 are damaged to form cracks. Take the thermal node T1 as an example. There are two triangular computational grids connected around it, and the mapping relationship is T1:{F1→F4}. Based on the mapping relationship represented by the linked list, the heat flowing into the T1 node through the two triangular computational grids can be calculated. The calculation formula is as follows:

[0092]

[0093] In the formula, is the heat flowing into T1 through the computational grid; Q ΔF1F2F3→F1 is the heat flowing into F1 through ΔF1F2F3; Q ΔF4F5F6→F4 is the heat that flows into F4 through ΔF4F5F6.

[0094] Q ΔF1F2F3→F1 For example, the calculation formula is as follows:

[0095] Q ΔF1F2F3→F1 =q x L y +q y L x (4)

[0096] In the formula, q x and q y are the heat flow rates in the x and y directions respectively; L x and L y The x- and y-direction components of the interpolated coordinates of the midpoints of sides F1F2 and F1F3. x and q y The calculation formula is as follows:

[0097]

[0098] Where k is the thermal conductivity coefficient; and They are the x and y coordinates of the corresponding nodes; is the corresponding node temperature, a=1,2,3.

[0099] At the same time, the angular computational grids on both sides of the crack will also exchange heat. The heat flowing into the T1 node through the crack is calculated according to the following formula:

[0100]

[0101] In the above formula, is the heat flowing into T1 through the crack; QΔF10F11F12→F1 is the heat flowing into F1 through ΔF10F11F12; Q ΔF7F8F9→F4 is the heat that flows into F4 through ΔF7F8F9.

[0102] Q ΔF10F11F12→F1 For example, the calculation formula is as follows:

[0103] Q ΔF10F11F12→F1 =0.5*k J (T F10 -T F1 )L (7)

[0104] In the formula, k J is the crack heat conductivity; T F1 and T F10 are the temperatures of nodes F1 and F10 respectively; L is the length of the crack.

[0105] Step S4.2. Calculated according to step S4.1 Substitute into the rock mass heat conduction formula (8) and calculate the temperature of each heat calculation node at this time;

[0106]

[0107] In the above formula, and are the temperatures of the i-th thermal calculation node Ti in the model at time t and t+Δt respectively; is the heat flowing into the i-th node; c p is the specific heat capacity of the rock mass; M i is the node mass corresponding to the i-th thermal calculation node, which is equal to the sum of the masses of all associated force calculation nodes.

[0108] Step S4.3. According to the temperature of each thermal calculation node at this moment calculated in step S4.2, the current temperature field distribution of the rock mass model is obtained, and the load caused by the temperature change is calculated; for each triangle, the grid ΔF is calculated. j F j+1 F j+2 , calculate the load f caused by temperature change according to formula (9) T :

[0109]

[0110] In the above formula, is the total temperature difference of the corresponding node compared to the beginning of the simulation, j is a positive integer, a=0,1,2; K * is the coefficient.

[0111] Step S4.4. According to the latest temperature field distribution of the rock model obtained in step S4.3, the load caused by the temperature change is substituted into formula (10) to obtain the acceleration of each force calculation node, and then the velocity and deformation of each force calculation node of the rock model are updated to perform rock mechanics simulation;

[0112] m k a k =f ext -f int -f T (10)

[0113] In the above formula, a k is the acceleration vector of the kth force calculation node, m k is the mass of the kth force calculation node, which is associated with the ith thermal calculation node; f ext is the external force on the node, including body force, external load, contact force and joint force, etc.; f int is the internal force on the node, f T is the temperature load.

[0114] The velocity and deformation calculation formulas of the kth force calculation node are as follows (11) and (12):

[0115]

[0116] In the above formula, and are the velocity vectors of the kth force calculation node of the rock mass model at time t and t+Δt respectively; and are the displacement vectors of the kth force calculation node of the calculation model at time t and t+Δt respectively.

[0117] Step S5. After each calculation step is completed, determine whether the joint unit is damaged based on the deformation and stress state of the joint unit, and update the mapping relationship between the thermal calculation node system and the force calculation node system. The mode of joint unit damage is as follows:

[0118] 1) Tensile failure: d n ≥d nc And d s <d s0 ;

[0119] 2) Shear failure: d s ≥d sc And d n <d n0 ;

[0120] 3) Mixed destruction: d s0 <d s <d sc And dn0 <d n <d n0 ,satisfy

[0121] Among them, d n and d s are the normal deformation and tangential deformation of the joint element respectively; the other variables satisfy the following form:

[0122] d n0 =f t / k n (13)

[0123] d nc =d n0 +2G I / f t (14)

[0124] d s0 =f s / k s (15)

[0125] d sc =d s0 +2G II / f s (16)

[0126] In the formula, f t and f s are respectively tensile strength and shear strength, which can be measured by experiment; k n and k s are the normal and tangential joint stiffnesses, respectively, which are microscopic calculation parameters and are usually taken as multiples of the elastic modulus; G I and G II are mode I and mode II fracture energies, respectively.

[0127] After the joint unit is destroyed, the mapping relationship between the thermal node system and the force node system needs to be updated. Figure 4 (b)-(c). When the joint unit connecting ΔF4F5F6 and ΔF7F8F9 is damaged, the linked list (T1:{F1→F4→F7→F10→F1}) storing the mapping relationship needs to be updated. The updated linked list structure still needs to satisfy the counterclockwise order, that is, T1:{F7→F10→F1→F4}. As the crack further expands, the joint unit connecting ΔF1F2F3 and ΔF10F11F12 is damaged. At this time, there is no joint unit connection between the triangular meshes on the upper and lower sides. The model is divided into two parts. A new thermal calculation node T2 is required. T1 and T2 are respectively associated with the force calculation nodes on the upper and lower sides. The original T1 linked list is updated to T1:{F1→F4}, and the newly generated T2 linked list also satisfies the counterclockwise order, which is T2:{F7→F10}.

[0128] Step S6. Repeat the above step S5 until the last calculation step is completed, and the thermal cracking process simulation of the rock mass is completed. The following is an explanation in conjunction with a specific application example:

[0129] Embodiment 1:

[0130] In this embodiment, taking tunnel engineering as an example, strategies are applied to the division of potential damage areas and indestructible areas and joint units. Figure 2 As shown in the figure, in order to avoid the influence of boundary effect on simulation accuracy, the rock modeling area usually needs to be much larger than the tunnel size. However, actual engineering mainly focuses on the surrounding rock near the tunnel boundary, and the surrounding rock area that can be affected by external disturbance is limited. Referring to the research of prior art, the disturbance area of ​​a circular tunnel usually does not exceed 3 times the diameter of the tunnel (see the paper "X.Du, P.Zhang, L.Jin, D.Lu, 2019, A multi-scale analysis method for the simulation of tunnel excavation in sandy cobble stratum, Tunnelling and Underground Space Technology"). Therefore, when setting the potential damage zone, its size can be set to a value slightly larger than the range of the disturbance zone. Similarly, for tunnels of any shape, the multiple of its maximum size is taken as the size of the potential damage zone. After discretizing the model using a triangular computational grid, inserting joint units inside and on the boundary of the potential damage zone, the complete computational model is obtained.

[0131] Embodiment 2:

[0132] In this embodiment, the simulation of thermal conduction of rock columns containing joint units is taken as an example to divide the potential damage area and the indestructible area and apply strategies to the joint units. The generation strategy of the thermal calculation node system and the force calculation node system in the present invention is explained, and the advantages of the method are reflected. Figure 6 (a) shows a rectangular rock column, 12m long and 1m wide. The model is discretized using triangular meshes to obtain 386 triangular meshes, as shown in Figure 6 As shown in (b). Before the joint unit is inserted, there are common edges between adjacent triangular meshes, and the nodes on the edges are shared by the triangular meshes on both sides, with a total of 246 nodes. The 246 nodes are set as thermal calculation nodes, and then the joint units are inserted into the model. At this time, the nodes of all triangular meshes are independent of each other, and the total number of points is three times the number of meshes, that is, 386*3=1158 nodes, which are force calculation nodes. The force calculation nodes with the same coordinates belong to the same thermal calculation node, constituting the mapping relationship described in step S3 of the invention content, as shown in FIG. Figure 6 (c) as shown.

[0133] For comparison, two methods are used to simulate the heat transfer of this model: the first method is the method proposed by Fang Ming, referred to as Method 1; the second method introduces a joint heat transfer model, which is the most common heat conduction model in the current discontinuous method and hybrid method. This model assumes that heat is exchanged between adjacent units through joint units. By introducing an artificially assumed joint heat transfer coefficient, the heat conduction simulation between models can be realized without cracking. This method is hereinafter referred to as Method 2. According to the research of the paper "C. Yan, Y. Jiao, 2020, A 2D discrete heat transfer model considering the thermal resistance effect of fractures for simulating the thermal cracking of brittle materials, Acta Geotechnica", the value of the joint heat transfer coefficient should be much larger than the true heat transfer coefficient of the material, and its calculation formula is as follows:

[0134] k J =nk / L e (17)

[0135] Where, L e is the joint length, the coefficient n should be greater than 100, and is taken as 100 here.

[0136] The initial temperature of the rock column is 0°C, and a constant temperature heat source is applied to its bottom to maintain its temperature at T B = 100°C, while the top temperature is fixed at T T =0℃. A horizontal measuring line is arranged in the middle of the rock column to monitor its temperature change. The temperature at any point in the model has an analytical solution, which satisfies the following formula:

[0137]

[0138] κ= k / ρc p (19)

[0139] Where y is the distance from the point to the bottom boundary; L is the length of the model; thermal conductivity k = 10 W / (m·℃); density ρ = 2000 kg / m 3 Specific heat capacity c p =1J / (kg·℃).

[0140] Under different calculation steps Δt, the temperature evolution curves of the measuring line simulated by method 1 and method 2 are as follows: Figure 7 and Figure 8 As shown. Among them, Figure 7(a), 7(b), 7(c), and 7(d) represent the temperature evolution curves of the measuring line simulated using method 1 when Δt is 4s, 2.5s, 2s, and 1s, respectively. Figure 8 (a) and (b) respectively show the results of method 2 when Δt is 2×10 -2 s and 1×10 -2 s, the simulated temperature evolution curve of the measuring line.

[0141] In method 1, when the value of Δt is less than 2.5s, Figure 7 As shown in (b), the simulated temperature evolution law of the measuring line is completely consistent with the analytical solution. In method 2, the requirement for the calculation step size is very strict: when Δt = 2×10 -2 s, the simulation results are completely unreliable, such as Figure 8 (a). Further reducing the value of the calculation step, when Δt = 1×10 -2 s, the simulation results are consistent with the analytical solution, such as Figure 8 (b) is shown. It is worth noting that, under the same simulation time, the value of the step size is inversely proportional to the total number of calculation steps. A smaller calculation step size means that more calculation steps are required. Taking this embodiment as an example, the number of calculation steps of method 1 and method 2 can differ by 250 times. The calculation time of method 1 is only 4s, while the calculation time of method 2 is 406s.

[0142] Embodiment 3:

[0143] In this embodiment, the numerical simulation method of the present invention is described by taking the simulation of heat transfer in fractured rock mass as an example and combining with specific examples. Fig. 9 The figure shows a square rock sample with 9 natural cracks inside. The initial temperature of the rock is 0℃, and a constant temperature heat source is applied to the top of the rock to maintain its temperature at 100℃. Fig.10 To simulate the temperature evolution process, Fig.10 (a)-(d), it can be seen that a heat flow from top to bottom is formed due to the temperature difference. The temperature of the bottom rock mass gradually increases with the increase of heat transfer time. Since the thermal conductivity of the crack is weaker than that of the intact rock mass, there is an obvious difference in temperature on both sides of the crack, which slows down the speed of heat flow.

[0144] Embodiment 4:

[0145] In this embodiment, the thermal fracture simulation of geothermal mining in conglomerate strata is taken as an example, and the numerical simulation method of the present invention is explained in combination with specific examples for the division of potential damage areas and indestructible areas and the application strategy of joint units. Fig.11As shown in the figure, it is a deeply buried geothermal reservoir with conglomerate as the lithology. In order to extract geothermal energy, an energy transmission pipeline is laid in the rock layer, and heat exchange with the surrounding rock is achieved by injecting cold water to bring geothermal energy back to the surface. Assume that the initial temperature of the surrounding rock is 200℃ and the initial water temperature in the pipeline is 20℃. Conglomerate consists of two parts: block stone and matrix. Considering that the strength of block stone is usually much greater than that of matrix structure and it is difficult to damage, the block stone is set as an indestructible zone and the matrix is ​​set as a potential damage zone in the modeling process to simulate the thermal cracking phenomenon that may occur during the heat exchange process.

[0146] Fig.12 In order to simulate the crack growth and temperature evolution process, Fig.12 (a)-(d), as geothermal energy is extracted, the temperature of the surrounding rock near the pipeline gradually decreases, generating additional temperature stress. When the temperature stress reaches the strength limit, the surrounding rock near the pipeline is first damaged, and the cracks extend away from the pipeline. Due to the initiation and expansion of the cracks, the temperature field inside the rock mass becomes discontinuous, and the temperature on both sides of the crack jumps.

[0147] The above descriptions are merely embodiments of the present invention and are not intended to limit the patent scope of the present invention. Any equivalent structure or equivalent process transformation made using the contents of the present invention specification and drawings, or directly or indirectly applied in other related technical fields, are also included in the patent protection scope of the present invention.

Claims

1. A continuous-discontinuous numerical simulation method considering the thermal-mechanical coupling problem of rock mass, characterized in that: The following steps are involved: Step S1. According to the stress characteristics of the rock mass, the simulation area of ​​the rock mass model is divided into a potential damage area and an indestructible area, and the potential damage area and the indestructible area are meshed; Step S2. setting a joint unit at each grid boundary of the potential damage zone, the joint unit connects the grids on both sides thereof, and simulates the rock mass fracture process through the destruction of the joint unit; Step S3. Based on the model mesh division and joint unit insertion results, a thermal calculation node system and a force calculation node system for simulation calculation are generated, and there is a mapping relationship between the two. The mapping relationship is continuously and dynamically updated during the process of simulating rock mass fracture; Step S4. Decompose the simulation process of rock mass fracture into several calculation steps, and solve the heat conduction process and the mechanical calculation process in each calculation step in turn; Step S5. After each calculation step is completed, determine whether the joint unit is damaged based on its deformation and stress state, and update the mapping relationship between the thermal calculation node system and the force calculation node system; Step S6. Repeat the above steps S4 and S5 until the last calculation step is completed, and the thermal cracking process simulation of the rock mass is completed.

2. The continuous-discontinuous numerical simulation method considering the thermal-mechanical coupling problem of rock mass according to claim 1 is characterized in that: In the step S1, the stress concentration area of ​​the rock mass caused by its internal defects or external disturbance is set as a potential damage area, and the remaining areas are set as indestructible areas; For the two-dimensional rock mass model, triangular computational grid is used to discretize the model.

3. The continuous-discontinuous numerical simulation method considering the thermal-mechanical coupling problem of rock mass according to claim 2 is characterized in that: In step S2, the computational grids in the potential failure zone are independent of each other, and adjacent computational grids are connected together by inserting joint units. The failure of the joint units indicates the initiation and expansion of real cracks in the rock mass.

4. The continuous-discontinuous numerical simulation method considering the thermal-mechanical coupling problem of rock mass according to claim 3 is characterized in that: In step S3, based on the model mesh division and joint unit insertion results, a thermal calculation node system and a force calculation node system for simulation calculation are generated, specifically: When meshing the model in step S1, a large number of nodes are first inserted inside and on the boundary of the model, and three adjacent nodes are connected to each other to form a triangular computational mesh; All nodes are divided into two sets of node systems: thermal calculation and force calculation. The nodes in the thermal calculation system are called thermal calculation nodes, and the nodes in the force calculation system are called force calculation nodes. In the indestructible zone, the thermal calculation nodes and the force calculation nodes are exactly the same, and there are multiple adjacent triangular calculation meshes that share the same node; inside the potential damage zone, due to the insertion of joint units, the number of force calculation nodes is greater than the number of thermal calculation nodes, and one thermal calculation node is associated with at least one or more force calculation nodes.

5. The continuous-discontinuous numerical simulation method considering the thermal-mechanical coupling problem of rock mass according to claim 4 is characterized in that: In the step S4, solving the heat conduction process and the mechanical calculation process in each calculation step specifically includes the following steps: Step S4.

1. Calculate the heat flowing into each heat calculation node at this moment, specifically: 1) When the joint unit is not destroyed, the crack does not exist, and the node corresponding to the thermal calculation node Ti is shared by multiple triangular calculation grids. This is equivalent to calculating the heat flowing into Ti through multiple triangles sharing the node. sum; 2) When the joint unit is destroyed and there is a crack, the crack divides the multiple triangular computational meshes that share the node into two parts, one on each side of the crack. Calculated by formula (1): In the above formula, The heat flowing into Ti for the triangles located on the same side of the thermal calculation node Ti is calculated. The heat that flows into node Ti through the crack through heat exchange for the triangle calculation mesh on the other side of the crack; Step S4.

2. Calculated according to step S4.1 Substitute into the rock mass heat conduction formula (2) and calculate the temperature of each heat calculation node at this time; In the above formula, and are the temperatures of the i-th thermal calculation node Ti in the model at time t and t+Δt respectively; is the heat flowing into the i-th node; c p is the specific heat capacity of the rock mass; M i is the node mass corresponding to the i-th thermal calculation node, which is equal to the sum of the masses of all associated force calculation nodes; Step S4.

3. According to the temperature of each thermal calculation node at this moment calculated in step S4.2, obtain the current temperature field distribution of the rock mass model, and calculate the load caused by the temperature change; Calculate the mesh ΔF for each triangle j F j+1 F j+2 , calculate the load f caused by temperature change according to formula (3) T : In the above formula, is the total temperature difference of the corresponding node compared to the beginning of the simulation, j is a positive integer, a=0, 1, 2; K* is the coefficient; Step S4.

4. According to the latest temperature field distribution of the rock model obtained in step S4.3, the load caused by the temperature change is substituted into formula (4) to obtain the acceleration of each force calculation node, and then the velocity and deformation of each force calculation node of the rock model are updated to perform rock mechanics simulation; m k a k =f ext -f int -f T (4) In the above formula, a k is the acceleration vector of the kth force calculation node, m k is the mass of the kth force calculation node, which is associated with the ith thermal calculation node; f ext is the external force on the node; f int is the internal force on the node, f T is the temperature load; The velocity and deformation calculation formulas of the kth force calculation node are as follows (5) and (6): In the above formula, and are the velocity vectors of the kth force calculation node of the rock mass model at time t and t+Δt respectively; and are the displacement vectors of the kth force calculation node of the calculation model at time t and t+Δt respectively.

6. The continuous-discontinuous numerical simulation method considering the thermal-mechanical coupling problem of rock mass according to claim 5, characterized in that: In the step S4.1, Through the common node F j The sum of the heat flowing into the computing node Ti from several triangular computing grids, where each triangular computing grid ΔF j F j+1 F j+2 The heat flowing into Ti is calculated using formula (4.1): In the above formula, q x and q y are the heat flow rates in the x and y directions respectively; L x and L y Calculate mesh ΔF for triangles j F j+1 F j+2 Edge F j F j+1 and edge F j F j+2 The x- and y-direction components of the midpoint coordinate interpolation; q x and q y The calculation formula (4.2) is as follows: In the above formula, k is the thermal conductivity coefficient, which is known; and They are the x and y coordinates of the corresponding nodes respectively; is the corresponding node temperature, a=0,1,2.

7. The continuous-discontinuous numerical simulation method considering the thermal-mechanical coupling problem of rock mass according to claim 6, characterized in that: In step S4.1, the heat flowing into the heat calculation node Ti through the crack is Equal to the common node F on the other side of the crack p The sum of the heat flowing into Ti through heat exchange in the triangular computational grid of the side; among them, the common node F p Calculate the mesh ΔF for each triangle p F p+1 F p+2 The heat flowing into Ti through the cracks for heat exchange is calculated using formula (4.3): In the above formula, k J is the crack heat conductivity; and Node F p and F j temperature; L is the length of the crack; the node F j and F p They overlap when no joint unit is inserted, and p is a positive integer.

8. The continuous-discontinuous numerical simulation method considering the thermal-mechanical coupling problem of rock mass according to claim 7, characterized in that: In step S4.3, the coefficient K * The value of is related to the plane assumption. Under the plane stress assumption, In the above formula, E is the elastic modulus, α is the thermal expansion coefficient, and μ is the Poisson's ratio, which can all be measured experimentally; Under the plane strain assumption, 9. The continuous-discontinuous numerical simulation method considering the thermal-mechanical coupling problem of rock mass according to claim 8, characterized in that: In step S5, after each calculation step is completed, the damaged joint units are obtained and removed according to the deformation and stress state of the joint units. The joint unit damage mode is as follows: 1) Tensile failure: d n ≥d nc And d s <d s0 ; 2) Shear failure: d s ≥d sc And d n <d n0 ; 3) Mixed destruction: d s0 <d s <d sc And d n0 <d n <d n0 ,satisfy Among them, d n and d s are the normal deformation and tangential deformation of the joint element respectively; the other variables satisfy the following form: d n0 =f t / k n (5.1) d nc =d n0 +2G I / f t (5.2) d s0 =f s / k s (5.3) d sc =d s0 +2G II / f s (5.4) In the formula, f t and f s are the tensile and shear strengths of the rock mass, respectively; k n and k s are the normal and tangential joint stiffnesses respectively; G I and G II are mode I and mode II fracture energies, respectively.

Images