Adaptive local grid encryption method based on embedded discrete crack

Abstract

The invention relates to a self-adaptive local mesh encryption method based on embedded discrete fractures, which comprises the following steps: aiming at a complex fracture system formed by horizontal well fracturing, establishing a two-dimensional tight reservoir physical model, carrying out initial mesh division, and adopting an EDFM method, two sets of independent mesh systems are adopted for the fractures and a matrix; determining a self-adaptive local grid encryption area; the method comprises the following steps: establishing a self-adaptive local grid encryption criterion based on crack distribution characteristics and physical flow rules, and automatically identifying a key region needing grid encryption: realizing a multi-stage grid encryption strategy, and establishing a tree-shaped hierarchical structure to manage a multi-stage encryption grid system; calculating a connectivity index CI after grid encryption; calculating conductivity, and constructing a conductivity matrix; and solving the pressure and the yield to realize fixed-production pressure solving and fixed-pressure production solving. According to the method, the depiction capability of the embedded discrete fracture model on the flow characteristics near the fracture is improved, and the yield change condition of the fractured reservoir is accurately predicted.

Description

TECHNICAL FIELD

[0001] The present application relates to a numerical simulation method for fractured horizontal well reservoirs in tight oil reservoirs, in particular to an adaptive local grid refinement method based on embedded discrete fractures. BACKGROUND

[0002] Fractured reservoirs play an important role in the world's oil and gas resources. According to statistics, about 60% of the global oil and gas reserves are stored in fractured reservoirs. The fractured reservoirs of carbonate rocks and tight sandstones in Tarim, Ordos and Sichuan basins in China have huge reserves. As the main fluid flow channel, the distribution characteristics, connectivity and percolation characteristics of fractures directly affect the development effect of oil and gas reservoirs. Therefore, accurately simulating the fluid flow behavior in fractured reservoirs is of great significance for formulating a reasonable development plan and improving oil and gas recovery.

[0003] The traditional numerical simulation methods for fractured reservoirs mainly include equivalent continuous medium model, dual medium model and discrete fracture model. The equivalent continuous medium model equivalent the permeability of fractures to the grid block, which cannot accurately describe the geometry and flow characteristics of fractures. The dual medium model (Warren J.E., Root P.J., 1963, 3 (3): 245-255) considers the matrix and fracture systems, but still uses the continuous medium assumption, which is difficult to depict the heterogeneity of complex fracture network. The discrete fracture model (Karimi-Fard M., Durlofsky L.J., Aziz K., 2004, 7 (3): 241-251) can explicitly represent each fracture, but requires grid partitioning for fractures, resulting in a sharp increase in the number of grids and high computational cost. SPE Journal SPE Reservoir Evaluation & Engineering In recent years, the embedded discrete fracture model (EDFM) as a new type of fracture modeling method has attracted widespread attention (Li L., Lee S.H., 2008, 227 (3): 2082-2098). This method embeds fractures as low-dimensional units into the background grid without explicit grid partitioning for fractures, greatly reducing the number of grids. However, when dealing with dense fracture areas or fracture and wellbore intersection areas, EDFM has difficulty in accurately capturing local flow details due to the large scale of the background grid, which affects the simulation accuracy. However, using uniformly refined grids in the entire simulation area will greatly increase the computational load, losing the computational advantage of EDFM.

[0004] In recent years, the embedded discrete fracture model (EDFM) as a new type of fracture modeling method has attracted widespread attention (Li L., Lee S.H., 2008, 227 (3): 2082-2098). This method embeds fractures as low-dimensional units into the background grid without explicit grid partitioning for fractures, greatly reducing the number of grids. However, when dealing with dense fracture areas or fracture and wellbore intersection areas, EDFM has difficulty in accurately capturing local flow details due to the large scale of the background grid, which affects the simulation accuracy. However, using uniformly refined grids in the entire simulation area will greatly increase the computational load, losing the computational advantage of EDFM. Journal of Computational Physics

[0005] ​​Local grid refinement (LGR) is an adaptive grid technique that refines the grid in key areas, which can control the total number of grids while ensuring the calculation accuracy. However, the combination of LGR and EDFM in the prior art still has the following technical problems: (1) how to accurately determine the area that needs to be locally encrypted, both to ensure the capture of key flow characteristics and to avoid excessive refinement; (2) how to handle the grid continuity problem between grids of different scales to ensure flow conservation and pressure continuity; (3) how to accurately and efficiently calculate the connectivity index of fracture-matrix, fracture-fracture and matrix-matrix in the encrypted grid, and then calculate the conductivity.

[0006] Therefore, the application provides an adaptive local grid refinement method based on embedded discrete fractures, which solves the contradiction between calculation accuracy and calculation efficiency in numerical simulation of dense oil reservoirs or shale oil and gas reservoirs, accurately depicts the flow characteristics of complex fracture networks, and accurately predicts the pressure and flow changes of the fracture system. The method has important value for forming a high-efficiency numerical simulation theory method for multi-scale fracture systems and building a reliable dense oil reservoir development effect prediction system, further improves the accuracy and calculation efficiency of numerical simulation of dense oil reservoirs, and provides important theoretical support for efficient development of complex fracture reservoirs such as dense oil reservoirs. SUMMARY

[0007] The application aims to provide an adaptive local grid refinement method based on embedded discrete fractures, which solves the contradiction between calculation accuracy and calculation efficiency in existing numerical simulation of dense oil reservoirs, accurately depicts the flow characteristics of complex fracture networks, and accurately predicts the pressure and flow changes of the fracture system.

[0008] The technical scheme adopted by the application to solve the technical problems is that the adaptive local grid refinement method based on embedded discrete fractures comprises the following steps: Step 1, for a complex fracture system formed by horizontal well fracturing, a two-dimensional dense oil reservoir physical model is established, and initial grid division is performed. The EDFM method is used, and the fracture and the matrix use two independent grid systems; Step 2, determine the adaptive local grid refinement area; based on the fracture distribution characteristics and physical flow law, establish an adaptive local grid refinement criterion to automatically identify the key area that needs to be grid-encrypted: Step 3, implementation of multi-level grid refinement strategy: before grid refinement, the storage space needs to be expanded; a tree-like hierarchical structure is established to manage the multi-level encrypted grid system, the “parent-child” grid relationship, grid level identification and global connection rules are clearly defined, and the “inactive” grid ID is reserved; Step 4, calculate the connectivity index CI after grid refinement: Step 5: Calculate the conductivity and construct the conductivity matrix. The matrix dimensions are... N total × N total ,in N total = N m + N f This is the sum of the number of active matrix grids and the number of crack segments; Step 6: Solve for pressure and output; solve the linear equations of pressure and source / sink terms in combination with boundary conditions to obtain the pressure solution; given different boundary conditions, realize pressure calculation for constant output and output calculation for constant pressure.

