A continuous-discontinuous numerical simulation method considering thermal coupling of rock mass

Abstract

The application discloses a continuous-noncontinuous numerical simulation method considering rock mass thermal coupling problems, divides a simulation area into a potential damage area and an undamageable area, and carries out grid division on the simulation area; then generates a thermal calculation node system and a force calculation node system for simulation calculation, and mapping relationship exists between the two; a joint element is arranged on the grid boundary of the potential damage area, and the rock mass fracture process is simulated through the damage of the joint element; the simulation process is decomposed into a plurality of calculation steps, the thermal conduction process and the mechanical calculation process are sequentially solved in each calculation step, whether the joint element is damaged is judged according to the deformation and stress state of the joint element, and the mapping relationship is updated, until the last calculation step is finished, and the simulation of the thermal fracture process of the rock mass is completed. The application can avoid introducing artificial heat transfer parameters in the thermal conduction simulation process, can pertinently study the fracture evolution law of the local area of the rock mass, can reduce the calculation cost of the model, and can effectively simulate the dynamic fracture process in the rock mass thermal coupling problems.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of geotechnical engineering numerical simulation, and particularly relates to a continuous-discontinuous numerical simulation method considering rock mass thermal-mechanical coupling problems. BACKGROUND

[0002] Thermal stress caused by temperature change can lead to the initiation and propagation of cracks in rock mass, further causing the instability and damage of engineering facilities. Although physical tests can intuitively show the thermal cracking process of rock, 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 damage behavior of deep geotechnical facilities, which can overcome the difficulties encountered in laboratory experiments. However, to ensure the accuracy of numerical simulation, the coverage area of the calculation model is usually much larger than the area of interest in engineering design. To simplify the modeling cost and facilitate program writing, existing thermal-mechanical 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 by temperature stresses, Computers and Geotechnics” and “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”, 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, to simplify the heat conduction calculation process, the additional introduction of artificial heat transfer parameters leads to the need for additional calibration of the algorithm when simulating rock mass heat transfer problems. SUMMARY

[0004] The purpose of the present application is to provide a continuous-discontinuous numerical simulation method considering rock mass thermal-mechanical coupling problems, which solves the technical problems in the prior art that the study area is usually not processed additionally when simulating the thermal cracking process, resulting in unnecessary computational consumption; at the same time, to simplify the heat conduction calculation process, the additional introduction of artificial heat transfer parameters leads to the need for additional calibration of the algorithm when simulating rock mass heat transfer problems.

[0005] To solve the above technical problems, the application adopts the following technical solutions to achieve the purpose:

[0006] A continuous-discontinuous numerical simulation method considering rock mass thermal coupling problems, comprising 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 zone and an undamageable zone, and the potential damage zone and the undamageable zone are meshed.

[0008] Step S2. A joint element is set at the grid boundary of the potential damage zone, the joint element connects the grids on both sides thereof, and the rock mass fracture process is simulated through the destruction of the joint element.

[0009] Step S3. Based on the model meshing and joint element insertion results, a thermal calculation node system and a force calculation node system for simulation calculation are generated, and the two systems have a mapping relationship, which is dynamically updated during the simulation of the rock mass fracture process.

[0010] Step S4. The simulation process of the rock mass fracture is divided into several calculation steps, and the heat conduction process and the mechanical calculation process are sequentially solved in each calculation step.

[0011] Step S5. After each calculation step is completed, it is judged whether the joint element is damaged according to the deformation and stress state of the joint element, and the mapping relationship between the thermal calculation node system and the force calculation node system is updated.

[0012] Step S6. Repeat the above step S5 to the end of the last calculation step, and then the simulation of the thermal fracture process of the rock mass is completed.

[0013] Further optimization, in the step S1, the stress concentration area caused by the internal defects or external disturbances of the rock mass is set as the potential damage zone, and the remaining area is set as the undamageable zone; for a two-dimensional rock mass model, triangular calculation grids are used to discretely divide the model.

