A numerical simulation method for full three-dimensional joint-joint and joint-hole cross mechanical behavior

Abstract

The application discloses a numerical simulation method for full three-dimensional seam-seam and seam-hole cross mechanical behavior, comprising the following steps: obtaining the geometric information of natural fractures and solution cavities in a seam-hole type reservoir, and constructing a full three-dimensional geometric model of the fractures and solution cavities based on the geometric information; performing mesh division on the full three-dimensional geometric model, and establishing a full three-dimensional numerical calculation grid containing a matrix, a fracture network and a solution cavity; determining a control equation for simulating the mechanical behavior of the fractures and solution cavities; performing fracture evolution analysis based on the numerical calculation grid and the control equation, simulating the cracking, expansion, turning of hydraulic fractures and the cross mechanical behavior of the hydraulic fractures with natural fractures and solution cavities; and outputting the fracture shape, opening distribution and flow field distribution results in the fracture evolution process. The method considers the heat and mass transfer process in the matrix and fractures, matrix deformation, full three-dimensional fracture expansion and fracturing fluid filtration, and can provide a prediction and optimization tool for seam-hole type oil and gas reservoir reconstruction design.

Description

Technical Field

[0001] This invention relates to the field of oil and gas reservoir development technology, and in particular to a three-dimensional numerical simulation method for the mechanical behavior of fracture-fracture and fracture-cavity intersections. Background Technology

[0002] Carbonate oil and gas reservoirs are characterized by ancient strata and tight formations, generally facing the challenge of no natural production or low and rapidly declining production after drilling. The vast majority of oil and gas wells require acid fracturing ("acid fracturing") to achieve large-scale development of these reservoirs. How to effectively implement acid fracturing to fully connect fracture-cavity reservoirs has become a key to achieving efficient development of carbonate oil and gas reservoirs. Clarifying the natural fracture / cavity connection mechanism in multi-field, multi-scale pore-fracture-cavity composite media is of great significance for improving the volumetric acid fracturing effect of ultra-deep carbonate oil and gas reservoirs and accelerating the development of deep and ultra-deep oil and gas resources. The interaction mechanism between artificial fractures and natural fractures or other discontinuous mechanical weak surfaces in the stimulation of deep and ultra-deep tight oil and gas reservoirs has always been a research hotspot.

[0003] With the discovery of large and super-large fractured-vuggy carbonate reservoirs in ultra-deep formations, existing technologies for understanding the interaction mechanism between hydraulic fractures and caverns mainly focus on three aspects: (i) the stress distribution induced by different fracture-vuggy combinations; (ii) the control effect of cavern pressure, size, and distribution on fracture propagation paths; and (iii) the fracture-vuggy communication patterns after acid fracturing. However, existing technologies lack numerical simulation methods for predicting the full three-dimensional mechanical behavior of fracture-fracture and fracture-vuggy intersections in fractured-vuggy reservoirs. The main challenges currently lie in the following two aspects:

[0004] (1) The mechanical communication mode between hydraulic fractures, natural fractures and karst caves in three dimensions is unclear.

[0005] During the stimulation of fractured-vuggy oil and gas reservoirs, hydraulic fractures extend into the reservoir from the perforation point under the drive of fracturing fluid. Since these reservoirs typically contain natural fractures and caves, the hydraulic fractures are significantly deflected during their extension due to stress induced by these fractures and caves. The mechanical interactions between hydraulic fractures, natural fractures, and caves in three dimensions are highly complex, and the underlying mechanical communication patterns remain unclear.

[0006] (2) Numerical calculation models for the complex mechanical behavior of full three-dimensional seam-seam and seam-hole have not yet been formed.

[0007] The three-dimensional mechanical behavior of hydraulic fractures, natural fractures, and karst caves is limited by factors such as stress conditions, fluid pressure, injection rate, fluid viscosity, cave size, and fracture mechanical properties. Currently, the most common research method is laboratory fracturing experiments. However, fracturing experiments face significant challenges in areas such as similar experimental design, fracture-cavity sample preparation, and quantitative characterization of experimental results. Therefore, conducting related research using numerical simulation methods is particularly important. However, current models used to address fracture propagation in fracture-cavity oil and gas reservoirs are still two-dimensional or pseudo-three-dimensional, and a fully three-dimensional numerical calculation model for the complex mechanical behavior of fracture-fracture and fracture-cavity interactions has not yet been developed. There is an urgent need to develop corresponding numerical calculation methods and numerical calculation models capable of simulating the complex mechanical behavior of fracture-cavity intersections. Summary of the Invention

[0008] To address the shortcomings of unclear communication patterns between hydraulic fractures, natural fractures, and karst caves during fracturing of fractured-vuggy oil and gas reservoirs, and the immaturity of numerical simulation methods and technologies, this invention provides a full three-dimensional numerical simulation method for the mechanical behavior of fracture-fracture and fracture-vuggy intersections. The method is based on a numerical calculation model using the indirect boundary element method, which comprehensively considers the superimposed induced stress fields of fractures and karst caves. This model can accurately calculate the stress concentration effect near the cave wall, providing an effective predictive tool for exploring the interaction patterns between hydraulic fractures, natural fractures, and karst caves.

[0009] This invention provides a numerical simulation method for the mechanical behavior of cross-slit and cross-hole joints in full three dimensions, the method comprising:

[0010] S1: Obtain the geometric information of natural fractures and caverns in fracture-vuggy reservoirs, and construct a full three-dimensional geometric model of fractures and caverns based on the geometric information;

[0011] S2: Mesh the full three-dimensional geometric model to create a full three-dimensional numerical computation mesh that includes the matrix, crack network, and cavern body;

[0012] S3: Determine the governing equations for simulating the mechanical behavior of cracks and caverns. The governing equations shall include at least the seepage-stress-heat transfer coupling equations and the crack constitutive model.

[0013] S4: Based on a full three-dimensional numerical computation grid and governing equations, crack evolution analysis is performed to simulate the initiation, propagation, and direction of hydraulic cracks, as well as their cross-mechanical behavior with natural cracks and karst caves.

[0014] S5: Output the results of crack morphology, aperture distribution, and flow field distribution during the crack evolution process.

[0015] Furthermore, the geometric information is multi-source heterogeneous geological data, including points, lines, surfaces, and volumes; wherein, points are core observation and downhole imaging data; lines are logging and tracer monitoring data along the wellbore; surfaces are fiber optic monitoring and seismic profile data; and volumes are three-dimensional seismic inversion and comprehensive interpretation data.

[0016] Furthermore, the specific process of S1 is as follows:

[0017] S11: Acquire multi-source heterogeneous geological data of fractured-vuggy reservoirs, perform data fusion processing, and extract the boundary point set of natural fractures and the spatial location and morphological parameters of karst bodies;

[0018] S12: Based on the extracted set of natural crack boundary points, the complex-shaped natural cracks are simplified into polygonal crack surfaces, and the vertex coordinates and spatial orientation of each polygon are determined; based on the simplified polygonal crack surfaces, a three-dimensional surface model of the natural cracks is constructed using parametric surface fitting technology.

[0019] S13: Based on the extracted spatial location and morphological parameters of the karst cave, the complex karst cave is simplified into an ellipsoid, and the coordinates of the center of each ellipsoid, the length of the semi-axis in three directions, and the normal vector of the major axis symmetry plane are determined; based on the simplified ellipsoid parameters, a three-dimensional closed ellipsoid model of the karst cave is generated using closed surface construction technology.

[0020] S14: Integrate the three-dimensional curved surface models of all natural cracks and the three-dimensional closed ellipsoidal models of the cave bodies into the same spatial coordinate system to form a full three-dimensional geometric model that includes the crack network and the cave system.

[0021] Furthermore, the specific process of S2 is as follows:

[0022] S21: Treat the matrix rock mass of the fractured-vuggy reservoir as a continuous porous medium, and use a Cartesian structured grid to divide the matrix region to generate a background grid; the background grid covers the entire computational domain, providing a basic grid framework for subsequent embedding of fractures and vuggies, and records the number, vertex coordinates and volume information of each matrix grid cell.

[0023] S22: Based on the three-dimensional surface model of the natural crack constructed by S12, the boundary closed contour line of each crack is extracted; cubic spline curves or B-spline curves are used to interpolate and fit the crack boundary point set to reconstruct the accurate surface morphology of the crack; the built-in functions of the mesh generation tool (such as the addSurfaceFilling function of GMSH) are used to triangularly mesh each crack surface to generate a crack mesh composed of triangular elements.

[0024] Preferably, all triangular units belonging to the same crack are categorized and stored as a triangular mesh with a half-side data structure;

[0025] S23: Based on the three-dimensional closed ellipsoidal model of the cave constructed in S13, the geometric parameters of the cave are determined according to the given coordinates of the sphere's center, the lengths of the semi-axis in three directions, and the normal vector of the major axis symmetry plane. A symmetric meshing strategy is adopted. First, the ellipsoid is uniformly divided into 12 symmetrical arc segments using built-in functions of the mesh generation tool (such as GMSH's addEllipseArc function). Based on these symmetrical arc segments, the boundary lines of 8 sub-surfaces are constructed, and each sub-surface is triangulated to generate a closed triangular mesh of the cave. Then, the built-in functions of the mesh generation tool are used to triangulate each cave surface, generating a cave mesh composed of triangular elements. This ensures the uniformity and smoothness of the cave surface mesh, providing a mesh foundation for accurately calculating the stress concentration effect of the cave walls.

