Simulation method and system for rough fracture invasion-percolation two-phase flow

Abstract

The invention discloses a rough fracture invasion-percolation two-phase flow simulation method, which comprises the following steps of: dividing a rough fracture opening field by adopting a plane grid unit, establishing an invasion-percolation two-phase flow numerical model, and simulating a displacement process as a quasi-static interface propulsion process. The unit with the minimum capillary pressure threshold value in the unit time step is invaded and displaced; a Young-Laplace equation is adopted, the curvature radius in the plane and the curvature radius in the opening direction are calculated, and the capillary pressure threshold value of each unit to be displaced is determined; a binary tree data structure is adopted, displacement process information of each unit is recorded, and in a single time step of displacement, the changed binary tree data structure is updated; and in combination with a deep search algorithm and a width search algorithm, unit position information of capillary capture is discriminated and found. According to the method, the influence of the curvature in the plane on the intrusion displacement process is considered, and an efficient data structure and a search algorithm are combined, so that the accurate and efficient simulation of the intrusion-percolation displacement process in the rough fracture is realized.

Description

Technical Field

[0001] The present invention relates to the numerical field of multiphase flow in fractured media, and in particular to a simulation method and system for rough fracture invasion-percolation two-phase flow. Background Art

[0002] Under the influence of long-term geological movements and engineering activities such as hydraulic fracturing, discontinuous structural surfaces such as cracks and joints are commonly developed within rock masses. These structural surfaces often serve as the main channels for fluid flow and material migration within the rock mass. Therefore, revealing the multiphase seepage characteristics of fractured media has become a key research topic in areas such as efficient oil and gas extraction from fractured reservoirs, geological storage of carbon dioxide, deep nuclear waste storage, and contaminant migration and remediation.

[0003] The macroscopic characteristics of gas / liquid-liquid two-phase flow in fractured media are closely related to the morphological features of the displacement interface at the mesoscale. Currently, numerical calculations are often used as an important means to reveal the two-phase flow displacement process and the evolution of interface morphology. Among them, the invasion-percolation two-phase flow model can effectively simulate the two-phase flow displacement process caused by capillary forces under quasi-static conditions. This model assumes that the location with the minimum capillary pressure threshold near the displacement interface is first displaced by invasion. However, previously proposed invasion-percolation two-phase flow models ignore the contribution of in-plane curvature to capillary pressure, which makes the calculated interface morphology more loose than the experimental results. Moreover, these models often incur a huge amount of computational effort when considering capillary capture phenomena and updating displacement interface information, resulting in low computational efficiency. Therefore, it is necessary to improve the traditional invasion-percolation model and propose an efficient numerical simulation method to achieve accurate and rapid simulation of the invasion-percolation two-phase flow displacement process in rough fractures. Summary of the Invention

[0004] To overcome the above-mentioned deficiencies of the prior art, the present invention considers the influence of in-plane curvature on capillary pressure and the capillary capture phenomenon, and proposes a fast and accurate simulation method and system for rough fracture invasion-percolation two-phase flow.

[0005] According to one aspect of the present invention, a method for simulating invasion-percolation two-phase flow in rough fractures is provided. The method comprises: dividing the rough fracture aperture field using plane grid units, establishing an invasion-percolation two-phase flow numerical model, simulating the displacement process as a quasi-static interface advancement process, and displacing the unit with the minimum capillary pressure threshold within a unit time step by invasion; determining the capillary pressure threshold of each unit to be displaced by calculating the in-plane curvature radius and the curvature radius in the aperture direction using the Young-Laplace equation;

[0006] A binary tree data structure is used to record the displacement process information of each unit, and the changing binary tree data structure is updated within a single time step of displacement; the depth search algorithm and the width search algorithm are combined to identify and find the unit position information of capillary capture.

[0007] Furthermore, the Young-Laplace equation is used to determine the capillary pressure threshold of each unit to be displaced by calculating the in-plane curvature radius and the curvature radius in the opening direction, including: using a boundary tracking algorithm to track the coordinates of several adjacent interface units of each unit to be displaced; fitting the coordinates of the tracked several interface units to the curvature circle equation to determine a fitting coefficient based on the curvature circle equation; calculating the corresponding in-plane curvature radius based on the fitting coefficient, calculating the curvature radius in the opening direction based on the static contact angle and the opening of each unit to be displaced, and calculating the capillary pressure threshold of each unit to be displaced in combination with the Young-Laplace equation.

