Structural network fracture modeling method based on multi-scale factor constraint

Abstract

The invention relates to the technical field of oil and gas reservoir development and geological modeling, and discloses a multi-scale factor constraint-based tectonic network fracture modeling method, which comprises the following steps of: synthesizing a Riemannian metric tensor field on the basis of earthquake, logging and geomechanics data, and defining non-Euclidean distance cost of a fracture expanded in an anisotropic medium; initial seed points are screened according to the elastic strain energy density, and initial growth potential energy in a limited range is distributed; solving the eikonal equation by using an anisotropic fast marching algorithm to carry out wavefront competitive growth, and dividing a grid region into generalized Voronoi units; identifying a wavefront contact interface, and extracting gradient features to judge a fusion or truncation type so as to establish fracture topological connection; and finally, tracking a geodesic line path along an anti-gradient direction to generate a three-dimensional discrete fracture network. According to the method, macro and micro constraints are unified through Riemannian geometry, clear physical significance is given to the fracture by utilizing an energy mechanism, and automatic and accurate construction of the complex fracture network topology structure is realized.

Description

Technical Field

[0001] This invention relates to the field of oil and gas reservoir development and geological modeling technology, specifically to a structural network fracture modeling method based on multi-scale factor constraints. Background Technology

[0002] Fractures serve as the primary reservoir space and seepage channels in unconventional oil and gas reservoirs such as tight sandstone, shale, and carbonate rocks. Accurate description of their three-dimensional spatial distribution characteristics is crucial for oil and gas resource evaluation and development planning. Currently, discrete fracture network modeling is the main method for characterizing fracture systems; however, existing technologies still have certain limitations when dealing with fracture generation and propagation under complex geological conditions.

[0003] First, in terms of integrating fracture morphology with geological properties, existing modeling methods struggle to effectively balance the multi-scale constraints of macroscopic structure and micromechanics. Common approaches treat fault information and rock mechanics parameters as independent control variables, or rely solely on simple geometric interpolation to constrain fracture orientation. This approach ignores the pervasive anisotropy of the subsurface medium, resulting in fractures that are typically regular planar or disk-shaped. It fails to capture the complex bending deformation and non-planar propagation characteristics of fractures under the combined effects of fault disturbance, geostress deflection, and lithological changes, thus reducing the model's geometric realism.

[0004] Secondly, traditional stochastic fracture modeling techniques rely on geometric simulations based on probability statistics to determine the location and size of fractures. The location of the fracture center point is usually randomly generated by a Poisson process or fractal distribution, while the fracture length and pore size are sampled according to statistical laws (such as power-law distribution). This purely geometric and statistical method lacks a clear geomechanical physical driving mechanism, fails to fully consider the energy dissipation principle of rock fracture, resulting in deviations between the simulated fracture distribution and the actual stress concentration area or high strain energy area, and making it difficult to explain the limiting factors of fracture size from a physical perspective.

[0005] Furthermore, in constructing fracture network topologies, existing techniques mostly rely on spatial intersection operations of geometric objects to determine connectivity. This approach, based on geometric Boolean operations, is not only computationally complex and inefficient, but also struggles to simulate the dynamic interactions during fracture growth. When determining the morphology of intersecting fractures, it often simply treats them as continuous intersections, lacking physical criteria based on wavefront propagation dynamics to distinguish between fracture fusion, branching, or truncation (T-type, X-type nodes), thus affecting the reliability of fracture network connectivity analysis and subsequent fluid flow simulation. Summary of the Invention

[0006] To address the shortcomings of existing technologies, this invention provides a method for constructing network fracture modeling based on multi-scale factor constraints. This method solves the problems of existing technologies, such as the difficulty in unifying and integrating macroscopic structural and microscopic mechanical multi-scale constraints to accurately simulate the non-planar propagation path of fractures in anisotropic media, the lack of means for controlling fracture nucleation distribution and size based on physical energy mechanisms, and the difficulty in efficiently and automatically constructing complex fracture network topology connections that conform to physical laws based on dynamic characteristics.

[0007] To achieve the above objectives, the present invention provides the following technical solution: This invention provides a structural network fracture modeling method based on multi-scale factor constraints. First, a unified mapping mechanism for multi-scale geological data is established. Based on seismic attributes, well logging inversion, and geomechanical simulation data, a structural steering tensor reflecting macroscopic fault control and a rock mechanics constitutive tensor reflecting microscopic geostress and rock brittleness characteristics are constructed, respectively. These tensors are then fused using the principle of linear weighted superposition, synthesizing a Riemannian metric tensor field at the three-dimensional grid nodes. This tensor field and the diffusion tensor field obtained from its inverse operation define the non-Euclidean distance cost of fracture propagation in anisotropic media, providing an anisotropic geometric background for fracture growth.

[0008] In terms of fracture nucleation and growth control, this invention introduces an energy-driven mechanism. The elastic strain energy density is calculated based on the stress and strain tensors obtained from geomechanical simulations, thus selecting initial seed points for fracture initiation. Initial growth potential energy is then allocated based on a critical fracture energy threshold. Subsequently, an anisotropic fast marching algorithm is used to solve the equations, simulating the competitive growth process of wavefronts in the Riemannian metric space. Wavefront propagation is controlled by the diffusion tensor field and naturally stops when the accumulated non-Euclidean distance cost exceeds the initial growth potential energy, thereby providing a clear physical constraint on fracture size.

[0009] To construct the fracture network topology, this invention employs a dynamic determination technique based on generalized Voronoi elements. The wavefront propagation process divides the grid region into elements controlled by different seed points, and the common boundary between these elements is the wavefront contact interface. By analyzing the spatial gradient characteristics of the wavefront arrival time field on both sides of the contact interface, the gradient angle and positive hindrance coefficient are calculated to quantitatively determine whether fractures merge or terminate, thereby automatically establishing fracture topology relationships that conform to physical laws.

[0010] Ultimately, this invention starts from the contact interface that establishes topological connections, traces the geodesic path along the inverse gradient direction of the wavefront arrival time field to the initial seed point, and materializes the backtracking streamline into a three-dimensional discrete fracture network model with hydraulic gap width properties. This method utilizes Riemannian manifold theory to unify multi-scale geological constraints, ensuring that the generated fracture path is the path with minimum energy consumption, and achieving a realistic simulation of fracture networks in complex underground structures.

[0011] This invention provides a method for constructing network crack modeling based on multi-scale factor constraints. It has the following beneficial effects: 1. This invention integrates macroscopic fault structure properties with microscopic rock mechanics parameters into a unified Riemannian metric tensor field, defining anisotropic non-Euclidean distance costs. This makes the generated fracture paths no longer simple geometric connections, but geodesics that can adaptively respond to fault disturbances, geostress deflection, and changes in rock brittleness. This ensures that the fracture model conforms to the anisotropic physical characteristics of the subsurface medium in both macroscopic distribution trends and local geometry.

[0012] 2. This invention abandons the traditional approach of relying on probability statistical distributions to determine crack location and size in stochastic modeling. Instead, it introduces elastic strain energy density as the nucleation driving force and uses initial growth potential energy to limit wavefront propagation range. This energy-mechanism-based modeling method ensures that cracks preferentially form in high strain energy regions, and that crack size is controlled by the local energy dissipation limit, thus improving the accuracy of the model's prediction of stress concentration zones and densely developed crack areas.

[0013] 3. This invention utilizes an anisotropic fast-moving algorithm to simulate the competitive growth of multi-source wavefronts, and automatically determines the fusion or truncation relationship between fractures by using the gradient characteristics and hindrance coefficients of the wavefront contact interface. This method avoids complex geometric Boolean operations and can efficiently and robustly automatically generate three-dimensional discrete fracture networks containing complex topological features such as T-shaped truncation and X-shaped intersections, thereby improving the reliability of fracture connectivity analysis and fluid flow simulation. Attached Figure Description

