Numerical simulation methods, systems, media, and terminals for three-dimensional crack network dissolution processes

Abstract

This invention discloses a numerical simulation method, system, medium, and terminal for the dissolution process of a three-dimensional fracture network. The method includes: acquiring geological parameters of the target reservoir and constructing a geometric model of the three-dimensional fracture network; meshing the geometric model and embedding the fracture network into the bedrock mesh using an embedded discrete fracture method; determining the multi-field coupled control equations of solute transport, fluid flow, and water-rock reaction; performing numerical simulation of the fracture dissolution evolution process based on the control equations; and outputting the geometric morphology, solute concentration distribution, and conductivity changes of the fracture network after dissolution. This method maintains the true geometric morphology of the fractures, avoids the technical challenges of unstructured meshes, and balances computational accuracy and efficiency. It can effectively simulate water-rock reaction, solute transport, and convection-diffusion processes, and is suitable for predicting the dissolution process of a three-dimensional fracture network in soluble oil and gas reservoirs.

Description

Technical Field

[0001] This invention relates to the field of oil and gas reservoir development, and in particular to a numerical simulation method and system for the dissolution process of a three-dimensional fracture network. Background Technology

[0002] Fractured-vuggy carbonate reservoirs are typical high-quality reservoirs and a research hotspot in the field of oil and gas extraction. Numerical simulation of fractured-vuggy reservoirs, as a core method for quantitatively analyzing the movement of oil and water within the reservoir, provides important theoretical support for optimizing the extraction technology of fractured-vuggy reservoirs. However, existing technologies have the following shortcomings in the numerical simulation of dissolution in a full three-dimensional fracture network:

[0003] (1) The challenge of dynamic prediction of conductivity after dissolution of three-dimensional fracture network

[0004] The conductivity of the three-dimensional fracture network after dissolution is a core evaluation indicator of the effectiveness of acid fracturing, and existing models struggle to achieve accurate dynamic prediction of this conductivity. On one hand, dissolution causes non-uniform changes in the roughness and aperture distribution of the fracture walls, resulting in spatial differences in the conductivity of each fracture within the three-dimensional network. Existing models often extend the conductivity formula for a single fracture (such as the cubic law), neglecting the synergistic effect of network topology on conductivity. On the other hand, acid loss and reaction product precipitation during dissolution lead to partial fracture closure, making it impossible for existing models to accurately predict the attenuation of network conductivity during long-term production.

[0005] (2) The problem of efficient numerical solution and verification for three-dimensional crack network dissolution simulation

[0006] Simulating the dissolution of three-dimensional crack networks involves complex geometric topology, multi-field coupled equations, and dynamic evolution processes. Existing numerical methods face a trade-off between solution efficiency and simulation accuracy, and lack an effective experimental verification system. From a numerical solution perspective, to ensure simulation accuracy, fine-grid discretization of the crack network and dissolution region is required. However, in three-dimensional scenarios, fine-grid methods lead to an exponential increase in computational cost. Existing traditional numerical methods such as the finite volume method and the finite element method have extremely low solution efficiency, making it difficult to support long-term dissolution simulations of large-scale crack networks.

[0007] Therefore, there is an urgent need to develop an efficient numerical simulation method that can faithfully reproduce the geometry of cracks and accurately simulate the dissolution evolution process. Summary of the Invention

[0008] To address the shortcomings of existing technologies, this invention provides a numerical simulation method and system for the dissolution process of a three-dimensional crack network. The method aims to solve the technical problem of the lack of a quantifiable and reproducible dynamic characterization mechanism for the non-uniform dissolution process of a three-dimensional crack network in open or closed systems.

[0009] In a first aspect, the present invention provides a numerical simulation method for the dissolution process of a three-dimensional crack network, comprising:

[0010] S1: Obtain the geological parameters of the target reservoir and construct a geometric model of the three-dimensional fracture network;

[0011] S2: The geometric model is meshed, and the fracture network is embedded into the bedrock mesh using the embedded discrete fracture method;

[0012] S3: Determine the multi-field coupled control equations of solute transport-fluid flow-water-rock reaction; wherein, the multi-field coupled control equations of solute transport-fluid flow-water-rock reaction include the control equations of solute transport and water-rock reaction in the bedrock system and the control equations of solute transport and water-rock reaction in the fracture system;

[0013] S4: Numerical simulation of the crack dissolution evolution process based on the aforementioned governing equations;

[0014] S5: Outputs the geometric morphology, solute concentration distribution, and changes in conductivity of the fracture network after dissolution.

[0015] Furthermore, the specific process of S1 is as follows:

[0016] S11: Based on the seismic interpretation, well logging data and core observation results of the target reservoir, statistically analyze the probability distribution characteristics of key geometric parameters of fractures and establish the probability density function of each geometric parameter; wherein, the key geometric parameters of fractures include strike, dip angle, length, aperture and spatial density, etc.

[0017] S12: Based on the probability density function of each geometric parameter, a Monte Carlo simulation method is used to generate three-dimensional crack elements in batches to obtain cracks, and generation constraints are set; wherein, each crack is generated by parameters such as spatial starting point coordinates, direction, dip angle, length and aperture; the constraints include minimum crack spacing, size of the simulation region boundary where the crack is located, crack length range, cracks cannot completely overlap, crack aperture range and spatial distribution density.

[0018] S13: A spatial geometric intersection detection algorithm is used to verify the topological relationship between the generated cracks, remove cracks that do not meet the constraints, split the nodes of the intersecting cracks, clarify the spatial coordinates and connectivity of the crack intersection points, and obtain the crack database. The topological relationship between cracks includes intersection, inclusion, penetration and other cases. The data in the crack database includes crack number, spatial coordinates, geometric dimensions, strike and dip angle, and intersection node information.

