Hydraulic and electrical power optimal path parallel computing method and system for rock pore network

Abstract

The present application belongs to the technical field of digital rock and porous medium percolation, and provides a rock pore network hydraulic and electric optimal path parallel computing method and system, a micro-scale three-dimensional pore network model is constructed according to a pre-processed two-dimensional gray image slice; a boundary condition is added to the micro-scale three-dimensional pore network model, including a hydraulic field inlet pressure, an outlet pressure, an electric field inlet potential and an outlet grounding; a single pipe hydraulic and electric conduction coefficient is solved, and the water pressure and the electric potential of each node and the hydraulic flux and the electric flux of each pipe are obtained; the hydraulic flux and the electric flux of each pipe are taken as the weight of the corresponding pipe, the hydraulic and electric flux path with the maximum path weight is selected as the hydraulic and electric optimal path of the rock sample, and the present application realizes the synchronous numerical simulation of hydraulic and electric transport in the three-dimensional pore network for the first time.

Description

Technical Field

[0001] This invention belongs to the field of digital rock and porous media seepage technology, specifically relating to a parallel calculation method and system for optimal hydraulic and electrical paths in rock pore networks. Background Technology

[0002] The statements in this section are merely background information related to the present invention and do not necessarily constitute prior art.

[0003] In the study of porous rock media, accurately obtaining the optimal hydraulic and electrical paths and tortuosity is crucial for understanding seepage and electrical transport characteristics. Traditional methods are often limited to two-dimensional analysis or cannot simultaneously consider the two transport processes, resulting in insufficient understanding of the hydraulic and electrical transport laws at the three-dimensional pore level. Summary of the Invention

[0004] To address the aforementioned problems, this invention proposes a parallel calculation method and system for optimal hydraulic and electrical paths in rock pore networks. This invention solves the problem of difficulty in obtaining optimal hydraulic and electrical paths and tortuosity of existing rock samples at the three-dimensional pore level, and for the first time realizes synchronous numerical simulation of hydraulic and electrical transport in a three-dimensional pore network.

[0005] According to some embodiments, the present invention adopts the following technical solution:

[0006] A parallel computational method for optimal hydraulic and electrical paths in rock pore networks includes the following steps:

[0007] Two-dimensional grayscale image slices of rock samples are obtained and preprocessed.

[0008] A microscale three-dimensional pore network model is constructed based on the preprocessed two-dimensional grayscale image slices.

[0009] Boundary conditions are added to the microscale three-dimensional porous network model, including inlet and outlet pressures of the hydraulic field, and inlet potential and outlet grounding of the electric field.

[0010] Based on the boundary conditions, the hydraulic and electrical conduction coefficients of a single pipe are solved, and the water pressure and electric potential at each node, as well as the hydraulic and electrical flux of each pipe, are obtained.

[0011] The hydraulic and electrical flux of each pipe are used as the weight of the corresponding pipe. The path with the largest hydraulic and electrical flux weight is selected as the optimal hydraulic and electrical path for the rock sample. The path weight of the optimal path is the sum of the weights of all the pipes it passes through.

[0012] As an alternative implementation, the process of obtaining two-dimensional grayscale image slices of rock samples and preprocessing them includes: performing image denoising based on Gaussian filtering in filtering and denoising, performing image segmentation, and calculating the maximum inter-class variance of the optimal segmentation threshold.

[0013] As a further defined implementation, the process of performing image segmentation and calculating the optimal segmentation threshold to maximize the inter-class variance includes performing image segmentation using the Otsu algorithm, calculating the optimal threshold T to maximize the inter-class variance, and the calculation formula is as follows:

[0014] ;

[0015] Where T represents the segmentation threshold; This represents the inter-class variance when the threshold is T. This represents the proportion of the foreground class in the image when the threshold is T. This represents the proportion of the background class in the image when the threshold is T. This represents the average grayscale value of the foreground class when the threshold is T. This represents the average grayscale value of the background class when the threshold is T.

[0016] As an alternative implementation, the process of constructing a microscale three-dimensional pore network model based on preprocessed two-dimensional grayscale image slices includes: reading the preprocessed two-dimensional grayscale image slices to form a three-dimensional data volume;

[0017] The two-dimensional grayscale image slices after preprocessing and image segmentation are stacked sequentially to form a three-dimensional pore network model.

[0018] The grayscale image is subjected to clear binarization to obtain the final binary 3D data;

