Three-dimensional rough fractured rock mass seepage heat transfer numerical simulation method and device

Abstract

The invention discloses a numerical simulation method and device for seepage heat transfer of a three-dimensional rough fractured rock mass, and relates to the technical field of rock mechanics engineering.The dynamic opening degree of the fractured rock mass is obtained by combining the initial opening degree and the equivalent opening degree of the fractured rock mass, and in the process of obtaining the dynamic opening degree of the fractured rock mass, the dynamic opening degree of the fractured rock mass is obtained. According to the method, based on the initial opening and the equivalent opening, the coupling influence of the shear expansion effect and roughness degradation on the fracture opening under the action of shear stress is quantified, and then based on the dynamic opening of the fractured rock mass, the point cloud data of the upper and lower fracture surfaces of the fractured rock mass after shearing is updated; the dynamic updating of the three-dimensional geometrical morphology of the real fracture in the shearing process is realized, and the bidirectional coupling simulation of the seepage field and the temperature field is realized based on the three-dimensional geometrical morphology data of the real fracture, so that the precision of the seepage heat transfer numerical simulation of the fractured rock mass is improved.

Description

Technical Field

[0001] The present invention relates to the field of rock mechanics engineering technology, and in particular to a three-dimensional rough fractured rock mass seepage heat transfer numerical simulation method, device, equipment and medium. Background Art

[0002] In engineering fields such as geothermal energy development and geological disposal of nuclear waste, the seepage-heat transfer coupling behavior of fractured rock masses under shear stress has a significant impact on system performance. Studies have shown that the shear strength of fractures increases from room temperature to a critical value of 200°C and then decreases with increasing temperature. The thermal hardening of fracture shear strength may be caused by the removal of adsorbed water, crystals, or thermal expansion of the rock matrix. The thermal weakening of fracture shear strength is mainly caused by the generation of microcracks in the rock matrix and the degradation of fracture roughness. In addition, a few studies have conducted shear behavior of rock fractures after heat treatment and cooling in air or water. With increasing heating temperature, the hardness, basic friction angle, shear stiffness, and shear strength of rock fractures change significantly.

[0003] At present, the numerical simulation of the seepage-heat transfer coupling process of fractured rock mass under shear stress mainly focuses on establishing a simplified theoretical model of shear-seepage-heat transfer coupling; researchers have successfully implemented geological forces in geothermal energy systems by applying water pressure below the minimum principal stress to achieve large-scale water injection in EGS to hydraulically stimulate existing fractures, induce the generation of shear expansion, and perform numerical simulation in this process. The main simulation methods are equivalent continuous medium model, dual-medium model and discrete fracture network model; the equivalent continuous medium model regards the fractured porous medium as a continuous medium, and characterizes the influence of fractures on fluid flow, heat transfer and mechanical response through equivalent parameters; the dual-medium model divides the fractured rock mass into fracture area and rock matrix area, and calculates the mass and energy transfer between the two superimposed areas; the discrete fracture network model can clearly characterize the distribution of fractures, and can capture the thermal-hydro-mechanical coupling behavior of fractures in the formation and their influence on fluid flow and heat transfer.

[0004] However, the three simulation methods, namely the equivalent continuous medium model, the dual-medium model and the discrete fracture network model, normalize the size and geometry of the fracture system in the rock mass on a two-dimensional level. It is difficult to characterize the influence of the three-dimensional geometric morphology of the real fractures on the seepage-heat transfer coupling process during the shear process, resulting in low accuracy in the numerical simulation of seepage heat transfer in fractured rock masses. Summary of the Invention

[0005] The embodiment of the present invention provides a three-dimensional rough fractured rock mass seepage heat transfer numerical simulation method and device, which can solve the problem in the prior art that the current numerical simulation method is to normalize the size and geometric shape of the fracture system in the rock mass on a two-dimensional level, making it difficult to characterize the influence of the three-dimensional geometric morphology of the real fracture during the shear process on the seepage-heat transfer coupling process, resulting in low accuracy in the numerical simulation of the seepage heat transfer in the fractured rock mass.

