A simulation method and system for phase change-hydrothermal force coupled crack propagation of fractured rock mass for low-temperature underground engineering

Abstract

This invention discloses a method and system for simulating crack propagation in fractured rock masses using a phase transformation-hydrothermal coupling approach for cryogenic underground engineering, belonging to the fields of rock mechanics and underground engineering safety evaluation. The method includes: firstly, establishing an engineering stage model of the fractured rock mass, simultaneously treating fractured elements as mechanical, hydraulic, thermal, and phase transformation state units; then identifying the aqueous phase, transitional phase, and icy phase, and calculating the equivalent phase transformation normal load; establishing a temperature-direction dual-variable parameter library, and jointly correcting strength and fracture toughness; inheriting displacement discontinuities, fracture aperture, temperature, fluid pressure, and phase state identifiers during engineering stage transitions, thereby simulating crack initiation, propagation, and seepage channel reconstruction throughout the entire process of excavation, cooling, water injection / drainage, and operation. This invention can output crack propagation paths, phase transformation damage zones, ice-induced crack segments, and seepage channel reconstruction zones, and can be used for stability evaluation of cryogenic underground storage facilities, cold-region tunnels, and deep cryogenic disturbance engineering projects.

Description

Technical Field

[0001] This invention relates to the fields of rock mechanics, underground engineering safety assessment and numerical simulation technology, specifically to a method and system for simulating crack propagation in fractured rock masses using phase transformation-hydrothermal coupling for low-temperature underground engineering. Background Technology

[0002] In underground engineering, fractured rock masses are often simultaneously affected by in-situ stress, excavation unloading, fracture water pressure, temperature changes, and material heterogeneity. Numerical methods are available to calculate fracture propagation, hydromechanical coupling, hydrothermal coupling, or thermomechanical coupling. For example, the displacement discontinuity method and the boundary element method can be used to simulate the mechanical response and crack propagation of explicit fractures; fracture flow can be calculated based on fracture aperture and pressure difference; and thermal boundary conditions and heat sources can be used to solve for temperature fields and thermal stresses.

[0003] However, fissure water in low-temperature underground engineering undergoes phase change. The ice expansion caused by phase change is not a simple constant pressure, but changes with fissure aperture, fissure normal deformation, ice layer stiffness, and surrounding rock contact stiffness. Simultaneously, temperature decreases or increases alter rock tensile strength, cohesion, internal friction angle, and fracture toughness; rock masses with well-developed bedding, foliation, or joints also exhibit significant direction-dependent strength. If conventional water-mechanical or water-thermal coupling models are still used, the effects of phase change pressure feedback, temperature-directional strength correction, and boundary reconstruction during construction on crack propagation paths are easily overlooked. Currently, when overcoming the shortcomings of water-mechanical or water-thermal coupling models, a common inherent technical bias among those skilled in the art is that ice expansion in fissure rock masses can be simplified to a constant pressure or equivalent thermal load acting on the fissure walls. This bias leads technicians to ignore the physical nature of the dynamic feedback evolution of ice expansion pressure with fracture aperture, ice layer stiffness, and surrounding rock contact state, and also hinders attempts to incorporate fracture water phase change, hydrothermal coupling, and crack propagation processes into a unified computational framework. Influenced by this bias, existing simulation methods typically treat hydro-mechanical coupling, hydro-thermal coupling, or thermo-mechanical coupling separately, or calculate engineering stages in isolation, making it difficult to accurately reflect the crack initiation, propagation, and seepage channel remodeling behavior of fractured rock masses throughout the entire process of excavation, cooling, water injection / drainage, and operation of low-temperature underground engineering projects.

[0004] Furthermore, underground engineering projects typically evolve in stages such as excavation, support, cooling, drainage, and operation. The boundary conditions and fracture networks of subsequent stages depend on the residual deformation, fracture aperture, phase state, and fluid pressure of the previous stage. Existing calculation processes often treat engineering stages in isolation, making it difficult to directly provide results that can be used for engineering decision-making, such as phase transformation damage zones, ice expansion-controlled cracks, and seepage channel reconstruction zones.

[0005] This shows that the existing technology needs further improvement. Summary of the Invention

[0006] One of the objectives of this invention is to provide a simulation method for phase transformation-hydrothermal coupling crack propagation in fractured rock mass for low-temperature underground engineering, which can more accurately simulate the initiation, propagation, and penetration of cracks in fractured rock mass throughout the entire process of excavation, cooling, water injection and drainage, and operation.

[0007] To achieve the above objectives, the present invention adopts the following technical solution:

[0008] A method for simulating crack propagation in fractured rock masses using a phase transformation-hydrothermal coupling approach for cryogenic underground engineering includes the following steps:

[0009] S1. Establish a fractured rock mass engineering stage model that includes surrounding rock boundary, excavation boundary, existing fractures, candidate crack ends, thermal boundary, hydraulic boundary, mechanical boundary, and at least one engineering stage sequence.

[0010] S2. Discretize the fractures, boundaries and crack ends in the fractured rock mass engineering stage model into computational units, and establish state variables for each computational unit. Use each computational unit as a common carrier for mechanical units, hydraulic units, thermal and phase change state units.

[0011] S3. At each engineering stage and time step, phase transition identification is performed on the computational unit according to the temperature field and water content, and the unit is divided into aqueous phase unit, transition phase unit or ice phase unit.