[0009] Step one of the above scheme is specifically as follows: (1) Define the computational domain parameters; Define the computational domain in a two-dimensional plane, considering the complex fracture mesh formed by hydraulic fracturing: Computational domain setting: reservoir length L x Determined based on well control area and fracturing scale; reservoir width L y The extent of the fracturing fracture should be considered in determining its location. Boundary condition type: The outer boundary is a closed boundary with no flow; the wellbore location of the inner boundary is treated using the equivalent wellbore radius. Coordinate system establishment: x The axis is along the horizontal well direction. y Axis perpendicular to x Axis, along the direction of the main crack; (2) Create the basic grid system; Determining the number of grid cells: , ; in, N x , N y They are the initial ones. x and y Number of directional grids; ∆ x 0, ∆ y 0 represents the initial value. x and y Orientation grid step size; Grid numbering and storage: A row-major natural numbering method is used, with each grid ( i , j The global number of ) is: ; (3) Establish a geometric model of the artificial crack; For fracture systems formed by hydraulic fracturing, each fracture is discretized into multiple continuous fracture segments, and a two-level numbering system is established: Level 1: Crack Number ; Level 2: Crack segment numbering, each crack i Discretized into N i s part: ; Global numbering of crack segments: ; Among them, L i f The length of each crack is in meters (m); Δl f The length of each crack segment, in meters; j For this section in the 1st i The number in the crack.

[0010] Step two in the above scheme is specifically as follows: A step-by-step encryption strategy is adopted, starting from level 1 and continuing until the maximum encryption level (max_level) is reached. For each encryption level level ϵ{1,2,...,max_level}, perform the following steps: (1) Candidate grid filtering; For each active cell, three conditions must be met simultaneously to be marked as to be encrypted: Condition 1: Level matching, only the grid of the previous level can be encrypted at the current level (grid.level=level−1). Condition 2: Size constraint, min(∆) x , ∆ y )>∆min, where Δmin is the minimum allowable mesh size; Condition 3: Crack intersection judgment; (2) Crack-mesh intersection detection; The intersection detection uses a line segment-rectangle intersection algorithm, and the judgment logic is as follows: Let the endpoint of the crack segment be ( x 1, y 1) and ( x 2, y 2), the grid boundary is [ xmin , xmax ]×[ ymin , ymax ] (a) Endpoints are within the grid: ; (b) Cracks horizontally penetrate the grid: ; (c) Cracks penetrate the grid perpendicularly: ; (d) Cracks diagonally penetrate the grid: For inclined cracks, parametric line segment representation and boundary intersection tests are used. The crack segment parametric equations are as follows: ; Calculate the intersection parameters of the crack segment with each of the four edges of the mesh. t If it exists t If the intersection points are on the grid boundary and are ∈[0,1], then they are considered to be intersecting; If any of the above conditions are met, the two objects are considered to be intersecting. (3) Encryption terminated; After each level of encryption, check if any grids have been marked. If no grids need to be encrypted, terminate the encryption process. If grids need to be encrypted, perform quarter-encryption and continue to the next level.