[0019] Based on the obtained binary three-dimensional data, each pore phase voxel is assigned the distance to the nearest solid phase boundary, and the skeleton representing the center line of the pore space is extracted and calculated using the three-dimensional Barthalli distance transformation.

[0020] Based on the obtained solid boundary distance and pore space centerline skeleton, the maximum sphere algorithm is used for identification and positioning to identify all non-overlapping spheres with the largest possible radius in the three-dimensional pore space.

[0021] Based on the identified spheres, spheres in contact with each other are defined as belonging to the same pore; and spheres from two different pore clusters are defined as a throat, and the pore radius and length in the three-dimensional pore network model are calculated.

[0022] As a further defined implementation, the process of stacking preprocessed and segmented two-dimensional grayscale image slices sequentially to form a three-dimensional porous network model includes: layer-to-layer registration based on mutual information maximization, feature point matching, and phase correlation methods, with the specific calculation formula as follows:

[0023] ;

[0024] in, It is the mutual information between slices A and B, that is, it measures the degree of correlation between the two. It is the joint probability distribution of gray value a in slice A and gray value b in slice B; and It is a marginal probability distribution.

[0025] As a further defined implementation, based on the obtained binary three-dimensional data, each pore phase voxel is assigned its distance to the nearest solid phase boundary, and a skeleton representing the centerline of the pore space is extracted. The calculation process using three-dimensional Bartsarian distance transformation includes the following formula:

[0026] ;

[0027] in, This is a Euclidean distance transformation of coordinate points in three-dimensional space. Represents a three-dimensional coordinate point on the surface of an object; This represents the square of the Euclidean distance, measuring the current point. With a point on the boundary The square of the straight-line distance between them.

[0028] As a further defined implementation, the process of identifying and locating the maximum sphere algorithm based on the obtained solid boundary distance and pore space centerline skeleton, and identifying all non-overlapping spheres with the largest possible radius in the three-dimensional pore space, includes: first creating an empty maximum heap, sorting them by radius from largest to smallest, and taking all pore phase voxels as potential sphere center candidates, but with an initial radius of 0.

[0029] Remove the candidate sphere with the largest radius from the priority queue. , inspection Center, radius Whether the sphere overlaps with any existing sphere that has been placed, and whether it is completely contained within the aperture phase;

[0030] For the existing center of the ball and radius New ball center and radius satisfy Then it is determined to be overlapping;

[0031] If there is no overlap, add this sphere to the network and update the radius of the candidate points in its neighborhood. Value and re-enqueue;

[0032] The algorithm stops when the priority queue is empty or when the radius of all remaining candidate balls is less than a very small threshold.

[0033] As a further defined implementation, the process of calculating the pore radius and length in the three-dimensional pore network model includes: based on the identified spheres, defining spheres that are in contact with each other as the same pore; and defining spheres from two different pore clusters as a throat, the geometric properties of which are determined by the position and radius of the sphere pair they connect.

[0034] Independent geometric parameters are calculated for each pore and throat. The formula for calculating pore parameters is as follows:

[0035] ;

[0036] in, For the equivalent radius, It is the union of the volumes of all spheres in the pore;

[0037] The throat length is determined by the shortest path length between the centers of the two pore clusters, and then the average radius, standard deviation of radius, and average length of all pores are calculated.

[0038] As an alternative implementation method, the process of solving the hydraulic and electrical conduction coefficients of a single pipe based on boundary conditions includes: hydraulic conduction coefficient The calculation formula is as follows:

[0039] ;

[0040] in, Where the pipe radius is; For fluid dynamic viscosity; This refers to the length of the pipe.

[0041] Electricity conduction coefficient The calculation formula is as follows:

[0042] ;

[0043] in, Where the pipe radius is; The fluid conductivity; This refers to the length of the pipe.

[0044] Based on the obtained hydraulic and electrical conduction coefficients, two sparse matrices are constructed in parallel to achieve the goal of simultaneously solving for hydraulic and electrical pressure.

[0045] A hydraulic sparse matrix is ​​constructed based on the hydraulic conduction coefficient and the connectivity between adjacent nodes. For each internal node in the sparse matrix, the flow balance equation is:

[0046] ;

[0047] in, Let ik be the hydraulic conductivity at point ik; and The pressure at points i and k;

[0048] Based on the boundary conditions at the inlet and outlet nodes, the LU decomposition method is used to efficiently solve the sparse matrix and obtain the hydraulic pressure of each node.

[0049] The electric potential at each node is also calculated by constructing the conductance matrix based on the electric conduction coefficient and the connectivity between adjacent nodes.