[0012] S4. Calculate the equivalent phase transformation normal load that changes with the fracture normal deformation feedback for the water phase element and the ice phase element, and input the equivalent phase transformation normal load together with water pressure, in-situ stress and external boundary load into the fracture mechanics solution.

[0013] S5. Based on the temperature and crack orientation of the calculation unit, tensile strength, cohesion, internal friction angle, type I fracture toughness and type II fracture toughness are obtained or interpolated from the temperature-direction bivariate parameter library.

[0014] S6. Based on the equivalent phase transformation normal load, fluid pressure and temperature-corrected strength parameters, calculate crack deformation, crack initiation index and crack propagation index, and update crack aperture and fluid channel.

[0015] S7. When the stage transition conditions are met, reconstruct the excavation boundary or thermal / hydraulic boundary, and map the displacement discontinuity, fluid pressure, temperature, phase identification and crack aperture of the previous engineering stage to the next engineering stage.

[0016] S8. Repeat S3 to S7 until the set engineering time or convergence condition is reached, and output the crack propagation path, phase transformation damage zone and seepage channel reconstruction zone.

[0017] In the above-mentioned simulation method for phase transformation-hydrothermal coupling crack propagation in fractured rock mass for low-temperature underground engineering, step S2 uses the displacement discontinuity method to discretize fractures of arbitrary shapes; in the fracture element, the tangential stress and normal stress are obtained by superimposing the boundary influence coefficient and the displacement discontinuity.

[0018] In the above-mentioned simulation method for phase transformation-hydrothermal coupling crack propagation in fractured rock mass for low-temperature underground engineering, in step S3, the water phase unit has a temperature > 0℃; the transition phase unit has a temperature of -2.5℃ ≤ temperature ≤ 0℃; and the ice phase unit has a temperature < -2.5℃.

[0019] The above-mentioned simulation method for phase transformation-hydrothermal coupling crack propagation in fractured rock mass for low-temperature underground engineering, in step S4, the equivalent phase transformation normal load is determined by the initial phase transformation pressure, the equivalent stiffness of the ice layer, the contact stiffness of the fracture, and the increment of the fracture normal displacement, so that the phase transformation pressure is corrected by feedback as the fracture opens or closes; for the j-th fracture element, the calculation formula for the ice expansion pressure is shown in equation (1):

[0020] (1);

[0021] In formula (1): The initial ice expansion pressure is for unjointed deformation; For the pressure of ice expansion, The equivalent stiffness of the ice layer; This represents the increment of the normal displacement of the crack.

[0022] The above-mentioned simulation method for phase transformation-hydrothermal coupling crack propagation in fractured rock mass for low-temperature underground engineering uses the following formulas: the calculation formula for ice expansion pressure is shown in Equation (2), and the calculation formula for equivalent stiffness of ice layer is shown in Equation (3).

[0023] (2);

[0024] (3);

[0025] In equation (2), This represents the elastic modulus during the initial stage of freezing.

[0026] In equation (3), e is the hydraulic opening;

[0027] The formula for calculating the contact stress in fracture mechanics is shown in equation (4):

[0028] (4);

[0029] In equation (4): The normal contact stiffness of the crack; This represents the initial displacement in the crack normal direction.

[0030] The formula for calculating the total fracture ice expansion pressure is shown in equation (5):

[0031] = (5).

[0032] The above-mentioned simulation method for phase transformation-hydrothermal coupling crack propagation in fractured rock mass for low-temperature underground engineering determines the equivalent stiffness of the ice layer based on the elastic modulus of the ice, the hydraulic aperture of the fracture, the phase transformation participation coefficient, and the minimum aperture limit.

[0033] The above-mentioned simulation method for phase transformation-hydrothermal coupling crack propagation in fractured rock mass for low-temperature underground engineering includes the following steps: In step S5, the temperature-direction bivariate parameter library includes tensile strength, cohesion, internal friction angle, Type I fracture toughness, and Type II fracture toughness recorded by temperature nodes and azimuth nodes, and the calculation parameters are generated by interpolation for non-node temperature and azimuth; In step S6, the crack propagation index includes Type I propagation components and Type II propagation components, which are normalized by temperature-direction corrected Type I fracture toughness and Type II fracture toughness, respectively.

[0034] The aforementioned simulation method for phase change-hydrothermal coupling crack propagation in fractured rock masses for low-temperature underground engineering includes, in step S6, a fracture hydraulic update step, a fracture hydrothermal exchange step, and a damping update step for fracture aperture changes. The fracture hydraulic update step calculates the inter-domain flow rate based on the average aperture, pressure difference, fluid viscosity, and element length of adjacent fracture domains, updates the fluid pressure and hydraulic connectivity based on the inter-domain flow rate, and calculates the equivalent leakage from the fracture to the surrounding rock based on the surrounding rock permeability coefficient, effective leakage distance, fracture domain pressure, and initial pore water pressure. The fracture hydrothermal exchange step calculates the heat carried by the fluid based on the flow rate entering the fracture domain, fluid density, fluid specific heat, upstream fracture domain temperature, and current fracture domain temperature, and uses this heat as an equivalent heat or cold source for the fracture wall surface in the temperature field update. The damping update step for fracture aperture changes includes updating the fracture aperture, fluid pressure, temperature, and phase change pressure step by step according to the damping coefficient when the fracture aperture change caused by a single iteration exceeds a preset proportion.