[0011] Step three in the above plan specifically refers to: (1) Before performing grid encryption, storage space needs to be expanded; the number of new grids = the number of grids to be encrypted × 4, and the coordinates, size, ID, and level are expanded; (2) The "parent" mesh is deactivated (active[parent_idx]=0) and the "child" mesh is generated (active[new_idx]=1); for each "parent" mesh to be encrypted, the following steps are performed: Step 1: Deactivate the "parent" grid. The "parent" grid will no longer participate in the calculation, but its information will be retained for hierarchical management. Step 2: Calculate the geometric parameters of the "child" mesh; the range of the "parent" mesh is [ xmin , xmax ]×[ ymin , ymax Calculate half-size: , ; Step 3: Generate 4 "sub" grid coordinates. Subgrid 1 (bottom left):[ xmin , xmin +d x ]×[ ymin , ymin +d y ] Subgrid 2 (bottom right):[ xmin +d x , xmax ]×[ ymin , ymin +d y ] Subgrid 3 (top left):[xmin , xmin +d x ]×[ ymin +d y , ymax ] Subgrid 4 (top right):[ xmin +d x , xmax ]×[ ymin +d y , ymax ] (3) Assigning sub-grid attributes: For each sub-grid, set attributes such as coordinates, size, global ID, and level; (4) The encrypted mesh hierarchy is a tree-like hierarchy: After encryption, the total number of active grids is: ; for m Level 1 encryption: ; in, N initial This is the initial number of grid cells; N refined The number of encrypted grid cells; N level1 The number of grid cells for the first level of encryption; N level2 The number of grid cells for the second level of encryption; (5) Grid number management: new grid IDs are continuously assigned starting from the current maximum ID + 1.

[0012] Step four in the above scheme is specifically as follows: (1) Crack-matrix connectivity index; For the connection between the crack segment f and the matrix mesh m CI The calculation formula is: ; in, CI fm The connectivity index between the crack and the matrix mesh; l eff This is the effective contact length of the crack segment; ω f The crack opening; d cm The characteristic distance is the vertical distance from the center of the crack segment to the center of the matrix grid. Let the coordinates of the starting point of the crack segment be ( x 1, y 1), The coordinates of the endpoint are ( x 2,y 2), the coordinates of the center point are ( x c,f , y c,f The coordinates of the matrix center point are ( x c,m , y c,m ),but: ; ; (2) Matrix-matrix connectivity index; ; ; in, CI mm The connectivity index between the matrix and the matrix grid; A mm The area of ​​the common contact surface between the two matrix meshes; d ij This is the equivalent flow distance; d i matrix mesh i The distance from the center point to the common contact surface; d j matrix mesh j The distance from the center point to the common contact surface; (3) Crack-crack intersection detection; Set crack section i The starting coordinates are P1 = ( x 1, y 1), The endpoint coordinates are P2 = ( x 2, y 2) Cracked section j The starting coordinates are Q1 = ( x 3, y 3), The endpoint coordinates are Q2 = ( x 4, y 4); For the crack segment i : ; For the crack segment j : ; in, ; make ,when D When =0, the two crack segments are parallel; when hour, ; ; When 0 < s <1 and 0< t <1 indicates that two crack segments intersect; Wherein, s1 and s2 are the crack segments. i , j The direction vector; D It is the cross product of two direction vectors; s , t These are the crack segments. i , j Position parameters on, s ∈[0,1], t ∈[0,1].

[0013] Step five in the above scheme is specifically as follows: (1) Calculate the conductivity between cracks and matrix, matrix and matrix, and cracks and cracks; Crack-matrix conductivity: ; in, T fm The conductivity between the crack and the matrix; K m For matrix permeability; K f The permeability of the crack; μ For fluid viscosity; Matrix-matrix conductivity: ; in, T mm The conductivity between the matrix and the matrix; Crack-crack conductivity: The conductivity of adjacent segments within the same crack, for two adjacent crack segments in the same crack. i and i +1, conductivity reflects one-dimensional linear flow within the crack: ; in, T ff,i,i+1 The conductivity between two adjacent crack segments within the same crack; h ∆ represents the oil layer thickness. x f This is the distance between the centers of the two crack segments; The conductivity of the intersecting segments of different cracks, when the segments of two different cracks... i and j When they intersect, the fluids exchange at the intersection point, and the conductivity is calculated using a harmonic average. ; The unilateral conductivity of each fracture segment is: ; in, T ff,i,j The conductivity of the intersection segment of different cracks; α i For crack segment i conductivity; α j For crack segment j Conductivity; 0.5⋅ l i This represents the distance from the center of the crack segment to the intersection point; (2) Construct the conductivity matrix.

[0014] The method for constructing the conductivity matrix in the above scheme is as follows: First, the diagonal elements: for the matrix mesh m, the diagonal elements are the sum of the time accumulation term and the conductivity of all connections. ; For crack segment f (global number) N m +f), the diagonal elements are: ; in, V m The volume of the matrix mesh; V f The volume of the crack segment; φ m For matrix porosity; φ f For crack porosity; C t,m The compressibility coefficient of the matrix; C t,f ∆ is the compressibility coefficient of the crack; t For time intervals; N m (m) is the set of matrix grids adjacent to matrix grid m; N f (m) is the set of crack segments connected to the matrix mesh m; N m (f) is the matrix mesh set connected to the fracture segment f; N f (f) is the set of crack segments connected to crack segment f; Then there are the off-diagonal elements, for any connected units. i and j Off-diagonal elements are negative values ​​of conductivity: ; The negative sign reflects the physical property of fluid flowing from a high-pressure unit to a low-pressure unit, when the conductivity satisfies the reciprocity T. ij =T ji When the matrix exhibits symmetry.

[0015] The method for defining computational domain parameters in the above scheme is as follows: Define the computational domain in a two-dimensional plane, considering the complex fracture mesh formed by hydraulic fracturing: Computational domain setting: reservoir length L x Determined based on well control area and fracturing scale; reservoir width L y The extent of the fracturing fracture should be considered in determining its location. Boundary condition type: The outer boundary is a closed boundary with no flow; the wellbore location of the inner boundary is treated using the equivalent wellbore radius. Coordinate system establishment: x The axis is along the horizontal well direction. y Axis perpendicular to x The axis is along the direction of the main crack.

[0016] Beneficial effects: (i) This invention establishes an adaptive local mesh refinement method based on the geometric intersection judgment of fracture-matrix mesh. Through the triple criteria of level matching, size restriction and intersection judgment, it realizes the automatic identification and progressive refinement of the refinement area, avoiding the subjectivity of relying on manual experience to specify the refinement area in the traditional method. This enables the formation of a multi-level mesh structure from coarse to fine around the complex fracture network, which significantly improves the accuracy of local flow characterization in fracture-dense areas and near-wellbore zones. At the same time, it controls the overall mesh number within a reasonable range, achieving the optimal balance between computational accuracy and computational efficiency.

[0017] (II) This invention establishes a complete method system for managing multi-level encrypted grids using a tree-like hierarchical structure, clarifies the "parent-child" grid relationship, grid level identifiers, and global numbering rules, and ensures the integrity and data consistency of the grid hierarchy tree by retaining the IDs of "deactivated" grids. It solves the problems of chaotic multi-scale grid management and low query efficiency in the prior art, and provides an efficient data structure foundation for rapid access, dynamic updates, and visualization of grid information. At the same time, it lays the technical framework for adaptive grid optimization and parallel computing.

[0018] (III) This invention proposes an improved connectivity index calculation method for the connection of matrix grids of different sizes after multi-level densification. It adopts accurate overlap area calculation and introduces equivalent flow distance, which effectively improves the flow conservation problem caused by grid size difference and ensures the accuracy of pressure and flow calculation.

[0019] (iv) This invention establishes a unified conductivity calculation framework for three types of connection relationships: crack-matrix, crack-crack, and matrix-matrix. It proposes a crack intersection detection algorithm based on line segment parameterized equations. Through harmonic average and series resistance models, it ensures the physical rationality of fluid exchange between different types of units. It systematically solves the problem of calculating connection relationships when embedded discrete crack models are coupled with multi-level densified meshes, and lays a mathematical foundation for constructing sparse symmetric positive definite conductivity matrices and using efficient iterative solvers.

[0020] (V) This invention forms a complete technical process from physical model establishment, adaptive local mesh generation, connection relationship calculation to conductivity matrix assembly, filling the research gap of local mesh refinement and embedded discrete fracture model coupling method, providing a reliable numerical simulation tool for the efficient development of unconventional oil and gas resources such as tight oil and gas reservoirs, deep carbonate reservoirs, and shale oil and gas, and promoting the development of fractured reservoir numerical simulation technology towards intelligence and automation, which has important theoretical value and broad application prospects.

[0021] (vi) This invention solves the following key technical problems existing in the numerical simulation of tight oil reservoirs: Traditional uniform grid methods face a trade-off between accuracy and efficiency when dealing with multi-scale fracture systems, requiring optimized allocation of computational resources. Embedded discrete fracture models suffer from insufficient accuracy in fracture-dense regions and near-wellbore areas, necessitating improved characterization of local flow features. After mesh refinement, the issue of mesh continuity between multi-scale grids needs to be addressed, ensuring flow conservation and pressure continuity. Efficient computation of fracture-matrix, fracture-fracture, and matrix-matrix connectivity relationships in complex fracture networks is addressed, providing a unified framework for conductivity calculation. Attached Figure Description

[0022] Figure 1 This is a schematic diagram of the multi-level encryption of the present invention; Figure 2 This is a schematic diagram of the initial mesh generation and fracture distribution of the reservoir physical model of the present invention; Figure 3 This is a flowchart illustrating the technical process of the present invention. Figure 4 This is a schematic diagram of the 4-level local mesh encryption of the present invention; Figure 5 This is a schematic diagram of the grid numbering management system of the present invention; Figure 6 This invention provides a comparison of the pressure fields during 100 days of production using unencrypted and Level 4 encrypted methods. Figure 7 This is a comparison chart of daily output and cumulative output under different grid configuration schemes of the present invention; Figure 8This is a comparison chart of daily and cumulative output under different encryption levels of the present invention; Figure 9 This is a schematic diagram of a 4-level grid encryption method in an example application of the present invention. Detailed Implementation

[0023] The invention will be further described below with reference to the accompanying drawings: This adaptive local mesh refinement method based on embedded discrete cracks includes the following steps: (i) Establishing a two-dimensional physical model of a tight oil reservoir and performing initial mesh generation. A two-dimensional planar reservoir model was established to address the complex fracture systems formed by horizontal well fracturing. The EDFM method was employed, with two independent mesh systems for the fractures and the matrix.

[0024] (1) Define computational domain parameters Define the computational domain in a two-dimensional plane, considering the complex fracture mesh formed by hydraulic fracturing: Computational domain setting: reservoir length ( L x The reservoir width is determined based on the well control area and fracturing scale. L y The extent of the fracturing fracture should be considered in determining the extent of the fracture.

[0025] Boundary condition type: The outer boundary is a closed boundary (no flow); the wellbore location of the inner boundary is treated using the equivalent wellbore radius.

[0026] Coordinate system establishment: x The axis is along the horizontal well direction. y Axis perpendicular to x Axis (along the direction of the main crack).

[0027] (2) Create a basic grid system Determining the number of grid cells: ,

[0028] in, N x , N y They are the initial ones. x and y Number of directional grids, ∆ x 0, ∆ y 0 represents the initial value. x and y The directional grid step size, in meters (m). Setting the initial grid step size requires a balance between computational accuracy and efficiency.

[0029] Grid numbering and storage: A row-major natural numbering method is used, with each grid ( i ,j The global number of ) is:

[0030] (3) Establishing a geometric model of artificial cracks The core idea of ​​the embedded discrete fracture model is to directly embed the fracture as a one-dimensional linear unit into the existing background mesh, without explicit meshing of the fracture. This method maintains the regularity of the background mesh while accurately characterizing the geometric features and flow properties of the fracture. For fracture systems formed by hydraulic fracturing, each fracture is discretized into multiple continuous fracture segments, establishing a two-level numbering system: Level 1: Crack Number

[0031] Level 2: Crack segment numbering, each crack i Discretized into N i s part:

[0032] Global numbering of crack segments: ; Among them, L i f The length of each crack is in meters (m); Δl f The length of each crack segment, in meters; j For this section in the 1st i The number in the crack.

[0033] (II) Determination of Adaptive Local Mesh Refinement Region Based on crack distribution characteristics and physical flow laws, an adaptive local mesh refinement criterion is established to automatically identify key areas requiring mesh refinement. A step-by-step encryption strategy is adopted, starting from level 1 and continuing until the maximum encryption level (max_level) is reached. For each encryption level For levelϵ{1,2,...,max_level}, perform the following steps: (1) Candidate grid selection For each active cell, three conditions must be met simultaneously to be marked as pending encryption: Condition 1: Level matching, only the grid of the previous level can be refined at the current level (grid.level=level−1), ensuring a hierarchical structure of progressive refinement.

[0034] Condition 2: Size constraint, min(∆) x , ∆ y >∆min, to prevent over-refinement.

[0035] Where Δmin is the minimum allowable mesh size, in meters.

[0036] Condition 3: Crack intersection determination. This is the core step and requires inspection of all crack segments.

[0037] (2) Crack-mesh intersection detection The intersection detection uses a line segment-rectangle intersection algorithm, and the judgment logic is as follows: Let the endpoint of the crack segment be ( x 1, y 1) and ( x 2, y 2), the grid boundary is [ xmin , xmax ]×[ ymin , ymax ] (a) Endpoints are within the grid:

[0038] (b) Cracks horizontally penetrate the grid:

[0039] (c) Cracks penetrate the grid perpendicularly:

[0040] (d) Cracks diagonally penetrate the grid: For general inclined cracks, parametric line segment representation and boundary intersection testing are used. The crack segment parametric equations are as follows:

[0041] Calculate the intersection parameters of the crack segment with each of the four edges of the mesh. t If it exists t If the intersection point is on the grid boundary and the intersection point is ∈[0,1], then the intersection is determined to be intersecting.

[0042] If any of the above conditions are met, the two objects are considered to be intersecting.

[0043] (3) Encryption Termination After each level of encryption, check if any grids have been marked. If no grids need to be encrypted, the encryption process terminates. If a grid needs to be encrypted, perform a four-part encryption and continue to the next level.

[0044] (III) Implementation of Multi-level Grid Encryption Strategy (1) Before performing grid encryption, the storage space needs to be expanded. The number of new grids = the number of grids to be encrypted × 4, and the arrays of coordinates, dimensions, IDs, and levels are expanded.

[0045] (2) The "parent" grid is deactivated (active[parent_idx]=0) and the "child" grid is generated (active[new_idx]=1). During program development, to facilitate calculation and grid searching, the active array of all deactivated grids is 0, and the active array of all active grids is 1. For each parent grid to be encrypted, the following steps are performed: Step 1: Deactivate the "parent" grid. The "parent" grid no longer participates in the calculation, but its information is retained for hierarchical management.

[0046] Step 2: Calculate the geometric parameters of the "child" mesh. The range of the "parent" mesh is [ xmin , xmax ]×[ ymin , ymax ] Calculate the half-size: ,

[0047] Step 3: Generate 4 "sub" grid coordinates. Subgrid 1 (bottom left):[ xmin , xmin +d x ]×[ ymin , ymin +d y ] Subgrid 2 (bottom right):[ xmin +d x , xmax ]×[ ymin , ymin +d y ] Subgrid 3 (top left):[ xmin , xmin +d x ]×[ ymin +d y , ymax ] Subgrid 4 (top right):[ xmin +d x , xmax ]×[ ymin +d y , ymax ] (3) Assigning values ​​to the "sub" grid attributes. Set the coordinates, size, global ID, level and other attributes for each "sub" grid.

[0048] (4) The encrypted mesh hierarchy This ultimately forms a tree-like hierarchical structure (taking two-level encryption as an example, where level represents the encryption level): Level 0: Initial mesh (partially "deactivated") Level 1: First-level "sub" mesh (partially "deactivated") Level 2: Second-level "sub" grid (all activities) Level 1: Other first-level "sub" grids (activities) After encryption, the total number of active grids is:

[0049] Generally speaking, for m Level 1 encryption:

[0050] in, N initial This is the initial number of grid cells; N refined The number of encrypted grid cells; N level1 The number of grid cells for the first level of encryption; N level2 The number of grid cells for the second level of encryption.

[0051] (5) Grid number management The ID allocation principle for new grids is as follows: new grid IDs are continuously allocated starting from the current maximum ID + 1.

[0052]

[0053] Among them, ID new The new grid number is max(ID), which is the maximum grid number.

[0054] The ID of the "deactivated" parent grid is retained and will not be reused. The IDs of all grids (active and "deactivated") remain unchanged throughout the simulation. For a detailed explanation of the grid numbering process, please refer to the grid numbering diagram.

[0055] (iv) Calculation of CI after mesh refinement The connectivity index (CI) is a key parameter in the EDFM method that quantitatively describes the fluid exchange capacity between fractures and matrix, fractures and fractures, and matrix and matrix, directly affecting the accuracy of conductivity calculation. In multi-level refined mesh systems, CI The calculation is more complex and requires consideration of the geometric characteristics of different mesh levels. CI This characterizes the geometric connection between two computational units, reflecting the ease with which fluid flows from one unit to another. Its definition is based on Darcy's law and fundamental principles of fluid mechanics, taking into account geometric factors such as contact area and flow distance.

[0056] (1) Crack-matrix connectivity index For the connection between the crack segment f and the matrix mesh m CI The calculation formula is:

[0057] in, CI fm The connectivity index between the crack and the matrix mesh is m; l eff The effective contact length of the crack segment is in meters (m). ω f Let the crack aperture be m; d cm The characteristic distance is the vertical distance from the center of the crack segment to the center of the matrix mesh, in meters. Let the coordinates of the starting point of the crack segment be ( x 1, y 1), The coordinates of the endpoint are ( x 2, y 2) The coordinates of the center point are ( x c,f , y c,f The coordinates of the matrix center point are ( x c,m , y c,m ),but:

[0058]

[0059] (2) Matrix-matrix connectivity index Considering that after local mesh refinement, there may be a situation where a large matrix mesh is adjacent to multiple small matrix meshes, the formula for calculating the connectivity index between different matrices is:

[0060]

[0061] in, CI mm The connectivity index between the matrix and the matrix mesh is m; A mm m is the area of ​​the common contact surface between the two matrix meshes. 2 ; d ij The equivalent flow distance is in meters (m). d i matrix mesh i The distance from the center point to the common contact surface, in meters (m). d j matrix meshj The distance from the center point to the common contact surface, in meters (m).

[0062] (3) Crack-crack intersection detection Intersection detection of cracks lays the foundation for establishing the conductivity matrix between cracks. Let crack segments... i The starting coordinates are P1 = ( x 1, y 1), The endpoint coordinates are P2 = ( x 2, y 2) Cracked section j The starting coordinates are Q1 = ( x 3, y 3), The endpoint coordinates are Q2 = ( x 4, y 4).

[0063] For the crack segment i :

[0064] For the crack segment j :

[0065] in,

[0066] make ,when D When =0, the two crack segments are parallel; when hour,

[0067]

[0068] When 0 < s <1 and 0< t If the value is less than 1, it can be determined that two crack segments intersect.

[0069] Wherein, s1 and s2 are the crack segments. i , j The direction vector; D It is the cross product of two direction vectors; s , t These are the crack segments. i , j Position parameters on, s ∈[0,1], t ∈[0,1].

[0070] (V) Calculation of conductivity and construction of matrix (1) Calculation of conductivity between crack and matrix, matrix and matrix, and crack and crack. Crack-matrix conductivity: The conductivity between the crack segment f and the matrix mesh m was determined using a harmonic averaging method based on the connectivity index.

[0071] in, T fm m represents the conductivity between the crack and the matrix. 3 / (Pa·s); K m m is the matrix permeability. 2 ; K f m is the fracture permeability. 2 ; μ Let be the fluid viscosity, Pa·s. This formula reflects the series flow characteristics of the fluid from the matrix to the fracture. Harmonic averaging ensures that when the matrix permeability is much smaller than the fracture permeability, the total conductivity is mainly controlled by the matrix, which is consistent with physical reality.

[0072] Matrix-matrix conductivity:

[0073] in, T mm m is the conductivity between matrices. 3 / (Pa·s).

[0074] Crack-crack conductivity: Conductivity of adjacent segments within the same crack. For two adjacent crack segments within the same crack. i and i +1, conductivity reflects one-dimensional linear flow within the crack:

[0075] in, T ff,i,i+1 m is the conductivity between two adjacent crack segments within the same crack. 3 / (Pa·s); h Here is the oil layer thickness, in meters (m); ∆ x f The distance between the centers of the two crack segments is in meters (m).

[0076] Conductivity of the intersecting segments of different cracks. When two different cracks intersect... i and j When they intersect, fluids can exchange at the intersection point. Conductivity is calculated using a harmonic average.

[0077] The unilateral conductivity of each fracture segment is:

[0078] in, T ff,i,j m represents the conductivity of the intersection segments of different cracks. 3 / (Pa·s); α i For crack segment i conductivity, m 3 / (Pa·s); α j For crack segment j conductivity, m 3 / (Pa·s); 0.5⋅ l i The distance from the center of the crack segment to the intersection point is represented in meters (m).

[0079] (2) Construction of the conductivity matrix The conductivity matrix (coefficient matrix) A is the core component that transforms discrete finite volume equations into a system of linear equations. The matrix dimension is... N total × N total ,in N total = N m + N f This represents the sum of the number of active matrix grids and the number of fracture segments. The construction of the conductivity matrix is ​​explained from both diagonal and off-diagonal elements.

[0080] First, the diagonal elements: for the matrix mesh m, the diagonal elements are the sum of the time accumulation term and the conductivity of all connections.

[0081] For crack segment f (global number) N m +f), the diagonal elements are:

[0082] in, V m m is the volume of the matrix mesh. 3 ; V f Let m be the volume of the crack segment. 3 ; φ m For matrix porosity; φ f For crack porosity; Ct,m Pa is the compressibility coefficient of the matrix. -1 ; C t,f Let be the compressibility coefficient of the crack. Pa -1 ;∆ t The time interval is s; N m (m) is the set of matrix grids adjacent to matrix grid m; N f (m) is the set of crack segments connected to the matrix mesh m; N m (f) is the matrix mesh set connected to the fracture segment f; N f (f) is the set of crack segments connected to crack segment f.

[0083] Next are the off-diagonal elements, for any connected units. i and j (Matrix-matrix, crack-matrix, or crack-crack), off-diagonal elements are negative values ​​of conductivity:

[0084] The negative sign reflects the physical property of fluid flowing from a high-pressure unit to a low-pressure unit. When the conductivity satisfies reciprocity (T... ij =T ji When ), the matrix exhibits symmetry.

[0085] (vi) Solving for pressure and output Based on the conductivity matrix obtained in technical solution (V), the pressure matrix is ​​constructed, and the linear equations between pressure and source / sink terms are solved in conjunction with boundary conditions to obtain the pressure solution. Given different boundary conditions, pressure calculation for constant production and production calculation for constant pressure can be achieved. Through the above conductivity calculation and matrix construction methods, this invention realizes an accurate mathematical description of the multi-scale crack-matrix coupled system, laying the foundation for numerical solutions.

[0086] The following is a detailed description in conjunction with the accompanying drawings: (I) Establishing an initial model of a tight oil reservoir and performing mesh generation using this invention. This invention establishes a two-dimensional planar numerical simulation model of the oil reservoir for the "cross"-shaped fracture system formed by vertical well fracturing. Considering the vertical well fracturing region of a low-permeability oil reservoir, this region is a homogeneous isotropic reservoir with a square planar shape. This regular geometric shape facilitates mesh generation and boundary condition handling. Embedded discrete fractures are used for mesh generation. Reservoir parameters and reservoir side lengths are given. L oil layer thickness h , wellbore radius r w Matrix permeabilityK m Original formation pressure P i fluid viscosity μ Half a length of crack L f Crack width w f Porosity φ Crack permeability K f . Figure 2 This diagram illustrates the initial uniform grid division and fracture distribution based on the reservoir characteristics of a cross-shaped fractured well. The figure clearly shows a 400m × 400m square area divided into 400 uniform grid cells (20 × 20 cells each). Grid numbers are row-major, numbered sequentially from left to right and bottom to top; therefore, the bottom left grid is numbered 1, and the top right grid is numbered 400. The Embedded Discrete Fracture Model (EDFM) method is employed, embedding fractures as one-dimensional linear units within a two-dimensional background grid. The advantages of this method are: fractures are independent of the matrix grid, fractures can pass through the background grid at any angle without generating complex unstructured grids, and computational efficiency and geometric flexibility are maintained. After hydraulic fracturing in a vertical well, two mutually perpendicular main fractures are formed, exhibiting a typical cross-shaped distribution. Each fracture segment is considered an independent computational unit, possessing its own volume, porosity, permeability, and other physical properties, and is connected to adjacent matrix grids and fracture segments through a connectivity index. The intersection of the two fractures marks the location of the wellbore. Figure 2 This initial mesh generation scheme provides the foundation for subsequent adaptive refinement. It is foreseeable that the mesh around the crack will be progressively refined to capture flow details near the crack, while corner regions far from the crack will retain a coarse mesh to conserve computational resources. This initial uniform mesh provides a unified starting point and reference benchmark for adaptive refinement.

[0087] (II) Technical process for implementing adaptive local mesh refinement using this invention. Based on the established initial uniform mesh model, this invention designs an adaptive local mesh refinement algorithm for the system, setting the maximum refinement level max_level=4, the minimum allowable mesh size ∆min=1.25m, and the geometric factor α=0.2. Figure 2This is a flowchart of the adaptive mesh densification algorithm of the present invention, which details the complete computation process from the initial mesh to the multi-level densification mesh. The specific steps include: first, starting the algorithm and initializing the computational environment, reading input files such as reservoir parameters, initial mesh information, and fracture geometry data; then, performing initial mesh generation, establishing a uniform Cartesian background mesh, marking the densification level of all meshes as level=0, and the active state as active=1; next, embedding the discretized fracture segments into the background mesh, traversing all matrix meshes for each fracture segment to determine whether the fracture segment crosses the mesh, and establishing the connection relationship between the fracture segment and the mesh using parameterized line segment representation and boundary intersection testing methods; setting densification parameters and initializing the densification loop control variable level=1; detecting whether the mesh intersects with the fracture, and performing a triple judgment by traversing all currently active meshes. The process includes checking if the grid level matches (grid.level=level-1), if the grid size meets the limit (min(∆x,∆y)>∆min), and if the grid intersects with the crack. If all three conditions are met, the grid is added to the list of grids to be encrypted. A four-part encryption operation is performed, deactivating the parent grid (active=0) and generating four child grids arranged in the bottom left, bottom right, top left, and top right positions. Each child grid is assigned a consecutive global ID and its level, active, and parent_id attributes are set. The intersection relationship between the crack and the child grid is recalculated, and the connection table is updated. The maximum level is checked: if the current encryption level is less than max_level, it returns to the next round of encryption if not, otherwise the process ends and the final multi-level encrypted grid system is output. This flowchart uses a progressive strategy to ensure a smooth transition in grid size. Based on geometry-driven criteria, it automatically identifies complex regions near cracks to achieve intelligent encryption, only refining areas that require refinement to avoid a surge in computation caused by global uniform refinement. The algorithm framework has good scalability and can be easily extended to 3D cases and complex crack networks.

[0088] (III) Achieving adaptive local mesh encryption effect using this invention. Based on Figure 3 The algorithm flow shown performs multi-level encryption. First, it performs level 1 encryption, traversing 400 initial grids and performing a triple check. 44 grids are then marked as needing encryption. These grids are distributed in a cross shape, covering the two main cracks and their adjacent areas. After being divided into quarters, 176 sub-grids are generated, each 10m × 10m in size. The total number of active grids after encryption is... N active =532. Figure 3This is a schematic diagram of the mesh distribution after Level 1 local mesh refinement. The diagram clearly shows two perpendicular band-shaped refinement zones forming around the "+" shaped crack. The mesh lines in these refinement zones are densely packed at 10m × 10m, while the four corner areas away from the crack retain the original 20m × 20m mesh. The spatial resolution along the crack line is improved from 20m to 10m, achieving optimized allocation of computational resources. Based on Level 1 refinement, Levels 2, 3, and 4 refinements are performed, refining 84, 164, and 324 meshes respectively to generate sub-mesh sizes of 5m × 5m, 2.5m × 2.5m, and 1.25m × 1.25m, resulting in a final total of 2248 active meshes. Figure 4 This is a comprehensive visualization of the multi-level refinement. From the center of the wellbore outwards, the mesh size exhibits a clear hierarchical structure: 1.25m → 2.5m → 5m → 10m → 20m, a total of five levels. The size ratio between adjacent levels is 2:1, forming a gradual mesh distribution that avoids mesh distortion. The innermost 1.25m mesh, located adjacent to the fracture, accurately captures the flow characteristics at the fracture surface. Compared to the 102,400 meshes required for global uniform refinement to 1.25m, this method achieves comparable computational accuracy using only 2,248 meshes, reducing the mesh count by 97.8%.

[0089] (iv) Using this invention for grid numbering management and hierarchical structure maintenance. During multi-level encryption, a large number of "parent-child" relationship grids and "deactivated" grids are generated. This invention employs a tree-like hierarchical structure management method to effectively manage this grid information. Figure 6 This diagram illustrates the grid numbering management process of the present invention. A 5×5 grid area is selected as an example to demonstrate the numbering process of the encrypted grid. The grid ID (G) uses globally unique identifiers G1-G25, assigned in the order of generation, and the ID never changes once assigned. The red line represents a crack segment located in grid G13. This grid is encrypted according to the grid encryption principle, and the four generated "sub" grids are further numbered as G26, G27, G28, and G29. However, for the continuity of the subsequent conductivity matrix construction, since G13 has been "deactivated," its matrix index is extended to G14 (i.e., the matrix index of grid number 14 is 13), and so on. Through this grid numbering and information management method, the present invention can ensure the uniqueness of the grid number and the continuity of the matrix index, while completely recording all grid information, providing basic data support for subsequent construction of the conductivity matrix, solving the pressure field, and calculating yield.

[0090] (v) Based on the aforementioned multi-level grid encryption system, the present invention can solve for CI and conductivity matrix, and then solve for pressure and output. Figure 6 As can be clearly seen from the pressure distribution of this invention, after level 4 encryption, the fluid flow capacity is enhanced, the pressure propagation speed is accelerated, and the pressure wave range is increased. Figure 7This is a comparison chart of daily and cumulative output for three different grid configuration schemes of the present invention. Figure 7 (a) shows the daily oil production calculated under three mesh configurations: local refinement (using the method of this invention), full refinement (using a fine mesh across the entire domain), and no refinement (using the original coarse mesh). As can be seen from the figure, the initial daily oil production shows a rapid downward trend under all three mesh configurations. The curves for the fully and locally refined mesh configurations almost overlap, while the daily oil production for the no-refinement configuration is significantly lower. This is because the coarse mesh size is too large, failing to accurately capture the pressure changes at the interface between the matrix and the fracture. This results in an underestimation of the connectivity index between the fracture and the matrix, leading to a lower conductivity and ultimately a lower production prediction. Figure 7 (b) shows the cumulative oil production comparison curves. The cumulative production of the fully densified and partially densified schemes is close and significantly higher than that of the undensified scheme, which verifies that the local mesh densification can achieve the same computational accuracy as the global fine mesh, while avoiding the problem that the coarse mesh cannot adequately characterize the flow characteristics near the crack.

[0091] Figure 8 This is a comparison chart of daily and cumulative output under different encryption levels of the present invention. Figure 8 (a) shows the daily oil production curves for different encryption levels. The four curves correspond to encryption levels 4, 3, 2, and 1, respectively. Initially, the four curves showed significant differences, with higher encryption levels resulting in higher initial production. However, the curves gradually converged over time. Figure 8 (b) shows the cumulative oil production comparison curves. Level 4 infill drilling has the highest cumulative production, decreasing sequentially. The difference in cumulative production between Level 4 and Level 3 infill drilling is small, indicating that the Level 3 to Level 4 infill drilling scheme can achieve a balance between accuracy and efficiency while ensuring precision. Based on the above research, this invention follows the above steps to study a fractured horizontal well (e.g., Figure 8 This study analyzes the computational speed and accuracy of output under different grid configuration schemes. Example

[0092] This embodiment uses numerical simulation to demonstrate the effectiveness of the method of the present invention in actual reservoir development through multi-stage fracturing of a horizontal well in a tight oil reservoir. Figure 9As shown, a two-dimensional planar model of a fractured horizontal well in a tight oil reservoir was established. The reservoir length is 800m, width is 400m, oil layer thickness is 10m, matrix permeability is 0.5mD, original formation pressure is 20MPa, and fluid viscosity is 2.0mPa·s. The horizontal well section is 640m long, with 8 fractures fractured, a half-length of 120m, a fracture aperture of 0.001m, and a fracture permeability of 10000mD. The porosity is 0.2, the bottomhole flowing pressure is 15MPa, and the simulation time is 1000 days. The initial mesh consists of 800 uniform 40×20 grids with an initial grid size of 20m×20m. The maximum mesh refinement level is set to 4, and the minimum allowable grid size is 1.25m. The adaptive encryption process was implemented as follows: Level 1 encryption divided 164 grids into four sub-grids, generating 656 10m×10m sub-grids, resulting in 1292 active grids; Level 2 encryption encrypted 324 grids, generating 1296 5m×5m sub-grids, increasing the number of active grids to 2264; Level 3 encryption encrypted 644 grids, generating 2576 2.5m×2.5m sub-grids, increasing the number of active grids to 4196; Level 4 encryption encrypted 1284 grids, generating 5136 1.25m×1.25m sub-grids, resulting in a final number of 8048 active grids. The grid size exhibited a five-layer gradient structure from the wellbore center outwards: 1.25m→2.5m→5m→10m→20m, achieving fine characterization around the fractures and efficient far-field calculations. To verify the effectiveness of the method, three grid schemes were compared, and the results of the comparison of computational efficiency and accuracy are shown in Table 1.

[0093] Table 1 Comparison of computational efficiency and accuracy of different mesh schemes Scheme type Grid configuration Total grid number Simulation time (s) 1000 days cumulative yield (m 3 ) No encryption 20m*20m coarse grid 800 101 4340 Full encryption 1.25m*1.25m full domain uniform fine grid 204800 1162 5720 The method of the present invention 4-level local encryption minimum grid 1.25m 8048 120 5650 As shown in Table 1, compared with the fully encrypted scheme, the method of the present invention reduces the number of grids by 96.07%, shortens the calculation time by 89.7%, and the cumulative production accuracy is only 1.2% lower; while the unencrypted scheme, although having a calculation time of 101s, has a cumulative production accuracy reduced by 24.1%. This embodiment fully verifies the automatic identification capability of the encrypted area, the advantages of improved calculation efficiency, and the accuracy guarantee effect of the method of the present invention, proving that the method has been successfully applied to the numerical simulation of multi-stage fracturing horizontal wells.

[0094] This invention effectively solves the problem of balancing computational accuracy and efficiency in numerical simulation of tight oil reservoirs, particularly addressing issues such as insufficient accuracy of embedded discrete fracture models in handling fracture-dense regions and near-wellbore areas, chaotic multi-scale grid management, and difficulties in calculating connectivity relationships between grids of different scales. It employs a method coupling local grid refinement with embedded discrete fractures, establishing an adaptive refinement criterion based on fracture-grid geometric intersection judgment. Through a triple criterion of level matching, size constraints, and intersection judgment, it achieves automatic identification and progressive refinement of the refinement region, establishing a mathematical model of tight oil reservoirs considering multi-level refinement grids. A tree-like hierarchical structure is used to manage "parent-child" grid relationships, and an improved connectivity index calculation method is introduced, forming a unified conductivity calculation framework for three types of connectivity relationships: fracture-matrix, fracture-fracture, and matrix-matrix, thereby obtaining numerical solutions for pressure and production. Comparison of different grid configuration schemes verifies that the local refinement method maintains computational accuracy comparable to global fine grids while reducing computational costs. The principle is clear and feasible, the method is scientific and reliable, and the process is clear and straightforward. It can further expand the refined numerical simulation method for complex fractured reservoirs and provide important technical support for the efficient development of unconventional oil and gas resources such as deep carbonate reservoirs, tight oil and gas reservoirs, and shale oil and gas.

Claims

1. An adaptive local mesh refinement method based on embedded discrete j unctions, characterized in that The method comprises the following steps: Step one, for the complex fracture system formed by horizontal well fracturing, a two-dimensional tight reservoir physical model is established, and initial grid division is performed, and an EDFM method is used, and fractures and matrix use two independent grid systems; Step two, determine the adaptive local grid refinement area; based on the fracture distribution characteristics and physical flow law, an adaptive local grid refinement criterion is established, and the key area needing grid refinement is automatically identified; Step three, implementation of multi-level grid refinement strategy: before grid refinement, the storage space needs to be expanded; a tree structure is established to manage the multi-level refined grid system, the "parent-child" grid relationship, grid level identification and global connection rules are clear, and the "inactive" grid ID is reserved; Step four, calculate the connectivity index CI after grid refinement: Step five, calculate the conductivity, build the conductivity matrix, the matrix dimension is N total × N total where N total = N m + N f is the sum of the active matrix grid number and the fracture segment number; Step six, solve the pressure and production; the linear equation set of pressure and source and sink items is solved in combination with the boundary conditions to obtain the pressure solution; different boundary conditions are given to realize the pressure solution with fixed production and the production solution with fixed pressure.

2. The method of claim 1, wherein: The step one is specifically: Define the calculation domain parameters; (2) Create a basic grid system; Determination of the number of meshes: , ; wherein, N x , N y are initial x and y direction grid numbers; Δ x 0, Δ y 0are initial x and y direction grid steps; Grid numbering and storage: A row-major natural numbering method is used, with each grid ( i , j The global number of ) is: ; (3) Establish an artificial fracture geometric model; For the fracture system formed by hydraulic fracturing, each fracture is discretized into multiple continuous fracture segments, and a two-level numbering system is established: First level: crack number, ; Second level: Crack segment number, each crack i Discretized into N i s Segment: ; Fracture segment global number: ; where L i f is the length of each crack segment, m; f is the length of each crack segment, m; j is the number of the segment in the i crack.

3. The embedded discrete-junction based self-adaptive local mesh refinement method of claim 2, wherein: The step two is specifically: A step-by-step refinement strategy is adopted, starting from the first level until the maximum refinement level (max_level) is reached, For each refinement level level∈{1,2,...,max_level}, the following steps are performed: (1) Candidate grid screening; for each active cell, three conditions need to be met to be marked for refinement: Condition 1: level matching, only the grid of the previous level can be refined at the current level (grid.level=level−1); Condition 2: Size constraint, min(∆) x , ∆ y )>∆min, where Δmin is the minimum allowable mesh size; Condition 3: fracture intersection judgment; (2) Fracture-grid intersection detection; The intersection detection adopts a line segment-rectangle intersection algorithm, and the judgment logic is as follows: Let the endpoint of the crack segment be ( x 1, y 1) and ( x 2, y 2), the grid boundary is [ xmin , xmax ]×[ ymin , ymax ] (a) Endpoints within the grid: ; (b) Fracture horizontally across the grid: ; (c) Fracture vertically across the grid: ; (d) Fracture obliquely across the grid: For inclined fractures, a parameterized line segment representation and boundary intersection test are used, and the fracture segment parameter equation is: ; Calculate the intersection point parameters of the crack segment and the four edges of the grid respectively t , if there exists t ∈[0, 1] such that the intersection point is on the grid boundary, then it is determined to be intersected; If any of the above conditions is met, it is determined that there is an intersection; (3) Refinement termination; After each level of refinement, it is checked whether there are grids that need to be refined, if there are no grids that need to be refined, the refinement process is terminated, if there are grids that need to be refined, four-part refinement is performed, and the next level is continued.

4. The method of claim 3, wherein: The step three is specifically: (1) Before grid refinement, the storage space needs to be expanded; the number of new grids = the number of grids to be refined × 4, and the coordinates, size, ID and level are expanded; (2) "Parent" grid "deactivation" (active[parent_idx]=0) and "child" grid generation (active[new_idx]=1); for each "parent" grid to be refined, the following steps are performed: Step 1: deactivate the "parent" grid, the "parent" grid no longer participates in the calculation, and its information is reserved for level management; Step 2: Calculate "child" grid geometry parameters; "parent" grid range is [ xmin , xmax ] x [ ymin , ymax Calculate half-size: , ; Step 3: generate 4 "sub" grid coordinates, Subgrid 1 (lower left):[ xmin , xmin +d x ] x [ ymin , ymin +d y ] Subgrid 2 (lower right): xmin +d x , xmax ] x [ ymin , ymin +d y ] Subgrid 3 (top left): xmin , xmin +d x ]×[ ymin +d y , ymax ] Subgrid 4 (top right): xmin + d x , xmax ] x [ ymin + d y , ymax ] (3) "Sub" grid attribute assignment: set the coordinates, size, global ID, level, etc. for each "sub" grid; (4) The encrypted grid hierarchy is a tree hierarchy: After encryption, the total number of active grids: ; For m Level encryption: ; wherein, N initial is the initial number of grids; N refined is the number of encrypted grids; N level1 is the number of grids of the first level of encryption; N level2 is the number of grids of the second level of encryption; (5) Grid number management, new grid ID starts from the current maximum ID + 1 continuous allocation.