[0019] S14: Based on core experiments, well logging interpretation, and field data, set the physical and mechanical parameters for each fracture; the physical and mechanical parameters include initial permeability, wall roughness, and mineral reactivity;

[0020] S15: Integrate the information from S11 to S14 to obtain the geometric model of the three-dimensional crack network.

[0021] Furthermore, the specific process of S2 is as follows:

[0022] S21: The background matrix is ​​divided using a structured mesh, the mesh size is set, each bedrock mesh cell is numbered and its parameters are stored. The stored parameters include physical and mechanical parameters such as coordinate range, volume, matrix permeability, porosity, and elastic modulus. In practice, it also includes adjusting meshes with excessive distortion and supplementing meshes with missing boundaries to ensure that the quality of the bedrock mesh meets the requirements of numerical solution.

[0023] S22: Determine the intersection relationship between the crack and the bedrock grid unit, and divide the crack unit;

[0024] S23: Construct the topological association between bedrock mesh cells and fracture cells, and store geometric and physical parameters to obtain a fracture network embedded in the bedrock mesh.

[0025] Furthermore, the specific process of S22 is as follows:

[0026] S221: For each crack, traverse all bedrock grid cells, calculate the intersection area between the crack surface and the grid cell, and if the area of ​​the intersection area is greater than zero, then the grid cell is determined to be a crack-containing grid cell, and the correspondence between the crack-containing grid cell number and the crack number is recorded.

[0027] S222: Based on the intersection of the crack and the bedrock grid, the crack is divided into several crack units, and each crack unit corresponds to a crack-containing grid unit.

[0028] S223: The matrix mesh portion embedded with the crack is mapped to a weighted graph using a continuous projection algorithm. The Dijkstra algorithm is then applied to find the shortest path between the two endpoints of the crack. This path is mapped back to the mesh to determine the surface receiving the projection, thereby completing the crack element division.

[0029] Furthermore, the solute transport and water-rock reaction governing equations of the bedrock system in the multi-field coupled governing equations of solute transport-fluid flow-water-rock reaction shown in S3 are as follows:

[0030] ;

[0031] in, Matrix porosity (dimensionless); This represents the concentration of the solute in the matrix. For simulating time; It is the divergence operator; The seepage velocity in the matrix; The dispersion coefficient of the solute in the matrix. ,in The molecular diffusion coefficient is... For longitudinal dispersion, I represents the lateral dispersion, and I represents the unit tensor. This refers to the dry density of the matrix. This represents the amount of solute adsorbed in the matrix. ,in This represents the maximum adsorption capacity. It is the adsorption equilibrium constant; For the source and sink terms of solutes in the matrix; This is the water-rock reaction rate term for the solute in the matrix. ,in The reaction rate constant is... The surface area of ​​the bedrock reaction. This represents the equilibrium concentration of the solute. This is the exchange term for solute between the matrix and the crack. ,in is the solute exchange coefficient.

[0032] Furthermore, the solute transport and water-rock reaction governing equations of the fracture system in the multi-field coupled solute transport-fluid flow-water-rock reaction equations shown in S3 are as follows:

[0033] ;

[0034] in, The crack opening; The concentration of the solute in the crack; The seepage velocity in the crack; The dispersion coefficient of the solute in the crack; The dry density of the crack wall; This represents the amount of solute adsorbed on the crack wall surface; These are the source and sink terms of solutes in the cracks; This is the water-rock reaction rate term for solutes in the fracture; This refers to the exchange of solute between connected fracture units. , The solute exchange coefficient between the cracks. For adjacent connected cracks Solute concentration in each unit.

[0035] Furthermore, the specific process of S4 is as follows:

[0036] S41: Initialize simulation parameters, set the total simulation duration and initial time step, adopt implicit time discretization scheme, and initialize the solute concentration field, pressure field, velocity field and aperture field of each grid element in the matrix system and fracture system;

[0037] S42: The solute transport and water-rock reaction coupling equation is split into convection operator, dispersion operator and adsorption-reaction operator by operator splitting method, and spatial discretization is performed by integral finite difference method respectively. The solute concentration field is updated by iterative coupling.

[0038] S43: Constructing a system of linear algebraic equations based on discrete results The coefficient matrix Based on the geometric parameters and coupling conditions of the matrix mesh and crack elements, the solution vector x includes fluid pressure and solute concentration. The right-hand term vector b is constructed based on boundary conditions and source-sink terms. The sparse matrix is ​​solved using the generalized minimum residual method to obtain the pressure field and concentration field at the current time step.

[0039] S44: Calculate the water-rock reaction rate of solute in each matrix grid and fracture unit based on the concentration field, update the bedrock porosity, fracture aperture and permeability according to the reaction rate, and feed the updated physical property parameters back to the governing equations of the next time step;

[0040] S45: Determine whether the solution at the current time step has converged. If it has converged, proceed to the next time step. If it has not converged, reduce the time step size and recalculate S42 to S44. Repeat the iteration until the preset total simulation time is reached.

[0041] Secondly, the present invention provides a numerical simulation system for a three-dimensional crack network dissolution process, the system being used to perform the steps of the method described above, including:

[0042] Geometric model construction module: Obtains geological parameters of the target reservoir and constructs a geometric model of the three-dimensional fracture network;

[0043] Mesh generation and embedding module: The geometric model is meshed, and the fracture network is embedded into the bedrock mesh using the embedded discrete fracture method;

[0044] The governing equation determination module determines the multi-field coupled governing equations of solute transport, fluid flow, and water-rock reaction; wherein, the multi-field coupled governing equations of solute transport, fluid flow, and water-rock reaction include the governing equations of solute transport and water-rock reaction in the bedrock system and the governing equations of solute transport and water-rock reaction in the fracture system;

