A SPH method and device for calculating the permeability of porous matrix digital core

Abstract

The present invention belongs to the field of mining and geotechnical engineering technology, and relates to an SPH method and device for calculating the permeability of a porous matrix digital core. The method comprises: obtaining a three-dimensional grayscale data body of a core by micrometer CT scanning, dividing the voxels of the three-dimensional grayscale data body into three types of media: pores, matrix, and skeleton by setting a threshold and dividing the voxels of the three-dimensional grayscale data body into three types of media: pores, matrix, and skeleton according to the grayscale value from low to high; calculating the matrix porosity according to the grayscale value of the matrix voxel, and calculating the sub-resolution permeability of the matrix voxel according to the Kozeny-Carman formula; under the inlet and outlet pressure boundary conditions and the no-slip boundary conditions of the skeleton wall, using SPH to solve the Navier-Stokes equation at the pore and the Brinkman equation at the matrix for a fluid of a given viscosity, and obtaining the velocity field of the digital core; calculating the permeability of the digital core according to Darcy's law. The method provided by the present invention can calculate the overall velocity field of a porous matrix digital core, and expands the scope of application of the calculation of the permeability of the digital core.

Description

Technical Field

[0001] The invention belongs to the technical field of mining and geotechnical engineering, and relates to an SPH method and device for calculating porous matrix type digital core permeability. Background Art

[0002] The three-dimensional pore skeleton data of the core obtained through micron CT scanning and image recognition is called digital core. The porous media seepage numerical simulation technology based on digital core can quickly provide important core physical parameters such as porosity, permeability, saturation, etc., providing a key reference for oil and gas exploration and development decisions.

[0003] However, digital cores based on pore skeleton binary data bodies are almost impossible to use to calculate the permeability of common low-porosity and low-permeability cores. This is because the current scanning range is the entire core, and the scanning range (i.e., field of view) is determined by the resolution. The higher the resolution, the smaller the field of view. To scan small holes, high resolution is required, and the field of view is also small, and the three-dimensional data of the entire core cannot be obtained. To scan the entire core, the scanning resolution must be sacrificed, and the price is that the small holes cannot be seen. Therefore, under the current scanning resolution of micron CT, only incompletely connected large pores can be obtained, and the unrecognized tiny pores are classified as rock skeletons. In the case of non-connected pores, the permeability is zero.

[0004] Therefore, a new digital core model is needed for low-porosity and low-permeability cores, which can include the matrix area where the pores and skeletons cannot be distinguished under the current scanning resolution, so that the pores of the digital core are connected as a whole in the sub-resolution sense; and a numerical method for calculating the permeability of pore-matrix digital cores is needed, which can quickly simulate the seepage problem of this cross-scale porous medium. Summary of the invention

[0005] In view of the above technical problems, the present invention proposes an SPH method and device for calculating the permeability of porous matrix digital cores, so as to achieve the purpose of establishing a porous matrix digital core model that maintains connectivity and quickly simulate cross-scale porous media seepage, thereby solving the permeability calculation problem of low-porosity and low-permeability digital cores.

[0006] The technical solution adopted by the present invention is:

[0007] A SPH (Smoothed Particle Hydrodynamics) method for calculating porous matrix digital core permeability, the method comprising:

[0008] (1) A three-dimensional grayscale data volume of the core is obtained by micrometer CT scanning, and voxels of the three-dimensional grayscale data volume are divided into three types of media: pores, matrix, and skeleton by setting a threshold and following the grayscale values ​​from low to high;

[0009] (2) Calculate the matrix porosity according to the gray value of the matrix voxel, and calculate the sub-resolution permeability of the matrix voxel according to the Kozeny-Carman formula;

[0010] (3) Under the inlet and outlet pressure boundary conditions and the no-slip boundary conditions on the skeleton wall, the SPH is used to solve the Navier-Stokes equation at the pores and the Brinkman equation at the matrix for a fluid with a given viscosity to obtain the velocity field of the digital core;

[0011] (4) Calculate the permeability of the digital core according to Darcy's law based on the velocity field, inlet and outlet pressure difference and fluid viscosity of the digital core.

[0012] Furthermore, step (1) includes:

[0013] The core was scanned by micrometer CT at a scanning resolution of D, and the three-dimensional grayscale data volume of the core was obtained by data reconstruction. The multi-level Otsu algorithm was used to arrange the three-dimensional grayscale data volume from low to high according to the grayscale value of the voxels, and the threshold was set to divide the voxels of the three-dimensional grayscale data volume into three types of media: pores, matrix, and skeleton; the grayscale threshold between pores and matrix was recorded as , the gray threshold between matrix and skeleton is recorded as ;

[0014] The three-dimensional grayscale data volume is transformed, and the grayscale values ​​of the pore voxels are set to , the grayscale values ​​of the skeleton voxels are set to ; The gray value of the matrix voxel remains unchanged. and between.

[0015] Further,

[0016] Step (2) is as follows:

[0017] Constructing gray value g to porosity Functional relationship , so that the porosity of pore voxels is 1, the porosity of skeleton voxels is 0, and the porosity of matrix voxels is between the two. Preset parameters Freely adjustable, when When the porosity Gray value It is a linear relationship. For most real cores, the matrix porosity is very low and its properties are close to those of the skeleton. The value is above 4. The porosity of pore voxels is 1, the porosity of skeleton voxels is 0, and the porosity of matrix voxels is between 0 and 1.

[0018] The pores in the matrix are not visible at the current scanning resolution D of micron CT, so they are called sub-resolution pores, and the corresponding permeability is called sub-resolution permeability, denoted by k.

[0019] The sub-resolution permeability k within the matrix is ​​calculated using the Kozeny-Carman formula:

[0020] ;

[0021] Where d is the effective sub-resolution particle size of the matrix, and d is taken as the current scanning resolution, that is, d=D; is the matrix porosity; .

[0022] Furthermore, in step (3), the inlet and outlet pressure boundary conditions include: setting a buffer zone with a thickness of a smooth length at the inlet and outlet;

[0023] For the inlet, the fluid particle pressure in the inlet buffer is set to the inlet pressure, and the velocity component of the fluid particles in the inlet buffer parallel to the inlet cross section is set to 0, so that the fluid particles in the inlet buffer move in a direction perpendicular to the inlet cross section, and only the velocity is updated after the fluid particles move, but the pressure is not updated. When a fluid particle in the inlet buffer leaves the inlet buffer and enters the flow field, the fluid particle that leaves is regarded as a regular fluid particle. If the fluid particle leaves the buffer for the first time, a new fluid particle is added to the inlet buffer at a position that is a smooth length behind the current position of the fluid particle that leaves the buffer for the first time along the direction of movement, and the velocity and pressure of the fluid particle that leaves the buffer for the first time are assigned to the new fluid particle.

[0024] For the outlet, the fluid particle pressure in the outlet buffer is set as the outlet pressure, the movement direction of the fluid particles in the outlet buffer is limited to the direction perpendicular to the outlet cross section, and after the fluid particles move, only the velocity is updated without updating the pressure. When a fluid particle in the outlet buffer leaves the buffer and moves away from the flow field, the corresponding fluid particle is deleted;

[0025] The skeleton wall no-slip boundary condition is specifically as follows: the relative velocity between the skeleton wall particles and the nearby fluid particles is calculated at each time step; the shortest distance from each skeleton wall particle to the fluid particle is searched once every 10 time steps;

[0026] Will be located at The distance from the i-th skeleton wall particle to the current nearest fluid particle is recorded as , then the current i-th skeleton wall particle and the nearby The position speed is The relative velocity of the fluid particles is:

[0027] ;

[0028] In the formula is the scanning resolution of micrometer CT.

[0029] Furthermore, step (3) is specifically as follows:

[0030] The pore voxels, matrix voxels and skeleton surface voxels are regarded as fluid particles at the pores, fluid particles at the matrix and solid particles on the wall, respectively, that is, the particle positions of the SPH system are initialized, and the initial velocities of all particles are set to 0;

[0031] The momentum equation of fluid particles includes two cases: the pore region and the matrix region:

[0032] Moment, pore area The Navier-Stokes equation for the particle at is:

[0033] ;

[0034] Moment, matrix area The Brinkman equation for the particle at is:

[0035] ;

[0036] ;

[0037] ;

[0038] In the formula, is the resultant acceleration of particles in the pore region; is the resultant acceleration of the particles in the matrix region, where is the pressure acceleration, is the viscous acceleration, is the Brinkman drag acceleration;

[0039] , They are The position of the i-th and j-th fluid particles in the pore region at time instant; , They are The position of the i-th and j-th fluid particles in the matrix region at time instant;

[0040] , , , The meanings are The velocity, pressure, number density and mass of the i-th fluid particle at the moment; , , The meanings are The velocity, pressure, and number density of the jth fluid particle at time instant; is the fluid viscosity; is the matrix porosity; k is the sub-resolution permeability at the matrix; subscript j represents the other particles except the i-th fluid particle within the support set of the smooth kernel function of the i-th fluid particle; is the gradient of the smooth kernel function between the ith fluid particle and the jth particle with respect to the position of the ith fluid particle;

[0041] The conventional time-marching scheme is used for the Navier-Stokes equations in the pore region, and the Strang splitting scheme is used for the Brinkman equations in the matrix region.

[0042] Furthermore, due to the large coefficient of the Brinkman drag term in the matrix region, the momentum equation is extremely rigid. If the conventional SPH time-marching format is used to update the flow field, a very small time step must be used due to the influence of the Courant-Friedrichs-Lewy stability condition. To avoid a time step that is too small, the conventional time-marching format is used only in the pore region for the Navier-Stokes equation, and the Strang splitting format is used for the time marching of the Brinkman equation in the matrix region.

[0043] Furthermore, in step (3), the Navier-Stokes equation in the pore region adopts a conventional time-marching format, specifically:

[0044] remember The time step length is , the time-marching format of the Navier-Stokes equation in the pore region is:

[0045] ;

[0046] ;

[0047] ;

[0048] in, , The meanings are respectively the predicted position of the first 1 / 2 time step forward and the predicted position according to and Calculated acceleration;

[0049] yes The velocity of the i-th fluid particle at the moment; yes The position of the i-th fluid particle in the pore region at the moment;

[0050] In the matrix region, the time advancement of the Brinkman equation uses the Strang splitting format:

[0051] ;

[0052] ;

[0053] ;

[0054] ;

[0055] ;

[0056] ;

[0057] in, , , , , The meanings are respectively the predicted position of the first half time step forward, the first predicted velocity obtained by advancing the pressure acceleration by one time step, the second predicted velocity obtained by advancing the Brinkman resistance acceleration by the first half time step after entering the Strang splitting format, the third predicted velocity obtained by advancing the viscosity acceleration by one time step, and the third predicted velocity obtained by advancing the Brinkman resistance acceleration by the last half time step after completing the Strang splitting format. The particle speed at the moment; yes Moment i The position of each fluid particle in the matrix region;

[0058] Step size of the time marching format The Courant-Friedrichs-Lewy stability condition is satisfied, which is independent of the Brinkman resistance term in the matrix region:

[0059] ;

[0060] ;

[0061] ;

[0062] In the formula, h, , , , The meanings are smooth length, speed of sound, mass of the i-th fluid particle, The pressure acceleration and number density of the i-th fluid particle at the moment; is the fluid viscosity; is the porosity; k is the sub-resolution permeability at the matrix.

[0063] Furthermore, step (4) is specifically as follows:

[0064] Calculate Darcy velocity of digital core based on velocity field of digital core , that is, calculate the ratio of the sum of the velocity components of the velocity field of fluid particles in the pore area and the matrix area along the overall seepage direction to the total number of fluid particles;

[0065] According to Darcy's law, the overall permeability of the digital core is:

[0066] ;

[0067] In the formula, is the fluid viscosity, L is the length of the digital core along the overall seepage direction, , are the inlet pressure and outlet pressure respectively.

[0068] A SPH device for calculating porous matrix type digital core permeability, the device comprising:

[0069] Digital core 3D data acquisition module: Use micron CT scanning to obtain the 3D grayscale data of the core, and divide the 3D data into three types of media: pores, matrix, and skeleton by threshold division according to the grayscale value from low to high;

[0070] Matrix sub-resolution porosity and permeability configuration module: used to set the sub-resolution porosity of matrix voxels according to the gray value of the matrix voxels, and to set the sub-resolution permeability of the matrix voxels according to the Kozeny-Carman formula;

[0071] SPH initialization module: used to convert pore voxels, matrix voxels and skeleton surface voxels into fluid particles at pores, fluid particles at matrix and solid particles on the wall according to their positions, and set the initial velocity of all particles to 0;

[0072] SPH simulation and digital core overall permeability calculation module: used to perform pressure-driven SPH numerical simulation and count the Darcy velocity of the flow field to calculate the overall permeability of the digital core.

[0073] An electronic device, comprising:

[0074] a memory storing computer-readable instructions;

[0075] The processor reads the computer-readable instructions stored in the memory to execute the SPH method for calculating the permeability of porous matrix digital core.

[0076] A computer-readable storage medium stores computer-readable instructions. When the computer-readable instructions are executed by a processor of a computer, the computer is caused to execute the SPH method for calculating the permeability of a porous matrix digital core.

[0077] The beneficial technical effects of the present invention are:

[0078] The SPH method for calculating the permeability of porous matrix digital core provided by the present invention establishes a porous matrix digital core with connectivity maintained in a sub-resolution sense, and then uniformly solves the flow fields at the pores and the matrix, so that even if the sub-resolution micropores at the matrix are invisible, the overall velocity field of the porous matrix digital core can still be calculated, thereby expanding the scope of application of the digital core permeability calculation.

[0079] In addition, if the conventional time-advancing format is used for calculation, the equation will be too rigid, resulting in a time that is too short and inefficient. To this end, the present invention uses the Strang splitting technique to decouple the Brinkman resistance term from other terms in the momentum equation, and analytically processes the resistance term, so that the SPH method advances time at a conventional time step, thereby ensuring high efficiency of flow field simulation. BRIEF DESCRIPTION OF THE DRAWINGS

[0080] Figure 1 A flow chart of an SPH method for calculating porous matrix digital core permeability provided by one embodiment of the present invention.

[0081] Figure 2 A structural block diagram of an SPH device for calculating porous matrix digital core permeability provided by one embodiment of the present invention. DETAILED DESCRIPTION

[0082] Example embodiments will now be described more fully with reference to the accompanying drawings. However, example embodiments can be implemented in a variety of forms and should not be construed as limited to the examples set forth herein; rather, these example embodiments are provided so that the description of the present invention will be more comprehensive and complete and the concepts of the example embodiments will be fully conveyed to those skilled in the art. The accompanying drawings are only schematic illustrations of the present invention and are not necessarily drawn to scale. The same reference numerals in the figures represent the same or similar parts, and thus their repeated description will be omitted.

