Method for obtaining representative volume element of rock

Abstract

The invention relates to a method for obtaining a representative volume element of a rock, which comprises the following steps of: converting a three-dimensional grey-scale map of the rock into a three-dimensional binary image, inwards taking eight cube-shaped sub-blocks from eight vertexes of the three-dimensional binary image, starting from one voxel of the side length of each sub-block, adding 1 every time, respectively calculating the porosity, and when the difference value of the porosity of any two sub-blocks is smaller than a set difference value, calculating the representative volume element of the rock. Performing mercury injection simulation on the eight sub-blocks: performing corrosion and expansion operation on the three-dimensional binary images of the sub-blocks by using a spherical structure element with the radius of R, and finding pores which are communicated with the outside and are larger than R in size; subtracting 1 from R every time, and repeating mercury injection simulation to obtain eight capillary pressure-mercury injection saturation curves; when the eight curves converge into one, the size of the sub-block is the size of the representative volume unit of the rock; otherwise, continuing mercury injection simulation after adding one voxel to the side length of the sub-block every time; according to the method, the method for obtaining the representative volume element of the rock is optimized, so that the obtained representative volume element of the rock is more accurate.

Description

(1) Technical Field:

[0001] The present invention relates to a method for rock detection, and particularly to a method for obtaining a representative volume element of a rock. (2) Background Art:

[0002] A representative elementary volume (REV for short) refers to the smallest volume element that can represent the macroscopic properties of a porous medium when describing the properties of the porous medium. Within this volume, the microscopic structural characteristics of the porous medium, including the size, shape, and distribution of pores, can be fully reflected, so that the properties of this volume element can represent the properties of the entire porous medium. This volume element should be large enough so that the measured properties (such as porosity, permeability, fractal dimension, etc.) are statistically representative within this volume element and do not fluctuate too much due to being too small. Porosity is closely related to the REV. The selection of the size of the REV should be much larger than that of a single pore or particle, so as to contain enough pores and matrix to average out the microscopic inhomogeneity. When measuring porosity, if the selected volume element is too small, for example, just selecting a place with only solid or only pores, the resulting porosity will be 0 or 1, which cannot contain enough pore information and leads to inaccurate measurement results of porosity. Gradually increasing the selected volume element, the measured value of porosity will tend to be stable, and the volume element at this time is the REV; after that, even if the volume is further increased, the porosity will not change significantly.

[0003] The prior art uses the porosity calculation method for rock REV analysis; the rock porosity calculation method is based on the three-dimensional gray-scale image of the rock obtained by CT scanning. Starting from the 8 vertices of the image, the size of the sub-blocks is gradually increased inward. Each time the side length of the 8 sub-blocks is increased, the ratio of the voxel of the pores in each sub-block to the total voxel of each sub-block is calculated to obtain the porosity. When the porosity tends to be stable, the volume element at this time is the REV of the rock. This rock porosity calculation method for obtaining the REV of the rock has a simple calculation process and high speed, but the calculation result is relatively rough, unable to distinguish the distribution of pores of different sizes inside the rock, and it is difficult to obtain the true REV. (3) Summary of the Invention:

[0004] The technical problem to be solved by the present invention is: to provide a method for obtaining a representative volume element of a rock, which overcomes the problems of relatively rough and inaccurate calculation results in the prior art porosity calculation method, adds a mercury intrusion simulation operation, and uses the permeability of the rock as a measurement index for obtaining the representative volume element of the rock, so that the obtained representative volume element of the rock is more accurate.

[0005] The technical solution of the present invention:

[0006] A method for obtaining a representative volume element of a rock, comprising the following steps:

[0007] Step 1: Scan the rock to obtain a three-dimensional gray-scale image of the rock;

[0008] Step 2: Convert the three-dimensional gray-scale image into a three-dimensional binary image according to a threshold. In the three-dimensional binary image, void voxels are labeled as 1, matrix voxels are labeled as 0, and the three-dimensional binary image is cube-shaped;

[0009] The CT scan image of the rock is a three-dimensional gray-scale image, and its mathematical essence is a three-dimensional matrix; the value of each voxel in the three-dimensional gray-scale image corresponds to the value of the element at the corresponding position in the three-dimensional matrix; taking an 8-bit gray-scale image as an example, the range of any voxel value in the three-dimensional gray-scale image is 0 to 2 8 -1, that is, (0 to 255), and the specific value of the voxel is determined by the composition of the corresponding position of the physical entity; the density of the void part of the physical entity is the smallest, and the voxel value at the corresponding position in the CT scan gray-scale image is close to 0 and appears black in the image; the density of the matrix part of the physical entity is relatively high, and the voxel value at the corresponding position in the CT scan gray-scale image is close to 255 and appears bright in the image; based on the above physical facts, a critical value d can be determined. In the three-dimensional gray-scale image, all voxels with values in the range of (0, d] are identified as voids, and each void voxel is labeled with 1; while all voxels with values in the range of (d, 255] are identified as the matrix, and each matrix voxel is labeled with 0; in this way, the three-dimensional gray-scale image is transformed into a three-dimensional binary image through threshold segmentation. In the three-dimensional binary image, there are only two colors, black for the matrix and white for the pores;