[0050] Based on the boundary conditions at the inlet and outlet nodes, the LU decomposition method is used to solve the linear equation system and obtain the potential of each node.

[0051] As an alternative implementation method, the process of selecting the hydraulic and electrical flux paths with the largest path weights as the optimal hydraulic and electrical flux paths for rock samples, using the hydraulic and electrical flux of each pipe as the corresponding pipe weights, includes: using the Dijkstra algorithm to perform optimal path analysis on the three-dimensional pore network model, and simultaneously obtaining the optimal hydraulic and electrical paths through parallel computing.

[0052] Let the set of nodes in the three-dimensional porous network model be . The pipeline assembly is Source The starting point and ending point of the path For the endpoint of the path, the weight of the pipeline between adjacent nodes u and v in each path. Defined as the absolute value of the water pressure difference and electric potential at both ends of the pipe connecting node u and node v, i.e.:

[0053] ;

[0054] in, , These are the water pressure and potential values ​​at node u at one end of the pipeline and node v at the other end of the pipeline, respectively. The greater the pressure difference between the two ends of the pipeline, the greater the flow rate in the pipeline, and the higher the weight.

[0055] In Dijkstra's algorithm, a distance array d[v] is introduced, which represents the current known maximum cumulative path weight from the source point s to the node v. The initialization condition is set to d[s]=0.

[0056] During the algorithm iteration, for the current node u, traverse all its adjacent nodes v, if the following conditions are met... If the sum of the path weights from the source point through u to v is greater than the currently recorded d[v], then update d[v] and record the predecessor node for subsequent path backtracking.

[0057] A parallel computational system for optimal hydraulic and electrical paths in rock pore networks includes:

[0058] The image acquisition module is configured to acquire two-dimensional grayscale image slices of rock samples and preprocess them;

[0059] The 3D pore network model construction module is configured to construct a microscale 3D pore network model based on preprocessed 2D grayscale image slices.

[0060] The boundary condition addition module is configured to add boundary conditions to the microscale three-dimensional porous network model, including inlet and outlet pressure of the hydraulic field, and inlet potential and outlet grounding of the electric field.

[0061] The hydraulic and electric conduction coefficient calculation module is configured to solve the hydraulic and electric conduction coefficients of a single pipe based on boundary conditions, and obtain the water pressure and electric potential at each node, as well as the hydraulic flux and electric flux of each pipe.

[0062] The optimal path determination module is configured to use the hydraulic flux and electrical flux of each pipe as the weight of the corresponding pipe, and select the hydraulic and electrical flux path with the largest path weight as the optimal hydraulic and electrical path of the rock sample. The path weight of the optimal path is the sum of the weights of all pipes it passes through.

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

[0064] This invention solves the problem of difficulty in obtaining the optimal hydraulic and electrical paths and tortuosity of existing rock samples from the perspective of three-dimensional pores, and for the first time realizes the synchronous numerical simulation of hydraulic and electrical transport in a three-dimensional pore network.

[0065] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, preferred embodiments are described below in detail with reference to the accompanying drawings. Attached Figure Description

[0066] The accompanying drawings, which form part of this invention, are used to provide a further understanding of the invention. The illustrative embodiments of the invention and their descriptions are used to explain the invention and do not constitute an improper limitation of the invention.

[0067] Figure 1This application provides a parallel calculation method and system for hydraulic and electrodynamic optimal paths and tortuosity of rock pore networks based on a three-dimensional pore network model, as provided in Embodiment 1 of this application.

[0068] Figure 2 The histogram of frequency distribution of pore radius in shale sample;

[0069] Figure 3 This is a diagram showing the water pressure distribution at the nodes of a three-dimensional pore network model.

[0070] Figure 4 This is a diagram showing the potential distribution of nodes in a three-dimensional porous network model.

[0071] Figure 5 A diagram showing the hydraulic flux distribution in a pipe within a three-dimensional pore network model;

[0072] Figure 6 This is a diagram showing the power flux distribution in a three-dimensional pore network model of a pipeline.

[0073] Figure 7 The hydraulically optimal path for the three-dimensional pore network model;

[0074] Figure 8 The optimal power path for a three-dimensional porous network model;

[0075] Figure 9 The curves are fitted to the tortuosity of the optimal hydraulic path and the optimal power path for shale. Detailed Implementation

