A method for restricted random walk simulation to calculate flow tortuosity in porous media

Abstract

The present invention discloses a method for simulating restricted migration to calculate the flow tortuosity in a porous medium, including simplifying and correcting a sample digital core based on the watershed algorithm; including the collision probability of microscopic-scale continuous-phase particles characterized by the gray value of the digital core in different regions and different migration directions; including the migration direction of the continuous-phase operator and the digital core mirror extension mapping method for different digital core resolutions; quantifying the tortuosity in different states by simulating the stable diffusion coefficient under infinite mirror mapping and combining the operator migration mode in the real situation in the above text. Based on the digital core, the tortuosity data obtained by the method of simulating random migration in the porous medium in the present invention, compared with the existing methods of ideal mathematical models such as koponen, sousa, and barrande, introduces the complex pore structure and force influence in migration. The tortuosity obtained by the new method is simulated and calculated, and the coincidence rate of the result is increased by 11% compared with the tortuosity of the traditional method.

Description

Technical Field

[0001] The present invention relates to the technical field of gas reservoir development, particularly to the field of unconventional oil and gas development, and specifically to a method for simulating restricted migration to calculate the flow tortuosity in porous media. Background Art

[0002] As an important parameter for describing the seepage channel, the pore tortuosity is defined as the ratio of the actual length of a specified migration during the seepage process to the macroscopic length of the seepage channel. With the development of unconventional oil and gas reservoirs in recent years, the oil and gas seepage process has become more complex. There are large errors in using the tortuosity obtained by the conventional method of simply calculating through the pore structure to explain the seepage process.

[0003] In recent years, experts and scholars have tended to use high-precision CT inversion to make digital cores to simulate the migration in complex unconventional oil and gas reservoirs. However, the interpretation of tortuosity only stays at the quantification of different average particle sizes of the pore medium minerals, ignoring the influence of the inter-particle force and the anisotropy of the reservoir rock during the actual seepage process. Summary of the Invention

[0004] In view of the above problems, the present invention aims to provide a technique for accurately characterizing the tortuosity method of complex porous media, making its results closer to the actual seepage situation.

[0005] The technical solution of the present invention is as follows:

[0006] A method for simulating restricted migration to calculate the flow tortuosity in porous media, characterized by including the following steps:

[0007] Step1: Perform a 360-degree scan on the rock sample and conduct three-dimensional correction processing based on the multi-angle images of the rock sample.

[0008] Step2: Separate the pores and the rock skeleton of the rock sample from the three-dimensional gray data of the processed rock sample through image threshold segmentation. Compare the matrix statistics of the pores with the actual measured porosity to adjust the threshold data and obtain the data characterizing the real rock skeleton and pores.

[0009] Step3: Combine with on-site experiments, use the watershed algorithm to segment the Euclidean distance matrix of the void phase, convert the complex matrix into a relatively simple set of corresponding sample space connection points, and set this matrix area as the collision-free area in the random walk.

[0010] Step4: Continuously remove the boundary voxels of the rock void phase, and combine with on-site experiments to obtain the collision and rebound probabilities of each voxel of the corresponding rock sample at the corresponding scan resolution.

[0011] Step 5: Select any point inside the pores in the sample digital core matrix and perform a random walk simulation. During the walk, record the corresponding walk distance for the corresponding number of steps, and perform mirror mapping on the matrix boundary during the walk. Different measures are taken to limit the number of steps or the walk distance when the random walk simulation coordinates touch the digital core skeleton.

[0012] Step 6: Characterize the tortuosity of the sample reservoir through the relationship between the random walk steps and the distance.

[0013] Further, the specific process of the said Step 1 is as follows:

[0014] Use a computer to traverse the reconstructed three-dimensional digital core matrix. Taking the gray data volume as the standard, divide the gray data into two frequency bands. For the data after the frequency band division, in the parts that are difficult to distinguish, based on the reconstruction algorithm of the corresponding gray value and density, convert the gray value into density data and replace its gray value data.

[0015] Further, the specific process of the said Step 2 is as follows:

[0016] Step 201: Measure the porosity of the sample through rock sample porosity measurement methods such as the liquid saturation drainage method or the helium method.

[0017] Step 202: Based on the corrected three-dimensional digital core data, divide the skeleton and pores by setting a gray or density threshold. Traverse the binarized pores after the division, perform the bwlabel algorithm on the coordinates marked as pores in the matrix to calculate their independence, mark and calculate the volume occupied by the non-connected pores after adjusting the threshold. Compare the porosity φ1 obtained from the experiment in Step 201, the porosity φ2 and the non-connected porosity φ after adjusting the threshold. u By continuously adjusting the threshold, make φ1 = φ2 - φ u , so as to determine the data that can reflect the true pore and skeleton structure of the core.