5. The embedded discrete-junction based self-adaptive local mesh refinement method of claim 4, wherein: The step four is specifically: (1) Fracture-matrix connectivity index; For the connection between the fracture segment f and the matrix grid m, CI The calculation formula is: ; wherein, CI fm is the connectivity index between the fracture and the matrix grid; l eff is the effective contact length of the fracture segment; ω f is the fracture opening; d cm is the characteristic distance, which is the vertical distance from the center of the fracture segment to the center of the matrix grid; Let the starting coordinate of the crack segment be (x x 1, y 1), the ending coordinate be (x x 2, y 2), the center coordinate be (x x c,f , y c,f ), and the matrix center coordinate be (x x c,m , y c,m ). Then: ; ; (2) Matrix-matrix connectivity index; ; ; wherein, CI mm is an index of connectivity between the matrix and the matrix grid; A mm is an area of a common contact surface of two matrix grids; d ij is an equivalent flow distance; d i is a matrix grid i distance of a center point to a common contact surface; d j is a matrix grid j distance of a center point to a common contact surface; (3) Fracture-fracture intersection detection; Set the coordinates of the start of the crack segment i P1= ( x 1, y 1) , the end coordinates P2= ( x 2, y 2) , the coordinates of the start of the crack segment j Q1= ( x 3, y 3) , the end coordinates Q2= ( x 4, y 4) ; For the crack segment i : ; For the crack segment j : ; Wherein, ; Let , when D =0, the two fracture segments are parallel, when , ; ; When 0 < s <1 and 0< t <1 indicates that two crack segments intersect; where s1, s2are position parameters on the crack segment i , j are direction vectors of the two crack segments; D is the cross product of the two direction vectors; s , t are position parameters on the crack segment i , j , s ∈ [0, 1], t ∈ [0, 1].

