A numerical discretization method and device for complex fractures in shale reservoir fracturing

Abstract

The present invention provides a method and apparatus for numerically discretizing complex fractures in shale reservoir fracturing. The method comprises loading an original numerical simulation background matrix grid and fracture polygons, parsing the fracture polygons, recording the fracture polygon coordinates, removing duplicate points from the fracture polygon coordinates, and characterizing the fracture polygons as a topologically consistent geometric structure based on the fracture polygon coordinates without duplicate points. Based on the topologically consistent geometric structure, the method calculates the fracture polygons intercepted by the background matrix grid when each fracture polygon is embedded in the background matrix grid. Connected fracture polygons within the background matrix grid are identified and grouped, and the conductivity between the background matrix grid and connected fracture polygons, and between connected fracture polygons and connected fracture polygons, is calculated. When calculating the conductivity parameters of embedded fractures, the method merges and divides fractures based on their connectivity, thereby improving the accuracy and efficiency of the calculation of non-connected fractures within the background grid.

Description

Technical Field

[0001] The present invention belongs to the technical field of oil reservoir description and oil reservoir numerical simulation, and in particular relates to a method and device for numerically discretizing complex fractures in shale reservoir fracturing. Background Art

[0002] Shale matrix permeability is extremely low (hundreds of nanodarcies to microdarcies). Horizontal wells and hydraulic fracturing are the core technologies that enable industrial shale oil and gas production. The fracture network created by hydraulic fracturing provides high-permeability pathways for fluid flow, increasing formation permeability while expanding the contact area between the reservoir and the wellbore, thereby increasing production. Techniques such as staged horizontal well fracturing can further increase reservoir volume and boost shale oil and gas production.

[0003] The morphology of hydraulic fractures in shale reservoirs is more complex. Warpinski et al., through field experiments and microseismic monitoring, found that hydraulic fractures in shale reservoirs are not traditional planar, symmetrical, double-wing fractures, but rather a complex, interlocking fracture system. Blanton, Pater, Gale, Chen Mian et al., through laboratory experiments, found that complex hydraulic fracture networks often form in formations with well-developed natural fractures.

[0004] Thanks to the advancement of fracture extension simulation and monitoring technology, such as Figure 1 As shown in Figure 2, the characterization of shale reservoir fracture networks is becoming increasingly refined. Due to the strong heterogeneity of shale reservoir physical properties and natural fractures, the fracture networks formed by fracturing exhibit differentiated characteristics, with significant differences in fracture morphology and parameters between layers, segments, and clusters.

[0005] In order to more accurately describe the mechanism of complex fractures, discrete fracture models are often used for reservoir numerical simulation. However, existing discrete fracture numerical models are all based on matching grids, that is, the fracture is used as the inner boundary and the grid is divided based on it as the constraint surface. Figure 2 As shown in the figure on the left, a simple fracture is embedded, while the figure on the right shows a complex fracture. Due to the complex fracture geometry, unstructured meshing techniques are required, which is a very complex and tedious process. This is especially true when the distance or angle between fractures is very small, often resulting in computational difficulties due to poor meshing quality. To address this issue, the Lee group and Moinfar et al. proposed an embedded discrete fracture model. This model embeds the fracture network directly into the bedrock structured mesh system, avoiding the complex unstructured meshing process. Because the embedded discrete fracture model does not consider the fracture morphology within the reservoir, it only requires a simple meshing of the bedrock system, significantly reducing meshing complexity and improving computational efficiency. In recent years, it has become increasingly widely accepted and used. Furthermore, this model organically combines existing mature reservoir numerical simulation techniques with the discrete fracture network model, enabling precise simulation of fluid flow in fractured reservoirs.

[0006] Complex fracture network polygons contain a large number of coincident points, which wastes storage space and makes it difficult to process fracture connectivity. It is necessary to identify and remove all the duplicate points and construct a valid point list.

[0007] There are many algorithms for finding duplicate points. The conventional pairwise comparison method is the most direct and effective method. i and point p j The method of whether the overlap is based on the formula dist(p i , p j )=||p i -p j ||<ε, where ε is a decimal (1E-2 in conventional shale reservoirs). When the above formula holds, that is, the distance between the two points is less than ε, then p i and point p j The pairwise comparison method is very accurate and effective for dealing with small-scale problems, but its time complexity is O(n 2 ), it is inefficient for large-scale problems. For example, for a complex crack with 10,000 points, of which 2,000 are random and repeated, it takes 2-5 minutes to determine the repeated points on a typical PC.

[0008] In addition to the pairwise comparison method, the most commonly used duplicate point search algorithm is the search algorithm based on point sorting (such as point moment sorting and three-axis block sorting algorithm). The computational complexity of this type of algorithm is O(nlog2n), which is still not efficient for very large and complex fracture network models.

[0009] However, many technical challenges remain when applying the embedded discrete fracture method to complex fracture networks in shale reservoirs. A major challenge is that complex shale fracture networks are typically composed of numerous small grid cells. Directly applying the embedded discrete fracture method results in a large number of computational nodes and long simulation times. Summary of the Invention

[0010] In response to the above problems, the present invention proposes a method and device for numerically discretizing complex fractures in shale reservoir fracturing.

[0011] In order to achieve the above object, the present invention adopts the following technical solutions:

[0012] A numerical discretization method for complex fractures in shale reservoir fracturing includes the following steps:

[0013] Load the original numerical simulation background matrix grid and fracture polygons, parse the fracture polygons, and record the fracture polygon coordinates;

[0014] Eliminate duplicate points in the crack polygon coordinates;

[0015] The fracture multi-deformation is characterized as a geometric structure with consistent topology according to the fracture polygon coordinates after removing duplicate points;

[0016] Based on the geometric structure of consistent topology, the fracture polygons intercepted by the background matrix grid are calculated when each fracture polygon is embedded in the background matrix grid;

[0017] Identify connected fracture polygons within the background matrix grid and group them;

[0018] The conductivity is calculated between adjacent background matrix grids, between background matrix grids and connected fracture polygons, and between connected fracture polygons.