[0014] Further optimization, in the step S2, the calculation grids of the potential damage zone are independent of each other, and adjacent calculation grids are connected together by inserting joint elements, and the destruction of the joint element represents the initiation and expansion of the real rock mass crack.

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

[0016] In the step S1 of meshing the model, a large number of nodes are inserted in the model and on the boundary, and three adjacent nodes are connected to form a triangular calculation grid.

[0017] All nodes are divided into two sets of nodes 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 non-destructive zone, the thermal calculation nodes and the force calculation nodes are completely consistent, and there are multiple adjacent triangular calculation grids sharing the same node. In the potential damage zone, because joint elements are inserted, the number of force calculation nodes is greater than the number of thermal calculation nodes, and one thermal calculation node is associated with one or more force calculation nodes.

[0018] Further optimization, in step S4, the heat conduction process and the mechanical calculation process are solved in each calculation step, specifically including the following steps.

[0019] Step S4.1. Calculate the inflow heat of each thermal calculation node at this moment, specifically:

[0020] 1) When the joint element is not damaged, the crack does not exist, the node corresponding to the thermal calculation node Ti is shared by multiple triangular calculation grids, then is equal to the heat flowing into Ti through the multiple triangular calculation grids sharing the node Total;

[0021] 2) When the joint element is damaged, the crack exists, the crack divides the multiple triangular calculation grids sharing the node into two parts, which are located on the two sides of the crack, then is calculated by formula (1):

[0022]

[0023] In the above formula, is the heat flowing into Ti through the triangular calculation grid on the same side of the thermal calculation node Ti, is the heat flowing into node Ti through the crack by heat exchange of the triangular calculation grid on the other side of the crack.

[0024] Step S4.2. According to the calculated in step S4.1, the temperature of each thermal calculation node at this moment is calculated by bringing the rock mass heat conduction formula (2) into the model;

[0025]

[0026] In the above formula, and are the temperatures of the i-th thermal calculation node Ti in the model at 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 mass of the i-th thermal calculation node corresponding to the node, which is equal to the sum of the masses of all associated force calculation nodes.

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

[0028]

[0029] In the above formula, is the total temperature difference of the corresponding node compared with 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 mass model obtained in step S4.3, the load caused by temperature change is brought 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 mass 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, the kth force calculation node is associated with the ith thermal calculation node; f ext is the external force received by the node, including body force, external load, contact force and joint force, etc.; f int is the internal force received by the node, f T is the temperature load.

[0033] The velocity and deformation calculation formula of the kth force calculation node is as follows:

[0034]

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

[0036] Further optimization, in step S4.1, The sum of the heat flowing into the calculation node Ti is calculated by several triangular calculation meshes of the shared node F j F j F j+1 F j+2 The heat flowing into Ti is calculated by 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 are the x and y direction components of the midpoint coordinate interpolation of the triangular calculation mesh ΔF j F j+1 F j+2 F j F j+1 F j F j+2 are the x and y direction components of the midpoint coordinate interpolation of the edge F x F y

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

[0040]

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

[0042] Further optimization, in the step S4.1, the heat flowing into the heat calculation node Ti through heat exchange by the crack is equal to the sum of the heat flowing into Ti through heat exchange by the triangular calculation meshes of the shared node F p located on the other side of the crack; wherein each triangular calculation mesh ΔF p F p F p+1 F p+2 The heat flowing into Ti through heat exchange by the crack is calculated by formula (4.3):

[0043]

[0044] In the above formula, k J is the thermal conductivity coefficient of the crack; and are the nodes F p and F jtemperature; L is the length of the fracture; the node F j and F p coincide when the joint element is not inserted, and are located on both sides of the fracture after the joint element is destroyed, and correspond. p is a positive integer.

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

[0046]

[0047] In the above formula, E is the elastic modulus, a is the thermal expansion coefficient, and μ is the Poisson's ratio, which can be measured by experiment.

[0048] Under the plane strain assumption,

[0049]

[0050] Further optimization, in the step S5, after each calculation step is completed, the damaged joint element is obtained and removed according to the deformation and stress state of the joint element, and the mode of joint element destruction 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 failure: d s0 < d s < d sc and d n0 < d n < d n0 , which satisfies