[0026] Preferably, all triangular units belonging to the same cave are categorized and stored as a triangular mesh with a half-side data structure.

[0027] S24: Treat the fracture triangular mesh and the cave triangular mesh as independent embedded objects; using a cyclic search algorithm, embed each fracture unit and cave unit into the matrix background mesh generated in S21 one by one, and establish spatial correspondences between fractures and matrix, fractures and fractures, fractures and caves, and caves and matrix; based on the spatial correspondences, establish a connectivity table for all mesh units, and calculate the conductivity coefficients between adjacent units, including matrix-fracture conductivity, fracture-fracture conductivity, and fracture-cavity conductivity;

[0028] S25: Integrate the matrix background mesh, crack mesh, and cave mesh to obtain a full three-dimensional numerical calculation mesh.

[0029] Preferably, the integrated mesh is subjected to quality checks and optimizations to ensure that the quality indicators of the mesh cells meet the requirements of numerical computation; a complete mesh data file containing all cell numbers, vertex coordinates, cell types and connectivity relationships is output to provide a computational mesh for solving the governing equations in S3.

[0030] Furthermore, the seepage-stress-heat transfer coupling equations include:

[0031] The governing equations of saturated rock mass mechanical equilibrium are:

[0032] ;

[0033] in, It is the divergence operator; For the fourth-order elastic tensor under drainage conditions; Biot coefficient; Pore ​​fluid pressure; For the Kronecker tensor; It is the linear thermal expansion coefficient; It is the bulk modulus of the drainage. The change in temperature , Indicates the current temperature. Indicates the reference temperature; For saturated rock mass density, , For fluid density, For the density of solid skeleton particles, It is porosity; It is the vector of gravitational acceleration; For strain tensor, , For fractional Laplace operators, For displacement;

[0034] mass conservation equation:

[0035] ;

[0036] in, For the volumetric strain of the rock, , , and They represent , and Displacement in three directions; Biot modulus; For time; This is the total thermal expansion coefficient of the saturated rock mass. , The coefficient of thermal expansion of the fluid; For fluid viscosity; This represents the volume factor of the fluid, with 0 indicating a reference state. , For permeability tensor; For fluid sources / sinks;

[0037] Energy conservation equation:

[0038] ;

[0039] in, The bulk modulus of saturated rock mass; For the total volumetric heat capacity, , The density of the rock matrix particles. For the heat capacity of the solid framework, For fluid heat capacity; This is the thermal conductivity tensor; Enthalpy; For thermal energy / sink.

[0040] Furthermore, the crack constitutive model includes a closed crack constitutive model and an open crack constitutive model, and the crack stiffness... Represented as:

[0041] ;

[0042] in, The effective stress is in the normal direction; The crack aperture of the crack element; For fluid pressure;

[0043] For natural fractures in a closed state during fracturing, their mechanical behavior is controlled by the effective normal stress and the surface roughness of the fracture. The constitutive model of a closed fracture is a mathematical expression used to describe the mechanical response behavior of a natural fracture under compressive closure, specifically characterizing the nonlinear relationship between the normal stress and the amount of fracture closure. Therefore, the constitutive model of a closed fracture is transformed into a description of the relationship between the normal stiffness and deformation of the fracture. The constitutive model of a closed fracture is as follows:

[0044] ;

[0045] ;

[0046] ;

[0047] in, The initial normal stiffness; The maximum allowable degree of closure; This is the joint roughness coefficient; Joint compressive strength;

[0048] The constitutive model of an opening fracture is a mathematical expression used to describe the mechanical response behavior of hydraulic fractures in an opening state under fluid pressure. Specifically, it characterizes the dynamic relationship between the fluid pressure within the fracture and the fracture aperture, as well as the fracture's resistance to deformation (i.e., fracture stiffness). For a fracture in an opening state driven by high fluid pressure, its stiffness depends on the fluid pressure within the fracture, the matrix mechanical properties (Young's modulus and Poisson's ratio), and the fracture geometry (size and surface morphology). The constitutive model of an opening fracture is as follows:

[0049] ;

[0050] ;

[0051] ;

[0052] ;

[0053] in, A vector consisting of the stiffnesses of N crack elements; The disturbance stress of the crack element; The deformation of the crack under disturbed stress conditions; This is a vector composed of the crack apertures of N crack elements after the addition of perturbation stress; It is a vector composed of the crack apertures of N crack elements under undisturbed stress. The discontinuous normal displacement vector of N crack elements after adding disturbance stress; For N crack elements, the normal displacement is discontinuous when there is no disturbance stress. The number of crack elements; This is the influence coefficient matrix; This is a vector representation of the stress borne by N crack elements; For disturbance stress, .

[0054] Furthermore, the specific process of S4 is as follows:

[0055] S401: Determining the initiation condition of hydraulic fractures: Based on the seepage-stress-heat transfer coupling equation, calculate the stress state at each potential initiation point in the current time step; use the maximum principal stress criterion to determine the tension initiation condition: when the maximum effective stress on the rock mass... Greater than or equal to the tensile strength of the rock At that time, the crack initiation occurred; simultaneously, considering the thermal stress effect caused by the injection of cryogenic fracturing fluid, when the deviatoric stress level... When the value is greater than or equal to 1.0, the natural crack is determined to have initiated shear slip cracking.

[0056] ;

[0057] in, This represents the deviatoric stress under the current condition; This represents the ultimate deviatoric stress at which the rock fails. This represents the current maximum effective principal stress; This represents the current minimum effective principal stress; This is the ultimate deviatoric stress at failure.

[0058] S402. Calculation of stress intensity factor at the crack tip: Based on displacement discontinuity, crack geometry, and rock properties, calculate the stress intensity factor at the crack tip using the displacement discontinuity method. :

[0059] ;

[0060] in, These are empirical constants; For parameters related to rock properties, , , It is Young's modulus. It is Poisson's ratio; For discontinuities in displacement; subscript These represent the pure opening mode (i.e.) Type), sliding mode (i.e.) (type), scissor mode (i.e.) Type); Subscript , , These represent the directions of normal opening, strike-slip shear, and dip-slip shear, respectively. It is the distance from the centroid of the triangular unit to the crack propagation front;

[0061] S403. Determine the crack propagation direction and turning angle:

[0062] The maximum circumferential stress criterion is used to handle pure opening mode extension and mixed mode extension, when the equivalent stress intensity factor Achieving rock toughness ( When the stress is constant, the crack will propagate; according to the Richard criterion, the crack orientation angle is calculated based on three stress intensity factors. The hybrid modes include I+II, I+III, and I+II+III.

[0063] Equivalent stress intensity factor for:

[0064] ;

[0065] ;

[0066] in, The crack deflection angle under triaxial loading; As a substitute variable;

[0067] Crack turning angle The calculation formula is:

[0068] ;

[0069] Among them, when When the crack turning angle is positive, and When the crack turning angle is negative;

[0070] S404. Calculate the crack propagation rate:

[0071] The crack propagation rate was calculated using the Paris scaling law model, and the propagation rate at the crack tip was correlated with the stress intensity factor amplitude to determine the crack propagation length per unit time step.

[0072] S405. Simulated stress interference between cracks and natural cracks / cavities:

[0073] The induced stress field is obtained by calculating the induced stress on the target mesh element by all crack and cavern elements using the indirect boundary element method, where the target mesh element is either a crack element or a cavern element. The induced stress is expressed as:

[0074] ;

[0075] in, , and These represent all crack and cavern elements within the target mesh element. In the local coordinate system, the normal induced stress, the shear induced stress along the strike direction, and the shear induced stress along the dip direction; The influence coefficient is as follows: Indicates a crack or cave unit Discontinuous displacement of the crack normal or virtual stress of the cavity normal in the target mesh element The influence coefficient of the normal induced stress generated at the location; Indicates a crack or cave unit Discontinuous displacement of the crack normal or virtual stress of the cavity normal in the target mesh element The influence coefficient of shear-induced stress generated at the location; Indicates a crack or cave unit Discontinuous displacement of the crack normal or virtual stress of the cavity normal in the target mesh element The influence coefficient of the tendency shear-induced stress generated at the location; Indicates a crack or cave unit The direction of cracks, discontinuous displacement, or cavern direction should be virtually reflected in the target mesh element. The influence coefficient of the normal induced stress generated at the location; Indicates a crack or cave unit The direction of cracks, discontinuous displacement, or cavern direction at the location are virtual stresses in the target mesh element. The influence coefficient of shear-induced stress generated at the location; Indicates a crack or cave unit The direction of cracks, discontinuous displacement, or cavern direction at the location are virtual stresses in the target mesh element. The influence coefficient of the tendency shear-induced stress generated at the location; Indicates a crack or cave unit Cracks at the location tend to exhibit discontinuous displacement or cavitation tend to cause virtual stress in the target mesh element. The influence coefficient of the normal induced stress generated at the location; Indicates a crack or cave unit Cracks at the location tend to exhibit discontinuous displacement or cavitation tend to cause virtual stress in the target mesh element. The influence coefficient of shear-induced stress generated at the location; Indicates a crack or cave unit The cracks at the location tend to exhibit discontinuous displacement or cavern tendency, resulting in virtual stress within the crack / cavity unit. The influence coefficient of the tendency shear-induced stress generated at the location; for crack elements, , and Representing units respectively The displacement is discontinuous along the normal, strike, and dip directions; for a cave unit, , and Representing units respectively Virtual stress along the normal, strike, and dip directions; This represents the total number of all fracture and cavern cell grids;

[0076] For the target mesh element, force balance requires that the induced stress equals the net stress, which is the internal fluid pressure minus the in-situ stress, i.e.:

[0077] ;

[0078] in, Target mesh cell Fluid pressure in the middle; It is the pore pressure in the matrix mesh into which the target mesh element is embedded. Background stress of the matrix mesh;

[0079] S406. Based on the induced stress field calculated in S405, determine the cross-behavior pattern when hydraulic fractures approach natural fractures or karst caves; wherein, the cross-behavior pattern includes penetration mode, fracture arrest mode, shear slip mode, fusion mode, and attraction or repulsion mode.

[0080] Furthermore, the judgment process is as follows: when the hydraulic fracture tip meets the restart condition on the other side of the natural fracture and the natural fracture does not undergo shear slip, it is judged as a penetration mode; when the natural fracture undergoes shear slip, the fracturing fluid flows into the natural fracture, causing pressure dissipation and the hydraulic fracture tip becomes blunt, it is judged as a termination mode; when the shear stress on the natural fracture exceeds its shear strength, shear slip occurs, activating the permeability of the natural fracture, it is judged as a shear slip mode; when two coplanar fractures approach each other, it is judged as fusion, and the fracture tips merge to form a new extension front, it is judged as a fusion mode; based on the stress concentration characteristics around the cave, the interaction type between the hydraulic fracture and the cave is judged, including attraction followed by gradual approach, attraction followed by repulsion, and repulsion followed by moving away from the cave body, it is judged as an attraction or repulsion mode.

[0081] Furthermore, S4 also includes using an adaptive mesh method to track the full three-dimensional crack propagation front and its geometry, the specific process of which is as follows:

[0082] S407: Based on the crack propagation direction determined by S403 and the crack propagation velocity calculated by S404, generate the propagation vector set at the crack tip at the current time step; connect the endpoints of the propagation vectors with the endpoints of the crack leading edge in the previous propagation step to obtain the new propagation leading edge contour line of the crack; record the spatial coordinates of each node of the new leading edge and the number of the crack element to which it belongs.

[0083] S408. Based on the geometric characteristics of the new leading edge, the quality of the extended leading edge mesh is adjusted in three cases:

[0084] Case a: When the included angle between two adjacent sides is greater than the preset angle, connect the vertex of the included angle and the midpoint of the opposite side to split the triangular unit into two calculation units, so as to avoid the unit being too long and narrow and affecting the calculation efficiency.

[0085] Case b: When the length of the expansion vector is less than the first preset threshold of the maximum expansion step size, the endpoints of each endpoint of the previous step are moved directly to the corresponding new endpoints according to the expansion vector, and the original endpoints are deleted to avoid the excessively dense grid causing a rapid increase in degrees of freedom and a significant decrease in computational efficiency.

[0086] Case c: When the expansion vector exceeds the second preset threshold of the maximum expansion step size, find the new leading edge endpoint according to the expansion vector, and add a new endpoint at the midpoint of the new expansion vector. Treat the single-step expansion as two-step uniform expansion to avoid the crack being too narrow and affecting the computational efficiency.

[0087] S409. Smoothing the leading edge of the crack: The Savitzky-Golay filter fitting method is used to smooth the latest leading edge node; the sawtooth fluctuation of the leading edge node is eliminated by local polynomial fitting to maintain the smooth and continuous characteristics of the crack leading edge; distorted mesh elements are avoided to ensure the computational stability of subsequent expansion steps.

[0088] S410: Handling penetration, crack arrest, and fusion of intersecting cracks:

[0089] Based on the cross-behavior patterns determined by S406, mesh processing is performed on the interactions of different types of cracks:

[0090] Penetration mode: When a hydraulic fracture restarts on the other side of a natural fracture, no additional processing of the mesh is required, allowing the fracture to continue to propagate.

[0091] Crack arrest mode: When a hydraulic crack stops at a natural crack, first determine the crack intersection line, and then perform grid truncation treatment with the intersection line as a constraint to stop the further propagation of the hydraulic crack.

[0092] Fusion mode: When two coplanar cracks approach and merge, the overlapping mesh is deleted to merge the leading edges of the two cracks and form a new continuous crack surface.

[0093] Shear slip mode: When a natural crack undergoes shear slip, new crack elements are generated at the slip surface, updating the topology of the crack network.

[0094] S411: Update crack network topology and computational mesh: Integrate the mesh cells adjusted in S408, the leading edge nodes smoothed in S409, and the cross crack mesh processed in S410; generate the complete crack network triangular mesh after the current expansion step; update the numbering, node coordinates, and connectivity of all crack cells to provide an updated computational mesh for the S401-S406 analysis in the next expansion step.

[0095] This invention proposes a numerical simulation method for the cross-mechanical behavior of fracture-fracture and fracture-vuggy structures in three dimensions. The model constructed by this method is a closed / open three-dimensional constitutive model of fractures, which can be used to calculate fracture stiffness under different stress states before, during, and after compression. A "fracture network + cavern" mesh generation technique is established for three-dimensional mesh modeling of fracture-vuggy oil and gas reservoirs. An adaptive mesh generation algorithm is proposed to track the three-dimensional fracture propagation front and its geometry. A numerical calculation model for the cross-mechanical behavior of fracture-fracture and fracture-vuggy structures in three dimensions is formed, providing prediction and optimization tools for the design of fracture-vuggy oil and gas reservoir stimulation. This method is not only applicable to the optimization design of fracture-vuggy oil and gas reservoir fracturing, but also to wellbore stability prediction, airtightness and stability analysis of compressed air energy storage artificial chambers, stability analysis of CO2 underground buried storage caprock, and stability and integrity assessment of salt cavern oil and gas reservoirs. It has broad application prospects in oil and gas resource development and carbon, waste, and energy storage fields. Attached Figure Description

[0096] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0097] Figure 1 This is a flowchart illustrating a numerical simulation method for the mechanical behavior of a fully three-dimensional seam-seam and seam-hole intersection provided in an embodiment of the present invention;

[0098] Figure 2 This is a schematic diagram of the full three-dimensional numerical computation grid provided in an embodiment of the present invention;

[0099] Figure 3 This is a schematic diagram of the constitutive models of closed and open three-dimensional curved surface cracks provided in the embodiments of the present invention;

[0100] Figure 4 This is a schematic diagram of the leading edge tracking of the extended crack provided in an embodiment of the present invention: Figure 4 (a) is a schematic diagram of the crack propagation vector and the new crack leading edge obtained from the propagation vector; Figure 4 (b) is a schematic diagram of cell splitting, node merging and cell addition operations for the crack propagation front mesh;

[0101] Figure 5 This is a schematic diagram of mesh processing when cracks intersect: Figure 5 (a) is a schematic diagram of the fusion of coplanar cracks when they meet; Figure 5 (b) is a schematic diagram of mesh generation when non-coplanar cracks meet and pass through or capture each other;

[0102] Figure 6 This is a schematic diagram of the three-dimensional morphology and crack aperture distribution of crack groups under different preset crack network morphologies provided in the embodiments of the present invention: Figure 6 (a) is a historical curve of the total area change of the crack group; Figure 6 (b) is a schematic diagram of five typical interaction modes of three-dimensional non-planar cracks;

[0103] Figure 7 This is a schematic diagram illustrating the competitive and collaborative evolution characteristics of three-dimensional non-planar crack groups under different initial crack geometries, as provided in this embodiment of the invention. Figure 7 (a) is a diagram showing the competitive propagation results of the initial crack at a dip angle of 0°; Figure 7 (b) is a diagram showing the competitive propagation results of the initial crack at a dip angle of 30°; Figure 7 (c) is a diagram showing the competitive propagation results of the initial crack at a 45° dip angle; Figure 7 (d) is a diagram showing the competitive propagation results of the initial crack at a dip angle of 60°; Figure 7(e) is a diagram showing the competitive propagation results of the initial crack at a 90° dip angle; Figure 7 (f) is a diagram showing the competitive propagation results of randomly distributed initial cracks with varying dip angles;

[0104] Figure 8 This is a schematic diagram of the karst cave-fracture-wellbore grid partitioning results provided in an embodiment of the present invention;

[0105] Figure 9 This is a diagram illustrating the simulation results of the mechanical behavior of crack-crack and crack-cavity intersections provided in this embodiment of the invention: Figure 9 (a) is a schematic diagram of the fracture propagation results after 5 steps during hydraulic fracturing; Figure 9 (b) is a schematic diagram of the fracture propagation results after 10 steps during hydraulic fracturing; Figure 9 (c) is a schematic diagram of the fracture propagation results after 15 steps during hydraulic fracturing; Figure 9 (d) is a schematic diagram of the crack propagation results after 20 steps during hydraulic fracturing. Detailed Implementation

[0106] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be described in detail below. Obviously, the described embodiments are merely some embodiments of this invention, and not all embodiments. Based on the embodiments of this invention, all other implementation methods obtained by those skilled in the art without creative effort are within the scope of protection of this invention.

[0107] Example 1

[0108] like Figure 1 As shown in the figure, this embodiment provides a numerical simulation method for the mechanical behavior of cross-slits and cross-holes in full three dimensions, the method comprising:

[0109] S1: Obtain the geometric information of natural fractures and caves in fracture-vuggy reservoirs, and construct a full three-dimensional geometric model of fractures and caves based on the geometric information.

[0110] Specifically, the geometric information is multi-source heterogeneous geological data, including points, lines, surfaces, and volumes; where points are core observation and downhole imaging data; lines are logging and tracer monitoring data along the wellbore; surfaces are fiber optic monitoring and seismic profile data; and volumes are three-dimensional seismic inversion and comprehensive interpretation data.

[0111] More specifically, the specific process of S1 is as follows:

[0112] S11: Acquire multi-source heterogeneous geological data of fractured-vuggy reservoirs, perform data fusion processing, and extract the boundary point set of natural fractures and the spatial location and morphological parameters of karst bodies;

[0113] S12: Based on the extracted set of natural crack boundary points, the complex-shaped natural cracks are simplified into polygonal crack surfaces, and the vertex coordinates and spatial orientation of each polygon are determined; based on the simplified polygonal crack surfaces, a three-dimensional surface model of the natural cracks is constructed using parametric surface fitting technology.

[0114] S13: Based on the extracted spatial location and morphological parameters of the karst cave, the complex karst cave is simplified into an ellipsoid, and the coordinates of the center of each ellipsoid, the length of the semi-axis in three directions, and the normal vector of the major axis symmetry plane are determined; based on the simplified ellipsoid parameters, a three-dimensional closed ellipsoid model of the karst cave is generated using closed surface construction technology.

[0115] S14: Integrate the three-dimensional curved surface models of all natural cracks and the three-dimensional closed ellipsoidal models of the cave bodies into the same spatial coordinate system to form a full three-dimensional geometric model that includes the crack network and the cave system.

[0116] S2: Mesh the full 3D geometric model to create a full 3D numerical computation mesh that includes the matrix, crack network, and cavern body.

[0117] S21: The reservoir rock mass is considered as a three-dimensional medium, namely, the matrix, fractures, and vugs. The matrix rock mass of the fractured-vug reservoir is treated as a continuous porous medium. A Cartesian structured grid is used to partition the matrix region, generating a background grid. This background grid covers the entire computational domain, providing a basic grid framework for subsequent embedding of fractures and vugs, and records the number, vertex coordinates, and volume information of each matrix grid cell. Specifically, the vertex coordinates of the triangulated cells (…) The corresponding number You can use the GMSH built-in functions getNodes and getNodesByElementType to obtain them.

[0118] S22: Based on the 3D surface model of the natural crack constructed in S12, the closed contour line of the boundary of each crack is extracted; cubic spline curves or B-spline curves are used to interpolate and fit the crack boundary point set to reconstruct the accurate surface morphology of the crack; the built-in functions of the mesh generation tool (such as the addSurfaceFilling function of GMSH) are used to triangularly mesh each crack surface to generate a crack mesh composed of triangular elements; all triangular elements belonging to the same crack are classified and stored as triangular meshes with half-side data structures through OpenMesh. In specific implementation, for the crack part, based on a series of crack boundary point sets, such as crack F#1{B1, B2, ..., Bn} as parametric surface constraints, cubic spline curves or B-spline curves are used to interpolate and fit discrete points to obtain closed contour lines. Based on these contour lines, the GMSH built-in function addSurfaceFilling is used to triangularly mesh each crack, and the same processing method is used for each crack to obtain the triangular mesh of the crack network.

[0119] S23: Based on the three-dimensional closed ellipsoidal model of the cave constructed in S13, the geometric parameters of the cave are determined according to the given coordinates of the sphere's center, the lengths of the semi-axis in three directions, and the normal vector of the major axis symmetry plane. A symmetric meshing strategy is adopted. First, the ellipsoid is uniformly divided into multiple symmetrical arc segments using built-in functions of mesh generation tools (such as GMSH's addEllipseArc function). Based on these symmetrical arc segments, the boundary lines of multiple sub-surfaces are constructed, and each sub-surface is triangulated to generate a closed triangular mesh of the cave. All triangular elements belonging to the same cave are categorized and stored as triangular meshes with half-edge data structures using OpenMesh. This ensures the uniformity and smoothness of the mesh on the cave surface, providing a mesh foundation for accurately calculating the stress concentration effect of the cave walls. In practice, for the cave section, based on the given ellipsoid center, minor axis, major axis, and major axis symmetry plane normal vector, the ellipsoid is first divided into 12 symmetrical arc segments using the GMSH built-in function addEllipseArc to ensure the uniformity of the mesh. Then, the function addCurveLoop is used to construct the boundary lines of 8 sub-surfaces. Next, the addSurfaceFilling function is used to generate 8 sub-surfaces. Finally, each surface is triangulated to obtain a series of triangulated meshes of the cave body.

[0120] S24: Treat the fracture triangular mesh and the cave triangular mesh as independent embedded objects; using a cyclic search algorithm, embed each fracture unit and cave unit into the matrix background mesh generated in S21 one by one, and establish spatial correspondences between fractures and matrix, fractures and fractures, fractures and caves, and caves and matrix; based on the spatial correspondences, establish a connectivity table for all mesh units, and calculate the conductivity coefficients between adjacent units, including matrix-fracture conductivity, fracture-fracture conductivity, and fracture-cavity conductivity;

[0121] S25: Integrating the matrix background mesh, crack mesh, and cavern mesh to obtain, as shown... Figure 2 The full three-dimensional numerical computation grid is shown.

[0122] Preferably, the integrated mesh is subjected to quality checks and optimizations to ensure that the quality indicators of the mesh cells meet the requirements of numerical computation; a complete mesh data file containing all cell numbers, vertex coordinates, cell types and connectivity relationships is output to provide a computational mesh for solving the governing equations in S3.

[0123] S3: Determine the governing equations for simulating the mechanical behavior of cracks and caverns. The governing equations shall include at least the seepage-stress-heat transfer coupling equations and the crack constitutive model.

[0124] In practical implementation, under the assumption of small deformation, the energy and mass conservation equations and stress balance equations can be obtained through the thermoelastic porous media theory. Applying Newton's second law to a unit mass of saturated porous media (ignoring the inertial forces of the rock mass), the mechanical governing equations of the system can be expressed as:

[0125] ;

[0126] in, It is the divergence operator; Let Cauchy's total stress tensor be the total stress tensor. It is the vector of gravitational acceleration; It is the density of saturated rock mass. For fluid density, For the density of solid skeleton particles, This refers to porosity. Note that the stress sign is defined as follows: tensile stress is positive, and compressive stress is negative. The effective stress principle indicates that the total external force on a rock mass is borne jointly by the rock skeleton and pore fluids, and the Biot coefficient... The weights of the fluid distribution on the total stress are given, mathematically expressed as:

[0127] ;

[0128] Furthermore, rocks shrink or expand under different temperature conditions, and the resulting change in rock volume (volume strain) can be equivalent to a stress load, commonly referred to as thermal stress. Therefore, the complete expression of the constitutive equation of a rock mass considering temperature can be written as:

[0129] ;

[0130] in, It is the thermal stress tensor;

[0131] This leads to the governing equations of saturated rock mass mechanics:

[0132] ;

[0133] in, It is the divergence operator; For the fourth-order elastic tensor under drainage conditions; For strain tensor; Biot coefficient; Pore ​​fluid pressure; For the Kronecker tensor; It is the linear thermal expansion coefficient; It is the bulk modulus of the drainage. The change in temperature , Indicates the current temperature. Indicates the reference temperature; For saturated rock mass density, , For fluid density, For the density of solid skeleton particles, It is porosity; It is the vector of gravitational acceleration; The strain tensor, under the assumption of small deformation, can be expressed as: , For fractional Laplace operators, For displacement;

[0134] The volumetric strain of rock can be expressed as the divergence of the displacement vector:

[0135] ;

[0136] in, , and They represent , and Displacement in three directions.

[0137] The fluid-thermal flow equations include the mass and energy conservation equations, and the equations of motion are Darcy's law and Fourier's law, respectively. The mass and energy conservation equations can be uniformly expressed as:

[0138] ;

[0139] in, and They represent fluid and heat, respectively. Indicates mass or energy; For source and sink items; For traffic; Let time be the constant. According to the theory of thermoelastic porous media, the changes in mass and specific heat within the system can be expressed in differential form as follows:

[0140] ;

[0141] ;

[0142] in, For Biot modulus, The coefficient of thermal expansion is given by the subscript. , and These represent solid skeleton, fluid, and saturated rock mass, respectively. The coefficient of thermal expansion of the fluid. For the total volumetric heat capacity, and These represent the solid skeleton heat capacity and the fluid heat capacity, respectively. and These are the specific heats of the saturated rock mass and the fluid, respectively. Additionally, some supplementary equations are as follows:

[0143] ;

[0144] in, , and These are the bulk moduli of saturated rock mass, rock skeleton, and fluid, respectively. Darcy's law applies:

[0145] ;

[0146] Thus, the mass conservation equation can be obtained:

[0147] ;

[0148] in, For the volumetric strain of the rock, , , and They represent , and Displacement in three directions; Biot modulus; For time; This is the total thermal expansion coefficient of the saturated rock mass. , The coefficient of thermal expansion of the fluid; For fluid viscosity; This represents the volume factor of the fluid, with 0 indicating a reference state. , For fluid density; For permeability tensor; For fluid sources / sinks.