[0019] Preferably, the analyzing the crack polygon and recording the coordinates of the crack polygon include the following steps:

[0020] Parse the crack FAB file, traverse each crack polygon recorded in it, and record its coordinates;

[0021] All points are grouped into point sets and numbered sequentially, and the point number corresponding to each crack polygon vertex is recorded.

[0022] Preferably, removing duplicate points in the crack polygon comprises the following steps:

[0023] Define the hash function. The hash function of any point P(x, y, z) is:

[0024] h(p)=floor{[D×10 -3 -floor(D×10 -3 )×10 8 ]};

[0025] in, is the Euclidean norm of point P, (x, y, z) are the coordinates of point P;

[0026] Traverse all points and apply the above formula to calculate the hash function value of each point and store it in vector H;

[0027] Define an initially empty set of non-repeated points S, traverse all points again, and for point P, find all points in the point set S with the same hash function value according to its hash function value h(p), and store them in the point set S';

[0028] If the point set S' is an empty set, then insert point P into the point set S. If the point set S' is not empty, then traverse any point P' in the point set S' and apply dist(p i -p j )=||pi -p j The ||<ε formula is used to determine whether P and P' are duplicate points. In the formula, dist represents the distance function, pi and pj represent the two points to be determined as duplicates, and ε represents the distance threshold value, which can be 0.01 for general shale reservoir models.

[0029] If there is no duplicate point, insert point P into the point set S. If it is a duplicate point, skip it until all points in the point set S are traversed to form a point set S without duplicates of the origin point set.

[0030] Preferably, the method further comprises the following steps:

[0031] After finding and removing all duplicate points, regenerate the point set list and record the number of each polygon vertex.

[0032] Preferably, the calculation of the crack polygons intercepted by the background matrix grid when each crack polygon is embedded in the background matrix grid comprises the following steps:

[0033] Construct a corner point grid;

[0034] Determine the position of the corner point grid and obtain the background matrix grid;

[0035] The background matrix grid intercepts the discrete fracture polygons;

[0036] Based on the geometric structure with consistent topology, the fracture polygons intercepted by the background matrix grid are calculated when each fracture polygon is embedded in the background matrix grid.

[0037] Preferably, determining the connected fracture polygons within the background matrix grid and grouping them comprises the following steps:

[0038] Traverse all the fracture polygons in any background matrix grid block and randomly select two fracture polygons for connectivity calculation;

[0039] Calculate the equation of each edge in each crack polygon;

[0040] Calculate whether there is an intersection point for each crack polygon. If so, the two crack polygons intersect. Otherwise, if there is no intersection point on all edges, they are not connected.

[0041] The crack polygons with intersection points are divided into one group, and the crack polygons without intersection points are divided into another group;

[0042] The set of connected fracture polygons is denoted as F;

[0043] F={F1,F2,…F i …F j …F n}, where i, j, and n represent the i-th, j-th, and n-th connected fracture polygons, respectively.

[0044] A numerical discretization device for complex fractures in shale reservoir fracturing, comprising a model loading module, a duplicate point elimination module, a fracture geometry topology characterization module, a discrete fracture embedding calculation module, a fracture connectivity calculation module, and a conductivity calculation module;

[0045] Model loading module, used to load the original numerical simulation background matrix grid and fracture polygons, parse the fracture polygons, and record the fracture polygon coordinates;

[0046] Duplicate point elimination module, used to eliminate duplicate points in the crack polygon coordinates;

[0047] A crack geometry topology characterization module is used to characterize the crack polygon into a geometric structure with consistent topology based on the crack polygon coordinates with duplicate points removed;

[0048] A discrete fracture embedding calculation module is used to calculate the fracture polygons intercepted by the background matrix grid when each fracture polygon is embedded in the background matrix grid based on the geometric structure with consistent topology;

[0049] a fracture connectivity calculation module for determining connected fracture polygons within the background matrix grid and grouping them;

[0050] The conductivity calculation module is used to calculate the conductivity between adjacent background matrix grids, between background matrix grids and connected fracture polygons, and between connected fracture polygons.

[0051] Preferably, the model loading module analyzes the cracks and records the polygonal coordinates of the cracks, including the following steps:

[0052] Parse the crack FAB file, traverse each crack polygon recorded in it, and record its coordinates;

[0053] All points are grouped into point sets and numbered sequentially, and the point number corresponding to each crack polygon vertex is recorded.

[0054] Preferably, the duplicate point elimination module eliminates duplicate points in the crack polygon, comprising the following steps:

[0055] Define the hash function. The hash function of any point P(x, y, z) is:

[0056] h(p)=floor{[D×10 -3 -floor(D×10 -3 )×10 8 ]};

[0057] in, is the Euclidean norm of point P, (x, y, z) are the coordinates of point P;

[0058] Traverse all points and apply the above formula to calculate the hash function value of each point and store it in vector H;

[0059] Define an initially empty set of non-repeated points S, traverse all points again, and for point P, find all points in the point set S with the same hash function value according to its hash function value h(p), and store them in the point set S';

[0060] If the point set S' is an empty set, then insert point P into the point set S. If the point set S' is not empty, then traverse any point P' in the point set S' and apply dist(p i -p j )=||p i -p j The ||<ε formula determines whether P and P' are duplicate points. If there are no duplicate points, point P is inserted into the point set S. If it is a duplicate point, it is skipped until all points in the point set S are traversed to form a non-duplicate point set S of the origin set.

[0061] Preferably, after the duplicate point elimination module finds and eliminates all duplicate points, it regenerates the point set list and records the number of each polygon vertex.

[0062] Preferably, the discrete crack embedding calculation module is used to calculate the crack polygons intercepted by the background matrix grid when each crack polygon is embedded in the background matrix grid, comprising the following steps:

[0063] Construct a grid of corner points;

[0064] Determine the position of the corner point grid and obtain the background matrix grid;

[0065] The background matrix grid intercepts the discrete fracture polygons;