[0014] Figure 1 The overall flowchart of the network fracture modeling method based on multi-scale factor constraints provided in the embodiments of the present invention is shown below. Figure 2 This is a flowchart illustrating the synthesis of the Riemannian metric tensor field and the diffusion tensor field in an embodiment of the present invention. Figure 3 This is a detailed flowchart of wavefront competitive growth and dynamic topology connection construction in an embodiment of the present invention. Detailed Implementation

[0015] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0016] See attached document Figure 1 One embodiment of the present invention can run on a computer device including a processor and a memory. The memory stores computer-executable instructions, which the processor executes to perform the following steps. The computer device needs to access a storage device for reading seismic data volumes, well logging data, and geomechanical simulation results data. The method includes the following steps: Step S100: Construct a multi-scale anisotropic Riemannian metric tensor field The computer equipment first establishes a three-dimensional geological mesh model of the target work area. Let the mesh region of this three-dimensional geological mesh model be... The position of any discrete node in the grid is determined by the grid nodes. express.

[0017] In this step, the computer equipment, based on the input multi-scale geological data, at each grid node... The calculation and synthesis of a second-order symmetric positive definite tensor, namely the Riemannian metric tensor field, is performed. Multi-scale geological data includes: seismic attribute data reflecting large-scale tectonic features, well logging inversion data reflecting small- and medium-scale petrophysical features, and geomechanical numerical simulation data reflecting geostress states. Riemannian metric tensor field. Used to define the non-Euclidean distance between two points in space, which characterizes the energy cost of a crack extending in a specific direction within that region.

[0018] Step S200: Optimization of crack nucleation points based on strain energy density The computer equipment calculates the elastic strain energy density of all grid nodes based on the stress tensor field and strain tensor field obtained from geomechanical simulation.

[0019] The computer equipment, based on a preset nucleation energy threshold and structural integrity coefficient, performs calculations within the grid region. Nodes with high strain energy density were selected as the initial seed set for crack growth. Initial seed point set Each seed point in It is assigned an initial growth potential that is positively correlated with its local excess strain energy.

[0020] Step S300: Synchronous Wavefront Competitive Growth Based on Generalized Voronoi Diagram Computer devices with an initial set of seed points Using all seed points as the starting source, the anisotropic fast traversal algorithm is used to solve the anisotropic functional equation.

[0021] In this process, the computer equipment simulates the energy wavefronts simultaneously excited at all seed points within a Riemannian metric tensor field. Propagation in a defined non-uniform Riemannian manifold space. The propagation velocity and direction of the wavefront are governed by the local metric tensor. Control. Through the synchronous advancement of multi-source wavefronts, the computational grid region is... The wavefront is divided into multiple generalized Voronoi units. Each generalized Voronoi unit contains a seed point, and the Riemann geodesic distance from any point within the unit to the seed point is less than the Riemann geodesic distance to any other seed point. Simultaneously, the initial growth potential energy is deducted from the wavefront's propagation energy based on the path integral during propagation; when the remaining energy is exhausted, the wavefront stops propagating.

[0022] Step S400: Wavefront Collision Detection and Dynamic Topology Connection The computer device detects contact events of wavefronts from different sources in real time during wavefront propagation. When wavefronts from different seed points meet at the boundary of a generalized Voronoi cell, the computer device calculates the wavefront gradient angle at the contact location and the local metric tensor features.

[0023] Computer equipment determines the contact type based on preset topology connection criteria: When the contact angle between wavefronts is less than the preset parallel threshold and the local rock brittleness meets the conditions, it is determined to be fusion, and the connection relationship between wavefronts is established. When the contact angle between the wavefronts is greater than the preset vertical threshold, it is determined to be truncation, and the wavefront with lower energy terminates on the path of the wavefront with higher energy, forming a branch structure. When the contact point is located in a high-stress barrier zone, it is determined to be repulsive, and no connection is established between the wavefronts.

[0024] Step S500: Geodesic backtracking and fracture geometry solidification The computer device performs a gradient descent backtracking operation based on the wavefront arrival time field generated in step S300 and the topological connectivity determined in step S400.

[0025] For each stopping boundary point or topological connection point of the wavefront, the computer device traces the minimum energy path, i.e., the geodesic, in the Riemannian manifold space along the inverse gradient direction of the wavefront arrival time field. The computer device connects the series of discrete points obtained by the tracing to generate a three-dimensional fracture skeleton surface, and assigns fracture pore size and permeability parameters according to the mesh properties along the path, ultimately generating a three-dimensional discrete fracture network model.

[0026] See attached document Figure 2 Before constructing the Riemannian metric tensor field, it is necessary to first establish a unified spatial reference system and complete the mapping and standardization of multi-source data. This step involves discretizing geological data with different resolutions and physical meanings into a unified grid region. middle.

[0027] The computer equipment first constructs a 3D corner grid or Cartesian grid for the target work area. This grid covers the entire geological modeling area, and each cell in the grid is considered a volume element, with its geometric center point defined as a grid node. For each grid node The system initializes a set of attribute channels to store scalar and tensor field data required for subsequent calculations.

[0028] For processing large-scale structural constraint data, the computer equipment reads the interpretation horizon data and fault polygon data. The computer equipment discretizes the interpretation horizon, calculates the plane normal vector and principal curvature (including the maximum and minimum principal curvature) at each grid node, and generates a structural curvature attribute field. Simultaneously, based on the fault polygon data, the computer equipment uses the Euclidean distance transformation algorithm to calculate the vertical distance from all grid nodes to the nearest fault plane, generating a fault distance field. The fault distance field is used to subsequently define the control effect of large-scale tectonic structures on fracture development intensity; that is, the closer to the fault, the lower the resistance to fracture propagation along the fault strike, and the more obvious the anisotropic characteristics.

[0029] For processing small- to medium-scale rock physical constraint data, computer equipment calculates Young's modulus using a rock physical volumetric model based on well logging data (such as P-wave velocity, S-wave velocity, and density logging) and pre-stack seismic inversion data. Compared to Poisson The three-dimensional distribution of the rock. Computer equipment uses the aforementioned elastic parameters to calculate the rock brittleness index. In one implementation, the rock brittleness index The rock brittleness index was obtained by weighted averaging of normalized Young's modulus and normalized Poisson's ratio, and further normalized to the (0,1] interval. The higher the value, the more easily the rock is broken, and the shorter the corresponding geodesic metric length in the subsequent metric tensor construction.

[0030] For processing geostress constraint data, the computer equipment reads the stress field data output from the geomechanical numerical simulation. Since geomechanical simulations typically use tetrahedral unstructured meshes, the computer equipment performs three-dimensional spatial interpolation operations to convert the stress tensor on the unstructured mesh into... Mapped to the grid area in this embodiment Grid nodes Above. Each mapped node stores the complete stress tensor components, including three normal stress components and three shear stress components. The computer further performs eigenvalue decomposition on the stress tensor of each node to extract the direction vector of the maximum principal stress. , intermediate principal stress direction vector and the direction vector of minimum principal stress And the corresponding principal stress scalar values.

[0031] To eliminate the influence of different physical dimensions on the construction of the metric tensor, the computer equipment processed the generated fault distance field. , structural curvature property field, rock brittleness index Furthermore, the principal stress magnitude field is dimensionlessly standardized. After the above processing, the mesh region... Each grid node in Each of them carries all the geometric and physical components required to construct the anisotropic Riemannian metric tensor, forming an initialized multi-physical property field.

[0032] The steps involved in constructing the macrostructure-guided tensor aim to transform large-scale geological structural features into mathematical tensor representations, which are then used to impose macroscopic directional constraints on a Riemannian metric field. Macrostructure-guided tensor It is mainly determined by the geometric characteristics of faults and their associated fracture zones, and is also regulated by the curvature of strata folds.

