Quantitative characterization method for microscopic characteristics of CO2 flooding in gas reservoirs based on CO2 water-rock reaction and diffusion effect

Abstract

The present application belongs to the technical field of oil and gas field development and carbon dioxide geological storage, and particularly relates to a gas reservoir injection and drive micro-characteristic quantitative characterization method based on water-rock reaction and diffusion effect. First, three-dimensional CT images before and after reaction are collected to construct a real pore tetrahedron calculation grid, a Navier-Stokes flow and Maxwell-Stefan convection-diffusion coupling model is established, a viscosity polynomial and a deviation factor density correction are introduced, a mass transfer dominant zone is divided by a characteristic cross-section flux, a mixed layer front edge proportion and a diffusion influence proportion are calculated, and the quantitative coupling influence of parameters before and after reaction is compared to provide micro-quantitative support for injection, production and storage evaluation. The present application can solve the problem of inaccurate characterization of the real displacement mechanism in the existing micro-pore scale simulation characterization method.

Description

Technical Field

[0001] This invention belongs to the field of oil and gas field development and carbon dioxide geological storage technology, and particularly relates to... Gas reservoir injection based on water-rock reaction and diffusion effects A quantitative characterization method for microscopic features. Background Technology

[0002] As global climate change becomes increasingly severe, carbon capture, utilization, and storage (CCUS) technology has become a key pathway to achieve large-scale greenhouse gas emission reduction. Among numerous storage options, depleted oil and gas reservoirs are considered the most suitable due to their complete sealing structure, good sealing performance, and detailed geological exploration foundation. Early geological sites that have been sealed off.

[0003] The technology can simultaneously achieve the dual benefits of long-term geological sequestration of greenhouse gases and increased natural gas production, and has already attracted widespread attention globally, leading to the launch of several demonstration projects. To effectively guide... In engineering practice, existing research often uses multiphysics coupling to explore the impact of parameters such as injection duration, matrix-fracture ratio, and permeability heterogeneity on extraction efficiency from a macroscopic numerical simulation perspective. However, Darcy's law often characterizes complex pore geometry as permeability, neglecting the complex flow characteristics and mass transfer mechanisms at the microscopic pore scale. In microscopic numerical simulations, traditional methods are mostly based on Fick's law to describe gas displacement processes. This law is prone to significant errors when dealing with real non-ideal (high temperature and high pressure) mixtures, making it difficult to accurately characterize the complex diffusion coupling effects between multi-component systems. In contrast, diffusion models based on the Maxwell-Stefan framework can more realistically reflect the interdiffusion behavior between components in complex microscopic pore mass transfer.

[0004] Furthermore, existing microscopic pore-scale simulation and characterization methods still have two significant limitations: Firstly, in the treatment of fluid properties, current studies often use simple linear weighting of fluid density and viscosity, failing to fully consider the nonlinear evolution of fluid properties under high temperature and pressure conditions in reservoirs as a result of dynamic mixing of multiple components, leading to deviations in the prediction of convection-diffusion mechanisms; secondly, actual reservoirs such as carbonate rocks experience long-term contact... After formation water is saturated, significant dynamic water-rock reactions such as dissolution and precipitation occur, leading to changes in the topology of the pore network. However, existing methods rarely integrate this real dynamic evolution of rock microstructure with complex multi-component mass transfer processes.

[0005] In summary, existing methods lack a high-precision quantitative characterization of the impact of non-ideal fluid property changes and dynamic evolution of micropore structure on gas reservoir extraction and production, resulting in inaccurate characterization of the true displacement mechanism. Summary of the Invention

[0006] The technical problem solved by this invention is to provide a method based on Gas reservoir injection based on water-rock reaction and diffusion effects A quantitative characterization method for displacement microstructures is proposed to address the inaccuracy of existing micro-pore scale simulation methods in representing the true displacement mechanism.

[0007] The basic solution provided by this invention: based on Gas reservoir injection based on water-rock reaction and diffusion effects Methods for quantitative characterization of microscopic features include: S1: Collect target core samples Two sets of real three-dimensional CT image sequences were obtained before and after the water-rock reaction. The three-dimensional CT image sequence before the reaction was binarized and segmented to distinguish the rock skeleton and pore space, thus obtaining three-dimensional pore volume data. Based on the three-dimensional pore volume data, the pore network was extracted, and the pore space isosurface was extracted to generate an initial surface mesh. S2: Perform topology repair and mesh reconstruction operations on the initial surface mesh to construct a tetrahedral computational mesh; and use the tetrahedral computational mesh as the computational domain to construct a multiphysics coupling model that couples the Navier-Stokes flow equations with the Maxwell-Stefan multi-component convection-diffusion equations. S3: Dynamically correct the viscosity and density of the mixed fluid in the multiphysics coupling model, and solve the problem based on the corrected multiphysics coupling model. Obtain the changes in mole fraction of different fluid components inside the target core at different injection time points, the convection flux and diffusion flux at different cross sections, and the thickness ratio of the mixing layer front. Set up a no-diffusion control condition and calculate the proportion of diffusion mechanism influence by comparing the component mole fractions of the baseline condition and the no-diffusion condition. S4: For the acquired post-reaction 3D CT image sequences, steps S1 to S3 are processed independently, and the changes in quantitative characteristic parameters such as component mole fraction, convective flux, diffusion flux, mixing layer front thickness ratio, and diffusion mechanism influence ratio before and after the reaction are compared to quantify the changes. Water-rock reaction and diffusion effects on gas reservoir injection The extent to which the extraction mechanism is affected.

[0008] Furthermore, S1 includes: S1-1: Micron-level industrial CT was used to perform separate operations on the target core. Two whole-core scans were performed before and after the saturated formation water-rock reaction to obtain three-dimensional CT image sequences before and after the reaction. For the three-dimensional CT image sequence before the reaction, a cube-shaped computational subdomain without edge artifacts was cropped from the original cylindrical scan image, and invalid black areas at the scan boundary were removed. S1-2: The cropped 3D image sequence is subjected to median filtering noise reduction using a 3×3×3 sliding window. The grayscale values ​​of the filtered image satisfy the following:

[0009] in, The original 3D CT image in voxel coordinates The grayscale value at that location; This represents the grayscale value at the corresponding location after noise reduction. The median value of the neighborhood grayscale values ​​is calculated. These represent the spatial coordinate translations of the sliding window along the x, y, and z axes, respectively. S1-3: Based on the globally optimal segmentation threshold t, the denoised image is binarized to distinguish between the rock skeleton and the pore space, resulting in binarized three-dimensional pore volume data, expressed as:

[0010] in, The voxel label values ​​of the binarized image. Indicates the rock skeleton region. This refers to the porous space region; S1-4: Based on the binarized 3D porous volume data, extract the pore network and extract the isosurfaces of the pore network to generate an initial surface mesh.

[0011] Furthermore, S2 includes: S2-1: Traverse all initial surface meshes generated by S1, identify defects including non-manifold edges and vertices, self-intersecting triangular faces, open boundary holes, and isolated porous surfaces not connected to the main channel, and mark the location and type of defect elements; S2-2: Delete all surface mesh elements corresponding to isolated pores that are not connected to the main channel of the injection end-extraction end, and seal the open boundary voids on the surface of the remaining effective pores. S2-3: For the triangular elements with self-intersecting triangular faces in the initial surface mesh, the face intersection is eliminated by fine-tuning the local vertex coordinates and reconnecting with the triangular elements. The vertex position is optimized by using the Laplacian smoothing algorithm to eliminate the step-like pixelation error and obtain the smoothed initial surface mesh. S2-4: Using the smoothed surface mesh of the initial surface mesh after smoothing as the boundary constraint, a volume mesh is generated. This mesh is then imported into the finite element simulation software to refine the mesh locally in the narrow channel region, resulting in a tetrahedral computational mesh for the entire pore domain.