[0035] The above-mentioned simulation method for phase transformation-hydrothermal coupling crack propagation in fractured rock mass for low-temperature underground engineering includes the following steps: in step S7, the stage transition conditions include the completion of the excavation step, change of water injection pressure, change of drainage boundary, change of temperature boundary, change of insulation layer parameters, or reaching the set crack propagation length; in step S8, the output phase transformation damage zone includes the continuous distribution area of ​​ice phase units, the concentrated area of ​​transition phase units, the fracture segment where the phase transformation normal load exceeds the threshold, and the candidate crack tip where the crack propagation index exceeds the threshold.

[0036] Another objective of this invention is to provide a simulation system for phase transformation-hydrothermal coupling crack propagation in fractured rock masses for low-temperature underground engineering. The system executes the above-described method and includes a fractured rock mass engineering stage model module, a temperature field update module, a phase transformation identification module, an equivalent phase transformation load module, a hydraulic update module, a hydrothermal exchange module, a temperature-direction parameter correction module, a crack propagation determination module, an engineering stage update module, and a result output module.

[0037] Compared with the prior art, the present invention brings the following beneficial technical effects:

[0038] ① The feedback processing of phase change load improves the rationality of frost heave force calculation.

[0039] Existing methods typically simplify frost heave in fractured rock masses as a constant pressure or equivalent thermal load, neglecting the physical fact that frost heave force dynamically changes with fracture opening and closing and the stress state of the ice layer. This invention establishes joint governing equations (Equations 1-5) between the equivalent stiffness of the ice layer, the contact stiffness of the fracture, and the increment of the fracture normal displacement. This allows the phase change pressure to self-correct as the fracture opens or closes, avoiding calculation errors caused by the constant pressure assumption. This approach provides a more physically accurate numerical implementation path for assessing frost heave force in low-temperature fractured rock masses.

[0040] ② The dual-variable correction of the strength parameters by temperature and direction enriches the applicable conditions of the fracture criterion.

[0041] In practical engineering, the tensile strength, cohesion, internal friction angle, and fracture toughness of fractured rock masses are affected by both temperature and exhibit significant direction dependence. The temperature-direction bivariate parameter library (M(T,β)) constructed in this invention can interpolate the corresponding strength and fracture toughness parameters based on the fracture tip temperature and candidate propagation directions, making the determination of crack initiation and propagation closer to the actual mechanical response of the material. For isotropic rock masses or cases where temperature influence is insignificant, this parameter library can also degenerate into a univariate or constant form, providing a degree of flexibility.

[0042] ③ The inheritance mechanism of state variables between engineering stages improves the initial mismatch problem caused by stage separation.

[0043] Underground engineering projects are typically constructed in stages. The crack propagation behavior in subsequent stages largely depends on the residual deformation, crack aperture, temperature distribution, fluid pressure, and phase distribution formed in the earlier stages. This invention inherits and maps these state variables during stage transitions, rather than re-initializing them. This helps maintain the continuity of the simulation process, reduces the initial condition deviations introduced by stage segmentation, and provides a more consistent computational framework for the analysis of the entire multi-stage construction process.

[0044] ④ The simultaneous output of evaluation indicators for multiple types of projects enhances the engineering reference value of the simulation results.

[0045] In addition to outputting crack propagation paths, this invention can also simultaneously provide information on the continuous distribution area of ​​ice phase units, the concentrated area of ​​transition phase units, the crack segment where the phase transformation normal load exceeds the threshold, and the seepage channel reconstruction area. These indicators correspond to the frost heave influence range, the phase transformation transition region, the main controlling ice heave crack segment, and the potential seepage path, respectively. This provides multiple reference bases for the stability assessment of surrounding rock in cryogenic reservoirs, the thermal insulation design of tunnels in cold regions, and the optimization of drainage schemes, and has certain practical engineering value.

[0046] In summary, this invention has made improvements over existing technologies in the parameterized description of the phase change-water thermodynamic coupling mechanism, the strength correction method, and the stage connection strategy, providing a technical solution for simulating crack propagation in fractured rock masses in low-temperature underground engineering. Attached Figure Description

[0047] The present invention will be further described below with reference to the accompanying drawings:

[0048] Figure 1 This is a flowchart of the simulation method of the present invention.

[0049] Figure 2 A schematic diagram of a crack is shown, representing the discontinuous displacement of N elements.

[0050] Figure 3 To identify the low-temperature aqueous phase-transition phase-ice phase and the phase transition pressure gradient.

[0051] Figure 4 This is a schematic diagram for calculating the equivalent phase transition normal load.

[0052] Figure 5 Delineation of fissure water basins.

[0053] Figure 6 Region partitioning for hydraulically coupled fluid flow simulation.

[0054] Figure 7 The fluid flow within the hydrothermal coupling region j. Detailed Implementation

[0055] This invention proposes a method and system for simulating crack propagation in fractured rock masses using a phase transformation-hydrothermal coupling approach for low-temperature underground engineering. To make the advantages and technical solutions of this invention clearer and more explicit, the invention will be further described below with reference to specific embodiments.