[0149] Similarly, the total differential form of energy is:

[0150] ;

[0151] and Fourier's law (heat conduction equation) and heat convection equation After substituting the values, we can obtain the energy conservation equation:

[0152] ;

[0153] in, The bulk modulus of saturated rock mass; For the total volumetric heat capacity, , The density of the rock matrix particles. For the heat capacity of the solid framework, For fluid heat capacity; This is the thermal conductivity tensor; Enthalpy; It is a heat energy sink.

[0154] like Figure 3 As shown, the fracturing process of fractured-vuggy oil and gas reservoirs involves two main types of fractures: open and closed. The fracture stiffness of closed and open fractures... Defined as:

[0155] ;

[0156] in, The effective stress is in the normal direction; The crack aperture of the crack element; This refers to fluid pressure.

[0157] For closed cracks, crack stiffness depends on the effective normal stress and the surface roughness characteristics of the crack (i.e., the joint roughness coefficient JRC and the joint compressive strength JCS). The shear strength of rough closed cracks... With shear displacement The relationship is non-linear; when the peak shear strength is reached... The preceding peak shear displacement is (corresponding to the peak shear displacement) ), - The shear strength initially exhibits a near-linear relationship, but then decreases non-linearly with increasing shear displacement, eventually approaching the residual shear strength. The corresponding shear displacement becomes the residual shear displacement. For closed cracks, the classic Barton-Bandis model is used to calculate crack deformation, resulting in the following constitutive model for closed cracks:

[0158] ;

[0159] ;

[0160] ;

[0161] in, The initial normal stiffness; The maximum allowable degree of closure; This is the joint roughness coefficient; This refers to the joint compressive strength.

[0162] The constitutive model of an opening crack is a mathematical expression used to describe the mechanical response behavior of hydraulic cracks in an opening state under fluid pressure. Specifically, it characterizes the dynamic relationship between fluid pressure within the crack and crack aperture, as well as the crack's resistance to deformation (i.e., crack stiffness). When a crack opens under high fluid pressure, its stiffness is related to the fluid pressure within the crack. The size of the crack opening may depend on the deformation of the matrix around the crack under hydraulic loads on the crack surface. According to the Kelvin solution, the displacement at any point is closely related to the hydraulic load applied to another point. Therefore, the stiffness of an opening crack depends on the fluid pressure within the crack, the mechanical properties of the matrix (i.e., Young's modulus and Poisson's ratio), and the crack geometry (i.e., crack shape and size). For opening cracks, their stiffness is affected by crack size, geometry, and the magnitude of the stress. Currently, there is no effective calculation method for the stiffness of curved surface cracks. To accurately predict the deformation resistance of dynamically opening cracks, starting from the original definition of crack stiffness, a small perturbation stress is introduced into the curved surface crack. Then, the deformation of the crack under disturbed stress conditions was calculated using the discontinuous displacement method. This leads to the determination of the stiffness of the open crack. The project uses a numerical method for implicit solution. The relationship between the implicit crack stiffness and fluid pressure of the open crack is as follows:

[0163] ;

[0164] According to the definition of the derivative, crack stiffness is expressed as follows:

[0165] The fluid pressure p is obtained by solving the problem. f and stress Subsequently, the crack aperture vector of all crack elements was calculated using the discontinuous displacement method. for:

[0166] ;

[0167] ;

[0168] The influence coefficient matrix [C] is a function of the crack mesh geometry and is independent of stress and pressure. To be taken from A vector consisting of N discontinuous normal displacements. The crack aperture w of the crack element. f Equal to the discontinuity of normal displacement D n Introduce a small fluid pressure disturbance. (100 Pa used in this model) to generate another crack element crack aperture vector By solving the equation again:

[0169] ;

[0170] ;

[0171] in Then obtain the secant crack stiffness for the current iteration step:

[0172] ;

[0173] Therefore, the constitutive model of the opening crack is:

[0174] ;

[0175] ;

[0176] ;

[0177] ;

[0178] in, A vector consisting of the stiffnesses of N crack elements; The disturbance stress of the crack element; The deformation of the crack under disturbed stress conditions; This is a vector composed of the crack apertures of N crack elements after the addition of perturbation stress; It is a vector composed of the crack apertures of N crack elements under undisturbed stress. The discontinuous normal displacement vector of N crack elements after adding disturbance stress; For N crack elements, the normal displacement is discontinuous when there is no disturbance stress. The number of crack elements; This is the influence coefficient matrix; This is a vector representation of the stress borne by N crack elements; For disturbance stress, .

[0179] S4: Based on numerical computation grids and governing equations, fracture evolution analysis is performed to simulate the initiation, propagation, and reversal of hydraulic fractures, as well as their interaction with natural fractures and karst caves. Full three-dimensional fracture evolution includes fracture initiation, propagation (shear displacement), and reversal; different criteria are used to describe these processes.

[0180] S401. Determining the initiation conditions of hydraulic fractures: Based on the seepage-stress-heat transfer coupling equation, calculate the stress state at each potential initiation point in the current time step; use the maximum principal stress criterion to determine the tension initiation conditions: when the fluid pressure is sufficient to overcome the minimum in-situ principal stress and the tensile strength of the rock mass itself, or when the rock mass is subjected to the maximum effective stress. Greater than or equal to the tensile strength of the rock At that time, the cracks began to initiate:

[0181] ;

[0182] Simultaneously, considering the thermal stress effect caused by low-temperature fracturing fluid injection (i.e., when fracturing fluid at a lower temperature is injected into the reservoir, on the one hand, the fluid pressure increases due to volume expansion caused by the rise in fluid temperature, and on the other hand, the cooling stress is generated due to the decrease in rock temperature near the fracture wall, thus causing the initial Mohr's circle to gradually approach the failure envelope, and the rock undergoes shear failure when the stress Mohr's circle contacts the failure envelope), when the deviatoric stress level When the value is greater than or equal to 1.0, the natural crack is determined to have initiated shear slip cracking.

[0183] ;

[0184] in, This represents the deviatoric stress under the current condition; This represents the ultimate deviatoric stress at which the rock fails. This represents the current maximum effective principal stress; This represents the current minimum effective principal stress; This is the ultimate deviatoric stress at failure.

[0185] S402. Calculation of the stress intensity factor at the crack tip: The stress intensity factor reflects the magnitude of the stress singularity at the crack tip and can be expressed as a function of displacement discontinuity, crack geometry, and rock properties. The stress intensity factor at the crack tip is calculated using the displacement discontinuity method. :

[0186] ;

[0187] in, As an empirical constant, this embodiment takes... =0.806; For parameters related to rock properties, , , It is Young's modulus. It is Poisson's ratio; For discontinuities in displacement; subscript These represent the pure opening mode (i.e.) Type), sliding mode (i.e.) (type), scissor mode (i.e.) Type); Subscript , , These represent the directions of normal opening, strike-slip shear, and dip-slip shear, respectively. It is the distance from the centroid of the triangular unit to the crack propagation front;

[0188] S403. Determine the crack propagation direction and turning angle:

[0189] The maximum principal stress criterion is used to handle pure opening mode extension and mixed mode extension, when the equivalent stress intensity factor Achieving rock toughness hour( (As a fixed value), the crack will propagate; according to the Richard criterion, the crack orientation angle is calculated based on three stress intensity factors. The hybrid modes include I+II, I+III, and I+II+III.

[0190] Equivalent stress intensity factor for:

[0191] ;

[0192] ;

[0193] in, The crack deflection angle under triaxial loading; As a substitute variable; , and The stress intensity factor is related to the pure opening model, the sliding mode (i.e., type II), and the shear mode (i.e., type III), respectively.

[0194] Crack turning angle The calculation formula is:

[0195] ;

[0196] Among them, when When the crack turning angle is positive, and When the crack turning angle is negative;

[0197] S404. Calculate the crack propagation rate:

[0198] The crack propagation rate was calculated using the Paris scaling law model, and the propagation rate at the crack tip was correlated with the stress intensity factor amplitude to determine the crack propagation length per unit time step.

[0199] S405. Simulated stress interference between cracks and natural cracks / cavities:

[0200] The induced stress field is obtained by calculating the induced stress on the target mesh element by all crack and cavern elements using the indirect boundary element method, where the target mesh element is either a crack element or a cavern element. The induced stress is expressed as:

[0201] ;

[0202] in, , and These represent all crack and cavern elements within the target mesh element. In the local coordinate system, the normal induced stress, the shear induced stress along the strike direction, and the shear induced stress along the dip direction; This is the influence coefficient; Indicates a crack or cave unit Discontinuous displacement of the crack normal or virtual stress of the cavity normal in the target mesh element The influence coefficient of the normal induced stress generated at the location; Indicates a crack or cave unit Discontinuous displacement of the crack normal or virtual stress of the cavity normal in the target mesh element The influence coefficient of shear-induced stress generated at the location; Indicates a crack or cave unit Discontinuous displacement of the crack normal or virtual stress of the cavity normal in the target mesh element The influence coefficient of the tendency shear-induced stress generated at the location; Indicates a crack or cave unit The direction of cracks, discontinuous displacement, or cavern direction at the location are virtual stresses in the target mesh element. The influence coefficient of the normal induced stress generated at the location; Indicates a crack or cave unit The direction of cracks, discontinuous displacement, or cavern direction at the location are virtual stresses in the target mesh element. The influence coefficient of shear-induced stress generated at the location; Indicates a crack or cave unit The direction of cracks, discontinuous displacement, or cavern direction at the location are virtual stresses in the target mesh element. The influence coefficient of the tendency shear-induced stress generated at the location; Indicates a crack or cave unit Cracks at the location tend to exhibit discontinuous displacement or cavitation tend to cause virtual stress in the target mesh element. The influence coefficient of the normal induced stress generated at the location; Indicates a crack or cave unit Cracks at the location tend to exhibit discontinuous displacement or cavitation tend to cause virtual stress in the target mesh element. The influence coefficient of shear-induced stress generated at the location; Indicates a crack or cave unit Cracks at the location tend to exhibit discontinuous displacement or cavitation tend to cause virtual stress in the target mesh element. The influence coefficient of the tendency shear-induced stress generated at the location; for crack elements, , and Representing units respectively The displacement is discontinuous along the normal, strike, and dip directions; for a cave unit, , and Representing units respectively Virtual stress along the normal, strike, and dip directions; This represents the total number of mesh elements for fracture and cavern elements;

[0203] For the target mesh element, force balance requires that the induced stress equals the net stress, which is the internal fluid pressure minus the in-situ stress, i.e.:

[0204] ;

[0205] in, Target mesh cell Fluid pressure in the middle; It is the pore pressure in the matrix mesh into which the target mesh element is embedded. Background stress of the matrix mesh;

[0206] S406. Based on the induced stress field calculated in S405, determine the cross-behavior pattern when hydraulic fractures approach natural fractures or karst caves; wherein, the cross-behavior pattern includes penetration mode, fracture arrest mode, shear slip mode, fusion mode, and attraction or repulsion mode. The specific determination process is as follows: when the hydraulic fracture tip satisfies the restart condition on the other side of the natural fracture, and the natural fracture does not undergo shear slip, it is determined to be a penetration mode; when the natural fracture undergoes shear slip, fracturing fluid flows into the natural fracture, causing pressure dissipation, and the hydraulic fracture tip becomes blunt, it is determined to be a fracture arrest mode; when the shear stress on the natural fracture exceeds its shear strength, shear slip occurs, activating the permeability of the natural fracture, it is determined to be a shear slip mode; when two coplanar fractures approach each other, it is determined to be a fusion mode, with the fracture tips merging to form a new extension front, it is determined to be a fusion mode; based on the stress concentration characteristics around the karst cave, determine the interaction type between the hydraulic fracture and the karst cave, including attraction followed by gradual approach, attraction followed by repulsion, and repulsion followed by moving away from the karst cave body, it is determined to be an attraction or repulsion mode.

[0207] More specifically, such as Figure 4 As shown, S4 also includes using an adaptive mesh method to track the full three-dimensional crack propagation front and its geometry. The specific process is as follows:

[0208] S407. For example Figure 4 As shown in (a), a new propagation leading edge is determined: based on the crack propagation direction determined in S403 and the crack propagation velocity calculated in S404, a propagation vector set at the crack tip at the current time step is generated. , , , , , Connect the endpoints of the expansion vector to the previous expansion step. At each endpoint of the crack leading edge, obtain new information about the crack. The leading edge contour line is continuously expanded; the spatial coordinates of each node of the new leading edge and the number of the crack element to which it belongs are recorded;

[0209] S408. (e.g.) Figure 4 As shown in (b), the quality of the extended leading edge mesh is adjusted:

[0210] Based on the geometric characteristics of the new leading edge, the quality of the extended leading edge mesh is adjusted in three cases:

[0211] Case a: When the included angle between two adjacent sides When the angle is greater than 120°, connect the vertex of the included angle to the midpoint of the opposite side. The triangular unit is split into two computational units to avoid the unit being too narrow and long, which would affect the computational efficiency.

[0212] Case b: When the length of the expansion vector is less than 20% of the maximum expansion step size, move each endpoint of the previous step directly to the corresponding new endpoint according to the expansion vector, and delete the original endpoints to avoid excessively dense grids that would cause a rapid increase in degrees of freedom and a significant decrease in computational efficiency.

[0213] Case c: When the expansion vector exceeds 50% of the maximum expansion step size, find the new leading edge endpoint according to the expansion vector, and add a new endpoint at the midpoint of the new expansion vector. and ), for The midpoint of the crack is used to treat a single-step expansion as a two-step uniform expansion, thus avoiding excessively narrow cracks that could affect computational efficiency.

[0214] S409. Smoothing the leading edge of the crack: The Savitzky-Golay filter fitting method is used to smooth the latest leading edge node; the sawtooth fluctuation of the leading edge node is eliminated by local polynomial fitting to maintain the smooth and continuous characteristics of the crack leading edge; distorted mesh elements are avoided to ensure the computational stability of subsequent expansion steps.

[0215] S410: As Figure 5 (a) and Figure 5 As shown in (b), the penetration, arrest, and fusion of intersecting cracks are handled by truncating the natural crack surface and re-meshing the local mesh along the intersection line:

[0216] Based on the cross-behavior patterns determined by S406, mesh processing is performed on the interactions of different types of cracks:

[0217] Penetration mode: When a hydraulic fracture restarts on the other side of a natural fracture, no additional processing of the mesh is required, allowing the fracture to continue to propagate.

[0218] Crack arrest mode: When a hydraulic crack stops at a natural crack, first determine the crack intersection line, and then perform grid truncation treatment with the intersection line as a constraint to stop the further propagation of the hydraulic crack.

[0219] Fusion mode: When two coplanar cracks approach and merge, the overlapping mesh is deleted to merge the leading edges of the two cracks and form a new continuous crack surface.

[0220] Shear slip mode: When natural cracks undergo shear slip, new crack elements are generated at the slip surface, updating the topology of the crack network;

[0221] S411: Update crack network topology and computational mesh: Integrate the mesh cells adjusted in S408, the leading edge nodes smoothed in S409, and the cross crack mesh processed in S410; generate the complete crack network triangular mesh after the current expansion step; update the numbering, node coordinates, and connectivity of all crack cells to provide an updated computational mesh for the S401-S406 analysis in the next expansion step.

[0222] S5: Output the results of crack morphology, aperture distribution, and flow field distribution during the crack evolution process.

[0223] In practice, S3 calculates the flow field distribution within the crack, while S4 calculates the crack propagation morphology and crack aperture distribution. The calculation results obtained from S3 and S4 are stored in the computer in the binary format required by TECPLOT, and can then be used for 3D visualization and analysis of the calculation results via TECPLOT.

[0224] It is understood that the same or similar parts in the above embodiments can be referred to each other, and the contents not described in detail in some embodiments can be referred to the same or similar contents in other embodiments.

[0225] Although embodiments of the present invention have been shown and described above, it is understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those skilled in the art can make changes, modifications, substitutions and variations to the above embodiments within the scope of the present invention.

[0226] To further illustrate the technical solution of this invention, several application examples are provided, as follows:

[0227] Application Example 1

[0228] This study focuses on the propagation of fractures induced by hydrocarbon generation in source rocks. Based on the proposed three-dimensional numerical simulation method for the mechanical behavior of fracture-fracture and fracture-cavity intersections, the propagation process of multiple fractures is simulated. The evolution of a fracture swarm consisting of 50 initial fractures is studied. The study area measures 1.8m × 1.8m × 1.8m and is discretized into 3375 rectangular grids of size 0.12m × 0.12m × 0.12m. The shale source rock is located at a relatively shallow depth, and its in-situ stress state is under a reverse fault stress mechanism. The main calculation parameters are set as follows: maximum horizontal principal stress... Minimum horizontal principal stress Vertical stress Assuming the initial pore pressure and the initial fluid pressure within the fracture are the same, both being 1.0 MPa, the Young's modulus is 39 GPa, the Poisson's ratio is 0.25, the Biot coefficient is 0.85, and the permeability is 10... -4 mD, density is 2200 kg / m³ 3The fluid has a bulk modulus of 45 GPa, a porosity of 0.09, a viscosity of 1.0 cP, and a density of 850 kg / m³. 3 The bulk modulus is 2.0 GPa, the initial crack inclination is 0.01 mm, the crack normal stiffness is 20 GPa / m, and the rock fracture toughness is 1.0 MPa. The injection rate for each fracture is 4.0 × 10⁻⁶. -7 m 3 / s.

[0229] like Figure 6 As shown, the three-dimensional morphology and crack aperture distribution of crack clusters under different preset crack network morphologies are illustrated. Figure 6 As shown in (a), the propagation status of each crack after 2665 seconds was recorded using crack surface area as an indicator. All cases underwent more than 35 crack propagation steps. The histogram of individual crack surface area reflects the competitive growth of crack swarms, with a few large cracks dominating. As the crack swarm density increases, the preferential propagation of a few cracks becomes more pronounced. Under the condition of the same crack density in the study area, the initial crack morphology is another reason for the differential growth of crack swarms. For fluid-driven cracks, the crack aperture at the crack tip determines the degree of stress concentration, thus determining the crack propagation potential. When cracks intersect, the crack tip becomes blunted, causing the crack propagation rate to gradually slow down. Isolated cracks will stop propagating due to the strong stress shielding effect caused by other widely distributed neighboring cracks. This calculation answers why a three-dimensional non-planar intersecting crack network is generated during hydrocarbon generation, determines five typical interaction modes of mechanical interaction of three-dimensional curved surface cracks, and clarifies the mechanical mechanism behind the alternating, differential, and gradual propagation rate behaviors of multiple three-dimensional cracks during propagation. Figure 6 As shown in (b), the competitive and collaborative evolution characteristics of three-dimensional nonplanar crack groups under different initial crack geometries are as follows: Figure 7 As shown, where, Figure 7 (a) shows the competitive propagation results of the initial crack at a dip angle of 0°; Figure 7 (b) shows the competitive propagation results of the initial crack at a dip angle of 30°; Figure 7 (c) shows the competitive propagation results of the initial crack at a dip angle of 45°; Figure 7 (d) shows the competitive propagation results of the initial crack at a dip angle of 60°; Figure 7 (e) shows the competitive propagation results of the initial crack at a 90° dip angle; Figure 7 (f) shows the competitive propagation results of the initial cracks with randomly distributed dip angles.

