Geological reservoir fracture modeling grid cell determination method based on oscillator dynamic energy optimization and high-order interaction

Abstract

The invention provides a geological reservoir fracture modeling grid cell determination method based on oscillator dynamic energy optimization and high-order interaction, and relates to the field of geology. According to the method, a dynamic energy optimization framework of an Ising machine is introduced into geomechanical modeling through fusion of amplitude dynamics and a high-order interaction model, the defects that in a traditional method, the physical significance of an energy function is fuzzy, and grid size effect correction is insufficient are overcome, and based on the real-value energy function of a complex oscillator variable and gradient descent, the dynamic energy optimization framework of the Ising machine is introduced into geomechanical modeling. Global convergence of mechanical parameters of the fracture network is realized, and the problem of local extremum of a traditional phase model is avoided; the fracture rigidity parameter is converted into a third-order polynomial energy item, and nonlinear characteristics of normal / shear deformation of the fracture surface are accurately described; a three-cycle method and a dynamic energy threshold value are combined, the grid size is synchronously optimized in the position-scale-azimuth dimension, the characterization error of fracture surface density and mechanical parameter anisotropy is reduced, and a standard grid format compatible with commercial software can be directly output.

Description

Technical Field

[0001] This invention relates to the field of geology, and more specifically, to a method for determining the grid cells for modeling fractures in geological reservoirs based on oscillator dynamic energy optimization and higher-order interactions. Background Technology

[0002] The exploration and development of fractured reservoirs (such as shale gas and tight oil reservoirs) heavily relies on geomechanical modeling techniques. The core objective is to determine the reservoir's "engineering sweet spot" and fracturing strategies through stress field simulation and fracture prediction. However, these reservoirs typically develop multi-scale, multi-directional natural fracture networks, leading to significant size effects and anisotropy in rock mechanical parameters. Traditional mesh generation methods (such as uniform meshes or empirical thresholding methods) struggle to accurately characterize the complex mechanical behavior of fracture networks, specifically manifested in the following ways:

[0003] 1. Size effect error: Small-scale mesh elements cannot cover the crack characterization element (REV), resulting in drastic fluctuations in mechanical parameters;

[0004] 2. Anisotropic distortion: Single-direction meshing ignores the asymmetric influence of crack orientation on Young's modulus and Poisson's ratio;

[0005] 3. Low computational efficiency: Existing finite element methods require iterative adjustments to mesh size and orientation, and computational complexity increases exponentially with crack density.

[0006] The main technologies for solving the above problems include the three-cycle statistical method, the phase oscillator model, and crack stiffness embedding. The three-cycle statistical method calculates equivalent mechanical parameters through position-scale-orientation cycles to determine the minimum mesh size that meets the accuracy threshold. The phase oscillator model is based on the Kuramoto model to simulate the synchronous behavior of the crack network and optimizes the mesh distribution using an energy function. Crack stiffness embedding uses the normal / shear stiffness coefficient as a constraint term to improve the physical consistency of crack surface deformation.

[0007] While existing methods have partially alleviated the challenges of crack modeling, the following bottlenecks remain:

[0008] Traditional phase models (such as Kuramoto) ignore amplitude dynamics and cannot support higher-order fracture interactions, causing mesh optimization to get stuck in local extrema; the three-cycle method relies on static fracture parameters and does not combine dynamic energy gradients to adjust the mesh size in real time, making it difficult to adapt to the multi-scale characteristics of heterogeneous reservoirs.

[0009] In summary, current technologies struggle to balance accuracy, efficiency, and dynamic adaptability in modeling fractured reservoirs, necessitating a mesh generation method that integrates higher-order interactions, dynamic energy optimization, and hardware acceleration. Summary of the Invention

[0010] To address the needs mentioned in the background art, this invention provides a method for determining the grid cells for geological reservoir fracture modeling based on oscillator dynamic energy optimization and higher-order interaction. This method aims to solve the technical problems of low grid division accuracy and poor computational efficiency caused by size effects, anisotropy, and complex fracture networks in traditional fractured reservoir mechanical modeling.

[0011] A method for determining the grid cells for geological reservoir fracture modeling based on oscillator dynamic energy optimization and higher-order interactions, comprising the following steps:

[0012] Step 1: Input the laboratory triaxial experimental data and the well logging P / S wave transit time and density data into the parameter mapping model to obtain the isostatic Young's modulus and Poisson's ratio; at the same time, based on the mapping model, generate a static mechanical parameter database covering various lithologies and stress conditions.