[0076] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0077] It should be noted that the following detailed description is illustrative and intended to provide further explanation of the invention. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0078] It should be noted that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit the scope of exemplary embodiments according to the invention. As used herein, the singular form is intended to include the plural form as well, unless the context clearly indicates otherwise. Furthermore, it should be understood that when the terms "comprising" and / or "including" are used in this specification, they indicate the presence of features, steps, operations, devices, components, and / or combinations thereof.

[0079] Where there is no conflict, the embodiments and features described in this application may be combined with each other.

[0080] Example 1

[0081] This embodiment provides a parallel calculation method for the optimal hydraulic and electrical paths and tortuosity of rock pore networks based on a three-dimensional pore network model. It solves the problem of the difficulty in obtaining the optimal hydraulic and electrical paths and tortuosity of existing rock samples from the three-dimensional pore level, and realizes the synchronous numerical simulation of hydraulic and electrical transport in a three-dimensional pore network for the first time.

[0082] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0083] Example 1

[0084] like Figure 1 As shown, this invention provides a parallel calculation method and system for the optimal hydraulic and electrodynamic paths and tortuosity of rock pore networks based on a three-dimensional pore network model. The method includes the following steps:

[0085] Step S1: Obtain original rock sample slices based on CT scans. The specific process is as follows:

[0086] In this embodiment, shale core samples were used, and a series of high-resolution two-dimensional grayscale image slices of the shale core samples were obtained by means of X-ray computed tomography (micron CT).

[0087] The preprocessing of the two-dimensional grayscale image slices is based on Gaussian filtering in the filtering and noise reduction process to reduce noise. While smoothing the noise, it can better preserve the edge information between the pores and the skeleton.

[0088] The general calculation formula for filtering and noise reduction is as follows:

[0089] ;

[0090] in For the input image in grayscale value; To output the image in coordinates The grayscale value after noise reduction; The weights of the convolution kernel (which need to be normalized): (to avoid brightness shift) and Then the convolution kernel is in and Radius in the direction;

[0091] The formula for calculating a two-dimensional Gaussian filter is as follows:

[0092] ;

[0093] in, It is the mathematical expression of a two-dimensional Gaussian function, which defines the weight distribution of neighboring pixels during the filtering process; The relative coordinates of the current pixel with respect to the center of its neighborhood; The standard deviation of the Gaussian distribution controls the decay rate of the weights; This is a normalization constant, ensuring that the sum of the integrals of the Gaussian function over the entire plane is 1;

[0094] For the preprocessed 2D grayscale image, the Otsu algorithm can be used for image segmentation to automatically calculate the optimal threshold T that maximizes the inter-class variance. The calculation formula is as follows:

[0095] ;

[0096] Where T represents the segmentation threshold; This represents the inter-class variance when the threshold is T. This represents the proportion of the foreground class in the image when the threshold is T. This represents the proportion of the background class in the image when the threshold is T. This represents the average grayscale value of the foreground class when the threshold is T. This represents the average grayscale value of the background class when the threshold is T.

[0097] Based on the Otsu algorithm described above, the difference between the image before and after segmentation is most obvious.

[0098] Step S2: Reshape the 2D slices and construct a 3D pore network model

[0099] All processed two-dimensional grayscale image slices from step S1 are read into memory to form a three-dimensional data volume;

[0100] The two-dimensional grayscale image slices after preprocessing and image segmentation are stacked sequentially to form a three-dimensional pore network model.

[0101] The registration between layers is based on mutual information maximization, feature point matching, and phase correlation methods. The specific calculation formula is as follows:

[0102] ;

[0103] in, It is the mutual information between slices A and B, that is, it measures the degree of correlation between the two. It is the joint probability distribution of gray value a in slice A and gray value b in slice B; and It is a marginal probability distribution.

[0104] To facilitate the acquisition of various parameters of the subsequent three-dimensional porous network model, Avizo software was used to perform clear binarization processing on the grayscale image to obtain the final binary three-dimensional data.

[0105] Based on the obtained binary 3D data, each pore phase voxel needs to be assigned its distance to the nearest solid boundary, and the skeleton representing the centerline of the pore space needs to be extracted. The calculation is performed using the 3D Barthalya distance transformation, and the calculation formula is as follows:

[0106] ;

[0107] in, This is a Euclidean distance transformation of coordinate points in three-dimensional space. Represents a three-dimensional coordinate point on the surface of an object; This represents the square of the Euclidean distance, which measures the distance from the current point. With a point on the boundary The square of the straight-line distance between them;

[0108] Based on the obtained solid boundary distance and pore space centerline skeleton, the maximum sphere algorithm is used for identification and positioning to identify all non-overlapping spheres with the largest possible radius in the three-dimensional pore space.