[0054] where d n and d s are the normal deformation and shear deformation of the joint element, respectively; the remaining 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 stiffness respectively; G I and G II are the I and II type fracture energy respectively, which are known.

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

[0061] The present application divides the calculation model into potential damage area and non-damage area, only considers the crack evolution process in the potential damage area, reduces the degree of freedom of the calculation system, and saves the calculation time; meanwhile, two different heat and force calculation nodes are introduced, which avoids the introduction of artificial heat transfer parameters, and can simulate the whole process of thermal and force coupling of complex rock mass system. BRIEF DESCRIPTION OF DRAWINGS

[0062] Figure 1 It is a flow chart of the continuous-discontinuous numerical simulation method considering thermal and solid coupling problems according to the present application;

[0063] Figure 2 It is a schematic diagram of the division of the potential damage area and the non-damage area according to the present application;

[0064] Figure 3 It is a schematic diagram of the setting of joint elements on the grid boundary of the potential damage area according to the present application;

[0065] Figure 4 It is a schematic diagram of the dynamic updating of the mapping relationship between the heat calculation node system and the force calculation node system according to the present application;

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

[0067] Figure 6 It is a schematic diagram of the rock column model containing joint elements according to the second embodiment of the present application;

[0068] Figure 7 It is a temperature evolution curve diagram of the measurement line simulated by the method 1 in the second embodiment of the present application;

[0069] Figure 8 The temperature evolution curve diagram of the measuring line simulated by the method 2 in the embodiment two of the present application is shown in the figure;

[0070] Figure 9 The model diagram of the complex fissure rock mass in the embodiment three of the present application is shown in the figure;

[0071] Figure 10 The diagram of the temperature evolution of the complex fissure rock mass simulated in the embodiment three of the present application is shown in the figure;

[0072] Figure 11 The gravel stratum diagram in the embodiment four of the present application is shown in the figure;

[0073] Figure 12 The diagram of the temperature and crack evolution of the gravel stratum simulated in the embodiment four of the present application is shown in the figure. DETAILED DESCRIPTION

[0074] In order to make the purpose, technical scheme and advantages of the embodiments of the present application more clear, the technical scheme of the present application will be described clearly and completely below in combination with the drawings. Obviously, the described embodiments are some embodiments of the present application, but not all the embodiments. Based on the embodiments in the present application, all the other embodiments obtained by those skilled in the art without creative labor fall within the protection scope of the present application.

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