[0018] Further, the specific process of the said Step 3 is as follows:

[0019] Perform Gaussian smoothing on the gray values in the areas divided as pores in the digital core according to the process of Step 2, erase the minimum values, then adjust the pore gray threshold in the watershed algorithm, divide the digital core into multiple similar small intervals according to the gray values and number them according to their gray value sizes. The voxel radius R and its corresponding threshold Th are stored in the matrix at the corresponding [x, y, z] positions of the corresponding digital core.

[0020] Based on the influence of the continuous phase velocity on the collision phase in the Boltzmann equation:

[0021]

[0022] Where: f is a dimensionless external force, X α is the corresponding position, ξ α is the dimensionless particle diffusion velocity

[0023] Considering that the tortuosity calculation in porous media ignores the external force, and combining the equations for the velocity field and the collision phase in the Boltzmann equation, the diffusion distribution probability of a pore point in the digital core can be equivalent to:

[0024]

[0025] Considering that the pore network in the digital core can be equivalent to an elliptical tube, so combining the Poiseuille equation in the flow process and the independent pores separated by the watershed algorithm in Step 3, the distribution probability of each voxel in the rock sample during the flow process can be expressed by the following equation:

[0026]

[0027] After simplification, the probability in each direction when a single voxel in the core pores diffuses freely at the microscale is:

[0028]

[0029] Furthermore, the specific process of Step 5 is as follows:

[0030] Step501: Based on the pore network model obtained in Step 2, randomly generate an operator in the pore network, read the data matrix obtained in Step 3, and combine the calculation in Step 4 to calculate the migration probability of the N neighboring voxels of the operator at the next iteration time step. If the resolution of the digital core is low, that is, the size of a single voxel is large, due to precision limitations, the model cannot accurately represent the pore structure, so there may be pore structures below the precision. Consider 27 neighboring voxels as the migration directions; if the resolution of the digital core is high and can represent the fine pore structure, then take 14 or 6 directions around a single voxel. The corresponding voxel migration probability is:

[0031]

[0032] When the operator migrates to the skeleton position in the digital core matrix during the random walk process, the position of the operator remains unchanged, and the probability of the current migration position is ignored in the next iteration migration process.

[0033] Step502: After obtaining the different direction probabilities in the digital core matrix through the Step501 process, call the built-in random number library in the programming software to make the operator perform a random walk with the corresponding direction probability within the corresponding pore position in the digital core matrix; perform horizontal and vertical mirror mappings on the digital core matrix to expand the boundary, and the corresponding process of the boundary mirror mapping is:

[0034]

[0035] Wherein: x i+1 , y i+1 represent the coordinates within the matrix where the operator is located before mapping; REF(x, y) represents the coordinates of the digital core matrix after mapping; b is the size of the matrix boundary.

[0036] Step503: Use programming software to repeatedly iterate the process of Step501 - Step502. At the same time, perform a migration simulation in an unrestricted space once, that is, all within the matrix are pores; design and record the positions of the operator at different time steps according to different needs.

[0037] Furthermore, the sub - process of selecting different step sizes according to different needs in the Step503 process to obtain the specific positions of different operators can be explained as follows: If the research process is two - dimensional planar migration, then record the distance r at the corresponding two - dimensional coordinates at different iteration steps according to the iteration process; if considering the overall tortuosity, then record the distance r at the corresponding three - dimensional coordinates at different iteration steps according to the iteration process. If considering the tortuosity at different migration scales, then combine the corresponding digital core resolution to record the corresponding steps when the operator reaches the specified migration distance.

[0038] On the basis of ensuring easy implementation, the present invention is more practical than the existing methods in screening the main controlling factors affecting the tight gas production, and has extremely profound significance for the subsequent research on tight gas production prediction. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the following drawings are only some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can be obtained based on these drawings.

[0040] Figure 1 It is a schematic flowchart of a method for restricted random walk simulation to calculate the flow tortuosity in a porous medium according to the present invention;

[0041] Figure 2 It is the digital core matrix after visualization by Avizo calculated in Example 1;

[0042] Figure 3 It is the matrix data (Avizo visualization) calculated in Example 1 based on the calculation and processing by the watershed algorithm;

[0043] Figure 4The diffusion coefficients at different time steps were calculated for Example 1; Detailed implementation manners

[0044] The present invention will be further described below in conjunction with the accompanying drawings and embodiments. It should be noted that, without conflict, the embodiments in the present application and the technical features in the embodiments may be combined with each other. It should be pointed out that unless otherwise specified, all technical and scientific terms used in the present application have the same meaning as commonly understood by those of ordinary skill in the technical field to which the present application belongs. As used in the present disclosure, the words such as "including" or "comprising" and the like are intended to mean that the elements or items appearing before the word cover the elements or items listed after the word and their equivalents, without excluding other elements or items.