[0109] First, create an empty max-heap, sort it by radius from largest to smallest, and treat all pore phase voxels as potential sphere center candidates, but with an initial radius of 0;

[0110] Remove the candidate sphere with the largest radius from the priority queue. , inspection Center, radius Whether the sphere overlaps with any existing sphere that has been placed, and whether it is completely contained within the aperture phase;

[0111] For the existing center of the ball and radius New ball center and radius satisfy Then it is determined to be overlapping;

[0112] If there is no overlap, add this sphere to the network and update the radius of the candidate points in its neighborhood. Value and re-enqueue;

[0113] The algorithm stops when the priority queue is empty or the radius of all remaining candidate balls is less than a very small threshold.

[0114] Based on the identified spheres, we define spheres that are in contact with each other as belonging to the same pore; and we define spheres from two different pore clusters as a throat, the geometric properties of which are determined by the position and radius of the sphere pair they connect.

[0115] Independent geometric parameters are calculated for each pore and throat. The formula for calculating pore parameters is as follows:

[0116] ;

[0117] in, For the equivalent radius, It is the union of the volumes of all spheres in the pore;

[0118] The throat length is determined by the shortest path length between the centers of the two pore clusters.

[0119] Based on the pore radius and length in the obtained three-dimensional pore network model, the average radius, standard deviation of radius, and average length of all pores are obtained.

[0120] like Figure 2 As shown, the average radius of the shale core samples taken in the example is 25 micrometers, the standard deviation of the radius is 0.50372, and the average length is 500 micrometers.

[0121] Step S3: Set boundary conditions

[0122] Boundary conditions were set for the microscale three-dimensional porous network model in the hydraulic and electric fields. The nodal pressure of the hydraulic field inlet pressure (k=1 layer) was set to 1, and the nodal pressure of the hydraulic field outlet pressure (k=Nk layer) was set to 0. The nodal potential of the electric field inlet potential (k=1 layer) was set to 1, and the nodal potential of the electric field outlet potential (k=Nk layer) was set to 0. There was no flux in other directions.

[0123] Step S4: Calculate the hydraulic conductivity coefficient of each pipe based on the pipe radius. and electrical conductivity It describes the laws governing water flow and electric current in laminar and steady-state motion, and calculates the pressure at each node of the model and the flow rate in the pipeline.

[0124] Calculate the hydraulic and electrical conduction coefficients of each shale pipeline based on the obtained radius distribution of the shale pipelines;

[0125] hydraulic conductivity The calculation formula is as follows:

[0126] ;

[0127] in, Where the pipe radius is; For fluid dynamic viscosity; This refers to the length of the pipe.

[0128] Electricity conduction coefficient The calculation formula is as follows:

[0129] ;

[0130] in, Where the pipe radius is; The fluid conductivity; This refers to the length of the pipe.

[0131] Based on the obtained hydraulic and electrical conduction coefficients, two sparse matrices are constructed in parallel to achieve the goal of simultaneously solving for hydraulic and electrical pressure.

[0132] A hydraulic sparse matrix is ​​constructed based on the hydraulic conduction coefficient and the connectivity between adjacent nodes. For each internal node in the sparse matrix, the flow balance equation is:

[0133] ;

[0134] in, Let ik be the hydraulic conductivity at point ik; and The pressure at points i and k;

[0135] Based on the boundary conditions at the inlet and outlet nodes, the LU decomposition method is used to efficiently solve the sparse matrix and obtain the hydraulic pressure of each node.

[0136] The electric potential at each node is also calculated by constructing the conductance matrix based on the electric conduction coefficient and the connectivity between adjacent nodes.

[0137] Similarly, based on the boundary conditions at the inlet and outlet nodes, the LU decomposition method is used to solve the linear equation system to obtain the potential of each node.

[0138] Figure 3 and Figure 4 Here are the water pressure distribution diagram and the electric potential distribution diagram of the nodes in the three-dimensional pore network model of Example 1;

[0139] The hydraulic pressure and electric potential of the nodes are calculated synchronously and in parallel.

[0140] The hydraulic flux is;

[0141] ;

[0142] Electricity flux is;

[0143] ;

[0144] Figure 5 and Figure 6 The hydraulic flux distribution and electrical flux distribution of the pipe in the three-dimensional pore network model of Example 1 are shown.

[0145] Step S5: Use an improved path algorithm to search in parallel for the path with the largest average weight, which will be the optimal path for both hydraulic and electric power.