[0010] Step 3: Take 8 cube-shaped sub-blocks from the 8 vertices of the three-dimensional binary image of the rock towards the inside of the three-dimensional binary image. The side lengths of the 8 sub-blocks are all n, and let n = 1 voxel;

[0011] Step 4: Calculate the porosity of each sub-block. The calculation method of the porosity is as follows:

[0012]

[0013] In the formula, K i is the porosity of the i-th sub-block, and i is a natural number from 1 to 8;

[0014] Step 5: When the difference in porosity between any two of the 8 sub-blocks is less than the set difference, it is considered that the porosity converges. At this time, the side length n of the sub-block = n s , and step 6 is executed; otherwise, the side length of the sub-block is increased by 1 voxel, that is, n = n + 1, and step 4 is executed;

[0015] The side length of the sub-block starts from n = 1 and increases by 1 each time. The porosity K of each sub-block is calculated respectively. The (K, n) of 8 sub-blocks are plotted in a rectangular coordinate system, and the points of (K, n) of each sub-block are connected to draw the porosity-side length curve of each sub-block;

[0016] On the porosity-side length curve, as the side length of the sub-block increases, the porosity-side length curves of 8 sub-blocks gradually approach until they entangle. When they converge to a single line, the side length position at this time is used as the rough convergence point of the porosity, and the side length of the sub-block at the convergence point is used as the starting size for mercury intrusion simulation; The porosity at this position is only used as a basis for the starting size of mercury intrusion simulation. For example, when the porosity converges, then n s ×n s ×n s is the starting size for mercury intrusion simulation of the sub-block. Starting from the 8 vertices of the three-dimensional binary image, starting from n s ×n s ×n s At this position, mercury intrusion simulation operations are performed on each sub-block;

[0017] Step 6: Perform mercury intrusion simulation operations on 8 sub-blocks respectively. The mercury intrusion simulation operation method for each sub-block is as follows:

[0018] In the mercury intrusion simulation operation, erosion operation and dilation operation are used. The basic operation units of erosion operation and dilation operation are spherical structuring elements (SSE); SSE is a discrete sphere assembled by a group of voxels. The coordinates of the central voxel of the discrete sphere are defined as (0, 0, 0), and the coordinates of any other voxel t are (x t , y t , z t ). The distance of any voxel t from the central voxel is:

[0019]

[0020] dist(t, 0) satisfies the relationship dist(t, 0) ≤ R, where R is the integer radius of SSE; SSE is stored in the set Se;

[0021] Step 6.1: Set a parameter R, and the initial value of R is a non-zero natural number;

[0022] Step 6.2: Use a solid spherical structuring element with a radius of R to perform erosion operation on the three-dimensional binary image of the sub-block. The purpose of the erosion operation is to convert the void voxels in the three-dimensional binary image of the sub-block into solid voxels; All solid spherical structuring elements and the three-dimensional binary image of the sub-block perform a logical AND operation to convert the void voxels in the three-dimensional binary image of the sub-block whose distance from the solid spherical structuring element is less than or equal to R into solid voxels;

[0023] Step 6.3: If all the void voxels have become solid voxels, it indicates that there are no pores with a size larger than R, then go to Step 6.6; otherwise, go to Step 6.4;

[0024] Step 6.4: Use a void spherical structural element with a radius of R to perform a dilation operation on the three-dimensional binary image that has been eroded in Step 6.2. The purpose of the dilation operation is to convert the solid voxels in the three-dimensional binary image into void voxels; perform a logical OR operation on all the void spherical structural elements and the three-dimensional binary image that has been eroded in Step 6.2, and convert the solid voxels within a distance less than or equal to R from the void spherical structural element in the three-dimensional binary image that has been eroded in Step 6.2 into void voxels;

[0025] Since there are still some residual voids in the pore space with a size larger than R after the erosion operation, therefore, after performing the dilation operation on the remaining pore space, the pore space with a size larger than R returns to its original form;

[0026] After performing the erosion operation and the dilation operation on the pore space in the three-dimensional binary image through a spherical structural element with a radius R, the pore space with a size less than or equal to R is converted into a solid matrix, thus finding the pores with a size larger than R in the sub-block;