[0033] Computer equipment is primarily based on fault distance fields. For each grid node Determine the fault influence region. The computer equipment sets a threshold for the width of the fault damage zone. For distances satisfying The nodes were determined to be within the fault control domain; for The node is determined to be in the background domain.

[0034] For a node located within a fault-controlled domain, the computer searches for the nearest fault mesh patch and extracts the unit normal vector of that discrete fault patch. The unit normal vector This represents the direction perpendicular to the fault plane, typically the direction in which a fracture must overcome the greatest resistance to penetrate the fault. Simultaneously, the tangential plane of the fault plane is spanned by two orthogonal vectors: the fault strike vector and the fault strike vector. and fault dip vector These two directions represent the dominant directions in which fractures extend along the fault fracture zone.

[0035] The computer device constructs an anisotropic weighting function based on distance attenuation. The function reaches its maximum value at the fault plane, and its value increases with the distance from the fault. The nonlinear decay of the weighting function decreases to zero. Using this weighting function, the computer device constructs the macroscopic construction steering tensor according to the following formula. ; in, It is a third-order unit tensor, representing an isotropic background metric; To construct anisotropic strength coefficients, which are used to adjust the strength of the constraint of the fault on the fracture orientation; This is the dyadic tensor (extraproduct) along the fault normal, which increases the metric eigenvalue along the normal direction; This represents the transpose operation of a vector.

[0036] Furthermore, the computer equipment incorporates a correction function for the constructed curvature. For high curvature regions (such as the anticline hinge zone), the computer equipment extracts the direction vector of the maximum principal curvature. In regions far from faults, computer equipment will use the tensor construction formula described above. Replace with perpendicular to The direction vector is such that the dominant orientation of the generated fracture is parallel to the axis of maximum curvature, which conforms to the development law of tension fractures in fold deformation. The final result is... It is a symmetric positive definite tensor that not only encodes the dominant orientation of fracture development, but also quantifies the degree of this dominance through the difference in eigenvalues.

[0037] The steps for constructing the microscopic rock mechanics constitutive tensor are based on the rock physical properties and geostress state of the grid nodes, constructing a microscopic rock mechanics constitutive tensor that reflects the ease of local medium fracture and stress orientation. This tensor physically describes the distribution of a rock's ability to resist the propagation of fractures in different directions under geological stress conditions.

[0038] The computer device first calls the stress tensor For each grid node The stress tensor at a given location is subjected to eigenvalue decomposition. The computer extracts three mutually orthogonal unit eigenvectors, which are the vectors representing the directions of the maximum principal compressive stress. , Vector of intermediate principal compressive stress and the direction vector of minimum principal compressive stress Simultaneously, the computer equipment acquires the corresponding three principal stress scalar values, which are denoted as follows: , and In this embodiment, the compressive stress is defined as a positive value, and the magnitude relationship satisfies... .

[0039] According to the propagation criterion in fracture mechanics, in a compressive stress field, the dominant propagation plane of a tensile fracture typically includes the direction of the maximum principal stress and is perpendicular to the direction of the minimum principal stress. This means that the fracture tip extends along the direction of the maximum principal stress. The energy required for forward propulsion is minimal, and along the direction of minimum principal stress... The energy resistance required for expansion (i.e., attempting to increase the fracture aperture or change the fracture surface normal) is the greatest.

[0040] Based on the aforementioned physical mechanism, computer equipment constructs the microscopic rock mechanics constitutive tensor using spectral decomposition. This tensor is composed of the superposition of resistance components in three orthogonal directions, and its mathematical expression is as follows: ; in, Represents a grid node; The values ​​are 1, 2, and 3, which correspond to the three principal stress directions, respectively. For the first Eigenvectors of the principal stress directions; Indicates along the eigenvector The directional dyadic tensor (outer product matrix) is used to define the directional skeleton of the tensor; For along the feature vector The Riemannian metric eigenvalue of the direction physically represents the cost of a crack extending per unit length along that direction; This represents the transpose operation of a vector.

[0041] To reflect the combined control of rock brittleness and stress anisotropy on propagation resistance, the computer equipment calculates the Riemannian metric eigenvalues ​​according to the following logic. : Computer equipment first introduced the rock brittleness index (Values ​​range from 0 to 1), serving as a basic inverse proportional factor. The higher the rock brittleness, the lower the foundation cost. Subsequently, the computer equipment anisotropically modulates the foundation cost using stress state.