[0045] Numerical solution module: Performs numerical simulation of the crack dissolution evolution process based on the aforementioned governing equations;

[0046] Results output module: Outputs the geometry of the fracture network after dissolution, the distribution of solute concentration, and the changes in conductivity.

[0047] Thirdly, the present invention provides a readable storage medium storing a computer program that, when invoked by a processor, performs the steps of the method described above.

[0048] Fourthly, the present invention provides an electronic terminal comprising a processor and a memory, the memory storing a computer program, the processor invoking the computer program to perform the steps of the method described above.

[0049] This invention proposes a numerical simulation method, system, medium, and terminal for the dissolution process of a three-dimensional fracture network. The method integrates three-dimensional fracture topology generation based on geostatistical parameters with the embedded discrete fracture method (EDFM) to achieve high-fidelity reconstruction of fracture geometry and spatial distribution, avoiding the technical challenges of unstructured meshes and balancing computational accuracy and efficiency. It constructs a multi-field coupled control equation set of solute transport, fluid flow, and water-rock reaction, and establishes a non-uniform dissolution prediction model driven by mineral composition, fluid chemical potential gradient, and local velocity field. The generalized minimum residual method (GMRES) is used to achieve efficient solution of the model, making up for the time-consuming and labor-intensive shortcomings of three-dimensional physical simulation experiments. This enables accurate and efficient simulation of the dissolution process of a three-dimensional fracture network, providing reliable technical support for related engineering design and optimization. Attached Figure Description

[0050] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0051] Figure 1 This is a flowchart illustrating a numerical simulation method for a three-dimensional crack network dissolution process according to an embodiment of the present invention.

[0052] Figure 2 This is a flowchart of the three-dimensional crack network generation and parameter characterization provided in the embodiments of the present invention;

[0053] Figure 3 This is a schematic diagram illustrating the construction principle of the embedded discrete crack model (EDFM) provided in this embodiment of the invention;

[0054] Figure 4 This is a schematic diagram of the multi-order water-rock reaction kinetic model provided in an embodiment of the present invention; wherein, Figure 4 (a) shows the dissolution rate of CaCO3 as a function of Ca. 2+ Schematic diagram of the changes; Figure 4 (b) represents the calcium ion equilibrium concentration C. eq Schematic diagram showing the change with temperature; Figure 4(c) is a schematic diagram showing the change of CO2 partial pressure with temperature;

[0055] Figure 5 This is a simulation diagram of the dissolution effect of a three-dimensional crack network provided in an embodiment of the present invention; Figure 5 (a) is a model diagram of the coupled mesh generation of cracks and matrix; Figure 5 (b) is a graph showing the change in the width of the dissolution fracture in the 60th year; Figure 5 (c) is a graph showing the calcium ion concentration distribution in the 60th year; Figure 5 (d) is a chemical reaction rate diagram; Figure 5 (e) is a graph showing the change in the width of the dissolution fracture in the 240th year; Figure 5 (f) is a graph showing the calcium ion concentration distribution in the 240th year.

[0056] Figure 6 This is a schematic diagram of the dissolution morphology of a salt rock fracture surface provided in an embodiment of the present invention;

[0057] Figure 7 This is a schematic diagram comparing simulation results with physical experimental results provided in the embodiments of the present invention;

[0058] Figure 8 This is a numerical model for the dissolution of a three-dimensional crack network in a soluble carbonate layer, provided in an embodiment of the present invention.

[0059] Figure 9 This image shows the calculation results of the three-dimensional fracture network dissolution of the carbonate layer. Figure 9 (a) The change of calcium ion flow rate over time at different replenishment sites; Figure 9 (b) shows the variation of the width of the dissolution fracture with the dissolution distance; Figure 9 (c) The calcium ion concentration distribution in the dissolution fracture network 100 years later; Figure 9 (d) shows the distribution of the fracture width of the dissolution fracture network after 100 years. Detailed Implementation

[0060] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be described in detail below. Obviously, the described embodiments are merely some embodiments of this invention, and not all embodiments. Based on the embodiments of this invention, all other implementation methods obtained by those skilled in the art without creative effort are within the scope of protection of this invention.

[0061] Example 1

[0062] like Figure 1 As shown in the figure, this embodiment provides a numerical simulation method for the dissolution process of a three-dimensional crack network, the method comprising:

[0063] S1: As Figure 2As shown, the geological parameters of the target reservoir are obtained, and a geometric model of the three-dimensional fracture network is constructed. Specifically, the construction of the geometric model of the three-dimensional fracture network includes the following steps:

[0064] S11. Extract statistical patterns of cracks and construct a probabilistic model:

[0065] Based on seismic interpretation, well logging data, and core observation results of the target reservoir, the probability distribution characteristics of key geometric parameters of fractures are statistically analyzed, and probability density functions for each geometric parameter are established. These key geometric parameters include strike (e.g., normal / uniform distribution), dip angle (0-90° range), length (log-normal distribution), aperture (power-law distribution), and spatial density (Poisson distribution). This ensures that the generated fracture network conforms to actual geological patterns.

[0066] S12. Monte Carlo Random Crack Generation and Constraint Control:

[0067] Based on the probability density function of each geometric parameter, a Monte Carlo simulation method is used to generate three-dimensional crack elements in batches to obtain cracks, and generation constraints are set. Each crack is generated by parameters such as spatial starting point coordinates, direction, dip angle, length, and aperture. The constraints include minimum crack spacing (≥0.3m), size of the simulation region boundary where the crack is located (specifically, the generated crack does not exceed the simulation region boundary), crack length range (0.5-50m), avoiding complete overlap with existing cracks, crack aperture range (10-1000μm), and spatial distribution density (0.1-5 cracks / m³).

[0068] S13. Crack topology verification and validity screening:

[0069] A spatial geometric intersection detection algorithm is used to verify the topological relationships between generated cracks, including intersections, inclusions, and penetrations. Cracks that penetrate each other, exceed boundaries, or whose spacing does not meet constraints are removed, while valid independent crack units are retained. Intersecting cracks are split into nodes to clarify the spatial coordinates and connectivity of the crack intersection points, forming a complete crack database containing crack numbers, spatial coordinates, geometric dimensions, orientation and dip angles, and intersection node information.

[0070] S14. Personalized assignment of crack heterogeneity parameters:

[0071] Based on core experiments, well logging interpretation, and field data, personalized physical and mechanical parameters are assigned to each effective fracture. Among them, the initial permeability is calculated based on the fracture aperture and wall roughness (corrected by the cubic law), the wall roughness is characterized by the fractal dimension (1.0-1.5), and the mineral reactivity is determined by the mineral composition of the fracture wall (the ratio of limestone / dolomite) to determine the reaction rate constant, which fully reflects the heterogeneous characteristics of the fracture network.

[0072] S15: Visual verification and optimization of the crack network:

[0073] The generated three-dimensional random fracture network is imported into a visualization platform to intuitively display the spatial distribution, connectivity, and positional relationship of the fractures with the matrix model. Based on the experience of geological experts, local adjustments are made to areas where the fracture density is too high or too low (such as adding or deleting some fractures) to ensure that the fracture network conforms to statistical laws and reflects the particularity of local fracture development in the reservoir. Finally, the information from S11 and S14 is integrated to obtain the geometric model of the three-dimensional fracture network.

[0074] S2: The geometric model is meshed, and the fracture network is embedded into the bedrock mesh using the embedded discrete fracture method. Specifically, as follows... Figure 3 As shown, a bedrock mesh is constructed, and the geometric model of the established three-dimensional fracture network is embedded into the bedrock mesh to complete the fracture element division and geometric information association, laying the foundation for the subsequent discretization and solution of the coupled equations. This mainly includes the following three steps:

[0075] S21. Mesh type selection and parameter settings:

[0076] In this embodiment, based on the simulation scenario of oil and gas reservoir seepage, a Cartesian structured grid (not limited to it, other structured grid methods, such as corner grids, can also be used) is selected as the background matrix partitioning method; the grid size is set according to the simulation accuracy and computational efficiency. To ensure that the grid size is less than a preset proportion of the minimum effective crack length (in this embodiment, it is set to 1 / 3 of the minimum effective crack length), and that the number of grid cells is controlled within a reasonable range, each grid cell is then numbered and its parameters are stored, including coordinate range, volume, matrix permeability, porosity, elastic modulus, and other physical and mechanical parameters.

[0077] Preferably, it also includes: adjusting the mesh size or mesh generation method for meshes with excessive distortion; supplementing mesh cells for meshes with missing boundaries to ensure that the quality of the bedrock mesh meets the requirements of numerical solution.

[0078] S22. 3D fracture network embedding: A spatial geometric intersection algorithm is used to determine whether the fracture surface of each fracture in the fracture network intersects with the bedrock mesh unit, and then the fracture unit is divided.

[0079] S221: For each crack, traverse all bedrock mesh elements, calculate the intersection area between the crack surface and the mesh element. If the area of ​​the intersection area is greater than zero, the mesh element is determined to be a cracked mesh element, and the correspondence between the cracked mesh element number and the crack number is recorded.

[0080] S222: For each fracture, based on its intersection with the bedrock grid, the fracture is divided into several fracture units, and each fracture unit corresponds to a fracture-containing grid unit.

[0081] S223: A continuous projection algorithm is used to ensure the continuity of crack element projection and avoid the accuracy reduction problem caused by discontinuous projection. That is, the matrix mesh part embedded in the crack is mapped as a weighted graph (nodes correspond to mesh nodes, and edges correspond to mesh surfaces). Dijkstra's algorithm is applied to find the shortest path between the two endpoints of the crack. This path is mapped back to the mesh to determine the surface that receives the projection, thereby completing the crack element division.

[0082] S23. Crack geometric information association and storage:

[0083] Construct topological relationships between bedrock mesh cells and fracture cells, including fracture-to-fracture relationships. ), cracks and matrix ( There are three types of non-adjacent connections between cracks, including crack element connections across cracks. ) and the connection between crack elements in the same crack ( Based on the fracture network of S22, the correspondence between intersecting fracture elements is determined. Then, the bedrock mesh data and fracture element data are stored in the integrated platform's database using an efficient data format to ensure rapid data retrieval and provide data support for subsequent coupled equation construction and numerical solutions. Finally, the core parameters of each fracture element are calculated, including the fracture element's length, width, area of ​​the intersecting region, aperture, and permeability; the number and core parameters of each matrix and fracture element are stored.

[0084] S3: Determine the multi-field coupled control equations for solute transport, fluid flow, and water-rock reaction.

[0085] To accurately characterize the dissolution evolution of three-dimensional fracture networks and overcome the limitation of traditional single-reaction kinetic models in reflecting the stage differences in the dissolution process, thereby improving the accuracy and practicality of fractured reservoir dissolution simulation, this invention introduces a multi-order chemical reaction kinetic model to describe the fracture network dissolution process, such as... Figure 4 As shown, the core is based on the dynamic evolution of calcium ion concentration during reservoir dissolution, constructing a functional relationship between chemical reaction rate and calcium ion concentration that is staged and continuously correlated. At the same time, combined with the connectivity characteristics between the reservoir and the atmosphere, the reservoir dissolution system is clearly divided into two categories: open system and closed system, so as to achieve accurate simulation of the fracture network dissolution process under different geological conditions.

[0086] Based on the connectivity between the reservoir and the atmosphere, reservoir dissolution systems are further divided into two categories: open systems and closed systems. Open systems are defined as reservoir dissolution scenarios with good connectivity to the atmosphere. These reservoirs have continuous atmospheric recharge channels (such as faults connecting to the surface or high-permeability fracture networks connecting to the atmosphere). During dissolution, the oxygen and carbon dioxide content of the fluid remains stable, and the calcium ions generated during dissolution can be continuously discharged from the reservoir through fluid flow, preventing calcium ion concentration saturation. The dissolution process mainly goes through an initial dissolution stage, a rapid dissolution stage, and finally stabilizes in a steady dissolution stage. The multi-stage kinetic model mainly focuses on parameter calibration for the first three stages and needs to introduce a fluid convection term to correct for the migration and discharge effect of calcium ions, ensuring that the coupling relationship between the reaction rate and calcium ion concentration conforms to the dynamic characteristics of an open system.

[0087] The closed system is defined as a reservoir dissolution scenario with poor atmospheric connectivity. Such reservoirs have no obvious atmospheric recharge channels, and the fluid is in a relatively closed environment during the dissolution process. The carbon dioxide content is gradually consumed, and the calcium ions generated by dissolution cannot be effectively discharged and can only accumulate in the fluid inside the reservoir, easily reaching a saturation state. The dissolution process goes through four stages: initial, rapid, stable, and decay. The multi-order kinetic model needs to fully calibrate the parameters of the four stages, focusing on optimizing the exponential decay parameters of the decay stage. At the same time, an ion diffusion term is introduced to correct the internal accumulation effect of calcium ions. Combined with the connectivity of the fracture network, the cumulative differences in calcium ion concentration and the spatial heterogeneity of dissolution rate in different fracture units are characterized.

[0088] Combining the solute transport mechanisms (convection, dispersion, adsorption-desorption) and water-rock reaction characteristics (dissolution-precipitation reaction) of fractured media, governing equations for solute transport and water-rock reaction are constructed for both bedrock and fractured systems. Fracture-matrix solute exchange terms and reaction coupling terms are introduced to couple solute transport and water-rock reaction. Solute transport in bedrock is mainly controlled by convection, molecular diffusion, and mechanical dispersion. Simultaneously, adsorption-desorption of solutes with the bedrock medium and water-rock dissolution-precipitation reactions occur. These reaction processes affect bedrock porosity and permeability, which in turn influence the seepage field and solute transport. Therefore, the governing equations for solute transport and water-rock reaction in the bedrock system are as follows:

[0089] ;

[0090] in, The matrix porosity (dimensionless) is dynamically evolved by the water-rock reaction. The initial value is determined based on field measurement data and is updated during the reaction process according to the amount of dissolution and precipitation. This represents the concentration of the solute in the matrix (mol / m³). The simulation time is in seconds. It is the divergence operator; The seepage velocity in the matrix is ​​(m / s). denoted as the dispersion coefficient (m² / s) of the solute in the matrix. ,in The molecular diffusion coefficient is... For longitudinal dispersion, I represents the lateral dispersion, and I represents the unit tensor. The dry density of the matrix is ​​(kg / m³). The amount of solute adsorbed in the matrix (mol / kg), which follows the Langmuir adsorption isotherm. , This represents the maximum adsorption capacity. It is the adsorption equilibrium constant; For the source and sink terms of solutes in the matrix (mol / (m³·s)); The water-rock reaction rate term for the solute in the matrix is ​​(mol / (m³·s)), with the dissolution reaction being positive and the precipitation reaction being negative, according to the reaction kinetic equation. calculate, The reaction rate constant is... The surface area of ​​the bedrock reaction. This represents the equilibrium concentration of the solute. This is the exchange term for solutes between the matrix and fractures (mol / (m³·s)). A positive value indicates the migration of solutes from the fractures to the bedrock, and a negative value indicates the migration from the bedrock to the fractures. calculate, It is the solute exchange coefficient, which is related to the crack-matrix contact area and seepage velocity.

[0091] Solute transport in fractures is primarily convection-driven, with relatively weak dispersion. Simultaneously, solute adsorption-desorption and water-rock dissolution-precipitation reactions occur (at higher rates than in bedrock due to the larger fracture surface area and stronger fluid flow). Fracture aperture is influenced by water-rock reactions and mechanical forces, which in turn affect seepage velocity and solute transport. Therefore, the governing equations for solute transport and water-rock reactions in the fracture system are as follows:

[0092] ;

[0093] in, Crack aperture (m); The concentration of the solute in the crack is (mol / m³). The seepage velocity in the crack is (m / s). The dispersion coefficient (m² / s) of the solute in the fracture is calculated in the same way as that of the bedrock. Since convection dominates in the fracture, the dispersion value is less than that of the bedrock. The dry density of the crack wall (kg / m³). The adsorption capacity (mol / kg) of solute on the fracture wall is given. The adsorption model is consistent with that of the bedrock, and the adsorption parameters are adjusted according to the mineral composition of the fracture wall. This represents the source and sink terms of solutes in the fracture (mol / (m³·s)), corresponding to the source and sink terms of the bedrock; The water-rock reaction rate term for solutes in the fracture (mol / (m³·s)) is calculated in the same way as for the bedrock, and the reaction rate constant is used. Larger than bedrock, reaction surface area The surface area per unit volume of the crack wall; The exchange term of solute between connected fracture units (mol / (m³·s)) is given by calculate, The solute exchange coefficient between the cracks. For adjacent connected cracks Solute concentration in each unit.

[0094] S4: Numerical simulation of the crack dissolution evolution process based on the aforementioned governing equations.

[0095] In practical implementation, considering the strong coupling characteristics of solute transport and water-rock reaction, a combination of operator splitting and integral difference methods is used to solve the coupled equations. The coupled equations are split into convection operators, dispersion operators, and adsorption-reaction operators, solved separately, and then coupled iteratively. This avoids the problems of singular coefficient matrices and convergence difficulties caused by direct coupling, while ensuring solution accuracy. The integral difference method is used for discretization of the convection-dispersion equations, overcoming the shortcomings of traditional finite difference methods in complex meshes (EDFM matrix mesh + fracture elements) where the dispersion term discretization accuracy is low and the convection term is prone to numerical dispersion. It balances local conservation and numerical stability, and is suitable for the complex geometry of three-dimensional random fracture networks. Based on the above methods, the partial differential equations are transformed into a system of linear algebraic equations. Where A is the coefficient matrix (constructed based on the parameters and coupling conditions of the matrix mesh and fracture elements), x is the solution vector (including fluid pressure and solute concentration), and b is the right-hand side vector (constructed based on boundary conditions and source / sink terms). Based on the characteristics of the coefficient matrix (sparse matrix), the Generalized Minimal Residual Method (GMRES) is used to efficiently solve the sparse matrix. An implicit time discretization scheme is employed, with a set time step, to progressively solve the linear algebraic equations at each time step, ensuring the stability of the calculation results. Therefore, the specific process of numerical simulation of the fracture dissolution evolution is as follows:

[0096] S41: Initialize simulation parameters, set the total simulation duration and initial time step, adopt implicit time discretization scheme, and initialize the solute concentration field, pressure field, velocity field and aperture field of each grid element in the matrix system and fracture system;

[0097] S42: The solute transport and water-rock reaction coupling equation is split into convection operator, dispersion operator and adsorption-reaction operator by operator splitting method, and spatial discretization is performed by integral finite difference method respectively. The solute concentration field is updated by iterative coupling.

[0098] S43: Constructing a system of linear algebraic equations based on discrete results The coefficient matrix Based on the geometric parameters and coupling conditions of the matrix mesh and crack elements, the solution vector x includes fluid pressure and solute concentration. The right-hand term vector b is constructed based on boundary conditions and source-sink terms. The sparse matrix is ​​solved using the generalized minimum residual method to obtain the pressure field and concentration field at the current time step.

[0099] S44: Calculate the water-rock reaction rate of solute in each matrix grid and fracture unit based on the concentration field, update the bedrock porosity, fracture aperture and permeability according to the reaction rate, and feed the updated physical property parameters back to the governing equations of the next time step;

[0100] S45: Determine whether the solution at the current time step has converged. If it has converged, proceed to the next time step. If it has not converged, reduce the time step size and recalculate S42 to S44. Repeat the iteration until the preset total simulation time is reached.

[0101] S5: Outputs the geometric morphology, solute concentration distribution, and changes in conductivity of the fracture network after dissolution.

[0102] Numerical solutions for the multi-field coupled model of solute transport, fluid flow, and water-rock reaction can be obtained through S3 and S4, including parameters such as fluid pressure, solute concentration, dissolution fracture width, and fracture conductivity coefficient. By outputting these parameters at each fracture grid and matrix grid in the binary format required by the TECPLOT software, a three-dimensional display of the fracture network geometry, solute concentration distribution, and conductivity variation can be obtained.

[0103] To verify the method, this implementation constructed a model with dimensions of 200m × 200m × 50m, embedding 100 cracks. Using the embedded discrete crack method, the matrix portion of the model was discretized into a 10m × 10m × 5m mesh, as shown below. Figure 5 As shown in (a). The bottom of the model is assumed to be an impermeable boundary, and the water supply is all from atmospheric precipitation, using a constant flow recharge method. The model is an open system. The main parameters of the model are set as follows: CO2 equilibrium concentration of 3.3 mol / m³; CO2 partial pressure of 0.5 MPa; initial temperature of 10℃; initial fracture width of 0.2 mm; precipitation recharge of 600 mm / year; initial fracture permeability of 50 mD; initial matrix permeability of 0.1 mD; porosity of 0.24; and local mass exchange coefficient of 2.0 × 10⁻⁴. -3m / s; molecular diffusion coefficient is 5.0 × 10⁻⁶ m / s; -9 m 2 / s; fluid viscosity is 1.0 cP; fluid density is 1000 km / m 3 Initial specific surface area 8000 m² 2 / m 3 Chemical reaction rate diagram as shown below Figure 5 As shown in (d), the calcium ion concentration distribution and dissolution fracture width obtained from simulation calculations after 60 years and 24 years are as follows. Figure 5 (b) Figure 5 (c) Figure 5 (e) Figure 5 As shown in (f).

[0104] Example 2

[0105] This embodiment provides a numerical simulation system for a three-dimensional crack network dissolution process. The system is used to perform the steps of the method described above, including:

[0106] Geometric model construction module: Obtains geological parameters of the target reservoir and constructs a geometric model of the three-dimensional fracture network;

[0107] Mesh generation and embedding module: The geometric model is meshed, and the fracture network is embedded into the bedrock mesh using the embedded discrete fracture method;

[0108] The governing equation determination module determines the multi-field coupled governing equations of solute transport, fluid flow, and water-rock reaction; wherein, the multi-field coupled governing equations of solute transport, fluid flow, and water-rock reaction include the governing equations of solute transport and water-rock reaction in the bedrock system and the governing equations of solute transport and water-rock reaction in the fracture system;

[0109] Numerical solution module: Performs numerical simulation of the crack dissolution evolution process based on the aforementioned governing equations;

[0110] Results output module: Outputs the geometry of the fracture network after dissolution, the distribution of solute concentration, and the changes in conductivity.

[0111] Example 3

[0112] This embodiment provides a readable storage medium storing a computer program that, when invoked by a processor, performs the steps of the method described above.

[0113] Example 4

[0114] This embodiment provides an electronic terminal including a processor and a memory, wherein the memory stores a computer program, and the processor calls the computer program to perform the steps of the method described above.

[0115] It should be understood that, in the embodiments of the present invention, the processor may be a Central Processing Unit (CPU), or it may be other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. The general-purpose processor may be a microprocessor or any conventional processor. The memory may include read-only memory and random access memory, and provides instructions and data to the processor. A portion of the memory may also include non-volatile random access memory. For example, the memory may also store device type information.

[0116] The readable storage medium is a computer-readable storage medium, which can be an internal storage unit of the controller described in any of the foregoing embodiments, such as the controller's hard drive or memory. The readable storage medium can also be an external storage device of the controller, such as a plug-in hard drive, Smart Media Card (SMC), Secure Digital (SD) card, or Flash Card equipped on the controller. Further, the readable storage medium can include both the controller's internal storage unit and external storage devices. The readable storage medium is used to store the computer program and other programs and data required by the controller. The readable storage medium can also be used to temporarily store data that has been output or will be output.

[0117] Based on this understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or all or part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present invention. The aforementioned readable storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0118] It is understood that the same or similar parts in the above embodiments can be referred to each other, and the contents not described in detail in some embodiments can be referred to the same or similar contents in other embodiments.