[0008] Furthermore, a binary tree data structure is used to record the displacement process information of each unit, and the changing binary tree data structure is updated within a single time step of displacement, including: defining a one-dimensional array of unit position, invasion status, and capillary pressure threshold; establishing an initial tree according to the layout principle of the binary tree data structure; finding the unit to be displaced with the minimum capillary pressure threshold for invasion within a single time step of displacement; deleting the nodes corresponding to the invaded units and adding the nodes corresponding to the new units to be displaced, and updating the binary tree data structure.

[0009] Furthermore, a one-dimensional array of cell position, invasion state, and capillary pressure threshold is defined, including: converting the two-dimensional position coordinates of each cell into a one-dimensional index; marking the invasion state for each cell and storing it in a one-dimensional array; calculating the capillary pressure threshold of each cell to be displaced and storing it in another dimensional array.

[0010] Furthermore, according to the arrangement principle of the binary tree data structure, an initial tree is established, including: sorting all the units to be displaced from small to large according to the capillary pressure threshold value; inserting all the units to be displaced into the binary tree according to the capillary pressure threshold value; the arrangement principle of the binary tree data structure is that the capillary pressure threshold value of each node in the binary tree is greater than the capillary pressure threshold value of the parent node and less than the capillary pressure threshold value of the child node.

[0011] Furthermore, the nodes corresponding to the invaded units are deleted, and the nodes corresponding to the units to be displaced are added, and the binary tree data structure is updated, including: deleting the original root node from the binary tree; checking the two child nodes corresponding to the root node, and selecting the child node with a smaller capillary pressure threshold as the new root node; updating the positions of each parent node and child node until the last node has no child nodes; moving the end node to the position where the child node is missing due to the node update at the front end; checking whether the updated binary tree meets the layout principle, and if not, swapping the positions; adding the new node corresponding to the unit to be displaced, and checking again whether the updated binary tree meets the layout principle, and if not, swapping the positions.

[0012] Furthermore, the rough fracture aperture field is divided using plane grid units, including: performing plane grid division on the rough fracture aperture field according to the measured rough fracture structural characteristics, so as to generalize the entire rough fracture into a collection of several tiny local flat plate units; the local flat plate unit is a parallelepiped, the side length of the local flat plate unit is determined by the measurement accuracy of the fracture aperture field, and the thickness of the local flat plate unit is equal to the local average aperture of the fracture at the location.

[0013] Furthermore, the depth search algorithm and the width search algorithm are combined to identify and find the unit position information of capillary capture, including: obtaining the position information of the interface unit; based on the position information of the interface unit, using the depth search algorithm to search for the area where capillary capture does not occur; excluding the area where capillary capture does not occur, using the width search algorithm to search for the unit position information where capillary capture occurs.

[0014] According to one aspect of the present invention, the present invention provides a simulation system for rough fracture invasion-percolation two-phase flow, comprising: a grid division module, which uses plane grid units to divide the rough fracture aperture field, establishes an invasion-percolation two-phase flow numerical model, simulates the displacement process as a quasi-static interface advancement process, and the unit with the minimum capillary pressure threshold is invaded and displaced within a unit time step; a threshold calculation module, which uses the Young-Laplace equation to determine the capillary pressure threshold of each unit to be displaced by calculating the in-plane curvature radius and the curvature radius in the opening direction; a data management module, which uses a binary tree data structure to record the displacement process information of each unit and updates the changed binary tree data structure within a single time step of the displacement; and a capillary capture module, which combines a depth search algorithm and a width search algorithm to identify and find the unit position information of the capillary capture.

[0015] According to one aspect of the present disclosure, a non-transitory computer-readable storage medium is provided. The non-transitory computer-readable storage medium stores computer instructions, which enable the computer to execute the method for simulating rough fracture invasion and percolation two-phase flow.

[0016] The above technical solution proposes a rough fracture invasion and percolation two-phase flow simulation method that takes into account in-plane curvature and capillary capture. That is, an in-plane curvature calculation method is proposed in the capillary pressure calculation, and a binary tree data structure is used to realize the recording and rapid update of displacement process information. The depth search algorithm and the width search algorithm are combined to accelerate the judgment of capillary capture, thereby realizing efficient and accurate simulation of the rough fracture invasion and percolation displacement process.

[0017] Compared with the prior art, the present invention has the following beneficial effects:

[0018] (1) This paper proposes a method for calculating the in-plane curvature of the fracture displacement interface, and further reveals the contribution of in-plane curvature to capillary pressure, thereby achieving an accurate solution to the capillary pressure threshold and an accurate simulation of the invasion-percolation displacement process.

[0019] (2) The present invention adopts a binary tree data structure and records the displacement process information by defining a one-dimensional array of unit position, invasion state and capillary pressure threshold. Within a single time step of displacement, the characteristics of the binary tree data structure are utilized to quickly find the unit with the minimum capillary pressure threshold, and to quickly delete the nodes corresponding to the invaded unit and quickly add the nodes corresponding to the unit to be displaced, thereby quickly updating the displacement process information and greatly improving the calculation speed of the invasion-percolation numerical model.

[0020] (3) The present invention combines a depth search algorithm and a width search algorithm to accelerate the identification of capillary capture. That is, a fast depth search algorithm is first used to exclude areas where no capillary capture occurs, and then a width search algorithm is used in the remaining areas to identify the occurrence and specific location of capillary capture. BRIEF DESCRIPTION OF THE DRAWINGS

[0021] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, a brief introduction will be given below to the drawings used 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.

[0022] Figure 1 A flowchart of a method for simulating rough fracture invasion and percolation two-phase flow provided in an embodiment of the present invention.

[0023] Figure 2 This is a first sub-flowchart of the method for simulating rough fracture invasion-percolation two-phase flow provided by an embodiment of the present invention.

[0024] Figure 3This is a second sub-flowchart of the method for simulating rough fracture invasion-percolation two-phase flow provided by an embodiment of the present invention.

[0025] Figure 4 This is a third sub-flowchart of the method for simulating rough fracture invasion-percolation two-phase flow provided by an embodiment of the present invention.

[0026] Figure 5 This is a calculation flow chart of a method for simulating rough fracture invasion and percolation two-phase flow provided by an embodiment of the present invention.

[0027] Figure 6 Schematic diagram of plane mesh division of rough cracks provided by an embodiment of the present invention.

[0028] Figure 7 Schematic diagram of the opening direction curvature and the in-plane curvature provided by an embodiment of the present invention.

[0029] FIG8( a ) is a diagram showing the distribution of rough crack openings according to an embodiment of the present invention.

[0030] FIG8( b ) is a diagram of the intrusion process provided by an embodiment of the present invention.

[0031] FIG8( c ) is a diagram showing the final intrusion result provided by an embodiment of the present invention.

[0032] Figure 9 A schematic diagram of a method for calculating in-plane curvature provided by an embodiment of the present invention.

[0033] FIG10( a ) is a schematic diagram illustrating the association between unit positions and intrusion states according to an embodiment of the present invention.

[0034] FIG10( b ) is a schematic diagram illustrating the association between the capillary pressure threshold and the binary tree data structure provided by an embodiment of the present invention.

[0035] FIG11( a ) is a schematic diagram of a first update of a binary tree data structure provided by an embodiment of the present invention.

[0036] FIG11( b ) is a schematic diagram of a second update of the binary tree data structure provided by an embodiment of the present invention.

[0037] FIG11( c ) is a schematic diagram of a third update of the binary tree data structure provided by an embodiment of the present invention.

[0038] Figure 12 Schematic diagram of combining a depth search algorithm and a width search algorithm for determining capillary capture according to an embodiment of the present invention.

[0039] Figure 13 Schematic diagram of a simulation system for rough fracture intrusion-percolation two-phase flow provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0040] It should be noted that:

[0041] The terms "including" and "having" and any variations thereof in the description and claims of the present invention and the above-mentioned drawings are intended to cover non-exclusive inclusions, for example, a process, method, system, product or apparatus that includes a series of steps or units is not necessarily limited to the steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to these processes, methods, products or apparatuses.

[0042] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are 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 work are within the scope of protection of the present invention. In addition, the technical features in the various embodiments or single embodiments provided by the present invention are arbitrarily combined with each other to form a new technical solution. This combination is not restricted by the sequence of steps and / or structural composition mode, but must be based on the ability of ordinary technicians in this field to implement it. When the combination of technical solutions is contradictory or cannot be implemented, it should be deemed that this combination of technical solutions does not exist and is not within the scope of protection required by the present invention.

[0043] Please refer to the attached Figure 1 and Figure 5 The present invention provides a simulation method for rough fracture invasion-percolation two-phase flow (including steps S101-S107).