[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 zone and an undamaged zone, and the potential damage zone and the undamaged zone are meshed.

[0077] Specifically, the stress concentration area caused by the internal defects of the rock mass or external disturbance is set as the potential damage zone, and the remaining area is set as the undamaged zone. In the two-dimensional calculation, triangular calculation grid is used to discretize the model, which is the most simple calculation grid and is convenient for program implementation. For reference, see the diagram Figure 2 .

[0078] Among them, Figure 2 (a) is a diagram of dividing the potential simulation area of the circular tunnel into a potential damage zone and an undamaged zone; Figure 2 (c) is a diagram of the grid division result of the potential simulation area of the circular tunnel. Among them, Figure 2 (b) is a diagram of dividing the potential simulation area of the arch tunnel into a potential damage zone and an undamaged zone; Figure 2 ​(d) the result map of meshing the potential failure zone of the arch-shaped tunnel.

[0079] Step S2. Set joint elements at the grid boundaries of the potential failure zone, the joint elements connecting the grids on both sides thereof, and the failure of the joint elements simulating the rock mass cracking process.

[0080] The calculation grids of the potential failure zone are independent of each other, and adjacent calculation grids are connected by inserting a plurality of joint elements, and the failure of the joint elements represents the initiation and propagation of the real cracks, which can be referred to the schematic Figure 3 , and the blue part in the figure is the joint element, which is distributed inside and on the boundary of the potential failure zone. Among them, Figure 3 (a) a schematic diagram of inserting joint elements into the potential failure zone of the potential model of the circular tunnel; Figure 3 (b) a schematic diagram of inserting joint elements into the potential failure zone of the potential model of the arch-shaped tunnel.

[0081] Step S3. Based on the results of the model meshing and the joint element insertion, generate a thermal calculation node system and a force calculation node system for simulation calculation, and the two systems have a mapping relationship, which is dynamically updated during the simulation of the rock mass cracking process. Specifically, when the model is meshed in step S1, a large number of nodes are inserted inside and on the boundary of the model, and three adjacent nodes are connected to form a triangular calculation grid; all nodes are divided into two 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 non-failure zone, the thermal calculation nodes and the force calculation nodes are completely consistent, and there are multiple triangular calculation grids adjacent to each other sharing the same node; inside the potential failure zone, due to the insertion of joint elements, the number of force calculation nodes is greater than that of thermal calculation nodes, and one thermal calculation node is associated with one or more force calculation nodes.

[0082] When the joint is not broken, the adjacent triangular calculation grids share one node, and one thermal calculation node is associated with 4 or 3 or 2 force calculation nodes; when the joint is broken, the adjacent triangular calculation grids no longer share one node, and one thermal calculation node is associated with one corresponding calculation node.

[0083] The above node system setting can be referred to the schematic Figure 4(a), four triangular calculation grids exist in the potential damage zone, joint elements are inserted at the grid interfaces, for the convenience of display, the joint elements are enlarged in the figure, before the simulation starts, the joint elements actually have no thickness. All the vertices of the four triangular elements are numbered in turn to obtain 12 mechanical calculation nodes F1 to F12, wherein the coordinates of F1, F4, F7 and F10 are the same, a thermal calculation node T1 is set at the coordinates for thermal conduction calculation, T1 is associated with the four mechanical calculation nodes F1, F4, F7 and F10, the mapping relationship can be expressed as a counterclockwise ordered linked list T1:{F1→F4→F7→F10→F1}, the linked list can be used for subsequent thermal conduction calculation, and the temperatures of the points on the linked list are the same, which satisfies the following form:

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

[0085] Step S4. The simulation process of rock mass fracture is divided into several calculation steps, and the thermal conduction process and the mechanical calculation process are solved in each calculation step. Specifically, the following steps are included:

[0086] Step S4.1. Calculate the inflow heat of each thermal calculation node at this moment, specifically:

[0087] 1) When the joint element is not damaged, the fracture does not exist, such as Figure 4 (a), the node corresponding to the thermal calculation node Ti is shared by four triangular calculation grids, then is equal to the heat flowing into the Ti node through the four triangular calculation grids sharing the node .

[0088] 2) When the joint element is damaged, the fracture exists, the fracture divides the four shared node triangular calculation grids into two parts, which are located on the two sides of the fracture, as shown in Figure 4 (c), then is calculated by formula (2):

[0089]

[0090] In the above formula, is the heat flowing into Ti from the triangular calculation grids on the same side of the thermal calculation node Ti, is the heat flowing into the node Ti through the fracture from the triangular calculation grids on the other side of the fracture.

[0091] Referring to the schematic Figure 5 , the calculation process of thermal conduction is further illustrated, wherein Figure 5(a) An example diagram of heat flow into F1 in the triangular mesh ΔF1F2F3; Figure 5 (b) A schematic diagram of heat flow through a fracture. Assuming that the joint elements connecting ΔF1F2F3 and ΔF10F11F12 and connecting ΔF4F5F6 and ΔF7F8F9 are broken to form a fracture, and taking heat node T1 as an example, two triangular meshes are 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 meshes can be calculated, and the calculation formula is as follows:

[0092]

[0093] In the formula, is the heat flowing into T1 through the triangular mesh; Q ΔF1F2F3→F1 is the heat flowing into F1 through ΔF1F2F3; Q ΔF4F5F6→F4 is the heat flowing into F4 through ΔF4F5F6.

[0094] Taking Q ΔF1F2F3→F1 as an 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 are the x-direction components and y-direction components of the coordinate interpolation of the edges F1F2 and F1F3. q x and q y are calculated as follows:

[0097]

[0098] In the formula, k is the thermal conductivity; and are the x and y coordinates of the corresponding nodes; are the temperatures of the corresponding nodes, a = 1, 2, 3.

[0099] At the same time, the triangular meshes on both sides of the fracture also exchange heat, and the heat flowing into the T1 node through the fracture is calculated according to the following formula:

[0100]

[0101] In the above formula, is the heat flowing into T1 through the fracture; QΔF10F11F12→F1 Q = QF1F2F3 + QF1F4F5 + QF1F6F7 + QF1F8F9 + QF1F10F11 + QF1F12F13 (1) ΔF7F8F9→F4 Q = QF1F2F3 + QF1F4F5 + QF1F6F7 + QF1F8F9 + QF1F10F11 + QF1F12F13 (1)

[0102] Q = QF1F2F3 + QF1F4F5 + QF1F6F7 + QF1F8F9 + QF1F10F11 + QF1F12F13 (1) ΔF10F11F12→F1 For example, the calculation formula is as follows:

[0103] Q = QF1F2F3 + QF1F4F5 + QF1F6F7 + QF1F8F9 + QF1F10F11 + QF1F12F13 (1) ΔF10F11F12→F1 = 0.5 * k J (T F10 -T F1 )L (7)

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

[0105] Step S4.2. According to the temperature of each thermal calculation node calculated in step S4.1 , the temperature of each thermal calculation node at this time is calculated by substituting the temperature of each thermal calculation node calculated in step S4.1

[0106]

[0107] In the above formula, and are the temperatures of the i-th thermal calculation node Ti in the model at times 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; and M i is the mass of 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 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 triangular calculation grid ΔF j F j+1 F j+2 , the load f T caused by the temperature change is calculated according to formula (9):

[0109]

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

[0111] Step S4.4. According to the latest temperature field distribution of the rock mass model obtained in step S4.3, the load caused by temperature change is introduced 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 mass 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, the kth force calculation node is associated with the ith thermal calculation node; f ext is the external force received by the node, including body force, external load, contact force and joint force, etc.; f int is the internal force received by the node, and f T is the temperature load.

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

[0115]

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

[0117] Step S5. After each calculation step is completed, it is determined whether the joint element is damaged according to the deformation and stress state of the joint element, and the mapping relationship between the thermal calculation node system and the force calculation node system is updated. The modes of joint element damage are 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 failure: d s0 <d s <d sc and dn0 <d n <d n0 , satisfy

[0121] where d n and d s are the normal and shear deformation of the joint element, respectively; the rest of the 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 the tensile and shear strength, which can be measured by experiments; k n and k s are the normal and shear joint stiffness, which are the meso-scale calculation parameters, usually take the multiple of the elastic modulus; G I and G II are the I and II type fracture energy.

[0127] After the joint element is destroyed, the mapping relationship between the thermal node system and the force node system needs to be updated, as shown in the schematic Figure 4 (b)-(c). When the joint element connecting ΔF4F5F6 and ΔF7F8F9 is destroyed, the linked list (T1:{F1→F4→F7→F10→F1}) storing the mapping relationship needs to be updated, and the updated linked list structure still needs to satisfy the counterclockwise order, that is, T1:{F7→F10→F1→F4}. With the further expansion of the crack, the joint element connecting ΔF1F2F3 and ΔF10F11F12 is destroyed, at this time there is no joint element connection between the upper and lower triangular grid, the model is divided into two parts, a new thermal calculation node T2 needs to be added, T1 and T2 are associated with the upper and lower force calculation nodes, respectively, the original T1 linked list is updated to T1:{F1→F4}, and the newly generated T2 linked list also satisfies the counterclockwise order, that is, T2:{F7→F10}.

[0128] Step S6. Repeat the above steps S5 to the end of the last calculation step, then the thermal cracking process simulation of the rock mass is completed. The following will be described in combination with a specific application example:

[0129] Example 1:

[0130] In this embodiment, the tunnel engineering is taken as an example to divide the potential damage zone and the non-damage zone and to apply the strategy of joint element. As shown in Figure 2 to avoid the influence of boundary effect on the simulation accuracy, the rock mass modeling area usually needs to be much larger than the size of the tunnel, however, the actual engineering mainly concerns the surrounding rock near the tunnel boundary, and the surrounding rock area affected by external disturbance is limited. According to the research of the 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, the size thereof can be set to a value slightly larger than the range of the disturbance zone. Similarly, for a tunnel of any shape, a multiple of the maximum size thereof is taken as the size of the potential damage zone. After the model is discretized by using a triangular calculation grid, joint elements are inserted inside the potential damage zone and on the boundary thereof, and then a complete calculation model is obtained.

[0131] Example 2:

[0132] In this embodiment, the thermal conduction simulation of a rock column containing joint elements is taken as an example to divide the potential damage zone and the non-damage zone and to apply the strategy of joint element. The generation strategy of the thermal calculation node system and the force calculation node system in the present application is illustrated, and the advantages of the present method are embodied. As shown in Figure 6 (a), it is a rectangular rock column with a length of 12 m and a width of 1 m. The model is discretized by using a triangular grid, and 386 triangular grids are obtained, as shown in Figure 6 (b). Before the joint elements are inserted, adjacent triangular grids have a common edge, and the nodes on the edge are shared by the two triangular grids, and there are 246 nodes. The 246 nodes are set as thermal calculation nodes, and then the joint elements are inserted inside the model, at this time, the nodes of all triangular grids are independent of each other, and the total number of nodes is three times the number of grids, that is, 386*3=1158 nodes, which are the force calculation nodes. The force calculation nodes with the same coordinates belong to the same thermal calculation node, and form the mapping relationship described in step S3 of the present application, as shown in Figure 6 (c).

[0133] As a comparison, two methods are used to simulate the heat transfer in the model: the first method is the method proposed by the author, which is referred to as Method 1; the second method introduces a joint heat transfer model, which is the most common heat transfer model in current non-continuous methods and hybrid methods. This model assumes that heat exchange occurs between adjacent elements through joint elements. By introducing an artificially assumed joint heat transfer coefficient, heat transfer simulation between models under unbroken conditions can be achieved. This method is referred to as Method 2 hereinafter. According to the research in 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 real 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 length of the joint, and the coefficient n should be greater than 100, which is taken as 100 here.

[0136] The initial temperature of the rock pillar is 0°C, a constant temperature heat source is applied at the bottom, keeping the temperature at T B = 100°C, while the temperature at the top is fixed at T T = 0°C. A horizontal measuring line is placed in the middle of the rock pillar to monitor the temperature changes. 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 of the point from the bottom boundary; L is the length of the model; the thermal conductivity k = 10 W / (m·°C); the density p = 2000 kg / m 3 ; and the specific heat capacity c p = 1 J / (kg·°C).

[0140] The temperature evolution curves of the measuring line simulated by Method 1 and Method 2 under different calculation steps At are shown in Figure 7 and Figure 8 respectively. Among them, Figure 7(a), 7(b), 7(c), and 7(d) represent the temperature evolution curves of the measuring line obtained by simulating the method 1 when Δt is 4s, 2.5s, 2s, and 1s, respectively. Figure 8 (a) and (b) respectively represent 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, as 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 calculation step size, when Δt=1×10 -2 s, the simulation results are consistent with the analytical solution, such as Figure 8 (b) It's worth noting that, given the same simulation time, the step size is inversely proportional to the total number of computational steps; a smaller step size means more computational steps are required. Taking this example, the number of computational steps between Method 1 and Method 2 can differ by a factor of 250. Method 1 takes only 4 seconds, while Method 2 takes 406 seconds.

[0142] Example 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. Figure 9 The figure shows a square rock sample with nine natural cracks. The initial temperature of the rock is 0°C, and a constant temperature heat source is applied to the top of the rock to maintain its temperature at 100°C. Figure 10 To simulate the temperature evolution process, Figure 10 (a)-(d) It can be seen that a top-down heat flow 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 a significant temperature difference on both sides of the crack, which slows down the speed of heat flow.