[0006] An embodiment of the present invention provides a method for numerically simulating seepage heat transfer in a three-dimensional rough fractured rock mass, comprising the following steps: Obtaining point cloud data of the fracture surface of the fractured rock mass before and after shearing; wherein the point cloud data of the fracture surface includes point cloud data of the upper fracture surface and the lower fracture surface; Match the point cloud data of the upper and lower fracture surfaces after shearing, discretize the matched point cloud data of the upper and lower fracture surfaces into grid units of the same size, and obtain the initial opening of the corresponding grid unit through the elevation difference between the upper and lower fracture surfaces within the grid unit, as well as the initial opening of the fractured rock mass to represent the shear expansion effect; Using the elevation data of the upper and lower fracture surfaces of the fractured rock mass before and after shearing, the elevation changes of the upper and lower fracture surfaces before and after shearing are respectively obtained; based on the elevation changes of the upper and lower fracture surfaces before and after shearing, the equivalent opening of the fractured rock mass formed by roughness degradation is obtained; The dynamic opening of the fractured rock mass is obtained based on the initial opening and equivalent opening of the fractured rock mass; the point cloud data of the upper and lower fracture surfaces of the fractured rock mass after shearing are updated based on the dynamic opening of the fractured rock mass, and the geometric morphology data of the fractured rock mass after shearing are obtained using the updated point cloud data of the upper and lower fracture surfaces; The seepage and heat transfer coupling process of fractured rock mass is simulated by using the geometric morphology data of fractured rock mass after shearing and based on the seepage field and temperature field of fractured rock mass.

[0007] Preferably, the step of obtaining point cloud data of the fracture surface of the fractured rock mass before and after shearing includes: Granite specimens with a single crack were prepared using the splitting method, and the initial crack surface roughness was controlled to be in the range of JRC 5-15; A VTOP 300T 3D laser scanner was used to obtain 3D point cloud data of the crack surface before shearing with a resolution of ≤10μm. The granite samples were subjected to a three-stage heat treatment in a high-temperature furnace, followed by heating, insulation, and natural cooling to room temperature to simulate a real geothermal environment; Compression shear tests were carried out on heat-treated granite specimens under the conditions of normal stress of 5-20 MPa and shear rate of 0.1 mm / min. The fracture surface after shearing was scanned by a VTOP 300T three-dimensional laser scanner to obtain the point cloud data of the fracture surface of the sheared fractured rock mass.

[0008] Preferably, obtaining the initial opening of the fractured rock mass to characterize the shear expansion effect includes: The interpolation function is used to perform spatial difference calculation on the point cloud data of the upper and lower crack surfaces after shearing, and the point cloud data of the upper and lower crack surfaces after shearing are matched based on the spatial difference of the point cloud data of the upper and lower crack surfaces; The point cloud data of the upper and lower crack surfaces after matching are discretized into grid cells of the same size. The initial opening of the grid cell is obtained by the elevation difference between the upper and lower crack surfaces in the corresponding grid cell, which is expressed as: ; in: e I Represents grid points I The opening of the position; and Respectively represent the upper and lower crack surfaces on the grid I The elevation of the place; According to the initial openings of all grid cells, the initial openings of the fractured rock mass representing the shear dilation effect are obtained.

[0009] Preferably, obtaining the equivalent opening of the fractured rock mass formed by roughness degradation includes: The point cloud data of the fracture surface of the fractured rock mass before and after shearing are transformed into the same coordinate system by coordinate transformation. The elevation data of the upper and lower fracture surfaces of the fractured rock mass before and after shearing are used to obtain the elevation changes of the upper and lower fracture surfaces before and after shearing, respectively, which can be expressed as: ; ; According to the elevation changes of the upper and lower fracture surfaces before and after shearing, the equivalent opening of the fractured rock mass formed by roughness degradation is obtained. e d , expressed as: ; in: 、 、 、 Represent the mesh points of the upper and lower crack surfaces before and after shearing, respectively I Elevation changes at the site.

[0010] Preferably, the dynamic aperture of the fractured rock mass is obtained by: ; in: e m Indicates the dynamic opening of the fractured rock mass; Δ un Indicates the change in mechanical opening caused by the closure of the normal load; Δ u s Represents the change in mechanical opening caused by shear expansion.

[0011] Preferably, the method of simulating seepage heat transfer in fractured rock mass includes: Import the XYZ coordinate data of the fracture surface in the geometric morphology data of the fractured rock mass after shearing into COMSOL™ software to generate an unstructured mesh; The fluid velocity field of fractured rock mass is simulated based on the Navier-Stokes equation, and the equivalent aperture correction factor is introduced. e d Dynamically adjust local permeability; Based on the local heat balance equation and combined with the correction coefficient of thermal conductivity in the contact area, the thermal resistance effect of the contact surface of the fractured rock mass is calculated to simulate the temperature field of the fractured rock mass. Through the multi-physics coupling of COMSOL™ software, based on the simulated fluid velocity and temperature fields, the velocity field data is transferred to the temperature field equation in real time to simulate the seepage heat transfer in fractured rock mass.