[0044] In step S101, the rough fracture aperture field is divided using plane grid units, and an invasion-percolation two-phase flow numerical model is established. The displacement process is simulated as a quasi-static interface advancement process, and the unit with the minimum capillary pressure threshold in a unit time step is displaced by invasion.

[0045] In step S101, the known rough fracture aperture field is divided into a plane grid according to the measured rough fracture structural characteristics, and the entire rough fracture is generalized into a grid composed of several small local flat plate units (i.e. parallelepipeds, as shown in the attached figure). Figure 6The local flat cells are composed of several elements (shown as ). These elements include invaded cells, displaced cells, and uninvaded cells. Invaded cells represent local flat cells in the fracture that have been completely occupied by the invading phase (i.e., the displacing phase, such as the injected fluid). Undisplaced cells are located at the front of the displacement interface and are local flat cells that may soon be invaded by the invading phase. They are also adjacent to invaded cells. Uninvaded cells represent local flat cells in the fracture that are still completely occupied by the displaced phase (i.e., the original fluid). Uninvaded cells are not adjacent to invaded cells. The displacement interface is essentially the contact surface between two fluids (i.e., the invading phase and the displaced phase). Interface cells are defined as invaded cells at the displacement interface. The invading and displaced phases are either gases or liquids. For example, if the displaced phase is crude oil and the invading phase is CO2, CO2 injection can be used to displace the reservoir crude oil to improve the efficiency of oil and gas field development.

[0046] It should be noted that the rough crack aperture refers to the vertical distance between the two walls of the rough crack, and the aperture of each point of the rough crack is different. The rough crack aperture field refers to the aperture distribution of each point in the plane space. Therefore, after the rough crack aperture field is plane meshed, each local flat plate unit corresponds to an aperture. In this embodiment, the thickness of each local flat plate unit is equal to the local average aperture of the crack at its location. In addition, the side length of each local flat plate unit in this embodiment is determined by the measurement accuracy of the rough crack aperture field, so as to avoid invalid mesh refinement while meeting the accuracy requirements and balance the simulation efficiency and reliability.

[0047] Furthermore, after constructing the invasion-percolation two-phase flow numerical model, it is also necessary to initialize the basic parameters, convergence conditions, and interface advancement process of the invasion-percolation two-phase flow numerical model according to the simulation conditions. The core idea of ​​the invasion-percolation two-phase flow numerical model is to study the invasion law of the invading phase in the displaced phase through a displacement process dominated by capillary forces. The simulation conditions are the initial conditions that need to be set, which are determined according to the properties of the fluid and the accuracy of the simulation and are not limited here. The basic parameters include but are not limited to the contact angle, the number of adjacent interface units, and the total volume of the invading phase. The contact angle is usually in the range of 0°~90°, and the number of adjacent interface units is usually set between 3 and 10, both of which are constants. The total volume of the invading phase refers to the total volume of the invading phase that needs to be injected into or occupy the rough fracture space during the displacement process.

[0048] In this embodiment, the entire displacement process is simulated as a quasi-static gas / liquid-liquid interface advancement process, in which the interface advancement is the complete invasion of the unit within each time step. The boundary conditions include the inlet boundary condition, the outlet boundary condition, and the impermeable boundary conditions on both sides. This embodiment also defines the initial interface as the inlet boundary, and the inlet boundary refers to the entry boundary of the invading phase fluid (gas / liquid). The outlet boundary refers to the outflow boundary of the fluid (gas / liquid). The impermeable boundary conditions on both sides refer to the constraints on the remaining two side boundaries of the simulation area, so that the fluids of the invading phase and the displaced phase cannot pass through. The convergence condition refers to the end condition of the simulation of the invasion-percolation two-phase flow numerical model. In this embodiment, the convergence condition is that the displacement interface breaks through (the invading phase reaches the outlet boundary) or the local flat plate unit volume occupied by the invading phase is equal to the total volume of the invading phase initially set.

[0049] Step S103 : Using the Young-Laplace equation, the capillary pressure threshold of each unit to be displaced is determined by calculating the in-plane curvature radius and the opening direction curvature radius.

[0050] In step S103, the contribution of the in-plane curvature to the capillary pressure is considered to achieve an accurate solution to the capillary pressure threshold, thereby achieving an accurate simulation of the invasion-percolation displacement process. Figure 2 , step S103 (including steps S1031 - S1035 ) will be further introduced below.

