Method and system for calculating rough crack slip tendency considering stress shadow effect
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- CENT SOUTH UNIV
- Filing Date
- 2026-03-09
- Publication Date
- 2026-05-12
Smart Images
Figure CN121787201B_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of oil and gas reservoir development and is mainly applied to the accurate calculation and assessment of fracture network evolution and its induced secondary geological hazards during the stimulation of unconventional energy (such as deep geothermal, shale gas, and coalbed methane) reservoirs. This application can provide core calculation methods for parameter optimization of multi-cluster fracturing in shale gas horizontal wells, seismic risk assessment of geothermal dry hot rock reservoir stimulation, and fault activity monitoring. Background Technology
[0002] Under deep geological conditions, the geostress environment is extremely complex, and fracture development is significantly affected by tectonic movements. In multi-stage and multi-cluster fracturing processes, the mutual interference between fractures (i.e., the stress shadow effect) and the fault activation behavior of rough natural fractures have become key factors determining the success or failure of fracturing.
[0003] Rock fissures in nature are not ideally smooth planes. Due to differences in mineral composition and irregular fracture during shear failure, fissure surfaces often exhibit significant roughness and fractal characteristics. Current techniques show that roughness significantly alters fissure closure characteristics. When a fissure closes under pressure, surface micro-protrusions first come into contact, generating non-uniform contact stress. This microscopic contact mechanism leads to a highly non-uniform effective normal stress on the fissure surface, resulting in a non-linear spatial evolution of the stress field. The classic Barton-Bandis model laid the foundation for rough joint mechanics, but this model primarily focuses on normal closure. Current computational methods cannot quantitatively address how roughness induces a "dilatation effect" during shearing and how this shear activation affects the surrounding geostress field.
[0004] When fracturing is performed on oil and gas reservoirs, the opening of a single fracture under fluid pressure generates an additional compressive stress increment in the fracture normal direction. This induced stress increment leads to an increase in the initiation and propagation pressures of adjacent fractures, and may even force subsequent fractures to deflect, forming a "fishbone" or "dendritic" fracture network. Existing research mainly focuses on stress superposition in isotropic media, assuming that the intensity of stress shadowing increases with decreasing fracture spacing and increases with increasing fracture length. However, existing industrial software (such as simulators based on the plane fracture assumption) often underestimates the significant impact of stress shadowing, especially in considering residual stress interference during fracture closure, lacking effective computational methods.
[0005] For calculating stress shadowing and slip trends, current numerical methods mainly include the finite element method (FEM), discrete element method (DEM), and displacement discontinuity method, each with its own advantages and disadvantages. While the FEM can simulate heterogeneous strata and complex boundary conditions, the dramatic increase in mesh size makes computational time unacceptable in engineering practice for reservoirs containing hundreds or thousands of fractures. The DEM can naturally handle contact, slip, and failure problems, but its biggest disadvantage lies in the lack of a standard calibration process between macroscopic mechanical parameters and microscopic particle parameters, and it suffers from significant "boundary effects" and low computational efficiency in large-scale stress field simulations. The displacement discontinuity method simplifies fractures into discontinuous surfaces in elastic space, resulting in extremely high computational efficiency and making it the preferred choice for handling stress disturbances caused by multiple fractures with roughness.
[0006] During hydraulic fracturing, shear slip of natural fractures can induce earthquakes while effectively improving fracture conductivity. Furthermore, stress shadowing not only determines the fracture propagation geometry but also directly affects the shear instability and slip risk of the fracture surface. Therefore, accurately describing the stress field evolution and slip trend of rough fractures under high stress conditions is a pressing scientific challenge in petroleum engineering, geothermal engineering, and rock mechanics. Currently, commonly used methods include slip trend analysis and expansion trend analysis. However, these methods all assume a smooth fracture surface and a constant fracture friction coefficient, while neglecting the stress shadowing effect of propagating fractures. It is evident that the industry urgently needs a new method that can simultaneously satisfy computational efficiency and physical accuracy. This method should fully consider the stress shadowing effect and the shear slip mechanics of rough fractures, accurately predicting the slip initiation conditions and evolution trends under multi-fracture interference.
[0007] Current techniques for assessing the slip tendency of rough natural cracks, taking into account the stress shadowing effect generated during crack propagation, are still immature. The main difficulties lie in the following three aspects:
[0008] (1) Constitutive coupling mechanism of nonlinear closure and shear expansion of rough walls: Rough cracks exhibit significant nonlinear closure characteristics under compression, and shear dilatation occurs during shearing due to the "climbing" of the wall protrusions. This deformation process acts inversely on the surrounding rock mass, changing the local stress distribution. How to construct a unified constitutive equation that can simultaneously couple the effective aperture changes caused by normal nonlinear compression and shear is the core obstacle to achieving accurate stress shadow solutions.
[0009] (2) Dynamic full-field tensor superposition of stress shadow effect in multi-fracture system: In the process of multi-cluster fracturing, the induced stress fields generated by each fracture are highly overlapping and interfere with each other in three-dimensional space. When the fracture has a rough morphology, the stress shadow is no longer a uniform pressure barrier, but exhibits a very strong non-uniform distribution characteristic. How to calculate the "shielding" or "enhancing" effect of this non-uniform stress field on adjacent fractures, and quantitatively describe the real-time evolution of stress components on each fracture surface, puts extremely high demands on the algorithm.
[0010] (3) Nonlinear Iterative Convergence Control for Large-Scale Rough Crack Network Computation: After introducing the roughness constitutive model, the displacement discontinuity equations evolve into a highly nonlinear implicit equation system. When dealing with complex networks of hundreds or thousands of cracks, the frequent switching of contact states at the crack surfaces (cracking, contact, slippage) can easily lead to singular stiffness matrices or computational oscillations. How to design an efficient nonlinear iterative control strategy that ensures computational accuracy while also considering the convergence speed of large-scale computation is the core technical bottleneck for the practical application of this technology. Summary of the Invention
[0011] The purpose of this application is to overcome the shortcomings of existing technologies, such as the deviation in stress field calculation and inaccurate assessment of slip risk caused by simplifying hydraulic cracks as smooth planes. It provides a rough crack slip trend calculation method and system that can deeply couple the micro-morphological characteristics of the crack surface with the macro-mechanical response and consider the stress shadow effect of the rough crack morphology.
[0012] The technical solution provided in this application is as follows:
[0013] In a first aspect, this application provides a method for calculating the slip tendency of rough cracks considering the stress shadow effect, including:
[0014] S1. Three-dimensional reconstruction and feature parameter extraction of rough cracks: The original point cloud data of the rock fracture surface is obtained by non-contact three-dimensional laser scanning. Based on the original point cloud data, continuous analytical surfaces of the upper and lower walls of the crack are obtained. The crack aperture field is calculated through the continuous analytical surfaces of the upper and lower walls. The continuous analytical surfaces are discretized into displacement discontinuous elements. The crack aperture field parameters are assigned to each displacement discontinuous element to obtain the initial crack width distribution. Crack feature parameters are extracted based on the initial crack width distribution.
[0015] S2. Construction of the shear dilatation constitutive model for rough cracks and calculation of the opening: Combining the crack characteristic parameters obtained in step S1, a rough crack shear dilatation constitutive model integrating the normal nonlinear closure and shear dilatation effects is constructed to calculate the normal closure amount and shear expansion displacement of the crack; the mechanical opening of the crack is obtained by combining the normal closure amount and the shear expansion displacement and converted into the hydraulic opening.
[0016] S3. Construction of thermal-fluid-structure interaction crack model considering stress shadow and solution of induced stress: Combining the hydraulic opening obtained in step S2, a thermal-fluid-structure interaction crack model considering stress shadow effect is constructed. The thermal-fluid-structure interaction crack model adopts a two-way coupling method of finite element method and boundary element method, and an embedded discrete crack model is introduced to iteratively solve the multi-physics coupling process in porous elastic medium, quantify the stress interference effect between multiple cracks, and obtain the crack network induced stress considering stress shadow effect.
[0017] S4. Crack slip trend assessment: Based on the crack network induced stress and initial far-field stress tensor obtained in step S3, solve for the current geostress tensor of the crack element; combine the current geostress tensor to calculate the rough crack slip trend value considering the stress shadow effect.
[0018] In some possible implementations, step S1, obtaining the continuous analytical surfaces of the upper and lower walls of the crack based on the original point cloud data, includes: performing data preprocessing, quantitative characterization and analysis on the original point cloud data in sequence to obtain the continuous analytical surfaces of the upper and lower walls of the crack; the data preprocessing includes: using a median filtering algorithm to remove noise data from the original point cloud data, and using the least squares method to perform plane leveling on the denoised point cloud data so that the overall direction of the crack is parallel to the calculation coordinate axis plane; the quantitative characterization and analysis includes: introducing fractal dimension and roughness amplitude parameters to characterize the plane-leveled point cloud data, calculating the surface power spectral density and establishing the power-law relationship between spatial frequency and amplitude to form the continuous analytical surfaces of the upper and lower walls of the crack.
[0019] In some possible implementations, the calculation of the crack aperture field by the continuous analytical surfaces of the upper and lower walls in step S1 specifically involves: aligning the continuous analytical surfaces of the upper and lower walls, calculating the projection distance between them along the local normal, and dividing the projection distance into a grid according to a preset resolution to obtain the crack aperture field.
[0020] In some possible implementations, the discretization of the continuous analytical surfaces of the upper and lower walls into displacement discontinuous elements in step S1 specifically involves using GMSH to discretize the continuous analytical surfaces of the upper and lower walls of the crack into triangular or quadrilateral elements respectively; the initial crack width distribution is the rough crack initial width distribution data formed after each displacement discontinuous element corresponds one-to-one with the assigned aperture field parameters.
[0021] In some possible implementations, the rough crack shear dilatation constitutive model in step S2 includes a crack normal nonlinear closure model, used to calculate the normal closure amount of the crack under compressive conditions, as shown in the formula:
[0022] ;
[0023] in, It is a normal closed quantity. For effective normal stress, The initial positive stiffness of the crack, It is the maximum normal closure quantity;
[0024] and The calculation formula is:
[0025] ;
[0026] ;
[0027] in, For the crack wall strength, This refers to the surface roughness of the crack wall. This represents the initial crack aperture.
[0028] In some possible implementations, the rough crack shear dilatation constitutive model in step S2 includes a shear dilatation displacement model for calculating the shear dilatation displacement. The shear dilatation displacement model divides the crack deformation under shear displacement into a pre-peak stage and a post-peak stage. The pre-peak stage is the stage where the shear displacement is not greater than the peak shear displacement, and the post-peak stage is the stage where the shear displacement is greater than the peak shear displacement. The shear dilatation displacement is calculated for each pre-peak and post-peak stage using the following formula:
[0029] ;
[0030] in, This is the shear expansion displacement. This is the shear displacement. This represents the peak shear displacement. The peak shear expansion displacement is expressed as follows:
[0031] ;
[0032] ;
[0033] in, The characteristic length of the jointed rock block.
[0034] In some possible implementations, the finite element method and boundary element method are bidirectionally coupled in step S3, specifically as follows: the finite element method is used to solve for the pore pressure in the matrix, the fluid pressure in the crack, and the stress state of the matrix, which are then used as the mechanical boundary conditions of the boundary element method; the boundary element method is used to solve for the induced stress and induced displacement of the crack element based on the displacement discontinuity method, obtaining the crack width and crack stiffness, which are then fed back to the finite element method, iteratively correcting the fluid pressure field and stress field until the calculation results meet the preset convergence conditions; after each time step iteration converges, the crack propagation state is checked according to the crack propagation criterion; if the propagation conditions are met, the crack propagation distance and angle are calculated, and propagation crack elements are added at the crack tip; at the same time, the physical quantities calculated at the end of the current time step are used as new initial guess values, and the bidirectional coupling iterative calculation process of the finite element method and boundary element method is repeated until convergence; if the propagation conditions are not met, the calculation directly proceeds to the next time step.
[0035] In some possible implementations, in step S3, the formula for calculating the crack network-induced stress is:
[0036] ;
[0037] in, , and These represent all crack elements within the crack element. Shear induced stress along the strike direction, shear induced stress along the dip direction, and normal induced stress generated in the local coordinate system; Represents crack element Discontinuous displacement at the location in the crack element The influence coefficient of the tendency shear-induced stress generated at the location; Represents crack element Discontinuous displacement at the location in the crack element The influence coefficient of the tendency shear-induced stress generated at the location; Represents crack element Discontinuous displacement at the location in the crack element The influence coefficient of the normal induced stress generated at the location; Represents crack element Discontinuous displacement at the crack element The influence coefficient of shear-induced stress generated at the location; Represents crack element Discontinuous displacement at the crack element The influence coefficient of the tendency shear-induced stress generated at the location; Represents crack element Discontinuous displacement at the crack element The influence coefficient of the normal induced stress generated at the location; Represents crack element Discontinuous normal displacement at the crack element The influence coefficient of shear-induced stress generated at the location; Represents crack element Discontinuous normal displacement at the crack element The influence coefficient of the tendency shear-induced stress generated at the location; Represents crack element Discontinuous normal displacement at the crack element The influence coefficient of the normal induced stress generated at the location; and Different representation unit The discontinuous displacement along the dip, strike, and normal directions; This represents the total number of all crack elements.
[0038] In some possible implementations, the solution to obtain the current geostress tensor of the natural fracture element in step S4 is specifically: the initial far-field stress tensor is superimposed with the strike-induced shear stress, dip-induced shear stress, and normal-induced stress of the fracture element calculated in step S3 to obtain the current geostress tensor of the natural fracture element in the local coordinate system.
[0039] In some possible implementations, the slip tendency value in step S4 is the ratio of the absolute value of the shear stress of the crack element to the product of the absolute value of the effective normal stress and the crack friction coefficient.
[0040] Secondly, this application provides a rough crack slip tendency calculation system considering the stress shadow effect, for implementing the above-mentioned method, including:
[0041] A rough crack 3D reconstruction and feature parameter extraction module is used to acquire raw point cloud data of rock fracture surface using non-contact 3D laser scanning, and to obtain continuous analytical surfaces of the upper and lower walls of the crack based on the raw point cloud data; the crack aperture field is calculated from the continuous analytical surfaces of the upper and lower walls; the continuous analytical surfaces are discretized into displacement discontinuous elements, and the crack aperture field parameters are assigned to each displacement discontinuous element to obtain the initial crack width distribution; and crack feature parameters are extracted from the initial crack width distribution.
[0042] The module for constructing a constitutive model of rough crack shear dilatation and calculating its aperture is used to construct a constitutive model of rough crack shear dilatation that integrates normal nonlinear closure and shear dilatation effects by combining crack characteristic parameters, so as to calculate the normal closure amount and shear dilatation displacement of the crack; the mechanical aperture of the crack is obtained by combining the normal closure amount and shear dilatation displacement and converted into hydraulic aperture.
[0043] A module for constructing a thermal-fluid-structure interaction (TFI) crack model and solving for induced stress considering stress shadowing is used to construct a TFI crack model considering stress shadowing effect by combining hydraulic opening. The TFI crack model adopts a bidirectional coupling method of finite element method and boundary element method, and introduces an embedded discrete crack model to iteratively solve the multi-physics coupling process in porous elastic medium, quantify the stress interference effect between multiple cracks, and obtain the crack network induced stress considering stress shadowing effect.
[0044] The crack slip trend assessment module is used to solve the current geostress tensor of the crack element based on the crack network induced stress and the initial far-field stress tensor; and to calculate the rough crack slip trend value considering the stress shadow effect by combining the current geostress tensor.
[0045] Thirdly, this application provides an electronic device, including: a memory and a processor;
[0046] The memory is used to store computer programs;
[0047] The processor is used to invoke the computer program to execute the method described above.
[0048] Fourthly, this application provides a computer-readable storage medium storing a computer program that, when executed on an electronic device, causes the electronic device to perform the method described above.
[0049] Fifthly, this application provides a computer program product, including a computer program that, when run on an electronic device, causes the electronic device to perform the method described above.
[0050] The specific implementation methods of the second to fifth aspects of this application can refer to the implementation methods of the first aspect, and will not be elaborated here.
[0051] Beneficial effects:
[0052] This application constructs continuous analytical surfaces on the upper and lower walls of the crack and obtains the crack aperture field. The aperture field parameters are then assigned to the displacement discontinuous elements to form the initial crack width distribution. This can truly reflect the spatial geometry and aperture distribution characteristics of natural rough cracks, providing an accurate geometric basis for subsequent mechanical constitutive and multi-field coupling calculations, and improving the reliability of crack deformation and slip trend calculations.
[0053] Based on the initial crack width distribution, crack characteristic parameters such as crack wall roughness and wall strength are extracted. A crack normal nonlinear closure model and a shear expansion model are established respectively to achieve accurate quantification of normal closure amount and shear expansion displacement. Through staged conversion and linear interpolation of mechanical opening and hydraulic opening, the hydraulic opening is ensured to be continuous and smooth throughout the entire shear process, which is more in line with the actual crack seepage law.
[0054] A thermo-fluid-structure interaction crack model is constructed by bidirectional coupling of the finite element method and the boundary element method, combined with an embedded discrete crack model. After each time step iteration converges, the crack propagation state is checked. An extended crack element is added at the crack tip and the physical quantities of the current time step are iterated again until convergence. This can dynamically and continuously simulate the entire crack propagation process.
[0055] By calculating the induced stress field generated by the interaction between fractures in a multi-fracture system, the current true stress of the fracture is obtained by superimposing the induced stress with the invariant initial far-field stress tensor. This fully reflects and quantifies the stress shadowing effect between fractures. Then, the slip trend value is calculated based on the true stress, which significantly improves the accuracy and engineering applicability of slip trend assessment in multi-fracture systems. This can provide a more reliable theoretical basis for reservoir stimulation and engineering risk prevention and control. Attached Figure Description
[0056] Figure 1 This is a flowchart of a method in one embodiment of this application.
[0057] Figure 2 This is a diagram illustrating the distribution of rough crack (network) width in one embodiment of this application. Figure 2 (a) is a schematic diagram of the crack width distribution when the crack roughness coefficient is 0.1; Figure 2 (b) is a schematic diagram of the crack width distribution when the crack roughness coefficient is 0.5; Figure 2 (c) is a schematic diagram of the geometric morphology of a three-dimensional rough crack network; Figure 2 (d) is a schematic diagram of the crack width distribution of the three-dimensional rough crack network.
[0058] Figure 3 This is a schematic diagram of a constitutive model for rough crack shear dilatation in one embodiment of this application.
[0059] Figure 4 This is a schematic diagram comparing the results of this application with those of existing research in one embodiment, wherein... Figure 4 (a) is a schematic diagram showing the change of shear expansion displacement with shear displacement; Figure 4 (b) is a schematic diagram showing the changes in the mechanical aperture, hydraulic aperture and their ratio of the crack as a function of shear displacement.
[0060] Figure 5 This is a schematic diagram of the principle of a thermal-fluid-structure interaction crack propagation model considering stress shadowing in one embodiment of this application.
[0061] Figure 6 This is a schematic diagram showing the construction of the fracturing model for a coal-rock gas reservoir and the simulation results of the natural fracture slip trend in a case study of coal-rock gas reservoir fracturing.
[0062] Figure 7 This is a schematic diagram showing the post-fracturing matrix pore pressure and fracture fluid pressure distribution in a simulated case of hydraulic fracturing in a coal-rock gas reservoir; where... Figure 7 (a) is a top view; Figure 7 (b) is a three-dimensional view.
[0063] Figure 8 These are microseismic cloud images obtained from in-situ monitoring from different perspectives in a geothermal reservoir fracturing simulation case, where different colors represent different fracturing sections; Figure 8 (a) is a front view; Figure 8 (b) is a side view; Figure 8 (c) is a top view; Figure 8 (d) is a three-dimensional view.
[0064] Figure 9 This is a schematic diagram illustrating the relative spatial relationship between microseismic cloud maps and fracture networks in a geothermal reservoir fracturing simulation case; where... Figure 9 (a) is a geometric diagram of the three-dimensional crack network; Figure 9 (b) is an overlay diagram of the distribution of microseismic events and the crack network.
[0065] Figure 10 This diagram illustrates the comparison of natural fracture slip trends before and after hydraulic fracturing in a geothermal reservoir simulation case, as well as the influence of stress shadowing on slip trends. Figure 10 (a) is a diagram showing the fracture slip trend before fracturing; Figure 10 (b) is a diagram showing the slip trend of the first post-compression crack without considering stress shadowing; Figure 10 (c) is a diagram showing the slip trend of the second-stage post-compression crack without considering stress shadowing; Figure 10 (d) is a diagram showing the slip trend of the second-stage post-compression crack when considering stress shadowing. Detailed Implementation
[0066] To enable those skilled in the art to better understand the present application, the technical solution of the present application will be further described in detail below with reference to the embodiments and accompanying drawings.
[0067] This application employs a three-dimensional discontinuous displacement method and an embedded discrete fracture model to achieve high-precision capture of the spatial evolution of stress shadowing effects under multi-fracture interference conditions. It reveals the mechanical mechanism of the reaction between roughness-induced dilatation behavior and frictional strengthening mechanisms on stress field reconstruction, thereby enabling dynamic monitoring and quantitative evaluation of the slip trend of complex fracture networks throughout the entire construction process. This method provides an efficient quantitative calculation tool for multi-cluster fracturing parameter design in unconventional oil and gas development, reservoir stimulation volume optimization, and induced seismic risk prevention in deep geothermal development, significantly improving the reliability of fracture evolution prediction in complex geological environments.
[0068] The specific embodiments of this application will now be described with reference to the accompanying drawings.
[0069] Example 1:
[0070] like Figure 1 As shown in the embodiments of this application, a method for calculating the slip trend of rough cracks considering the stress shadow effect is disclosed, including:
[0071] S1. Three-dimensional reconstruction and feature parameter extraction of rough cracks: The original point cloud data of the rock fracture surface is obtained by non-contact three-dimensional laser scanning. Based on the original point cloud data, continuous analytical surfaces of the upper and lower walls of the crack are obtained. The crack aperture field is calculated through the continuous analytical surfaces of the upper and lower walls. The continuous analytical surfaces are discretized into displacement discontinuous elements. The crack aperture field parameters are assigned to each displacement discontinuous element to obtain the initial crack width distribution. Crack feature parameters are extracted based on the initial crack width distribution.
[0072] S2. Construction of the shear dilatation constitutive model for rough cracks and calculation of the opening: Combining the crack characteristic parameters obtained in step S1, a rough crack shear dilatation constitutive model integrating the normal nonlinear closure and shear dilatation effects is constructed to calculate the normal closure amount and shear expansion displacement of the crack; the mechanical opening of the crack is obtained by combining the normal closure amount and the shear expansion displacement and converted into the hydraulic opening.
[0073] S3. Construction of thermal-fluid-structure interaction crack model considering stress shadow and solution of induced stress: Combining the hydraulic opening obtained in step S2, a thermal-fluid-structure interaction crack model considering stress shadow effect is constructed. The thermal-fluid-structure interaction crack model adopts a two-way coupling method of finite element method and boundary element method, and an embedded discrete crack model is introduced to iteratively solve the multi-physics coupling process in porous elastic medium, quantify the stress interference effect between multiple cracks, and obtain the crack network induced stress considering stress shadow effect.
[0074] S4. Crack slip trend assessment: Based on the crack network induced stress and initial far-field stress tensor obtained in step S3, solve for the current geostress tensor of the crack element; combine the current geostress tensor to calculate the rough crack slip trend value considering the stress shadow effect.
[0075] This application employs a three-dimensional discontinuous displacement method and an embedded discrete fracture model to solve the multi-physics coupling processes involved in hydraulic fracturing, including fluid flow, heat transfer, rock mass deformation, natural fracture activation, and fracture propagation. This includes three-dimensional reconstruction of rough fractures and extraction of characteristic parameters; construction of a rough fracture shear dilatation constitutive model and calculation of its aperture; construction of a thermo-fluid-structure interaction fracture model considering stress shadowing and solution of induced stresses; and assessment of fracture slip trends. These parts are described in detail below.
[0076] (1) Three-dimensional reconstruction and feature parameter extraction of rough cracks:
[0077] To accurately reproduce the true rough morphology of natural cracks and provide a realistic geometric basis for subsequent mechanical calculations, firstly, a high-precision non-contact 3D laser scanner was used to perform a full-field scan of the rock fracture surface (upper and lower walls) to obtain raw point cloud data; the upper and lower walls of the rock cracks.
[0078] Then, a median filtering algorithm is used to remove noisy data from the original point cloud data, and the least squares method is used for plane leveling to ensure that the overall direction of the crack is parallel to the calculation coordinate axis plane.
[0079] We then introduce two parameters, fractal dimension (describing the irregularity of the crack surface) and roughness amplitude (describing the amplitude of protrusions / concavities on the crack surface), to characterize the denoised point cloud. By calculating the surface power spectral density, we establish the power-law relationship between spatial frequency and amplitude.
[0080] ;
[0081] in, For power spectral density, For spatial frequency, A constant characterizing the strength of a structure. For fractal dimensions.
[0082] This step allows discrete point cloud data to be transformed into crack surfaces with continuous analytical properties.
[0083] Based on the same scanning, processing, and analytical methods, continuous analytical surfaces of the upper and lower walls of the crack are constructed respectively. Then, the projected distances of the upper and lower walls along the local normal (the direction perpendicular to the crack wall) are found by alignment. The distances are then meshed according to an appropriate resolution to obtain the crack aperture field. The opening field It refers to different locations on the surface of the crack. The crack width varies at different locations due to the roughness of the crack. The aperture field is a quantitative and spatial description of this difference.
[0084] Finally, using 3D geometric modeling and mesh generation software (such as GMSH software), the reconstructed crack surface is discretized into a series of displacement discontinuous triangular or quadrilateral elements, and the aperture field is... By assigning these tasks to each computing unit, a system can be reconstructed, such as... Figure 2 The initial crack width distribution of different roughness cracks (network) is shown. The renderings provide core geometric parameters and mesh foundations for the subsequent mechanical constitutive modeling of rough cracks.
[0085] Fracture characteristic parameters are extracted from the initial fracture width distribution, including fracture wall roughness, fracture wall strength, initial fracture aperture, characteristic length of jointed rock blocks, and fracture element normal vector, which serve as inputs for subsequent constitutive modeling and stress calculation.
[0086] (2) Construction of constitutive model for shear dilatation of rough crack and calculation of aperture:
[0087] To address the nonlinear deformation characteristics of rough cracks under compression, a nonlinear normal closure model for the crack is established to calculate the normal closure amount of the crack under compression conditions. This model can be expressed as:
[0088] ;
[0089] in, This is the normal closure amount (unit: mm). Effective normal stress (unit: MPa). The initial normal stiffness of the crack (i.e., the normal stiffness of the crack when the effective normal stress approaches zero, which characterizes the initial compressive strength of the crack). It is the maximum normal closure amount (the ultimate compressive deformation of the crack when the effective normal stress approaches infinity).
[0090] According to rock mechanics experiments, and A relationship can be established with the crack wall strength and crack wall roughness:
[0091] ;
[0092] ;
[0093] in, For the crack wall strength, This refers to the surface roughness of the crack wall. The initial crack aperture, i.e., the original crack aperture when the normal stress is 0, is given by... The calculation yielded: .
[0094] Natural cracks are actually subjected to a combination of compression and shear, making their deformation behavior more complex than that under simple normal compression. As shear displacement increases, the protruding portions of the crack wall are gradually eroded, thus reducing the surface roughness of the crack wall. This reduces the crack friction coefficient, ultimately leading to a decrease in the crack's shear resistance. As the shear displacement continuously increases, the shear stress and shear dilatation exhibit different characteristics, which can be roughly divided into four stages (…). Figure 3 ): Linear elastic stage Hardening stage softening stage and residual stages The first two stages are collectively referred to as the pre-peak stage, and the latter two stages are collectively referred to as the post-peak stage. Based on the results of numerous direct shear tests on rocks, existing research has summarized the dynamic wall friction coefficient. The expression differs before and after the peak, establishing its dynamic relationship with shear displacement. From this, the formula for calculating the shear expansion displacement caused by shear displacement can be derived, as follows:
[0095] ;
[0096] in, This is the shear expansion displacement. This is the shear displacement. This represents the peak shear displacement. The peak shear expansion displacement is expressed as follows:
[0097] ;
[0098] ;
[0099] in, The characteristic length (in meters) of the jointed rock block.
[0100] Combined with the calculated normal closure quantity and shear expansion displacement Taking into account the combined effects of normal compression and shear expansion on crack aperture, the mechanical aperture of the rough crack is obtained:
[0101] ;
[0102] in, For mechanical opening.
[0103] Due to the roughness of the fracture surface, the fluid inside the fracture is not a smooth, flat plate flow, but rather a tortuous channel flow. Therefore, the mechanical aperture needs to be equivalent to the hydraulic aperture, which allows direct substitution into the cubic law to calculate the fracture permeability. The piecewise conversion formula between mechanical aperture and hydraulic aperture is as follows:
[0104] ;
[0105] in, This refers to the hydraulic opening.
[0106] for In such cases, linear interpolation is used for calculation. and Given two endpoints, calculate the hydraulic opening values corresponding to these two endpoints using the above piecewise conversion formula. and ,exist Within the transition range, the hydraulic opening value corresponding to the two endpoints is determined based on the actual shear displacement ratio. and The linear weighted calculation is performed between them, and the formula is as follows: ,in This is the actual shear displacement ratio. This ensures the hydraulic opening. The continuity and rationality of the calculations throughout the shearing process avoid abrupt changes at the boundaries of the piecewise transformation formula.
[0107] In the above formula, the mechanical opening degree and hydraulic opening All units are mm.
[0108] Figure 4 The changes in crack shear expansion displacement, mechanical aperture, and hydraulic aperture with shear displacement under compressive and shear stress conditions are presented. Under the same parameter conditions, the results are in complete agreement with existing research results, indicating that the rough crack shear expansion constitutive model proposed in this application is reliable.
[0109] This coarse crack shear dilatation constitutive model can realize the quantitative coupling calculation of crack normal and shear deformation, and the output hydraulic aperture provides key seepage parameters for subsequent multiphysics coupling calculations.
[0110] (3) Construction of a thermal-fluid-structure interaction crack model considering stress shadows and solution of induced stress:
[0111] To accurately capture the spatial evolution of stress shadowing effects in multi-fracture systems, this application proposes a novel thermo-fluid-structure interaction (TFI) simulation framework, combining core parameters such as hydraulic aperture and fracture stiffness output from the aforementioned coarse fracture dilatation constitutive model. This framework uses an iterative coupling method based on time / scale-dependent fracture stiffness to simulate full-three-dimensional fracture propagation in porous elastic media, enabling collaborative calculation of the multi-physics coupling process of fluid flow, heat transfer, rock deformation, and fracture propagation. Within this framework, triangular elements, independent of the matrix discretized from hexahedral meshes, are used to explicitly track and fit the morphology of the propagating fracture in each propagation step. Based on an embedded discrete fracture model, the finite element method (FEM) is used to solve the fluid dynamics system. The calculated fracture pressure and stress state of the embedded fracture master mesh provide accurate mechanical boundary conditions for the Boundary Element Method (BEM) solution. Conversely, the BEM module, by solving for the mechanical response of the fracture, provides the finite element method module with continuously changing fracture stiffness and aperture, which are crucial parameters for the former to close the mass balance equation related to the fracture elements, achieving bidirectional data feedback between the two. Finally, at the end of each time step, the total stress and crack tip displacement are calculated to estimate the velocity and direction of the newly formed crack in front of the crack tip.
[0112] Specifically, such as Figure 5 As shown, the finite element model and the boundary element model are coupled bidirectionally to solve the problem using the stress consistency condition at the crack boundary. First, the finite element method is used to calculate the pore pressure in the matrix, the fluid pressure in the crack, and the stress state of the matrix. The current stress state can then be calculated using the stress superposition principle. ,in For far-field stress, For temperature-induced stress, This represents the change in pore pressure. The stress induced by crack deformation is calculated. The calculated pressure and stress are directly used as boundary conditions for the discontinuous displacement method (indirect boundary element method). The crack is discretized into a set of triangular meshes. After applying the specified boundary stress to the surface of the crack element, the induced stress and induced displacement of the crack element can be solved in the local coordinate system of the crack element using the discontinuous displacement method, forming a two-way coupled solution scheme of finite element and boundary element methods.
[0113] The coupling between the finite element model and the boundary element crack propagation model is mainly achieved through iterative solutions between the fluid-structure interaction model and the crack propagation model, using the conservation of fluid mass within the crack and the interchange of stress boundaries. For example... Figure 5As shown, firstly, based on the results of the above-mentioned three-dimensional reconstruction and dilatation constitutive model of the rough crack, initial guesses are given for the crack width, initial crack stiffness, fluid pressure, and stress distribution. Then, the fluid pressure in the crack and the effective stress on the crack wall are solved using the finite element method. In each iteration step, the set of crack boundary stress conditions constructed from the physical site obtained by the finite element model is used as the stress boundary conditions of the boundary element model to further calculate the crack width (crack permeability), crack stiffness, and crack-induced stress. Subsequently, the crack parameters obtained by the boundary element model are fed back to the finite element model to correct the fluid pressure field and stress field. The above process is repeated iteratively until the changes in crack parameters, fluid pressure, and stress meet the pre-given convergence conditions, and then the calculation proceeds to the next time step.
[0114] Before proceeding to the next time step, the crack propagation state needs to be checked according to the crack propagation criterion: if the propagation conditions are met, the crack propagation distance and angle are given according to the aforementioned criteria, and a displacement discontinuity element (propagating crack element) is added at the crack tip. Simultaneously, the physical quantities calculated at the end of the current time step (such as crack pressure, stress distribution, crack width, stiffness, etc.) are used as initial guesses for new iterative calculations until convergence. If the propagation conditions are not met, the calculation proceeds directly to the next time step. This process proceeds step by step according to the two major calculation cycles of iterative convergence and crack propagation check until the calculation time reaches the given requirement, at which point the calculation ends. Ultimately, the full-field induced stress field and crack propagation state considering the stress shadow effect in the multi-crack system are obtained, providing core stress data support for crack slip trend assessment (natural crack slip trend assessment).
[0115] (4) Crack slip trend assessment:
[0116] Based on the bidirectional coupling solution results of the thermo-fluid-structure interaction crack model considering crack stress shadows, the crack network-induced stress can be predicted using the following formula when the crack fluid pressure and stress conditions are given:
[0117] ;
[0118] in, , and These represent all crack elements within the crack element. Shear induced stress along the strike direction, shear induced stress along the dip direction, and normal induced stress generated in the local coordinate system; Represents crack element Discontinuous displacement at the location in the crack element The influence coefficient of the tendency shear-induced stress generated at the location; Represents crack element Discontinuous displacement at the location in the crack element The influence coefficient of the tendency shear-induced stress generated at the location; Represents crack element Discontinuous displacement at the location in the crack element The influence coefficient of the normal induced stress generated at the location; Represents crack element Discontinuous displacement at the crack element The influence coefficient of shear-induced stress generated at the location; Represents crack element Discontinuous displacement at the crack element The influence coefficient of the tendency shear-induced stress generated at the location; Represents crack element Discontinuous displacement at the crack element The influence coefficient of the normal induced stress generated at the location; Represents crack element Discontinuous normal displacement at the crack element The influence coefficient of shear-induced stress generated at the location; Represents crack element Discontinuous normal displacement at the crack element The influence coefficient of the tendency shear-induced stress generated at the location; Represents crack element Discontinuous normal displacement at the crack element The influence coefficient of the normal induced stress generated at the location; and Different representation unit The discontinuous displacement along the dip, strike, and normal directions; This represents the total number of all crack elements.
[0119] Based on the predicted crack network-induced stress and combined with the initial far-field stress tensor measured in engineering, the actual stress conditions at any spatial location of a natural crack are obtained:
[0120] ;
[0121] in, For the current geostress tensor, For the initial far-field stress tensor, This is the normal vector of the crack element.
[0122] To quantitatively assess the shear slip risk of natural cracks considering the stress shadowing effect, a friction constitutive relation is introduced to define the shear slip tendency value of the crack. Its expression is:
[0123] ;
[0124] in, The coefficient of friction of the crack. and Let these represent the shear stress and normal stress acting on the fracture element, respectively, which can be derived from the current geostress tensor and the fracture element normal vector:
[0125] ;
[0126] ;
[0127] Therefore, the final calculation formula for the rough crack slip tendency value considering the stress shadow effect is obtained by integration:
[0128] .
[0129] This formula can be used to calculate the slip trend value of each fracture unit, enabling a quantitative and spatial assessment of the slip trend of rough fractures in multi-fracture systems. This provides a direct decision-making basis for seismic risk prevention and fracturing parameter optimization in the transformation of unconventional energy reservoirs such as shale oil and gas and geothermal development.
[0130] The effectiveness of this application will be verified through simulation experiments in two real-world scenarios: fracturing of coal-rock gas reservoirs and fracturing of geothermal reservoirs.
[0131] Coal-gas reservoir fracturing simulation:
[0132] In the development of coal-rock gas, the mechanical interaction between hydraulic fractures and natural fractures leads to the shear activation of local natural fractures. Similarly, the pressure increase in coal matrix pores and the stress shadowing effect of propagating fractures significantly alter the stability of natural fractures. This application's method is used to simulate the fracturing process of a coal-rock gas reservoir. The main parameters of the coal fracturing model are set as follows: the modeling area is 500 m × 500 m × 500 m, containing a total of 247 natural fractures with a length range of 25-50 m, an average fracture length of 40 m, a uniform fracture height of 65 m, a Young's modulus of coal of 8.0 GPa, a Poisson's ratio of 0.35, a vertical reservoir stress of 53 MPa, a maximum horizontal principal stress of 48 MPa, a minimum horizontal principal stress of 43 MPa, a coal matrix permeability of 0.1 mD, a fracturing fluid viscosity of 45 cP, an initial pore pressure of 25 MPa, and a pumping rate of 0.233 m³ / h. 3 / s, injection time of 90 min, top burial depth of coal reservoir of 2790 m, rock fracture toughness of 0.5 MPa .
[0133] Figure 6This paper presents a coal and rock fracturing grid model and a distribution map of the slip trend of natural fractures. The map demonstrates the modeling capability of the proposed method for complex coal and rock reservoirs containing 247 natural fractures, verifying the applicability of the scheme to large-scale fracture network calculations in engineering practice and solving the problem that existing methods struggle to handle complex fracture systems. The map also shows the numerical distribution of the slip trend of natural fractures, with specific slip trend values ranging from "0" to "1". When tensile or shear failure occurs, the slip trend value is "1". The value (0, 1) represents the distance from the fracture to failure, verifying that this application can achieve quantitative and visual assessment of slip trend. Furthermore, the slip trend judgment criteria (a value of 1 indicates failure) are consistent with engineering practice and have practical value. Simulation results show that during the fracturing process of coal and rock, the pressurization range is mainly located near the main hydraulic fracturing fracture and near natural fractures that are hydraulically connected to the main fracture. The main driving force for shear slip in the fracture comes from three parts: pore pressurization caused by the loss of fracturing fluid, the increase in fluid pressure within the fracture due to the connection between the hydraulic fracture and the natural fracture, and the stress shadowing effect generated after the fracture opens or shear slip occurs. This proves the necessity and rationality of incorporating the stress shadowing effect into the slip trend calculation in this application, demonstrating that the scheme can accurately capture the influence of stress shadowing on slip.
[0134] Figure 7 This paper presents the simulation results of the thermo-fluid-structure interaction model of this application on the pore pressure distribution in the post-fracturing matrix and the fluid pressure distribution in the fracture. It verifies that the scheme can accurately characterize the non-uniform distribution characteristics of fluid pressure in the matrix and fracture, and that fluid pressure is a core parameter for calculating induced stress and stress shadowing effects, thus proving the accuracy of the physical field calculations in this application. The spatial correlation characteristics between the pressure-increasing region and hydraulic and natural fractures are consistent with the actual fracturing pressure propagation laws in engineering, verifying that the embedded discrete fracture model + finite element-boundary element bidirectional coupling calculation method of this application can effectively realize the multi-physics coupling and transmission between the fracture and the matrix, solving the problems of underestimating stress shadowing and inaccurate pressure distribution calculations in existing industrial software.
[0135] Geothermal reservoir fracturing simulation:
[0136] The method described in this application was used to simulate the hydraulic fracturing process of a naturally fractured geothermal reservoir. The studied geothermal reservoir spanned 50 m × 50 m × 50 m and contained a total of 101 natural fractures. All natural fractures had an initial width of 0.13 mm, and their permeability could be estimated using the cubic law. The fracture and matrix meshes were completely independent, allowing the hydraulic fractures to expand flexibly within the three-dimensional matrix domain. The fracture network was divided into 13,029 triangular meshes, and the matrix was divided into 50,000 Cartesian structured meshes. The connection between the fractures and the matrix was established using an embedded discrete fracture method. The injection fracturing fluid temperature was 11℃, the reservoir temperature was 35℃, the Young's modulus was 71 GPa, the Poisson's ratio was 0.22, the vertical in-situ stress was 41.8 MPa, the maximum horizontal principal stress was 34 MPa, the minimum horizontal principal stress was 21.7 MPa, the matrix permeability was 0.1 mD, the fracturing fluid viscosity was automatically calculated based on temperature fracturing, the initial pore pressure was 8.3 MPa, the pumping rate was 400 mL / min, and the rock density was 2900 kg / m³. 3 The porosity is 0.01, and the rock fracture toughness is 2.0 MPa. The thermal conductivity is 5.0 W / (m·K), the friction coefficient of the natural crack is 0.6, the initial stiffness of the crack is 10 GPa, the heat capacity of the rock is 805 J / (kg·K), the Biot coefficient is 0.75, the compressibility modulus of the matrix particles is 650 GPa, and the maximum crack propagation step is 0.5 m.
[0137] Figure 8 The in-situ microseismic monitoring results for geothermal reservoir fracturing recorded 1743 microseismic events. These events were clustered and mapped onto the three-dimensional space of the test rig, with each microseismic event corresponding to the energy release process of the rock mass during fracturing. It is worth noting that all injection operations, including hydraulic fracturing, short-term flow testing, and long-term fluid circulation, can trigger microseismic events of varying intensities. Besides velocity-enhanced non-seismic slip, microseismic monitoring can capture all activities, including the opening and expansion of tensile fractures and the shearing and activation of natural fractures. Therefore, microseismic contour maps can effectively reflect the three-dimensional trajectory and expansion range of the hydraulic fracturing network, and to a certain extent, represent the reservoir stimulation volume range achieved by hydraulic fracturing.
[0138] Figure 9This paper presents a spatial overlay of the geometric distribution of hydraulic fractures and natural fractures after two stages of hydraulic fracturing, along with microseismic event cloud maps obtained from in-situ monitoring. This image spatially overlays the three-dimensional fracture network simulated by the method described in this application with the field microseismic cloud maps, verifying a high degree of match between the simulated fracture propagation trajectory and reservoir stimulation range and the distribution of microseismic events measured in the field. Sporadic microseismic events far from the hydraulic fractures also correspond to the simulated induced stress distribution. The image shows that microseismic events are mostly concentrated near the hydraulic fractures and the natural fractures connected to them. Some scattered microseismic events can also be detected in areas far from the hydraulic fractures. These microseismic events consist of two parts: one part is caused by the induced stress generated by the pressurized fractures, and the other part is generated by unknown (unresolved) natural fractures altering the stress boundary conditions during fracturing. Both major causes are covered by the stress shadow effect calculation and induced stress tensor superposition method of this application, proving that this application can accurately capture the core causes of crack slip / microseismic events in actual engineering, and solve the problem that existing methods ignore stress shadow and cannot explain complex microseismic distributions; at the same time, it verifies the accuracy of the rough crack shear dilatation constitutive model + thermo-fluid-structure interaction model of this application, and the spatial correspondence between the simulated crack activation region and the on-site microseismic region proves that the physical accuracy of the scheme meets the engineering requirements.
[0139] Figure 10 The evolution of shear slip trends in a natural fracture network is presented. Before fracturing, the natural fractures are generally in a stable state, with most fractures exhibiting slip trend values less than 0.75. After fracturing, due to factors such as matrix pressurization, fracture fluid pressurization, stress shadowing effects, and temperature-induced stress, shear failure occurs in natural fractures near hydraulic fractures and those directly communicating with them. Comparing the slip trend evolution before and after fracturing verifies that this application can dynamically monitor fracture slip trends throughout the entire fracturing process, solving the problem that existing methods cannot achieve dynamic assessment. Notably, comparing numerical simulation results with and without considering stress shadowing effects reveals that stress shadowing caused by fracture propagation is a crucial mechanism for triggering slip / microseismic events. Ignoring this mechanism leads to inaccurate predictions of the magnitude, scale, and geometry of microseismic events, and also underestimates the effectiveness of hydraulic fracturing in reservoir stimulation. This application incorporates stress shadowing effects into the calculation, enabling precise quantification of the impact of stress shadowing on fracture slip and achieving high-precision prediction of slip trends.
[0140] Example 2:
[0141] This application provides a rough crack slip trend calculation system considering the stress shadow effect, used to implement the method described in Embodiment 1, including:
[0142] A rough crack 3D reconstruction and feature parameter extraction module is used to acquire raw point cloud data of rock fracture surface using non-contact 3D laser scanning, and to obtain continuous analytical surfaces of the upper and lower walls of the crack based on the raw point cloud data; the crack aperture field is calculated from the continuous analytical surfaces of the upper and lower walls; the continuous analytical surfaces are discretized into displacement discontinuous elements, and the crack aperture field parameters are assigned to each displacement discontinuous element to obtain the initial crack width distribution; and crack feature parameters are extracted from the initial crack width distribution.
[0143] The module for constructing a constitutive model of rough crack shear dilatation and calculating its aperture is used to construct a constitutive model of rough crack shear dilatation that integrates normal nonlinear closure and shear dilatation effects by combining crack characteristic parameters, so as to calculate the normal closure amount and shear dilatation displacement of the crack; the mechanical aperture of the crack is obtained by combining the normal closure amount and shear dilatation displacement and converted into hydraulic aperture.
[0144] A module for constructing a thermal-fluid-structure interaction (TFI) crack model and solving for induced stress considering stress shadowing is used to construct a TFI crack model considering stress shadowing effect by combining hydraulic opening. The TFI crack model adopts a bidirectional coupling method of finite element method and boundary element method, and introduces an embedded discrete crack model to iteratively solve the multi-physics coupling process in porous elastic medium, quantify the stress interference effect between multiple cracks, and obtain the crack network induced stress considering stress shadowing effect.
[0145] The crack slip trend assessment module is used to solve the current geostress tensor of the crack element based on the crack network induced stress and the initial far-field stress tensor; and to calculate the rough crack slip trend value considering the stress shadow effect by combining the current geostress tensor.
[0146] Example 3:
[0147] This embodiment provides an electronic device, including: a memory and a processor;
[0148] The memory is used to store computer programs;
[0149] The processor is configured to invoke the computer program to execute the method as described in Embodiment 1.
[0150] Example 4:
[0151] This embodiment provides a computer-readable storage medium storing a computer program. When the computer program is run on an electronic device, it causes the electronic device to perform the method described in Embodiment 1.
[0152] Example 5:
[0153] This embodiment provides a computer program product, including a computer program that, when run on an electronic device, causes the electronic device to perform the method described in Embodiment 1.
[0154] The specific implementation of the system, electronic device, computer-readable storage medium, and computer program product provided in this application can be referred to the specific embodiments of the above methods, and will not be repeated here.
[0155] Obviously, those skilled in the art should understand that the various units or steps of this application described above can be implemented using general-purpose computing devices. They can be centralized on a single computing device or distributed across a network of multiple computing devices. Optionally, they can be implemented using computer-executable program code, thereby storing them in a storage device for execution by a computing device, or fabricating them separately as individual integrated circuit modules, or fabricating multiple modules or steps into a single integrated circuit module. Thus, this application is not limited to any particular combination of hardware and software.
[0156] The above description is merely a preferred embodiment of this application and is not intended to limit this application. Various modifications and variations can be made to this application by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the protection scope of this application.
Claims
1. A method for calculating the slip tendency of rough cracks considering the stress shadow effect, characterized in that, include: S1. Three-dimensional reconstruction and feature parameter extraction of rough cracks: The original point cloud data of the rock fracture surface is obtained by non-contact three-dimensional laser scanning. Based on the original point cloud data, continuous analytical surfaces of the upper and lower walls of the crack are obtained. The crack aperture field is calculated through the continuous analytical surfaces of the upper and lower walls. The continuous analytical surfaces are discretized into displacement discontinuous elements. The crack aperture field parameters are assigned to each displacement discontinuous element to obtain the initial crack width distribution. Crack feature parameters are extracted based on the initial crack width distribution. S2. Construction of the shear dilatation constitutive model for rough cracks and calculation of the opening: Combining the crack characteristic parameters obtained in step S1, a rough crack shear dilatation constitutive model integrating the normal nonlinear closure and shear dilatation effects is constructed to calculate the normal closure amount and shear expansion displacement of the crack; the mechanical opening of the crack is obtained by combining the normal closure amount and the shear expansion displacement and converted into the hydraulic opening. S3. Construction of thermal-fluid-structure interaction crack model considering stress shadow and solution of induced stress: Combining the hydraulic opening obtained in step S2, a thermal-fluid-structure interaction crack model considering stress shadow effect is constructed. The thermal-fluid-structure interaction crack model adopts a two-way coupling method of finite element method and boundary element method, and an embedded discrete crack model is introduced to iteratively solve the multi-physics coupling process in porous elastic medium, quantify the stress interference effect between multiple cracks, and obtain the crack network induced stress considering stress shadow effect. S4. Crack slip trend assessment: Based on the crack network induced stress and initial far-field stress tensor obtained in step S3, solve for the current geostress tensor of the crack element; combine the current geostress tensor to calculate the rough crack slip trend value considering the stress shadow effect.
2. The method according to claim 1, characterized in that, In step S1, obtaining continuous analytical surfaces of the upper and lower walls of the crack based on the original point cloud data includes: performing data preprocessing, quantitative characterization and analysis of the morphology on the original point cloud data in sequence to obtain continuous analytical surfaces of the upper and lower walls of the crack; the data preprocessing includes: using a median filtering algorithm to remove noise data from the original point cloud data, and using the least squares method to perform plane leveling on the denoised point cloud data so that the overall direction of the crack is parallel to the calculation coordinate axis plane; the quantitative characterization and analysis of the morphology includes: introducing fractal dimension and roughness amplitude parameters to characterize the plane-leveled point cloud data, calculating the surface power spectral density and establishing the power law relationship between spatial frequency and amplitude to form continuous analytical surfaces of the upper and lower walls of the crack.
3. The method according to claim 1, characterized in that, The step S1, which calculates the crack aperture field by means of the continuous analytical surfaces of the upper and lower walls, specifically involves aligning the continuous analytical surfaces of the upper and lower walls, calculating the projection distance between them along the local normal, and dividing the projection distance into a grid according to a preset resolution to obtain the crack aperture field.
4. The method according to claim 1, characterized in that, The rough crack shear dilatation constitutive model in step S2 includes a crack normal nonlinear closure model, used to calculate the normal closure amount of the crack under compressive conditions. The formula is: ； in, It is a normal closed quantity. For effective normal stress, The initial positive stiffness of the crack, It is the maximum normal closure quantity; and The calculation formula is: ； ； in, For the crack wall strength, This refers to the surface roughness of the crack wall. This represents the initial crack aperture.
5. The method according to claim 4, characterized in that, In step S2, the rough crack shear dilatation constitutive model includes a shear dilatation displacement model used to calculate the shear dilatation displacement. This shear dilatation displacement model divides the crack deformation under shear displacement into a pre-peak stage and a post-peak stage. The pre-peak stage is where the shear displacement is no greater than the peak shear displacement, and the post-peak stage is where the shear displacement is greater than the peak shear displacement. The shear dilatation displacement is calculated for each of the pre-peak and post-peak stages using the following formula: ； in, This is the shear expansion displacement. This is the shear displacement. This represents the peak shear displacement. The peak shear expansion displacement is expressed as follows: ； ； in, The characteristic length of the jointed rock block.
6. The method according to claim 1, characterized in that, The bidirectional coupling of the finite element method and the boundary element method in step S3 is as follows: The finite element method is used to solve for the pore pressure in the matrix, the fluid pressure in the crack, and the stress state of the matrix, which are then used as the mechanical boundary conditions of the boundary element method. The boundary element method is used to solve for the induced stress and induced displacement of the crack element based on the displacement discontinuity method, obtaining the crack width and crack stiffness, which are then fed back to the finite element method. The fluid pressure field and stress field are iteratively corrected until the calculation results meet the preset convergence conditions. After each time step of iteration convergence, the crack propagation state is checked according to the crack propagation criterion. If the propagation conditions are met, the crack propagation distance and angle are calculated, and a propagating crack element is added at the crack tip. At the same time, the physical quantities calculated at the end of the current time step are used as new initial guess values, and the bidirectional coupling iterative calculation process of the finite element method and the boundary element method is repeated until convergence. If the propagation conditions are not met, the calculation proceeds directly to the next time step.
7. The method according to claim 1, characterized in that, In step S3, the formula for calculating the induced stress in the crack network is: ； in, , and These represent all crack elements within the crack element. Shear induced stress along the strike direction, shear induced stress along the dip direction, and normal induced stress generated in the local coordinate system; Represents crack element Discontinuous displacement at the location in the crack element The influence coefficient of the tendency shear-induced stress generated at the location; Represents crack element Discontinuous displacement at the location in the crack element The influence coefficient of the tendency shear-induced stress generated at the location; Represents crack element Discontinuous displacement at the location in the crack element The influence coefficient of the normal induced stress generated at the location; Represents crack element Discontinuous displacement at the crack element The influence coefficient of shear-induced stress generated at the location; Represents crack element Discontinuous displacement at the crack element The influence coefficient of the tendency shear-induced stress generated at the location; Represents crack element Discontinuous displacement at the crack element The influence coefficient of the normal induced stress generated at the location; Represents crack element Discontinuous normal displacement at the crack element The influence coefficient of shear-induced stress generated at the location; Represents crack element Discontinuous normal displacement at the crack element The influence coefficient of the tendency shear-induced stress generated at the location; Represents crack element Discontinuous normal displacement at the crack element The influence coefficient of the normal induced stress generated at the location; and Different representation unit The discontinuous displacement along the dip, strike, and normal directions; This represents the total number of all crack elements.
8. The method according to claim 1, characterized in that, In step S4, the current geostress tensor of the fracture element is solved by superimposing the initial far-field stress tensor with the strike-induced shear stress, dip-induced shear stress, and normal-induced stress of the fracture element calculated in step S3 to obtain the current geostress tensor of the fracture element in the local coordinate system.
9. The method according to claim 1, characterized in that, The slip tendency value mentioned in step S4 is the ratio of the absolute value of the shear stress of the crack element to the product of the absolute value of the effective normal stress and the crack friction coefficient.
10. A rough crack slip tendency calculation system considering stress shadowing effect, characterized in that, To implement the method of any one of claims 1 to 9, comprising: A rough crack 3D reconstruction and feature parameter extraction module is used to acquire raw point cloud data of rock fracture surface using non-contact 3D laser scanning, and to obtain continuous analytical surfaces of the upper and lower walls of the crack based on the raw point cloud data; the crack aperture field is calculated from the continuous analytical surfaces of the upper and lower walls; the continuous analytical surfaces are discretized into displacement discontinuous elements, and the crack aperture field parameters are assigned to each displacement discontinuous element to obtain the initial crack width distribution; crack feature parameters are extracted from the initial crack width distribution. The module for constructing a constitutive model of rough crack shear dilatation and calculating its aperture is used to construct a constitutive model of rough crack shear dilatation that integrates normal nonlinear closure and shear dilatation effects by combining crack characteristic parameters, so as to calculate the normal closure amount and shear dilatation displacement of the crack; the mechanical aperture of the crack is obtained by combining the normal closure amount and shear dilatation displacement and converted into hydraulic aperture. A module for constructing a thermal-fluid-structure interaction (TFI) crack model and solving for induced stress considering stress shadowing is used to construct a TFI crack model considering stress shadowing effect by combining hydraulic opening. The TFI crack model adopts a bidirectional coupling method of finite element method and boundary element method, and introduces an embedded discrete crack model to iteratively solve the multi-physics coupling process in porous elastic medium, quantify the stress interference effect between multiple cracks, and obtain the crack network induced stress considering stress shadowing effect. The crack slip trend assessment module is used to solve the current geostress tensor of the crack element based on the crack network induced stress and the initial far-field stress tensor; and to calculate the rough crack slip trend value considering the stress shadow effect by combining the current geostress tensor.