6. The method of claim 5, wherein: The step five is specifically: (1) Calculate the fracture-matrix, matrix-matrix, and fracture-fracture conductivities; Fracture-matrix conductivity: ; wherein, T fm Kf is the conductivity between the fracture and the matrix; K m Km is the matrix permeability; K f Kf is the fracture permeability; μ Kf is the fracture permeability; Matrix-matrix conductivity: ; wherein T mm is the conductivity between matrix and matrix. Fracture-fracture conductivity: Conductivity of adjacent segments within the same fracture, for two adjacent fracture segments in the same fracture i and i +1, conductivity reflects one-dimensional linear flow within the fracture: ; wherein, T ff,i,i+1 K is the conductivity between two fracture segments adjacent in the same fracture; h h is the reservoir thickness; Δ x f L is the distance between the centers of two fracture segments. Conductivity of different fracture intersections, when two different fracture segments i and j intersect, fluid exchange at the intersection, conductivity takes the harmonic average: ; Where the unilateral conductivity of each fracture segment is: ; wherein, T ff,i,j is the conductivity of the different fracture intersection segment; α i is the conductivity of the fracture segment i ; α j is the conductivity of the fracture segment j ;0.5⋅ l i represents the distance from the center of the fracture segment to the intersection point; (2) Construct the conductivity matrix.