[0051] Step S1031, using a boundary tracking algorithm to track the coordinates of several adjacent interface units of each unit to be displaced. In this embodiment, the number of adjacent interface units is a parameter set in the initialization of the invasion-percolation two-phase flow numerical model. Figure 9 , select a unit to be displaced (gray diamond point) for which the capillary threshold needs to be calculated, and track the coordinates of n adjacent interface units in the vicinity around this point (black triangle points). This can be achieved through the boundary tracing algorithm in image processing (for example, using the bwtraceboundary function in MATLAB), which is not limited here.

[0052] Step S1033, fitting the coordinates of the tracked adjacent interface units to the curvature circle equation to determine the fitting coefficient according to the curvature circle equation. In this embodiment, the coordinates (x i ,y i ) is fitted to the defined in-plane curvature circle equation using the least squares method, which can be expressed as: , where a1, a2, and a3 are the coefficients to be fitted; n is the number of interface elements to be tracked. The least squares solution of the vector a=[a1, a2, a3] can be obtained using the following formula: .

[0053] Step S1035: Calculate the corresponding in-plane curvature radius based on the fitting coefficient, calculate the opening direction curvature radius based on the static contact angle and the opening of each unit to be displaced, and calculate the capillary pressure threshold of each unit to be displaced using the Young-Laplace equation. In this embodiment, the curvature radius of each unit to be displaced can be calculated based on the a1, a2, and a3 obtained by fitting. The formula is: It should be noted that the sign of the curvature is determined by the position of the in-plane curvature circle. When the in-plane curvature circle is located inside the invading phase, the local interface is convex and the curvature value is positive. When the in-plane curvature circle is located inside the displaced phase, the local interface is concave and the curvature value is negative.

[0054] Please refer to the attached Figure 7 In this embodiment, the capillary pressure threshold P c The capillary pressure threshold is calculated based on the Young-Laplace equation, which depends on the interfacial tension and the curvature of the interface in the opening direction and the in-plane direction: , where γ represents the interfacial tension; k1 and k2 represent the curvature in the opening direction and the in-plane curvature of the unit to be displaced, respectively; r1 and r2 are the corresponding curvature radii in the opening direction and the in-plane curvature radii, respectively; θ c is the static contact angle; b represents the opening of the unit to be displaced. It should be noted that the interfacial tension γ and the static contact angle θ c It is related to the properties of the fluid, so it is a known constant and is the parameter set for the initialization of the numerical model of intrusion-percolation two-phase flow. The opening b is also a known constant after the plane mesh is divided. Therefore, according to the static contact angle θ c The curvature radius r1 in the opening direction is calculated by combining the interfacial tension γ and the in-plane curvature radius r2 calculated in step S1035 to obtain the capillary pressure threshold P. c .

[0055] Step S105 , using a binary tree data structure to record the displacement process information of each unit, and updating the changed binary tree data structure within a single time step of displacement.

[0056] In step S105, a binary tree data structure is used to sort all the units to be displaced according to the capillary pressure threshold value, so as to efficiently manage and quickly retrieve the units to be displaced. Figure 3 , step S105 (including steps S1051 - S1057 ) will be further introduced below.

[0057] Step S1051 defines a one-dimensional array of cell positions, invasion states, and capillary pressure thresholds. In this embodiment, all local flat cell units, except for their boundaries, are connected to four adjacent cells. The sub2ind function in MATLAB is used to convert the two-dimensional position coordinates (x, y) of all local flat cell units into a one-dimensional index i, thereby assigning a unique index i to each local flat cell. Referring to FIG. 10(a), for example, the index corresponding to the local cell (2, 3) in the second row and third column is i = 8. The invasion state of each local cell is labeled and stored in a one-dimensional array s(i). Where s = 2 indicates that the local flat cell is in the invaded state, corresponding to an invaded cell; s = 1 indicates that the local flat cell is in the pending invasion state, corresponding to a pending displacement cell. s = 0 indicates that the local flat cell is in the uninvaded state, corresponding to an uninvaded cell. The capillary pressure threshold of the local flat cell with the invasion state s=1 is calculated by the Young-Laplace equation mentioned above, and the capillary pressure threshold is stored in a one-dimensional array r(i), which represents the difficulty of invasion (the smaller the invasion capillary pressure threshold, the easier the cell is to be invaded).