[0146] The path algorithm uses Dijkstra's algorithm to perform optimal path analysis on the three-dimensional pore network model, and obtains the hydraulic optimal path and the electric optimal path simultaneously through parallel computing.

[0147] Let the set of nodes in the three-dimensional porous network model be . The pipeline assembly is Source The starting point and ending point of the path The endpoint of the path. The weights of the pipelines between adjacent nodes u and v in each path. Defined as the absolute value of the water pressure difference and electric potential at both ends of the pipe connecting node u and node v, i.e.:

[0148] ;

[0149] in, , These are the water pressure and potential values ​​at node u at one end of the pipeline and node v at the other end of the pipeline, respectively. The greater the pressure difference between the two ends of the pipeline, the greater the flow rate in the pipeline, and therefore the higher the weight.

[0150] In Dijkstra's algorithm, a distance array d[v] is introduced, representing the current known maximum cumulative path weight from the source node s to node v. The initialization condition is set to d[s] = 0;

[0151] During the algorithm iteration, for the current node u, traverse all its adjacent nodes v, if the following conditions are met... If the sum of the path weights from the source point through u to v is greater than the currently recorded d[v], then d[v] is updated, and the predecessor node is recorded for subsequent path backtracking.

[0152] like Figure 7 , Figure 8 As shown, the optimal hydraulic path and the optimal electric path for shale gas are obtained in this embodiment;

[0153] Step S6: Calculate hydraulic tortuosity based on the optimal path and power tortuosity It also outputs comparative analysis results.

[0154] The tortuosity of the optimal hydraulic path in sandstone is calculated using the following formula:

[0155] .

[0156] in, This represents the actual length of the optimal path for shale hydraulic and electric power generation. This represents the straight-line distance between the starting and ending points of the optimal path for shale hydraulic and electric power operations.

[0157] In the shale used in the example, the tortuosity of the hydraulic optimal path can reach 2.1379, while the tortuosity of the electric optimal path is 1.7586. Obviously, the tortuosity of the hydraulic optimal path is greater than that of the electric optimal path.

[0158] After 200 iterations of obtaining the optimal paths for shale hydraulics and electricity, an empirical calculation formula between the two can be fitted:

[0159] ;

[0160] Where x is the tortuosity of the optimal hydraulic path for shale, and y is the tortuosity of the optimal power path for shale.

[0161] like Figure 9 As shown, the tortuosity of the optimal path for shale hydraulics is 0.2 times higher than that of the optimal path for shale power.

[0162] Example 2

[0163] A parallel computational system for optimal hydraulic and electrical paths in rock pore networks includes:

[0164] The image acquisition module is configured to acquire two-dimensional grayscale image slices of rock samples and preprocess them;

[0165] The 3D pore network model construction module is configured to construct a microscale 3D pore network model based on preprocessed 2D grayscale image slices.

[0166] The boundary condition addition module is configured to add boundary conditions to the microscale three-dimensional porous network model, including inlet and outlet pressure of the hydraulic field, and inlet potential and outlet grounding of the electric field.

[0167] The hydraulic and electric conduction coefficient calculation module is configured to solve the hydraulic and electric conduction coefficients of a single pipe based on boundary conditions, and obtain the water pressure and electric potential at each node, as well as the hydraulic flux and electric flux of each pipe.

[0168] The optimal path determination module is configured to use the hydraulic flux and electrical flux of each pipe as the weight of the corresponding pipe, and select the hydraulic and electrical flux path with the largest path weight as the optimal hydraulic and electrical path of the rock sample. The path weight of the optimal path is the sum of the weights of all pipes it passes through.

[0169] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can be implemented in one or more computer-usable storage media (including, but not limited to, disk storage, etc.) containing computer-usable program code. CD - ROM It takes the form of a computer program product implemented on (such as optical memory, etc.).

[0170] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, as well as combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0171] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0172] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0173] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made by those skilled in the art without creative effort within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A parallel computation method for optimal hydraulic and electrical paths in rock pore networks, characterized in that, Includes the following steps: Two-dimensional grayscale image slices of rock samples are obtained and preprocessed. A microscale three-dimensional pore network model is constructed based on the preprocessed two-dimensional grayscale image slices. Boundary conditions are added to the microscale three-dimensional porous network model, including inlet and outlet pressures of the hydraulic field, and inlet potential and outlet grounding of the electric field. Based on the boundary conditions, the hydraulic and electrical conduction coefficients of a single pipe are solved, and the water pressure and electric potential at each node, as well as the hydraulic and electrical flux of each pipe, are obtained. The hydraulic and electrical flux of each pipe are used as the weight of the corresponding pipe. The path with the largest hydraulic and electrical flux weight is selected as the optimal hydraulic and electrical path for the rock sample. The path weight of the optimal path is the sum of the weights of all the pipes it passes through.