[0083] In addition, the described features, structures or characteristics may be combined in one or more example embodiments in any suitable manner. In the following description, many specific details are provided to provide a full understanding of the example embodiments of the present invention. However, those skilled in the art will appreciate that the technical solution of the present invention may be practiced while omitting one or more of the specific details, or other methods, components, steps, etc. may be adopted. In other cases, known structures, methods, implementations or operations are not shown or described in detail to avoid obscuring various aspects of the present invention.

[0084] Some of the blocks shown in the accompanying drawings are functional entities that do not necessarily correspond to physically or logically independent entities. These functional entities can be implemented in software form, or in one or more hardware modules or integrated circuits, or in different networks and / or processor devices and / or microcontroller devices.

[0085] Example 1: An SPH method for calculating the permeability of a porous matrix digital core, such as Figure 1 As shown, the method includes:

[0086] (1) Using micrometer CT scanning to obtain a three-dimensional grayscale data volume of the core, and dividing the voxels of the three-dimensional grayscale data volume into three types of media: pores, matrix, and skeleton by setting a threshold;

[0087] (2) The matrix porosity is set according to the gray value of the matrix voxel, and the sub-resolution permeability of the matrix voxel is set according to the Kozeny-Carman formula;

[0088] (3) Under the inlet and outlet pressure boundary conditions and the no-slip boundary conditions of the skeleton wall, the velocity field of the digital core is obtained by using SPH to solve the Navier-Stokes equation at the pores and the Brinkman equation at the matrix for a fluid with a given viscosity;

[0089] (4) Based on the velocity field, inlet and outlet pressure difference and fluid viscosity of the digital core, the permeability of the digital core is calculated according to Darcy's law.

[0090] In the present invention, step (1) is specifically:

[0091] The multilevel Otsu algorithm is used to arrange the three-dimensional grayscale data volume from low to high according to the grayscale value of the voxels, and the threshold is set to divide the voxels of the three-dimensional grayscale data volume into three types of media: pores, matrix, and skeleton; the grayscale threshold between pores and matrix is ​​recorded as , the gray threshold between matrix and skeleton is recorded as ; Among them, the Multilevel Otsu algorithm is a conventional technology, see the literature: Ping-Sung Liao, Tse-Sheng Chen, Pau-Choo Chung. A Fast Algorithm for Multilevel Thresholding[J]. Journal of Information Science and Engineering , 2001, 17(5):713-727. DOI:10.1016 / S1352-0237(02)00379-9.

[0092] The three-dimensional grayscale data volume is transformed, and the grayscale values ​​of the pore voxels are set to , the grayscale values ​​of the skeleton voxels are set to , the grayscale value of the matrix voxel remains unchanged (the grayscale value of the matrix voxel is Hewei between).

[0093] In the present invention, step (2) is specifically:

[0094] Constructing gray value g to porosity Functional relationship :

[0095] ;

[0096] in, is the gray value of the pore voxel; is the gray value of the skeleton voxel; preset parameters , Used to adjust the mapping intensity. Matrix porosity Voxel gray value It is a linear relationship. The larger the value, the smaller the matrix porosity, and the matrix properties are closer to the skeleton. The value of is adjusted according to the actual situation of the core and the experience of the algorithm user. For most real cores, the matrix porosity is very low and its properties are close to those of the skeleton, so the value is 4 or above; the porosity of the pore voxel is 1, the porosity of the skeleton voxel is 0, and the porosity of the matrix voxel is between the two;

[0097] The pores in the matrix are scanned at the current resolution of micrometer CT. It is invisible under the microscope, so it is called sub-resolution pores, and the corresponding permeability is called sub-resolution permeability, denoted by k.

[0098] The sub-resolution permeability k corresponding to the sub-resolution pores in the matrix is ​​calculated using the Kozeny-Carman formula:

[0099] ;

[0100] Where d is the effective sub-resolution particle size of the matrix, and d is taken as the current scanning resolution, that is, d=D; is the matrix porosity; .

[0101] In this embodiment, step (3) includes initialization of the SPH system, the momentum equation form solved by the SPH method, boundary conditions, and the time advancement format of the solution process;

[0102] In step (3) of the present invention, the inlet and outlet pressure boundary conditions include: setting a buffer zone with a thickness of a smooth length at both the inlet and the outlet;

[0103] For the inlet, the fluid particle pressure in the inlet buffer is set to the inlet pressure, and the velocity component of the fluid particles in the inlet buffer parallel to the inlet cross section is set to 0, so that the fluid particles in the inlet buffer move in a direction perpendicular to the inlet cross section, and only the velocity is updated after the fluid particles move, but not the pressure. When a fluid particle in the inlet buffer leaves the buffer and enters the flow field, the corresponding fluid particle is regarded as a regular fluid particle. If the fluid particle leaves the buffer for the first time, a new fluid particle is added to the inlet buffer at a position that is a smooth length behind the current position of the fluid particle that leaves the buffer for the first time along the direction of movement, and the velocity and pressure of the fluid particle that leaves the buffer for the first time are assigned to the new fluid particle.