[0045] Example 1

[0046] As shown in the attached Figure 1 figures, a method for accurately measuring the dynamic tortuosity of seepage flow in porous media based on random walk includes the following steps:

[0047] Step1: Perform a 360-degree scan on the rock sample and perform three-dimensional correction processing based on the multi-angle images of the rock sample; the converted three-dimensional digital core matrix is shown in Table 1, and the data read by Avizo visualization is as shown in the attached drawings of the specification. Figure 2 figures.

[0048] Table 1 Three-dimensional digital core data matrix after conversion

[0049]

[0050]

[0051] Step2: Separate the pore space and the rock sample skeleton of the processed rock sample three-dimensional gray data through image threshold segmentation. Compare the matrix statistics of the pores with the actual measured porosity to adjust the threshold data and obtain the data characterizing the real rock sample skeleton and pores;

[0052] After adjusting the threshold, the porosity of 0.143 obtained with the threshold set at 4817 is the closest to the test result. Therefore, the pores are divided with this threshold, Figure 2 which is the result after Avizo visualization of the processed matrix data.

[0053] Step3: In combination with on-site experiments, use the watershed algorithm to segment the Euclidean distance matrix of the void phase, convert the complex matrix into a relatively simple set of corresponding sample space connection points, and set this matrix area as the collision-free area in random walk.

[0054] The collision-free area of the random walk after the watershed algorithm and correction is as shown in Figure 3 the three blue positions.

[0055] Step 4: By continuously removing the voxel at the boundary of the rock pore phase and combining with on-site experiments, the collision and rebound probabilities of individual voxels of the corresponding rock sample at the corresponding scanning resolution can be obtained.

[0056] Taking the digital core matrix coordinates [77, 77, 77] as an example, based on its gray value, the migration probabilities in 9 directions are calculated as 163 / 814, 5 / 906, 16 / 103, 49 / 370, 5 / 56, 25 / 222, 24 / 299, 50 / 841, 35 / 213.

[0057] Step 5: Select any point within the pores of the digital core matrix of the sample and perform a random walk simulation. During the walk, record the corresponding number of steps and the corresponding migration distance, and perform mirror mapping on the matrix boundary during the walk. Different measures are taken to limit the number of steps or the migration distance when the random walk simulation coordinates touch the digital core skeleton.

[0058] There are different particle sizes in the digital core sample. The digital core sample is divided into large particle size, small particle size populations, and 10,000 * 100,000-step walks are carried out separately according to the unconstrained space, and the migration distances in different cases are recorded.

[0059] Step 6: Characterize the tortuosity of the sample reservoir through the relationship between the random walk steps and the distance.

[0060] As Figure 4 shown, when the number of walk steps is higher than 70,000 times, the diffusion coefficients in four different cases all tend to be constant values. Therefore, the tortuosities of the large particle size region, small particle size region, and the whole of this sample are 1.7, 2.1, and 1.9 respectively. Substitute this value into the Poiseuille equation for flow simulation, and the simulation result differs from the real result by 8%. Compared with the tortuosity value calculated by porosity, the simulation error is reduced by 11%.

[0061] The above is only a preferred embodiment of the present invention, and it does not impose any form of limitation on the present invention. Although the present invention has been disclosed above with the preferred embodiment, it is not intended to limit the present invention. Any person skilled in the art can make some changes or modifications to it as equivalent embodiments within the scope of the technical solution of the present invention. However, as long as it does not depart from the content of the technical solution of the present invention, any simple modification, equivalent change, and modification made to the above embodiments based on the technical essence of the present invention still fall within the scope of the technical solution of the present invention.

Claims

1. A method for simulating restricted random walks to calculate the flow tortuosity in porous media, characterized by the following steps: Step1: Perform a 360-degree scan on the rock sample and conduct three-dimensional correction processing based on multi-angle images of the rock sample; Step2: Separate the pores and the rock skeleton of the rock sample from the three-dimensional gray-scale data of the processed rock sample through image threshold segmentation, compare the matrix statistics of the pores with the actual measured porosity to adjust the threshold data, and obtain the data characterizing the real rock skeleton and pores; Step3: Combine field experiments, use the watershed algorithm to segment the Euclidean distance matrix of the void phase, convert the complex matrix into a relatively simple set of connected points in the corresponding sample space, and set this matrix area as the collision-free area in the random walk; Step4: Continuously remove the boundary voxels of the rock void phase, and combine field experiments to obtain the collision and rebound probabilities of individual voxels of the corresponding rock sample at the corresponding scanning resolution; Step5: Select any point in the pores of the sample digital core matrix and perform a random walk simulation. Record the corresponding number of steps and the corresponding walking distance during the walking process. Perform mirror mapping on the matrix boundary during the random walk simulation. Different measures are taken to limit the number of steps or the walking distance when the random walk simulation coordinates contact the digital core skeleton; Step6: Characterize the tortuosity of the sample reservoir through the relationship between the random walk steps and the distance. The tortuosity is determined by the slope of the image of the variance of the average migration distance and the number of steps. The corresponding tortuosity is: Where: d diff is the diffusion parameter for free migration, d lmtd is the diffusion coefficient for restricted migration, is the corresponding image slope.