2. The parallel calculation method for optimal hydraulic and electrical paths in rock pore networks as described in claim 1, characterized in that, The process of obtaining two-dimensional grayscale image slices of rock samples and preprocessing them includes: denoising the image based on Gaussian filtering in filtering and denoising, performing image segmentation, and calculating the maximum inter-class variance of the optimal segmentation threshold. Alternatively, image segmentation processing, and the process of calculating the optimal segmentation threshold to maximize the inter-class variance, includes using the Otsu algorithm for image segmentation processing, calculating the optimal threshold T to maximize the inter-class variance, and the calculation formula is: Where T represents the segmentation threshold; This represents the inter-class variance when the threshold is T. This represents the proportion of the foreground class in the image when the threshold is T. This represents the proportion of the background class in the image when the threshold is T. This represents the average grayscale value of the foreground class when the threshold is T. This represents the average grayscale value of the background class when the threshold is T.

3. The parallel calculation method for optimal hydraulic and electrical paths in rock pore networks as described in claim 1, characterized in that, The process of constructing a microscale three-dimensional pore network model based on preprocessed two-dimensional grayscale image slices includes: reading the preprocessed two-dimensional grayscale image slices to form a three-dimensional data volume; The two-dimensional grayscale image slices after preprocessing and image segmentation are stacked sequentially to form a three-dimensional pore network model. The grayscale image is subjected to clear binarization to obtain the final binary 3D data; Based on the obtained binary three-dimensional data, each pore phase voxel is assigned the distance to the nearest solid phase boundary, and the skeleton representing the center line of the pore space is extracted and calculated using the three-dimensional Barthalli distance transformation. Based on the obtained solid boundary distance and pore space centerline skeleton, the maximum sphere algorithm is used for identification and positioning to identify all non-overlapping spheres with the largest possible radius in the three-dimensional pore space. Based on the identified spheres, spheres in contact with each other are defined as belonging to the same pore; and spheres from two different pore clusters are defined as a throat, and the pore radius and length in the three-dimensional pore network model are calculated.

4. The parallel calculation method for optimal hydraulic and electrical paths in rock pore networks as described in claim 3, characterized in that, The process of stacking preprocessed and segmented 2D grayscale image slices sequentially to form a 3D porous network model includes: layer-to-layer registration based on mutual information maximization, feature point matching, and phase correlation methods. The specific calculation formula is as follows: in, It is the mutual information between slices A and B, that is, it measures the degree of correlation between the two. It is the joint probability distribution of gray value a in slice A and gray value b in slice B; and It is a marginal probability distribution.

5. The parallel calculation method for optimal hydraulic and electrical paths in rock pore networks as described in claim 3, characterized in that, Based on the obtained binary 3D data, each pore phase voxel is assigned its distance to the nearest solid boundary, and a skeleton representing the centerline of the pore space is extracted. The calculation process using the 3D Bartschariah distance transformation includes the following formulas: in, This is a Euclidean distance transformation of coordinate points in three-dimensional space. Represents a three-dimensional coordinate point on the surface of an object; This represents the square of the Euclidean distance, measuring the current point. With a point on the boundary The square of the straight-line distance between them.

6. The parallel calculation method for optimal hydraulic and electrical paths in rock pore networks as described in claim 3, characterized in that, The process of identifying and locating all non-overlapping spheres with the largest possible radius in the three-dimensional pore space based on the obtained solid boundary distance and pore space centerline skeleton includes: first, creating an empty max-heap, sorting them by radius from largest to smallest, and using all pore phase voxels as potential sphere center candidates, but with an initial radius of 0; Remove the candidate sphere with the largest radius from the priority queue. , inspection Center, radius Whether the sphere overlaps with any existing sphere that has been placed, and whether it is completely contained within the aperture phase; For the existing center of the ball and radius New ball center and radius satisfy Then it is determined to be overlapping; If there is no overlap, add this sphere to the network and update the radius of the candidate points in its neighborhood. Value and re-enqueue; The algorithm stops when the priority queue is empty or when the radius of all remaining candidate balls is less than a very small threshold.