[0027] Step 6.5: When mercury enters the pore space of the core as a non-wetting phase, the capillary pressure is a resistance, and an external mercury injection pressure needs to be applied to overcome the capillary pressure; according to the Young-Laplace equation, the capillary pressure P C (R) of the pores during the mercury injection simulation operation of the sub-block and the void radius (i.e., the radius R of the spherical structural element) have the following relationship:

[0028]

[0029] Since the pore wall is rough, this will affect the intrusion of mercury. According to the Wenzel model,

[0030] cosθ = fcosθ *

[0031] Finally, calculate the capillary pressure P C (R) of the pores in the sub-block:

[0032]

[0033] In the formula, f is the geometric roughness coefficient of the pore wall of the sub-block, f is the ratio of the actual area of the rough surface to its apparent area, σ is the surface tension of mercury, θ is the actual contact angle of mercury on the rough surface, and θ * is the ideal contact angle of mercury;

[0034] Assume that the surface tension and contact angle of mercury are constant during the mercury injection process, then there is a negative correlation between the capillary pressure and the pore radius; when the mercury injection pressure is low, only a small capillary pressure can be overcome to press mercury into pores with a large radius; when the mercury injection pressure increases, a larger capillary pressure can be overcome to press mercury into pores with a small radius.

[0035] Iterate from large to small for R to simulate the mercury injection process with the injection pressure increasing from low to high. A larger pore radius corresponds to a smaller capillary pressure, and a smaller pore radius corresponds to a larger capillary pressure.

[0036] Calculate the mercury injection saturation S(R): The mercury injection saturation S(R) is the ratio of the total mercury injection volume V to the number of voxels in the entire pore space of the three-dimensional binary image of the sub-block. The total mercury injection volume V is the number of voxels of the pore space where mercury invades the sub-block, that is, the number of voxels of the pore space with a size larger than R and connected to the outside of the sub-block in the sub-block.

[0037] Step 6.6: Let R = R - 1. If R = 0, then end the operation; otherwise, go to Step 6.2 and perform mercury injection simulation using a solid spherical structure element with a radius reduced by 1.

[0038] Step 7: First, plot the relationship curves of the capillary pressure P C (R) and the mercury injection saturation S(R) for 8 sub-blocks in the rectangular coordinate system respectively. Then, uniformly take M mercury injection saturation values on the mercury injection saturation S(R) number axis, where M is a natural number greater than or equal to 2. Each mercury injection saturation value corresponds to 8 capillary pressure values on the capillary pressure P C (R) number axis; Arbitrarily select two relationship curves from the 8 relationship curves, calculate the differences between the 2 capillary pressure values corresponding to each mercury injection saturation value to obtain M differences, and take the average of these M differences.

[0039] Step 8: When the average value obtained regardless of which two relationship curves are arbitrarily selected is less than the set average value, it is considered that the relationship curves of the capillary pressure P C (R) and the mercury injection saturation S(R) for the 8 sub-blocks converge, and the size n×n×n of the sub-block at this time is the size of the representative elementary volume REV of the rock; otherwise, let the side length n of the sub-block be n = n + 1, then go to Step 6, and repeat the mercury injection simulation operation on a larger sub-block until the relationship curves of the capillary pressure P C (R) and the mercury injection saturation S(R) converge together.

[0040] In Step 1, a three-dimensional gray-scale image of the rock is obtained using CT scanning.

[0041] In Step 5, the set difference is 0.001.

[0042] In Step 6.1, the method for determining the initial value of the parameter R is as follows:

[0043] Step 6.1.1: Set a parameter R′, where R′ is a non - zero natural number and R′ is greater than the size of the largest pore in the three - dimensional binary image of the sub - block;

[0044] Step 6.1.2: Use a solid spherical structural element with a radius of R′ to perform an erosion operation on the three - dimensional binary image of the sub - block. Perform a logical AND operation on all solid spherical structural elements and the three - dimensional binary image of the sub - block, and convert the void voxels in the three - dimensional binary image of the sub - block whose distance from the solid spherical structural element is less than or equal to R′ into solid voxels;

[0045] Step 6.1.3: If all void voxels have become solid voxels, let R′ = R′−1, and execute Step 6.1.2; otherwise, determine the initial value of the parameter R as R = R′ + 1, and then end.

[0046] In Step 6.5, the method for determining the total mercury intrusion volume V is as follows:

[0047] Step 6.5.1: First, define the boundary void voxels and internal void voxels of the three - dimensional binary image after the dilation operation in Step 6.4, and search for the boundary void voxels;

[0048] Step 6.5.2: Starting from each boundary void voxel in turn, find all the internally connected void voxels connected to each boundary void voxel; in the three - dimensional binary image after the dilation operation in Step 6.4, the total number of voxels of all boundary void voxels and all internally connected void voxels is the total mercury intrusion volume V;

[0049] The method for finding all the internally connected void voxels connected to a boundary void voxel is as follows:

[0050] Step 6.5.2.1: First, search for the internal void voxels connected to the upper, lower, left, right, front, and back six faces of this boundary void voxel; if not found, end the operation, otherwise, mark the found internal void voxels as internally connected void voxels, and use these found internally connected void voxels as the starting voxels for the next search, and then execute Step 6.5.2.2;

[0051] Step 6.5.2.2: Search for the internal void voxels connected to the upper, lower, left, right, front, and back six faces of the starting voxel; if not found, end the operation, otherwise, mark the found internal void voxels as internally connected void voxels, and use these found internally connected void voxels as the starting voxels for the next search, and then execute Step 6.5.2.2.

[0052] In Step 8, set the average value to 0.001.

[0053] Advantages of the present invention:

[0054] 1. The present invention optimizes the method for obtaining a representative volume element of rock. First, the convergence point of the porosity of the sub-blocks is determined by the porosity calculation method, and then mercury intrusion simulation is performed on the sub-blocks starting from this convergence point position, so as to obtain the representative volume element of the rock. The present invention overcomes the problem that the calculation results of the porosity calculation method in the prior art are relatively rough and inaccurate, adds mercury intrusion simulation operations, and uses the permeability of the rock as a measurement index for obtaining the representative volume element of the rock through mercury intrusion simulation operations, making the obtained representative volume element of the rock more accurate. (IV) Description of the drawings:

[0055] Figure 1 Schematic diagram of the conversion of a three-dimensional gray-scale image of rock into a three-dimensional binary image;

[0056] Figure 2 Schematic diagram of taking 8 cube-shaped sub-blocks in the three-dimensional binary image of rock;

[0057] Figure 3 Schematic diagram of the porosity K-side length n curve of 8 sub-blocks;

[0058] Figure 4 (a) Schematic diagram of a spherical structuring element (SSE) with R = 1;

[0059] Figure 4 (b) Schematic diagram of a spherical structuring element (SSE) with R = 2;

[0060] Figure 4 (c) Schematic diagram of an XY cross-section of the set Se;

[0061] Figure 5 (a) Schematic diagram of erosion operation;

[0062] Figure 5 (b) Schematic diagram of the result of erosion operation;

[0063] Figure 5 (c) Schematic diagram of dilation operation;

[0064] Figure 5 (d) Schematic diagram of the result of dilation operation;

[0065] Figure 6 Schematic diagram of the detailed process from erosion operation to result;

[0066] Figure 7 Schematic diagram of the detailed process from dilation operation to result;

[0067] Figure 8 Schematic diagram of the method for determining the total mercury injection volume V;

[0068] Figure 9 Schematic diagram of the pore space of the mercury intrusion sub-block;

[0069] Figure 10 is the capillary pressure P of 8 sub-blocks C Schematic diagram of the (R)-mercury intrusion saturation S(R) curve;

[0070] Figure 11 is the capillary pressure P of 8 sub-blocks C Schematic diagram of the convergence of the (R)-mercury intrusion saturation S(R) curve.

[0071] In the figure: 3 is the sub-block, 4 is the three-dimensional gray-scale image of the rock, 5 is the three-dimensional binary image of the rock, 6 is the matrix of the rock, 7 is the pore of the rock; P1, P2, P3 and P4 are voids, T1, T2 and T3 are throats; a is the boundary void voxel, b is the internal connected void voxel. (V). Specific implementation method:

[0072] The method for obtaining the representative volume element of the rock includes the following steps:

[0073] Step 1: Use CT to scan the rock to obtain the three-dimensional gray-scale image 4 of the rock, as Figure 1 shown on the left;

[0074] Step 2: Convert the three-dimensional gray-scale image 4 into a three-dimensional binary image 5 according to the threshold, as Figure 1 shown on the right. In the three-dimensional binary image 5, the void voxels are marked as 1, and the matrix voxels are marked as 0. The three-dimensional binary image 5 is cube-shaped;

[0075] The rock CT scan image is the three-dimensional gray-scale image 4, and its mathematical essence is a three-dimensional matrix; the value of each voxel in the three-dimensional gray-scale image 4 corresponds to the value of the corresponding element in the three-dimensional matrix; taking an 8-bit gray-scale image as an example, the range of any voxel value in the three-dimensional gray-scale image 4 is 0 to 2 8 -1, that is, (0 to 255), and the specific value of the voxel is determined by the composition of the corresponding position of the physical entity; the density of the void part of the physical entity is the smallest, and the voxel value at the corresponding position in the CT scan gray-scale image is close to 0 and is black in the image; the density of the matrix 6 part of the physical entity is relatively high, and the voxel value at the corresponding position in the CT scan gray-scale image is close to 255 and is bright in the image; based on the above physical facts, a critical value d can be determined. In the three-dimensional gray-scale image 4, all voxels with values in the range of (0, d] are identified as voids, and each void voxel is marked with 1; while all voxels with values in the range of (d, 255] are identified as the matrix 6, and each matrix voxel is marked with 0; in this way, the three-dimensional gray-scale image 4 is transformed into a three-dimensional binary image 5 through threshold segmentation. In the three-dimensional binary image 5, there are only two colors, black is the matrix 6 of the rock, and white is the pore 7 of the rock;