[0058] Step S1053: Establish an initial tree according to the layout principles of the binary tree data structure. Please refer to FIG10(b). In this embodiment, cells with an intrusion status of s=1 are sorted from smallest to largest according to the value of r(i). Based on the layout principles of the binary tree data structure, cells with s=1 are added to the binary tree to establish the initial tree. The specific layout principles of the binary tree data structure include: ① Each node in the binary tree has a parent node and two child nodes (with the exception that, if the total number of nodes in the binary tree is even, the parent node of the last node in the binary tree has only one child node). ② The intrusion capillary pressure threshold of the cell corresponding to any node in the binary tree is always greater than the intrusion capillary pressure threshold of the parent node and less than the intrusion capillary pressure threshold of the child node, meaning that the parent node is more susceptible to intrusion, which also makes the root node (the node at the top of the tree) most susceptible to intrusion and displacement.

[0059] Step S1055: Within a single displacement time step, a displacement cell with the minimum capillary pressure threshold is searched for invasion. Specifically, within a unit time step, the displacement cell with the lowest current capillary pressure threshold is selected for invasion, and the displacement cell information is updated. In the next time step, a new displacement cell with the lowest capillary pressure threshold is selected for invasion. This invasion process is repeated until the displacement interface is breached (i.e., the cell at the outlet is invaded) or the set total volume of the invading phase is reached. At this point, the numerical model for the invasion-percolation two-phase flow has reached convergence, and the simulation is terminated. Based on the simulation results, the spatiotemporal evolution of the displacement interface in the rough fracture, the distribution characteristics of the capillary capture phase, and macroscopic characteristics such as saturation and interface analysis dimension are analyzed to reveal the fracture invasion-percolation two-phase flow patterns under the influence of different roughness structures (see the rough fracture aperture distribution diagram in Figure 8(a), the invasion process diagram in Figure 8(b), and the final invasion result diagram in Figure 8(c)).

[0060] Step S1057: Delete the node corresponding to the invaded unit, add the node corresponding to the new unit to be displaced, and update the binary tree data structure. Specifically, please refer to Figures 11(a)-(c). The binary tree update process includes: ① Delete the root node (Operation 1 in Figure 11(a)); ② Check the two child nodes corresponding to the root node and select the child node with the smaller capillary pressure threshold as the new root node (Operation 2 in Figure 11(a)); ③ Repeat the same operation and continue to update the parent node and child nodes (Operations 3 and 4 in Figure 11(a)) until the last node has no child nodes; ④ In order to maintain the integrity of the tree, move the last node to the position where the child node is missing due to the node update (Figure 11(a)). ) in the figure; ⑤ Check whether the updated binary tree meets the layout principle, that is, check whether the r(i) value of the child node after the move is greater than that of the parent node. If not, swap the positions (operations 6 and 7 in Figure 11(b)); ⑥ Add a new node corresponding to the unit to be displaced whose invasion state becomes s(i) = 1 (operation 8 in Figure 11(c)); ⑦ Check again whether the updated binary tree meets the layout principle, that is, check whether the r(i) value of the new node after the addition is greater than that of its parent node. If not, swap the positions (operations 9 to 11 in Figure 11(c)).

[0061] In the above embodiment, a binary tree data structure is used to record displacement process information by defining a one-dimensional array of cell locations, invasion states, and capillary pressure thresholds. Within a single displacement time step, the binary tree data structure leverages the characteristics of the structure to rapidly locate cells with the minimum capillary pressure threshold, quickly delete nodes corresponding to invaded cells, and quickly add nodes corresponding to cells to be displaced. This allows for rapid updates of displacement process information, significantly improving the computational speed of the invasion-percolation numerical model.

[0062] Step S107 , combining the depth search algorithm and the width search algorithm to identify and find the unit position information of the capillary capture.

[0063] In step S107, a fast depth search algorithm is first used to exclude areas where no capillary capture occurs, and then a width search algorithm is used in the remaining areas to identify the occurrence and specific location of capillary capture, thereby accelerating the identification of capillary capture. Figure 4 and Figure 12 , step S107 (including steps S1071 - S1075 ) will be further introduced below.

[0064] Step S1071: Obtain the position information of the interface cells. In this embodiment, based on the stored binary tree data, the position information of the most advanced invaded cells in the flow direction (the flow direction refers to the propagation direction of the invading phase, i.e., the x-direction) and the sidewall direction (the sidewall direction refers to the direction perpendicular to the flow direction and pointing toward the fracture boundary, i.e., the y-direction) is extracted. This information, i.e., the position information of the invaded cells at the displacement interface, is stored in a one-dimensional array.