[0012] The embodiment of the present invention further provides a three-dimensional rough fractured rock mass seepage heat transfer numerical simulation device, comprising: A data module is used to obtain point cloud data of the fracture surface of the fractured rock mass before and after shearing; wherein the point cloud data of the fracture surface includes point cloud data of the upper fracture surface and the lower fracture surface; The characterization module is used to match the point cloud data of the upper and lower fracture surfaces after shearing, discretize the matched point cloud data of the upper and lower fracture surfaces into grid units of the same size, obtain the initial opening of the corresponding grid unit through the elevation difference between the upper and lower fracture surfaces within the grid unit, and obtain the initial opening of the fractured rock mass to characterize the shear expansion effect; Using the elevation data of the upper and lower fracture surfaces of the fractured rock mass before and after shearing, the elevation changes of the upper and lower fracture surfaces before and after shearing are respectively obtained; based on the elevation changes of the upper and lower fracture surfaces before and after shearing, the equivalent opening of the fractured rock mass formed by roughness degradation is obtained; The simulation module is used to obtain the dynamic aperture of the fractured rock mass based on the initial aperture and equivalent aperture of the fractured rock mass; based on the dynamic aperture of the fractured rock mass, the point cloud data of the upper and lower fracture surfaces of the fractured rock mass after shearing are updated, and the geometric morphology data of the fractured rock mass after shearing are obtained using the updated point cloud data of the upper and lower fracture surfaces; The seepage and heat transfer coupling process of fractured rock mass is simulated by using the geometric morphology data of fractured rock mass after shearing and based on the seepage field and temperature field of fractured rock mass.

[0013] An embodiment of the present invention further provides an electronic device, including a memory and a processor; The memory is used to store computer programs; The processor is configured to implement the steps of the above-mentioned method for numerical simulation of seepage heat transfer in a three-dimensional rough fractured rock mass when executing the computer program stored in the memory.

[0014] An embodiment of the present invention further provides a computer-readable storage medium for storing a computer program, which, when executed by a processor, implements the steps of the above-mentioned method for numerical simulation of seepage heat transfer in a three-dimensional rough fractured rock mass.

[0015] The embodiments of the present invention provide a method and apparatus for numerically simulating seepage heat transfer in a three-dimensional rough fractured rock mass. Compared with the prior art, the method and apparatus have the following beneficial effects: The present invention discretizes the point cloud data of the upper and lower fracture surfaces after matching into grid units of the same size, obtains the initial opening of the corresponding grid unit and the initial opening of the fractured rock mass through the elevation difference between the upper and lower fracture surfaces in the grid unit; obtains the equivalent opening of the fractured rock mass according to the elevation change of the upper and lower fracture surfaces before and after shearing, and obtains the dynamic opening of the fractured rock mass by combining the initial opening and equivalent opening of the fractured rock mass. In the process of obtaining the dynamic opening of the fractured rock mass, the present invention quantifies the coupling influence of the shear expansion effect and roughness degradation on the fracture opening under the action of shear stress based on the initial opening and the equivalent opening, and then updates the point cloud data of the upper and lower fracture surfaces of the fractured rock mass after shearing based on the dynamic opening of the fractured rock mass, so as to realize the dynamic update of the three-dimensional geometric morphology of the real fracture during the shearing process, and realizes the two-way coupling simulation of the seepage field and the temperature field based on the three-dimensional geometric morphology data of the real fracture, thereby improving the accuracy of the numerical simulation of seepage heat transfer in the fractured rock mass. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 A schematic diagram of the overall process of a three-dimensional rough fractured rock mass seepage heat transfer numerical simulation method provided by an embodiment of the present invention; Figure 2 A schematic diagram of an experimental process for the compression-shear characteristics of granite fractures after high-temperature treatment, according to a method for numerical simulation of seepage heat transfer in a three-dimensional rough fractured rock mass provided by an embodiment of the present invention; Figure 3 A schematic diagram of the calculation process of crack opening change and roughness degradation during compression and shearing of a three-dimensional rough fractured rock mass seepage heat transfer numerical simulation method provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0017] To make the above-mentioned objects, features, and advantages of the present invention more readily apparent, specific embodiments of the present invention are described in detail below with reference to the accompanying drawings. The following description sets forth numerous specific details to facilitate a full understanding of the present invention. However, the present invention can be implemented in many other ways than those described herein, and those skilled in the art may make similar modifications without departing from the scope of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0018] See also Figure 1 The embodiment of the present invention provides a three-dimensional numerical simulation method for seepage and heat transfer in rough fractured rock mass. At present, there are mainly two methods for numerical simulation of seepage and heat transfer coupling in fractured rock mass under shear stress.

