A simulation method and system of rough fracture invasion-percolation two-phase flow

Abstract

The application discloses a simulation method of rough fissure invasion-percolation two-phase flow, adopts plane grid units to divide a rough fissure opening field, establishes an invasion-percolation two-phase flow numerical model, simulates a displacement process as a quasi-static interface advancing process, and units with a minimum capillary pressure threshold value are invaded and displaced in a unit time step; a Young-Laplace equation is adopted, a curvature radius in a plane and an opening direction curvature radius are calculated, a 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 a changed binary tree data structure is updated in a single time step of displacement; deep search algorithm and width search algorithm are combined to identify and find unit position information of capillary trapping. The application realizes accurate and efficient simulation of the invasion-percolation displacement process in the rough fissure by considering the influence of the curvature in the plane on the invasion displacement process and combining efficient data structure and search algorithm.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation of multiphase flow in fractured media, and specifically to a simulation method and system for rough fracture intrusion-percolation two-phase flow. Background Technology

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

[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 computation is often used as an important means to reveal the two-phase flow displacement process and interface morphology evolution. Among them, the intrusion-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 will be intruded and displaced first. However, the previously proposed intrusion-percolation two-phase flow models ignore the contribution of in-plane curvature to capillary pressure, which leads to a looser interface morphology compared with experimental results. Furthermore, the above models often generate huge computational loads when considering capillary trapping phenomena and updating displacement interface information, resulting in low computational efficiency. Therefore, it is necessary to improve the traditional intrusion-percolation model and propose an efficient numerical simulation method to achieve accurate and rapid simulation of the intrusion-percolation two-phase flow displacement process in rough fractures. Summary of the Invention

[0004] To overcome the shortcomings of the prior art, this invention considers the influence of in-plane curvature on capillary pressure and the capillary trapping phenomenon, and proposes a fast and accurate simulation method and system for rough fracture intrusion-percolation two-phase flow.

[0005] According to one aspect of the present invention, the present invention provides a simulation method for intrusion-percolation two-phase flow in rough fractures. The simulation method includes: dividing the rough fracture aperture field using planar mesh elements to establish a numerical model of intrusion-percolation two-phase flow; simulating the displacement process as a quasi-static interface propagation process; and displacing the element with the minimum capillary pressure threshold within a unit time step by intrusion. The method also employs the Young-Laplace equation to determine the capillary pressure threshold of each element to be displaced by calculating the in-plane radius of curvature and the radius of curvature in the aperture direction.

[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; a depth search algorithm and a width search algorithm are combined to identify and find the unit location 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 radius of curvature and the radius of curvature in the opening direction. This includes: using a boundary tracing algorithm to track the coordinates of several neighboring interface units of each unit to be displaced; fitting the tracked coordinates of the interface units to the curvature circle equation to determine the fitting coefficients based on the curvature circle equation; calculating the corresponding in-plane radius of curvature based on the fitting coefficients; calculating the radius of curvature in the opening direction based on the static contact angle and the opening of each unit to be displaced; and combining this with the Young-Laplace equation to calculate the capillary pressure threshold of each unit to be displaced.

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

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

[0010] Furthermore, in accordance with the arrangement principle of the binary tree data structure, an initial tree is established, including: sorting all the units to be replaced according to the size of the capillary pressure threshold in ascending order; inserting all the units to be replaced into the binary tree according to the size of the capillary pressure threshold; the arrangement principle of the binary tree data structure is that the capillary pressure threshold of each node in the binary tree is greater than the capillary pressure threshold of the parent node and less than the capillary pressure threshold of the child node.

[0011] Furthermore, the nodes corresponding to the invaded units are deleted, and the nodes corresponding to the units to be replaced are added. 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 the 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 last node to the position where the child node is missing due to the node update at the front; checking whether the updated binary tree meets the layout principle, and if not, swapping the positions; adding the nodes corresponding to the new units to be replaced, 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 planar mesh elements, including: dividing the rough fracture aperture field into planar meshes based on the measured rough fracture structural characteristics, so as to generalize the entire rough fracture into an assembly of several tiny local flat plate elements; the local flat plate element is a parallelepiped, the side length of the local flat plate element is determined by the measurement accuracy of the fracture aperture field, and the thickness of the local flat plate element is equal to the local average aperture of the fracture at its location.

[0013] Furthermore, by combining depth search and width search algorithms, the unit location information of capillary capture is determined and found, including: obtaining the location information of the interface unit; based on the location information of the interface unit, using a depth search algorithm to search for areas where capillary capture has not occurred; excluding the areas where capillary capture has not occurred, using a width search algorithm to search for the unit location information where capillary capture has occurred.

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

[0015] According to one aspect of the present invention, a non-transitory computer-readable storage medium is provided that stores computer instructions that cause the computer to perform the simulation method of rough fracture intrusion-percolation two-phase flow.

[0016] The above technical solution proposes a simulation method for rough fracture intrusion-percolation two-phase flow that considers in-plane curvature and capillary trapping. Specifically, it proposes an in-plane curvature calculation method in capillary pressure calculation, uses a binary tree data structure to record and quickly update displacement process information, and combines depth search and width search algorithms to accelerate the identification of capillary trapping, thereby achieving efficient and accurate simulation of rough fracture intrusion-percolation displacement process.

[0017] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0018] (1) This invention proposes a method for calculating the in-plane curvature of the fracture displacement interface, thereby revealing the contribution of in-plane curvature to capillary pressure, realizing the accurate solution of the capillary pressure threshold, and the 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 cell position, invasion state and capillary pressure threshold. Within a single time step of displacement, the characteristics of the binary tree data structure are used to quickly find the cell with the minimum capillary pressure threshold and realize the quick deletion of the corresponding node of the invaded cell and the quick addition of the corresponding node of the cell 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 depth search algorithm and width search algorithm to accelerate the identification of capillary capture. That is, firstly, a fast depth search algorithm is used to exclude the area where no capillary capture occurs, and then the width search algorithm is used in the remaining area to identify the occurrence of capillary capture and its specific location. Attached Figure Description

[0021] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the accompanying drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the accompanying drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

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

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

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

[0025] Figure 4 The third sub-flowchart of the simulation method for rough fracture intrusion-percolation two-phase flow provided in the embodiments of the present invention.

[0026] Figure 5 The calculation flowchart is shown for the simulation method of rough fracture intrusion-percolation two-phase flow provided in the embodiment of the present invention.

[0027] Figure 6 This is a schematic diagram of planar mesh division for a rough crack provided in an embodiment of the present invention.

[0028] Figure 7 This is a schematic diagram of the opening direction curvature and in-plane curvature provided in an embodiment of the present invention.

[0029] Figure 8(a) is a roughness crack aperture distribution diagram provided in the embodiment of the present invention.

[0030] Figure 8(b) is a diagram of the intrusion process provided in an embodiment of the present invention.

[0031] Figure 8(c) shows the final intrusion result provided by the embodiment of the present invention.

[0032] Figure 9 This is a schematic diagram of the in-plane curvature calculation method provided in an embodiment of the present invention.

[0033] Figure 10(a) is a schematic diagram showing the association between unit location and intrusion state provided in an embodiment of the present invention.

[0034] Figure 10(b) is a schematic diagram showing the relationship between capillary pressure threshold and binary tree data structure provided in the embodiment of the present invention.

[0035] Figure 11(a) is a schematic diagram of the first update of the binary tree data structure provided in the embodiment of the present invention.

[0036] Figure 11(b) is a second update diagram of the binary tree data structure provided in the embodiment of the present invention.

[0037] Figure 11(c) is a schematic diagram of the third update of the binary tree data structure provided in the embodiment of the present invention.

[0038] Figure 12 A schematic diagram of capillary capture determination combining depth search algorithm and width search algorithm provided in this embodiment of the invention.

[0039] Figure 13 A schematic diagram of a simulation system for rough fracture intrusion-percolation two-phase flow provided in this embodiment of the invention. Detailed Implementation

[0040] It should be noted that:

[0041] The terms “comprising” and “having”, and any variations thereof, in the specification, claims, and accompanying drawings of this invention are intended to cover a non-exclusive inclusion, such as a process, method, system, product, or apparatus that includes a series of steps or units, not necessarily limited to those explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0042] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention. In addition, the technical features of the various embodiments or individual embodiments provided by the present invention can be arbitrarily combined to form new technical solutions. Such combinations are not bound by the order of steps and / or structural composition patterns, but must be based on the ability of those skilled in the art to implement them. When the combination of technical solutions is contradictory or cannot be implemented, it should be considered that such a combination of technical solutions does not exist and is not within the scope of protection claimed by the present invention.

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