[0104] For the outlet, the fluid particle pressure in the outlet buffer is set as the outlet pressure. Similar to the fluid particles in the inlet buffer, the movement direction of the fluid particles in the outlet buffer is restricted to the direction perpendicular to the outlet cross section, and after the movement, only the speed is updated without updating the pressure. When a fluid particle in the outlet buffer leaves the buffer and moves away from the flow field, the corresponding fluid particle is deleted.

[0105] The no-slip boundary condition for the skeleton wall is specifically:

[0106] The relative speed between the skeleton wall particles and the nearby fluid particles is calculated at each time step; the shortest distance from each skeleton wall particle to the fluid particle is searched once every 10 time steps;

[0107] Will be located at The distance from the i-th skeleton wall particle to the current nearest fluid particle is recorded as , then the current i-th skeleton wall particle and the nearby The position speed is The relative velocity of the fluid particles is:

[0108] ;

[0109] In the formula is the scanning resolution of micrometer CT.

[0110] The SPH system initialization in step (3) is specifically as follows:

[0111] The pore voxels, matrix voxels and skeleton surface voxels are regarded as fluid particles at the pores, fluid particles at the matrix and solid particles on the wall, respectively, that is, the particle positions of the SPH system are initialized, and the initial velocities of all particles are set to 0;

[0112] The momentum equation solved by the SPH method in step (3) is specifically in the form of:

[0113] The momentum equation of fluid particles includes two cases: the pore region and the matrix region:

[0114] Moment, pore area The Navier-Stokes equation for the particle at is:

[0115] ;

[0116] Moment, matrix area The Brinkman equation for the particle at is:

[0117] ;

[0118] ;

[0119] ;

[0120] In the formula, is the resultant acceleration of particles in the pore region; is the resultant acceleration of the particles in the matrix region, where is the pressure acceleration, is the viscous acceleration, is the Brinkman drag acceleration;

[0121] , They are The position of the i-th and j-th fluid particles in the pore region at time instant; , They are The position of the i-th and j-th fluid particles in the matrix region at time instant;

[0122] , , , The meanings are The velocity, pressure, number density and mass of the i-th fluid particle at the moment; , , The meanings are The velocity, pressure, and number density of the jth fluid particle at time instant; is the fluid viscosity; is the matrix porosity; k is the sub-resolution permeability at the matrix; subscript j represents the other particles except the i-th fluid particle within the support set of the smooth kernel function of the i-th fluid particle; is the gradient of the smooth kernel function between the ith fluid particle and the jth particle with respect to the position of the ith fluid particle;

[0123] Since the coefficient of the Brinkman resistance term in the matrix region is too large, the momentum equation is extremely rigid. If the conventional SPH time-marching format is used to update the flow field, a very small time step must be used due to the influence of the Courant-Friedrichs-Lewy stability condition. To avoid too small a time step, the conventional time-marching format is used only in the pore region for the Navier-Stokes equation, and the Strang splitting format is used for the time marching of the Brinkman equation in the matrix region.

[0124] In step (3) of the present invention, the Navier-Stokes equation in the pore region adopts a conventional time-marching format, specifically:

[0125] remember The time step length is , the time-marching format of the Navier-Stokes equation in the pore region is:

[0126] ;

[0127] ;

[0128] ;

[0129] in, , The meanings are respectively the predicted position of the first 1 / 2 time step forward and the predicted position according to and Calculated acceleration;

[0130] yes The velocity of the i-th fluid particle at the moment; yes The position of the i-th fluid particle in the pore region at the moment;

[0131] Since the coefficient of the Brinkman drag term in the matrix region is too large, the momentum equation is extremely rigid. If the conventional SPH time advancement format is used to update the flow field, a very small time step must be used due to the influence of the Courant-Friedrichs-Lewy stability condition. To avoid a time step that is too small, the Strang splitting format is used for the time advancement of the Brinkman equation in the matrix region. The time advancement format of the Brinkman equation in the matrix region is:

[0132] ;

[0133] ;

[0134] ;

[0135] ;

[0136] ;

[0137] ;

[0138] in, , , , , The meanings are respectively the predicted position of the first half time step forward, the first predicted velocity obtained by advancing the pressure acceleration by one time step, the second predicted velocity obtained by advancing the Brinkman resistance acceleration by the first half time step after entering the Strang splitting format, the third predicted velocity obtained by advancing the viscosity acceleration by one time step, and the third predicted velocity obtained by advancing the Brinkman resistance acceleration by the last half time step after completing the Strang splitting format. The particle speed at the moment; yes Moment i The position of each fluid particle in the matrix region;

[0139] The step size of the above time marching format The only requirement is the Courant-Friedrichs-Lewy stability condition, which is independent of the Brinkman resistance term in the matrix region:

[0140] ;

[0141] ;

[0142] ;

[0143] In the formula, h, , , , The meanings are smooth length, speed of sound, mass of the i-th fluid particle, The pressure acceleration and number density of the i-th fluid particle at the moment; is the fluid viscosity; is the porosity; k is the sub-resolution permeability at the matrix.