[0019] The first approach is to establish a simplified theoretical model of shear-seepage-heat transfer coupling; with geothermal energy extraction as the background, scholars have successfully implemented geological forces in geothermal energy systems by applying water pressure below the minimum principal stress to achieve large-scale water injection in EGS to stimulate existing fractures through hydraulics, thereby inducing the generation of shear expansion; the numerical methods used to simulate these problems are mainly divided into equivalent continuum models, dual-medium models and discrete fracture network models; the equivalent continuum model regards fractured porous media as continuous media and characterizes the effects of fractures on fluid flow, heat transfer and mechanical response through equivalent parameters; the continuum model has the advantage of high computational efficiency , but cannot fully consider the influence of fracture connectivity and other fracture properties on thermal performance; the dual-medium model divides the fractured rock mass into a fracture region and a rock matrix region, and calculates the mass and energy transfer between the two superimposed regions; however, these methods over-normalize the size and geometry of the fracture system and cannot fully represent the anisotropy and geometric characteristics of the fractures; the discrete fracture network model can clearly characterize the distribution of fractures and capture the thermal-hydro-mechanical coupling behavior of fractures in the formation and their influence on fluid flow and heat transfer, but it still cannot embed the influence of the three-dimensional geometric morphology of real fractures on the seepage-heat transfer coupling process during shear.

[0020] The second method is to decompose the research steps and conduct a numerical study on the seepage-material transport characteristics inside the fracture under shear stress; at present, there are few reports on the research on the seepage-heat transfer characteristics inside rough fractures under shear stress. The most representative one that can be found is the numerical study on the influence of fracture aperture anisotropy on the thermal properties of shear-induced fractures conducted by Esuru Rita Okoroafor in 2022. The research focus is on the influence of the relationship between flow direction and shear direction on the heat transfer inside the fracture; therefore, this method is summarized by drawing on the research ideas of the influence of shear displacement on fluid flow and solute transport in three-dimensional rough fractures; since the solute transport process in the fracture is highly dependent on the fluid flow in the pore space (ie, opening), the solute transport process in fractures with different shear displacements may change with significant changes in the fluid flow path; Zhou determined the shear induced fracture by using the normal displacement measured in the indoor shear test and the digitized fracture surface. The changes in the crack opening caused by the fracture were studied, and the fluid flow in the cracks under different shear displacements was simulated by solving the Reynolds equation. The results showed that for matched cracks, due to shear expansion, the increase in shear displacement will enhance the flow and transport; Zhao's study found that shear-induced roughness degradation significantly increased the value of the solute retardation coefficient by providing more adsorption surfaces in the crack pores; Vilarrasa et al. found that shear displacement increased the flow velocity perpendicular to the shear direction and promoted the solute migration in the cracks; it is not difficult to find that a series of studies on seepage transport inside cracks under shear action conducted by the second method have attracted widespread attention, and also proved the feasibility of using this method to study the multi-field coupling process of rough cracks.

[0021] Based on the above analysis, it can be seen that the current modeling of fracture seepage-heat transfer coupling process has the following defects: One is that the current two-dimensional modeling cannot reflect the dynamic evolution of the actual crack morphology; specifically, the three-dimensional rough cracks are simplified into two-dimensional planes or equivalent openings, resulting in excessive smoothing of three-dimensional geometric changes such as convex fractures and particle migration during the shear process; the reason is that the two-dimensional model cannot represent the asymmetric wear of the crack surface along the shear direction, the three-dimensional movement path of particles (such as rolling and embedding) and the local opening heterogeneity, which will lead to significant deviations in the prediction of seepage channels.

[0022] One is that the current thermal-fluid-solid coupling accuracy is low; specifically, when the heat transfer equation is coupled with the seepage field, the local thermal resistance effect of the roughness degradation of the crack surface and the anisotropic heat conduction caused by particle accumulation are not considered; the reason is that the two-dimensional geometric simplification leads to the weakening of the spatial correlation between the temperature field and the seepage field, and it is impossible to capture the true distribution of the thermal boundary layer in the three-dimensional crack.

[0023] One is that the current means are highly dependent on empirical constitutive models and lack universality; specifically, a large number of experiments are required to calibrate the shear-openness relationship parameters, which makes it difficult to extend to complex geological conditions (such as high-temperature dissolution and chemical-mechanical coupling); the reason is that the two-dimensional model lacks physical driving modeling of the microscopic mechanical behavior of the fracture surface (such as particle crushing and dissolution) and only relies on macroscopic statistical corrections.

[0024] To address the problem that two-dimensional models cannot reflect the three-dimensional asymmetric wear, particle migration, and local aperture heterogeneity of the crack surface during shearing, the present invention accurately models the dynamic evolution of the three-dimensional rough crack geometry. The specific process is as follows: ①Introduction of three-dimensional crack morphology dynamic reconstruction technology: Based on high temperature compression shear test and 3D scanning (VTOP 300T), the shear dilatancy effect (Δ u s ) and roughness degradation ( e d ) coupled mathematical model to quantify the normal closure (Δ u n ) and shear expansion ( e d ) on the joint effect of crack aperture.

[0025] Through point cloud meshing and contact wear criterion ( e I <0 is the contact area, e I ≥0 for seepage channels), and dynamically update the geometric morphology of the fracture surface.