[0013] Step 2: Using occurrence, density, and combination pattern data, construct a fracture attribute set. Under the multi-scale discrete network framework, use Monte Carlo sampling combined with lithological parameters in the database to generate a three-dimensional discrete fracture set. Then, convert the discrete fracture set into a mesh geometric unit that can be recognized by the finite element / discrete element modeling platform through a general geometric interface.

[0014] Step 3: Map each mesh element to a complex oscillator node and construct a comprehensive potential energy model to characterize the energy distribution between mesh elements; introduce a stiffness constraint term into the potential energy model to achieve a unified expression of the crack surface normal and shear nonlinear characteristics;

[0015] Step 4: Overlay fracture node attributes onto the global reservoir grid to form an attribute mapping layer; based on the gradient information of the potential energy model, initiate energy-driven optimization to refine the high-energy gradient region and coarsen the low-gradient region; through a multi-dimensional cyclic correction framework, sequentially verify and adjust the grid size to meet the characterization unit accuracy and anisotropic error threshold.

[0016] Step 5: Adjust the weights of the high-frequency crack orientation elements according to the orientation term in the potential energy model to balance the anisotropic effect; convert the oscillator phase information into a binary mesh of crack and matrix through a tag mapping mechanism, and output a mesh format file that conforms to commercial software standards.

[0017] Step Six: Import the binarized mesh into the simulation platform to simulate working conditions such as water injection and fracturing, and compare it with the measured in-situ stress data; through three sub-modules of error statistics, efficiency comparison and stability verification, quantify the improvement of model accuracy and computational performance; automatically generate a closed-loop evaluation report containing error distribution, efficiency improvement rate and optimization suggestions to provide a basis for decision-making.

[0018] Furthermore, the node attributes of the complex oscillator node are jointly constituted by the database parameters stiffness and crack orientation.

[0019] Furthermore, the integrated potential energy model includes self-energy, second-order coupling, higher-order interactions, and stiffness constraints.

[0020] The beneficial effects of this invention are as follows: By integrating amplitude dynamics and higher-order interaction models, this invention introduces the dynamic energy optimization framework of the Ising machine into geomechanical modeling, overcoming the shortcomings of traditional methods such as the ambiguity of the physical meaning of the energy function and insufficient correction of the mesh size effect. Based on the real-valued energy function of the complex oscillator variable and gradient descent, this invention achieves global convergence of the mechanical parameters of the fracture network, avoiding the local extremum problem of the traditional phase model; it transforms the fracture stiffness parameter into a third-order polynomial energy term, accurately characterizing the nonlinear features of the fracture surface normal / shear deformation; and by combining the three-cycle method and dynamic energy threshold, it simultaneously optimizes the mesh size in the location-scale-orientation dimensions, reducing the characterization error of the fracture surface density and the anisotropy of mechanical parameters, and can directly output a standard mesh format compatible with commercial software. Attached Figure Description

[0021] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the embodiments will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as a limitation on the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0022] Figure 1 A flowchart of the method of the present invention is shown. Detailed Implementation

[0023] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. It should be understood that the accompanying drawings in the present invention are for illustrative and descriptive purposes only and are not intended to limit the scope of protection of the present invention. Furthermore, it should be understood that the schematic drawings are not drawn to scale. The flowcharts used in this invention illustrate operations implemented according to some embodiments of the present invention. It should be understood that the operations in the flowcharts may not be implemented in sequence, and steps without logical contextual relationships may be reversed or implemented simultaneously. In addition, those skilled in the art, guided by the content of this invention, may add one or more other operations to the flowcharts, or remove one or more operations from the flowcharts.

[0024] Furthermore, the embodiments described herein are merely some, not all, of the embodiments of the invention. The components of the embodiments of the invention described and illustrated herein can typically be arranged and designed in various different configurations. Therefore, the following detailed description of the embodiments of the invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the invention without inventive effort are within the scope of protection of the invention.

[0025] It should be noted that the term "comprising" will be used in the embodiments of the present invention to indicate the presence of a feature subsequently declared, but does not preclude the addition of other features. It should also be noted that similar reference numerals and letters in the following figures indicate similar items; therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures. In the description of the present invention, it should also be noted that the terms "first," "second," "third," etc., are used only for distinguishing descriptions and should not be construed as indicating or implying relative importance.

[0026] The following is a detailed description of this case, in conjunction with the relevant accompanying drawings in the instruction manual.