[0012] Furthermore, S2 also includes: S2-5: Using a tetrahedral computational mesh as the computational domain, a multiphysics coupled model is established in finite element simulation software, which combines the Navier-Stokes flow equations with the Maxwell-Stefan convection-diffusion equations. The Navier-Stokes flow equations include: The equation for the conservation of mass is expressed as follows:

[0013] The equation for the conservation of momentum is expressed as follows:

[0014] in, The density of the gas mixture, The velocity vector of the fluid. For fluid pressure, The dynamic viscosity of the gas mixture. For unit tensors, For volume forces, It is a vector differential operator; The Maxwell-Stefan convection-diffusion equations include: The equation for the conservation of mass of the components is expressed as follows:

[0015] The equation for the relative diffusion flux of the components is expressed as follows:

[0016] The diffusion driving force equation is expressed as follows:

[0017]

[0018] in, Let t be the mass fraction of component i, and t be the injection time. Source terms generated for components, Let be the relative diffusion flux of component i. The Maxwell-Stefan diffusion coefficient is... Where is the thermal diffusivity, and T is the system temperature. As a diffusion driving force, Let k be the mole fraction of component k. Total pressure; Let k be the molar mass of component k. The average molar mass of the gas mixture; Let be the molar mass of the i-th component; S2-6: Define the two end faces of the core along the flow direction as the injection end and the production end, including: Rock solid wall: no-slip boundary and zero diffusion flux boundary , where n is the out-of-plane normal unit vector; Injection end: constant velocity boundary and constant injection component boundary ; For a constant inlet flow rate, The components are mole fraction; Production end: Constant pressure boundary and zero diffusion flux boundary , This creates relative pressure on exports.

[0019] Furthermore, S3 includes: S3-1: Based on the fundamental viscosity calculated using Wilke's mixing rule, a polynomial correction factor related to the mole fraction of components is introduced to dynamically correct the non-ideal mixing viscosity under high temperature and high pressure. The expression is as follows:

[0020]

[0021] in, To correct the viscosity of the fluid mixture, The viscosity of the fluid mixture calculated using the Wilke formula. These are the mole fractions of component i, component j, and component k, respectively. It is a constant correction factor; S3-2: The deviation factor is constructed as a function of temperature, pressure, and component mole fraction. Substituting this into the ideal gas law completes the density correction of the mixed fluid. The expression is:

[0022] in, For the density of the fluid mixture, For pressure, The average molar mass of the fluid mixture. The gas deviation factor is a function of temperature, pressure, and mole fraction of components. This is the universal gas constant. This refers to the system temperature.

[0023] Furthermore, S3 also includes: S3-3: Using the viscosity and density of the mixed fluid as modified in S3-1 and S3-2 as parameters, the Navier-Stokes flow equation of the multiphysics coupling model is input for solution to obtain the steady-state velocity field and pressure field; then, using the steady-state velocity field as a fixed flow field, the Maxwell-Stefan convection-diffusion equation is substituted into the solution to output the spatial distribution of the mass fraction and mole fraction of each component in the whole pore domain at different injection times; S3-4: Based on the results output from S3-3, extract the characteristic cross-sectional data perpendicular to the flow direction, and calculate the convective flux and diffusion flux at this cross-section at different injection times. The calculation expression is:

[0024] diffusion flux The calculation expression is:

[0025] Based on the convective and diffusion fluxes at this cross section at different injection times, flux-time evolution curves are generated. S3-5: Based on the flux time evolution curve, identify the flux boundary point where the convective flux and diffusion flux are equal in value, and determine the position and length of the boundary interface along the z-axis corresponding to this boundary point. The spatial division is based on the interface, with the upstream region designated as the convection-dominant region and the downstream region as the diffusion-dominant region. The output zoning results are as follows: parameter; S3-6: Based on the defined convection-dominant and diffusion-dominant regions, the overall mixing layer is divided into a rear mixing layer and a leading mixing layer. The proportion of the leading mixing layer is obtained by calculating the ratio of the pore volume at the leading mixing layer within the diffusion-dominant region to the total pore volume of the overall mixing layer. The expression is:

[0026]

[0027]

[0028] in, The volume of the mixed layer front is defined by taking the core injection end section as the origin of the z-axis and the flow direction as the positive z-axis direction. Let z be the z-axis coordinate of the interface where the flux and diffusion flux are equal. It is the length of the distance from the z-axis. Let i be the mole fraction of component i. The total volume of the mixed layer, The percentage of the leading edge of the hybrid layer; S3-7: Based on the displacement simulation results in S3-3 as the baseline simulation results, a diffusion-free control condition was set up to obtain the whole core sample without considering diffusion. The curve showing the change of average mole fraction over time; and the location of the whole core sample in the baseline simulation results considering diffusion. The point at which the average mole fraction drops to 5% is considered the point of complete displacement. Residual core samples from the whole core at the same time point are extracted for the diffusion-free control. Average mole fraction, calculated under the baseline simulation conditions and the diffusion-free control conditions. The absolute difference in mole fractions is used to quantify the proportion of the diffusion mechanism's influence on the displacement and evacuation process. The expression is:

[0029]

[0030] in, For diffusion-free control conditions mole fraction To consider the diffusion scenario mole fraction for The time point when the mole fraction is 5%, The proportion of the influence of diffusion mechanism.

[0031] Furthermore, S4 includes: S4-1: Obtain the value from S1 The three-dimensional CT image sequence after the water-rock reaction is used as input, and the entire process from S1 to S3 is executed independently and repeatedly to obtain the complete set of results for the component mole fraction field, convection flux and diffusion flux evolution curves, mixing layer front ratio, and diffusion mechanism influence ratio of the core after the reaction. S4-2: Select characteristic sections at the same injection time and spatial location, compare the differences in three types of field parameters—component mole fraction, convection flux, and diffusion flux—of the cores before and after the reaction, and calculate the relative change rate of the corresponding parameters; S4-3: Extract the core quantitative characteristic parameters from the simulation results before and after the reaction, and calculate the changes in the proportion of the mixing layer front and the proportion of the diffusion mechanism caused by the water-rock reaction. S4-4: By combining the relative rate of change of corresponding production, the change in the proportion of the mixing layer front, and the change in the proportion of diffusion mechanism influence, the individual contribution weights of water-rock reaction to convection displacement and diffusion mass transfer through pore structure evolution are decomposed to obtain the coupling mechanism of water-rock reaction and diffusion effect, ultimately completing... Water-rock reaction and diffusion effects on gas reservoir injection Quantitative characterization of the impact of the extraction mechanism.

[0032] The technical principle of this invention lies in: targeting existing gas reservoir injection... To address the technical shortcomings of traditional microscopic characterization methods, such as insufficient accuracy in depicting multi-component diffusion, lack of consideration for the non-ideal properties of high-temperature and high-pressure fluids, and missing coupling between the evolution of pore structure and mass transfer processes in water-rock reactions, this paper uses 3D CT scan images of real rock cores as the geometric data source. Through image preprocessing, pore network extraction, surface mesh topology repair, and tetrahedral mesh reconstruction, a numerical computational domain highly consistent with the real pore morphology is constructed. A modeling approach coupling the steady-state Navier-Stokes flow equation and the Maxwell-Stefan multi-component convection-diffusion equation is adopted to fully incorporate the diffusion driving force resulting from the combined effects of concentration and pressure gradients, accurately reproducing the diffusion at the pore scale. and The interdiffusion mass transfer mechanism was investigated; the properties of the mixed fluid under high temperature and high pressure conditions in the reservoir were dynamically corrected using polynomial viscosity correction and deviation factor density correction methods based on component mole fraction, matching the nonlinear evolution characteristics of properties during multi-component mixing; based on transient simulation results, convection-diffusion dominant zones were identified through flux evolution at characteristic sections, the proportion of the mixing layer front was calculated by combining the spatial partitioning of the mixing layer, and the displacement contribution of the diffusion mechanism was quantified through diffusion-free control simulation; finally, the interdiffusion mass transfer mechanism was investigated. Independent full-process calculations and multi-dimensional parameter comparisons were conducted on two sets of core samples before and after the water-rock reaction to quantitatively characterize the coupling effect of water-rock reaction on pore structure modification and diffusion. The influence mechanism of the displacement process was investigated, and high-precision quantitative characterization of micro-displacement characteristics was achieved.