[0044] Step S101: The rough fracture aperture field is divided using planar mesh elements to establish an intrusion-percolation two-phase flow numerical model. The displacement process is simulated as a quasi-static interface propagation process, in which the element with the minimum capillary pressure threshold within a unit time step is intruded and displaced.

[0045] In step S101, based on the measured rough fracture structural characteristics, the known rough fracture aperture field is divided into planar meshes, and the entire rough fracture is generalized into several tiny local flat plate elements (i.e., parallelepipeds, as shown in the attached figure). Figure 6The structure consists of several local plate units, including invaded units, displaced units, and uninvaded units. Invaded units represent local plate units within a fracture that are completely occupied by the invading phase (i.e., the displacing phase, such as the injected fluid). Displaced units represent local plate units located at the leading edge of the displacement interface, where the invading phase may soon invade; these are adjacent units to the invaded units. Uninvaded units represent local plate units within a fracture that are still completely occupied by the displaced phase (such as the original fluid). Uninvaded units are not adjacent to invaded units. The displacement interface is essentially the contact surface between two fluids (such as the invading phase and the displaced phase), and the interface unit is defined as the invaded unit at the displacement interface. The invading phase and the displaced phase can be gas or liquid. For example, the displaced phase is crude oil, and the invading phase is CO2. Injecting CO2 to displace the crude oil in the reservoir can 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 is different at every point in the rough crack. The rough crack aperture field refers to the aperture distribution at every point in planar space. Therefore, after dividing the rough crack aperture field into a planar mesh, each local flat plate element corresponds to an aperture. In this embodiment, the thickness of each local flat plate element is equal to the local average aperture of the crack at its location. Furthermore, in this embodiment, the side length of each local flat plate element is determined by the measurement accuracy of the rough crack aperture field to avoid invalid mesh refinement while meeting accuracy requirements, thus balancing simulation efficiency and reliability.