[0026] ②High temperature-shear coupling experimental data drive: Through the three-stage high-temperature treatment of "heating-insulation-cooling" to simulate the real geothermal environment, combined with compression-shear experiments and three-dimensional morphology scanning after shearing, the correlation data of temperature-mechanics-geometry evolution are obtained, providing a physical verification basis for the model.

[0027] Aiming at the problem that the two-dimensional simplified model ignores the local thermal resistance effect and seepage anisotropy caused by crack roughness degradation, the present invention conducts a refined modeling of heat-fluid-solid multi-field coupling; the specific process is as follows: ①Multi-field coupling algorithm based on real three-dimensional morphology: The dynamically updated fracture geometry model (XYZ coordinate data) is imported into COMSOL™, and the fluid mechanics (Navier-Stokes equations) and heat transfer equations (local thermal equilibrium theory) are coupled through the secondary development module to achieve two-way interactive calculations of the seepage field and temperature field.

[0028] An equivalent opening correction factor is embedded in the fluid module to reflect the blockage or expansion of the seepage path caused by particle migration; dynamic calculation of contact thermal resistance is introduced in the heat transfer module to capture the heat conduction inhibition effect in the rough surface contact area.

[0029] ② Chemical coupling expansion reserved interface: A dissolution / precipitation module interface is reserved in the model to support future integrated chemical processes (such as the dissolution rate equation), avoiding the lack of scalability of existing technologies due to their closed architecture.

[0030] In order to solve the problem that the current two-dimensional experimental calibration cannot verify the internal dynamic behavior of three-dimensional cracks, the present invention conducts model verification; the specific process is as follows: ① Direct verification: geometric morphology comparison.

[0031] The contact distribution ( e I <0 area) to perform spatial matching analysis to verify the reliability of the geometric reconstruction algorithm.

[0032] ②Indirect verification: seepage-heat transfer coupling verification.

[0033] Combine shear-seepage experiments (flow rate, pressure gradient) with numerical simulation results to compare the seepage field distribution; verify the coupling accuracy of the temperature field and seepage field through infrared thermal imaging or distributed temperature sensing (DTS) data.

[0034] A multi-scale verification strategy was adopted: microscale (local contact thermal resistance), mesoscale (heterogeneity of seepage channels), and macroscale (overall seepage heat transfer efficiency).

[0035] Specific implementation steps include: Step 1: Dynamic reconstruction of 3D rough crack geometry, e.g. Figure 2 shown.

[0036] 1. High temperature compression and shear test and 3D topography data acquisition; including: Specimen preparation: Single-crack granite specimens were prepared by the splitting method to control the initial roughness (JRC 5-15).

[0037] Initial morphology scanning: Use a VTOP 300T scanner to obtain initial fracture 3D point cloud data (resolution ≤ 10 μm); the initial fracture 3D point cloud data includes the 3D point cloud data of the upper fracture surface and the lower fracture surface of the fractured rock mass.

[0038] High temperature treatment: The sample was treated in a high temperature furnace in three stages: heating (rate 5°C / min) → holding (target temperature 200-500°C, 2h) → cooling (natural cooling to room temperature)

[0039] Shear experiment: Shearing was performed using a compression shear testing machine (normal stress 5-20 MPa, shear rate 0.1 mm / min), and the shear displacement-stress curve was recorded.

[0040] Post-shear morphology scanning: Perform a final 3D scan of the fracture surface after shearing; wherein, the fractured rock mass after shearing includes an upper fracture surface and a lower fracture surface. After the final 3D scan of the fractured rock mass after shearing, 3D point cloud data of the upper fracture surface and the lower fracture surface of the fractured rock mass after shearing can be obtained.

[0041] 2. Crack opening and roughness degradation modeling, such as Figure 3 shown.

[0042] Establish a mathematical model, expressed as: .

[0043] in: represents the initial opening under the initial effective stress; Δ u n Indicates the change in mechanical opening caused by the closure of the normal load; Δ u s Represents the change in mechanical opening caused by shear expansion.

[0044] The dynamic calculation of opening is expressed as: .

[0045] in: e d Represents the equivalent openness resulting from roughness degradation, evaluated by assuming a uniform distribution of the total amount of variation among the void grid points.

[0046] 3. Point cloud data processing and geometric reconstruction; including: Meshing processing: discretize the upper and lower crack surface point cloud data into grid units of the same size (side length ≤ 0.1 mm).

[0047] Contact criteria: If the grid cell opening e I ≥0, marked as percolation channel.

[0048] like e I <0, marked as contact wear area.

[0049] Dynamic update: The crack morphology data at each stage of shear displacement is iteratively updated through the MATLAB program.

[0050] According to the above formula, first make the initial opening ei The solution.