[0033] The beneficial effects of this invention are as follows: 1. This invention solves the problem that the traditional Fick diffusion law has large errors in high-temperature and high-pressure multi-component systems and cannot accurately characterize the inter-diffusion coupling effect between components: The present invention uses the Maxwell-Stefan diffusion framework to construct a convection-diffusion coupling model, which fully considers the diffusion driving force composed of the concentration gradient and the pressure gradient, and can truly reflect the diffusion at the pore scale. and The multi-component interactive diffusion behavior is more rigorously characterized, and the physical mechanism of the mass transfer process is significantly better than that of the traditional Fick model in terms of the accuracy of microscopic mass transfer simulation.

[0034] 2. This invention addresses the problem of existing methods using simple linear weighting of fluid density and viscosity, which fails to consider the nonlinear evolution of reservoir properties under high temperature and pressure, leading to displacement prediction bias: Based on Wilke's mixing rule, this invention introduces a multi-component mole fraction polynomial correction factor to achieve dynamic viscosity correction. By correcting non-ideal gas density through deviation factors related to temperature, pressure, and composition, it can dynamically track the real-time changes in the properties of the mixed fluid during displacement, effectively reducing property calculation errors under non-ideal conditions and ensuring the reliability of convection-diffusion mechanism prediction.

[0035] 3. Resolved the issue of existing methods not incorporating [this technology / method]. The problem of the disconnect between the dynamic evolution of pore structure and displacement mechanism characterization caused by water-rock reaction and the actual long-term development process of reservoirs: Based on real core CT images before and after the reaction, this invention uses independent simulation of the entire process and multi-dimensional comparison of pore structure, field distribution and core quantitative characteristics to quantitatively separate the individual contribution weights of water-rock reaction to convection displacement and diffusion mass transfer, clearly revealing the enhancement mechanism of the coupling effect between the two, and the characterization results are more consistent with the development characteristics of actual reservoirs such as carbonate rocks.

[0036] 4. The quantitative indicators proposed in this invention, such as the proportion of the hybrid layer front and the proportion of diffusion mechanism influence, compensate for the shortcomings of macroscopic simulation in revealing microscopic mechanisms. They enable quantitative assessment of the contribution of diffusion effects to the displacement and cavitation process at the pore scale, providing a basis for injection... This provides quantitative, micro-level data for optimizing injection parameters and evaluating storage effectiveness in gas reservoirs, and can effectively guide [efforts / guidelines]. On-site engineering practice. Attached Figure Description

[0037] Figure 1 This is a flowchart of an embodiment of the present invention; Figure 2 This is a flowchart of CT image processing in an embodiment of the present invention; Figure 3 This is a flowchart of surface mesh processing in an embodiment of the present invention; Figure 4 This is a convection-diffusion flux diagram in an embodiment of the present invention; Figure 5 This is a diagram showing the proportion and thickness of the hybrid layer leading edge in an embodiment of the present invention; Figure 6 This is a graph showing the change of mole fraction at different injection rates over time under diffusion conditions in an embodiment of the present invention; Figure 7 This is a graph showing the change of mole fraction at different injection rates over time under non-diffusion conditions in an embodiment of the present invention. Figure 8 This is a comparison chart of the effect of different injection rates on diffusion, replacement completion time, and baseline group in embodiments of the present invention; Figure 9 This is a graph showing the change of mole fraction over time before and after the reaction, considering diffusion and without diffusion, in an embodiment of the present invention. Figure 10 This is a diagram showing the diffusion effect and displacement completion time before and after the reaction in an embodiment of the present invention. Detailed Implementation

[0038] The following detailed description illustrates the specific implementation method: The basic implementation examples are as follows: Figure 1 As shown: Based on Gas reservoir injection based on water-rock reaction and diffusion effects Methods for quantitative characterization of microscopic features include: S1: Collect target core samples Two sets of real 3D CT image sequences before and after the water-rock reaction were obtained. The 3D CT image sequence before the reaction was binarized and segmented to distinguish the rock skeleton and pore space, thus obtaining 3D pore volume data. Based on the 3D pore volume data, a pore network was extracted, and isosurfaces of the pore space were extracted to generate an initial surface mesh. S1 includes: S1-1: Micron-level industrial CT was used to perform separate operations on the target core. Two full-core scans were performed before and after the formation water-rock saturation reaction to obtain 3D CT image sequences. For the pre-reaction 3D CT image sequence, a cube-shaped computational subdomain free of edge artifacts was cropped from the original cylindrical scan image, and invalid black areas at the scan boundaries were removed. Specifically: Natural core plunger samples from the target gas reservoir, such as carbonate cores, with a diameter of 25 mm and a length of 50 mm, were selected. The samples were divided into two groups: one group served as a baseline sample before the reaction and underwent direct CT scanning; the other group was placed in a high-temperature, high-pressure reactor and introduced under target reservoir conditions, such as a temperature of 80°C and a pressure of 20 MPa. Saturated formation water is used to carry out water-rock reaction. The reaction time is set according to the reservoir development cycle, such as 30 days. After the reaction is completed, the core sample is taken out.

[0039] Micron-level X-ray micro-CT equipment was used to perform three-dimensional scanning of the core samples before and after the reaction. The scanning parameters were set as follows: tube voltage 140kV, tube current... The detector exposure time is 2 seconds, the scanning angle range is 0°-360°, the angle step is 0.25°, and the spatial resolution is set to... After scanning, a filtered back-projection algorithm is used for 3D reconstruction, resulting in two sets of 3D CT grayscale image sequences. The actual spatial size corresponding to a single voxel is [missing information]. .

[0040] For the pre-reaction 3D CT image sequence, an 800×800×800 voxel cube computational subdomain was cropped from the central region of the original cylindrical scan data, removing annular artifacts and invalid black boundary regions at the scan edges. Selecting the central subdomain ensures the statistical representativeness of the pore structure while avoiding interference from image distortion at the scan boundaries in subsequent calculations. The same cropping rules were used for subdomain extraction in the post-reaction 3D CT image sequence.

[0041] S1-2: The cropped 3D image sequence is subjected to median filtering for noise reduction using a 3×3×3 sliding window. Specifically: voxel coordinates Centered on the target voxel, take the x-axis offset. y-axis offset z-axis offset For all 27 neighboring voxels, extract the corresponding gray values ​​of these 27 voxels, and take the median of the gray values ​​as the gray value of the central voxel after filtering. The gray values ​​of the filtered image satisfy the following:

[0042] in, The original 3D CT image in voxel coordinates The grayscale value at that location; This represents the grayscale value at the corresponding location after noise reduction. The median value of the neighborhood grayscale values ​​is calculated. These represent the spatial coordinate translations of the sliding window along the x, y, and z axes, respectively. Three-dimensional median filtering can effectively eliminate salt-and-pepper noise and random noise introduced by CT scans, while preserving the edge features of the pore-rock skeleton interface and avoiding pore boundary blurring caused by filtering.