[0144] Example 4:

[0145] In this embodiment, the numerical simulation method of the present invention is described by taking the thermal fracture simulation of geothermal mining in conglomerate strata as an example, dividing the potential damage zone and the indestructible zone and applying the strategy of the joint unit with a specific example. Figure 11As shown, it is a deep buried geothermal reservoir, and the lithology is conglomerate. In order to extract geothermal energy, the energy pipeline is laid in the rock layer, and the heat exchange with the surrounding rock is realized by pouring cold water to bring the geothermal energy back to the ground. Assuming that the initial temperature of the surrounding rock is 200 DEG C, and the initial water temperature in the pipeline is 20 DEG C. The conglomerate is composed of two parts of block stone and matrix. Considering that the strength of the block stone is usually much larger than that of the matrix structure, it is difficult to be damaged, so the block stone is set as an inviolable area in the modeling process, and the matrix is set as a potential damage area to simulate the thermal cracking phenomenon that may occur in the heat exchange process.

[0146] Figure 12 In order to simulate the crack propagation and temperature evolution process, the temperature field of the surrounding rock is calculated from the heat conduction equation. Figure 12 (a)-(d), with the extraction of geothermal energy, the temperature of the surrounding rock near the pipeline gradually decreases, and additional temperature stress is generated. When the temperature stress reaches the strength limit, the surrounding rock near the pipeline is damaged first, and the crack extends away from the pipeline. Due to the initiation and propagation of the crack, the temperature field in the rock mass appears discontinuous, and the temperature on both sides of the crack jumps.

[0147] The above only describes the embodiments of the present application, and does not limit the patent range of the present application, and any equivalent structure or equivalent process transformation using the content of the specification and drawings, or direct or indirect application in other related technical fields, are also included in the patent protection range of the present application.

Claims

1. A continuous-discontinuous numerical simulation method considering the thermal-mechanical coupling problem of rock mass, characterized by: The following steps are involved: Step S1. Divide the simulation area of ​​the rock mass model into a potential damage area and an indestructible area according to the stress characteristics of the rock mass, and perform grid division on the potential damage area and the indestructible area; Step S2. Setting a joint unit at each grid boundary of the potential failure zone. The joint unit connects the grids on both sides of the joint unit, and the rock mass fracture process is simulated through the destruction of the joint unit. Step S3. Based on the model meshing and joint element insertion results, a thermal calculation node system and a force calculation node system are generated for simulation calculations. A mapping relationship exists between the two, and the mapping relationship is continuously and dynamically updated during the simulation of 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 sequence; 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, thus completing the simulation of the thermal cracking process of the rock mass.

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 step S1, the stress concentration area of ​​the rock mass caused by internal defects or external disturbances 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 propagation 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 the model and on the boundary, and three adjacent nodes are connected 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 completely consistent, and there is a situation where multiple adjacent triangular calculation meshes 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 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 flow 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, then Equivalent to the heat flowing into Ti through multiple triangles that share the node sum; 2) When the joint unit is damaged and cracks exist, the cracks divide the multiple triangular computational meshes that share the node into two parts, one on each side of the crack. Calculated using formula (1): In the above formula, The heat flowing into Ti for the triangle calculation mesh located on the same side of the thermal calculation node Ti is calculated. The heat flowing into node Ti through heat exchange through the crack is calculated for the triangle mesh on the other side of the crack; Step S4.

2. Calculated according to step S4.1 Substitute the rock mass heat conduction formula (2) to 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. Based on the current temperature of each thermal calculation node 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. Based on 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. 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 triangle computing grids, where each triangle 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 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 exchanged. 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; where 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 is calculated using formula (4.3): In the above formula, k J is the crack heat conductivity coefficient; and Node F p and F j The temperature of the crack; L is the length of the crack; the node F j and F p When no joint unit is inserted, they overlap 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, all of which can 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) Where, 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 mode I and II fracture energies, respectively.

Images