[0056] This invention proposes a simulation method for phase transformation-hydrothermal coupled crack propagation in fractured rock masses for low-temperature underground engineering, such as... Figure 1 As shown, the specific steps include:

[0057] Step 1: Establish the engineering stage model of fractured rock mass: Input the engineering boundary, existing fractures, excavation steps, thermal boundary, hydraulic boundary, mechanical boundary, rock mass elastic parameters, fracture contact parameters, fluid parameters, phase transformation parameters, temperature-related strength parameters, and direction-related strength parameters into the computer program. The processor executes at least one engineering stage sequence, including the engineering boundary and existing fractures, to simulate crack initiation and propagation in elastic and isotropic rock media.

[0058] Step 2, Discretizing fractures and boundaries: The fractures, boundaries, and crack ends in the fractured rock mass engineering stage model are discretized into computational units, and state variables are established for each computational unit. Each computational unit serves as a common carrier for mechanical units, hydraulic units, thermal and phase change state units.

[0059] (1) Discontinuous displacement: Discretize arbitrary crack shapes into N linear crack elements, such as Figure 2 As shown, each element has a local tangential direction s and a normal direction n, defining the discontinuous displacement of the element. Elements are chosen to occupy a specific portion of the s-axis, for example, |s|≤a, n=0. If this element is considered a line crack, it can be addressed by setting a surface on the positive side of n=0 (denoted as n=0). + The other is on the negative side (represented as n=0). - The two surfaces are distinguished by a coordinate system. When moving from one side of an element to the other, the displacement undergoes a continuous, specific change in the value Di = (Ds, Dn). In the diagram, 1, 2, 3, and N are the element numbers, and 2a is the element length. The element positions are determined by the x, y coordinate system. If the crack surfaces are subjected to stress, they will move relative to each other. The displacement discontinuity method is a discrete approximation for finding a smooth distribution of actual relative displacements; the discrete approximation is based on... Figure 2 The crack shown is found through N subdivisions, each subdivision being a boundary element representing a discontinuity displacement of the element. The displacement discontinuity of the j-th crack element is defined as the relative displacement difference between the upper and lower crack surfaces:

[0060] ;

[0061] ;

[0062] and This refers to the shear force and normal displacement of the j-th crack element. The superscripts "+" and "-" indicate the positive and negative surfaces of the crack relative to local coordinate n.

[0063] Local displacement and The two components of the vector are formed. Regardless of whether the crack's positive or negative surface is considered, they are always positive in the positive directions of s and n. Therefore, from the equation, it can be derived that if the two surfaces of the crack are displaced relative to each other, the normal component of the discontinuous displacement is... just.

[0064] (2) Tangential stress and normal stress: In a fractured element, tangential stress and normal stress can be obtained by superimposing the boundary influence coefficient and the displacement discontinuity:

[0065] ;

[0066] ;

[0067] In the formula The numbers are 1, 2, and up to N. and They represent the first The tangential and normal stresses of each element; , , and The boundary influence coefficient of stress is determined by the element geometry and elastic parameters; and Indicates the first Discontinuous displacements of individual elements; and It is the component of far-field stress or external boundary load in local coordinates.

[0068] Rock discontinuities can exist in three states: open, elastic contact, or sliding.

[0069] Cracks opened: ;

[0070] Elastic contact: and The size depends on the crack stiffness ( ) and discontinuous displacement ( );

[0071]

[0072]

[0073] slide:

[0074]

[0075]

[0076] in It is the friction angle of the crack surface. The sign depends on the direction of the slide.

[0077] Based on this, the present invention integrates water pressure, thermal stress and equivalent phase transformation normal load into the crack mechanics solution to obtain updated results of crack aperture and displacement discontinuity.

[0078] Step 3: Temperature field update: Solve or read the temperature field for the current time step, and map the temperature parameters to the crack calculation unit, crack tip, and potential crack initiation point.

[0079] Step 4: Phase Transition State Identification: Identify the water phase, transition phase, and ice phase based on temperature, saturation, and fracture aperture, and calculate the phase transition participation coefficients for the water and ice phases. Liquefied natural gas underground storage caverns operate at extremely low temperatures (-160°C). During normal storage periods, the surrounding rock will cool to sub-zero temperatures depending on time and distance from the cavern walls. If the rock mass is already saturated before cooling, ice will form in the initially water-filled pores and fractures. Ice formation generates significant expansion pressure. Under low-temperature conditions, three distinct stages of the water phase may exist within the rock mass, such as… Figure 3 As shown.

[0080] Aqueous phase unit: temperature > 0℃; Transitional phase unit: -2.5℃ ≤ temperature ≤ 0℃; Ice phase unit: temperature < -2.5℃.

[0081] Step 5: Calculation of Equivalent Phase Transformation Normal Load: For ice phase or transition phase elements, calculate the initial phase transformation pressure, equivalent ice layer stiffness, and fracture contact stiffness, and correct the phase transformation pressure based on the fracture normal displacement increment. In the transition phase, a temperature-dependent phase transformation pressure gradient can be introduced.

[0082] ;

[0083] In the formula, This is the phase transition pressure gradient that varies with temperature. The total expansion pressure during this stage can be calculated using the following formula:

[0084] ;

[0085] In the formula, The temperature range is -2.5 to 0℃.

[0086] In the ice phase, the transformation of water into ice can be described by the coefficient of volume expansion:

[0087] ;

[0088] To simplify the modeling process of ice rings, two stages were considered: a liquid phase and a solid phase, with the phase transition occurring at a temperature of 0°C. In rock joints, joints may be acted upon by two forces, such as... Figure 4 .

[0089] ①Contact stress on the joint surface ( The joint stiffness and joint deformation are determined by the joint contact stiffness and joint deformation.

[0090] ② The expansion pressure from ice ( ), is a function of ice stiffness and joint deformation.

[0091] The formula for calculating the ice expansion pressure for the j-th fracture element is as follows:

[0092] ;

[0093] In the formula The initial ice expansion pressure is the pressure at which joints are not deformed.

[0094] Let its hydraulic opening be e, and its elastic modulus at the initial stage of freezing be... The equivalent stiffness of the ice layer can then be written as:

[0095] ;

[0096] ;

[0097] In the formula This represents the equivalent stiffness of the ice layer.

[0098] The formula for calculating contact stress in fracture mechanics is as follows:

[0099] ;

[0100] In the formula For the normal contact stiffness of the crack, The initial displacement is the normal direction of the crack. This represents the increment of the normal displacement of the crack.

[0101] Therefore, the formula for calculating the total fracture ice expansion pressure is as follows:

[0102] ;

[0103] In the formula It is the initial stress. The total stiffness of the ice-filled fissures is considered. Therefore, ice expansion is no longer simply treated as a constant pressure, but rather evolves in conjunction with fissure opening and closing, the equivalent stiffness of the ice layer, and the contact stiffness of the fissures. When a fissure opens, the ice expansion pressure is released; when a fissure closes, the ice layer is compressed and intensifies within an upper limit.

[0104] Step 6: Water-thermal-mechanical iteration. Using water pressure, equivalent phase change normal load, thermal stress, and in-situ stress as boundary conditions, solve for crack deformation and update crack aperture, fluid pressure, temperature, and flow channels.

[0105] (1) The fracture hydraulic renewal module mainly focuses on fluid flow in rock fractures; however, leakage from the fracture fluid channel to the rock mass is also considered, such as... Figure 5 As shown. In this case, the fluid flows not only within the fracture channels, but also between the fracture channels and the surrounding rock.

[0106] When using the displacement discontinuity method for mechanical numerical calculations, the crack is divided into several displacement discontinuity elements. During the calculation of fluid flow, each displacement discontinuity element is considered a hydraulic region, with adjacent regions interconnected, such as... Figure 6 As shown, under the influence of pressure difference, fluid flows from one region to another.

[0107] For adjacent fracture domains i and j, the inter-domain flow can be calculated in the form of fracture channel flow:

[0108] ;

[0109] In the formula: -Average aperture of adjacent fracture domains;

[0110] -Equivalent length of adjacent domains;

[0111] – Inter-domain pressure difference;

[0112] – Fluid viscosity.

[0113] Calculated and It is used to update the pressure in the fracture domain and influence fracture deformation, aperture, and the generation of new fracture channels through pressure change feedback.

[0114] (2) Fracture water-heat exchange module: When fluid flows in the fracture and the fluid temperature is inconsistent with the surrounding rock temperature, the heat carried by the fluid will change the temperature of the fracture wall. For example Figure 7 As shown, fluid flows from fracture region i to region j, and the rock temperatures at units i and j are T, respectively. i and T j The fluid temperature entering region j from region i is T. i As the fluid moves from one end of region j to the other, it is "heated" by the fracture walls, and when it reaches the other end, the fluid temperature becomes T. j .

[0115] In this process, the fluid absorbs heat energy from the fracture wall through fracture region j, which acts as a negative heat source, causing the rock temperature to decrease in the next step. For fracture region j, the volume of fluid flowing in from region i within time step Δt is:

[0116] ;

[0117] The heat exchanged between the fluid and the fracture wall can be expressed as:

[0118] ;

[0119] In the formula: -Specific heat of fluid; – Fluid density; –Upstream fluid temperature; -Current fracture wall temperature; -Time step.

[0120] This is transmitted from the rock mass to the region. The thermal energy of the fluid. The temporary heat source is calculated over the time step using the following formula. The heat intensity generated internally:

[0121] ;

[0122] The heat or cold source is input into the temperature field update module to determine the fracture temperature and phase transition state for the next time step. Note the calculation time step. When doing so, the cumulative thermal intensity of all previous steps should be used.

[0123] Step 7: Temperature-direction bivariate strength correction: Based on the crack tip temperature and the candidate crack propagation direction, interpolation is performed to obtain the corresponding strength and fracture toughness parameters.

[0124] Experimental data show that many rock properties are temperature-dependent. Key rock property parameters known to exhibit significant temperature dependence include: ① uniaxial compressive strength decreases with increasing temperature; ② tensile strength decreases with increasing temperature; ③ fracture toughness decreases with increasing temperature; ④ thermal conductivity decreases with increasing temperature; and ⑤ the coefficient of thermal expansion increases with increasing temperature.

[0125] To accurately simulate the thermodynamic behavior of rocks, this invention can easily account for the temperature dependence of rock strength parameters, such as compressive strength, tensile strength, and fracture toughness. It includes the following steps:

[0126] ① Establish a parameter library M(T,β), where T is the temperature at the calculation point and β is the direction of the candidate crack surface or crack propagation direction, i.e., the azimuth angle. The parameter library includes strength parameters with several temperature values, such as tensile strength. Cohesion internal friction angle Type I fracture toughness Type II fracture toughness For any temperature between two given values, the intensity parameter can be calculated using a linear function.

[0127] Fracture occurs when the combination of the two principal stresses reaches a critical value. The critical value is the rock's initial fracture strength, which is usually proportional to the rock's overall strength. Using the Mohr-Coulomb criterion, the rock's compressive strength is given by the following formula:

[0128] ;

[0129] In the formula and These are the maximum and minimum principal stresses; φ is the angle of internal friction in the rock; c is the cohesion. The direction of the failure surface is perpendicular to... The included angle

[0130] ;

[0131] The tensile strength is given by the following formula:

[0132] ;

[0133] In the formula, For the tensile strength of intact rock, the failure surface is... In the direction of. The above two equations only apply to isotropic rocks.

[0134] ② Determine the temperature values ​​at points inside the rock and at the tips of the cracks.

[0135] ③ Use the obtained temperature calculation point or the intensity value of the tip.

[0136] ④ Determine whether fracture initiation and / or fracture propagation will occur.

[0137] Step 8: Crack Initiation and Propagation Determination: Calculate the propagation index of candidate crack directions and select the direction with the largest propagation index as the priority propagation direction; when the propagation index reaches the threshold, a new crack element is generated.

[0138] Use the F-criteria to determine the likelihood and direction of crack propagation:

[0139] ;

[0140] In the formula and For strain energy release rates of type I and type II, and This represents the critical strain energy release rate for Type I and Type II cracks. The crack propagation direction is the direction in which F reaches its maximum value.

[0141] For Type I and Type II, the F criterion can also be expressed using the fracture stress intensity factor ( , ) and fracture toughness ( and The following is an example:

[0142] ;

[0143] when When F ≥ 1, the crack propagates along the direction of maximum F; when When the value is less than 1, the crack tip remains stable. If the rock mass is approximately isotropic, the β dimension degenerates to a constant; if the temperature effect is insignificant, the T dimension degenerates to a constant.

[0144] Step 9, Engineering Stage Update: During stage transition, the boundary and fracture network are reconstructed, and the displacement discontinuities, fracture aperture, fluid pressure, temperature, phase state, and phase transition pressure history of the previous stage are inherited.

[0145] When the engineering phase moves from the excavation phase to the cooling phase, or from the cooling phase to the water injection / drainage phase, this invention does not re-initialize the model. Instead, it uses the fracture geometry, displacement discontinuity, fracture aperture, temperature field, fluid pressure, phase state identifier, and phase transition pressure already formed in the previous phase as the initial conditions for the next phase.

[0146] The boundary of the excavated borehole is divided into m DD elements. First, using the displacement discontinuities of each element, the stress and displacement at the boundary between the rock mass and the borehole can be obtained. Second, the subsequent excavation of the second borehole is considered. The boundary of the second borehole is divided into (km) elements, making the total number of elements for both boreholes k. Special treatment is needed for the boundary conditions in the dual-borehole model. The stress applied to the boundary of the first borehole is set to zero, while the stress applied to the boundary of the second borehole is the combined stress from the first step at the "future" boundary location of the second borehole. These are the sum of the in-situ stress and induced stress generated by the excavation of borehole 1. Finally, the sum of the two steps yields the displacement discontinuities for the next stage.

[0147] (i=m+1,k);

[0148] (i=m+1,k);

[0149] For a sequential excavation problem with n excavation steps, an equal number of computational steps are used for simulation. The model boundary is updated in each computational step, and the results of this step and the previous step (i.e., DD values) are stored and used in the next step.

[0150] ;

[0151] in, Represents the set of historical variables. and These represent the computational domain and boundary set of adjacent engineering stages, respectively. This represents variable mapping and new boundary initialization operations. Through the above processing, the present invention can simulate the continuous evolution of fractured rock mass from excavation disturbance, thermal disturbance, hydraulic disturbance to low-temperature phase transformation disturbance throughout the entire construction and operation process.

[0152] Step 10: Output Engineering Evaluation Results: Repeat the above steps until the set engineering time or convergence condition is reached, generating crack propagation paths, phase transformation damage zones, ice expansion-controlled crack segments, seepage channel reconstruction zones, and risk classification maps. When F ≥ 1, the crack propagates along the direction of maximum F; when When the value is less than 1, the crack tip remains stable. If the rock mass is approximately isotropic, the β dimension degenerates to a constant; if the temperature effect is insignificant, the T dimension degenerates to a constant.

[0153] Based on the above method, this invention also proposes a phase transformation-hydrothermal coupled crack propagation simulation system for fractured rock masses in low-temperature underground engineering. The system executes the above method and includes a fractured rock mass engineering stage model module, a temperature field update module, a phase transformation identification module, an equivalent phase transformation load module, a hydraulic update module, a hydrothermal exchange module, a temperature-direction parameter correction module, a crack propagation judgment module, an engineering stage update module, and a result output module.

[0154] The present invention will be further described below with reference to specific embodiments.

[0155] Example 1:

[0156] Stability Simulation of Fractured Rock Mass in Low-Temperature Underground Storage: Taking the surrounding rock of a low-temperature underground storage facility as an example, several initial fractures exist in the surrounding rock. The project goes through four stages: excavation, cooling, water injection and drainage, and operation. First, a two-dimensional or quasi-three-dimensional calculation model is established, and the elastic modulus of the surrounding rock, Poisson's ratio, fracture normal stiffness, fracture tangential stiffness, initial in-situ stress, initial temperature, boundary temperature, fracture water pressure, saturation, ice elastic modulus, and phase transition temperature range are input.

[0157] During the excavation phase, the opening or closing of fractures caused by excavation unloading is calculated, and the fracture aperture and displacement discontinuities are stored. During the cooling phase, the temperature field is calculated and the temperature is mapped to each fracture element. When the temperature of a fracture element is below the phase transition threshold, it is identified as a transitional phase or an ice phase, and the phase transition participation coefficient is calculated.

[0158] Within each time step, the equivalent stiffness and equivalent phase change pressure of the ice layer are calculated for the ice phase element. The equivalent phase change pressure and water pressure are used together as the fracture normal boundary conditions to solve for fracture deformation. After fracture deformation, the fracture aperture is updated, and the fracture aperture, in turn, updates the equivalent stiffness of the ice layer, fluid channels, and phase change pressure.

[0159] Subsequently, the propagation of multiple candidate directions at the crack tip is evaluated. The fracture toughness of each candidate direction is given by a temperature-direction bivariate parameter library. When the comprehensive propagation index of a candidate direction reaches a threshold, a new fracture element is added along that direction, and the fracture network and hydraulic connectivity are updated.

[0160] Finally, based on the crack propagation path, the continuous distribution of phase transformation units, and the changes in fluid channels, the low-temperature phase transformation damage zone, the ice expansion-controlled crack segment, and the recommended key monitoring areas are output.

[0161] Example 2:

[0162] Phased Excavation and Freeze-Thaw Disturbance Simulation of Tunnels in Cold Regions: In tunnel engineering in cold regions, multiple engineering phases are established according to the advancement sequence of the tunnel face. After each phase is completed, the crack aperture, displacement discontinuity, temperature, fluid pressure, and phase characteristics formed in the previous phase are mapped to the next phase. For newly excavated boundaries, their initial stress release state is calculated first, and then cryogenic and hydraulic boundaries are applied.

[0163] When groundwater freezes in the fissures, the fissure boundary load is updated according to the equivalent phase transition normal load model of this invention; as the temperature rises, the phase transition participation coefficient decreases, and the equivalent stiffness of the ice layer and the phase transition pressure are gradually released. This can simulate fissure propagation, local spalling, and the development of loosened surrounding rock zones under freeze-thaw cycles.

[0164] Example 3:

[0165] Parameter library construction method: The temperature-direction bivariate parameter library can be obtained from indoor experiments, field inversion, or calibration of initial parameter values ​​from literature. For each temperature node Tj and direction node β, record... , , , and During calculation, adjacent nodes are read and interpolated based on the crack tip temperature and candidate orientation. If the rock mass is approximately isotropic, the orientation dimension can be degenerated to a constant; if the temperature effect is not significant, the temperature dimension can be degenerated to a constant.

[0166] Any parts not mentioned in this invention can be achieved by referring to existing technologies.

[0167] Those skilled in the art should recognize that the above embodiments are merely illustrative of this application and are not intended to limit this application. Any appropriate changes and variations made to the above embodiments within the essential spirit and scope of this application should fall within the scope of protection of the claims of this application.

Claims

1. A method for simulating crack propagation in fractured rock masses using a phase transformation-hydrothermal coupling approach for low-temperature underground engineering, characterized in that... The steps are as follows: S1. Establish a fractured rock mass engineering stage model that includes surrounding rock boundary, excavation boundary, existing fractures, candidate crack ends, thermal boundary, hydraulic boundary, mechanical boundary, and at least one engineering stage sequence. S2. Discretize the fractures, boundaries and crack ends in the fractured rock mass engineering stage model into computational units, and establish state variables for each computational unit. Use each computational unit as a common carrier for mechanical units, hydraulic units, thermal and phase change state units. S3. At each engineering stage and time step, phase transition identification is performed on the computational unit according to the temperature field and water content, and the unit is divided into aqueous phase unit, transition phase unit or ice phase unit. S4. Calculate the equivalent phase transformation normal load that changes with the fracture normal deformation feedback for the water phase element and the ice phase element, and input the equivalent phase transformation normal load together with water pressure, in-situ stress and external boundary load into the fracture mechanics solution. S5. Based on the temperature and crack orientation of the calculation unit, tensile strength, cohesion, internal friction angle, type I fracture toughness and type II fracture toughness are obtained or interpolated from the temperature-direction bivariate parameter library. S6. Based on the equivalent phase transformation normal load, fluid pressure and temperature-corrected strength parameters, calculate crack deformation, crack initiation index and crack propagation index, and update crack aperture and fluid channel. S7. When the stage transition conditions are met, reconstruct the excavation boundary or thermal / hydraulic boundary, and map the displacement discontinuity, fluid pressure, temperature, phase identification and crack aperture of the previous engineering stage to the next engineering stage. S8. Repeat S3 to S7 until the set engineering time or convergence condition is reached, and output the crack propagation path, phase transformation damage zone and seepage channel reconstruction zone.

2. The method for simulating crack propagation in fractured rock mass coupled with hydrothermal coupling for low-temperature underground engineering as described in claim 1, characterized in that: In step S2, the displacement discontinuity method is used to discretize the cracks of arbitrary shapes; in the crack element, the tangential stress and normal stress are obtained by superimposing the boundary influence coefficient and the displacement discontinuity.

3. The method for simulating crack propagation in fractured rock mass coupled with hydrothermal coupling for low-temperature underground engineering as described in claim 1, characterized in that: In step S3, the aqueous phase unit has a temperature greater than 0°C; the transition phase unit has a temperature between -2.5°C and 0°C; and the ice phase unit has a temperature between -2.5°C and -2.5°C.

4. The method for simulating crack propagation in fractured rock mass coupled with hydrothermal coupling for low-temperature underground engineering as described in claim 1, characterized in that: In step S4, the equivalent phase transformation normal load is determined by the initial phase transformation pressure, the equivalent stiffness of the ice layer, the fracture contact stiffness, and the fracture normal displacement increment, so that the phase transformation pressure is corrected by feedback as the fracture opens or closes; for the j-th fracture element, the calculation formula for the ice expansion pressure is shown in equation (1): （1）； In formula (1): The initial ice expansion pressure is for unjointed deformation; For the pressure of ice expansion, The equivalent stiffness of the ice layer; This represents the increment of the normal displacement of the crack.

5. The method for simulating phase transformation-hydrothermal coupled crack propagation in fractured rock masses for low-temperature underground engineering as described in claim 4, characterized in that: The formula for calculating ice expansion pressure is shown in equation (2), and the formula for calculating the equivalent stiffness of the ice layer is shown in equation (3): （2）； （3）； In equation (2), This represents the elastic modulus during the initial stage of freezing. In equation (3), e is the hydraulic opening; The formula for calculating the contact stress in fracture mechanics is shown in equation (4): （4）； In equation (4): The normal contact stiffness of the crack; This represents the initial displacement in the crack normal direction; The formula for calculating the total fracture ice expansion pressure is shown in equation (5): = （5）。 6. The method for simulating crack propagation in fractured rock mass coupled with hydrothermal coupling, as described in claim 4, is characterized in that: The equivalent stiffness of the ice layer is determined based on the elastic modulus of the ice body, the hydraulic aperture of the fissures, the phase transformation participation coefficient, and the minimum aperture limit.

7. The method for simulating phase transformation-hydrothermal coupled crack propagation in fractured rock masses for low-temperature underground engineering as described in claim 1, characterized in that: In step S5, the temperature-direction bivariate parameter library includes tensile strength, cohesion, internal friction angle, Type I fracture toughness, and Type II fracture toughness recorded by temperature nodes and azimuth nodes, and the non-nodal temperature and azimuth are generated by interpolation to generate calculation parameters; in step S6, the crack propagation index includes Type I propagation component and Type II propagation component, which are normalized by temperature-direction corrected Type I fracture toughness and Type II fracture toughness, respectively.

8. The method for simulating crack propagation in fractured rock mass coupled with hydrothermal coupling, as described in claim 7, is characterized in that: Step S6 further includes a fracture hydraulic renewal step, a fracture water-heat exchange step, and a damping update step based on changes in fracture aperture. The fracture hydraulic renewal step calculates the inter-domain flow rate based on the average aperture, pressure difference, fluid viscosity, and unit length of adjacent fracture domains, updates the fluid pressure and hydraulic connectivity based on the inter-domain flow rate, and calculates the equivalent leakage from the fracture to the surrounding rock based on the surrounding rock permeability coefficient, effective leakage distance, fracture domain pressure, and initial pore water pressure. The fracture water-heat exchange step calculates the heat carried by the fluid based on the flow rate entering the fracture domain, fluid density, fluid specific heat, upstream fracture domain temperature, and current fracture domain temperature, and uses it as an equivalent heat source or cold source for the fracture wall to participate in the temperature field update. The steps for updating the damping of crack aperture changes include updating the crack aperture, fluid pressure, temperature and phase change pressure step by step according to the damping coefficient when the crack aperture change caused by a single iteration exceeds a preset ratio.

9. The method for simulating crack propagation in fractured rock mass coupled with hydrothermal coupling, as described in claim 7, is characterized in that: In step S7, the stage transition conditions include the completion of the excavation step, the change of water injection pressure, the change of drainage boundary, the change of temperature boundary, the change of insulation layer parameters, or the achievement of the set crack propagation length; in step S8, the output phase transformation damage zone includes the continuous distribution zone of ice phase units, the concentrated zone of transition phase units, the crack segment where the phase transformation normal load exceeds the threshold, and the candidate crack tip where the crack propagation index exceeds the threshold.

10. A simulation system for phase transformation-hydrothermal coupled crack propagation in fractured rock masses for low-temperature underground engineering, characterized in that: The system described herein executes the method described in any one of claims 1 to 9, and the system comprises a fractured rock mass engineering stage model module, a temperature field update module, a phase transformation identification module, an equivalent phase transformation load module, a hydraulic update module, a hydrothermal exchange module, a temperature-direction parameter correction module, a crack propagation determination module, an engineering stage update module, and a result output module.

Images