[0066] Based on the geometric structure with consistent topology, the fracture polygons intercepted by the background matrix grid are calculated when each fracture polygon is embedded in the background matrix grid.

[0067] Preferably, the fracture connectivity calculation module is used to determine the connected fracture polygons within the background matrix grid and group them, comprising the following steps:

[0068] Traverse all the fracture polygons in any background matrix grid block and randomly select two fracture polygons for connectivity calculation;

[0069] Calculate the equation of each edge in each crack polygon;

[0070] Calculate whether there is an intersection point for each crack polygon. If so, the two crack polygons intersect. Otherwise, if there is no intersection point on all edges, they are not connected.

[0071] The crack polygons with intersection points are divided into one group, and the crack polygons without intersection points are divided into another group;

[0072] The set of connected fracture polygons is denoted as F;

[0073] F={F1,F2,…F i …F j …F n}, where i, j, and n represent the i-th, j-th, and n-th connected fracture polygons, respectively.

[0074] Beneficial effects of the present invention:

[0075] 1. When calculating the embedded crack conduction parameters, the present invention merges and divides the cracks according to their connectivity, which can improve the calculation accuracy and efficiency of non-connected cracks in the background grid.

[0076] 2. The present invention reduces the problems of multiple computing nodes and simulation time field by eliminating duplicate points in the crack polygon coordinates and representing the crack polygon as a consistent topological structure.

[0077] Other features and advantages of the present invention will be described in the following description, and in part will become apparent from the description, or will be understood by practicing the present invention. The purpose and other advantages of the present invention can be realized and obtained by the structures pointed out in the description, claims and drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0078] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following is a brief introduction to the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0079] Figure 1 A schematic diagram of the complex fracture network formed by hydraulic fracturing is shown;

[0080] Figure 2 A schematic diagram of an embedded crack is shown;

[0081] Figure 3 A flowchart of a numerical discretization method for complex fractures in shale reservoir fracturing is shown;

[0082] Figure 4 Shown is a schematic diagram of a hash table data structure;

[0083] Figure 5 The topological characterization of irregular crack geometry is shown;

[0084] Figure 6 A geometric topological representation showing the existence of overlap of common edges;

[0085] Figure 7 A schematic diagram of a corner grid is shown;

[0086] Figure 8 A polygon stitching method is shown;

[0087] Figure 9 The polygons before and after stitching are shown;

[0088] Figure 10 A schematic diagram of the least squares fitting method is shown;

[0089] Figure 11 A geometric characterization result diagram is shown;

[0090] Figure 12 A background matrix grid block diagram is shown;

[0091] Figure 13 A schematic diagram showing the presence of unconnected cracks in the background matrix grid is shown;

[0092] Figure 14 A comparison diagram of oil and water saturation is shown;

[0093] Figure 15 A comparison chart of the results of daily oil production curves is shown;

[0094] Figure 16 The figure shows the comparison between the number of grids and the calculation time. DETAILED DESCRIPTION

[0095] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts shall fall within the scope of protection of the present invention.

[0096] A numerical discretization method for complex fractures in shale reservoirs, such as Figure 3 As shown, the following steps are included:

[0097] Load the original numerical simulation background matrix grid and fracture polygons, parse the fracture polygons, and record the fracture polygon coordinates;

[0098] Eliminate duplicate points in the crack polygon coordinates;

[0099] The fracture multi-deformation is characterized as a geometric structure with consistent topology according to the fracture polygon coordinates after removing duplicate points;

[0100] Based on the geometric structure of consistent topology, the fracture polygons intercepted by the background matrix grid are calculated when each fracture polygon is embedded in the background matrix grid;

[0101] Identify connected fracture polygons within the background matrix grid and group them;

[0102] The conductivity is calculated between adjacent background matrix grids, between background matrix grids and connected fracture polygons, and between connected fracture polygons.

[0103] Furthermore, analyzing the cracks and recording the coordinates of the crack polygons include the following steps:

[0104] Parse the crack FAB file, traverse each crack polygon recorded in it, and record its coordinates;

[0105] All points are grouped into point sets and numbered sequentially, and the point number corresponding to each crack polygon vertex is recorded.

[0106] The details are shown in the following table:

[0107] Table 1: Crack FAB file analysis table

[0108]

[0109] Furthermore, removing duplicate points in the crack polygon includes the following steps:

[0110] Define the hash function. The hash function of any point P(x, y, z) is:

[0111] h(p)=floor{[D×10 -3 -floor(D×10 -3 )×10 8 ]};

[0112] in, is the Euclidean norm of point P, (x, y, z) are the coordinates of point P;

[0113] Traverse all points and apply the above formula to calculate the hash function value of each point and store it in vector H;

[0114] Define an initially empty set of non-repeated points S, traverse all points again, and for point P, find all points in the point set S with the same hash function value according to its hash function value h(p), and store them in the point set S';

[0115] If the point set S' is an empty set, then insert point P into the point set S. If the point set S' is not empty, then traverse any point P' in the point set S' and apply dist(p i -p j )=||p i -p j The ||<ε formula determines whether P and P' are duplicate points. If there are no duplicate points, point P is inserted into the point set S. If it is a duplicate point, it is skipped until all points in the point set S are traversed to form a non-duplicate point set S of the origin set.

[0116] It should be noted that the present invention proposes an efficient algorithm for determining large-scale (more than 10,000 points) random repeated points using a hash table. Figure 4 As shown in the figure, a hash table, also known as a hash table, is a data structure that directly accesses records based on a key value. In other words, it accesses records by mapping a key value to a location in the table, speeding up lookups. This mapping function is called a hash function, and the array that stores the records is called a hash table. Given a table M, there exists a function that, for any given key value, returns the address of the record in the table containing that key. The table M is called a hash table, and the function is called a hash function.