[0144] Step (4) of the present invention is specifically as follows:

[0145] Statistical Darcy velocity of digital core ; That is, the ratio of the sum of the velocity components of the average velocity field of fluid particles in the pore area and the matrix area along the overall seepage direction to the total number of fluid particles;

[0146] According to Darcy's law, the overall permeability of the digital core is:

[0147] ;

[0148] In the formula, is the fluid viscosity, L is the length of the digital core along the overall seepage direction, , are the inlet pressure and outlet pressure respectively.

[0149] The SPH method for calculating the permeability of porous matrix digital core provided by the present invention first establishes a porous matrix digital core model that maintains connectivity in a sub-resolution sense, and then uses a set of numerical methods to uniformly solve the flow field at the pores and the matrix, so that even if the sub-resolution micropores at the matrix are invisible, the overall velocity field of the porous matrix digital core can still be calculated, expanding the scope of application of the digital core permeability calculation. In addition, if this numerical method is calculated using a conventional time-marching format, the equation will be too rigid, resulting in a time that is not too long or too small, and the efficiency is extremely low. To this end, the present invention uses the Strang splitting technique to decouple the Brinkman resistance term from other terms of the momentum equation, and analytically processes the resistance term, so that the SPH method advances time with a conventional time step, thereby ensuring the high efficiency of flow field simulation.

[0150] The method provided by the present invention uniformly solves the flow field at the pores and the matrix, so that even if the sub-resolution micropores at the matrix are invisible, the overall velocity field of the pore matrix type digital core can still be calculated, which expands the scope of application of the digital core permeability calculation. When solving the Brinkman equation at the matrix, the sub-resolution permeability of the matrix voxel is too low, resulting in a huge coefficient of the Brinkman resistance term, causing the rigidity of the SPH momentum equation to be greatly enhanced. In order to avoid an extremely small time step, the present invention uses the Strang splitting technique to decouple the Brinkman resistance term from other terms of the momentum equation, and analytically processes the resistance term, ensuring that the SPH method performs flow field simulation with a conventional time step.

[0151] Example 2: An SPH device for calculating the permeability of a porous matrix digital core, such as Figure 2 As shown, the device comprises:

[0152] Digital core 3D data acquisition module 310: obtains 3D grayscale data of the core using micrometer CT scanning, and divides the 3D data into three types of media: pores, matrix, and skeleton by threshold division according to grayscale values ​​from low to high;

[0153] Matrix region sub-resolution porosity and permeability configuration module 320: used to set the sub-resolution porosity of the matrix voxel according to the gray value of the matrix voxel, and to set the sub-resolution permeability of the matrix voxel according to the Kozeny-Carman formula;

[0154] SPH initialization module 330: used to convert pore voxels, matrix voxels and skeleton surface voxels into fluid particles at pores, fluid particles at matrix and solid particles on the wall according to their positions, and set the initial velocity of all particles to 0;

[0155] The SPH simulation and digital core overall permeability calculation module 340 is used to perform pressure-driven SPH numerical simulation and to calculate the Darcy velocity of the flow field to calculate the overall permeability of the digital core.

[0156] In this embodiment, the digital core three-dimensional data volume acquisition module 310 and the matrix sub-resolution porosity and permeability configuration module 320 establish a pore matrix digital core model that maintains connectivity in the sub-resolution sense, and provide sub-resolution porosity and permeability that vary with position in the matrix area, so that the pore area and matrix area of ​​the digital core are respectively adapted to the Navier-Stoke equation and the Brinkman equation. Then, the velocity field of the digital core is simulated under pressure drive through the SPH initialization module 330 and the SPH simulation and digital core overall permeability calculation module 340, and finally the overall permeability of the core is obtained by Darcy's law statistics.

[0157] When performing SPH flow field numerical simulation, the Strang splitting technique is used in the matrix region to decouple the Brinkman drag term from other terms of the momentum equation, and the drag term is analytically processed. This ensures that although the momentum equation in the matrix region is extremely rigid, the time marching format can still be calculated using a conventional time step, thus ensuring the high efficiency of the SPH numerical simulation.

[0158] The implementation process of the functions and effects of each module / unit in the above-mentioned device is specifically described in the implementation process of the corresponding steps in the above-mentioned method, which will not be repeated here.

[0159] It is understood that these modules can be implemented by hardware, software, or a combination of the two. When implemented in hardware, these modules can be implemented as one or more hardware modules, such as one or more application specific integrated circuits. When implemented in software, these modules can be implemented as one or more computer programs executed on one or more processors.

[0160] Embodiment 3: An electronic device, comprising:

[0161] a memory storing computer-readable instructions;

[0162] The processor reads the computer-readable instructions stored in the memory to execute the method of simulating external force-driven scene flow by smooth particle dynamics described in Example 1.

[0163] Embodiment 4: A computer-readable storage medium stores computer-readable instructions. When the computer-readable instructions are executed by a processor of a computer, the computer executes the method of simulating external force-driven scene flow by smooth particle dynamics described in Embodiment 1.