[0043] S1-3: Based on the globally optimal segmentation threshold t, perform binarization segmentation on the denoised image. Specifically: Traverse every candidate threshold within the entire grayscale value range of the image, divide the image voxels into rock skeleton class with grayscale greater than or equal to the candidate threshold and pore space class with grayscale less than the candidate threshold, calculate the inter-class variance of the two classes of voxels, and select the candidate threshold corresponding to the maximum value of the inter-class variance as the global optimal segmentation threshold t. Binarization of the image is performed based on the optimal segmentation threshold t to distinguish between the rock skeleton and pore space, resulting in binarized three-dimensional pore volume data, expressed as:

[0044] in, The voxel label values ​​of the binarized image. Indicates the rock skeleton region. This refers to the porous space region; After segmentation, the porosity of the binarized image is calculated and compared with the core porosity measured by the laboratory helium method. The relative error of segmentation is controlled within 5% to ensure segmentation accuracy. Small isolated noise voxels in the segmentation results are removed by morphological opening operation of 3 voxels, and finally binarized three-dimensional pore volume data is obtained.

[0045] S1-4: Based on the binarized 3D porous volume data, the pore network is extracted, and the isosurfaces of the pore network are extracted to generate an initial surface mesh; specifically: First, the watershed oversegmentation-label merging (SNOW) algorithm is used to extract the pore network from the binarized 3D pore volume data. First, the Euclidean distance from all pore voxels to the nearest rock skeleton voxel is calculated to generate a 3D distance map. Then, a watershed transformation is applied to the distance map to obtain oversegmented pore sub-regions. Using the local maxima of the distance map as pore center markers, adjacent oversegmented sub-regions are merged to eliminate false pores and oversegmentation errors. Finally, a topological pore network consisting of independent pore volumes and connecting throats is extracted, outputting the volume, equivalent sphere radius, and spatial coordinates of each pore, as well as the equivalent radius, length, and connected pore number of each throat.

[0046] Next, the Marching Cubes algorithm is used to extract the isosurfaces of the binarized pore space. The gray-level isosurfaces corresponding to the segmentation threshold t are used as the pore-skeleton interface. All voxel cube units are traversed, and the intersection points of the isosurfaces and voxel edges are calculated by interpolation to generate an initial surface mesh composed of triangular facet units. The coordinates of the mesh vertices correspond to the three-dimensional spatial position of the interface, and the normal vectors of the triangular units are uniformly pointed to the inside of the pores. Finally, the initial surface mesh file in STL format is output as the input for subsequent mesh topology repair and tetrahedral mesh construction.

[0047] The generated initial surface mesh is as follows Figure 2 As shown.

[0048] S2: Perform topology repair and mesh reconstruction operations on the initial surface mesh to construct a tetrahedral computational mesh; and use the tetrahedral computational mesh as the computational domain to construct a multiphysics coupled model that combines the Navier-Stokes flow equations with the Maxwell-Stefan multi-component convection-diffusion equations; the tetrahedral mesh processing flow is as follows. Figure 3 As shown, S2 includes: S2-1: Traverse all initial surface meshes generated by S1, identify defects including non-manifold edges and vertices, self-intersecting triangular faces, open boundary holes, and isolated porous surfaces not connected to the main channel, and mark the location and type of defect elements; where: Non-manifold edges and vertices: Count the number of triangle faces shared by each edge, and determine the edge that shares 3 or more triangle faces as a non-manifold edge; count the number of edge loops connected to each vertex, and determine the vertex that connects 2 or more independent edge loops as a non-manifold vertex; Triangle face self-intersection: By traversing all triangle faces using a spatial bounding box intersection detection algorithm, triangle faces that have spatial geometric intersections but do not belong to common edge / common vertex contact are judged as self-intersection defects; Open Boundary Holes: Extract all boundary edges that belong to only one triangle face, close them according to the connection relationship to form a boundary loop, and each independent boundary loop corresponds to one open hole; Isolated pore surfaces: Using the region growth method, starting from the pore grid at the core injection end and the extraction end, all connected surface units are traversed. Closed surface areas that are not traversed are identified as isolated pores that are not connected to the main channel.

[0049] S2-2: Delete all surface mesh elements corresponding to isolated pores that are not connected to the main channel between the injection and extraction ends, and seal the open boundary voids on the surface of the remaining effective pores; specifically: Isolated pore removal: Based on the isolated pore elements marked in the defect list, delete all triangular faces corresponding to isolated pores in batches to eliminate invalid calculation dead zones and retain only the effective pore surface mesh connected to the main channel of the injection end-extraction end. Open hole closure: For each open boundary hole on the effective pore surface, the point-by-point triangulation method is used for closure: three adjacent boundary vertices are selected sequentially along the hole boundary ring to generate triangular filling units, while ensuring that the normal vector of the newly added triangle is consistent with the normal vector of the original surface mesh; all boundary rings are traversed until all holes are closed, resulting in a closed surface mesh that initially meets the watertightness requirements.

[0050] S2-3: For triangular elements with self-intersecting triangular faces in the initial surface mesh, the face intersections are eliminated by fine-tuning the local vertex coordinates and reconnecting them to the triangular elements. The Laplacian smoothing algorithm is then used to optimize the vertex positions, eliminating step-like pixelation errors, resulting in a smoothed initial surface mesh. Specifically: Self-intersection defect repair: For the marked self-intersecting triangular faces, spatial intersection is eliminated by fine-tuning the local vertex coordinates, flipping the edges, and reconstructing the local faces; during the adjustment process, the maximum displacement of the vertices is controlled to not exceed 0.2 times the voxel size to avoid distortion of the pore geometry. After repair, there is no spatial self-intersection of any of the triangular faces. Laplacian Mesh Smoothing: The Laplacian smoothing algorithm is used to optimize the vertex positions of the surface mesh, eliminating the step-like boundary errors caused by pixelation in CT images. The smoothing iteration formula is as follows:

[0051] in, Let be the original coordinate vector of the i-th grid vertex. The smoothed vertex coordinate vector; Let i be the set of vertices directly adjacent to vertex i. The number of adjacent vertices; The smoothing relaxation factor is 0.2 in this embodiment.

[0052] In this embodiment, the smoothing iteration number is set to 5 times, which effectively smooths the pore wall surface and avoids excessive smoothing that could cause deviations in key geometric parameters such as pore radius and throat size, ultimately resulting in a smoothed watertight surface mesh.

[0053] S2-4: Using the smoothed surface mesh of the initial surface mesh as boundary constraints, a volume mesh is generated. This mesh is then imported into the finite element simulation software for local mesh refinement in the narrow channel region, resulting in a tetrahedral computational mesh for the entire porosity domain. Specifically: Using the smooth, closed surface mesh obtained from S2-3 as boundary constraints, an initial tetrahedral volume mesh was generated using the constrained Delaunay tetrahedral meshing algorithm. The generated volume mesh was then imported into the finite element numerical simulation software, where local mesh refinement was performed on the pore throat and narrow flow channel regions. The maximum side length of the elements in the throat region was set to not exceed [a certain value]. The maximum side length of the wide-pore region element does not exceed .

[0054] After the mesh is generated, a quality check is performed: all negative volume elements are removed, and the aspect ratio of all tetrahedral elements is controlled to be no greater than 5. Finally, a tetrahedral computational mesh covering the entire pore domain and meeting the numerical solution accuracy requirements is obtained, which serves as the computational domain for multiphysics coupling simulation.

[0055] S2-5: Using a tetrahedral computational mesh as the computational domain, a multiphysics coupled model is established in finite element simulation software, which combines the Navier-Stokes flow equations with the Maxwell-Stefan convection-diffusion equations. The Navier-Stokes flow equations include: The equation for the conservation of mass is expressed as follows:

[0056] The equation for the conservation of momentum is expressed as follows:

[0057] in, The density of the gas mixture, The velocity vector of the fluid. For fluid pressure, The dynamic viscosity of the gas mixture. For unit tensors, For volume forces, It is a vector differential operator; The Maxwell-Stefan convection-diffusion equations include: The equation for the conservation of mass of the components is expressed as follows:

[0058] The equation for the relative diffusion flux of the components is expressed as follows:

[0059] The diffusion driving force equation is expressed as follows:

[0060]

[0061] in, Let t be the mass fraction of component i, and t be the injection time. Source terms generated for components, Let be the relative diffusion flux of component i. The Maxwell-Stefan diffusion coefficient is... Where is the thermal diffusivity, and T is the system temperature. As a diffusion driving force, Let k be the mole fraction of component k. Total pressure; Let k be the molar mass of component k. The average molar mass of the gas mixture; Let be the molar mass of the i-th component; This embodiment is for Displacement A two-component system, components Corresponding to and Under isothermal simulation conditions, the thermal diffusion term can be turned off, and only molecular diffusion driven by concentration gradient and pressure gradient can be considered.

[0062] S2-6: Define the two end faces of the core along the flow direction as the injection end and the production end, including: Rock solid wall: no-slip boundary and zero diffusion flux boundary , where n is the out-of-plane normal unit vector; Injection end: constant velocity boundary and constant injection component boundary ; To maintain a constant inlet flow rate, the value used in this embodiment is... , The components are The mole fraction, i.e., the injected gas is pure ; Production end: Constant pressure boundary and zero diffusion flux boundary , This creates relative pressure on exports.

[0063] The initial conditions of the model are set as follows: the entire pore domain is filled at the initial moment. ,Right now The initial pressure across the entire domain is ; S3: Dynamically correct the viscosity and density of the mixed fluid in the multiphysics coupling model, and solve the model based on the corrected model. Obtain the mole fraction changes of different fluid components inside the target core at different injection time points, the convective and diffusion fluxes at different cross sections, and the thickness ratio of the mixing layer front. Set up a diffusion-free control condition, and calculate the proportion of diffusion mechanism influence by comparing the component mole fractions of the baseline condition and the diffusion-free condition. S3 includes: S3-1: Based on the fundamental viscosity calculated by Wilke's mixing rule, a polynomial correction factor with respect to the mole fraction of the components is introduced to dynamically correct the non-ideal mixing viscosity under high temperature and high pressure. Specifically, the calculation of atmospheric pressure is first based on Wilke's mixing rule. The basic viscosity of a binary gas mixture Then, a polynomial correction function related to the mole fraction of the components is introduced. The high-pressure correction for the base viscosity is expressed as follows:

[0064]

[0065] in, To correct the viscosity of the fluid mixture, The viscosity of the fluid mixture calculated using the Wilke formula. Let i, j, and k be the mole fractions of components i, j, and k, respectively. In this embodiment, ... Two-component system; The constant correction factor is obtained by fitting experimental data from a high-temperature and high-pressure falling ball viscometer under the temperature and pressure conditions of the target reservoir. During the numerical solution process, the mixing viscosity of the corresponding cell is dynamically calculated based on the real-time mole fraction of the components in each grid cell, thereby achieving dynamic correction of viscosity non-ideality throughout the displacement process.

[0066] S3-2: Construct the deviation factor as a function of temperature, pressure, and component mole fraction, and substitute it into the ideal gas law to complete the density correction of the mixed fluid; Specifically, the deviation factor Z is solved using the Peng-Robinson equation of state. The deviation factor is a multivariate function of system temperature T, total pressure p, and component mole fraction. Substituting the obtained deviation factor into the density calculation formula, the density of the mixed fluid under the corresponding temperature, pressure, and composition conditions is obtained, as expressed in the following expression:

[0067] in, For the density of the fluid mixture, For pressure, The average molar mass of the fluid mixture. The gas deviation factor is a function of temperature, pressure, and mole fraction of components. This is the universal gas constant, with values ​​ranging from 1 to 10. ; System temperature; In sync with viscosity correction, the density parameters are updated in real time at each time step based on the component distribution of each grid cell, matching the nonlinear evolution characteristics of the physical properties of multi-component mixing during the displacement process.

[0068] S3-3: Using the corrected viscosity and density of the mixed fluid from S3-1 and S3-2 as parameters, the multiphysics coupling model is input. A sequential coupling strategy of steady-state flow and transient mass transfer is employed. First, the steady-state flow field is solved. Then, based on the fixed flow field, the transient multi-component mass transfer process is solved, outputting the spatiotemporal distribution of components across the entire pore domain. Specifically: The steady-state Navier-Stokes flow equations are solved as follows: The initial viscosity and density parameters of the mixed fluid, modified in S3-1 and S3-2, are substituted into the steady-state Navier-Stokes governing equations constructed in S2-5. The velocity inlet, pressure outlet, and wall no-slip boundary conditions set in S2-6 are applied. Spatial discretization is performed using the finite element method, and the solution is iteratively obtained using the SIMPLE algorithm. A velocity field residual convergence threshold is set to [value missing]. The convergence threshold of the pressure field residual is After iterative convergence, the steady-state velocity field and pressure field of the entire pore domain are obtained, which serve as the fixed flow field for subsequent transient mass transfer calculations.

[0069] Transient Maxwell-Stefan mass transfer solution: Taking the steady-state velocity field as a fixed flow field, substituting it into the Maxwell-Stefan convection-diffusion equations, and setting the initial conditions: the entire pore domain is filled at the initial moment. ,Right now , The backward Euler method is used for time discretization, and the time step is set. The total simulation duration extends to the extraction end. After the breakthrough, the entire core The average mole fraction dropped below 5%.

[0070] During the solution process, after updating the spatial distribution of components at each time step, the viscosity and density of the corresponding grid cells are updated simultaneously by calling the property correction formulas of S3-1 and S3-2 to ensure dynamic matching between the property parameters and the component distribution. After the solution is completed, the spatial distribution data of the mass fraction and mole fraction of each component in the entire pore domain at all time steps are output.

[0071] S3-4: Based on the results output from S3-3, extract the characteristic cross-sectional data perpendicular to the flow direction, calculate the convective flux and diffusion flux at this cross-section at different injection times, and generate, based on the convective flux and diffusion flux at this cross-section at different injection times, as shown below. Figure 4 The flux-time evolution curve is shown below; Specifically, along the main flow direction of the core, i.e., the z-axis direction, a vertical section at the midpoint of the core's axial direction is selected as the characteristic section; for each injection time step, the values ​​at each point on the section are calculated. For flux density and diffusion flux density, for flux The calculation expression is:

[0072] diffusion flux The calculation expression is:

[0073] in, Let i be the convective flux density vector of component i. Let be the diffusion flux density vector of component i, and the other parameters are defined in the same way as the mass transfer equation parameters in S2-5; In the above steps, the component of the flux vector in the cross-section normal is taken, and the average convective flux density and average diffusion flux density of the cross-section are calculated. The calculation is repeated for all injection time steps to obtain the flux values ​​corresponding to the characteristic cross-section at different injection times. Evolution curves of convective and diffusion fluxes with injection time are plotted. For example... Figure 4 The flux-time evolution curve shown can be seen that, in During the initial injection phase (0-0.019 s), both the growth rate and absolute value of diffusion flux were significantly higher than those of convective flux, indicating that the transport of the mixed gas front was primarily driven by the concentration gradient, with interdiffusion playing a decisive role in the early breakthrough. As displacement progressed, the two curves intersected at 0.019 s, marking a critical shift in the mass transfer-dominant mechanism; thereafter, convective flux rapidly surpassed diffusion flux and became absolutely dominant, eventually stabilizing at approximately 0.44. The constant level indicates The influx of subjects achieved the Rapid overall displacement. It is noteworthy that, due to the rapid enhancement of convection which briefly exacerbated the local concentration gradient, the diffusion flux did not immediately decrease after 0.019 s, but instead reached approximately 0.13 around 0.02 s. The peak value, followed by the residual peak value. The gradual emptying and smooth decay towards zero.