[0230] Application Example 2

[0231] The numerical simulation method for the cross-mechanical behavior of cracks and fissures and cavities in full three dimensions, provided by this invention, is used to simulate the cross-mechanical behavior of cracks and fissures, and cracks and cavities in full three dimensions. For example... Figure 8 As shown, the simulation parameters used are set as follows: Considering segmented fracturing in a 50m×50m×50m three-dimensional space, the matrix is ​​divided into 20×20×20=8000 Cartesian structured grids; the considered horizontal well section is 30m in length; the matrix porosity is 0.25; the matrix permeability is 0.5mD; the initial pore water pressure of the matrix is ​​25MPa; the initial fluid pressure inside the wellbore is 27MPa; the wellbore radius is 0.25m; the Young's modulus is 38.8GPa; the Poisson's ratio is 0.15; and the rock fracture toughness is 3.0MPa. In-situ stress is MPa, injection rate is The fluid viscosity is 5.0 cP. The horizontal well section has two pre-fabricated cutting fractures spaced 6 m apart. On both sides of the horizontal well are four spherical and four ellipsoidal cavities, each with an internal pressure of 27.5 MPa. (The last sentence appears to be incomplete and unrelated to the preceding text.) Figure 8 As shown, the first cave is an ellipsoid, with its center located at... m, the semi-axis lengths in the three directions are m; the second cave is a sphere, with its center located at... The first cave is 2.5m in diameter and has a radius of 2.5m. The third cave is spherical, with its center located at... The fourth cave is an ellipsoid with a radius of 2.5m; its center is located at... m, the semi-axis lengths in the three directions are m.

[0232] like Figure 9 As shown, the evolution of fractures due to the mutual interference between caves, fractures, and wellbore in the presence of both caves and natural fissures is illustrated. Caves of different morphologies exhibit an attractive force on adjacent fractures, with spherical caves showing a stronger attraction than ellipsoidal caves, leading to a more pronounced deflection of the fractures. This phenomenon can be attributed to the stress concentration characteristics around caves of different geometries. The stress shielding effect between adjacent fractures also causes them to move away from each other during propagation, and pre-existing natural fissures influence the propagation process and final geometry of the fractures. The attractive force exerted by the caves on the propagating fractures is weakened. This is partly attributed to the slippage of pre-existing natural fissures, which induces additional stress disturbances, potentially hindering the propagation of the fractures towards the caves. Specifically, according to... Figure 9 (a) It can be seen that the influence of far-field stress on the initial crack propagation direction changes, and the attraction of the cavern to the crack is not significant; according to Figure 9 (b) It can be seen that as the cracks gradually grow, the influence of cavern-induced stress on the expanding cracks gradually increases; according to Figure 9(c) It can be seen that the cracks are attracted by the karst caves and undergo significant diversion and expansion; by Figure 9 (d) It can be seen that the cracks are affected by the stress induced by the caves, and are gradually moved away from the caves due to repulsion; the hydraulic cracks have a strong penetrating ability, and the result of their interaction with the natural cracks is penetration; the interaction behavior of the hydraulic cracks with the four cave bodies is different, mainly manifested as insignificant attraction, gradual approach after attraction, repulsion after attraction, and moving away from the cave body after being repelled.

Claims

1. A numerical simulation method for the mechanical behavior of cross-slit and cross-hole joints in full three dimensions, characterized in that, include: S1: Obtain the geometric information of natural fractures and caverns in fracture-vuggy reservoirs, and construct a full three-dimensional geometric model of fractures and caverns based on the geometric information; S2: Mesh the full three-dimensional geometric model to create a full three-dimensional numerical computation mesh that includes the matrix, crack network, and cavern body; S3: Determine the governing equations for simulating the mechanical behavior of cracks and caverns. The governing equations shall include at least the seepage-stress-heat transfer coupling equations and the crack constitutive model. The crack constitutive models include closed crack constitutive models and open crack constitutive models, and the crack stiffness... Represented as: ； in, The effective stress is in the normal direction; The crack aperture of the crack element; For fluid pressure; The constitutive model for a closed crack is: ； ； ； in, The initial normal stiffness; The maximum allowable degree of closure; This is the joint roughness coefficient; Joint compressive strength; The constitutive model for an open-type crack is: ； ； ； ； in, A vector consisting of the stiffnesses of N crack elements; The disturbance stress of the crack element; The deformation of the crack under disturbed stress conditions; This is a vector composed of the crack apertures of N crack elements after the addition of perturbation stress; It is a vector composed of the crack apertures of N crack elements under undisturbed stress. The discontinuous normal displacement vector of N crack elements after adding disturbance stress; For N crack elements, the normal displacement is discontinuous when there is no disturbance stress. The number of crack elements; This is the influence coefficient matrix; This is a vector representation of the stress borne by N crack elements; For disturbance stress, ; S4: Based on a full three-dimensional numerical computation grid and governing equations, crack evolution analysis is performed to simulate the initiation, propagation, and direction of hydraulic cracks, as well as their cross-mechanical behavior with natural cracks and karst caves. S5: Output the results of crack morphology, aperture distribution, and flow field distribution during the crack evolution process.

2. The method according to claim 1, characterized in that, The geometric information is multi-source heterogeneous geological data, including points, lines, surfaces, and volumes; among them, points are core observation and downhole imaging data; lines are logging and tracer monitoring data along the wellbore; surfaces are fiber optic monitoring and seismic profile data; and volumes are three-dimensional seismic inversion and comprehensive interpretation data.

3. The method according to claim 2, characterized in that, The specific process of S1 is as follows: S11: Acquire multi-source heterogeneous geological data of fractured-vuggy reservoirs, perform data fusion processing, and extract the boundary point set of natural fractures and the spatial location and morphological parameters of karst bodies; S12: Based on the extracted set of natural crack boundary points, the natural cracks are simplified into polygonal crack surfaces, and the vertex coordinates and spatial orientation of each polygon are determined. Then, a three-dimensional surface model of the natural cracks is constructed using parametric surface fitting technology. S13: Based on the extracted spatial location and morphological parameters of the karst cave, the karst cave is simplified into an ellipsoid, and the coordinates of the center of each ellipsoid, the length of the semi-axis in three directions, and the normal vector of the major axis symmetry plane are determined. A three-dimensional closed ellipsoid model of the karst cave is generated using closed surface construction technology. S14: Integrate the three-dimensional curved surface models of all natural cracks and the three-dimensional closed ellipsoidal models of the cave bodies into the same spatial coordinate system to form a full three-dimensional geometric model that includes the crack network and the cave system.

4. The method according to claim 3, characterized in that, The specific process of S2 is as follows: S21: Treat the matrix rock mass of fractured-vuggy reservoirs as a continuous porous medium, and use a Cartesian structured grid to divide the matrix region to generate a background grid; S22: Based on the three-dimensional surface model of the natural crack constructed by S12, the boundary closed contour line of each crack is extracted; the crack boundary point set is interpolated and fitted to reconstruct the surface morphology of the crack; the built-in function of the mesh generation tool is used to triangularly mesh each crack surface to generate a crack mesh composed of triangular elements. S23: A three-dimensional closed ellipsoidal model of the cave body constructed based on S13. The geometric parameters of the cave are determined according to the given coordinates of the center of the sphere, the lengths of the semi-axis in three directions, and the normal vector of the major axis symmetry plane. The ellipsoid is uniformly divided into a preset number of symmetrical arc segments using the built-in functions of the mesh generation tool. Based on symmetrical arc segments, construct the boundary lines of a preset number of sub-surfaces, and triangulate each sub-surface to generate a closed triangular mesh of the cave body; use the built-in function of the mesh generation tool to triangulate each cave surface to generate a cave mesh composed of triangular units. S24: Treat the fracture triangular mesh and the cave triangular mesh as independent embedded objects; through a cyclic search algorithm, embed each fracture unit and cave unit into the matrix background mesh generated in S21 one by one, and establish the spatial correspondence between fracture and matrix, fracture and fracture, fracture and cave, and cave and matrix. Based on the spatial correspondence, establish a connectivity table for all grid cells and calculate the conductivity coefficient between adjacent cells; S25: Integrate the matrix background mesh, crack mesh, and cave mesh to obtain a full three-dimensional numerical calculation mesh.