[0027] This scheme proposes a method for determining the grid cells for geological reservoir fracture modeling based on oscillator dynamic energy optimization and higher-order interactions, specifically including the following steps:

[0028] Step 1: Calibrating Rock Mechanical Parameters

[0029] The static mechanical parameters of rocks (such as Young's modulus and Poisson's ratio) directly reflect their actual deformation capacity, while dynamic parameters indirectly characterize rock stiffness through the propagation speed of sound waves. Therefore, experimental calibration can, to a certain extent, ensure that the rock mechanical behavior in the numerical model is consistent with that of the real reservoir.

[0030] 1.1 Based on the above objectives, stress-strain curves of rocks under different confining pressures were first obtained through laboratory triaxial mechanical experiments, combined with the P-wave transit time Δt from well logging data. p Transverse wave time difference Δt s and density data ρ b Calculate the dynamic Young's modulus E d and dynamic Poisson ratio μ d .

[0031] Longitudinal wave time difference Δt p It represents the time required for a longitudinal acoustic wave (P-wave) to travel a unit distance (usually 1 foot) in rock, measured in microseconds per foot (μs / ft), and its response reflects the compressibility of the rock.

[0032] Transverse wave time difference Δts The time it takes for a transverse acoustic wave (S-wave) to travel a unit distance is represented by microseconds per foot (μs / ft), and it reflects the rock's shear resistance.

[0033] Density data ρ b This represents the rock bulk density measured by the logging instrument, expressed in grams per cubic centimeter (g / cm³). 3 ).

[0034] By converting the transverse wave time difference and longitudinal wave time difference into transverse wave velocity and longitudinal wave velocity respectively through unit conversion, the dynamic Young's modulus E can be obtained using the relationship between elastic wave velocity and density. d The mathematical expression for it is shown in equation (1):

[0035]

[0036] By using the ratio of shear wave velocity to longitudinal wave velocity, combined with density data ρ b This comprehensively reflects the rock's resistance to deformation;

[0037] Dynamic Poisson ratio μ d The mathematical expression for it is shown in equation (2):

[0038]

[0039] Dynamic Poisson ratio μ d The ratio of transverse strain to longitudinal strain in a rock under stress reflects the compressibility of the material.

[0040] Calculating dynamic Young's modulus and Poisson's ratio using well logging data is a key step connecting geophysical measurements and rock mechanics analysis, providing fundamental parameters for reservoir geological modeling, fracturing design, and other applications.

[0041] 1.2 Dynamic parameters need to be converted into static values ​​through experimental calibration to more accurately reflect the mechanical behavior of rocks. Therefore, based on the experimental calibration formula, the dynamic parameters (dynamic Young's modulus E) are converted into static values. d and dynamic Poisson ratio μ d The parameters are converted into static parameters (Young's modulus E and Poisson's ratio μ), and a database of static mechanical parameters of reservoir rocks is established, covering the stiffness and deformation characteristics of different lithologies. The mathematical expression of Young's modulus E is shown in equation (3):

[0042] E = k1E d +b1(3)

[0043] The mathematical expression for Poisson's ratio μ is shown in equation (4):

[0044] μ=k2μ d +b2(4)

[0045] In equations (3) and (4), k1, k2, μ1, and μ2 are experimental calibration coefficients.

[0046] Step 2: Construction of the crack network

[0047] 2.1. Through field observation, the attitude, density, and combination patterns of fractures are statistically analyzed, generating a fracture attitude statistical table, a density distribution map, and a combination pattern classification report, including:

[0048] Occurrence measurement involves using a geological compass or a 3D laser scanner to record the strike (e.g., N30°E) and dip (e.g., 45°) of the fracture.

[0049] Density calculation involves measuring the number of cracks per unit area (e.g., 5 cracks / m²) on an outcrop or core sample. 2 );

[0050] Combination pattern recognition distinguishes crack arrangement patterns (such as parallel groups and mesh intersections) through statistical analysis.

[0051] 2.2 Construct a three-dimensional discrete crack network model to transform field observation data into a digital three-dimensional model that reflects the spatial distribution of cracks.

[0052] Specifically, based on the Monte Carlo method, a fracture surface is generated according to statistical parameters (such as strike normal distribution and density Poisson distribution), and non-penetration conditions for fracture length and spacing are set (such as fractures terminating at rock layer interfaces). Then, a three-dimensional DFN model containing fracture location, attitude, and size is generated.

[0053] 2.3. Use finite element software to generate the geometry of non-penetrating crack surfaces, transforming the discrete crack network into a computable geometric entity for mechanical analysis.