7. The parallel calculation method for optimal hydraulic and electrical paths in rock pore networks as described in claim 3, characterized in that, The process of calculating the pore radius and length in a three-dimensional pore network model includes: based on the identified spheres, defining spheres that are in contact with each other as the same pore; and defining spheres from two different pore clusters as a throat, the geometric properties of which are determined by the position and radius of the spheres they connect. Independent geometric parameters are calculated for each pore and throat. The formula for calculating pore parameters is as follows: in, For the equivalent radius, It is the union of the volumes of all spheres in the pore; The throat length is determined by the shortest path length between the centers of the two pore clusters, and then the average radius, standard deviation of radius, and average length of all pores are calculated.

8. The parallel calculation method for optimal hydraulic and electrical paths in rock pore networks as described in claim 1, characterized in that, The process of solving for the hydraulic and electrical conduction coefficients of a single pipe based on boundary conditions includes: hydraulic conduction coefficient The calculation formula is as follows: in, Where the pipe radius is; For fluid dynamic viscosity; This refers to the length of the pipe. Electricity conduction coefficient The calculation formula is as follows: in, Where the pipe radius is; The fluid conductivity; This refers to the length of the pipe. Based on the obtained hydraulic and electrical conduction coefficients, two sparse matrices are constructed in parallel to achieve the goal of simultaneously solving for hydraulic and electrical pressure. A hydraulic sparse matrix is ​​constructed based on the hydraulic conduction coefficient and the connectivity between adjacent nodes. For each internal node in the sparse matrix, the flow balance equation is: in, Let ik be the hydraulic conductivity at point ik; and The pressure at points i and k; Based on the boundary conditions at the inlet and outlet nodes, the LU decomposition method is used to efficiently solve the sparse matrix and obtain the hydraulic pressure of each node. The electric potential at each node is also calculated by constructing the conductance matrix based on the electric conduction coefficient and the connectivity between adjacent nodes. Based on the boundary conditions at the inlet and outlet nodes, the LU decomposition method is used to solve the linear equation system and obtain the potential of each node.

9. The parallel calculation method for optimal hydraulic and electrical paths in rock pore networks as described in claim 1, characterized in that, The process of using the hydraulic and electrical flux of each pipe as the weight of the corresponding pipe and selecting the hydraulic and electrical flux path with the largest path weight as the optimal hydraulic and electrical path for the rock sample includes: using the Dijkstra algorithm to perform optimal path analysis on the three-dimensional pore network model, and simultaneously obtaining the optimal hydraulic and electrical paths through parallel computing. Let the set of nodes in the three-dimensional porous network model be . The pipeline assembly is Source The starting point and ending point of the path For the endpoint of the path, the weight of the pipeline between adjacent nodes u and v in each path. Defined as the absolute value of the water pressure difference and electric potential at both ends of the pipe connecting node u and node v, i.e.: in, , These are the water pressure and potential values ​​at node u at one end of the pipeline and node v at the other end of the pipeline, respectively. The greater the pressure difference between the two ends of the pipeline, the greater the flow rate in the pipeline, and the higher the weight. In Dijkstra's algorithm, a distance array d[v] is introduced, which represents the current known maximum cumulative path weight from the source point s to the node v. The initialization condition is set to d[s]=0. During the algorithm iteration, for the current node u, traverse all its adjacent nodes v, if the following conditions are met... If the sum of the path weights from the source point through u to v is greater than the currently recorded d[v], then update d[v] and record the predecessor node for subsequent path backtracking.

10. A parallel computing system for optimal hydraulic and electrical paths in rock pore networks, characterized in that, include: The image acquisition module is configured to acquire two-dimensional grayscale image slices of rock samples and preprocess them; The 3D pore network model construction module is configured to construct a microscale 3D pore network model based on preprocessed 2D grayscale image slices. The boundary condition addition module is configured to add boundary conditions to the microscale three-dimensional porous network model, including inlet and outlet pressure of the hydraulic field, and inlet potential and outlet grounding of the electric field. The hydraulic and electric conduction coefficient calculation module is configured to solve the hydraulic and electric conduction coefficients of a single pipe based on boundary conditions, and obtain the water pressure and electric potential at each node, as well as the hydraulic flux and electric flux of each pipe. The optimal path determination module is configured to use the hydraulic flux and electrical flux of each pipe as the weight of the corresponding pipe, and select the hydraulic and electrical flux path with the largest path weight as the optimal hydraulic and electrical path of the rock sample. The path weight of the optimal path is the sum of the weights of all pipes it passes through.

Images