[0065] In step S1073, based on the position information of the interface cells, a depth search algorithm is used to search for areas where capillary trapping has not occurred. Specifically, based on the position of the invaded cells at the displacement interface, a forward search along the flow direction and the sidewall direction is performed, resulting in the gray area shown in Figure 8, which is the area where capillary trapping has not occurred.

[0066] In step S1075, excluding regions without capillary trapping, a width search algorithm is used to search for cell locations where capillary trapping occurs. Specifically, the black region outside the gray region in Figure 8 is searched to determine whether capillary trapping occurs. Capillary trapping is determined by determining whether a simply connected region disconnected from the main flow region exists after the displacement interface. This can be implemented using the bwlabeln and regionprops functions in MATLAB.

[0067] It should be noted that capillary capture refers to the phenomenon in which, when an invading phase displaces a displaced phase, the displaced phase is locally "captured" and isolated within the displacement interface due to the combined effects of capillary forces and interfacial curvature, forming isolated liquid masses or clusters disconnected from the main flow region (i.e., simply connected regions). It can be understood that by precisely locating the capillary capture locations, the invading-percolation two-phase flow simulation method provided by the present invention can more realistically simulate the complex flow in rough fractures.

[0068] Please refer to the attached Figure 13The present invention also provides a simulation system 100 for rough fracture invasion and percolation two-phase flow. The simulation system 100 includes a grid division module 1, a threshold calculation module 2, a data management module 3, and a capillary capture module 4. The grid division module 1 uses plane grid cells to divide the rough fracture aperture field, establishes an invasion and percolation two-phase flow numerical model, and simulates the displacement process as a quasi-static interface advancement process. Within a unit time step, the unit with the minimum capillary pressure threshold is invaded and displaced. The threshold calculation module 2 uses the Young-Laplace equation to determine the capillary pressure threshold of each unit to be displaced by calculating the in-plane curvature radius and the curvature radius in the opening direction. The data management module 3 uses a binary tree data structure to record the displacement process information of each unit and updates the changing binary tree data structure within a single displacement time step. The capillary capture module 4 combines the depth search algorithm and the width search algorithm to identify and find the unit location information of the capillary capture.

[0069] In summary, the present invention discloses an efficient simulation method for invasion-percolation two-phase flow in rough fractures taking into account in-plane curvature and capillary capture. Based on the traditional numerical model of invasion-percolation two-phase flow, this method takes into account the contribution of in-plane curvature to capillary pressure, proposes an in-plane curvature calculation method based on the curvature circle equation and the least squares method, and realizes the accurate solution of the capillary pressure threshold; in the invasion displacement process, a binary tree data structure is used to record and quickly update displacement process information such as the position of the unit, invasion status and capillary pressure threshold; at the same time, combined with the depth search algorithm and the width search algorithm, a rapid judgment of capillary capture is achieved. In other words, the present invention realizes the accurate and efficient simulation of the invasion-percolation displacement process in rough fractures by considering the influence of in-plane curvature on the invasion displacement process and combining efficient data structures and search algorithms, so as to reveal the laws of fracture invasion-percolation two-phase flow under the influence of different rough structures.

[0070] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. 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 replace some or all of the technical features therein with equivalents. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the technical solutions of the embodiments of the present invention.

Claims

1. A method for simulating two-phase flow in rough fracture intrusion and percolation, characterized in that: The simulation method comprises: The rough fracture aperture field is divided using plane grid cells, and a numerical model of invasion-percolation two-phase flow is established. The displacement process is simulated as a quasi-static interface advancement process, and the cell with the minimum capillary pressure threshold within a unit time step is displaced by invasion. The Young-Laplace equation is used to determine the capillary pressure threshold of each unit to be displaced by calculating the in-plane curvature radius and the curvature radius in the opening direction. A binary tree data structure is used to record the displacement process information of each unit, and the changed binary tree data structure is updated within a single time step of displacement; Combine the depth search algorithm and width search algorithm to identify and find the unit position information of capillary capture.

2. A method for simulating two-phase flow in rough fracture intrusion and percolation according to claim 1, characterized in that: The Young-Laplace equation is used to calculate the in-plane curvature radius and the opening direction curvature radius to determine the capillary pressure threshold of each unit to be displaced, including: The boundary tracking algorithm is used to track the coordinates of several adjacent interface units of each unit to be displaced; Fitting the coordinates of the tracked plurality of interface units to a curvature circle equation to determine a fitting coefficient according to the curvature circle equation; The corresponding in-plane curvature radius is calculated based on the fitting coefficient, the opening direction curvature radius is calculated based on the static contact angle and the opening of each unit to be displaced, and the capillary pressure threshold of each unit to be displaced is calculated in combination with the Young-Laplace equation.