[0051] In order to achieve precise matching of the upper and lower crack surface points, the interpolation function is used to perform spatial difference calculation on the three-dimensional point cloud data obtained by scanning; due to the irregularity and complexity of the crack surface, the three-dimensional point cloud data directly obtained by scanning cannot be directly used for subsequent calculations and analysis; the interpolation function can perform spatial interpolation calculation on the point cloud data, fill in the gaps in the data, and improve the density and accuracy of the data, thereby achieving precise matching of the upper and lower crack surface points, that is, unifying the positions of the point cloud data corresponding to the upper and lower crack surfaces.

[0052] Among them, first, based on the naturally existing geometric features of the crack surface (such as edges and grooves), the FPFH (Fast Point Feature Histogram) feature descriptor is extracted; the RANSAC algorithm is used to screen the internal points, complete the initial transformation estimation, and control the registration error within twice the voxel resolution; the weighted ICP (Iterative Closest Point) algorithm is introduced, and the weight function is constructed by combining the normal vector angle and the point spacing; the interpolation kernel function bandwidth is automatically adjusted according to the change of the curvature of the crack surface, and the inverse distance weighted interpolation (IDW) is used in flat areas, and the thin plate spline interpolation is switched to in the curvature mutation area. The optimal parameter combination is determined through cross-validation to perform accurate matching of the upper and lower crack surface points.

[0053] The point cloud data of the upper and lower crack surfaces are processed by high-precision grids. The upper and lower crack surfaces are discretized into grid units of the same size. The crack opening at any grid point is determined by the height difference between the upper and lower crack surfaces, which can be expressed as: .

[0054] in: e I Represents grid points I The opening of the position; and Respectively represent the upper and lower crack surfaces on the grid I When e I When ≥0, it means there is a gap in the crack. e I <0, indicating that contact wear occurs during the shearing process.

[0055] The shearing process causes wear on the crack surface, resulting in a decrease in surface roughness. Therefore, the roughness degradation is converted into the change in the roughness height of the upper and lower crack surfaces after shearing. The crack sample morphology obtained before and after shearing is transformed into the same coordinate system through coordinate transformation. Then, the best fit in Geomagic Qualing (3D Systems) is used to obtain the morphology of the crack sample before and after shearing, and the change in the roughness height of the crack surface is calculated. The change in the height of the upper and lower crack surfaces is expressed as: .

[0056] .

[0057] It can be obtained that the equivalent opening formed by roughness degradation is e d for: .

[0058] in: 、 、 、 Represent the mesh points on the upper and lower surfaces of the crack before and after shearing, respectively I Elevation changes at the site.

[0059] Through the above formula, different shear displacements (U s ) stage to update the upper and lower crack surface point cloud data to obtain the geometric morphology data of the crack under shear stress.

[0060] Step 2: Force-flow-heat coupling.

[0061] 1. Model construction and equation coupling; including: Geometry model import: Import the updated fracture surface XYZ coordinate data into COMSOL™ and generate an unstructured mesh (minimum element size 0.05 mm).

[0062] Multiphysics coupling settings: ① Seepage field: Calculate the fluid velocity field based on the Navier-Stokes equation, and embed the equivalent opening correction factor e d To reflect the particle migration effect.

[0063] ② Temperature field: Using the local heat balance equation, a dynamic calculation module for contact thermal resistance is introduced (the thermal conductivity of the contact area is reduced by 50-70%).

[0064] ③ Bidirectional coupling: Flow velocity and temperature data are transmitted in real time through the COMSOL™ “Multiphysics Coupling” interface.

[0065] 2. Chemical coupling expansion interface is reserved.

[0066] An API interface is reserved in the model architecture to support subsequent integration of dissolution / precipitation modules (such as the dissolution rate model based on the Arrhenius equation).

[0067] Step 3: Model validation.

[0068] 1. Direct verification: geometric shape comparison.

[0069] The scanned data of the crack surface after shearing (resolution 20 μm) were compared with the contact distribution predicted by the model ( e I <0 and e I <0 region) were subjected to spatial coincidence analysis (error threshold ≤ 15%).

[0070] 2. Indirect verification: seepage-heat transfer data verification.

[0071] The flow-pressure gradient data were obtained through shear-seepage experiments and fitted with the simulation results (R² ≥ 0.90).

[0072] During implementation, the experimental parameters included: High temperature treatment conditions: heating rate 5°C / min, holding time 2h, temperature range 200-500°C.

[0073] Scanning accuracy: VTOP 300T scanner resolution ≤ 10 μm.

[0074] Model parameters include: Mesh size: 0.05-0.1 mm, contact area thermal conductivity correction factor 0.3-0.5.

[0075] Validation thresholds include: The geometric matching error is ≤15%, the seepage data R² is ≥0.90, and the temperature error is ≤10%.