[0117] If the keyword is k, its value is stored in the storage location of f(k). Thus, the searched record can be directly retrieved without comparison. This correspondence f is called a hash function, and the table built based on this concept is called a hash table. Different keywords may obtain the same hash address, that is, k1≠k2, and f(k1)=f(k2). This phenomenon is called a collision. Keywords with the same function value are called synonyms for this hash function. In summary, based on the hash function f(k) and a method for handling collisions, a set of keywords is mapped to a finite, continuous set of addresses (intervals), and the "image" of the keyword in the address set is used as the storage location of the record in the table. This table is called a hash table. This mapping process is called hashing, and the resulting storage location is called a hash address. Hash functions can make access to a data sequence faster and more efficient. Through hash functions, data elements can be located more quickly.

[0118] In actual work, different hash functions need to be used depending on different situations. Commonly used methods are:

[0119] 1. Direct addressing: The hash address is calculated by taking the keyword or a linear function of the keyword. This is H(key) = key or H(key) = a·key + b, where a and b are constants (this type of hash function is called a self-function). If a value already exists in H(key), the next value is searched, and so on until there is no value left in H(key). Then, the value is added.

[0120] 2. Digital analysis: Analyzing a set of data, such as the birth dates of a group of employees, we find that the first few digits of the birth date are generally the same, which increases the probability of a conflict. However, we find that the last few digits of the birth date, representing the month and the specific day, vary greatly. If the latter digits are used to construct the hash address, the probability of a conflict is significantly reduced. Therefore, digital analysis is to identify numerical patterns and use these data to construct hash addresses with a lower probability of conflict.

[0121] 3. Squaring the middle: When it's unclear which bits of a keyword are most evenly distributed, you can first square the keyword and then use the middle bits of the squared value as the hash address. This is because the middle bits of the squared value are correlated with every bit in the keyword, so different keywords will likely generate different hash addresses.

[0122] The search process of the hash table is basically the same as the table creation process. Some key codes can be found directly through the address converted by the hash function, while other key codes conflict with the address obtained by the hash function and need to be searched according to the conflict handling method. Among the three conflict handling methods introduced, the search after the conflict is still the process of comparing the given value with the key code. Therefore, the measure of the hash table search efficiency is still measured by the average search length. During the search process, the number of key code comparisons depends on the number of conflicts. The fewer conflicts generated, the higher the search efficiency. The more conflicts generated, the lower the search efficiency. Therefore, the factors that affect the number of conflicts are also the factors that affect the search efficiency. For the repeated points in the complex crack network, it is proposed that the key to improving the efficiency of the hash algorithm lies in constructing a suitable hash function and reducing the loading factor β of the hash table. The loading factor of the hash table is defined as:

[0123] β = number of elements filled in the table / length of the hash table;

[0124] β is a factor that indicates how full the hash table is. Since the table length is a fixed value, β is proportional to the number of elements in the table. Therefore, the larger β is, the more elements in the table, and the greater the possibility of conflicts. The smaller β is, the fewer elements in the table, and the less likely conflicts are.

[0125] The hash function can significantly reduce hash collisions. Testing has shown that for a complex crack example with 10,000 points, including 2,000 random duplicates, the time required to find duplicate points has been reduced from 2-5 minutes to 5-10 seconds.

[0126] Furthermore, the method further comprises the following steps:

[0127] After finding and removing all duplicate points, regenerate the point set list and record the number of each polygon vertex.

[0128] Furthermore, calculating the fracture polygons intercepted by the background matrix grid when each fracture polygon is embedded in the background matrix grid includes the following steps:

[0129] Construct a corner point grid;

[0130] Determine the position of the corner point grid and obtain the background matrix grid;

[0131] The background matrix grid intercepts the discrete fracture polygons;

[0132] Based on the geometric structure with consistent topology, the fracture polygons intercepted by the background matrix grid are calculated when each fracture polygon is embedded in the background matrix grid.

[0133] It should be noted that if Figure 7 As shown in FIG, constructing a corner grid requires parsing the geometric information data of the corner grid file. The corner grid file consists of major parts such as SPECGRID, COORD, ZCORN, and ACTNUM.

[0134]

[0135]

[0136] First, let's assume the reservoir has N1, N2, and N3 grid cells in the X, Y, and Z directions, respectively. To construct the corner grid, we first need to generate pillars along the depth direction. Two points define a line, and each point corresponds to three coordinate values, so defining a pillar requires six data points.

[0137] Even after the corner grid's pillars are fixed, the grid's exact position is still uncertain. Therefore, we need another parameter: the Z coordinates of the grid's eight corner points. Because the grid's corner points can only move along the defined pillars, knowing the coordinates of a point in any direction allows us to determine its position on the pillars—the three-dimensional coordinates of the corner point. Since the contact surfaces of adjacent grids in a corner grid don't necessarily overlap, each corner point of each grid needs to be defined individually. For a hexahedral grid, each grid has eight corner points, each requiring a coordinate value for location. The total number of Z coordinates we need is 8*N1*N2*N3.

[0138] Then, the polygons intercepted by the background matrix grid are calculated when each discrete crack polygon is embedded in the background matrix grid. In this step, the method calculates the polygons intercepted by the background matrix grid for each crack polygon based on the 8 vertex information of the background matrix grid (hexahedron). In specific implementation, the intersection of each edge (line segment) of the crack polygon and the 6 faces of the grid (quadrilateral q) is calculated, and the intersection is surrounded to form the intercepted polygon. The intersection of line segment e and the crack polygon is calculated as follows:

[0139] First, determine whether the line l on which the line segment e lies is parallel to the plane on which the quadrilateral q lies. If the line is not parallel to the plane, there will be an intersection;

[0140] It is known that the straight line l passes through the point m and the direction vector is v1, then its equation is:

[0141] pm-v1·t=0;

[0142] The plane passes through point n and its normal direction vector is v n , its equation can be written as:

[0143] v n (pn) = 0;

[0144] Combining the above two formulas, we can find the coordinates of the intersection point p0; to determine whether it is inside the line segment e(e1, e2), the judgment formula is:

[0145] (p0-e1)·(p0-e2)<0;

[0146] Determine whether it is inside the quadrilateral q(q1,q2,q3,q4). The judgment formula is that for any i∈[1,4], it satisfies:

[0147] [(q i+1 -q i )×(p0-q i )]·[(q i+2 -q i+1 )×(p0-q i+1 )]>0;

[0148] If the intersection point is both inside the line segment e and the quadrilateral q, then it is considered an intersection point. Otherwise, it is considered an invalid intersection point and is discarded. The intercepted polygons come in many shapes, mainly triangles, quadrilaterals, pentagons, and hexagons.

[0149] like Figure 8As shown in the figure, since irregular complex fracture networks are composed of many small polygons, embedding them may result in a surface composed of multiple continuous polygons, increasing the number of computational nodes and the complexity of the conduction relationship in the embedded discrete fracture model. Therefore, the key and difficulty in embedding irregular complex fracture networks lies in splicing the polygons formed by embedding and intercepting the grid blocks at each corner point to reduce the geometric complexity of the final embedded discrete fracture.

[0150] like Figure 9 As shown in the figure, due to the use of a geometric representation of the fracture grid with consistent topology, the process of splicing polygons is relatively convenient: (1) Search the edges of all polygons in a single grid, find the common edges and mark them. The search for common edges can be done through the point set hashing algorithm, or the edge hashing algorithm can be reconstructed, which will not be described here; (2) Delete the marked edges (common edges) and only keep the unmarked edges; (3) Connect the retained unmarked edges end to end to form a spliced ​​polygon.

[0151] Note that the polygons to be joined here may not be coplanar, so plane fitting is required. The fitting algorithm uses the least squares method for fitting. For the above surface, Figure 10 As shown, the fitted plane is a plane passing through point p(x0, y0, z0) and with a normal vector of (a, b, c). The plane equation is:

[0152] a(x-x0)+b(y-y0)+c(z-z0)=0;

[0153] The fitting process minimizes the offset of all points. For n sets of observation data p i =(x i ,y i , z i )(i=1,2,3,…,n), then its least squares objective function is:

[0154]

[0155] Among them, f(xi,u) 2 =[a(xi-x0)+b(yi-y0)+c(zi-z0)] 2 represents the deviation of the i-th data point, ω i is the weight of the i-th data point, which can be taken as a uniform weight or a weighted fit can be performed according to the value of the data point to solve the extreme value problem and obtain the fitting solution.

[0156] So far, if Figure 11 As shown in Figure 3, we obtained the geometric representation of irregular complex fractures for numerical simulation of shale oil and gas reservoirs, which can be used as input data for the embedded discrete fracture numerical simulation model.

[0157] It's important to note that topology is the study of properties of geometric figures or spaces that remain constant despite continuous changes in shape. It considers only the positional relationships between objects, regardless of their shape or size. Topology, literally translated as "geography," originally referred to the study of topography and landforms. Geometric topology is a branch of mathematics that emerged in the 19th century and falls within the realm of geometry. Some aspects of topology appeared as early as the 18th century. Isolated problems discovered at that time played a crucial role in the later development of topology.

[0158] Induced topology refers to a method of constructing a topology, while uniform topology refers to a topology induced by a uniform structure. The mathematical definition of topology is as follows: Let X be a non-empty set, and a subset of X's power set (i.e., a family of sets consisting of certain subsets of X) T is called a topology of X. A topology is a topology if and only if: (1) X and the empty set {} both belong to T; (2) the union of any number of members in T is still in T; and (3) the intersection of a finite number of members in T is still in T. The set X together with its topology T is called a topological space, denoted as (X, T), and the members of T are called open sets of this topological space. The three conditions in the definition are called topological axioms. Generally speaking, many different topologies can be specified on a set, so when talking about a topological space, both the set and the specified topology should be specified. When there is no misunderstanding, a set is often used to refer to a topological space, such as topological space X, topological space Y, etc. At the same time, in the category of topology, we discuss continuous mappings. The definition is: f:(X, T1)→(Y, T2) is continuous if and only if the preimage of the open set is open. Two topological spaces are homeomorphic if and only if there exists a one-to-one corresponding continuous mapping that is the inverse of each other. The homotopy of a map and the homotopy of a space are also very useful definitions.

[0159] like Figure 5 As shown, a uniform topology is a topology induced by a uniform structure. Let (X, U) be a uniform space, T be a subset of X, satisfying: for any x∈T, there exists u∈U such that where U(x) = {y(x, y)∈U}. The set T of all such Ts is a topology on X, called a topology induced by a consistent structure U or a consistent topology. When the topological space X induced by a consistent structure U is compact, then the consistent structure consistent with the topology of X is uniquely determined. A consistent space is completely regular, and completely regular spaces are topologically consistent with it. This leads to the following result: a necessary and sufficient condition for a topology T on a set X to be a consistent topology of a consistent structure on X is that (X, T) is a completely regular space.

[0160] Returning to the geometric representation of irregular complex crack networks, e.g. Figure 6 As shown in Figure 2, the process of constructing a consistent topological representation of the crack is as follows:

[0161] On the original fracture network, the vertices of each fracture polygon are numbered sequentially and recorded as P1, P2, ..., P n Then a certain crack polygon can be expressed as F i ={P fi1 ,P fi2 ,...,P fin}, where P fi1 ,P fi2 ,...,P fin Number all vertices of the polygon;

[0162] After finding and eliminating all duplicate points, record the number correspondence between the original point set and the non-duplicate point set according to the regenerated non-duplicate point set, and update the number F' of each polygon vertex on this basis. i ={P f,i1 ,P f,i2 ,...,P f,in}, where P f,i1 ,P f,i2 ,...,P f,in The numbers of all vertices of the polygon in the set of non-duplicate points;

[0163] Processing of overlapping common edges is divided into:

[0164] 1. For all crack edges e(P e1 , P e2 ), find its unit direction vector v e ;

[0165] in,

[0166] Find all units with the same direction vector v l Since the representation of vectors and points in three-dimensional space is the same (xyz coordinates), the hash algorithm can be used to obtain them.