5. The method according to claim 1, characterized in that, In S3, the seepage-stress-heat transfer coupling equation includes: The governing equations of saturated rock mass mechanical equilibrium are: ； in, It is the divergence operator; For the fourth-order elastic tensor under drainage conditions; Biot coefficient; Pore ​​fluid pressure; For the Kronecker tensor; It is the linear thermal expansion coefficient; It is the bulk modulus of the drainage. The change in temperature , Indicates the current temperature. Indicates the reference temperature; For saturated rock mass density, , For fluid density, For the density of solid skeleton particles, It is porosity; It is the vector of gravitational acceleration; For strain tensor, , For fractional Laplace operators, For displacement; mass conservation equation: ； in, For the volumetric strain of the rock, , , and They represent , and Displacement in three directions; Biot modulus; For time; This is the total thermal expansion coefficient of the saturated rock mass. , The coefficient of thermal expansion of the fluid; For fluid viscosity; This represents the volume factor of the fluid, with 0 indicating a reference state. ; For permeability tensor; For fluid sources / sinks; Energy conservation equation: ； in, The bulk modulus of saturated rock mass; For the total volumetric heat capacity, , For the heat capacity of the solid framework, For fluid heat capacity; This is the thermal conductivity tensor; Enthalpy; For thermal energy / sink.

6. The method according to claim 1, characterized in that, The specific process of S4 is as follows: S401. Determining the initiation conditions of hydraulic fractures: Based on the seepage-stress-heat transfer coupling equation, the stress state at each potential crack initiation point is calculated at the current time step; the maximum principal stress criterion is used to determine the tension initiation condition: when the rock mass experiences the maximum effective stress... Greater than or equal to the tensile strength of the rock At that time, the crack initiation occurred; simultaneously, considering the thermal stress effect caused by the injection of cryogenic fracturing fluid, when the deviatoric stress level... A value greater than or equal to 1.0 indicates that a natural crack has initiated shear slip cracking; among which, the deviatoric stress level... The calculation formula is as follows: ； in, This represents the deviatoric stress under the current condition; This represents the ultimate deviatoric stress at which the rock fails. This represents the current maximum effective principal stress; This represents the current minimum effective principal stress; This represents the ultimate deviatoric stress at failure. S402. Calculate the stress intensity factor at the crack tip: Based on displacement discontinuity, crack geometry, and rock properties, the stress intensity factor at the crack tip is calculated using the displacement discontinuity method. : ； in, These are empirical constants; For parameters related to rock properties, , , It is Young's modulus. It is Poisson's ratio; For discontinuities in displacement; subscript These represent the open mode, slide mode, and scissor mode, respectively; subscript , , These represent the directions of normal opening, strike-slip shear, and dip-slip shear, respectively. It is the distance from the centroid of the triangular unit to the crack propagation front; S403. Determine the crack propagation direction and turning angle: The maximum principal stress criterion is used to handle pure opening mode extension and mixed mode extension, when the equivalent stress intensity factor Achieving rock toughness At this time, the crack will propagate; according to the Richard criterion, the crack orientation angle is calculated based on three stress intensity factors. The hybrid modes include I+II, I+III, and I+II+III. Equivalent stress intensity factor for: ； ； in, The crack deflection angle under triaxial loading; As a substitute variable; Crack turning angle The calculation formula is: ； Among them, when When the crack turning angle is positive, and When the crack turning angle is negative; S404. The crack propagation rate is calculated using the Paris scaling law model; S405. Calculate the induced stress on the target mesh element by all crack and cavity elements using the indirect boundary element method to obtain the induced stress field, where the target mesh element is either a crack element or a cavity element; the target mesh element... The induced stress is expressed as: ； in, , and These represent all crack and cavern elements within the target mesh element. In the local coordinate system, the normal induced stress, the shear induced stress along the strike direction, and the shear induced stress along the dip direction; The influence coefficient is as follows: Indicates a crack or cave unit The corresponding discontinuous displacement in the crack normal or the virtual stress in the cave normal at the target mesh element The influence coefficient of the normal induced stress generated at the location; Indicates a crack or cave unit The corresponding discontinuous displacement in the crack normal or the virtual stress in the cave normal at the target mesh element The influence coefficient of shear-induced stress generated at the location; Indicates a crack or cave unit The corresponding discontinuous displacement in the crack normal or the virtual stress in the cave normal at the target mesh element The influence coefficient of the tendency shear-induced stress generated at the location; Indicates a crack or cave unit The corresponding crack direction, discontinuous displacement, or cavern direction virtual stress in the target mesh element The influence coefficient of the normal induced stress generated at the location; Indicates a crack or cave unit The corresponding crack direction, discontinuous displacement, or cavern direction virtual stress in the target mesh element The influence coefficient of shear-induced stress generated at the location; Indicates a crack or cave unit The corresponding crack direction, discontinuous displacement, or cavern direction virtual stress in the target mesh element The influence coefficient of the tendency shear-induced stress generated at the location; Indicates a crack or cave unit The corresponding crack tends to discontinuous displacement or cavitation tendency, and the virtual stress is in the target mesh element. The influence coefficient of the normal induced stress generated at the location; Indicates a crack or cave unit The corresponding crack tends to discontinuous displacement or cavitation tendency, and the virtual stress is in the target mesh element. The influence coefficient of shear-induced stress generated at the location; Indicates a crack or cave unit The corresponding crack tends to discontinuous displacement or cavitation tendency, and the virtual stress is in the target mesh element. The influence coefficient of the tendency shear-induced stress generated at the location; for crack elements, , and Representing units respectively The displacement is discontinuous along the normal, strike, and dip directions; for a cave unit, , and Representing units respectively Virtual stress along the normal, strike, and dip directions; This represents the total number of mesh elements for fracture and cavern elements; For the target mesh element, force balance requires that the induced stress equals the net stress, which is the internal fluid pressure minus the in-situ stress, i.e.: ； in, Target mesh cell Fluid pressure in the middle; It is the pore pressure in the matrix mesh into which the target mesh element is embedded. Background stress of the matrix mesh; S406. Based on the induced stress field calculated in S405, determine the cross-behavior pattern when hydraulic fractures approach natural fractures or karst caves; wherein, the cross-behavior pattern includes penetration mode, fracture arrest mode, shear slip mode, fusion mode, and attraction or repulsion mode.

7. The method according to claim 6, characterized in that, The judgment process of S406 is as follows: When the tip of a hydraulic fracture satisfies the restart condition on the other side of a natural fracture, and the natural fracture does not undergo shear slip, it is determined to be a penetration mode. When shear slip occurs in a natural fracture, fracturing fluid flows into the natural fracture, causing pressure dissipation and blunting of the hydraulic fracture tip, which is determined to be the fracturing termination mode. When the shear stress on a natural crack exceeds its shear strength, shear slip occurs, activating the permeability of the natural crack, and this is identified as a shear slip mode. When two coplanar cracks approach each other, it is determined to be a fusion, and the crack tips merge to form a new extension leading edge, which is determined to be a fusion mode; Based on the stress concentration characteristics around the cave, the interaction type between the hydraulic fissure and the cave is determined, including attraction followed by gradual approach, attraction followed by repulsion, and repulsion followed by moving away from the cave body, and is determined to be an attraction or repulsion mode.

8. The method according to claim 6, characterized in that, The S4 also includes using an adaptive mesh method to track the full three-dimensional crack propagation front and its geometry. The specific process is as follows: S407. Based on the crack propagation direction determined in S403 and the crack propagation velocity calculated in S404, generate the propagation vector set at the crack tip at the current time step; Connect the endpoints of the expansion vector to the endpoints of the crack leading edge in the previous expansion step to obtain the new expansion leading edge contour of the crack; record the spatial coordinates of each node of the new leading edge and the number of the crack element to which it belongs; S408. Based on the geometric characteristics of the new leading edge, the quality of the extended leading edge mesh is adjusted, specifically as follows: When the included angle between two adjacent sides is greater than a preset angle, connect the vertex of the included angle and the midpoint of the opposite side to split the triangular unit into two calculation units; When the length of the expansion vector is less than the first preset threshold of the maximum expansion step, the endpoints of the previous step are moved directly to the corresponding new endpoints according to the expansion vector, and the original endpoints are deleted. When the expansion vector exceeds the second preset threshold of the maximum expansion step size, a new leading edge endpoint is found according to the expansion vector, and a new endpoint is added at the midpoint of the new expansion vector, treating the single-step expansion as a two-step uniform expansion. S409. Smoothing the extended leading edge: The Savitzky-Golay filter fitting method is used to smooth the latest extended leading edge node; S410. Handling Penetration, Arrest, and Fusion of Intersecting Cracks: Based on the intersection behavior patterns determined in S406, mesh processing is performed on the interactions of different types of cracks. Penetration mode: When a hydraulic fracture restarts on the other side of a natural fracture, no additional processing of the mesh is required, allowing the fracture to continue to propagate. Crack arrest mode: When a hydraulic crack stops at a natural crack, first determine the crack intersection line, and then perform grid truncation treatment with the intersection line as a constraint to stop the further propagation of the hydraulic crack. Fusion mode: When two coplanar cracks approach and merge, the overlapping mesh is deleted to merge the leading edges of the two cracks and form a new continuous crack surface. Shear slip mode: When natural cracks undergo shear slip, new crack elements are generated at the slip surface, updating the topology of the crack network; S411: Update the crack network topology and computational mesh: Integrate the mesh cells adjusted in S408, the leading edge nodes smoothed in S409, and the cross crack mesh processed in S410 to generate the complete crack network triangular mesh after the current expansion step; update the numbering, node coordinates, and connection relationships of all crack cells.

Images