[0076] The present invention is driven by the real morphology data of high temperature compression and shearing experiments, combined with the contact criterion ( e I <0) dynamically updates the fracture aperture to avoid over-smoothing of the two-dimensional equivalent model; through three-dimensional scanning and point cloud meshing, the modeling errors of geometric changes such as fracture surface roughness degradation and particle migration are reduced compared with the two-dimensional model, and the accuracy of seepage path prediction is improved.

[0077] The present invention directly couples the Navier-Stokes equation with the local thermal equilibrium theory, embeds the dynamic correction of contact thermal resistance (thermal conductivity reduced by 50-70%) and the equivalent opening factor ( e d ), truly reflecting the anisotropic heat flow interactions within rough fractures. This method can be expanded to chemical-mechanical coupling scenarios (with reserved dissolution / precipitation interfaces). This paper employs a dual-track verification system of "direct verification (geometric comparison) + indirect verification (seepage-heat transfer data)" to cover multi-scale accuracy requirements from microscopic contact to macroscopic heat transfer.

[0078] The present invention realizes the fully coupled dynamic simulation of shear, seepage and heat transfer in three-dimensional rough fractures for the first time, breaking through the geometric distortion bottleneck of two-dimensional models. The present invention is based on real experimental data and the direct solution of the Navier-Stokes equations, reducing the reliance on empirical constitutive models (such as Barton-Bandis) and improving universality. The present invention optimizes the multi-field coupling algorithm through COMSOL™ secondary development to meet the needs of real-time engineering analysis.

[0079] The above-described embodiments merely illustrate several implementations of the present invention, and while their descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the patent. It should be noted that a person skilled in the art would be able to make numerous variations and improvements without departing from the spirit of the present invention, all of which fall within the scope of protection of the present invention. Therefore, the scope of protection of the patent for this invention shall be determined by the appended claims.

Claims

1. A three-dimensional rough fractured rock mass seepage heat transfer numerical simulation method, characterized in that: The following steps are involved: Obtaining point cloud data of the fracture surface of the fractured rock mass before and after shearing; wherein the point cloud data of the fracture surface includes point cloud data of the upper fracture surface and the lower fracture surface; Match the point cloud data of the upper and lower fracture surfaces after shearing, discretize the matched point cloud data of the upper and lower fracture surfaces into grid units of the same size, and obtain the initial opening of the corresponding grid unit through the elevation difference between the upper and lower fracture surfaces within the grid unit, as well as the initial opening of the fractured rock mass to represent the shear expansion effect; Using the elevation data of the upper and lower fracture surfaces of the fractured rock mass before and after shearing, the elevation changes of the upper and lower fracture surfaces before and after shearing are respectively obtained; based on the elevation changes of the upper and lower fracture surfaces before and after shearing, the equivalent opening of the fractured rock mass formed by roughness degradation is obtained; The dynamic opening of the fractured rock mass is obtained based on the initial opening and equivalent opening of the fractured rock mass; the point cloud data of the upper and lower fracture surfaces of the fractured rock mass after shearing are updated based on the dynamic opening of the fractured rock mass, and the geometric morphology data of the fractured rock mass after shearing are obtained using the updated point cloud data of the upper and lower fracture surfaces; The seepage and heat transfer coupling process of fractured rock mass is simulated by using the geometric morphology data of fractured rock mass after shearing and based on the seepage field and temperature field of fractured rock mass.

2. A three-dimensional rough fractured rock mass seepage heat transfer numerical simulation method according to claim 1, characterized in that: The step of obtaining point cloud data of the fracture surface of the fractured rock mass before and after shearing includes: Granite specimens with a single crack were prepared using the splitting method, and the initial crack surface roughness was controlled to be in the range of JRC 5-15; A VTOP 300T 3D laser scanner was used to obtain 3D point cloud data of the crack surface before shearing with a resolution of ≤10μm. The granite samples were subjected to a three-stage heat treatment in a high-temperature furnace, followed by heating, insulation, and natural cooling to room temperature to simulate a real geothermal environment; Compression shear tests were carried out on heat-treated granite specimens under the conditions of normal stress of 5-20 MPa and shear rate of 0.1 mm / min. The fracture surface after shearing was scanned by a VTOP 300T three-dimensional laser scanner to obtain the point cloud data of the fracture surface of the sheared fractured rock mass.