[0167] 2. For any two polygon sides e1 and e2 with the same unit direction vector, find their overlapping parts. Specifically, use the judgment formula for whether the point falls within the line segment:

[0168] (P * -P e1 )·(P * -P e2 )<0;

[0169] Where Pe1 and Pe2 are points on edges e1 and e2 respectively.

[0170] 3. Find the length of the overlapping part. The formula is:

[0171] L=||P eji -P ej2 ||;

[0172] Where Pej1 and Pej2 are the two vertices of the overlapping line segment.

[0173] Furthermore, connected fracture polygons within the background matrix grid are determined and grouped, including the following steps:

[0174] Traverse all the fracture polygons in any background matrix grid block and randomly select two fracture polygons for connectivity calculation;

[0175] Calculate the equation of each edge in each crack polygon;

[0176] Calculate whether there is an intersection point for each crack polygon. If so, the two crack polygons intersect. Otherwise, if there is no intersection point on all edges, they are not connected.

[0177] The crack polygons with intersection points are divided into one group, and the crack polygons without intersection points are divided into another group;

[0178] The set of connected fracture polygons is denoted as F;

[0179] F={F1,F2,…F i …F j …F n}, where i, j, and n represent the i-th, j-th, and n-th connected fracture polygons, respectively.

[0180] It should be noted that if Figure 12 As shown in the figure, the background matrix grid block is denoted as G. For any pair of fracture polygons f1 and f2 or f3 and f4, calculate whether they are connected: group the fracture polygons according to the connectivity relationship. If there are n groups of connected fractures, then divide the set F of connected fracture polygons into F1, F2, ..., F n Next, we need to calculate the background matrix grid M and the crack F i , any connected crack F i and crack F j The conductivity of the background matrix grid M and the crack F i The conductivity calculation formula between the matrix grid M and the fracture Fi is the same as that of the embedded discrete fracture model proposed by Moinfar (2013), that is, the conductivity calculation formula between the matrix grid M and the fracture Fi is as follows:

[0181]

[0182] in, k M is the permeability of the bedrock grid M, k F is the permeability of the fracture grid F;

[0183] Refers to the equivalent distance between the bedrock grid M and the fracture grid F, that is, the average vertical distance from all points in the bedrock grid to the fracture segment;

[0184] A MF Refers to the contact area between the crack and the bedrock.

[0185] Any connected crack f inside the grid block i and crack f j The conductivity is calculated using the Karimi-Fard (2004) method, as follows:

[0186]

[0187] Among them, a i = kA / D is the geometric conductivity factor of fracture i, k is the fracture permeability, A is the fracture cross-sectional area, and D is the distance from the center point of the fracture polygon to the common edge of the fracture.

[0188] Using the mass conservation governing equation for micro-compressible flow, when solving the initial local pressure field, the linear-no-flow boundary condition is used (but not limited to) as the boundary condition. After solving the governing equation and obtaining the local pressure field, the target conductivity T can be calculated using the following formula. ij Well Index

[0189]

[0190]

[0191] Among them, q IJ is the total flow from coarse grid I to coarse grid J, p I and p J are the average pressures of coarse grids I and J, respectively. i and j represent the fine grid numbers, I and J represent the coarse grid numbers, and the coarse grid refers to the merged fracture grid.

[0192]

[0193]

[0194] Where p represents the fracture pressure, v represents the fracture volume, q is the flow rate between fractures, i and j represent the fine grid numbers, and I and J represent the coarse grid numbers.

[0195] A numerical discretization device for complex fractures in shale reservoir fracturing, comprising a model loading module, a duplicate point elimination module, a fracture geometry topology characterization module, a discrete fracture embedding calculation module, a fracture connectivity calculation module, and a conductivity calculation module;

[0196] Model loading module, used to load the original numerical simulation background matrix grid and fracture polygons, parse the fracture polygons, and record the fracture polygon coordinates;

[0197] Duplicate point elimination module, used to eliminate duplicate points in the crack polygon coordinates;

[0198] A crack geometry topology characterization module is used to characterize the crack polygon into a geometric structure with consistent topology based on the crack polygon coordinates with duplicate points removed;

[0199] A discrete crack embedding calculation module is used to calculate the polygons intercepted by the background matrix grid when each crack polygon is embedded in the background matrix grid based on the geometric structure with consistent topology;

[0200] a fracture connectivity calculation module for determining connected fracture polygons within the background matrix grid and grouping them;

[0201] The conductivity calculation module is used to calculate the conductivity between adjacent background matrix grids, between background matrix grids and connected fracture polygons, and between connected fracture polygons.

[0202] It should be noted that the following Figure 13 (Left) Discrete fracture model, Figure 13 (Right) shows a grid area containing unconnected cracks, and exemplifies a comparison between a discrete crack simulation model established using the results of calculations using the conductivity calculation method for discrete cracks provided in the present disclosure and a discrete crack simulation model established using the results of calculations using the conductive connectivity relationship of the dual-medium model in the related art.

[0203] For example, the area of ​​the discrete fracture simulation model is 7000×7000 ft 2 , including 15 fractures and two wells, one of which is an injection well and the other a production well. The 3D geological grid established by this discrete fracture simulation model has a total of 85×75×5=31,875 grid cells. To compare the results of the methods provided in the embodiments of this disclosure with those in related art, a discrete fracture simulation model based on an unstructured grid (44,773 grid cells) was also established as a baseline model.

[0204] Some parameters of the discrete crack model are shown in Table 2 below:

[0205] Table 2: Discrete crack model parameters

[0206]