[0119] Although embodiments of the present invention have been shown and described above, it is understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those skilled in the art can make changes, modifications, substitutions and variations to the above embodiments within the scope of the present invention.

[0120] To further illustrate the technical solution of this application, we will use single-crack dissolution and complex crack network dissolution as application examples, as follows:

[0121] Application Example 1: Dissolution of a Single Crack

[0122] This example uses the simulation method proposed in this invention to simulate the dissolution process of a single crack. The dissolution morphology of the salt rock fracture surface is as follows: Figure 6 As shown, the main parameters of the model are set as follows: diffusion coefficient D = 2.0 × 10⁻⁶. -5 cm 2 / s, the molar mass of the salt rock is M = 58.5 kg / mol, the saturation concentration of the salt rock is Cs = 5.4 mol / L, and the density of the salt rock is =2160kg / m 3 The head difference between the two ends of the rock sample was ∆P=4.7cm, the length of the sample was L=6cm, the initial width was b=0.01cm, the initial concentration of the solution was C0=5mol / L, and the duration of the experiment was 125h. Figure 7 The figure shows the dissolution morphology of salt rock fracture surfaces and a comparison of experimental and simulation results. The calculated curves in the figure show that the calculated dissolution thickness is in excellent agreement with the experimentally obtained dissolution thickness, with a consistency rate of 86%. This indicates that the fracture dissolution simulation method proposed in this invention can well reproduce the dissolution process of fractures in soluble rock masses and has high calculation accuracy.