[0054] Specifically, in ANSYS or COMSOL, the crack surface is discretized into triangular or quadrilateral elements, contact pairs between the crack surface and the matrix are set to avoid penetration (e.g., friction coefficient 0.6), and a finite element model file containing crack geometry and contact properties is generated.

[0055] 2.4. Import it into the discrete element simulation platform to provide a geometric and physical property basis for subsequent mechanical analysis.

[0056] Specifically, the finite element mesh is imported into 3DEC or PFC3D, the block and joint surfaces are defined, and mechanical parameters are assigned, ultimately resulting in a discrete element model file (.3dec, .dat format) and a parameter configuration table.

[0057] The three-dimensional fracture network model output in step two ensures that the initial conditions of the model are consistent with those of the real reservoir.

[0058] Step 3: Construction of Dynamic Energy Function and Embedding of Crack Stiffness

[0059] 3.1 Based on the static mechanical parameters of step one and the three-dimensional crack network model of step two, and within the Ising machine framework, each mesh element is abstracted as a complex oscillator, the mathematical expression of which is shown in equation (4):

[0060]

[0061] In the complex oscillator represented by equation (4), the mode length |z i | Characterizes the element stiffness (positively correlated with Young's modulus), phase θ i It reflects the azimuth angle of the crack.

[0062] 3.2. Subsequently, the dynamic energy function H′ is constructed, and its mathematical expression is shown in equation (5):

[0063]

[0064] The core terms of the dynamic energy function H′ include:

[0065] ① Linear term h i z i This describes the contribution of the element's own stiffness to the potential energy.

[0066] ② Second-order coupling term Simulate the mechanical synergistic effects between adjacent units (such as energy transfer during crack propagation);

[0067] ③ Third-order interaction term Depicting stress concentration at the intersection of multiple cracks;

[0068] ④ The fracture stiffness term, based on experimental data from fractured rock samples, includes normal stiffness. and shear stiffness The energy consumption of crack surface deformation is incorporated into the energy function.

[0069] This process ensures that the energy function is a real value through complex variable conjugate operations. Its physical meaning is the total potential energy of the reservoir fracture network, which encompasses the combined effects of rock stiffness, fracture interaction, and deformation energy consumption, providing an energy-driven basis for subsequent grid optimization.

[0070] Step 4: Dynamic Gradient Descent and Mesh Size Optimization

[0071] 4.1 First, the computational mesh of the reservoir mechanical model is initialized. This initialization process includes mesh generation and fracture attribute embedding. The mesh generation is based on the reservoir geometric boundaries (such as layer surfaces and fault lines) from step one, dividing the entire reservoir space into an initial computational mesh (such as hexahedral elements, with an initial size r0 = 5m). Subsequently, the three-dimensional fracture network model from step two is superimposed on the initial computational mesh, elements containing fractures are marked, and fracture stiffness parameters are assigned, thereby achieving fracture attribute embedding. The fracture stiffness parameters include the normal stiffness K. n and shear stiffness K s .

[0072] 4.2 Utilizing the dynamic energy function H′ energy gradient Drive mesh parameter adjustment to identify high-energy gradient regions (such as crack intersections and stress concentration areas);

[0073] Furthermore, mesh refinement and coarsening are performed for high-energy gradient regions and low-energy gradient regions, for example:

[0074] In the high energy gradient region ( The mesh size r = 5m was refined to r = 2m to capture local details of crack interactions;

[0075] In the low energy gradient region ( Coarse the mesh to r = 10m to reduce computational load.

[0076] 4.3 Finally, the mesh size is optimized iteratively using a three-loop method:

[0077] ① Position Cyclic: Move the statistical window with a step size Δr = 1m to calculate the equivalent Young's modulus at different positions. and equivalent Poisson ratio Analyze local mechanical parameter fluctuations;

[0078] ② Scale Cycling: Gradually increase the mesh size r and repeat the position cycle until the fluctuation of the equivalent parameter tends to stabilize (i.e., reach the REV size of the characterizing unit cell), thus eliminating the size effect error;

[0079] ③ Orientation Cyclic: Rotate the statistical window direction (e.g., 0°, 45°, 90°) to quantify the influence of anisotropy on mechanical parameters (e.g., the Young's modulus along the crack direction is 20% lower than that in the vertical direction).

[0080] By setting a discrimination index threshold (E) y <0.01 GPa, μ y >0.005), filter the smallest grid size r that meets the accuracy requirements. opt Output the optimal mesh cell size and its equivalent parameter distribution.

[0081] In the above-mentioned discrimination index threshold, Ey The mathematical expression for it is shown in equation (6):

[0082]

[0083] Discriminant threshold μ y The mathematical expression for it is shown in equation (7):

[0084]

[0085] Step 5: Anisotropy Correction and Binarized Mesh Generation

[0086] 5.1 Based on the azimuth correlation term in the dynamic energy function H′, adjust the azimuth weight of the mesh elements to suppress the mechanical parameter deviation in the crack dominant direction (e.g., NE40°); the mathematical expression of this azimuth correlation term is shown in equation (8):

[0087] E θ =∑J ij cos(θ i -θ j )+∑P ijk cos(θ i +θ j -θ k (8)

[0088] 5.2 Subsequently, the complex oscillator z is mapped using binarization. i phase θ i Convert to discrete labels i =sign(Re(z) i )), where s i =+1 indicates a crack development unit, s i =-1 indicates a matrix element; this process transforms the continuous energy distribution into an engineering-usable binary mesh, while simultaneously generating a normalized corner mesh file containing the element size (r). opt It includes equivalent mechanical parameters and fracture activation tags, with a format compatible with commercial software such as Petrel, and can be directly used for reservoir stress field simulation and fracturing design.

[0089] Step Six: Model Validation and Error Analysis

[0090] 6.1 Import the binary mesh model into finite element software (such as ABAQUS) to simulate the reservoir stress distribution during water injection or fracturing, and compare it with downhole stress test data (such as hydraulic fracturing fracture monitoring results) to achieve stress field simulation.

[0091] 6.2 Calculate the relative error between the predicted stress and the measured value, and verify the accuracy of the model.

[0092] 6.3. The statistical calculation time is compared with that of the traditional uniform grid method, and the quantization efficiency is improved by a certain percentage.

[0093] 6.4 Output verification report, optimized grid model and parameter configuration table to provide high-precision and adaptive technical support for oilfield development plan.

[0094] The above are merely specific embodiments of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A method for determining the grid cells for geological reservoir fracture modeling based on oscillator dynamic energy optimization and higher-order interactions, characterized in that, The specific steps include: Step 1: Input the laboratory triaxial experimental data and the well logging P / S wave transit time and density data into the parameter mapping model to obtain the isostatic Young's modulus and Poisson's ratio; at the same time, based on the mapping model, generate a static mechanical parameter database covering various lithologies and stress conditions. Step 2: Using occurrence, density, and combination pattern data, construct a fracture attribute set. Under the multi-scale discrete network framework, use Monte Carlo sampling combined with lithological parameters in the database to generate a three-dimensional discrete fracture set. Then, convert the discrete fracture set into a mesh geometric unit that can be recognized by the finite element / discrete element modeling platform through a general geometric interface. Step 3: Map each grid cell to a complex oscillator node and construct a comprehensive potential energy model to characterize the energy distribution between grid cells; Step 4: Overlay fracture node attributes onto the global reservoir grid to form an attribute mapping layer; based on the potential energy model gradient information, initiate energy-driven optimization to refine the high-energy gradient region and coarsen the low-gradient region. Step 5: Adjust the weights of the high-frequency crack orientation elements according to the orientation term in the potential energy model to balance the anisotropic effect; convert the oscillator phase information into a binary mesh of crack and matrix through a tag mapping mechanism, and output a mesh format file that conforms to commercial software standards. Step Six: Import the binarized mesh into the simulation platform to simulate working conditions such as water injection and fracturing, and compare it with the measured in-situ stress data; through three sub-modules of error statistics, efficiency comparison and stability verification, quantify the improvement of model accuracy and computational performance; automatically generate a closed-loop evaluation report containing error distribution, efficiency improvement rate and optimization suggestions to provide a basis for decision-making.

2. The method according to claim 1, characterized in that, The node properties of the complex oscillator node are composed of the database parameters stiffness and crack orientation.

3. The method according to claim 1, characterized in that, The integrated potential energy model includes self-energy, second-order coupling, higher-order interactions, and stiffness constraints.

4. The method according to claim 1, characterized in that, The potential energy model introduces a stiffness constraint term to achieve a unified expression of the crack surface normal and shear nonlinear characteristics.

5. The method according to claim 1, characterized in that, Step four also includes sequentially verifying and adjusting the mesh size through a multidimensional cyclic correction framework to meet the characterization unit cell accuracy and anisotropic error thresholds.

Images