7. The adaptive local grid refinement method based on embedded discrete fractures according to claim 6, wherein: The method for constructing the conductivity matrix is: First, the diagonal elements, for matrix grid m, the diagonal element is the sum of the time accumulation item and all connection conductivities: ; For the crack segment f (globally numbered N m +f), the diagonal element is: ; wherein, V m V is the volume of the matrix grid; V f Vf is the volume of the fracture segment; φ m φm is the matrix porosity; φ f φf is the fracture porosity; C t,m Cm is the compressibility of the matrix; C t,f Cf is the compressibility of the fracture; Δ t Δt is the time interval; N m (m) is the set of matrix grids adjacent to matrix grid m; N f (m) is the set of fracture segments connected to matrix grid m; N m (f) is the set of matrix grids connected to fracture segment f; N f (f) is the set of fracture segments connected to fracture segment f; Then the off-diagonal elements, for any connected cell i and j The off-diagonal elements are negative values of the conductance: ; The negative sign reflects the physical property that fluid flows from a high pressure unit to a low pressure unit, and the matrix exhibits symmetry when the conductivity satisfies reciprocity T ij = T ji ​ 8. The method of claim 7, wherein: The method for defining the calculation domain parameters is: Define the calculation domain in the two-dimensional plane, considering the complex fracture grid formed by hydraulic fracturing: Computational domain setup: reservoir length L x Determined from well control area and fracture size; reservoir width L y Determined considering the extension of the fracture Boundary condition type: the outer boundary is a closed boundary with no flow; The wellbore position of the inner boundary is treated with an equivalent wellbore radius; Coordinate system establishment: x axis along horizontal well direction, y axis perpendicular to x axis, along main fracture direction.

Images