[0207] Figure 14 The left figure is the oil-water saturation distribution diagram of the discrete fracture simulation model based on unstructured grid. Figure 14 The middle figure is the oil-water saturation distribution diagram of the discrete fracture simulation model established by calculating the conductive connectivity relationship of the dual medium model. Figure 14 The figure on the right is the oil-water saturation distribution diagram of the discrete fracture simulation model provided by the embodiment of the present disclosure. The oil-water saturation distribution of the discrete fracture simulation model based on an unstructured grid is closest to the oil-water saturation distribution in actual fractured oil and gas reservoirs, where light-colored areas represent areas with higher water saturation, and dark-colored areas represent areas with higher oil saturation. This oil-water saturation distribution diagram is used as a reference. After calculating the conductivity of discrete fractures using the discrete fracture conductivity calculation method provided by the embodiment of the present disclosure, the oil-water saturation distribution diagram corresponding to the discrete fracture simulation model established using this conductivity value is compared to the oil-water saturation distribution diagram corresponding to the discrete fracture simulation model established using the conductivity value calculated using the conductive connectivity relationship of the dual-medium model in the related art. The oil-water saturation distribution of the discrete fracture simulation model based on the unstructured grid is closest to the oil-water saturation distribution in actual fractured oil and gas reservoirs, and has higher prediction accuracy.

[0208] Figure 15 This is a comparison chart of daily oil production curves provided by an embodiment of the present disclosure. The solid line curve in the figure represents the actual daily oil production curve of the oil and gas reservoir and serves as a reference curve. The dashed-dotted line curve represents the daily oil production curve calculated using the discrete fracture conductivity calculation method provided by the embodiment of the present disclosure, and then using the discrete fracture simulation model established using this conductivity value. The dashed-dotted line curve represents the daily oil production curve calculated using the conductivity connectivity relationship of the dual-medium model in the related art, and then using this conductivity value to establish the discrete fracture simulation model. The comparison results show that the trajectory of the dashed-dotted line curve is closer to the trajectory of the actual curve than the dashed-dotted line curve. Because unconnected fractures can be forcibly connected using the methods in the related art and substituted into the calculation, the discrete fracture simulation model established based on the conductivity result has excessively high fracture conductivity, resulting in a rapid decline in daily oil production in the dashed line curve and reduced accuracy.

[0209] Figure 16This patent provides a graph comparing the number of grids and calculation time. The left graph provides a graph comparing the daily oil production and saturation error, provided by an embodiment of the present disclosure. From the comparison of the numerical simulation grid number and calculation time, because the discrete fracture conductivity calculation method provided by the embodiment of the present disclosure carefully considers the unconnected fractures within each grid block, the number of fractures is reduced from 3115 to 1050, and the numerical simulation calculation time is reduced from 2310 seconds to 1033 seconds. However, from the comprehensive comparison of daily oil production and protection degree error, the average daily oil production error is reduced from 5.5% to 1.2%, and the average saturation error is reduced from 9.3% to 1.1%. The discrete fracture conductivity calculation method provided by the embodiment of the present disclosure is a method that achieves a better balance between accuracy and efficiency.

[0210] Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein; and these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A numerical discretization method for complex fractures in shale reservoir fracturing, characterized by: The following steps are involved: Load the original numerical simulation background matrix grid and fracture polygons, parse the fracture polygons, and record the fracture polygon coordinates; Eliminate duplicate points in the crack polygon coordinates, including: Define the hash function. The hash function of any point P (x, y, z) is: ； in, is the Euclidean norm of point P, (x, y, z) is the coordinate of point P; Traverse all points and apply the above formula to calculate the hash function value of each point and store it in vector H; Define an initially empty point set S without duplication, traverse all points again, and for point P, find all points in the point set S with the same hash function value according to its hash function value h(p), and store them in the point set S'; If the point set S' is an empty set, then insert point P into the point set S. If the point set S' is not empty, then traverse any point P' in the point set S' and apply The formula determines whether P and P' are repeated points. In the formula, dist represents the distance function, pi and pj represent the two points to be judged as repeated, and ε represents the distance threshold value, which can be 0.01 for the shale reservoir model. If there is no duplicate point, insert point P into the point set S. If it is a duplicate point, skip it until all points in the point set S are traversed to form a point set S without duplicates in the original point set. The fracture polygon coordinates with duplicate points removed are used to characterize the fracture multi-deformation into a geometric structure with consistent topology, including: On the original fracture network, the vertices of each fracture polygon are numbered sequentially and recorded as P1, P2, ..., P n , then a crack polygon can be expressed as F i ={P fi1 ,P fi2 ,...,P fin }, where P fi1 ,P fi2 ,...,P fin Number all vertices of the crack polygon; After finding and eliminating all duplicate points, record the number correspondence between the original point set and the non-duplicate point set according to the regenerated non-duplicate point set, and update the number F' of each crack polygon vertex on this basis. i ={P f , i1 ,P f , i2 ,...,P f , in }, where P f , i1 ,P f , i2 ,...,P f , in The numbers of all vertices of the Lei Feng polygon in the non-duplicate point set; Processing of overlapping common edges is divided into: For all crack edges , find its unit direction vector ; in, ; Use hashing algorithm to calculate all unit direction vectors with the same edge; For any two polygon edges e1 and e2 with the same unit direction vector, find their overlapping parts: ; Where, P e1 and P e2 are points on edges e1 and e2 respectively; Calculate the length of the overlapping part using the formula: ； Where, P ej1 and P ej2 are two vertices of overlapping line segments; Based on the geometric structure of consistent topology, the fracture polygons intercepted by the background matrix grid are calculated when each fracture polygon is embedded in the background matrix grid; Identify connected fracture polygons within the background matrix grid and group them; The conductivity is calculated between adjacent background matrix grids, between background matrix grids and connected fracture polygons, and between connected fracture polygons.

2. The method for numerical discretization of complex fractures in shale reservoir fracturing according to claim 1, characterized in that: The method of analyzing the crack polygon and recording the coordinates of the crack polygon includes the following steps: Parse the crack FAB file, traverse each crack polygon recorded in it, and record its coordinates; All points are grouped into point sets and numbered sequentially, and the point number corresponding to each crack polygon vertex is recorded.

3. The method for numerical discretization of complex fractures in shale reservoir fracturing according to claim 1, characterized in that: The following steps are also included: After finding and removing all duplicate points, regenerate the point set list and record the number of each polygon vertex.