[0042] Specifically, for the direction of maximum principal stress (i.e. (At that time), the computer equipment is set to the minimum cost. Its calculations are subject to brittleness control or minimal stress suppression: ; in, To prevent tiny positive numbers with a denominator of zero.

[0043] For the direction of minimum principal stress (i.e. (At that time), the computer equipment is set to the maximum cost. The computer device introduces a penalty factor based on the differential stress ratio, making... Greater than This stretches the edges in the Riemannian metric space. The distance in a direction forces the calculated geodesic path to avoid that direction.

[0044] For the direction of the intermediate principal stress (i.e. (Time), intermediate costs of setting up computer equipment Its value is between and The values ​​between them are usually determined by linear interpolation.

[0045] Through the above steps, the computer device fuses discrete scalar properties and vector stress fields into a continuously varying second-order tensor field, which is the microscopic rock mechanics constitutive tensor. It can precisely guide the growth of fractures at the microscale along the trajectory where the rock is weakest (highly brittle) and where stress resistance is least (direction of maximum principal stress).

[0046] The steps of Riemannian metric tensor field synthesis and eigenvalue decomposition aim to mathematically couple the macroscopic structure-guiding tensor obtained from the aforementioned steps with the microscopic rock mechanics constitutive tensor to generate the final Riemannian metric tensor field that controls fracture network growth. And perform eigenvalue decomposition on it to obtain the geometric parameters required for anisotropic propagation.

[0047] The computer equipment first performs a tensor weighted composition operation. At each grid node... At this point, the computer equipment utilizes the principle of linear weighted superposition to fuse tensor information of different scales. The synthesis formula is as follows: ; in, The synthesized Riemannian metric tensor field is a third-order symmetric positive definite matrix, which physically represents the state of the grid nodes. The fracture growth resistance field integrates the tectonic setting and rock mechanical properties; This is a macro-scale weighting coefficient used to adjust the strength of the control of large-scale fault structures on fracture distribution; This is a microscale weighting coefficient used to adjust the degree to which local lithological brittleness and stress state control the detailed morphology of fractures. To construct a macroscopic guiding tensor; It is the constitutive tensor of microscopic rock mechanics.

[0048] Through the above synthesis Near the fault, the direction of fracture growth is mainly dominated by the macroscopic structural guidance tensor, showing a strong consistency in structural orientation; while in areas far from the fault, it is mainly dominated by the microscopic rock mechanics constitutive tensor, and the direction of fracture growth is more responsive to local stress deflection and lithological changes.

[0049] After synthesis is complete, in order to support the subsequent anisotropic fast traversal algorithm, the computer device performs a Riemannian metric tensor field analysis at each node. Perform eigenvalue decomposition. The decomposition process aims to extract the principal axis directions and axis lengths of the local metric ellipsoid. The computer solves the characteristic equation to obtain three eigenvalues. and the corresponding three orthogonal unit eigenvectors .

[0050] In the physical definition of this embodiment, eigenvalues This represents the eigenvectors along the corresponding orthogonal unit vectors. The Riemannian cost or growth resistance accumulated per unit Euclidean distance traveled in a direction. A larger eigenvalue indicates greater resistance in that direction and slower wavefront propagation.

[0051] Based on the above decomposition results, the computer device further calculates the dual diffusion tensor. This tensor is the inverse tensor of the metric tensor, i.e. In numerical implementation, computer devices directly construct the eigenvalues ​​using their reciprocals. At this point, the reciprocal of the smallest eigenvalue ( The maximum virtual propagation velocity is represented by , and the direction of its corresponding eigenvector is the optimal growth direction of the crack; while the reciprocal of the maximum eigenvalue ( ) represents the local minimum virtual propagation speed, corresponding to the most difficult direction for crack growth.

[0052] Finally, the computer equipment stores the calculated set of eigenvectors and the corresponding set of velocity eigenvalues ​​(i.e., the reciprocals of the eigenvalues) in the grid nodes, which serve as input parameters for subsequent solutions to the anisotropic equations. Geometrically, this set of parameters defines an ellipsoid centered at the node, with its major axis pointing in the direction where fractures are prone to develop and its minor axis pointing in the direction where they are suppressed, thereby precisely controlling the instantaneous propagation morphology of the wavefront at that location.

[0053] The purpose of the elastic strain energy density calculation step is to quantify the deformation energy accumulated at various points within the geological body, which is the physical essence driving fracture nucleation and propagation. The computer equipment uses the full-field stress and strain data obtained through geomechanical simulation in the preceding steps to calculate the elastic strain energy density of each grid node.

[0054] The computer device first traverses the three-dimensional mesh space. Each grid node in For the current node, the computer device reads the stress tensor corresponding to that point from memory. and strain tensor Among them, stress tensor It includes the three normal stress components and three shear stress components acting on that point, and the strain tensor. This includes three linear strain components and three engineering shear strain components (or tensor shear strain components, depending on the definition of the input data; in this embodiment, they are uniformly converted to tensor form) occurring at that point. Based on the assumption of linear elasticity, the computer device calculates the elastic potential energy stored per unit volume at that node by calculating the dot product of the stress tensor and the strain tensor. Elastic strain energy density The calculation formula is as follows: ; in, Represents grid nodes The elastic strain energy density at a point, whose physical unit is energy / volume (e.g., joules / cubic meter), is a non-negative scalar field; The trace operation represents the sum of the elements on the main diagonal of a matrix. For grid nodes Stress tensor at the point; For grid nodes The strain tensor at that point.

[0055] The above formula, in component form, is equivalent to calculating half the sum of the products of all stress components and their corresponding strain components. This calculation process compresses the complex tensor field state into a single elastic strain energy density. The elastic strain energy density It directly reflects the degree of energy accumulation in rocks under the influence of tectonic movements.

[0056] To facilitate subsequent numerical processing and threshold determination, the computer equipment normalizes the elastic strain energy density field after completing the full-field calculation. The computer equipment then calculates the maximum value of the energy density across the entire field. and minimum value and each node's Mapping to a dimensionless interval. After this step, each node in the three-dimensional geological grid is assigned an energy attribute value. The higher the value, the more unstable the rock is in a high-energy state, and the greater the possibility of fracturing and releasing energy, thus providing a physical basis for the spatial location of fracture seeds.

[0057] The purpose of the fracture nucleation point selection and screening step is to discretize the continuously distributed energy field into specific fracture growth initiation locations, i.e., seed points. Based on a physical driving mechanism, the computer equipment identifies the set of discrete nodes most prone to fracture in a three-dimensional mesh space.

[0058] The computer equipment first constructs a comprehensive nucleation potential energy index field. The index is a scalar field used to comprehensively characterize grid nodes. The probability of fracture nucleation occurring at a given location. Computer equipment combined with elastic strain energy density. and rock brittleness index Weighted calculations were performed. Physically, high strain energy density provides the driving force for fracture, while a high brittleness index means that the rock has a low plastic yield threshold; both factors contribute to the formation of fractures.

[0059] The computer equipment calculates the nucleation potential energy index for each grid node using the following formula: ; in, The elastic strain energy density of the node; and These represent the maximum and minimum values ​​of the total energy density, respectively. The energy sensitivity index is a constant greater than 1, used to nonlinearly stretch the differences in high-energy regions and highlight high-energy anomalies. The rock brittleness index; This represents the fault distance field. To construct the control weight coefficients, It is a very small positive number. The function of this term is to assign additional nucleation weights to the region near the fault, simulating the development characteristics of fault-associated fractures.

[0060] The comprehensive nucleation potential energy index field was calculated. Then, the computer equipment performs a dual screening mechanism to determine the final set of initial seed points. .

[0061] The first screening step is absolute threshold determination. The computer device sets a nucleation threshold. (For example, take) (85th percentile value). The computer device traverses all grid nodes and will satisfy... Nodes that meet the criteria are marked as background nodes and excluded from the candidate list; The nodes are retained as candidate seed points.

[0062] The second screening step is spatial repulsion (non-maximum suppression) judgment. To avoid generating overly dense clusters of seed points in local high-value regions (which would cause simulated fractures to overlap at the same point), the computer equipment introduces a minimum nucleation spacing parameter. The computer equipment processes all candidate seed points according to... Sort the data from largest to smallest. The computer then selects the node with the highest potential energy in the sequence as the seed node. and add to the initial seed point set. Each seed point is established Computer devices are those located in the remaining candidate list, and all those located in the list are selected. Center of the sphere, radius is Other candidate points within the spatial neighborhood are eliminated.

[0063] The computer device repeats the above sorting and elimination process until the candidate list is empty. The final output is the initial seed point set. Include These nodes are discrete grid node coordinates. They are located in positions that are optimal in terms of both energy and brittleness, and they also maintain a spatial distribution density that conforms to geological laws, serving as the source for subsequent wavefront competitive growth.

[0064] The purpose of the initial growth potential energy allocation step for fractures is to convert the local high strain energy at the fracture seed point into the total driving force for fracture propagation, i.e., to determine the maximum theoretical area or length that each fracture surface can generate. The computer equipment provides the initial seed point set. Each seed point in Calculate a scalar value, namely the initial growth potential. This potential energy value defines the maximum Riemann geodesic distance that a wavefront originating from this seed point can accumulate before it stops propagating.

[0065] The computer equipment first determines the critical breaking energy density threshold at the location of each seed point. This threshold characterizes the ultimate energy that a rock can withstand while maintaining an elastic state; exceeding this value results in irreversible brittle fracture. The computer calculates this critical threshold using the maximum normal stress criterion, based on the rock's tensile strength and Young's modulus at the seed point.

[0066] Subsequently, the computer equipment calculates the excess effective strain energy at the seed point and converts it into initial growth potential energy for crack propagation. The calculation formula is as follows: ; in, Seed point Elastic strain energy density at the location; Seed point The critical fracture energy density threshold at the location; The volume of the grid cell is used to convert energy density into total energy. The energy conversion efficiency coefficient is a pre-defined dimensionless constant. This coefficient takes into account non-surface energy losses such as heat dissipation and plastic deformation dissipation during rock fracturing, indicating that only a specific proportion of excess energy is effectively used to overcome the resistance defined by the Riemannian metric field to create the fracture surface; This is a function that maximizes the value of a function to ensure that when the current energy is less than the critical threshold, the growth potential energy is zero, meaning that no expansion occurs.

[0067] The computer equipment will calculate The numerical values ​​are stored in the data structure for each seed point. In subsequent wavefront competitive growth steps, this... This will serve as the cutoff threshold for wavefront propagation. Each step the wavefront takes on the Riemannian manifold consumes a certain amount of potential energy according to the local metric tensor; when the accumulated consumption exceeds... At this point, the growth process of the fracture immediately terminates. Through this mechanism, the size of the generated fracture is positively correlated with the degree of local energy accumulation in the geological body; that is, large fractures are generated in high-energy areas and small fractures are generated in low-energy areas, thus achieving physical constraints on the fracture size distribution.

[0068] See attached document Figure 3 The core of constructing the anisotropic equation and wavefront propagation lies in establishing a mathematical model describing the evolution of the fracture front in geological space. Instead of directly simulating the discrete geometric topology of the fracture, computer equipment uses the concept of level sets or rapid travel to transform the fracture growth process into a scalar field, namely the arrival time field. Solving boundary value problems.

[0069] In this physical model, each crack seed point is considered a wave source. The wavefront diffuses outward from the wave source, and its diffusion velocity and direction are strictly controlled by the diffusion tensor field (or velocity tensor field). To describe this complex anisotropic propagation behavior, the computer device establishes the static Hamilton-Jacobi equations, i.e., the anisotropic equations.

[0070] Computer equipment in grid area The internally defined anisotropic equation is as follows: ; in, This represents the arrival time field (or Riemann geodesic distance field). Physically, it represents the time-to-grid distance of the wavefront from the nearest seed point to the grid nodes. The minimum geological cost or generalized time consumed in the cumulative process; This represents the spatial gradient vector of the arrival time field. Geometrically, the direction of this vector is perpendicular to the local tangent plane (i.e., the normal direction) of the wavefront, and its magnitude reflects the slowness of the wavefront advance. This represents the transpose operation of a vector; The diffusion tensor is represented by the Riemannian metric tensor field. Inverse equation (i.e.) is obtained by reversing the equation. Physically, It represents the conductivity of the local medium for fracture propagation or a virtual velocity ellipsoid. Its principal axis direction determines the dominant direction of fracture growth, and the magnitude of the eigenvalue determines the growth rate in that direction.

[0071] The above formula reveals the non-Euclidean characteristics of wavefront propagation: in anisotropic geological media, the energy transmission direction of the wave (i.e., the group velocity direction, which is the ray path of the actual growth of the fracture) and the normal direction of the wavefront (i.e., the phase velocity direction, the gradient direction) usually do not coincide. The distribution of eigenvalues ​​determines the degree of this deviation: in the fault strike or in the direction of maximum principal stress, Having large eigenvalues ​​makes The projection component in this direction must be small to satisfy the equation (i.e. (Slow growth) leads to rapid forward propagation of the wavefront; conversely, in the direction perpendicular to the fault or the minimum principal stress, the wavefront is highly impeded and its propagation is significantly delayed.

[0072] To solve this partial differential equation, the computer device sets boundary conditions: for the seed point set, their arrival time... The arrival times of the background nodes are initialized to 0; the arrival times of the remaining background nodes are initialized to infinity. With this setup, the solution to the equation... This will be a monotonically increasing function, whose value precisely records the order in which the wavefront sweeps through each point in space. The final solution obtained... The field actually constructs a generalized Voronoi graph space based on Riemannian metrics, where each node... of The value is the minimum Riemannian energy cost from that point to its source point.

[0073] The purpose of the multi-source wavefront synchronous propagation and competitive growth algorithm is to simulate the process in a numerical grid where multiple fracture sources simultaneously release energy and compete for survival space. The computer equipment employs an anisotropic fast traverse algorithm based on a heap data structure to solve the aforementioned anisotropic equation, thereby obtaining the overall arrival time distribution and fracture attribution relationships. The computer equipment first initializes three core data fields: State marker field: used to record the computational state of grid nodes, divided into known points (values ​​are fixed), wavefront points (candidate values ​​are being calculated), and unknown points (not yet touched by the wavefront).

[0074] Arrival Time Field Initialized to infinity, only the initial seed point set. The nodes in the array are initialized to 0.

[0075] Source Index Field : Used to record which fracture seed point captured the node. During initialization, the seed point... place The remaining points are marked as invalid values.

[0076] The computer constructs a minimum priority queue and pushes all seed points into the queue. During algorithm iteration, the queue remains ordered by arrival time. The values ​​are arranged in ascending order. This simulates the shortest time principle or Fermat's principle in physics, which states that a wavefront always propagates preferentially to the position of least resistance (shortest time).

[0077] Subsequently, the computer device enters a loop-progression phase, repeating the following operations until the queue is empty: Node freezing and selection: Computer devices are popped from the top of the priority queue. The node with the smallest value is designated as the current master node. The state of this node is marked as a known point, meaning that the shortest anisotropic geodesic path from the source point to this point has been determined, and its arrival time is... and its source Permanently locked.

[0078] Neighborhood Update and Anisotropic Computation: Computer Device Traversal Neighbor nodes in three-dimensional space .like If the point is not known, the computer attempts to update its arrival time. At this point, the computer solves the equation numerically using a discretized difference scheme. To reflect anisotropy, the computer does not directly use Euclidean distance, but instead calls the local diffusion tensor at that location. .

[0079] The computer equipment uses the upwind difference scheme to solve for the candidate arrival times. The equation is a quadratic equation in one variable. This equation reflects the constraint that the wavefront must propagate along a characteristic direction. If the calculated... Less than the current record value of the neighbor node Then an update will occur: Update time value: ; Update ownership: This means that neighboring nodes inherit the fracture source ID from the master node. This step is crucial for generating the generalized Voronoi diagram, ensuring that the space is divided into regions controlled by different seed points; Update queue state: If If a point is initially unknown, add it to a priority queue and mark it as a wavefront point; if it is already in the queue, adjust its position in the heap to maintain order.

[0080] Energy depletion check (growth termination): During an attempt to update Previously, the computer equipment performed an energy check. The computer equipment read the initial growth potential energy of the source to which the wavefront belongs. If the calculated candidate arrival time This exceeds the initial growth potential of the source, that is: .

[0081] The computer then determines that the crack growth momentum in that branch has been exhausted and abandons the attempt. The update mechanism ensures that the final length of each fissure is strictly controlled by the local energy supply during its nucleation, thus preventing unlimited growth.

[0082] Through the aforementioned multi-source synchronous propagation, the wavefront advances rapidly in the low-resistance region (high brittleness, direction of maximum principal stress) and slowly in the high-resistance region. When two wavefronts from different sources meet, due to the irreversible state of known points, the later-arriving wavefront cannot cover the nodes occupied by the earlier-arriving wavefront. This natural collision and stopping mechanism mathematically automatically constructs the boundary of the weighted anisotropic Voronoi diagram, which corresponds to the mutual truncation and connectivity relationships in complex fracture networks.

[0083] The purpose of this step, Generalized Voronoi Partitioning and Fracture Geometry Extraction, is to analyze the final state of multi-source wavefront competitive growth, transforming the continuous arrival time field and source index field into a discrete fracture network model with clear geological significance. The computer equipment utilizes the source index field generated in the preceding steps. The three-dimensional geological space is divided into several non-overlapping volume units, which constitute the generalized anisotropic Voronoi partition.

[0084] Computer device definition seed points Controlled partitions A set of spatial points that satisfies the following conditions: ; in, This represents any grid node in a three-dimensional geological grid; This represents the mesh area of ​​the entire 3D geological model; Mathematical symbols indicating that an element belongs to a set; This indicates a condition qualifier, meaning that the following conditions must be met; Represents grid nodes The source index field at that location.

[0085] In a physical sense, partitioning Represents the origin of the first The volume of rock covered or conquered by the fissure at each seed point during growth. Unlike traditional Voronoi diagrams based on Euclidean distance, this embodiment generates... It exhibits obvious anisotropic characteristics: Shape distortion: subject to Riemannian metric tensor field Control The shape is not a regular polyhedron, but rather a flat or elongated irregular body stretched along the direction of the dominant cleavage or the direction of the maximum principal stress. This accurately reflects the geological fact that fractures are more likely to extend in areas with well-developed bedding or stress concentration.

[0086] Boundedness: Subject to initial growth potential Restrictions, partitions It doesn't necessarily fill the entire computational space. In regions far from the seed point or where energy is depleted, there exists a background region not covered by any partition. These boundaries are called free boundaries and represent the natural stopping points of the fissure tips.

[0087] Competition boundary: When two adjacent partitions and ( When they are in close contact, their contact surfaces form a competing boundary. This boundary corresponds to the arrival time field. The ridge line, or wavefront collision location, represents the point where two fractures connect or truncate each other in fracture network modeling. After completing the spatial partitioning, the computer equipment executes either a skeleton extraction algorithm or a mid-plane extraction algorithm to reconstruct the two-dimensional surface geometry of the crack from the volume data. For each partition... The computer device identifies the main axis direction (i.e. the longest extension direction) of its geometry and extracts the central patch in that direction as the discrete geometric representation of the crack.

[0088] The computer device specifically performs the following topology analysis operations: Independent fracture identification: if partitioned If the boundary of a crack consists entirely of free boundaries and does not contact other partitions, it is considered an isolated crack.

[0089] Crack intersection identification: If partitioning and A common contact surface exists, and the computer device detects the arrival time gradient directions on both sides of the contact surface. If the angle between the gradient directions on both sides exceeds a preset angle threshold (e.g., an obtuse angle), a T-shaped or X-shaped intersection is determined to have occurred. The computer device marks this contact surface as a connecting channel for fluid flow.

[0090] Through the above process, this invention freezes the complex physical-driven growth process into a final static geometric model. This model not only includes the location, size, and attitude of the fractures (determined by the anisotropic morphology of the partitions), but also inherently contains the topological connections between the fractures. Without the need for subsequent complex artificial geometric splicing, it directly generates a three-dimensional discrete fracture network (DFN) that conforms to the geomechanical mechanism.

[0091] The purpose of studying the energy dissipation and boundary formation steps of fracture growth is to introduce physical constraints in real time during wavefront propagation, simulating the mechanical phenomenon where the fracture tip naturally stops developing due to the decay of driving energy. Computer equipment precisely defines the geometrical free boundary of the fracture by comparing the cumulative generalized geological cost with the initially allocated energy budget.

[0092] When a computer device executes a multi-source wavefront competitive growth algorithm, for wavesfronts belonging to the first... seed points Arbitrary wavefronts arriving at grid nodes Perform the following energy dissipation determination: ; in, For grid nodes Arrival time field at the location; The initial growth potential energy of the fracture source represents the total work limit required for the fracture to overcome the resistance of the surrounding rock and expand.

[0093] Based on the above determination results, the computer equipment classifies the growth state of the crack into three cases and marks them in the mesh: Active growth zone: meets the requirements Within this region, the fracture has enough remaining energy to continue extending forward. The wavefront continues to advance towards neighboring nodes until it encounters the wavefront of another fracture (forming a competing boundary) or its energy is exhausted.

[0094] Natural termination boundary: satisfies The isosurface location is identified by the computer as the tip line of the fracture. Here, the energy driving the fracture to open is just balanced with the resistance to overcoming the rock's fracture toughness, and the fracture stops growing. This results in the formation of partitions. Instead of being an infinite region that fills the entire space, it is confined to a finite anisotropic geometry, thus avoiding the defect of traditional Voronoi diagrams that extend infinitely in non-competitive directions.

[0095] Background area: meets the requirements The area is marked as an intact, unbroken rock matrix.

[0096] It is worth noting that, due to The growth rate is affected by the diffusion tensor (i.e., Riemannian metric tensor field) The inverse of diffusion tensor control means that energy dissipation exhibits significant directionality in space. This is particularly evident in directions with large eigenvalues ​​of the diffusion tensor (such as the direction of maximum principal stress or bedding planes). The slower the increase in geometric distance, the less energy is consumed per unit distance, so the crack can extend further in this direction; conversely, in the direction with smaller eigenvalues ​​(such as high-stress retardation regions). As the distance increases rapidly, energy is depleted quickly, resulting in a short crack propagation distance.

[0097] Through this mechanism, the fractures generated by this invention automatically exhibit geometric characteristics adapted to the geological environment: nearly circular in isotropic media; elliptical or elongated in media with well-developed bedding, with the major axis along the bedding direction; and twisted, irregular curved surfaces in complex stress fields. The resulting fracture network model is not only accurate in topological connections but also strictly adheres to the laws of energy conservation and fracture mechanics in terms of the scale (aspect ratio, area) of individual fractures.

[0098] The purpose of wavefront collision interface identification and geometric localization steps is to accurately capture the geometric trajectories of fracture wavefronts from different nucleation points in contact with each other within a discretized three-dimensional grid space. These trajectories physically correspond to tangent lines or connectivity surfaces in the fracture network, determining the final connectivity topology of the fracture network.

[0099] The computer device iterates through each grid node and its neighboring nodes in the computational domain, based on the source index field generated in the previous steps. and arrival time field Identify wavefront collision interfaces.

[0100] The computer device first defines and searches for adjacent node pairs that satisfy the index jump condition. For any two spatially adjacent grid nodes... and If the following logical conditions are met, the computer device determines that there is a potential collision interface element between these two nodes: ; in, and Representing grid nodes respectively and The associated fracture source index ID; This indicates a background area not occupied by any wavefront; The AND operator means that all conditions within the parentheses must be met simultaneously for the condition to be true.

[0101] This condition indicates that, and They were captured by two different fracture sources, both of which were inside actively growing or fully grown fractures, rather than at the free boundary of the fracture tip.

[0102] To further determine the geometric properties and physical validity of the collision interface, the computer device calculates the arrival time gradient vectors at the aforementioned node pairs, i.e. and These two vectors indicate the local propagation direction of the wavefront at their respective locations. The computer device quantifies the severity of the collision by calculating the cosine of the angle between these two gradient vectors. ; in, Represents grid nodes Arrival time at the location The spatial gradient vector; Represents grid nodes Arrival time at the location The spatial gradient vector; Represents the spatial gradient vector The Euclidean norm; Represents the spatial gradient vector The Euclidean norm.

[0103] like The value is less than the preset negative threshold (i.e., the collision angle). A near 180-degree angle indicates that two wavefronts are propagating towards each other and meeting head-on. This corresponds to geological thrust cut-off, typically forming a fracture junction zone with high conductivity. Computer equipment labels such interfaces as ridges or skeletal boundaries.

[0104] Computer equipment utilizes geometric stitching algorithms or dual mesh techniques to... and The spatial positions between points are used to construct triangular facets, connecting discrete collision point pairs into a continuous curved surface. This surface constitutes the cell boundary of the generalized Voronoi diagram. Unlike the planar boundary of a traditional Voronoi diagram, the collision interface generated in this embodiment is a curved surface controlled by Riemannian metrics, and its normal direction is determined by the difference in propagation velocity between the wavefronts on both sides.

[0105] Through the above steps, the computer equipment not only identified the contact positions between fractures, but also eliminated artifacts through gradient analysis, ensuring that the identified interfaces are the physical intersections of wavefront energy transfer, thus providing a geometric basis for subsequent analysis of T-type, X-type, or Y-type topological relationships between fractures.

[0106] The purpose of the fracture collision geometric feature classification and topological evolution criterion steps is to further determine the fracture topology type generated by the collision from a geomechanical perspective, based on the identification of the wavefront collision interface (i.e., the cavity boundary of the generalized Voronoi diagram). The computer equipment no longer simply treats the boundary as a simple geometric dividing line, but instead, by introducing anisotropic metric parameters, it analyzes it as a fracture tangent zone with physical properties (such as conductivity and closure state), thereby determining whether the final generated discrete fracture network forms a connected, truncated, or traversing structure.

[0107] The computer device first determines the local geometric normals of the collision interface. For the mesh nodes... The wavefront collision occurring at a certain point, the normal vector of its interface. Defined by the difference vector of the arrival time gradients from both sides: ; in, and These are the two sides of the collision interface (belonging to the source respectively). Heyuan The arrival time-space gradient vector of the normal vector. It indicates the spatial orientation of the fracture tangent surface, that is, the direction perpendicular to the collision contact surface.

[0108] Subsequently, the computer equipment establishes a stress compatibility criterion based on anisotropy metrics to determine the topological connection form at the intersection. The computer equipment utilizes a Riemannian metric tensor field. Calculate the positive hysteresis coefficient This coefficient characterizes the degree to which the geological medium hinders fracture penetration in the direction perpendicular to the collision interface, and its calculation formula is as follows: ; In this formula, Physically, it represents the anisotropic geological cost of the local medium. The larger the value, the greater the geological resistance (such as high pressure stress or hard rock layers) in the direction perpendicular to the collision interface (i.e. the direction in which the fracture attempts to continue to expand through the interface), and the more difficult it is for the fracture to expand penetratingly.

[0109] Based on the calculated positive hysteresis coefficient The collision angle obtained in the previous step Computer devices execute dynamic topology evolution logic. When Less than the preset penetration threshold And the collision angle Within the effective shear range (e.g.) When the computer determines that the two fracture wavefronts have merged through a continuous flow at this point, the computer determines that the two fracture wavefronts have merged through a continuous flow. In the construction of the discrete fracture network model, the computer divides the two corresponding partitions... and The skeleton facets are geometrically stitched at this boundary to generate a connected edge, and the edge is marked as a high-conductivity channel for fluid flow, forming a Y-type or X-type connection.

[0110] Conversely, when Greater than the penetration threshold Or the collision is mainly manifested as a head-on collision ( When the system is in a high-impedance state, the computer equipment determines that the later-arriving wavefront will be blocked and intercepted by the earlier-arriving wavefront (or geological structure surface). In this situation, the computer equipment identifies the arrival time. On the larger side (the weaker energy side), its corresponding fracture skeleton terminates at this boundary, forming a T-node. In this case, the collision interface itself does not act as a fluid channel, or is marked as a low-permeability closed interface. Furthermore, in special regions with high brittleness and extreme anisotropy, if... The direction is exactly the same as The minimum eigenvalue direction (i.e., the minimum principal stress direction or the maximum permeability direction) is aligned, and even in the event of a collision, the computer equipment allows the fracture skeleton to cross each other, forming an X-shaped intersection.

[0111] By introducing the Riemannian metric tensor field By participating in collision determination, this invention overcomes the limitation that traditional geometric Voronoi diagrams can only generate straight or random boundaries, ensuring that the fracture intersection relationship (topology) and the stress field and medium properties underground remain physically consistent, thereby improving the realism of the fracture network model in subsequent fluid numerical simulations.

[0112] The purpose of the fracture network topology connection rule execution and entityization step is to transform the abstract collision events identified in the previous steps based on the stress criterion into concrete geometric entities and topological data structures in the Discrete Fracture Network (DFN) model. The computer device completes the construction of fracture network entities from wavefront contact by performing geometric Boolean operations and attribute mapping.

[0113] The computer device traverses the set of all interfaces marked as valid collisions. For each collision event, based on the topology type (connected, cut, or traverse), it performs a specific geometric reconstruction operation. When the determination result is a connectivity merge (X-type or Y-type), the computer device performs a mesh stitching operation. Specifically, the computer device stitches two adjacent partitions... and The internally generated fracture skeleton facets extend to the common collision interface (i.e., the ridge location). Along the contact line, the computer device forcibly aligns the vertex coordinates of the two facets, merging repeating nodes to form a physically continuous common intersection line. This common intersection line is defined as a highly conductive channel for fluid exchange, with an equivalent hydraulic gap width... It is usually set as a weighted function of the width of the fracture surface on both sides to simulate the expansion effect of the fracture zone at the intersection.

[0114] When the determination result is a blocking cutoff (T-shaped), the computer device performs a mesh clipping operation. The computer device identifies the cutoff side (i.e., the time of arrival). The larger, weaker energy side of the fracture is geometrically trimmed at the collision interface to prevent it from crossing the interface; while the truncation side (i.e., the arrival time) The smaller, higher-energy side of the fracture maintains its geometric integrity, serving as the termination boundary of the truncated fracture. This operation geometrically creates an asymmetric T-junction structure. At such intersections, the computer device does not establish a direct, fully conductive connection between the fluid nodes, but instead introduces a small interfacial conductivity multiplier (e.g., 0.01 to 0.1). It is used to simulate the fluid flow resistance caused by stress concentration, closure, or non-through properties at T-junctions.

[0115] After connecting and trimming the geometric entities, the computer equipment synchronously updates the graph data structure of the fracture network and assigns the fracture width attribute to the generated fracture surface entities using the diffusion tensor. For any point on the fracture surface... Its hydraulic gap width It is determined by the anisotropic diffusion coefficient at that point in the direction perpendicular to the fracture surface: ; in, This is the unit normal vector of the fracture surface at that point; The diffusion tensor at that location (i.e., the Riemannian metric tensor field) The inverse of the stress reflects the local stress state and brittle characteristics of the rock; This is a dimensionless scaling factor used to calibrate the model to match macroscopic geological statistics. This is an exponential parameter used to control anisotropic sensitivity. Through this formula, the present invention ensures that the crack has a larger aperture in the direction where it is prone to expansion (large diffusion coefficient), achieving integrated construction of geometric modeling and physical properties.

[0116] The purpose of the geodesic backtracking path integral and fracture geometry solidification steps is to utilize the gradient information of the arrival time field to accurately reconstruct the bending geometry of the fracture surface in anisotropic media through a reverse tracing algorithm. The computer equipment processes discrete volumetric data (partitions) It is transformed into a continuous and smooth curved surface geometry, ensuring that the generated fracture surface strictly follows the shortest path principle of wavefront propagation (i.e., Fermat's principle).

[0117] The computer device first extracts the first seed points Controlled partitions The set of all boundary nodes, including competing boundary nodes (collisions with other fractures) and free boundary nodes (energy-depleted apexes). For each boundary node in this set... The computer equipment performs a geodesic backtracking operation to find its path to the seed point. The optimal energy path.

[0118] The computer equipment uses a gradient-based path integral algorithm to calculate the backtracking trajectory. In the Riemannian manifold space, from the boundary node... Return to seed point The shortest path (i.e., geodesic) tangent direction and arrival time field The gradient directions are collinear, but modulated by the diffusion tensor. The computer generates the backtracking path point sequence by solving the following ordinary differential equation: ; in, The vector representing the spatial location on the backtracking path is related to the path length. The function; This represents the unit tangent vector of the path at the current point, indicating the instantaneous direction of backtracking; This is the spatial gradient vector of the arrival time field at the current location, pointing in the direction of wavefront propagation (i.e., the direction of the fastest time increase). For position The diffusion tensor (i.e., the Riemannian metric tensor field) at that location (The reverse of the original), which causes the backtracking path to deviate from a simple Euclidean straight line and instead deflect towards areas of high permeability or low stress, thus forming a curved trajectory that conforms to geological characteristics.

[0119] Computer equipment uses high-order numerical integration methods (such as the fourth-order Runge-Kutta method) to discretize and solve the above equations, thereby obtaining a series of streamlines converging from the boundary to the seed point. These streamlines physically represent the historical trajectory of fracture growth.

[0120] Subsequently, the computer device performs geometric skinning and triangulation operations. Using all the generated backflow lines as a skeleton, the computer device constructs triangular meshes between adjacent streamlines. Specifically, the computer device connects the streamlines based on isochronous nodes, reconstructing the one-dimensional set of lines into a two-dimensional fracture surface.

[0121] During the meshing process, the computer automatically performs surface smoothing. Due to the discreteness of the original mesh data, local noise may exist in the gradient calculation, and the generated initial surface may contain non-physical micro-vibrations. The computer uses Laplace smoothing or spline interpolation techniques to optimize the coordinates of the internal mesh nodes while keeping the overall fracture orientation and topological connectivity (i.e., boundary positions) unchanged, so that the fracture surface exhibits a smooth shape that conforms to the bending characteristics of rock mechanics.

[0122] Through this step, the present invention transforms the abstract scalar field (arrival time field) into a concrete, visualized three-dimensional geometric model. The generated fracture surface not only connects the seed point to the boundary, but its curvature change directly reflects the heterogeneity and anisotropy of the subsurface medium (such as bypassing high-stress zones and deflecting along bedding planes), thus providing high-fidelity geometric boundary conditions for subsequent fluid numerical simulations.

Claims

1. A method for modeling network fractures based on multi-scale factor constraints, characterized in that, The method comprises the following steps: The following calculation steps are performed by the computer device for the grid nodes of the three-dimensional geological grid model: Based on the multi-scale geological data, a Riemann metric tensor field is synthesized at the grid nodes, and a non-Euclidean distance cost required for a unit length of a fracture to expand in an anisotropic medium is defined by using the Riemann metric tensor field; The elastic strain energy density of the grid nodes is calculated, and an initial seed point is screened and an initial growth potential for limiting the growth range of the fracture is assigned according to the elastic strain energy density; An anisotropic eikonal equation is solved from the initial seed point, a whole-field wavefront arrival time field is calculated by spatially accumulating the non-Euclidean distance cost, and a grid region is dynamically divided into generalized Voronoi cells controlled by different initial seed points; A common boundary between the generalized Voronoi cells is identified as a wavefront contact interface, and a wavefront arrival time field gradient feature on both sides of the wavefront contact interface is extracted for determining a contact type and establishing a fracture topological connection; From the wavefront contact interface where the fracture topological connection is established, a geodesic line path is tracked to the corresponding initial seed point in the direction of the negative gradient of the wavefront arrival time field, and a three-dimensional discrete fracture network model is generated.

2. The method of claim 1, wherein, Before the Riemann metric tensor field is synthesized, the multi-scale geological data including seismic attribute data, logging inversion data and geomechanical simulation data are preprocessed, and the specific processing steps include: Based on the seismic attribute data, a structural curvature attribute field and a fault distance field to the nearest fault plane at the grid nodes are calculated as structural constraint data; Based on the logging inversion data, Young's modulus and Poisson's ratio of the grid nodes are calculated to obtain a rock brittleness index as rock physical constraint data; An unstructured stress tensor of the geomechanical simulation data is mapped to the grid nodes and subjected to eigenvalue decomposition, and a principal stress direction vector and a principal stress magnitude are extracted as geostress constraint data.

3. The method of claim 1, wherein the method is characterized by: The step of synthesizing the Riemann metric tensor field includes constructing a macroscopic structural guide tensor: For the grid nodes located in a fault control domain, a unit normal vector of the nearest fault plane is extracted; A structural anisotropy weighting function based on fault distance field attenuation is constructed; A macroscopic structural guide tensor is constructed by linearly superimposing a third-order unit tensor and a dyadic tensor along the fault unit normal vector direction, and the weight of the dyadic tensor is determined by the structural anisotropy weighting function and a structural anisotropy strength coefficient.

4. The method of claim 3, wherein, The step of synthesizing the Riemann metric tensor field includes constructing a microscopic rock mechanics constitutive tensor: Three principal stress direction characteristic vectors and corresponding three principal stress scalar values at the grid nodes are obtained; According to the rock brittleness index and the three principal stress scalar values at the grid nodes, Riemann metric eigenvalues along the three principal stress direction characteristic vectors are calculated respectively; The Riemann metric eigenvalue along the maximum principal stress direction is negatively correlated with the rock brittleness index, and the Riemann metric eigenvalue along the minimum principal stress direction is set to be greater than the Riemann metric eigenvalue along the maximum principal stress direction by introducing a differential stress ratio penalty factor. The micro rock mechanics constitutive tensor is constructed by weighted superposition of dyadic tensors in three principal stress directions and corresponding eigenvalues of Riemann metric.

5. The method of claim 4, wherein, The step of synthesizing the Riemann metric tensor field specifically comprises: The macro tectonic guiding tensor and the micro rock mechanics constitutive tensor are fused by using a linear weighted superposition principle to generate the Riemann metric tensor field; An inverse operation is performed on the tensors in the Riemann metric tensor field to obtain a diffusion tensor field, which is used to represent an anisotropic ellipsoid feature of wave front propagation speed.

6. The method of claim 1, wherein, The step of screening the initial seed points according to the elastic strain energy density specifically comprises: According to the stress tensor field and the strain tensor field obtained by geomechanical simulation, half of the double-point product of the stress tensor and the strain tensor at the grid node is calculated as the elastic strain energy density of the grid node; An integrated nucleation potential index field is constructed by combining the elastic strain energy density and a pre-acquired rock brittleness index; A nucleation threshold and a minimum nucleation spacing parameter are set, and nodes in the grid region whose integrated nucleation potential index is greater than the nucleation threshold and whose mutual distance is greater than the minimum nucleation spacing parameter are screened to form the initial seed point set.

7. The method of claim 1, wherein, The step of assigning the initial growth potential for limiting the growth range of the fracture specifically comprises: A critical fracture energy density threshold of the seed point position is determined; A difference between the elastic strain energy density at the seed point and the critical fracture energy density threshold is calculated; The difference is multiplied by the grid cell volume and an energy conversion efficiency coefficient to obtain the initial growth potential, which is equal in value to the maximum non-Euclidean distance cost that can be accumulated by the wave front from the seed point.

8. The method of claim 5, wherein the method further comprises: The step of solving the eikonal equation and calculating the full-field wave front arrival time field specifically comprises: The eikonal equation is constructed according to the numerical relationship that the quadratic form product of the spatial gradient vector of the arrival time field and the diffusion tensor field is equal to one; The anisotropic fast marching algorithm is used to solve the eikonal equation, specifically comprising: A minimum priority queue containing all seed points is used to manage the wave front propagation order, and the node with the minimum arrival time is processed preferentially; When the wave front advances to the neighbor node and the candidate arrival time is calculated by solving the locally discretized anisotropic eikonal equation, it is judged whether the candidate arrival time exceeds the initial growth potential of the corresponding seed point, and if so, the advancing is stopped to form a natural termination boundary of the fracture.

9. The method of claim 1, wherein the method is characterized by: The step of judging the contact type and establishing the fracture topology connection specifically comprises: The cosine value of the included angle between the spatial gradient vectors of the arrival time field on both sides of the wave front contact interface is calculated; The normal vector of the wave front contact interface is calculated by calculating the difference vector of the gradients of the arrival time field on both sides of the contact position, and the positive forward blocking coefficient is calculated according to the normal vector and the local Riemann metric tensor field; When the cosine value is less than a preset negative threshold and the positive forward blocking coefficient is less than a preset penetration threshold, it is judged that the contact type is fusion, and the wave front contact interface is marked as a connected channel; When the positive forward blocking coefficient is greater than the penetration threshold, it is judged that the contact type is truncation, and the fracture skeleton on the side with the larger arrival time is truncated at the wave front contact interface.

10. The method of claim 1, wherein the method is characterized by: tracking geodesic paths from the corresponding initial seed points in the direction of the negative gradient of the wavefront arrival time field comprises: for each node on the boundary of each generalized Voronoi cell, generating a backtracking streamline converging to the corresponding seed point by solving an ordinary differential equation containing the inverse of the local Riemannian metric tensor and the gradient of the arrival time field; using the backtracking streamlines as a skeleton, constructing a triangular mesh between adjacent streamlines and generating a fracture surface; assigning a hydraulic aperture attribute to the fracture surface according to the diffusion coefficient in the direction normal to the surface at each point on the fracture surface.

Images