[0076] Step 3: Take eight cube-shaped sub-blocks 3 from the eight vertices of the three-dimensional binary image 5 of the rock into the interior of the three-dimensional binary image 5, as Figure 2 shown. The side lengths of the eight sub-blocks 3 are all n, and let n = 1;

[0077] Step 4: Calculate the porosity of each sub-block 3. The calculation method of the porosity is as follows:

[0078]

[0079] In the formula, K i is the porosity of the i-th sub-block 3, and i is a natural number from 1 to 8;

[0080] Step 5: Among the eight sub-blocks 3, when the difference in porosity between any two sub-blocks 3 is less than 0.001, it is considered that the porosity converges. At this time, the side length n of the sub-block 3 is 200, and step 6 is executed; otherwise, the side length of the sub-block 3 is increased by 1 voxel, that is, n = n + 1, and step 4 is executed;

[0081] The side length of the sub-block 3 starts from n = 1 and increases by 1 each time. The porosity K of each sub-block 3 is calculated respectively. The (K, n) of the eight sub-blocks 3 are plotted in a rectangular coordinate system, and the points (K, n) of each sub-block 3 are connected to draw the porosity-side length curve of each sub-block 3, as Figure 3 shown;

[0082] On the porosity-side length curve, as the side length of the sub-block 3 increases, the porosity-side length curves of the eight sub-blocks 3 gradually approach until they are entangled and converge to a single line. At this time, the side length position is used as the rough convergence point of the porosity, as Figure 3 the position where the side length n is 200 in. The side length of the sub-block 3 at the convergence point is used as the starting size for mercury intrusion simulation; the porosity at this position is only used as a basis for the starting size of mercury intrusion simulation. For example, if the porosity converges, then 200×200×200 is the starting size for mercury intrusion simulation of the sub-block 3. Starting from the eight vertices of the three-dimensional binary image 5, mercury intrusion simulation operations are performed on each sub-block 3 starting from this position of 200×200×200;

[0083] Step 6: Perform mercury intrusion simulation operations on the eight sub-blocks 3 respectively. The mercury intrusion simulation operation method for each sub-block 3 is as follows:

[0084] The mercury intrusion simulation operation uses erosion operation and dilation operation. The basic operation unit of the erosion operation and dilation operation is the spherical structuring element (SSE); the SSE is a discrete sphere assembled by a group of voxels. The coordinates of the central voxel of the discrete sphere are defined as (0, 0, 0), and the coordinates of any other voxel t are (x t , y t, z t ), the distance of any voxel t from the central voxel is:

[0085]

[0086] dist(t, 0) satisfies the relation dist(t, 0) ≤ R, where R is the integer radius of the SSE; the SSEs with R = 1 and R = 2 are respectively as Figure 4 (a) and Figure 4 (b) shown, the SSEs are stored in the set Se, Figure 4 (c) shows an XY cross-section of Se, and this set stores an SSE with R = 2;

[0087] Step 6.1: Set a parameter R, and the initial value of R is a non-zero natural number;

[0088] Step 6.2: Use a solid spherical ball structuring element with a radius of R to perform an erosion operation on the three-dimensional binary image 5 of sub-block 3, as Figure 5 (a) shown. The purpose of the erosion operation is to convert the void voxels in the three-dimensional binary image 5 of sub-block 3 into solid voxels; all solid spherical ball structuring elements perform a logical AND operation with the three-dimensional binary image 5 of sub-block 3, and convert the void voxels in the three-dimensional binary image 5 of sub-block 3 whose distance from the solid spherical ball structuring element is less than or equal to R into solid voxels, as Figure 5 (b) shown;

[0089] The SSE of the erosion operation is composed of solid voxels, as Figure 6 shown in (a) of, and the purpose is to convert the void voxels in the binary image into solid voxels; as Figure 6 shown in (b) of, the central voxel of Se coincides with the solid voxel marked as j in the binary image, and at the same time Se performs a logical AND operation with the voxels in the corresponding area of the binary image, and converts the void voxels in the binary image whose distance from the solid voxel j is less than or equal to R into solid voxels; in the calculation process, in order to prevent Se from overflowing, a layer of solid voxels needs to be added outside the binary image boundary, and the number of added solid voxel layers is the same as the radius of the SSE; perform the erosion operation on each solid voxel in the binary image, and the result is as Figure 6 shown in (c) of;