[0164] Through the description of the above implementation methods, it is easy for those skilled in the art to understand that the example implementation methods described here can be implemented by software, or by software combined with necessary hardware. Therefore, the technical solution according to the implementation method of the present invention can be embodied in the form of a software product, which can be stored in a non-volatile storage medium (which can be a CD-ROM, a USB flash drive, a mobile hard disk, etc.) or on a network, and includes several instructions to enable a computing device (which can be a personal computer, a server, a mobile terminal, or a network device, etc.) to execute the method according to the implementation method of the present invention.

[0165] Those skilled in the art will readily appreciate other embodiments of the present invention after considering the specification and practicing the invention claimed herein. The present invention is intended to cover any variations, uses or adaptations of the present invention that follow the general principles of the present invention and include common knowledge or customary techniques in the art for which the present invention is not claimed. The specification and examples are to be considered exemplary only, and the true scope and spirit of the present invention are indicated by the appended claims.

Claims

1. A SPH method for calculating the permeability of porous matrix digital core, characterized in that: The method comprises: (1) obtaining a three-dimensional grayscale data volume of the core by micrometer CT scanning, and dividing the voxels of the three-dimensional grayscale data volume into three types of media: pores, matrix, and skeleton by setting a threshold and following the grayscale values ​​from low to high; (2) Calculate matrix porosity based on the grayscale value of matrix voxels, and calculate sub-resolution permeability of matrix voxels based on the Kozeny-Carman formula; (3) Under the inlet and outlet pressure boundary conditions and the no-slip boundary conditions on the skeleton wall, the SPH is used to solve the Navier-Stokes equation at the pores and the Brinkman equation at the matrix for a fluid with a given viscosity to obtain the velocity field of the digital core; In the matrix region, the time advancement of the Brinkman equation uses the Strang splitting format: in, The meanings are respectively the predicted position of the first half time step forward, the first predicted velocity obtained by advancing the pressure acceleration by one time step, the second predicted velocity obtained by advancing the Brinkman resistance acceleration by one half time step into the Strang splitting format, the third predicted velocity obtained by advancing the viscosity acceleration by one time step, and the t obtained by advancing the Brinkman resistance acceleration by one half time step after completing the Strang splitting format. n+1 The particle speed at the moment; Yes n+1 The position of the i-th fluid particle in the matrix region at the moment; (4) Calculate the permeability of the digital core according to the velocity field, inlet and outlet pressure difference and fluid viscosity of the digital core in accordance with Darcy's law.

2. The SPH method for calculating the permeability of porous matrix digital core according to claim 1, characterized in that: Step (1) is specifically as follows: The multi-level Otsu algorithm is used to arrange the three-dimensional grayscale data volume from low to high according to the grayscale value of the voxels, and the threshold is set to divide the voxels of the three-dimensional grayscale data volume into three types of media: pores, matrix, and skeleton; the grayscale threshold between the pores and the matrix is ​​recorded as g0, and the grayscale threshold between the matrix and the skeleton is recorded as g1; The three-dimensional grayscale data volume is transformed, the grayscale values ​​of the pore voxels are all set to g0; the grayscale values ​​of the skeleton voxels are all set to g1; the grayscale values ​​of the matrix voxels remain unchanged and are between g0 and g1.

3. The SPH method for calculating the permeability of porous matrix digital core according to claim 2, characterized in that: In step (2), the sub-resolution permeability k in the matrix is ​​calculated using the Kozeny-Carman formula: Wherein, d is the sub-resolution effective particle size of the matrix, and the value of d is the current scanning resolution, that is, d=D; ε is the matrix porosity; 0<ε<1.

4. The SPH method for calculating the permeability of porous matrix digital core according to claim 1, characterized in that: In step (3), the inlet and outlet pressure boundary conditions include: setting a buffer zone with a thickness of a smooth length at the inlet and outlet; For the inlet, the fluid particle pressure in the inlet buffer is set to the inlet pressure, and the velocity component of the fluid particles in the inlet buffer parallel to the inlet cross section is set to 0, so that the fluid particles in the inlet buffer move in a direction perpendicular to the inlet cross section, and only the velocity is updated after the fluid particles move, but the pressure is not updated. When a fluid particle in the inlet buffer leaves the inlet buffer and enters the flow field, the fluid particle that leaves is regarded as a regular fluid particle. If the fluid particle leaves the buffer for the first time, a new fluid particle is added to the inlet buffer at a position that is a smooth length behind the current position of the fluid particle that leaves the buffer for the first time along the direction of movement, and the velocity and pressure of the fluid particle that leaves the buffer for the first time are assigned to the new fluid particle. For the outlet, the fluid particle pressure in the outlet buffer is set as the outlet pressure, the movement direction of the fluid particles in the outlet buffer is limited to the direction perpendicular to the outlet cross section, and after the fluid particles move, only the velocity is updated without updating the pressure. When a fluid particle in the outlet buffer leaves the buffer and moves away from the flow field, the corresponding fluid particle is deleted; The skeleton wall no-slip boundary condition is specifically as follows: the relative velocity between the skeleton wall particles and the nearby fluid particles is calculated at each time step; the shortest distance from each skeleton wall particle to the fluid particle is searched once every 10 time steps; will be located at z i The distance from the i-th skeleton wall particle to the current closest fluid particle is recorded as d i , then the current i-th skeleton wall particle and the nearby particle at z j The position velocity is v j The relative velocity of the fluid particles is: Where D is the scanning resolution of micrometer CT.