3. The method for numerical simulation of seepage heat transfer in three-dimensional rough fractured rock mass according to claim 1, characterized in that: The acquisition of the initial opening of the fractured rock mass to characterize the shear expansion effect includes: The interpolation function is used to perform spatial difference calculation on the point cloud data of the upper and lower crack surfaces after shearing, and the point cloud data of the upper and lower crack surfaces after shearing are matched based on the spatial difference of the point cloud data of the upper and lower crack surfaces; The point cloud data of the upper and lower crack surfaces after matching are discretized into grid cells of the same size. The initial opening of the grid cell is obtained by the elevation difference between the upper and lower crack surfaces in the corresponding grid cell, which is expressed as: ； in: e I Represents grid points I The opening of the position; and Respectively represent the upper and lower crack surfaces on the grid I The elevation of the place; According to the initial openings of all grid cells, the initial openings of the fractured rock mass representing the shear dilation effect are obtained.

4. The method for numerical simulation of seepage heat transfer in three-dimensional rough fractured rock mass according to claim 3, characterized in that: The acquisition of the equivalent opening of the fractured rock mass formed by roughness degradation includes: The point cloud data of the fracture surface of the fractured rock mass before and after shearing are transformed into the same coordinate system by coordinate transformation. The elevation data of the upper and lower fracture surfaces of the fractured rock mass before and after shearing are used to obtain the elevation changes of the upper and lower fracture surfaces before and after shearing, respectively, which can be expressed as: ； ； According to the elevation changes of the upper and lower fracture surfaces before and after shearing, the equivalent opening of the fractured rock mass formed by roughness degradation is obtained. e d , expressed as: ； in: 、 、 、 Represent the mesh points of the upper and lower crack surfaces before and after shearing, respectively I Elevation changes at the site.

5. The method for numerical simulation of seepage heat transfer in three-dimensional rough fractured rock mass according to claim 4, characterized in that: The dynamic aperture of the fractured rock mass is obtained as follows: ； in: e m Indicates the dynamic opening of the fractured rock mass; Δ u n Indicates the change in mechanical opening caused by the closure of the normal load; Δ u s Represents the change in mechanical opening caused by shear expansion.

6. A three-dimensional rough fractured rock mass seepage heat transfer numerical simulation method according to claim 5, characterized in that: The method of simulating seepage heat transfer in fractured rock mass comprises: Import the XYZ coordinate data of the fracture surface in the geometric morphology data of the fractured rock mass after shearing into COMSOL™ software to generate an unstructured mesh; The fluid velocity field of fractured rock mass is simulated based on the Navier-Stokes equation, and the equivalent aperture correction factor is introduced. e d Dynamically adjust local permeability; Based on the local heat balance equation and combined with the correction coefficient of thermal conductivity in the contact area, the thermal resistance effect of the contact surface of the fractured rock mass is calculated to simulate the temperature field of the fractured rock mass. Through the multi-physics coupling of COMSOL™ software, based on the simulated fluid velocity and temperature fields, the velocity field data is transferred to the temperature field equation in real time to simulate the seepage heat transfer in fractured rock mass.

7. A three-dimensional rough fractured rock mass seepage heat transfer numerical simulation device, characterized in that: include: A data module is used to obtain point cloud data of the fracture surface of the fractured rock mass before and after shearing; wherein the point cloud data of the fracture surface includes point cloud data of the upper fracture surface and the lower fracture surface; The characterization module is used to match the point cloud data of the upper and lower fracture surfaces after shearing, discretize the matched point cloud data of the upper and lower fracture surfaces into grid units of the same size, obtain the initial opening of the corresponding grid unit through the elevation difference between the upper and lower fracture surfaces within the grid unit, and obtain the initial opening of the fractured rock mass to characterize the shear expansion effect; Using the elevation data of the upper and lower fracture surfaces of the fractured rock mass before and after shearing, the elevation changes of the upper and lower fracture surfaces before and after shearing are respectively obtained; based on the elevation changes of the upper and lower fracture surfaces before and after shearing, the equivalent opening of the fractured rock mass formed by roughness degradation is obtained; The simulation module is used to obtain the dynamic aperture of the fractured rock mass based on the initial aperture and equivalent aperture of the fractured rock mass; based on the dynamic aperture of the fractured rock mass, the point cloud data of the upper and lower fracture surfaces of the fractured rock mass after shearing are updated, and the geometric morphology data of the fractured rock mass after shearing are obtained using the updated point cloud data of the upper and lower fracture surfaces; The seepage and heat transfer coupling process of fractured rock mass is simulated by using the geometric morphology data of fractured rock mass after shearing and based on the seepage field and temperature field of fractured rock mass.

8. An electronic device, characterized in that: include: memory and processor; The memory is used to store computer programs; The processor is configured to implement the steps of a three-dimensional rough fractured rock mass seepage heat transfer numerical simulation method according to any one of claims 1 to 6 when executing the computer program stored in the memory.

9. A computer-readable storage medium, characterized in that Used to store a computer program, which, when executed by a processor, implements the steps of a three-dimensional rough fractured rock mass seepage heat transfer numerical simulation method as claimed in any one of claims 1 to 6.

Images