[0123] Application Example 2: Dissolution of Complex Crack Networks

[0124] This case study examines the dissolution process of a fracture network in a well-developed soluble carbonate layer. The thickness of the soluble layer is 90 meters. Other key simulation parameters are set as follows: water head difference between the top and bottom of the reservoir is 37 meters, reservoir temperature is 10℃, carbon dioxide partial pressure is 0.05 atmospheres, initial fracture aperture is 1.0 mm, initial fracture permeability is 100 mD (considering the roughness of the fractures); initial matrix permeability is 0.1 mD; porosity is 0.15; and local mass exchange coefficient is 1.5 × 10⁻⁶. -3 m / s; molecular diffusion coefficient is 3.6 × 10⁻⁶ m / s.-9 m 2 / s; fluid viscosity is 1.0 cP; fluid density is 1000 km / m 3 Initial specific surface area 5000 m² 2 / m 3 The simulation yielded the distribution of calcium ion concentration, the distribution of dissolution fracture width, the variation of fluid flow rate over time, and the variation of dissolution fracture width with dissolution distance after 100 years. Figure 8 A three-dimensional discrete crack model is presented. Figure 9 The distribution of calcium ions and the width distribution of dissolution fissures were shown. Figure 9 (a) It can be seen that different recharge locations exhibit different calcium ion flow rate history curves under the same head difference due to different frictional resistance; from Figure 9 (b) The width of the dissolution fracture increases significantly near the refill location, but not significantly away from the refill location; from Figure 9 (c) It can be seen that the calcium ion concentration varies significantly throughout the study area, with the lowest concentration near the recharge site; from Figure 9 (d) It can be seen that the most significant crack dissolution is limited to a small area near the replenishment site.

Claims

1. A numerical simulation method for the dissolution process of a three-dimensional crack network, characterized in that, include: S1: Obtain the geological parameters of the target reservoir and construct a geometric model of the three-dimensional fracture network; S2: The geometric model is meshed, and the fracture network is embedded into the bedrock mesh using the embedded discrete fracture method; S3: Determine the multi-field coupled control equations of solute transport-fluid flow-water-rock reaction; wherein, the multi-field coupled control equations of solute transport-fluid flow-water-rock reaction include the control equations of solute transport and water-rock reaction in the bedrock system and the control equations of solute transport and water-rock reaction in the fracture system; Among them, the governing equations for solute transport and water-rock reaction in the multi-field coupled governing equations of solute transport-fluid flow-water-rock reaction are: ； in, For matrix porosity; This represents the concentration of the solute in the matrix. For simulating time; It is the divergence operator; The seepage velocity in the matrix; The dispersion coefficient of the solute in the matrix This refers to the dry density of the matrix. This represents the amount of solute adsorbed in the matrix. ,in This represents the maximum adsorption capacity. It is the adsorption equilibrium constant; For the source and sink terms of solutes in the matrix; This is the water-rock reaction rate term for the solute in the matrix. ,in The reaction rate constant is... The surface area of ​​the bedrock reaction. This represents the equilibrium concentration of the solute. This is the exchange term for solute between the matrix and the crack. ,in The solute exchange coefficient; The governing equations for solute transport and water-rock reaction in the fracture system within the multi-field coupled solute transport-fluid flow-water-rock reaction equations are as follows: ； in, The crack opening; The concentration of the solute in the crack; The seepage velocity in the crack; The dispersion coefficient of the solute in the crack; The dry density of the crack wall; This represents the amount of solute adsorbed on the crack wall surface; These are the source and sink terms of solutes in the cracks; This is the water-rock reaction rate term for solutes in the fracture; This refers to the exchange of solute between connected fracture units. , The solute exchange coefficient between the cracks. For adjacent connected cracks Solute concentration in each unit; S4: Numerical simulation of the crack dissolution evolution process based on the aforementioned governing equations; S5: Outputs the geometric morphology, solute concentration distribution, and changes in conductivity of the fracture network after dissolution.