[0074] S3-5: Based on the flux time evolution curve, identify the flux boundary point where the convective flux and diffusion flux are equal in value, and determine the position and length of the boundary interface along the z-axis corresponding to this boundary point. The spatial division is based on the interface, with the upstream region designated as the convection-dominant region and the downstream region as the diffusion-dominant region. The output zoning results are as follows: parameter; S3-6: Based on the defined convection-dominant and diffusion-dominant regions, the overall mixing layer is divided into a rear mixing layer and a leading mixing layer. The ratio of the pore volume at the leading mixing layer within the diffusion-dominant region to the total pore volume of the overall mixing layer is calculated to obtain the following... Figure 5 The percentage of the leading edge of the hybrid layer is shown; specifically: Total Pore Domain Output by S3-3 mole fraction as input, with The molar fraction of 5%-95% is used as the defining threshold to ensure that the pore space of the entire core meets the following requirements. All porous regions are defined as a monolithic mixed layer. Indicates injected components mole fraction; Based on the dominant region division results of S3-5, the overall mixing layer is divided into two parts: the mixing layer located in the convection-dominant region is the rear part of the mixing layer, and the mixing layer located in the diffusion-dominant region is the front edge of the mixing layer. Next, volume integration is performed over the entire pore domain to calculate the volume at the mixing layer front and the total volume of the overall mixing layer, as expressed by:

[0075]

[0076]

[0077] in, For the leading edge volume of the mixing layer, It is the length of the distance from the z-axis. Let i be the mole fraction of component i. The total volume of the mixed layer, The value represents the proportion of the leading edge of the mixed layer; the values ​​within square brackets are the region determination conditions, which are 1 when all conditions are met and 0 when none are met.