[0047] Furthermore, after constructing the numerical model of the intrusion-percolation two-phase flow, it is necessary to initialize the basic parameters, convergence conditions, and interface propagation process of the numerical model according to the simulation conditions. The core idea of ​​the intrusion-percolation two-phase flow numerical model is to study the intrusion law of the intruding 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 elements, and the total volume of the intruding phase. The contact angle is usually in the range of 0° to 90°, and the number of adjacent interface elements is usually set between 3 and 10, both of which are constants. The total volume of the intruding phase refers to the total volume of the intruding phase that needs to be injected 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 propagation process, where interface propagation represents the complete invasion of a cell within each time step. Boundary conditions include inlet boundary conditions, outlet boundary conditions, and impermeable boundary conditions on both sides. This embodiment also defines the initial interface as the inlet boundary, which 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 constrain the remaining two boundaries of the simulation region, preventing both the invading and displaced phase fluids from passing through. The convergence condition refers to the termination condition of the numerical simulation of the invading-percolation two-phase flow. In this embodiment, the convergence condition is the breakthrough of the displacement interface (the invading phase reaches the outlet boundary) or the volume of the local flat plate cell occupied by the invading phase equals the initially set total volume of the invading phase.

[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 radius of curvature and the radius of curvature in the opening direction.

[0050] In step S103, the contribution of in-plane curvature to capillary pressure is considered to achieve an accurate solution for the capillary pressure threshold, thereby realizing an accurate simulation of the intrusion-percolation displacement process. Please refer to the appendix for further details. Figure 2 The following will further describe step S103 (including steps S1031-S1035).

[0051] Step S1031 involves using a boundary tracing algorithm to track the coordinates of several neighboring interface units for each unit to be displaced. In this embodiment, the number of neighboring interface units is a parameter set during the initialization of the intrusion-percolation two-phase flow numerical model. Please refer to the appendix for details. Figure 9 Select a cell (gray diamond point) whose capillary threshold needs to be calculated, and trace the coordinates (black triangle points) of n neighboring interface cells within the vicinity of this point. This can be achieved through boundary tracing algorithms in image processing (e.g., using the bwtraceboundary function in MATLAB), which is not limited here.

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

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

[0054] Please refer to the appendix for details. Figure 7 In this embodiment, the capillary pressure threshold P c The capillary pressure threshold, calculated using the Young-Laplace equation, depends on the interfacial tension and the curvature of the interface in the opening direction and in-plane direction. Where γ represents interfacial tension; k1 and k2 represent the opening direction curvature and in-plane curvature of the element to be displaced, respectively; r1 and r2 are the corresponding opening direction curvature radius and in-plane curvature radius, respectively; θ c θ represents the static contact angle; b represents the opening of the element to be replaced. It should be noted that the interfacial tension γ and the static contact angle θ... c The aperture b is related to the properties of the fluid, so it is a known constant value and is a parameter set during the initialization of the intrusion-percolation two-phase flow numerical model. After planar mesh generation, it is also a known constant value. Therefore, based on the static contact angle θ... c The radius of curvature r1 in the opening direction is calculated by combining the opening b of each unit to be displaced with the interfacial tension γ and the in-plane radius of curvature r2 calculated in step S1035 above, so that the capillary pressure threshold P can be obtained. c .

[0055] Step S105: Using a binary tree data structure, record the replacement process information of each unit, and update the changed binary tree data structure within a single time step of the replacement.

[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, thereby efficiently managing and quickly retrieving the units to be displaced. Please refer to the appendix for details. Figure 3 The following will further describe step S105 (including steps S1051-S1057).

[0057] Step S1051: Define a one-dimensional array of unit position, invasion state, and capillary pressure threshold. In this embodiment, each local plate unit, except for the boundary, is connected to four adjacent units. Using the sub2ind function in MATLAB, the two-dimensional position coordinates (x, y) of all local plate units are converted into a one-dimensional index i, thus assigning a unique index i to each local plate unit. Referring to Figure 10(a), for example, the index i=8 corresponds to the local unit (2,3) in the 2nd row and 3rd column. The invasion state of each local unit is marked and stored in a one-dimensional array s(i). Here, s=2 indicates that the local plate unit is in an invaded state, i.e., it corresponds to an invaded unit; s=1 indicates that the local plate unit is in a pending invasion state, i.e., it corresponds to a unit to be replaced; and s=0 indicates that the local plate unit is in a non-invaded state, i.e., it corresponds to a non-invaded unit. Using the Young-Laplace equations described above, the capillary pressure threshold of a local plate element with an invasion state of s=1 is calculated and stored in a one-dimensional array r(i), which represents the degree of difficulty of invasion (the smaller the invasion capillary pressure threshold, the easier the element is to be invaded).

[0058] Step S1053: Establish an initial tree according to the arrangement principles of the binary tree data structure. Referring to Figure 10(b), in this embodiment, for units with an intrusion state of s=1, sort them from smallest to largest according to the value of r(i). Based on the arrangement principles of the binary tree data structure, add the units under the condition s=1 to the binary tree to establish the initial tree. The specific arrangement principles of the binary tree data structure include: ① Each node in the binary tree has one parent node and two child nodes (an exception is that if the total number of nodes in the binary tree is even, the parent node of the last node in the binary tree will only have one child node). ② The intrusion capillary pressure threshold of 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 the parent node is more easily intruded, which also makes the root node (the node at the top of the tree) the most easily intruded and replaced.

[0059] Step S1055: Within a single time step of displacement, the unit to be displaced with the minimum capillary pressure threshold is searched for invasion. Specifically, within one unit time step, the unit to be displaced with the minimum current capillary pressure threshold is selected for invasion, and the information of the unit to be displaced is updated. In the next time step, a new unit to be displaced with the minimum capillary pressure threshold is selected for invasion, and the above invasion process is continuously executed until the displacement interface is breached (i.e., the unit at the outlet is invaded) or the set total volume of the invading phase is reached. At this point, the numerical model of invading-percolation two-phase flow has reached the convergence condition, and the simulation is stopped. Based on the simulation results, the spatiotemporal evolution law of the displacement interface in the rough fracture, the distribution characteristics of the capillary trapping phase, and macroscopic characteristic indicators such as saturation and interface analysis dimension are analyzed, thereby revealing the law of fracture invading-percolation two-phase flow under the influence of different rough structures (as shown in the rough fracture aperture distribution diagram in Figure 8(a), the invading process diagram in Figure 8(b), and the final invading result diagram in Figure 8(c)).

[0060] Step S1057: Delete the node corresponding to the invaded unit and add the node corresponding to the new unit to be replaced, updating the binary tree data structure. Specifically, please refer to Figures 11(a)-(c). The binary tree update process includes: ① Deleting the root node (operation 1 in Figure 11(a)); ② Checking the two child nodes corresponding to the root node and selecting the child node with the smaller capillary pressure threshold as the new root node (operation 2 in Figure 11(a)); ③ Repeating the same operation to continue updating the parent node and child nodes (operations 3 and 4 in Figure 11(a)) until the last node has no child nodes; ④ To maintain the integrity of the tree, move the last node to the position where the child nodes are missing due to the node update at the front (Figure 11(a)). Operation 5 in Figure 11(b)); ⑤ Check whether the updated binary tree satisfies the arrangement principle, that is, check whether the r(i) value of the moved child node is greater than that of the parent node. If not, swap the positions (operations 6 and 7 in Figure 11(b)); ⑥ Add the new node corresponding to the unit to be replaced whose intrusion state changes to s(i)=1 (operation 8 in Figure 11(c)); ⑦ Check again whether the updated binary tree satisfies the arrangement principle, that is, check whether the r(i) value of the added new node 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 embodiments, a binary tree data structure is used to record displacement process information by defining a one-dimensional array of cell location, 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 cell with the minimum capillary pressure threshold, and to quickly delete the corresponding node of the invaded cell and quickly add the corresponding node of the cell to be displaced, thereby rapidly updating the displacement process information and significantly improving the calculation speed of the invasion-percolation numerical model.

[0062] Step S107: Combine depth search algorithm and width search algorithm to identify and find the unit location information of capillary capture.

[0063] In step S107, a fast depth-first search algorithm is first used to exclude areas where capillary capture has not occurred. Then, a width-first search algorithm is used within the remaining areas to determine the occurrence and specific location of capillary capture, thereby accelerating the identification of capillary capture. Please refer to the appendix for further details. Figure 4 and Figure 12 The following will further describe step S107 (including steps S1071-S1075).

[0064] Step S1071: Obtain the position information of the interface unit. In this embodiment, based on the stored binary tree data, the position information of the foremost invading unit in the flow direction (the flow direction refers to the 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 towards the fracture boundary, i.e., the y direction) is extracted, that is, the position information of the invading unit at the displacement interface, and stored in a one-dimensional array.

[0065] Step S1073: Based on the location information of the interface units, a depth search algorithm is used to search for areas where capillary capture has not occurred. Specifically, based on the location of the intruded units at the displacement interface, the search proceeds forward along the flow direction and the sidewall direction to obtain the gray area in Figure 8, which is the area where capillary capture has not occurred.

[0066] Step S1075: Exclude areas where capillary capture did not occur, and use a width search algorithm to search for the location information of cells where capillary capture occurred. Specifically, continue searching the black areas outside the gray areas in Figure 8 and determine whether capillary capture has occurred. The method for determining capillary capture is to determine whether there are simply connected regions that are not connected to the main flow region after the displacement interface, which can be achieved using the bwlabeln and regionprops functions in MATLAB.

[0067] It should be noted that capillary capture refers to the phenomenon where, during the displacement of the displaced phase by the invading 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 clumps or clusters that are not connected to the main flow region (i.e., a single-connected region). Understandably, by precisely locating the capillary capture position, the simulation method for invading-percolation two-phase flow provided by this invention can more realistically simulate the complex flow in rough fractures.

[0068] Please refer to the appendix for details. Figure 13This invention also provides a simulation system 100 for rough fracture intrusion-percolation two-phase flow. The simulation system 100 includes a mesh generation module 1, a threshold calculation module 2, a data management module 3, and a capillary capture module 4. The mesh generation module 1 uses planar mesh elements to divide the rough fracture aperture field, establishing an intrusion-percolation two-phase flow numerical model. The displacement process is simulated as a quasi-static interface propagation process, where the element with the minimum capillary pressure threshold within a unit time step is intruded and displaced. The threshold calculation module 2 uses the Young-Laplace equation to determine the capillary pressure threshold of each element to be displaced by calculating the in-plane radius of curvature and the radius of curvature in the aperture direction. The data management module 3 uses a binary tree data structure to record the displacement process information of each element and updates the changing binary tree data structure within a single time step of displacement. The capillary capture module 4 combines depth search and width search algorithms to identify and locate the capillary capture element positions.

[0069] In summary, this invention discloses an efficient simulation method for intrusion-percolation two-phase flow in rough fractures, considering in-plane curvature and capillary trapping. Based on traditional numerical models of intrusion-percolation two-phase flow, this method considers the contribution of in-plane curvature to capillary pressure and proposes an in-plane curvature calculation method based on the curvature circle equation and the least squares method, achieving accurate solution for the capillary pressure threshold. During the intrusion displacement process, a binary tree data structure is used to record and quickly update displacement process information such as cell position, intrusion state, and capillary pressure threshold. Simultaneously, by combining depth search and width search algorithms, rapid identification of capillary trapping is achieved. In other words, by considering the influence of in-plane curvature on the intrusion displacement process and combining efficient data structures and search algorithms, this invention achieves accurate and efficient simulation of the intrusion-percolation displacement process in rough fractures, thereby revealing the laws governing intrusion-percolation two-phase flow in fractures 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, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the technical solutions of the embodiments of the present invention.

Claims

1. A simulation method for rough fracture intrusion-percolation two-phase flow, characterized in that, The simulation method includes: Planar mesh elements are used to divide the rough fracture aperture field, and an intrusion-percolation two-phase flow numerical model is established. The displacement process is simulated as a quasi-static interface propagation process, in which the element with the minimum capillary pressure threshold within a unit time step is intruded and displaced. Using the Young-Laplace equation, the capillary pressure threshold of each displaced element is determined by calculating the in-plane radius of curvature and the radius of curvature in the opening direction. This includes: using a boundary tracing algorithm to track the coordinates of several neighboring interface elements of each displaced element; fitting the tracked coordinates of the interface elements to the curvature circle equation to determine the fitting coefficients based on the curvature circle equation; and calculating the corresponding in-plane radius of curvature based on the fitting coefficients, using the following formula: r2 represents the radius of curvature in the plane. , and The fitting coefficients are represented; the radius of curvature in the opening direction is calculated based on the static contact angle and the opening of each unit to be displaced, using the following formula: r1 represents the radius of curvature in the opening direction, θ c Let represent the static contact angle, and b represent the opening of the element to be displaced; the capillary pressure threshold of each element to be displaced is calculated using the Young-Laplace equation, as follows: γ represents interfacial tension; k1 and k2 represent the opening direction curvature and in-plane curvature of the element to be displaced, respectively; P c Indicates the capillary pressure threshold; A binary tree data structure is used to record the displacement process information of each unit. Within a single time step of displacement, the changing binary tree data structure is updated, including: defining a one-dimensional array of unit position, invasion state, and capillary pressure threshold; building an initial tree according to the arrangement principle of the binary tree data structure; finding the unit to be displaced with the minimum capillary pressure threshold within a single time step of displacement and invading it; deleting the corresponding node of the invaded unit and adding the corresponding node of the new unit to be displaced, and updating the binary tree data structure. By combining depth search and width search algorithms, the location information of capillary capture units can be identified and located.

2. The simulation method for rough fracture intrusion-percolation two-phase flow as described in claim 1, characterized in that, Define a one-dimensional array of cell location, intrusion state, and capillary pressure threshold, including: Convert the two-dimensional position coordinates of each unit into a one-dimensional index; Each unit is marked with an intrusion status and stored in a one-dimensional array; Calculate the capillary pressure threshold for each unit to be displaced and store it in another dimensional array.

3. The simulation method for rough fracture intrusion-percolation two-phase flow as described in claim 1, characterized in that, Following the layout principles of a binary tree data structure, construct the initial tree, including: Sort all the units to be displaced in ascending order of capillary pressure threshold. All units to be replaced are inserted into the binary tree according to the capillary pressure threshold. The arrangement principle of the binary tree data structure is that the capillary pressure threshold of each node in the binary tree is greater than the capillary pressure threshold of the parent node and less than the capillary pressure threshold of the child node.

4. The simulation method for rough fracture intrusion-percolation two-phase flow as described in claim 1, characterized in that, Delete the node corresponding to the invaded unit, add the node corresponding to the unit to be replaced, and update the binary tree data structure, including: Remove 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 last node to the position where the child node is missing due to the node update at the front. Check if the updated binary tree satisfies the placement principle; if not, swap the positions. Add the node corresponding to the new unit to be replaced, and check again whether the updated binary tree meets the layout principle. If not, swap their positions.

5. The simulation method for rough fracture intrusion-percolation two-phase flow as described in claim 1, characterized in that, The rough fracture aperture field is divided using planar mesh elements, including: Based on the measured structural characteristics of the rough fracture, the rough fracture aperture field is divided into planar grids to generalize the entire rough fracture into an assembly of several tiny local flat plate units. The local flat plate unit is a parallelepiped, and the side length of the local flat plate unit is determined by the measurement accuracy of the fracture aperture field. The thickness of the local flat plate unit is equal to the local average aperture of the fracture at its location.

6. The simulation method for rough fracture intrusion-percolation two-phase flow as described in claim 1, characterized in that, Combining depth search and width search algorithms, the location information of capillary capture units is determined and located, 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 areas where capillary capture has not occurred; After excluding the areas where capillary capture did not occur, a width search algorithm is used to search for the location information of the units where capillary capture occurred.

7. A simulation system for rough fracture intrusion-percolation two-phase flow, characterized in that, include: The mesh generation module uses planar mesh elements to divide the rough fracture aperture field, establishes an intrusion-percolation two-phase flow numerical model, and simulates the displacement process as a quasi-static interface propagation process. The element with the minimum capillary pressure threshold within a unit time step is intruded and displaced. The threshold calculation module uses the Young-Laplace equation to determine the capillary pressure threshold of each displaced unit by calculating the in-plane radius of curvature and the radius of curvature in the opening direction. This includes: using a boundary tracing algorithm to track the coordinates of several neighboring interface units of each displaced unit; fitting the tracked coordinates of these interface units to the curvature circle equation to determine the fitting coefficients; and calculating the corresponding in-plane radius of curvature based on the fitting coefficients, using the following formula: r2 represents the radius of curvature in the plane. , and The fitting coefficients are represented; the radius of curvature in the opening direction is calculated based on the static contact angle and the opening of each unit to be displaced, using the following formula: r1 represents the radius of curvature in the opening direction, θ c Let represent the static contact angle, and b represent the opening of the element to be displaced; the capillary pressure threshold of each element to be displaced is calculated using the Young-Laplace equation, as follows: γ represents interfacial tension; k1 and k2 represent the opening direction curvature and in-plane curvature of the element to be displaced, respectively; P c Indicates the capillary pressure threshold; The data management module uses a binary tree data structure to record the displacement process information of each unit. Within a single time step of displacement, it updates the changing binary tree data structure, including: defining a one-dimensional array of unit location, invasion status, and capillary pressure threshold; establishing an initial tree according to the arrangement principle of the binary tree data structure; finding the unit to be displaced with the minimum capillary pressure threshold within a single time step of displacement and invading it; deleting the corresponding node of the invaded unit and adding the corresponding node of the new unit to be displaced, and updating the binary tree data structure. The capillary capture module combines depth search and width search algorithms to determine and locate the unit position information of capillary capture.

8. A non-transitory computer-readable storage medium, characterized in that, The non-transitory computer-readable storage medium stores computer instructions that cause the computer to perform the simulation method for rough fracture intrusion-percolation two-phase flow as described in any one of claims 1 to 6.

Images