[0090] Step 6.3: If all the void voxels have become solid voxels, indicating that there are no pores 7 with a size larger than R, then go to Step 6.6, otherwise go to Step 6.4;

[0091] Step 6.4: Use a void spherical ball structuring element with a radius of R to perform a dilation operation on the three-dimensional binary image 5 eroded in Step 6.2, as Figure 5As shown in (c), the purpose of the dilation operation is to convert the solid voxels in the three-dimensional binary image 5 into void voxels; all void spherical structuring elements are logically ORed with the three-dimensional binary image 5 that has been eroded in step 6.2, and the solid voxels within the three-dimensional binary image 5 that has been eroded in step 6.2 and are at a distance less than or equal to R from the void spherical structuring elements are converted into void voxels, as shown in Figure 5 (d);

[0092] Since there are still some residual voids in the pore space larger than R after the erosion operation, therefore, after performing the dilation operation on the remaining pore space, the pore space larger than R returns to its original form;

[0093] The SSE of the dilation operation consists of void voxels, as shown in Figure 7 (a), and the purpose is to convert the solid voxels in the binary image into void voxels; as shown in Figure 7 (b), the central voxel of Se coincides with the void voxel marked as k in the binary image, and at the same time, Se is logically ORed with the voxels in the corresponding region of the binary image, converting the solid voxels within the binary image that are at a distance less than or equal to R from the void voxel k into void voxels; similar to the erosion operation, a corresponding number of solid voxel layers need to be added outside the boundary of the binary image during the calculation process to prevent Se from overflowing; the dilation operation is performed on each void voxel in the binary image, and the result is as shown in Figure 7 (c);

[0094] After performing the erosion operation and the dilation operation on the pore space in the three-dimensional binary image 5 using the spherical structuring element with radius R, the pore space with a size less than or equal to R is converted into a solid matrix 6, thereby finding the pores 7 in sub-block 3 with a size larger than R;

[0095] Step 6.5. When mercury, as the non-wetting phase, enters the pore space of the core, the capillary pressure is a resistance, and an external mercury injection pressure needs to be applied to overcome the capillary pressure; according to the Young-Laplace equation, when performing the mercury injection simulation operation on sub-block 3, the relationship between the capillary pressure P C (R) of pore 7 and the pore radius (i.e., the radius R of the spherical structuring element) is:

[0096]

[0097] Since the pore wall is rough, which will affect the intrusion of mercury, according to the Wenzel model,

[0098] cosθ = fcosθ *

[0099] Finally, calculate the capillary pressure P C (R) of pore 7 in sub-block 3:

[0100]

[0101] In the formula, f is the geometric roughness coefficient of the pore wall of sub-block 3. f is the ratio of the actual area of the rough surface to its apparent area. σ is the surface tension of mercury, θ is the actual contact angle of mercury on the rough surface, and θ * is the ideal contact angle of mercury;

[0102] Assuming that the surface tension and contact angle of mercury are constant during the mercury injection process, there is a negative correlation between the capillary pressure and the pore radius. When the mercury injection pressure is low, only a small capillary pressure can be overcome to press mercury into pores with a large radius. When the mercury injection pressure increases, a larger capillary pressure can be overcome to press mercury into pores with a small radius;

[0103] R is iterated from large to small to simulate the mercury injection process with the injection pressure increasing from low to high. The larger the radius of pore 7, the smaller the corresponding capillary pressure, and the smaller the radius of pore 7, the larger the corresponding capillary pressure;

[0104] Calculate the mercury injection saturation S(R): The mercury injection saturation S(R) is the ratio of the total mercury injection volume V to the number of voxels in the entire pore 7 space of the three-dimensional binary image 5 of sub-block 3. The total mercury injection volume V is the number of voxels of mercury invading the pore 7 space of sub-block 3, that is, the number of voxels of the pore 7 space in sub-block 3 with a size larger than R and connected to the outside of sub-block 3;

[0105] Figure 9 The following shows the situations of mercury invading the pore 7 space of sub-block 3 at four different sizes of R (the red area in the figure is the mercury injection area): (a) No mercury invades the pore 7 space of sub-block 3; (b) Mercury only invades the pore P1 of sub-block 3. Due to the blockage at T1, P2 is not connected to the outside; (c) Mercury only invades the pores P1, T1, and P2 of sub-block 3. Due to the blockage at T2 and T3, P3 and P4 are not connected to the outside; (d) Mercury invades all the pores 7 of sub-block 3;

[0106] Step 6.6: Let R = R - 1. If R = 0, end the operation; otherwise, go to Step 6.2 and perform mercury injection simulation using a solid spherical structure element with a radius reduced by 1;

