Method for determining apparent permeability of carbon dioxide injected into shale gas by considering irreducible water distribution
Through scanning electron microscopy image reconstruction and combining with lattice Boltzmann simulation, the problem of low accuracy in determining the permeability of shale gas injection carbon dioxide is solved, and the accurate permeability calculation is achieved under the influence of bound water in shale. High flexibility and low cost.
Patent Information
- Application Number
- CN202411231031.8
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2024-09-04
- Publication Date
- 2024-12-13
AI Technical Summary
Existing experimental methods are difficult to accurately determine the permeability of shale gas injection carbon dioxide at ultra-low water saturation. The conventional plunger experiment has low accuracy and cannot effectively consider the impact of bound water on shale pore flow.
SEM images were used to obtain shale images, and the adsorption and desorption of bound water were simulated through digital core reconstruction and lattice Boltzmann model. Combined with D3Q19-LBGK and MRT-LBM models, the flow rate of carbon dioxide in aqueous shale was calculated, and finally used Darcy's formula calculates apparent penetration rate.
It realizes the accurate calculation of gas permeability in aqueous shale under different stratigraphic conditions, and has the advantages of high flexibility, good repeatability and low cost.
Smart Images
Figure CN119124962A_ABST
Abstract
Description
Technical Field
[0001] The invention relates to the technical field of microscopic numerical simulation of oil and gas reservoirs, and in particular to a method for determining the apparent permeability of shale gas carbon dioxide injection taking into account the distribution of irreducible water. Background Art
[0002] Shale gas is a green and environmentally friendly unconventional oil and gas resource. my country's shale gas resource development potential is huge, and it is one of the important successor resources for my country's natural gas resources to achieve increased reserves and production in the future. Under normal circumstances, there is a certain amount of primary water in the original shale reservoir, and the action mechanism in the shale pores is significantly different at different water saturations, which in turn affects the flow capacity of the shale pores. Existing experimental methods such as the equilibrium water experiment generally obtain shale cores that are saturated with bound water, while shale has the characteristics of ultra-low water saturation. Therefore, there is an urgent need for a method to obtain shale cores at ultra-low water saturation.
[0003] Shale is rich in micro-nano pores, including both organic and inorganic pores. These characteristics make the storage and transmission of fluids in shale nanopores extremely complex. However, the conventional plunger test method has low instrument accuracy and large errors in the permeability obtained by the test, making it difficult to determine the gas permeability of shale gas after carbon dioxide injection in the presence of bound water. Therefore, it is necessary to carry out relevant numerical simulation research to obtain the gas permeability of shale gas injected with carbon dioxide in the presence of bound water. Summary of the invention
[0004] The present invention mainly overcomes the deficiencies in the prior art and aims to provide a method for visualizing the distribution of original bound water in shale.
[0005] In order to achieve the above technical objectives, the present invention adopts the following technical solutions:
[0006] S1: Obtain shale SEM, porosity and pore size distribution;
[0007] S2: Based on the shale scanning electron microscope image, a shale digital core reconstruction method with porosity as a constraint is established to reconstruct the shale three-dimensional digital core;
[0008] S3: The thirteen-direction average algorithm is used to calculate the pore size distribution, and then the connectivity labeling algorithm is used to determine the connectivity of the pores, and the connected system is calibrated to the same calibration number;
[0009] S4: Establish an LBM model of bound water adsorption-desorption. In the digital core established in S2, the lattice Boltzmann (LBM) is used to simulate the adsorption and desorption of shale bound water to obtain a shale digital core containing original bound water.
[0010] S5: Establish a Lattice Boltzmann Method (LBM) model of CO2 flow under irreducible water conditions, using the digital core of the water-bearing shale obtained in S4 as the medium and using the newly established flow model to simulate the flow of CO2 in the water-bearing shale;
[0011] S6: Based on the average flow velocity at the outlet surface obtained by simulation in S5, Darcy's law is used to calculate the apparent permeability of carbon dioxide flow in shale.
[0012] Furthermore, step S2 uses a binarization algorithm based on the shale scanning electron microscope image to convert the scanning electron microscope grayscale image into a binary image through a segmentation threshold, and then reconstructs the shale three-dimensional digital core using a random algorithm based on the binary image.
[0013] Furthermore, the thirteen-direction algorithm in step S3 includes:
[0014] The algorithm scans in the x, y and z directions until it reaches the center of the pore (initial grid point); then, the algorithm scans forward and backward in thirteen directions along a given direction in the three-dimensional computational domain until it reaches a grid point occupied by a solid body, and calculates the average value of the distance from the pore center to the solid body in each direction as the pore length; thus, the pore size distribution in the core is obtained.
[0015] Furthermore, the connectivity marking algorithm in step S3 includes:
[0016] By scanning the entire digital domain, meshes containing connected pores are grouped together based on the connectivity between each mesh. The algorithm considers a total of 26 directions in three dimensions (3D) (1 and the other 26 points form 26 directions) as potential connection directions. The algorithm scans in the x, y, and z directions to obtain a mesh point occupied by a pore.
[0017] The algorithm then scans the 13 adjacent directions of the grid point in turn. The algorithm then marks the grid points according to the following logic:
[0018] 1. If the grid is not connected to any of the 13 adjacent grid points (no grid point is marked), a new disconnected calibration number is assigned to the grid point;
[0019] 2. If the grid point is connected to one or more of the 13 adjacent grid points (at least one adjacent grid point is labeled), find the point with the smallest number of labels among these connected grid points and label all grid points with the smallest number of labels. In the first scan, the algorithm identifies the connected pores and assigns labels (the label numbers of connected pores may be different from each other). After several scans (the number of scans of the domain depends on the size of the domain and the pore connectivity), the algorithm unifies these labels using the smallest label number. Eventually, the same minimum label belongs to a connected pore system, while each unconnected pore system receives a different label number.
[0020] Furthermore, the LBGK model of D3Q19 is selected to establish the model in step S4, and the lattice Boltzmann equation is expressed as:
[0021]
[0022] In the formula, f i (x,t) is the distribution function, is the density distribution function in equilibrium, Δt is the time step of each iteration; τ is the relaxation time, and υ is the kinematic viscosity. e About the speed of sound hour, The lattice velocities in the nineteen directions are:
[0023]
[0024] The equilibrium distribution function is defined as follows:
[0025]
[0026] where ω i Represents the weight of the speed in each direction, and the expression is:
[0027]
[0028] Macroscopic properties such as density ρ and velocity u are determined by the moment of the particle distribution function, which is expressed as:
[0029]
[0030] Furthermore, the step S4 uses a pseudo potential model to simulate the phase separation behavior of water, wherein the mixed form of the interaction between water is:
[0031]
[0032] where ψ is the field potential. G ffis the Green function, a parameter similar to temperature, used to determine the interaction force between fluids. If the value is negative, the magnitude of the negative value indicates the strength of the attraction between particles; if the value is positive, the magnitude of the positive value corresponds to the magnitude of the repulsion between particles. Thermodynamic inconsistency is suppressed by adjusting the coefficient β.
[0033] Furthermore, the inter-particle force between the solid and liquid phases in step S4 is expressed as:
[0034]
[0035] By adjusting the density ρ of the liquid on the solid wall w To obtain the interaction forces between different solid-liquid phases, s(x+e i ) is an index function, where 0 represents a fluid node and 1 represents a solid node. Note that ρ w It is just a free parameter used to set different wall properties.
[0036] Furthermore, the total interaction force F in step S4 is tatal The modified equilibrium distribution function and equilibrium speed are introduced into the LBM system, as shown in the formula:
[0037]
[0038] Through Taylor expansion, the momentum equation of the fluid can be restored, and its pressure equation is as follows:
[0039]
[0040] Furthermore, the step S4 uses the real state equation of water to replace the ideal state equation, which can not only improve the water-gas density ratio, but also reduce the false velocity near the phase interface, which is expressed as follows:
[0041]
[0042] Here G ff is no longer a constant like temperature, it just keeps the sign of the square root always positive. In the above formula, b=0.18727RT c / p. Here, the proposed comparison temperature (T r =T / T c ) or quasi-comparative pressure (P r =P / P c ) is the key to simulating the distribution of bound water.
[0043] Furthermore, in step S4, the left and right boundaries are set to a periodic boundary format, and the upper and lower boundaries are set to a rebound boundary format.
[0044] Furthermore, the model established in step S5 uses the D3Q19 multi-relaxation lattice Boltzmann model (MRT-LBM), and its collision migration process expression is as follows:
[0045]
[0046] In the formula, f i (x,t) is the particle distribution function, is the corresponding equilibrium distribution function. δt is the time increment, c i are the nineteen speeds in the D3Q19 model, expressed as:
[0047]
[0048] M is an orthogonal matrix, S is a relaxation time matrix, expressed as follows:
[0049]
[0050] S=diag[0,s e ,s ε ,0,s q ,0,s q ,0,s q ,s v ,s π ,s v ,s π ,s v ,s v ,s v ,s t ,s t ,s t ]
[0051] In the formula, s v is the relaxation time τ s Related parameters (s v =1 / τ s ), the relaxation time is a parameter related to the kinematic viscosity, and the expression is:
[0052]
[0053] Where μ and ρ represent gas viscosity and density respectively.
[0054] Furthermore, the equilibrium distribution function expression in step S5 is:
[0055]
[0056] In the formula, ω i is the weight in the i direction, c s is the speed of sound in the lattice, equal to
[0057] Furthermore, in step S5, macroscopic quantities such as density, velocity, pressure, etc. are obtained using the following formula:
[0058]
[0059] Furthermore, the step S5 adopts the DM boundary condition at the boundary, and its expression is:
[0060]
[0061] In the formula, represents the distribution function of all particles that flow toward the solid wall after collision, ∑ j ω j is the total weight of particles flowing toward the solid wall.
[0062] Furthermore, in step S5, a Bosanquet-type effective viscosity model is used to correct the influence of the Knudsen layer, which takes into account the average effective viscosity on the cross section, and the expression is as follows:
[0063]
[0064] In the formula, μ and μ e are the dynamic viscosity and effective dynamic viscosity of the gas respectively.
[0065] Furthermore, the relaxation time (τ s ) and the Knudsen number (Kn) are related as follows:
[0066]
[0067] Where H is the local characteristic length in the shale nanoporous medium.
[0068] Furthermore, in step S6, the average flow rate of carbon dioxide at the outlet cross section obtained by simulation in S5 is used to calculate the permeability K of carbon dioxide gas using Darcy's formula, and the expression is as follows:
[0069]
[0070] In the formula, is the average velocity of carbon dioxide at the cross section; μ is the viscosity of carbon dioxide gas; L is the length of the shale nanoporous medium; Δp is the pressure difference before and after the fluid passes through the core.
[0071] The present invention provides a method for determining the apparent permeability of carbon dioxide injected into shale gas considering the distribution of bound water. The method uses a Monte Carlo random algorithm to reconstruct a shale nano digital core. The reconstructed digital core is used as a platform. The LBM model for simulating the occurrence of bound water in shale is established based on the D3Q19-LBGK model, the pseudo-potential model, etc., to simulate the shale digital core containing original bound water. Then, the D3Q19-MRT model, the viscosity correction model, and the DM boundary are coupled to establish a set of LBM models considering the influence of the original bound water in the shale. The flow rate of carbon dioxide in water-bearing shale is simulated, and finally, the apparent permeability of carbon dioxide injected into the water-bearing shale is calculated according to the Darcy formula. The method can calculate the gas permeability in water-bearing shale under different formation conditions, and has the advantages of high flexibility, good repeatability, and low cost. BRIEF DESCRIPTION OF THE DRAWINGS
[0072] Figure 1 is a flow chart of a method according to an embodiment of the present invention;
[0073] Figure 2 This is a schematic diagram of a scanning electron microscope image of shale in an embodiment of the present invention;
[0074] Figure 3 This is a schematic diagram of a shale binarization image according to an embodiment of the present invention;
[0075] Figure 4 A schematic diagram of a three-dimensional digital core of shale reconstructed according to an embodiment of the present invention. DETAILED DESCRIPTION
[0076] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail in conjunction with the following embodiments. It should be understood that the specific embodiments described herein are only used to conclude the present invention and are not used to limit the present invention.
[0077] like Figure 1 A method for determining the apparent permeability of shale gas carbon dioxide injection considering the distribution of irreducible water comprises the following steps:
[0078] S1: Obtain shale scanning electron microscope images, porosity and pore size distribution. Shale scanning electron microscope images are measured using the industry standard "SYT 5162-1997 Rock Sample Scanning Electron Microscope Analysis Method";
[0079] S2: Based on the shale SEM image, a shale digital core reconstruction method with porosity as a constraint is established to reconstruct the shale three-dimensional digital core. The specific process is to use a binarization algorithm based on the shale SEM image, and convert the SEM grayscale image into a black and white binary image through segmentation thresholds to ensure that the porosity of the binary image is basically consistent with the core porosity. Then, based on the binary image, a Monte Carlo random algorithm is used to reconstruct the shale three-dimensional digital core;
[0080] S3: The thirteen-direction average algorithm is used to calculate the pore size distribution, and then the connectivity labeling algorithm is used to determine the connectivity of the pores, and the connected system is calibrated to the same calibration number; by scanning the entire digital domain, the grids containing connected pores are grouped together according to the connectivity between each grid. The algorithm considers a total of 26 directions (1 and the other 26 points form 26 directions) in three dimensions (3D) as potential connection directions; the algorithm scans in the x, y and z directions to obtain a grid point occupied by a pore; then, the algorithm scans the 13 adjacent directions of the grid point in turn; then, the algorithm marks the grid point according to the following logic: 1. If the grid is not connected to any of the 13 adjacent grid points (no grid point is marked), a new disconnected calibration number is assigned to the grid point; 2. If the grid point is connected to one or more of the 13 adjacent grid points (at least one adjacent grid point is marked), the point with the smallest number of labels among these connected grid points is found, and all grid points are marked with the minimum number of labels; in the first scan, the algorithm identifies connected pores and assigns labels (the label numbers of connected pores may be different from each other); after several scans (the number of scans of the domain depends on the size of the domain and the pore connectivity), the algorithm uses the smallest label number to unify these labels; eventually, the same minimum label belongs to a connected pore system, while each unconnected pore system receives a different label number.
[0081] S4: Establishing an LBM model of bound water adsorption-desorption, using the lattice Boltzmann (LBM) to simulate the adsorption and desorption of bound water in shale in the digital core established in S2, and obtaining a shale digital core containing original bound water; the specific process is as follows, the bound water adsorption-desorption LBM model of the present invention adopts the LBGK model of D3Q19, and the lattice Boltzmann equation is expressed as:
[0082]
[0083] In the formula, f i (x,t) is the distribution function, is the density distribution function in equilibrium, Δt is the time step of each iteration; τ is the relaxation time, and υ is the kinematic viscosity. e About the speed of sound hour, The lattice velocities in the nineteen directions are:
[0084]
[0085] The particle equilibrium distribution function is defined as follows:
[0086]
[0087] where ω i Represents the weight of the speed in each direction, and the expression is:
[0088]
[0089] Macroscopic properties such as density ρ and velocity u are determined by the moment of the particle distribution function, which is expressed as:
[0090]
[0091] The phase separation behavior of water is simulated using a pseudopotential model, where the interactions between water molecules take the form of a mixture:
[0092]
[0093] where ψ is the field potential. G ff is the Green function, a parameter similar to temperature, used to determine the interaction force between fluids; if the value is negative, the magnitude of the negative value indicates the strength of the attraction between particles; if the value is positive, the magnitude of the positive value corresponds to the magnitude of the repulsion between particles. Thermodynamic inconsistency is suppressed by adjusting the coefficient β;
[0094] The interparticle force between water molecules and the solid-liquid phase of the solid wall is expressed as:
[0095]
[0096] By adjusting the density ρ of the liquid on the solid wall w To obtain the interaction forces between different solid-liquid phases, s(x+e i ) is an index function, where 0 represents a fluid node and 1 represents a solid node; Note that ρ w It is just a free parameter used to set different wall properties;
[0097] The total interaction force F is obtained by superimposing the liquid-liquid force between water molecules and the solid-liquid force between the solid wall and water molecules. tatal , and the total interaction force F tatal The modified equilibrium distribution function and equilibrium speed are introduced into the LBM system, as shown in the formula:
[0098]
[0099] Through Taylor expansion, the momentum equation of the fluid can be restored, and its pressure equation is as follows:
[0100]
[0101] In order to make the simulation results more accurate, the Carnahan-Starling (CS) state equation is used to replace the ideal state equation used in the conventional model, which can not only improve the water-gas density ratio, but also reduce the false velocity near the phase interface. It is expressed as follows:
[0102]
[0103] Here G ff is no longer a constant like temperature, it just keeps the sign of the square root always positive; in the above formula, b=0.18727RT c / p;
[0104] The left and right boundaries are set to the periodic boundary format, and the upper and lower boundaries are set to the rebound boundary format. The LBM bound water adsorption-desorption simulation is carried out to obtain the shale digital core under the original water-containing condition.
[0105] S5: Establish a lattice Boltzmann method (LBM) model of carbon dioxide flow under bound water conditions. The water-bearing shale digital core obtained in S4 is used as the medium, and the newly established flow model is used to simulate the flow of carbon dioxide in the water-bearing shale. The specific process is as follows: The D3Q19 multi-relaxation lattice Boltzmann model (MRT-LBM) is selected to characterize the flow of carbon dioxide gas in the water-bearing shale core, and its collision migration process expression is as follows:
[0106]
[0107] In the formula, f i (x,t) is the particle distribution function, is the corresponding equilibrium distribution function in the i direction; δt is the time increment, c i are the nineteen speeds in the D3Q19 model, expressed as:
[0108]
[0109] M is an orthogonal matrix, S is a relaxation time matrix, expressed as follows:
[0110]
[0111] S=diag[0,s e ,s ε ,0,s q ,0,s q ,0,s q ,s v ,s π ,s v ,s π ,s v ,sv ,s v ,s t ,s t ,s t ]
[0112] In the formula, s v is the relaxation time τ s Related parameters (s v =1 / τ s ), the relaxation time is a parameter related to the kinematic viscosity, and the expression is:
[0113]
[0114] Where μ and ρ represent gas viscosity and density respectively.
[0115] The particle equilibrium distribution function expression is:
[0116]
[0117] In the formula, ω i is the weight in the i direction, c s is the speed of sound in the lattice, equal to
[0118] Macroscopic parameters such as density (ρ), velocity (u), pressure (p), etc. are obtained using the following formulas:
[0119]
[0120] For particles that flow toward the solid after collision, the DM boundary condition is used at the solid boundary, and its expression is:
[0121]
[0122] In the formula, represents the distribution function of all particles that flow toward the solid wall after collision, ∑ j ω j is the total weight of particles flowing toward the solid wall;
[0123] In order to ensure the accuracy of the simulation results, the Bosanquet-type effective viscosity model is introduced to correct the influence of the Knudsen layer, which takes into account the average effective viscosity on the cross section. The expression is as follows:
[0124]
[0125] In the formula, μ and μ e are the dynamic viscosity and effective dynamic viscosity of the gas respectively;
[0126] The relaxation time (τ) in the model after correcting the viscosity s) and the Knudsen number (Kn) are related as follows:
[0127]
[0128] Where H is the local characteristic length in shale nanoporous media;
[0129] LBM was used to carry out LBM simulation of carbon dioxide injection flow in water-bearing shale, and the fluid velocity in the shale outlet section was statistically obtained.
[0130] S6: Based on the average flow velocity at the outlet surface simulated in S5, the apparent permeability of carbon dioxide flow in shale is calculated using Darcy's law, and the apparent permeability of carbon dioxide flow in the shale core under the irreducible water saturation can be obtained; Darcy's law is expressed as follows:
[0131]
[0132] In the formula, is the average velocity of carbon dioxide at the cross section; μ is the viscosity of carbon dioxide gas; L is the length of the shale nanoporous medium; Δp is the pressure difference before and after the fluid passes through the core.
[0133] Embodiment 1:
[0134] Take a shale core from a certain block, and measure the porosity of the rock sample according to step S1 to be 9.08%. Through the scanning electron microscope image, the porosity of the binary image obtained at a threshold of 113 is 9.04%. The scanned image is as follows Figure 2 and two Valued pictures such as Figure 3 ;
[0135] According to step 2, by scanning the binary processed image, traversing the nodes in the x, y, and z directions, the Monte Carlo random algorithm is used to reconstruct the shale digital core. The size of the reconstructed core is 200×200×200, such as Figure 4 As shown, the porosity of the digital core is 9.11%, and the error of the constructed core porosity is 0.3%.
[0136] According to step 3, the pore size distribution is statistically analyzed using the 13-direction average algorithm; the connected pore algorithm is used to calibrate the connected pores to distinguish between isolated pores and connected pores in the reconstructed core.
[0137] According to step 4, a shale digital core containing original irreducible water is simulated, and the number of water-containing grids is counted and divided by the total number of grids to obtain the irreducible water saturation of 0.10.
[0138] According to step 5, a flow simulation of carbon dioxide gas injected into shale containing bound water is carried out to obtain the flow velocity of carbon dioxide at the outlet surface.
[0139] According to step 6, the apparent permeability of the injected carbon dioxide gas in the shale at the water saturation is calculated according to Darcy's formula to be 0.00071 mD.
[0140] The present invention provides a method for determining the apparent permeability of carbon dioxide injected into shale gas considering the distribution of bound water. The method uses a Monte Carlo random algorithm to reconstruct a shale nano digital core. The reconstructed digital core is used as a platform. The LBM model for simulating the occurrence of bound water in shale is established based on the D3Q19-LBGK model, the pseudo-potential model, etc., to simulate the shale digital core containing original bound water. Then, the D3Q19-MRT model, the viscosity correction model, and the DM boundary are coupled to establish a set of LBM models considering the influence of the original bound water in the shale. The flow rate of carbon dioxide in water-bearing shale is simulated, and finally, the apparent permeability of carbon dioxide injected into the water-bearing shale is calculated according to the Darcy formula. The method can calculate the gas permeability in water-bearing shale under different formation conditions, and has the advantages of high flexibility, good repeatability, and low cost.
[0141] The above description is not intended to impose any form of limitation on the present invention. Although the present invention has been disclosed through the above embodiments, it is not intended to limit the present invention. Any technician familiar with the profession can make some changes or modifications to equivalent embodiments of equivalent changes using the technical contents disclosed above without departing from the scope of the technical solution of the present invention. However, any simple modification, equivalent change and modification made to the above embodiments based on the technical essence of the present invention without departing from the content of the technical solution of the present invention still falls within the scope of the technical solution of the present invention.
Claims
1. A method for determining the apparent permeability of shale gas CO2 injection considering the distribution of bound water, characterized in that: The following steps are involved: S1: Obtain shale SEM images, porosity and pore size distribution; S2: Based on the shale scanning electron microscope image, a shale digital core reconstruction method with porosity as a constraint is established to reconstruct the shale three-dimensional digital core; S3: The thirteen-direction average algorithm is used to calculate the pore size distribution, and then the connectivity labeling algorithm is used to determine the connectivity of the pores, and the connected system is calibrated to the same calibration number; S4: Establish an LBM model of bound water adsorption-desorption. In the digital core established in S2, the lattice Boltzmann LBM is used to simulate the adsorption and desorption of bound water in shale to obtain a shale digital core containing original bound water. S5: Establish a Lattice Boltzmann Method (LBM) model of carbon dioxide flow under bound water conditions, using the digital core of the water-bearing shale obtained in S4 as the medium and using the newly established flow model to simulate the flow of carbon dioxide in the water-bearing shale; S6: Based on the average flow velocity at the outlet surface obtained by simulation in S5, Darcy's law is used to calculate the apparent permeability of carbon dioxide flow in shale.
2. The method for determining the apparent permeability of shale gas carbon dioxide injection considering the distribution of irreducible water according to claim 1, characterized in that: The step 2 is based on shale scanning electron microscopy, uses a binarization algorithm to set a threshold, converts the grayscale image into a binary image, and reconstructs the shale digital core using a Monte Carlo random algorithm on the premise of ensuring that the porosity of the two is consistent.
3. The method for determining apparent permeability of shale gas carbon dioxide injection considering irreducible water distribution according to claim 1, characterized in that: In step S3, the thirteen-direction algorithm is used to calculate the pore size distribution of the shale core and the interconnection marking algorithm is used to calibrate the pore connectivity; The Thirteen Directions Method specifically includes: The algorithm scans in the x, y and z directions until it reaches the center of the pore, i.e., the initial grid point. Then, the algorithm scans forward and backward in thirteen directions along a given direction in the three-dimensional computational domain until it reaches a grid point occupied by a solid. The distance from the pore center to the solid in each direction is calculated as the average value as the pore length, thereby obtaining the pore size distribution in the core. The connected labeling algorithm includes: By scanning the entire digital domain, the algorithm groups grids containing connected pores together based on the connectivity between each grid. The algorithm considers a total of three dimensions as potential connection directions. The algorithm scans in the x, y, and z directions to obtain a grid point occupied by a pore. The algorithm then scans the 13 neighboring directions of this grid point in turn. The algorithm then labels the grid points according to the following logic: (1) If the grid is not connected to any of the 13 adjacent grid points, that is, no grid point is marked, a new disconnected calibration number is assigned to the grid point; (2). If the grid point is connected to one or more of the 13 adjacent grid points, then find the point with the smallest number of labels among these connected grid points and label all grid points with the smallest number of labels. In the first scan, the algorithm identifies the connected pores and assigns labels. After several scans, the algorithm uses the smallest label number to unify these labels. Ultimately, the same minimum label belongs to a connected pore system, while each unconnected pore system receives a different label number.
4. The method for determining apparent permeability of shale gas carbon dioxide injection considering irreducible water distribution according to claim 1, characterized in that: In step S4, a single-component multiphase pseudo-potential LBM model is used to simulate the adsorption of original bound water in shale on the wall surface; specifically: First, the LBGK model of D3Q19 is selected to establish the model. The expression of the lattice Boltzmann equation is: In the formula, f i (x, t) is the distribution function, f i eq (x, t) is the density distribution function at equilibrium, Δt is the time step of each iteration; τ is the relaxation time, and υ is the kinematic viscosity. e About the speed of sound hour, The lattice velocities in the nineteen directions are: The equilibrium distribution function is defined as follows: where ω i Represents the weight of the speed in each direction, and the expression is: Macroscopic properties such as density ρ and velocity u are determined by the moment of the particle distribution function, which is expressed as: The phase separation behavior of water is simulated using a pseudopotential model, where the interactions between water molecules take the form of a mixture: where ψ is the field potential; G ff is the Green function, a parameter similar to temperature, used to determine the interaction force between fluids; if the value is negative, the magnitude of the negative value indicates the strength of the attraction between particles; if the value is positive, the magnitude of the positive value corresponds to the magnitude of the repulsion between particles; thermodynamic inconsistency is suppressed by adjusting the coefficient β; The expression of the interparticle force between the solid and liquid phases in the pseudo-potential model is: By adjusting the density ρ of the liquid on the solid wall w To obtain the interaction forces between different solid-liquid phases, s(x+e i ) is an index function, where 0 represents a fluid node and 1 represents a solid node; Note that ρ w It is just a free parameter used to set different wall properties; The total interaction force F tatal The modified equilibrium distribution function and equilibrium speed are introduced into the LBM system, as shown in the formula: Through Taylor expansion, the momentum equation of the fluid can be restored, and its pressure equation is as follows: In order to improve the simulation accuracy, the real state equation of water is used to replace the ideal state equation, which can not only improve the water-gas density ratio, but also reduce the false velocity near the phase interface. It is expressed as follows: Here G ff It is no longer a constant like temperature, it just keeps the sign of the square root always positive; in the above formula, a = 0.4963R 2 T c 2 / p c , b=0.18727RT c / p; Set the left and right boundaries to periodic boundary format and the upper and lower boundaries to bounce boundary format.
5. The method for determining apparent permeability of shale gas carbon dioxide injection considering irreducible water distribution according to claim 1, characterized in that: The step S5 uses a multi-relaxation lattice Boltzmann model to simulate the flow behavior of carbon dioxide in shale containing bound water; specifically: First, the established model uses the D3Q19 multi-relaxation lattice Boltzmann model MRT-LBM, and its collision migration process expression is as follows: f i (x+c i δt,t+δt)-f i (x,t)=-(M -1 SM)[f i (x,t)-f i (eq) (x,t)] In the formula, f i (x, t) is the particle distribution function, f i (eq) (x, t) is the corresponding equilibrium distribution function; δt is the time increment, c i are the nineteen speeds in the D3Q19 model, expressed as: M is an orthogonal matrix, S is a relaxation time matrix, expressed as follows: S=diag[0,s e ,s ε ,0,s q ,0,s q ,0,s q ,s v ,s π ,s v ,s π ,s v ,s v ,s v ,s t ,s t ,s t ] In the formula, s v is the relaxation time τ s The relevant parameters, s v =1 / τ s , relaxation time is a parameter related to kinematic viscosity, and its expression is: In the formula, μ and ρ represent gas viscosity and density respectively; The particle equilibrium distribution function expression in the MRT-LBM model is: In the formula, ω i is the weight in the i direction, c s is the speed of sound in the lattice, equal to In this model, the macroscopic quantities density ρ, velocity ρu, and pressure are obtained using the following formulas: The DM boundary condition is used at the solid boundary, and its expression is: In the formula, represents the distribution function of all particles that flow toward the solid wall after collision, ∑ j ω i is the total weight of particles flowing toward the solid wall; The Bosanquet-type effective viscosity model is used to correct the influence of the Knudsen layer, which takes into account the average effective viscosity on the cross section. The expression is as follows: In the formula, μ and μ e are the dynamic viscosity and effective dynamic viscosity of the gas respectively; The relaxation time τ in the model s The relationship with the Knudsen number Kn is as follows: Where H is the local characteristic length in the shale nanoporous medium.
6. The method for determining apparent permeability of shale gas carbon dioxide injection considering irreducible water distribution according to claim 1, characterized in that: The step S6 uses the average flow rate of carbon dioxide at the outlet cross section obtained by simulation in S5 and Darcy's formula to calculate the permeability K of carbon dioxide gas, which is expressed as follows: In the formula, is the average velocity of carbon dioxide at the cross section; μ is the viscosity of carbon dioxide gas; L is the length of the shale nanoporous medium; Δp is the pressure difference before and after the fluid passes through the core.