5. The SPH method for calculating the permeability of porous matrix digital core according to claim 2, characterized in that: Step (3) is specifically as follows: The pore voxels, matrix voxels and skeleton surface voxels are regarded as fluid particles at the pores, fluid particles at the matrix and solid particles on the wall, respectively, that is, the particle positions of the SPH system are initialized, and the initial velocities of all particles are set to 0; The momentum equation of fluid particles includes two cases: the pore region and the matrix region: t n Moment, pore area The Navier-Stokes equation for the particle at is: t n Moment, matrix area The Brinkman equation for the particle at is: In the formula, is the resultant acceleration of particles in the pore region; is the resultant acceleration of the particles in the matrix region, where is the pressure acceleration, is the viscous acceleration, is the Brinkman drag acceleration; They are t n The position of the i-th and j-th fluid particles in the pore region at time instant; They are t n The position of the i-th and j-th fluid particles in the matrix region at time instant; m i The meanings are t n The velocity, pressure, number density and mass of the i-th fluid particle at the moment; The meanings are t n The velocity, pressure, and number density of the jth fluid particle at the moment; μ is the fluid viscosity; ε is the matrix porosity; κ is the sub-resolution permeability at the matrix; the subscript j represents the other particles except the i-th fluid particle within the support set of the smooth kernel function of the i-th fluid particle; is the gradient of the smooth kernel function between the ith fluid particle and the jth particle with respect to the position of the ith fluid particle; The conventional time-marching scheme is used for the Navier-Stokes equations in the pore region, and the Strang splitting scheme is used for the Brinkman equations in the matrix region.

6. The SPH method for calculating the permeability of porous matrix digital core according to claim 5, characterized in that: In step (3), the Navier-Stokes equation in the pore region adopts a conventional time-marching format, specifically: Note n The time step length at time is Δt n , the time-marching format of the Navier-Stokes equation in the pore region is: in, The meanings are respectively the predicted position of the first 1 / 2 time step forward and the predicted position according to and Calculated acceleration; Yes n+1 The velocity of the i-th fluid particle at the moment; Yes n+1 The position of the i-th fluid particle in the pore region at the moment; Step size Δt of the time marching format n The Courant-Friedrichs-Lewy stability condition is satisfied, which is independent of the Brinkman resistance term in the matrix region: In the formula, h, c s 、m i , The meanings are smooth length, speed of sound, mass of the i-th fluid particle, t n The pressure acceleration and number density of the i-th fluid particle at time instant; μ is the fluid viscosity; ε is the porosity; ε is the sub-resolution permeability at the matrix.

7. The SPH method for calculating the permeability of porous matrix digital core according to claim 1, characterized in that: Step (4) is specifically as follows: Calculate Darcy velocity of digital core based on velocity field of digital core That is, the ratio of the sum of the velocity components of the velocity field of fluid particles in the pore area and the matrix area along the overall seepage direction to the total number of fluid particles is calculated; According to Darcy's law, the overall permeability of the digital core is: Where μ is the fluid viscosity, L is the length of the digital core along the overall seepage direction, pin and pout are the inlet pressure and outlet pressure, respectively.

8. An SPH device for calculating porous matrix digital core permeability, using the method according to any one of claims 1 to 7, characterized in that: The device comprises: Digital core 3D data acquisition module: Use micron CT scanning to obtain the 3D grayscale data of the core, and divide the 3D data into three types of media: pores, matrix, and skeleton by threshold division according to the grayscale value from low to high; Matrix area sub-resolution porosity and permeability configuration module: used to set the matrix voxel sub-resolution porosity according to the gray value of the matrix voxel, and to set the matrix voxel sub-resolution permeability according to the Kozeny-Carman formula; SPH initialization module: used to convert pore voxels, matrix voxels and skeleton surface voxels into fluid particles at pores, fluid particles at matrix and solid particles on the wall according to their positions, and set the initial velocity of all particles to 0; SPH simulation and digital core overall permeability calculation module: used to perform pressure-driven SPH numerical simulation and count the Darcy velocity of the flow field to calculate the overall permeability of the digital core.

9. An electronic device, characterized in that: include: a memory storing computer-readable instructions; A processor reads computer-readable instructions stored in a memory to execute the SPH method for calculating the permeability of a porous matrix digital core as described in any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that: Computer-readable instructions are stored, and when the computer-readable instructions are executed by a processor of a computer, the computer is caused to execute the SPH method for calculating the permeability of a porous matrix digital core as described in any one of claims 1 to 7.

Images