2. A method for simulating restricted random walks to calculate the flow tortuosity in porous media according to claim 1, characterized in that the specific process of Step1 is: Use a computer to traverse the reconstructed three-dimensional digital core matrix. Based on the gray-scale data volume, divide the gray-scale data into two frequency bands. For the data after frequency band division, in the parts that are difficult to distinguish, based on the reconstruction algorithm of the corresponding gray value and density, convert the gray value into density data and replace its gray value data.

3. A method for simulating restricted random walks to calculate the flow tortuosity in porous media according to claim 2, characterized in that the specific process of Step2 is: Step201: Measure the porosity of the sample through the liquid saturation drainage method or the helium method; Step202: Based on the corrected three-dimensional digital core data, divide the skeleton and pores by setting the gray level or density threshold. Traverse the binarized pores after division, perform the bwlabel algorithm on the coordinates marked as pores in the matrix to calculate their independence, mark and calculate the volume occupied by the non-connected pores after adjusting the threshold, compare the porosity φ1 obtained from the experiment in Step201, the porosity φ2 and the non-connected porosity φ after adjusting the threshold u , and continuously adjust the threshold so that φ1 = φ2 - φ u , so as to determine the data that can reflect the true pore and skeleton structure of the core.

4. A method for simulating restricted random walks to calculate the flow tortuosity in porous media according to claim 3, characterized in that the specific process of Step3 is: Perform Gaussian smoothing on the gray values in the areas divided into pores in the digital core according to the process of Step2, erase the minimum values, then adjust the pore gray threshold in the watershed algorithm, divide the digital core into multiple similar small intervals according to the gray values and number them according to their gray value sizes. The voxel radius R and its corresponding threshold Th are stored in the matrix at the corresponding [x, y, z] positions of the corresponding digital core.

5. A method for simulating restricted random walks to calculate the flow tortuosity in porous media according to claim 4, characterized in that the specific process of Step5 is: Step501: Based on the pore network model obtained in Step 2, a random operator is generated in the pore network. The data matrix obtained in Step 3 is read, and the migration probabilities of the surrounding N voxels of the random voxel in the next iteration time step are calculated in combination with the calculation in Step 4. If the resolution of the digital core is low, that is, the size of a single voxel is large, due to precision limitations, the model cannot finely represent the pore structure. Therefore, there may be pore structures below the precision. Consider the surrounding 27 voxels as the migration directions; if the resolution of the digital core is high and can characterize the fine pore structure, then take 14 or 6 directions around a single voxel, and the corresponding voxel migration probabilities are as follows: When the operator migrates to the skeleton position in the digital core matrix during the random walk process, the position of the operator remains unchanged, and the probability of the current migration position is ignored in the next iteration migration process; Step502: Obtain the probabilities in different directions obtained from the digital core matrix through the process of Step501, and call the built-in random number library in the programming software to make the operator perform a random walk with the corresponding direction probabilities within the corresponding pore positions in the digital core matrix; perform horizontal and vertical mirror mappings on the digital core matrix to expand the boundary. The process of the corresponding boundary mirror mapping is: Wherein: x i+1 , y i+1 represent the coordinates within the matrix where the operator is located before mapping; REF(x, y) represents the coordinates of the digital core matrix after mapping; b is the matrix boundary size; Step503: Use the programming software to repeatedly iterate the processes of Step501 - Step502, and at the same time run a migration simulation in an unrestricted space, that is, all within the matrix are pores; design and record the positions of the operator at different time steps according to different needs.

6. A method for restricted random walk simulation to calculate the flow tortuosity in a porous medium according to claim 5, characterized in that the sub - process of obtaining the specific positions of different operators by selecting different step sizes according to different needs in the process of Step503 can be explained as follows: If the research process is planar two - dimensional migration, then record the distance r at the corresponding two - dimensional coordinates at different iteration steps according to the iteration process; if considering the overall tortuosity, then record the distance r at the corresponding three - dimensional coordinates at different iteration steps according to the iteration process; if considering the tortuosity at different migration scales, then record the corresponding number of steps when the operator reaches the specified migration distance in combination with the corresponding digital core resolution.

Images