[0078] In the above steps, based on Figure 4 At the determined 0.019s alternation time point, this embodiment applies to the multi-component gas mixing layer ( The spatial distribution of the pore volume (with a mole fraction between 5% and 95%) was precisely quantified. The diffusion-dominated region in front of the cross-section was defined as the mixing layer front, and the convection-dominated region behind the cross-section was defined as the mixing layer rear. Volume integral calculations show that the mixing layer front accounts for 31.7% of the total pore volume of the mixing layer (e.g., ...). Figure 5 As shown in the three-dimensional microscopic distribution diagram, blue represents the leading edge of the mixing layer, red represents the rear part of the mixing layer, and yellow represents... The mole fraction is still higher than 95% in the unconducted region. This high percentage of 31.7% directly confirms that in complex micropore-throat structures (such as abrupt changes in pore size or blind-end pores), local velocity limitation greatly expands the diffusion-dominant region driven by concentration gradients, indicating that micro-diffusion is by no means a negligible local boundary effect.

[0079] S3-7: Based on the displacement simulation results in S3-3 as the baseline simulation results, a diffusion-free control condition was set up to obtain the whole core sample without considering diffusion. The curve of average mole fraction versus time; first, extract the results from the normal simulation and plot them. and The curve showing the change in average mole fraction over time, as shown below. Figure 6 As shown; next, a diffusion-free control condition is set up, the process of which is as follows: keep the model's geometric mesh, boundary conditions, and physical property parameters completely consistent with the baseline condition, and set the Maxwell-Stefan diffusion coefficient... Set to minimum value This makes the component diffusion flux approximately zero, completely eliminating the diffusion effect; the transient mass transfer solution process in S3-3 is repeated to obtain the following results. Figure 7 The whole core sample shown is without considering diffusion. The curve showing the change in average mole fraction over time.

[0080] Next, determine the time point of complete replacement: in Figure 6 In the baseline simulation results considering diffusion, the location of the whole core was determined. The point at which the average mole fraction drops to 5% is considered the point of complete displacement. ,satisfy .

[0081] Finally extract Figure 7 In the middle, the whole core residue under the non-diffusion control condition at the same time point Average mole fraction, calculated under the baseline simulation conditions and the diffusion-free control conditions. The absolute difference in mole fractions is used to obtain, for example... Figure 8 The percentage of influence of the diffusion mechanism on the overall displacement and evacuation process is shown below. This percentage is used to quantify the influence of the diffusion mechanism on the displacement and evacuation process, expressed as:

[0082]

[0083] in, For diffusion-free control conditions mole fraction To consider the diffusion scenario mole fraction for The time point when the mole fraction is 5%, The proportion of the influence of diffusion mechanism.

[0084] pass Figure 6 , Figure 7 and Figure 8 The results intuitively show that the injection rate has a significant positive control effect on the replacement efficiency; as the flow rate increases from 0.0006 m / s to 0.0018 m / s, The time to complete replacement when the mole fraction drops to 5% is significantly shortened. Meanwhile, a lateral comparison considering diffusion ( Figure 6 ) and exclusion diffusion effect ( Figure 7 The simulation results show that, at any constant flow velocity, after excluding diffusion... The attenuation curve became significantly slower, confirming that the concentration gradient-driven interdiffusion within the micropores can effectively accelerate the process. The stripping and removal of [materials] is an indispensable promoting mechanism in the extraction process. Further comprehensive quantitative results ( Figure 8 The study revealed a clear negative correlation between the diffusion effect and the injection rate. The diffusion effect gradually decreased from approximately 18% at low flow rates (0.0006 m / s) and stabilized at approximately 11% at high flow rates (0.0018 m / s). This phenomenon is attributed to the fact that high injection rates significantly enhance the convection dominance of the fluid and drastically compress the residence time of multi-component gases within the pore network, thus limiting the time window for interdiffusion. Conversely, at low flow rates, the longer residence time provides sufficient conditions for concentration gradient-driven diffusion to unfold, thereby significantly amplifying the mass transfer contribution of the diffusion mechanism in the mixing layer front and complex blind-end pores.

[0085] S4: For the acquired post-reaction 3D CT image sequences, steps S1 to S3 are processed independently, and the changes in quantitative characteristic parameters such as component mole fraction, convective flux, diffusion flux, mixing layer front thickness ratio, and diffusion mechanism influence ratio before and after the reaction are compared to quantify the changes. Water-rock reaction and diffusion effects on gas reservoir injection The degree of influence of the extraction mechanism; among which, S4 includes: S4-1: Obtain the value from S1 The three-dimensional CT image sequence following the water-rock reaction is used as input. The entire process from preprocessing (S1) to S3 is executed independently and repeatedly to obtain the complete set of results for the post-reaction core, including the component mole fraction field, convection and diffusion flux evolution curves, the proportion of the mixing layer front, and the proportion of diffusion mechanism influence. Specifically: Data collected using S1-1 The three-dimensional CT image sequence after the water-rock reaction is used as input, and the complete process of S1-1 image cropping, S1-2 median filtering and noise reduction, S1-3 threshold segmentation, S1-4 pore network extraction and initial surface mesh generation, S2-1 to S2-4 mesh topology repair and tetrahedral computational mesh construction, S2-5 to S2-6 multiphysics coupling model construction and boundary condition setting, and S3-1 to S3-7 dynamic property correction and transient solution and feature parameter extraction is executed independently and repeatedly.

[0086] During implementation, strict control was exercised over variable consistency: all image preprocessing parameters, mesh generation standards, physical property correction coefficients, boundary condition values, solver convergence thresholds and time step settings were completely consistent with the baseline simulation before the reaction. Only the pore geometry model was replaced with the pore structure of the core after the reaction, ensuring that the two sets of simulation results had unique variable comparability.

[0087] After the calculation is completed, the output will show the complete set of results, including the mole fraction field of the whole pore domain components, the evolution curves of convective and diffusion fluxes of the characteristic section, the proportion of the mixing layer front, and the proportion of the influence of diffusion mechanism, which will serve as the data basis for subsequent comparative analysis before and after the reaction.

[0088] S4-2: Select characteristic sections at the same injection time and spatial location, and compare the differences in three types of field parameters—component mole fraction, convection flux, and diffusion flux—of the core samples before and after the reaction, and calculate the relative change rate of the corresponding parameters; for example... Figure 9 As shown, it includes: Average relative change rate of mole fraction:

[0089] The calculation method for the average relative change rate of mole fraction is the same as described above. The calculation method for the average relative change rate of mole fraction is the same, and will not be elaborated further in this application.

[0090] The relative rate of change of convective flux at the characteristic cross section:

[0091] Relative rate of change of diffusion flux at characteristic cross section:

[0092] in, , These are the corresponding positions on the characteristic cross-sections before and after the reaction. mole fraction; , These are the average convective flux densities of the characteristic cross sections before and after the reaction at time n, respectively. , , , respectively, represent the average diffusion flux density of the characteristic cross section before and after the reaction at time n; A is the pore flow area of ​​the characteristic cross section; The total number of typical moments selected; S4-3: Extract the core quantitative characteristic parameters from the simulation results before and after the reaction, and calculate the changes in the proportion of the mixing layer front and the proportion of the diffusion mechanism's influence caused by the water-rock reaction; specifically including: Change in the proportion of the leading edge of the mixed layer :

[0093] The change in the proportion of diffusion mechanism influence :

[0094] in, , The percentage of the leading edge of the mixed layer before and after the reaction; , The percentage of diffusion mechanism influence before and after the reaction are respectively; S4-4: Combining the relative rate of change of corresponding production, the change in the proportion of the mixing layer front, the change in the proportion of diffusion mechanism influence, and the time required to complete replacement, the individual contribution weights of water-rock reaction to convection displacement and diffusion mass transfer through pore structure evolution are decomposed to obtain the coupling mechanism of water-rock reaction and diffusion effect, ultimately completing... Water-rock reaction and diffusion effects on gas reservoir injection Quantitative characterization of the impact of the extraction mechanism; diffusion effects before and after the reaction and the replacement completion time, etc. Figure 10 As shown.

[0095] Based on the above steps and Figure 9 , Figure 10 The results show that after the reaction... The time required for complete replacement to reduce the mole fraction to 5% was significantly prolonged compared to before the reaction (from approximately 1.25 s before the reaction to approximately 1.55 s after the reaction). This is because, although the prolonged physicochemical interactions increased the overall pore volume, local precipitation and uneven dissolution significantly exacerbated the heterogeneity and tortuosity of the pore network, leading to poorer effective connectivity for fluid transport and thus slowing down the overall replacement process. Simultaneously, by comparing data from control groups considering diffusion and excluding diffusion effects, and combining... Figure 10 The percentage histogram of diffusion effects shows that the absolute influence of the diffusion mechanism after the reaction increases compared to before the reaction (from approximately 13.5% to nearly 15%). This is because the water-rock reaction induces an increase in complex blind-end pores or micropore throats, which severely restricts convection in these local areas. Consequently, the gas reservoir relies more heavily on concentration gradient-driven interdiffusion to strip and extract deep residual gas in the later stages of displacement. .

[0096] The above are merely embodiments of the present invention. Commonly known structures and characteristics are not described in detail here. Those skilled in the art are aware of all common technical knowledge in the field prior to the application date or priority date, are aware of all existing technologies in that field, and have the ability to apply conventional experimental methods prior to that date. Those skilled in the art can, under the guidance of this application, improve and implement this solution in combination with their own capabilities. Some typical known structures or methods should not be obstacles for those skilled in the art to implement this application. It should be noted that those skilled in the art can make several modifications and improvements without departing from the structure of the present invention. These should also be considered within the scope of protection of the present invention, and will not affect the effectiveness of the implementation of the present invention or the practicality of the patent. The scope of protection claimed in this application should be determined by the content of its claims, and the specific embodiments described in the specification can be used to interpret the content of the claims.

Claims

1. Based on Gas reservoir injection based on water-rock reaction and diffusion effects The method for quantitative characterization of microscopic features is characterized by: include: S1: Collect target core samples Two sets of real three-dimensional CT image sequences were obtained before and after the water-rock reaction. The three-dimensional CT image sequence before the reaction was binarized and segmented to distinguish the rock skeleton and pore space, thus obtaining three-dimensional pore volume data. Based on the three-dimensional pore volume data, the pore network was extracted, and the pore space isosurface was extracted to generate an initial surface mesh. S2: Perform topology repair and mesh reconstruction operations on the initial surface mesh to construct a tetrahedral computational mesh; and use the tetrahedral computational mesh as the computational domain to construct a multiphysics coupling model that couples the Navier-Stokes flow equations with the Maxwell-Stefan multi-component convection-diffusion equations. S3: Dynamically correct the viscosity and density of the mixed fluid in the multiphysics coupling model, and solve the problem based on the corrected multiphysics coupling model. Obtain the changes in mole fraction of different fluid components inside the target core at different injection time points, the convection flux and diffusion flux at different cross sections, and the thickness ratio of the mixing layer front. Set up a no-diffusion control condition and calculate the proportion of diffusion mechanism influence by comparing the component mole fractions of the baseline condition and the no-diffusion condition. S4: For the acquired post-reaction 3D CT image sequences, steps S1 to S3 are processed independently, and the changes in quantitative characteristic parameters such as component mole fraction, convective flux, diffusion flux, mixing layer front thickness ratio, and diffusion mechanism influence ratio before and after the reaction are compared to quantify the changes. Water-rock reaction and diffusion effects on gas reservoir injection The extent to which the extraction mechanism is affected.

2. The method based on claim 1 Gas reservoir injection based on water-rock reaction and diffusion effects The method for quantitative characterization of microscopic features is characterized by: S1 includes: S1-1: Micron-level industrial CT was used to perform separate operations on the target core. Two whole-core scans were performed before and after the saturated formation water-rock reaction to obtain three-dimensional CT image sequences before and after the reaction. For the three-dimensional CT image sequence before the reaction, a cube-shaped computational subdomain without edge artifacts was cropped from the original cylindrical scan image, and invalid black areas at the scan boundary were removed. S1-2: The cropped 3D image sequence is subjected to median filtering noise reduction using a 3×3×3 sliding window. The grayscale values ​​of the filtered image satisfy the following: in, The original 3D CT image in voxel coordinates The grayscale value at that location; This represents the grayscale value at the corresponding location after noise reduction. The median value of the neighboring grayscale values ​​is calculated. These represent the spatial coordinate translations of the sliding window along the x, y, and z axes, respectively. S1-3: Based on the globally optimal segmentation threshold t, the denoised image is binarized to distinguish between the rock skeleton and the pore space, resulting in binarized three-dimensional pore volume data, expressed as: in, The voxel label values ​​of the binarized image. Indicates the rock skeleton region. This refers to the porous space region; S1-4: Based on the binarized 3D porous volume data, extract the pore network and extract the isosurfaces of the pore network to generate an initial surface mesh.

3. The method based on claim 2 Gas reservoir injection based on water-rock reaction and diffusion effects The method for quantitative characterization of microscopic features is characterized by: S2 includes: S2-1: Traverse all initial surface meshes generated by S1, identify defects including non-manifold edges and vertices, self-intersecting triangular faces, open boundary holes, and isolated porous surfaces not connected to the main channel, and mark the location and type of defect elements; S2-2: Delete all surface mesh elements corresponding to isolated pores that are not connected to the main channel of the injection end-extraction end, and seal the open boundary voids on the surface of the remaining effective pores. S2-3: For the triangular elements with self-intersecting triangular faces in the initial surface mesh, the face intersection is eliminated by fine-tuning the local vertex coordinates and reconnecting with the triangular elements. The vertex position is optimized by using the Laplacian smoothing algorithm to eliminate the step-like pixelation error and obtain the smoothed initial surface mesh. S2-4: Using the smoothed surface mesh of the initial surface mesh after smoothing as the boundary constraint, a volume mesh is generated. This mesh is then imported into the finite element simulation software to refine the mesh locally in the narrow channel region, resulting in a tetrahedral computational mesh for the entire pore domain.

4. The method based on claim 3 Gas reservoir injection based on water-rock reaction and diffusion effects The method for quantitative characterization of microscopic features is characterized by: S2 further includes: S2-5: Using a tetrahedral computational mesh as the computational domain, a multiphysics coupled model is established in finite element simulation software, combining the Navier-Stokes flow equations and the Maxwell-Stefan convection-diffusion equations. The Navier-Stokes flow equations include: The mass conservation equation is expressed as follows: The equation for the conservation of momentum is expressed as follows: in, The density of the gas mixture, The velocity vector of the fluid. For fluid pressure, The dynamic viscosity of the gas mixture. For unit tensors, For volume forces, It is a vector differential operator; The Maxwell-Stefan convection-diffusion equations include: The equation for the conservation of mass of the components is expressed as follows: The equation for the relative diffusion flux of the components is expressed as follows: The diffusion driving force equation is expressed as follows: in, Let t be the mass fraction of component i, and t be the injection time. Source terms generated for components, Let be the relative diffusion flux of component i. The Maxwell-Stefan diffusion coefficient is... Where is the thermal diffusivity, and T is the system temperature. As a diffusion driving force, Let k be the mole fraction of component k. Total pressure; Let k be the molar mass of component k. The average molar mass of the gas mixture; Let be the molar mass of the i-th component; S2-6: Define the two end faces of the core along the flow direction as the injection end and the production end, including: Rock solid wall: no-slip boundary and zero diffusion flux boundary , where n is the out-of-plane normal unit vector; Injection end: constant velocity boundary and constant injection component boundary ; For a constant inlet flow rate, The components are mole fraction; Production end: Constant pressure boundary and zero diffusion flux boundary , This creates relative pressure on exports.

5. The method based on claim 4 Gas reservoir injection based on water-rock reaction and diffusion effects The method for quantitative characterization of microscopic features is characterized by: S3 includes: S3-1: Based on the fundamental viscosity calculated using Wilke's mixing rule, a polynomial correction factor related to the mole fraction of components is introduced to dynamically correct the non-ideal mixing viscosity under high temperature and high pressure. The expression is as follows: in, To correct the viscosity of the fluid mixture, The viscosity of the fluid mixture calculated using the Wilke formula. These are the mole fractions of component i, component j, and component k, respectively. It is a constant correction factor; S3-2: The deviation factor is constructed as a function of temperature, pressure, and component mole fraction. Substituting this into the ideal gas law completes the density correction of the mixed fluid. The expression is: in, For the density of the fluid mixture, For pressure, The average molar mass of the fluid mixture. The gas deviation factor is a function of temperature, pressure, and mole fraction of components. This is the universal gas constant. This refers to the system temperature.

6. The method based on claim 5 Gas reservoir injection based on water-rock reaction and diffusion effects The method for quantitative characterization of microscopic features is characterized by: S3 further includes: S3-3: Using the viscosity and density of the mixed fluid as modified in S3-1 and S3-2 as parameters, the Navier-Stokes flow equation of the multiphysics coupling model is input for solution to obtain the steady-state velocity field and pressure field; then, using the steady-state velocity field as a fixed flow field, the Maxwell-Stefan convection-diffusion equation is substituted into the solution to output the spatial distribution of the mass fraction and mole fraction of each component in the whole pore domain at different injection times; S3-4: Based on the results output from S3-3, extract the characteristic cross-sectional data perpendicular to the flow direction, and calculate the convective flux and diffusion flux at this cross-section at different injection times. The calculation expression is: diffusion flux The calculation expression is: Based on the convective and diffusion fluxes at this cross section at different injection times, flux-time evolution curves are generated. S3-5: Based on the flux time evolution curve, identify the flux boundary point where the convective flux and diffusion flux are equal in value, and determine the position and length of the boundary interface along the z-axis corresponding to this boundary point. The spatial division is based on the interface, with the upstream region designated as the convection-dominant region and the downstream region as the diffusion-dominant region. The output zoning results are as follows: parameter; S3-6: Based on the defined convection-dominant and diffusion-dominant regions, the overall mixing layer is divided into a rear mixing layer and a leading mixing layer. The proportion of the leading mixing layer is obtained by calculating the ratio of the pore volume at the leading mixing layer within the diffusion-dominant region to the total pore volume of the overall mixing layer. The expression is: in, The volume of the mixed layer front is defined by taking the core injection end section as the origin of the z-axis and the flow direction as the positive z-axis direction. Let z be the z-axis coordinate of the interface where the flux and diffusion flux are equal. It is the length of the distance from the z-axis. Let i be the mole fraction of component i. The total volume of the hybrid layer, The percentage of the leading edge of the hybrid layer; S3-7: Based on the displacement simulation results in S3-3 as the baseline simulation results, a diffusion-free control condition was set up to obtain the whole core sample without considering diffusion. The curve showing the change of average mole fraction over time; and the location of the whole core sample in the baseline simulation results considering diffusion. The point at which the average mole fraction drops to 5% is considered the point of complete displacement. Residual core samples from the whole core at the same time point are extracted for the diffusion-free control. Average mole fraction, calculated under the baseline simulation conditions and the diffusion-free control conditions. The absolute difference in mole fractions is used to quantify the proportion of the diffusion mechanism's influence on the displacement and evacuation process. The expression is: in, For diffusion-free control conditions mole fraction To consider the diffusion scenario mole fraction for The time point when the mole fraction is 5%, The proportion of the influence of diffusion mechanism.

7. The method based on claim 6 Gas reservoir injection based on water-rock reaction and diffusion effects The method for quantitative characterization of microscopic features is characterized by: S4 includes: S4-1: Obtain the value from S1 The three-dimensional CT image sequence after the water-rock reaction is used as input, and the entire process from S1 to S3 is executed independently and repeatedly to obtain the complete set of results, including the mole fraction field of the components, the evolution curves of convection flux and diffusion flux, the proportion of the mixing layer front, and the proportion of diffusion mechanism influence, corresponding to the core after the reaction. S4-2: Select characteristic sections at the same injection time and spatial location, compare the differences in three types of field parameters—component mole fraction, convection flux, and diffusion flux—of the cores before and after the reaction, and calculate the relative change rate of the corresponding parameters; S4-3: Extract the core quantitative characteristic parameters from the simulation results before and after the reaction, and calculate the changes in the proportion of the mixing layer front and the proportion of the diffusion mechanism caused by the water-rock reaction. S4-4: By combining the relative rate of change of corresponding production, the change in the proportion of the mixing layer front, and the change in the proportion of diffusion mechanism influence, the individual contribution weights of water-rock reaction to convection displacement and diffusion mass transfer through pore structure evolution are decomposed, thus obtaining the coupling mechanism between water-rock reaction and diffusion effect, and finally completing... Water-rock reaction and diffusion effects on gas reservoir injection Quantitative characterization of the impact of the extraction mechanism.

Images