2. The method according to claim 1, characterized in that, The specific process of S1 is as follows: S11: Based on the seismic interpretation, well logging data and core observation results of the target reservoir, statistically analyze the probability distribution characteristics of key geometric parameters of fractures and establish the probability density function of each geometric parameter; S12: Based on the probability density function of each geometric parameter, the Monte Carlo simulation method is used to generate three-dimensional crack elements in batches to obtain cracks, and generation constraints are set. S13: Using a spatial geometric intersection detection algorithm, the topological relationship between the generated cracks is verified, cracks that do not meet the constraints are removed, nodes of the intersecting cracks are split, the spatial coordinates and connectivity of the crack intersection points are clarified, and a crack database is obtained. S14: Based on core experiments, well logging interpretation, and field data, set the physical and mechanical parameters for each fracture; S15: Integrate the information from S11 to S14 to obtain the geometric model of the three-dimensional crack network.

3. The method according to claim 1, characterized in that, The specific process of S2 is as follows: S21: Use structured grids to divide the background matrix, set the grid size, number each bedrock grid unit and store the corresponding physical and mechanical parameters; S22: Determine the intersection relationship between the crack and the bedrock grid unit, and divide the crack unit; S23: Construct the topological association between bedrock mesh cells and fracture cells, and store geometric and physical parameters to obtain a fracture network embedded in the bedrock mesh.

4. The method according to claim 3, characterized in that, The specific process of S22 is as follows: S221: For each crack, traverse all bedrock grid cells, calculate the intersection area between the crack surface and the grid cell, and if the area of ​​the intersection area is greater than zero, then the grid cell is determined to be a crack-containing grid cell, and the correspondence between the crack-containing grid cell number and the crack number is recorded. S222: Based on the intersection of the crack and the bedrock grid, the crack is divided into several crack units, and each crack unit corresponds to a crack-containing grid unit. S223: The matrix mesh portion embedded with the crack is mapped to a weighted graph using a continuous projection algorithm. The Dijkstra algorithm is then applied to find the shortest path between the two endpoints of the crack. This path is mapped back to the mesh to determine the surface receiving the projection, thereby completing the crack element division.

5. The method according to claim 1, characterized in that, The specific process of S4 is as follows: S41: Initialize simulation parameters, set the total simulation duration and initial time step, adopt implicit time discretization scheme, and initialize the solute concentration field, pressure field, velocity field and aperture field of each grid element in the matrix and crack; S42: The solute transport and water-rock reaction coupling equation is split into convection operator, dispersion operator and adsorption-reaction operator by operator splitting method, and spatial discretization is performed by integral finite difference method respectively. The solute concentration field is updated by iterative coupling. S43: Constructing a system of linear algebraic equations based on discrete results The coefficient matrix Based on the geometric parameters and coupling conditions of the matrix mesh and crack elements, the solution vector x includes fluid pressure and solute concentration. The right-hand term vector b is constructed based on boundary conditions and source-sink terms. The sparse matrix is ​​solved using the generalized minimum residual method to obtain the pressure field and concentration field at the current time step. S44: Calculate the water-rock reaction rate of solute in each matrix grid and fracture unit based on the concentration field, update the bedrock porosity, fracture aperture and permeability according to the reaction rate, and feed the updated physical property parameters back to the governing equations of the next time step; S45: Determine whether the solution at the current time step has converged. If it has converged, proceed to the next time step. If it has not converged, reduce the time step size and recalculate S42 to S44. Repeat the iteration until the preset total simulation time is reached.

6. A numerical simulation system for the dissolution process of a three-dimensional crack network, characterized in that, The system is used to perform the steps of the method according to any one of claims 1-5, including: Geometric model construction module: Obtains geological parameters of the target reservoir and constructs a geometric model of the three-dimensional fracture network; Mesh generation and embedding module: The geometric model is meshed, and the fracture network is embedded into the bedrock mesh using the embedded discrete fracture method; The governing equation determination module determines the multi-field coupled governing equations of solute transport, fluid flow, and water-rock reaction; wherein, the multi-field coupled governing equations of solute transport, fluid flow, and water-rock reaction include the governing equations of solute transport and water-rock reaction in the bedrock system and the governing equations of solute transport and water-rock reaction in the fracture system; Numerical solution module: Performs numerical simulation of the crack dissolution evolution process based on the aforementioned governing equations; Results output module: Outputs the geometry of the fracture network after dissolution, the distribution of solute concentration, and the changes in conductivity.

7. A readable storage medium, characterized in that: A computer program is stored, which, when invoked by a processor, performs the steps of the method according to any one of claims 1-5.

8. An electronic terminal, characterized in that: It includes a processor and a memory, the memory storing a computer program, the processor calling the computer program to perform the steps of the method according to any one of claims 1-5.

Images