4. A numerical discretization method for complex fractures in shale reservoir fracturing according to any one of claims 1 to 3, characterized in that: The method of calculating the crack polygons intercepted by the background matrix grid when each crack polygon is embedded in the background matrix grid comprises the following steps: Construct a corner point grid; Determine the position of the corner point grid and obtain the background matrix grid; The background matrix grid intercepts the discrete fracture polygons; Based on the geometric structure with consistent topology, the fracture polygons intercepted by the background matrix grid are calculated when each fracture polygon is embedded in the background matrix grid.

5. The method for numerical discretization of complex fractures in shale reservoir fracturing according to claim 1, characterized in that: The method of determining connected fracture polygons within the background matrix grid and grouping them comprises the following steps: Traverse all the fracture polygons in any background matrix grid block and randomly select two fracture polygons for connectivity calculation; Calculate the equation of each edge in each crack polygon; Calculate whether there is an intersection point for each crack polygon. If so, the two crack polygons intersect. Otherwise, if there is no intersection point on all edges, they are not connected. The crack polygons with intersection points are divided into one group, and the crack polygons without intersection points are divided into another group; The set of connected fracture polygons is denoted as F; , where i, j, and n represent the i-th, j-th, and n-th connected fracture polygons, respectively.

6. A numerical discretization device for complex fractures in shale reservoir fracturing, characterized by: It includes model loading module, duplicate point elimination module, crack geometry topology characterization module, discrete crack embedding calculation module, crack connectivity calculation module and conductivity calculation module; Model loading module, used to load the original numerical simulation background matrix grid and fracture polygons, parse the fracture polygons, and record the fracture polygon coordinates; The duplicate point removal module is used to remove duplicate points in the crack polygon coordinates, including: Define the hash function. The hash function of any point P (x, y, z) is: ； in, is the Euclidean norm of point P, (x, y, z) is the coordinate of point P; Traverse all points and apply the above formula to calculate the hash function value of each point and store it in vector H; Define an initially empty point set S without duplication, traverse all points again, and for point P, find all points in the point set S with the same hash function value according to its hash function value h(p), and store them in the point set S'; If the point set S' is an empty set, then insert point P into the point set S. If the point set S' is not empty, then traverse any point P' in the point set S' and apply The formula determines whether P and P' are duplicate points. If there are no duplicate points, point P is inserted into the point set S. If it is a duplicate point, it is skipped until all points in the point set S are traversed to form a non-duplicate point set S of the origin set. The fracture geometry topology characterization module is used to characterize the fracture polygon into a geometric structure with consistent topology based on the fracture polygon coordinates after removing duplicate points, including: On the original fracture network, the vertices of each fracture polygon are numbered sequentially and recorded as P1, P2, ..., P n , then a crack polygon can be expressed as F i ={P fi1 ,P fi2 ,...,P fin }, where P fi1 ,P fi2 ,...,P fin Number all vertices of the crack polygon; After finding and eliminating all duplicate points, record the number correspondence between the original point set and the non-duplicate point set according to the regenerated non-duplicate point set, and update the number F' of each crack polygon vertex on this basis. i ={P f , i1 ,P f , i2 ,...,P f , in }, where P f , i1 ,P f , i2 ,...,P f , in is the number of all vertices of the Lei Feng polygon in the non-duplicate point set; Processing of overlapping common edges is divided into: For all crack edges , find its unit direction vector ; in, ; Use hashing algorithm to calculate all unit direction vectors with the same edge; For any two polygon edges e1 and e2 with the same unit direction vector, find their overlapping parts: ; Where, P e1 and P e2 are points on edges e1 and e2 respectively; Calculate the length of the overlapping part using the formula: ； Where, P ej1 and P ej2 are two vertices of overlapping line segments; A discrete fracture embedding calculation module is used to calculate the fracture polygons intercepted by the background matrix grid when each fracture polygon is embedded in the background matrix grid based on the geometric structure with consistent topology; a fracture connectivity calculation module for determining connected fracture polygons within the background matrix grid and grouping them; The conductivity calculation module is used to calculate the conductivity between adjacent background matrix grids, between background matrix grids and connected fracture polygons, and between connected fracture polygons.

7. The device for numerically discretizing complex fractures in shale reservoir fracturing according to claim 6, characterized in that: The model loading module analyzes the cracks and records the coordinates of the crack polygons, including the following steps: Parse the crack FAB file, traverse each crack polygon recorded in it, and record its coordinates; All points are grouped into point sets and numbered sequentially, and the point number corresponding to each crack polygon vertex is recorded.

8. The device for numerically discretizing complex fractures in shale reservoir fracturing according to claim 6, characterized in that: After the duplicate point elimination module finds and eliminates all duplicate points, it regenerates the point set list and records the number of each polygon vertex.

9. The device for numerically discretizing complex fractures in shale reservoir fracturing according to claim 6, characterized in that: The discrete crack embedding calculation module is used to calculate the crack polygons intercepted by the background matrix grid when each crack polygon is embedded in the background matrix grid, and includes the following steps: Construct a corner point grid; Determine the position of the corner point grid and obtain the background matrix grid; The background matrix grid intercepts the discrete fracture polygons; Based on the geometric structure with consistent topology, the fracture polygons intercepted by the background matrix grid are calculated when each fracture polygon is embedded in the background matrix grid.

10. The device for numerically discretizing complex fractures in shale reservoir fracturing according to claim 6, characterized in that: The fracture connectivity calculation module is used to determine the connected fracture polygons within the background matrix grid and group them, including the following steps: Traverse all the fracture polygons in any background matrix grid block and randomly select two fracture polygons for connectivity calculation; Calculate the equation of each edge in each crack polygon; Calculate whether there is an intersection point for each crack polygon. If so, the two crack polygons intersect. Otherwise, if there is no intersection point on all edges, they are not connected. The crack polygons with intersection points are divided into one group, and the crack polygons without intersection points are divided into another group; The set of connected fracture polygons is denoted as F; , where i, j, and n represent the i-th, j-th, and n-th connected fracture polygons, respectively.

Images