[0107] Step 7: First, plot the relationship curves of the capillary pressure P C (R) and the mercury injection saturation S(R) of 8 sub-blocks 3 in a rectangular coordinate system, as Figure 10 shown. Then, evenly take 9 mercury injection saturation values (0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1) on the mercury injection saturation S(R) number axis S. Each mercury injection saturation value is on the capillary pressure P C (R) number axis P CEight capillary pressure values are correspondingly obtained above; Arbitrarily select two of the eight relationship curves, calculate the differences between the corresponding two capillary pressure values at each mercury intrusion saturation value to obtain nine differences, and take the average of these nine differences;

[0108] Step 8. When the average value obtained by arbitrarily selecting any two relationship curves is less than 0.001, it is considered that the relationship curve between the capillary pressure P C (R) and the mercury intrusion saturation S(R) converges. As Figure 11 shown, the size n×n×n of sub-block 3 at this time is the size of the representative volume element REV of the rock; Otherwise, make the side length n of sub-block 3 equal to n+1, then go to step 6, and repeat the mercury intrusion simulation operation on a larger sub-block 3. For example, continue to repeat the mercury intrusion simulation operation on 201×201×201 until the relationship curve between the capillary pressure P C (R) and the mercury intrusion saturation S(R) converges together.

[0109] In step 6.1, the method for determining the initial value of parameter R is as follows:

[0110] Step 6.1.1. Set a parameter R′, where R′ is a non-zero natural number and R′ is greater than the size of the largest pore 7 in the three-dimensional binary image 5 of sub-block 3, and R′ = 50;

[0111] Step 6.1.2. Use a solid spherical structural element with a radius of R′ to perform erosion operation on the three-dimensional binary image 5 of sub-block 3. Perform a logical AND operation on all solid spherical structural elements and the three-dimensional binary image 5 of sub-block 3, and convert the void voxels in the three-dimensional binary image 5 of sub-block 3 whose distance from the solid spherical structural element is less than or equal to R′ into solid voxels;

[0112] Step 6.1.3. If all void voxels have become solid voxels, it means that R′ is too large. Let R′ = R′-1 and execute step 6.1.2; Otherwise, determine the initial value of parameter R as R = R′+1, and then end.

[0113] In step 6.5, as Figure 8 shown, the method for determining the total mercury intrusion volume V is as follows:

[0114] Step 6.5.1. First, define the boundary void voxels a and the internal void voxels of the three-dimensional binary image 5 after the dilation operation in step 6.4, and search for the boundary void voxels a;

[0115] Step 6.5.2: Starting from each boundary void voxel a in sequence, find all the internally connected void voxels b that are connected to each boundary void voxel a; in the three-dimensional binary image 5 that has undergone the dilation operation in Step 6.4, the total number of voxels of all boundary void voxels a and all internally connected void voxels b is the total mercury intrusion volume V.

[0116] The method for finding all the internally connected void voxels b that are connected to a boundary void voxel a is as follows:

[0117] Step 6.5.2.1: First, search for the internal void voxels that are connected to the upper, lower, left, right, front, and rear six faces of this boundary void voxel a; if none are found, end the operation, otherwise, mark the found internal void voxels as internally connected void voxels b, and use these found internally connected void voxels b as the starting voxels for the next search, and then execute Step 6.5.2.2;

[0118] Step 6.5.2.2: Search for the internal void voxels that are connected to the upper, lower, left, right, front, and rear six faces of the starting voxel; if none are found, end the operation, otherwise, mark the found internal void voxels as internally connected void voxels b, and use these found internally connected void voxels b as the starting voxels for the next search, and then execute Step 6.5.2.2.

Claims