3. A method for simulating rough fracture invasion and percolation two-phase flow according to claim 1, characterized in that: A binary tree data structure is used to record the displacement process information of each unit, and the changed binary tree data structure is updated within a single time step of displacement, including: A one-dimensional array defining the cell position, invasion state, and capillary pressure threshold; According to the layout principle of binary tree data structure, establish the initial tree; In a single time step of displacement, find the unit to be displaced with the minimum capillary pressure threshold for invasion; Delete the nodes corresponding to the invaded units, add the nodes corresponding to the new units to be displaced, and update the binary tree data structure.

4. A method for simulating rough fracture invasion and percolation two-phase flow according to claim 3, characterized in that: A one-dimensional array defining the cell position, invasion state, and capillary pressure threshold, including: Convert the two-dimensional position coordinates of each unit into a one-dimensional index; Mark the invasion status of each unit and store it in a one-dimensional array; The capillary pressure threshold of each unit to be displaced is calculated and stored in another dimensional array.

5. A method for simulating rough fracture invasion and percolation two-phase flow according to claim 3, characterized in that: According to the layout principle of the binary tree data structure, the initial tree is established, including: Sort all units to be displaced from small to large according to the capillary pressure threshold; All units to be displaced are inserted into a binary tree according to the capillary pressure threshold value; the arrangement principle of the binary tree data structure is that the capillary pressure threshold value of each node in the binary tree is greater than the capillary pressure threshold value of the parent node and less than the capillary pressure threshold value of the child node.

6. A method for simulating rough fracture invasion and percolation two-phase flow according to claim 3, characterized in that: Delete the nodes corresponding to the invaded units, add the nodes corresponding to the units to be displaced, and update the binary tree data structure, including: Delete the original root node from the binary tree; Check the two child nodes corresponding to the root node and select the child node with the smaller capillary pressure threshold as the new root node; Update the positions of each parent node and child node until the last node has no child nodes; Move the end node to the position where the child node is missing due to node update at the front end; Check whether the updated binary tree meets the layout principle. If not, swap positions. Add the corresponding node of the new unit to be replaced, and check again whether the updated binary tree meets the layout principle. If not, swap the positions.

7. A method for simulating rough fracture invasion and percolation two-phase flow according to claim 1, characterized in that: Plane grid units are used to divide the rough fracture aperture field, including: The rough fracture aperture field is divided into a plane grid according to the measured structural characteristics of the rough fracture, so as to generalize the entire rough fracture into a collection of several tiny local flat plate units; the local flat plate unit is a parallelepiped, the side length of the local flat plate unit is determined by the measurement accuracy of the fracture aperture field, and the thickness of the local flat plate unit is equal to the local average aperture of the fracture at the location.

8. A method for simulating rough fracture invasion and percolation two-phase flow according to claim 1, characterized in that: Combining the depth search algorithm and the width search algorithm, the unit position information of the capillary capture is identified and found, including: Get the position information of the interface unit; Based on the position information of the interface unit, a depth search algorithm is used to search for an area where no capillary capture occurs; The area where no capillary capture occurs is excluded, and a width search algorithm is used to search for unit position information where capillary capture occurs.

9. A simulation system for rough fracture invasion-percolation two-phase flow, characterized in that: include: The grid division module uses plane grid units to divide the rough fracture aperture field, establishes an invasion-percolation two-phase flow numerical model, and simulates the displacement process as a quasi-static interface advancement process. The unit with the minimum capillary pressure threshold in a unit time step is displaced by invasion; The threshold calculation module uses the Young-Laplace equation to determine the capillary pressure threshold of each unit to be displaced by calculating the in-plane curvature radius and the curvature radius in the opening direction; The data management module uses a binary tree data structure to record the displacement process information of each unit and updates the changed binary tree data structure within a single time step of displacement; The capillary capture module combines the depth search algorithm and the width search algorithm to identify and find the unit position information of the capillary capture.

10. A non-transitory computer-readable storage medium, characterized in that The non-transitory computer-readable storage medium stores computer instructions, and the computer instructions enable the computer to execute the simulation method of rough fracture invasion and percolation two-phase flow according to any one of claims 1 to 8.

Images