1. A method for obtaining a representative volume element of a rock, characterized in that It includes the following steps: Step 1: Scan the rock to obtain a three-dimensional gray-scale image of the rock; Step 2: Convert the three-dimensional gray-scale image into a three-dimensional binary image. In the three-dimensional binary image, void voxels are labeled as 1 and matrix voxels are labeled as 0; Step 3: Take eight cubic sub-blocks from the eight vertices of the three-dimensional binary image of the rock towards the inside of the three-dimensional binary image. The side lengths of the eight sub-blocks are all n, and let n = 1 voxel; Step 4: Calculate the porosity of each sub-block. The calculation method of porosity is as follows: where K i is the porosity of the i-th sub-block, and i is a natural number from 1 to 8; Step 5: When the difference in porosity between any two of the eight sub-blocks is less than the set difference, execute Step 6; otherwise, increase the side length of the sub-block by 1 voxel, that is, n = n + 1, and execute Step 4; Step 6: Perform mercury intrusion simulation operations on the eight sub-blocks respectively. The mercury intrusion simulation operation method for each sub-block is as follows: Step 6.1: Set a parameter R. The initial value of R is a non-zero natural number; Step 6.2: Use a solid spherical structural element with a radius of R to erode the three-dimensional binary image of the sub-block, and convert the void voxels in the three-dimensional binary image of the sub-block whose distance from the solid spherical structural element is less than or equal to R into solid voxels; Step 6.3: If all void voxels have become solid voxels, go to Step 6.6; otherwise, go to Step 6.4; Step 6.4: Use a void spherical structural element with a radius of R to dilate the three-dimensional binary image eroded in Step 6.2, and convert the solid voxels in the three-dimensional binary image eroded in Step 6.2 whose distance from the void spherical structural element is less than or equal to R into void voxels; After the erosion operation and dilation operation, pores with a size larger than R in the sub-block are found; Step 6.5, calculate the capillary pressure P of the pores in the sub-block C (R): where f is the geometric roughness coefficient of the pore wall of the sub-block, σ is the surface tension of mercury, and θ * is the ideal contact angle of mercury; Calculate the mercury intrusion saturation S(R): The mercury intrusion saturation S(R) is the ratio of the total mercury intrusion volume V to the number of voxels in the entire pore space of the three-dimensional binary image of the sub-block. The total mercury intrusion volume V is the number of voxels of mercury invading the pore space of the sub-block, that is, the number of voxels of the pore space in the sub-block with a size larger than R and connected to the outside of the sub-block; Step 6.6: Let R = R - 1. If R = 0, end the operation; otherwise, go to Step 6.2; Step 7: First, plot the relationship curves of capillary pressure P C (R) and mercury injection saturation S(R) for 8 sub-blocks in a rectangular coordinate system respectively. Then, evenly select M mercury injection saturation values on the mercury injection saturation S(R) number axis, where M is a natural number greater than or equal to 2. Each mercury injection saturation value corresponds to 8 capillary pressure values on the capillary pressure P C (R) number axis. Arbitrarily select two relationship curves from the 8 relationship curves, calculate the differences between the corresponding 2 capillary pressure values at each mercury injection saturation value to obtain M differences, and take the average of these M differences; Step 8: When the average value obtained from any two selected relationship curves is less than the set average value, the size n×n×n of the sub-block at this time is the size of the representative volume element of the rock; Otherwise, let the side length n of the sub-block be n = n + 1, and then go to Step 6.

2. The method for obtaining a representative volume element of rock according to claim 1, wherein: In Step 1, a CT scan is used to obtain the three-dimensional gray-scale image of the rock.

3. The method for obtaining a representative volume element of a rock according to claim 1, characterized in that: In Step 5, the set difference is 0.

001.

4. The method for obtaining a representative volume element of rock according to claim 1, characterized in that: In Step 6.1, the method for determining the initial value of the parameter R is: Step 6.1.1: Set a parameter R'. R' is a non-zero natural number and R' is greater than the size of the largest pore in the three-dimensional binary image of the sub-block; Step 6.1.2: Use a solid spherical structural element with a radius of R' to erode the three-dimensional binary image of the sub-block, and convert the void voxels in the three-dimensional binary image of the sub-block whose distance from the solid spherical structural element is less than or equal to R' into solid voxels; Step 6.1.3: If all the void voxels have become solid voxels, let R′ = R′ - 1, and execute Step 6.1.2; otherwise, determine the initial value of parameter R as R = R′ + 1, and then end.

5. The method for obtaining a representative volume element of rock according to claim 1, characterized in that: In Step 6.5 mentioned above, the method for determining the total mercury intrusion volume V is as follows: Step 6.5.1: First, define the boundary void voxels and the internal void voxels of the three-dimensional binary image after the dilation operation in Step 6.4, and search for the boundary void voxels. Step 6.5.2: Starting from each boundary void voxel in turn, find all the internally connected void voxels connected to each boundary void voxel; in the three-dimensional binary image after the dilation operation in Step 6.4, the total number of voxels of all the boundary void voxels and all the internally connected void voxels is the total mercury intrusion volume V. The method for finding all the internally connected void voxels connected to a boundary void voxel is as follows: Step 6.5.2.1: First, search for the internal void voxels connected to the six faces of the upper, lower, left, right, front, and back of this boundary void voxel; if none are found, end the operation, otherwise, mark the found internal void voxels as internally connected void voxels, and use these found internally connected void voxels as the starting voxels for the next search, and then execute Step 6.5.2.

2. Step 6.5.2.2: Search for the internal void voxels connected to the six faces of the upper, lower, left, right, front, and back of the starting voxel; if none are found, end the operation, otherwise, mark the found internal void voxels as internally connected void voxels, and use these found internally connected void voxels as the starting voxels for the next search, and then execute Step 6.5.2.

2.

6. The method for obtaining a representative volume element of rock according to claim 1, characterized in that: In Step 8 mentioned above, set the average value to 0.001.

Images