A method for reservoir-scale reservoir acidification simulation

Abstract

The application discloses a kind of reservoir acidification simulation methods for oil reservoir scale, it is related to oil and gas field development technical field.The acid rock reaction flow system of reservoir is decoupled into three-dimensional matrix continuous domain and one-dimensional discrete wormhole topology net, constructs pressure and flow coupling equation, acid liquid transport and reaction equation and geometric shape evolution equation to carry out reservoir acidification simulation, the pressure field of three-dimensional matrix continuous domain and one-dimensional discrete wormhole topology net is fully implicit simultaneous solution, obtains the fluid pressure distribution in reservoir, determines the volume flow and flow direction in one-dimensional discrete wormhole topology net, constructs directed acyclic graph, executes topological sorting and generates downstream calculation sequence calculation to obtain the acid liquid concentration field in one-dimensional discrete wormhole topology net, in combination with geometric shape evolution equation, the geometric parameter of one-dimensional discrete wormhole topology net is updated explicitly, generates new tip node and new wormhole pipe section using geometric topology update, until tip node breaks through production boundary, obtain the simulation result of reservoir acidification.

Description

Technical Field

[0001] This invention relates to the field of oil and gas field development technology, and specifically to a reservoir acidification simulation method at the reservoir scale. Background Technology

[0002] Carbonate reservoir matrix acidizing is a key technology for achieving high-speed and efficient oil and gas field development, removing near-wellbore contamination, and increasing single-well productivity. During reservoir acidizing operations, the injected acid undergoes a violent nonlinear multi-field coupled dissolution reaction with the rock skeleton, spontaneously evolving into a dendritic macropore network with extremely high conductivity, forming vermiforms. Accurately predicting the spatial topology, propagation distance, and optimal acid dosage of vermiforms is crucial for optimizing construction parameters.

[0003] In recent years, theories such as the two-scale continuum model have matured, achieving unprecedented success in core-scale mechanism studies and optimal injection rate prediction by coupling Darcy's macroscopic flow with the reaction dynamics of microscopic pore interfaces. However, when attempting to extend the two-scale continuum model to the meter- or even hundred-meter-scale of real reservoirs using conventional reservoir numerical methods, it encounters limitations from both physical theory and computational power. First, reservoir numerical methods can essentially only resolve long-wavelength phenomena at the grid scale, but the absolute diameter of wormholes generated by reservoir acidification does not increase proportionally with the scaling of the system scale; instead, it is strictly locked by the characteristic correlation length of the intrinsic heterogeneity of the formation. Therefore, in order to accurately capture the dissolution front in three-dimensional space, the computational grid size of the dual-scale continuous medium model must not exceed this characteristic length, which usually limits the computational grid size to less than 3 cm. Secondly, the target treatment zone of actual reservoir acidization is usually tens of meters in size. This strict grid size rigidity limit forces the single-well near-field model to generate a massive computational grid of billions. Under this computational scale, if the traditional fully implicit method is used for solving, the assembly and solution of massive sparse Jacobian matrices at each time step during the simulation calculation will cause a cubic-level computational explosion, resulting in long machine time, memory crashes, and difficulty in applying it to actual engineering simulation. Conversely, if the traditional coarse grid parameter equivalent amplification method is used as a compromise, the geometric topological characteristics of the wormhole network will be completely lost due to the highly unstable dynamic heterogeneity of acid rock reaction, making it impossible to intuitively assess whether the acid has penetrated the contaminated zone.

[0004] Therefore, there is an urgent need to propose a reservoir acidification simulation method at the reservoir scale to break the theoretical constraints and computational disasters caused by the grid size limitation due to the characteristic length lock-in process in reservoir acidification simulation. This method would enable dynamic tracking and geometric growth simulation of wormhole height during reservoir acidification, and solve the engineering challenges of intelligent simulation of macroscale acidification. Summary of the Invention

[0005] This invention aims to solve the above-mentioned problems and proposes a reservoir acidification simulation method at the reservoir scale. By decoupling the spatial topology from the physical system, the complex reaction flow system is separated into a coarse three-dimensional matrix continuous domain that carries low-velocity seepage and a one-dimensional discrete vermiform topology that encapsulates high conductivity. Combined with the mapping and embedding at the mathematical integration level, the absolute conflict of the simulation grid scale is overcome, and a sub-grid is generated inside the three-dimensional coarse grid. This realistically restores the geometric growth and development of vermiforms in the reservoir and realizes the simulation of reservoir acidification at the reservoir scale.

[0006] The present invention adopts the following technical solution:

[0007] A reservoir acidification simulation method at the reservoir scale includes the following steps:

[0008] Step 1: Obtain the acid-rock reaction flow system of the reservoir to be simulated, and decouple it into a three-dimensional matrix continuous domain and a one-dimensional discrete wormhole topology network. Construct the pressure and flow coupling equation, acid transport and reaction equation, and geometric evolution equation to simulate reservoir acidification.

[0009] Step 2: Perform a fully implicit simultaneous solution to the pressure fields of the three-dimensional matrix continuous domain and the one-dimensional discrete wormhole topology network to obtain the fluid pressure distribution inside the three-dimensional matrix continuous domain and the one-dimensional discrete wormhole topology network in the reservoir, and determine the volumetric flow rate and flow direction inside the one-dimensional discrete wormhole topology network.

[0010] Step 3: Extract the flow direction of the fluid in the one-dimensional discrete wormhole topology network to construct a directed acyclic graph, perform topological sorting to generate a downstream calculation sequence, and solve the acid concentration at each node in the one-dimensional discrete wormhole topology network based on the first-order implicit upwind scheme along the downstream calculation sequence to obtain the acid concentration field in the one-dimensional discrete wormhole topology network.

[0011] Step 4: Based on the acid concentration field and geometric evolution equations within the one-dimensional discrete wormhole topology network, explicitly update the geometric parameters of the one-dimensional discrete wormhole topology network, including the radial expansion dimension of each pipe segment and the cumulative extension length of the tip node in the one-dimensional discrete wormhole topology network.

[0012] Step 5: Determine whether the preset cross-grid growth conditions are met based on the cumulative extension length of the tip node. Perform geometric topology update to generate new tip nodes and wormhole segments, update the time step, and repeat steps 2 to 4 until the maximum acid injection time is reached or the preset acidization breakthrough conditions are met. Then stop the reservoir acidization simulation and obtain the reservoir acidization simulation results.

[0013] Preferably, step 1 includes the following steps:

[0014] Step 1.1: Construct pressure-flow coupling equations, including three-dimensional matrix seepage equations, one-dimensional wormhole flow equations, and one-dimensional-three-dimensional fluid filtration coupling equations;

[0015] The fluid in the three-dimensional matrix flow equation is set as a slightly compressible single-phase flow. Combined with the reservoir pressure, the three-dimensional matrix flow equation is determined as follows:

[0016] ;

[0017] In the formula, For matrix porosity; This is the overall compression coefficient; For matrix pressure; For time; It is a divergence operator; For matrix permeability; The viscosity of the acid solution; For the three-dimensional Hamiltonian gradient operator; The filtration loss per unit length of wormhole into the substrate; Let be the spatial Dirac distribution function; It is a spatial position vector; This refers to the spatial location of the worm hole axis;

[0018] The fluid in the one-dimensional wormhole flow equation is set to be incompressible, the wormhole is set as a dynamically variable diameter circular pipe, and the fluid flows along the one-dimensional arc length of the wormhole using a Poiseuille flow. The one-dimensional wormhole flow equation is determined as follows:

[0019] ;

[0020] In the formula, For pipe flow volumetric flow rate, ; The dynamic radius of the worm pore; The arc length of the worm hole; This refers to the pressure inside the wormhole tube;

[0021] The one-dimensional-three-dimensional fluid filtration-loss coupling equation is established based on the Peaceman well index model, and its expression is:

[0022] ;

[0023] In the formula, The coupling conduction index, ; It is the natural logarithm function; The equivalent discharge radius of the three-dimensional matrix mesh. ,in, , , These represent the spatial dimensions of the three-dimensional matrix mesh. , , Dimensions;

[0024] Step 1.2: Ignoring longitudinal molecular diffusion and focusing solely on the flow-dominant principle, the acid concentration is calculated within a one-dimensional discrete worm-pore topology network. The acid transport and reaction equations are established as follows:

[0025] ;

[0026] In the formula, This refers to the acid concentration. The mass transfer coefficient is . ; The effective molecular diffusion coefficient; It is an asymptotic Sherwood number; The pore Reynolds number; It is a Schmidt number;

[0027] Step 1.3, based on the law of conservation of mass, convert the acid consumed on the side wall of the wormhole into an increase in the radius of the wormhole segment, resulting in:

[0028] ;

[0029] In the formula, The volumetric solubility of the acid solution;

[0030] The unreacted acid reaching the tip node is described using a tip extension equation, where the tip node is the node at the very front of a one-dimensional discrete wormhole topology. The extension velocity of the unreacted acid advancing through the borehole in front of the tip node is determined, and a geometric evolution equation is established to overcome the grid size limitations in reservoir acidization simulation, enabling the simulation of sub-grid scale dissolution growth. The results are as follows:

[0031] ;

[0032] In the formula, The rate of spread of the unreacted acid solution; The cumulative extension length of the tip node; The volumetric flow rate reaching the tip node; The acid concentration required to reach the tip node; The initial germination radius is the radius of the newly generated micropore channels.

[0033] Preferably, in step 2, a hidden-pressure, concentrated discretization scheme is used to connect the three-dimensional matrix mesh in the three-dimensional matrix continuous domain with its intersecting one-dimensional wormhole segments. Backward Euler time-difference discretization is then used to construct a large-scale sparse linear equation system. ,in, It is a large sparse coefficient matrix, composed of the coefficients of the three-dimensional matrix mesh discretization equation and the one-dimensional wormhole segment discretization equation. Let be the pressure field vector to be determined, containing The pressure of all 3D matrix meshes and wormhole segments at all times. It is a constant vector;

[0034] The discretization equation of the three-dimensional matrix mesh is:

[0035] ;

[0036] in,

[0037] ;

[0038] In the formula, , All are serial numbers of the three-dimensional matrix mesh; For the first The control volume of a three-dimensional matrix mesh; For the first Porosity of a three-dimensional matrix mesh; This is the overall compression coefficient; For the current time step Pressure values ​​of a three-dimensional matrix mesh; Each is the previous time step. Pressure values ​​of a three-dimensional matrix mesh; The discrete time step; A set of topologically adjacent grids or nodes; The fluid conductivity between adjacent three-dimensional matrix meshes; For the current time step Pressure values ​​of a three-dimensional matrix mesh; , These are all serial numbers of one-dimensional wormhole segments. ; For the first A collection of one-dimensional wormhole segments contained in a three-dimensional matrix mesh; For the first Discrete form well index of a one-dimensional wormhole segment; For the current time step Pressure value of a one-dimensional wormhole segment; For the first The coupling conduction index of a one-dimensional wormhole segment; For the first The physical length of a one-dimensional wormhole segment; For the first The matrix permeability of a three-dimensional matrix grid; The viscosity of the acid solution;

[0039] The discrete equation for the one-dimensional wormhole segment is:

[0040] ;

[0041] in,

[0042] ;

[0043] In the formula, The conductivity of the pipe section; For the first One-dimensional wormhole segment and the first The explicit radius between one-dimensional wormhole segments; For the first One-dimensional wormhole segment and the first Center distance between one-dimensional wormhole segments; For the current time step Pressure value of a one-dimensional wormhole segment.

[0044] Preferably, step 3 includes the following steps:

[0045] Step 3.1: Based on the fluid pressure distribution within the three-dimensional matrix continuum and the one-dimensional discrete wormhole topology network in the reservoir at the current time step, and combined with the volumetric flow rate and flow direction within each one-dimensional wormhole segment in the one-dimensional discrete wormhole topology network, and based on the physical property of the fluid flowing in a tree-like divergent manner, the one-dimensional discrete wormhole topology network is abstracted into a directed acyclic graph. A graph theory algorithm is used to perform topological sorting on the directed acyclic graph to generate a node substitution sequence. In the node substitution sequence, each node is arranged in order of the fluid flow direction from the injection well to the tip of each one-dimensional wormhole segment.

[0046] Step 3.2, based on the first-order implicit upwind scheme, the discrete mass conservation equation for the one-dimensional wormhole segment in the one-dimensional discrete wormhole topology is determined as follows:

[0047] ;

[0048] In the formula, For the node index in the directed acyclic graph; For the previous time step The control volume of the one-dimensional wormhole segment where each node is located; For the current time step Acid concentration at each node; For the previous time step Acid concentration at each node; outflow The set of flow rates for all pipe segments at each node; For outflow boundary flow; For the current time step Acid concentration at each node; Inflow is the first The set of flow rates for all pipe segments at each node; For inflow to the boundary flow; The acid concentration at the upstream node; For the previous time step The radius of the one-dimensional wormhole segment where each node is located; For the first The physical length of the one-dimensional wormhole segment where each node is located; This is the mass transfer coefficient of the previous time step; For the current time step Filtration at each node;

[0049] Step 3.3: Substitute each node in the generated node sequence into the discrete mass conservation equation of the one-dimensional wormhole pipe segment for solution. During the solution process, the acid concentration of the upstream node corresponding to each node in the directed acyclic graph is pre-extracted. This is then transformed into a known constant term, causing the discrete mass conservation equation for the one-dimensional wormhole pipe segment to degenerate into one containing only a single unknown. The acid concentration at each node in the current directed acyclic graph is calculated using a linear algebraic equation in one variable.

[0050] Preferably, in step 5, during the reservoir acidification simulation, the radius of each one-dimensional wormhole segment and the cumulative extension length of each tip node in the one-dimensional discrete wormhole topology network are updated at each time step. Then, a geometric topology evolution check is performed on the one-dimensional discrete wormhole topology network to determine whether the cumulative extension length of all tip nodes in the one-dimensional discrete wormhole topology network has reached a preset growth threshold. If any tip node has reached the growth threshold, new nodes and connected segments are added to the directed acyclic graph, and the cumulative extension distance corresponding to the tip node that has reached the growth threshold is deducted to generate a subgrid, thus completing the macroscopic cross-grid topology update and bifurcation evolution of the one-dimensional discrete wormhole topology network.

[0051] Preferably, the formula for updating the radius of the one-dimensional wormhole segment in the one-dimensional discrete wormhole topology network is as follows:

[0052] ;

[0053] In the formula, For the current time step The radius of a one-dimensional wormhole segment; For the previous time step The radius of a one-dimensional wormhole segment; The discrete time step; The volumetric solubility of the acid solution; The mass transfer coefficient of the previous time step; For the current time step The acid concentration within a one-dimensional wormhole segment; For matrix porosity;

[0054] The formula for calculating the cumulative extension length update of the tip node is as follows:

[0055] ;

[0056] In the formula, The cumulative extension length of the current time step's tip node; This is the cumulative extension length of the tip node in the previous time step; The rate of spread of the unreacted acid solution.

[0057] The present invention has the following beneficial effects:

[0058] (1) This invention proposes a reservoir acidification simulation method at the reservoir scale. It abandons the cumbersome fully implicit architecture and adopts a discrete format of implicit pressure and explicit concentration and geometry. Combining the closed-loop physical characteristic that fluids in underground porous media networks must flow unidirectionally according to the pressure gradient, it introduces directed acyclic graphs and topological sorting algorithms from computer graph theory to rearrange and organize the millions of discrete grids into a single linear operation sequence that strictly follows the causal flow direction of the fluid. This allows the algorithm engine to progress block by block along the streamline direction, independently solving the concentration field under a first-order upwind explicit format. This invention not only reduces the computational cost of inverting large matrix simultaneous equations with extremely low linear time complexity but also completely breaks free from the limitations of the CFL condition on small time steps, achieving rapid simulation of reservoir acidification.

[0059] (2) This invention proposes a reservoir acidification simulation method at the reservoir scale, which integrates the physical advantages of embedded dimensionality reduction model and the high-speed computing capability of computer graph theory topology sorting. It effectively overcomes the shortcomings of traditional numerical simulation methods in reservoir-scale acidification simulation, such as grid curse and inability to depict the real dissolution morphology of reservoirs. It provides a new idea for solving the problem of intelligent simulation engineering of reservoir acidification construction at the macro scale, effectively improves the design efficiency and economic benefits of oilfield acidification transformation, and has important engineering application value and broad industry promotion prospects. Attached Figure Description

[0060] Figure 1 This is a flowchart of a reservoir acidification simulation method at the reservoir scale according to the present invention.

[0061] Figure 2 This is a schematic diagram of the reservoir acidification simulation results obtained using the present invention. In the diagram, the black areas represent stagnant dissolution blind ends, and the red areas represent actively extending main trunk vermiform holes. Detailed Implementation

[0062] The specific embodiments of the present invention will be further described below with reference to the accompanying drawings.

[0063] Example 1

[0064] This invention proposes a reservoir acidification simulation method at the reservoir scale, such as... Figure 1 As shown, the specific steps include:

[0065] Step 1: Obtain the acid-rock reaction flow system of the reservoir to be simulated, and decouple it into a three-dimensional matrix continuum and a one-dimensional discrete wormhole topology network. Construct the pressure-flow coupling equation, acid transport and reaction equation, and geometric evolution equation to simulate reservoir acidification, including the following sub-steps:

[0066] Step 1.1: Construct pressure-flow coupling equations, which include three-dimensional matrix seepage equations, one-dimensional wormhole flow equations, and one-dimensional-three-dimensional fluid filtration coupling equations.

[0067] The fluid in the three-dimensional matrix flow equation is set as a slightly compressible single-phase flow. Combined with the reservoir pressure, the three-dimensional matrix flow equation is determined as follows:

[0068] ;

[0069] In the formula, For matrix porosity; This is the overall compression coefficient; For matrix pressure; For time; It is a divergence operator; For matrix permeability; The viscosity of the acid solution; This is a three-dimensional Hamiltonian gradient operator used for spatial derivative and divergence calculations; The filtration loss per unit length of wormhole into the substrate; Let be the spatial Dirac distribution function; It is a spatial position vector; This refers to the spatial location of the worm hole axis.

[0070] The fluid in the one-dimensional wormhole flow equation is set to be incompressible, the wormhole is set as a dynamically variable diameter circular pipe, and the fluid flows along the one-dimensional arc length of the wormhole using a Poiseuille flow. The one-dimensional wormhole flow equation is determined as follows:

[0071] ;

[0072] in,

[0073] ;

[0074] In the formula, This refers to the volumetric flow rate in the pipe. The dynamic radius of the worm pore; The arc length of the worm hole; This refers to the pressure inside the wormhole tube.

[0075] The one-dimensional-three-dimensional fluid filtration-loss coupling equation is established based on the Peaceman well index model, and the one-dimensional-three-dimensional fluid filtration-loss coupling equation is set as follows:

[0076] ;

[0077] in,

[0078] ;

[0079] In the formula, The coupling conduction index; It is the natural logarithm function; The equivalent discharge radius of the three-dimensional matrix mesh. ,in, , , These represent the spatial dimensions of the three-dimensional matrix mesh. , , The size.

[0080] Step 1.2: Ignoring longitudinal molecular diffusion and focusing solely on the flow pole, the acid concentration is calculated within a one-dimensional discrete wormhole topology network, and the acid transport and reaction equations are established.

[0081] The acid transport and reaction equations are set as follows:

[0082] ;

[0083] in,

[0084] ;

[0085] In the formula, This refers to the acid concentration. The mass transfer coefficient; The effective molecular diffusion coefficient; It is an asymptotic Sherwood number; The pore Reynolds number; It is a Schmitt number.

[0086] Step 1.3, based on the law of conservation of mass, convert the acid consumed on the side wall of the wormhole into an increase in the radius of the wormhole segment, resulting in:

[0087] ;

[0088] In the formula, This represents the volumetric solubility of an acid solution, used to express the volume of pure rock skeleton that can be dissolved by a unit mole of acid solution.

[0089] The unreacted acid reaching the tip node is described using a tip extension equation, where the tip node is the node at the very front of a one-dimensional discrete wormhole topology. The extension velocity of the unreacted acid advancing through the borehole in front of the tip node is determined, and a geometric evolution equation is established to overcome the grid size limitations in reservoir acidization simulation, enabling the simulation of sub-grid scale dissolution growth. The results are as follows:

[0090] ;

[0091] In the formula, The rate of spread of the unreacted acid solution; The cumulative extension length of the tip node; The volumetric flow rate reaching the tip node; The acid concentration required to reach the tip node; The initial germination radius is the radius of the newly generated micropore channels.

[0092] Step 2: Perform a fully implicit simultaneous solution to the pressure fields of the three-dimensional matrix continuous domain and the one-dimensional discrete wormhole topology network to obtain the fluid pressure distribution inside the three-dimensional matrix continuous domain and the one-dimensional discrete wormhole topology network in the reservoir, and determine the volumetric flow rate and flow direction inside the one-dimensional discrete wormhole topology network.

[0093] Specifically, a discretization scheme with implicit pressure and explicit density is adopted. Since the conductivity of the wormhole is proportional to the fourth power of the dynamic radius of the wormhole and has extremely strong nonlinearity, in order to ensure the stability of the calculation, the three-dimensional matrix mesh in the three-dimensional matrix continuous domain is connected with the intersecting one-dimensional wormhole line segments. Backward Euler time difference discretization is used to construct a large-scale sparse linear equation system. ,in, It is a large sparse coefficient matrix, composed of the coefficients of the three-dimensional matrix mesh discretization equation and the one-dimensional wormhole segment discretization equation. Let be the pressure field vector to be determined, containing The pressure of all 3D matrix meshes and wormhole segments at all times. It is a constant vector.

[0094] The discretization equation of the three-dimensional matrix mesh is:

[0095] ;

[0096] in,

[0097] ;

[0098] In the formula, , All are serial numbers of the three-dimensional matrix mesh; For the first The control volume of a three-dimensional matrix mesh; For the first Porosity of a three-dimensional matrix mesh; This is the overall compression coefficient; For the current time step Pressure values ​​of a three-dimensional matrix mesh; Each is the previous time step. Pressure values ​​of a three-dimensional matrix mesh; The discrete time step; A set of topologically adjacent grids or nodes; The fluid conductivity between adjacent three-dimensional matrix meshes; For the current time step Pressure values ​​of a three-dimensional matrix mesh; , These are all serial numbers of one-dimensional wormhole segments. ; For the first A collection of one-dimensional wormhole segments contained in a three-dimensional matrix mesh; For the first Discrete form well index of a one-dimensional wormhole segment; For the current time step Pressure value of a one-dimensional wormhole segment; For the first The coupling conduction index of a one-dimensional wormhole segment; For the first The physical length of a one-dimensional wormhole segment; For the first The matrix permeability of a three-dimensional matrix grid; This refers to the viscosity of the acid solution.

[0099] The discrete equation for the one-dimensional wormhole segment is:

[0100] ;

[0101] in,

[0102] ;

[0103] In the formula, The conductivity of the pipe section; For the first One-dimensional wormhole segment and the first The explicit radius between one-dimensional wormhole segments; For the first One-dimensional wormhole segment and the first Center distance between one-dimensional wormhole segments; For the current time step Pressure value of a one-dimensional wormhole segment.

[0104] Step 3: Extract the flow direction of the fluid in the one-dimensional discrete wormhole topology network to construct a directed acyclic graph. Execute the topology sorting algorithm to generate a downstream calculation sequence. Based on the first-order implicit upwind scheme, solve for the acid concentration at each node in the one-dimensional discrete wormhole topology network along the downstream calculation sequence to obtain the acid concentration field in the one-dimensional discrete wormhole topology network. This includes the following sub-steps:

[0105] Step 3.1: Based on the fluid pressure distribution within the three-dimensional matrix continuum and the one-dimensional discrete wormhole topology network in the reservoir at the current time step, and combined with the volumetric flow rate and flow direction within each one-dimensional wormhole segment in the one-dimensional discrete wormhole topology network, and based on the physical property of the fluid flowing in a tree-like divergent manner, the one-dimensional discrete wormhole topology network is abstracted into a directed acyclic graph. A graph theory algorithm is used to perform topological sorting on the directed acyclic graph to generate a node substitution sequence. In the node substitution sequence, each node is arranged in order of the fluid flow direction from the injection well to the tip of each one-dimensional wormhole segment.

[0106] Step 3.2, based on the first-order implicit upwind scheme, the discrete mass conservation equation for the one-dimensional wormhole segment in the one-dimensional discrete wormhole topology is determined as follows:

[0107] ;

[0108] In the formula, For the node index in the directed acyclic graph; For the previous time step The control volume of the one-dimensional wormhole segment where each node is located; For the current time step Acid concentration at each node; For the previous time step Acid concentration at each node; outflow The set of flow rates for all pipe segments at each node; For outflow boundary flow; For the current time step Acid concentration at each node; Inflow is the first The set of flow rates for all pipe segments at each node; For inflow to the boundary flow; The acid concentration at the upstream node; For the previous time step The radius of the one-dimensional wormhole segment where each node is located; For the first The physical length of the one-dimensional wormhole segment where each node is located; This is the mass transfer coefficient of the previous time step; For the current time step Filter loss at each node.

[0109] Step 3.3: Substitute each node in the node substitution sequence generated in Step 3.1 into the discrete mass conservation equation of the one-dimensional wormhole pipe segment in Step 3.2 for solution. During the solution process, the acid concentration of the upstream node corresponding to each node in the directed acyclic graph is extracted in advance. This is then transformed into a known constant term, causing the discrete mass conservation equation for the one-dimensional wormhole pipe segment to degenerate into one containing only a single unknown. The linear algebraic equation in one variable is used to calculate the acid concentration at each node in the current directed acyclic graph, thus ignoring the limitation of the CFL time step (i.e., the maximum safe time step in the explicit numerical simulation determined based on the Courant-Friedrichs-Lewy condition).

[0110] Step 4: Based on the acid concentration field and geometric evolution equations within the one-dimensional discrete wormhole topology, explicitly update the geometric parameters of the one-dimensional discrete wormhole topology, including the radial expansion dimension of each pipe segment and the cumulative extension length of the tip node in the one-dimensional discrete wormhole topology.

[0111] Step 5: Determine whether the preset cross-grid growth conditions are met based on the cumulative extension length of the tip node. Perform geometric topology update to generate new nodes and pipe segments, update the time step, and repeat steps 2 to 4 until the maximum acid injection time is reached or the preset acidization breakthrough conditions are met. Then stop the reservoir acidization simulation and obtain the reservoir acidization simulation results.

[0112] Furthermore, during the reservoir acidification simulation, at each time step, the radius of each one-dimensional wormhole segment and the cumulative extension length of each tip node in the one-dimensional discrete wormhole topology network are updated. Then, a geometric topology evolution check is performed on the one-dimensional discrete wormhole topology network to determine whether the cumulative extension length of all tip nodes in the one-dimensional discrete wormhole topology network has reached a preset growth threshold. If any tip node has reached the growth threshold, new nodes and connected segments are added to the directed acyclic graph, and the cumulative extension distance corresponding to the tip node that has reached the growth threshold is deducted to generate a sub-mesh, thus completing the macroscopic cross-mesh topology update and bifurcation evolution of the one-dimensional discrete wormhole topology network.

[0113] Specifically, the formula for updating the radius of the one-dimensional wormhole segment in the one-dimensional discrete wormhole topology network is as follows:

[0114] ;

[0115] In the formula, For the current time step The radius of a one-dimensional wormhole segment; For the previous time step The radius of a one-dimensional wormhole segment; The discrete time step; The volumetric solubility of the acid solution; The mass transfer coefficient of the previous time step; For the current time step The acid concentration within a one-dimensional wormhole segment; This refers to the matrix porosity.

[0116] The formula for calculating the cumulative extension length update of the tip node is as follows:

[0117] ;

[0118] In the formula, The cumulative extension length of the current time step's tip node; This is the cumulative extension length of the tip node in the previous time step; The rate of spread of the unreacted acid solution.

[0119] Example 2

[0120] To verify the feasibility and superiority of the reservoir acidization simulation method at the reservoir scale described in Example 1, this example applies the reservoir acidization simulation method at the reservoir scale described in Example 1 to a three-dimensional embedded discrete wormhole model established based on an actual reservoir to conduct hydrochloric acid acidization simulation.

[0121] In this embodiment, the size of the three-dimensional matrix mesh in the three-dimensional reservoir model is 0.4m × 0.4m × 0.4m, the total number of meshes is 64,000, and the three-dimensional matrix mesh dimensions in all three directions are set to 0.01m. The equivalent discharge radius of the three-dimensional matrix mesh for A three-dimensional heterogeneous permeability field with random variance is generated using a random seed, triggering competitive bifurcation in the wormhole space. The matrix permeability in the three-dimensional reservoir model... Set as matrix porosity Set to 0.15, overall compression ratio Set as ; Set up a simulated acid solution in the three-dimensional reservoir model, wherein the viscosity of the simulated acid solution is... Set as Effective molecular diffusion coefficient Set as kinematic viscosity Set as Schmidt number Set as Volume solubility of acid Set as Normalized concentration of injected acid solution Set as The initial germination radius of newly generated micropore channels Set as The growth threshold was set to 0.8 times the size of the 3D matrix mesh, and the injection end pressure of the 3D reservoir model was set to... Production end pressure set The hydrochloric acid acidification simulation terminated when the tip of the worm hole broke through the production boundary. Adaptive time step control was used, with the maximum time step limit set to 1.0s and the random seed set to 42 to ensure the experiment was repeatable.

[0122] Based on the reservoir acidification simulation method at the reservoir scale described in Example 1, this three-dimensional reservoir model is used to simulate reservoir acidification. The specific process is as follows:

[0123] Step 1: The acid-rock reaction flow system of the 3D reservoir model is decoupled into a 3D matrix continuum and a 1D discrete wormhole topology. Pressure-flow coupling equations, acid transport and reaction equations, and geometric evolution equations are constructed for reservoir acidification simulation. The pressure-flow coupling equations include a 3D matrix seepage equation, a 1D wormhole flow equation, and a 1D-3D fluid filtration coupling equation. The fluid in the 3D matrix seepage equation is set as a slightly compressible single-phase flow. The 1D wormhole flow equation follows Poiseuille flow laws. The 1D-3D fluid filtration coupling equation is based on the Peaceman well index model. The acid transport and reaction equation ignores longitudinal molecular diffusion and is only applied to solve for acid concentration within the 1D discrete wormhole topology. In the geometric evolution equation, based on mass conservation, the acid consumed on the lateral side of the wormhole wall is converted into an increase in the radius of the wormhole segment. A radial expansion equation is set to achieve sub-grid-scale expansion of the 3D wormhole segment. A tip extension equation is also set to ensure that wormhole growth exceeds grid size limitations and achieves sub-grid-scale dissolution.

[0124] Step 2: The pressure fields of the three-dimensional matrix continuous domain and the one-dimensional discrete wormhole topology are solved implicitly. A 7-point difference scheme is used to preassemble the three-dimensional matrix conductivity matrix. Combined with the superposition of one-dimensional wormhole pipe segment flow conduction and one-dimensional-three-dimensional filter loss coupling terms, a large sparse linear equation system is constructed. A constant pressure boundary is applied at the injection end. During the solution process, the inflow and outflow flow rates of each wormhole pipe segment are used to lock the acid flow direction to avoid acid backflow. The fluid pressure distribution inside the three-dimensional matrix continuous domain and the one-dimensional discrete wormhole topology in the reservoir is obtained, and the volumetric flow rate and flow direction inside the one-dimensional discrete wormhole topology are determined.

[0125] Step 3: Based on the volumetric flow rate and flow direction inside the one-dimensional discrete wormhole topology network, the one-dimensional discrete wormhole topology network is abstracted into a directed acyclic graph. Topological sorting is performed, and a downstream calculation sequence is generated in the order from the injection end to each tip. Based on the first-order implicit upwind scheme, the acid concentration at each node in the one-dimensional discrete wormhole topology network is solved along the downstream calculation sequence to obtain the acid concentration field in the one-dimensional discrete wormhole topology network.

[0126] Step 4: Based on the acid concentration field and geometric evolution equations within the one-dimensional discrete wormhole topology network, explicitly update the geometric parameters of the one-dimensional discrete wormhole topology network, including the radial expansion dimension of each pipe segment and the cumulative extension length of the tip node in the one-dimensional discrete wormhole topology network.

[0127] Step 5: Determine whether the preset cross-grid growth conditions are met based on the cumulative extension length of the tip node. Perform geometric topology update to generate new nodes and pipe segments, update the time step, and repeat steps 2 to 4 until acidization breakthrough occurs. Stop the reservoir acidization simulation. Otherwise, update the one-dimensional discrete wormhole topology network structure and cyclically advance the simulation until acidization breakthrough occurs. This completes the macroscopic cross-grid topology update and bifurcation evolution of the one-dimensional discrete wormhole topology network, and obtains the reservoir acidization simulation results.

[0128] In this embodiment, all active tip nodes are traversed. If the cumulative extension length of a tip node reaches 0.008m, a bifurcation is triggered. The basic forward direction is determined by extracting the three-dimensional pressure gradient of the mesh where the tip node is located. Then, a three-dimensional local basis is constructed based on Gram-Schmidt orthogonalization, generating two cones with random perturbations in the bifurcation direction. New tip nodes and new wormhole segments are generated, with the initial radius of the new wormhole segments set to 0.2mm. The wormhole topology and active tip list are updated, and the time step is advanced. Steps 2 to 4 are repeated until the wormhole tip breaks through the production boundary, at which point the reservoir acidization simulation is terminated, and the reservoir acidization simulation results are output. Figure 2 As shown.

[0129] Depend on Figure 2 It can be seen that the wormhole network simulated by the method of the present invention exhibits obvious fractal characteristics. The main wormholes extend along the direction of the maximum pressure gradient, and the branch wormholes branch in random directions in three-dimensional space, forming a regular three-dimensional tree structure. Moreover, all wormholes maintain dynamic diameter variation characteristics. Furthermore, by obtaining the pressure field slice image of the wormhole network generated by the hydrochloric acid acidification simulation, it was observed that there is an obvious pressure change in the contact area between the wormholes and the matrix, which clearly reflects the acid filtration effect of the wormholes to the matrix. This is highly consistent with the phenomenon of actual three-dimensional acidification experiments, verifying the accuracy of the reservoir acidification simulation results of the method of the present invention in actual scenarios.

[0130] In summary, the method of this invention decouples the acid-rock reaction flow system into a three-dimensional matrix and a one-dimensional wormhole, and combines implicit pressure solution, topological analytical upwind concentration calculation, and dynamic geometric bifurcation evolution. This not only ensures the accuracy of reservoir acidification simulation but also significantly reduces the computational complexity of acidification simulation, achieving accurate simulation of the dominant wormhole growth process. It provides technical support for reservoir acidification design at the reservoir scale and has broad practical engineering application value.

[0131] Of course, the above description is not intended to limit the present invention, and the present invention is not limited to the examples given above. Any changes, modifications, additions or substitutions made by those skilled in the art within the scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A reservoir acidification simulation method at the reservoir scale, characterized in that, Includes the following steps: Step 1: Obtain the acid-rock reaction flow system of the reservoir to be simulated, and decouple it into a three-dimensional matrix continuous domain and a one-dimensional discrete wormhole topology network. Construct the pressure and flow coupling equation, acid transport and reaction equation, and geometric evolution equation to simulate reservoir acidification. Step 2: Perform a fully implicit simultaneous solution to the pressure fields of the three-dimensional matrix continuous domain and the one-dimensional discrete wormhole topology network to obtain the fluid pressure distribution inside the three-dimensional matrix continuous domain and the one-dimensional discrete wormhole topology network in the reservoir, and determine the volumetric flow rate and flow direction inside the one-dimensional discrete wormhole topology network. Step 3: Extract the flow direction of the fluid in the one-dimensional discrete wormhole topology network to construct a directed acyclic graph, perform topological sorting to generate a downstream calculation sequence, and solve the acid concentration at each node in the one-dimensional discrete wormhole topology network based on the first-order implicit upwind scheme along the downstream calculation sequence to obtain the acid concentration field in the one-dimensional discrete wormhole topology network. Step 4: Based on the acid concentration field and geometric evolution equations within the one-dimensional discrete wormhole topology network, explicitly update the geometric parameters of the one-dimensional discrete wormhole topology network, including the radial expansion dimension of each pipe segment and the cumulative extension length of the tip node in the one-dimensional discrete wormhole topology network. Step 5: Determine whether the preset cross-grid growth conditions are met based on the cumulative extension length of the tip node. Perform geometric topology update to generate new tip nodes and wormhole segments, update the time step, and repeat steps 2 to 4 until the maximum acid injection time is reached or the preset acidization breakthrough conditions are met. Then stop the reservoir acidization simulation and obtain the reservoir acidization simulation results.

2. The reservoir acidization simulation method at the reservoir scale according to claim 1, characterized in that, Step 1 includes the following steps: Step 1.1: Construct pressure-flow coupling equations, including three-dimensional matrix seepage equations, one-dimensional wormhole flow equations, and one-dimensional-three-dimensional fluid filtration coupling equations; The fluid in the three-dimensional matrix flow equation is set as a slightly compressible single-phase flow. Combined with the reservoir pressure, the three-dimensional matrix flow equation is determined as follows: ； In the formula, For matrix porosity; The overall compression coefficient; For matrix pressure; For time; For divergence operators; For matrix permeability; The viscosity of the acid solution; For the three-dimensional Hamiltonian gradient operator; The filtration loss per unit length of wormhole into the substrate; Let be the spatial Dirac distribution function; It is a spatial position vector; This refers to the spatial location of the worm hole axis; The fluid in the one-dimensional wormhole flow equation is set to be incompressible, the wormhole is set as a dynamically variable diameter circular pipe, and the fluid flows along the one-dimensional arc length of the wormhole using a Poiseuille flow. The one-dimensional wormhole flow equation is determined as follows: ； In the formula, For pipe flow volumetric flow rate, ; The dynamic radius of the worm pore; The length of the worm hole arc; This refers to the pressure inside the wormhole tube; The one-dimensional-three-dimensional fluid filtration-loss coupling equation is established based on the Peaceman well index model, and its expression is: ； In the formula, The coupling conduction index, ; It is the natural logarithm function; The equivalent discharge radius of the three-dimensional matrix mesh. ,in, , , These represent the spatial dimensions of the three-dimensional matrix mesh. , , Dimensions; Step 1.2: Ignoring longitudinal molecular diffusion and focusing solely on the flow-dominant principle, the acid concentration is calculated within a one-dimensional discrete worm-pore topology network. The acid transport and reaction equations are established as follows: ； In the formula, This refers to the acid concentration. The mass transfer coefficient is . ; The effective molecular diffusion coefficient; It is an asymptotic Sherwood number; The pore Reynolds number; It is a Schmidt number; Step 1.3, based on the law of conservation of mass, convert the acid consumed on the side wall of the wormhole into an increase in the radius of the wormhole segment, resulting in: ； In the formula, The volumetric solubility of the acid solution; The unreacted acid reaching the tip node is described using a tip extension equation, where the tip node is the node at the very front of a one-dimensional discrete wormhole topology. The extension velocity of the unreacted acid advancing through the borehole in front of the tip node is determined, and a geometric evolution equation is established to overcome the grid size limitations in reservoir acidization simulation, enabling the simulation of sub-grid scale dissolution growth. The results are as follows: ； In the formula, The rate of spread of the unreacted acid solution; The cumulative extension length of the tip node; The volumetric flow rate reaching the tip node; The acid concentration required to reach the tip node; The initial germination radius is the radius of the newly generated micropore channels.

3. The reservoir acidification simulation method at the reservoir scale according to claim 1, characterized in that, In step 2, a discretization scheme with implicit pressure and apparent density is used to connect the three-dimensional matrix mesh in the three-dimensional matrix continuous domain with its intersecting one-dimensional wormhole segments. Backward Euler time-difference discretization is then used to construct a large-scale sparse linear equation system. ,in, It is a large sparse coefficient matrix, composed of the coefficients of the three-dimensional matrix mesh discretization equation and the one-dimensional wormhole segment discretization equation. Let be the pressure field vector to be determined, containing The pressure of all 3D matrix meshes and wormhole segments at all times. It is a constant vector; The discretization equation of the three-dimensional matrix mesh is: ； in, ； In the formula, , All are serial numbers of the three-dimensional matrix mesh; For the first The control volume of a three-dimensional matrix mesh; For the first Porosity of a three-dimensional matrix mesh; The overall compression coefficient; For the current time step Pressure values ​​of a three-dimensional matrix mesh; Each is the previous time step. Pressure values ​​of a three-dimensional matrix mesh; The discrete time step; A set of topologically adjacent grids or nodes; The fluid conductivity between adjacent three-dimensional matrix meshes; For the current time step Pressure values ​​of a three-dimensional matrix mesh; , These are all serial numbers of one-dimensional wormhole segments. ; For the first A collection of one-dimensional wormhole segments contained in a three-dimensional matrix mesh; For the first Discrete form well index of a one-dimensional wormhole segment; For the current time step Pressure value of a one-dimensional wormhole segment; For the first The coupling conduction index of a one-dimensional wormhole segment; For the first The physical length of a one-dimensional wormhole segment; For the first The matrix permeability of a three-dimensional matrix grid; The viscosity of the acid solution; The discrete equation for the one-dimensional wormhole segment is: ； in, ； In the formula, The conductivity of the pipe section; For the first One-dimensional wormhole segment and the first The explicit radius between one-dimensional wormhole segments; For the first One-dimensional wormhole segment and the first Center distance between one-dimensional wormhole segments; For the current time step Pressure value of a one-dimensional wormhole segment.

4. The reservoir acidification simulation method at the reservoir scale according to claim 1, characterized in that, Step 3 includes the following steps: Step 3.1: Based on the fluid pressure distribution within the three-dimensional matrix continuum and the one-dimensional discrete wormhole topology network in the reservoir at the current time step, and combined with the volumetric flow rate and flow direction within each one-dimensional wormhole segment in the one-dimensional discrete wormhole topology network, and based on the physical property of the fluid flowing in a tree-like divergent manner, the one-dimensional discrete wormhole topology network is abstracted into a directed acyclic graph. A graph theory algorithm is used to perform topological sorting on the directed acyclic graph to generate a node substitution sequence. In the node substitution sequence, each node is arranged in order of the fluid flow direction from the injection well to the tip of each one-dimensional wormhole segment. Step 3.2, based on the first-order implicit upwind scheme, the discrete mass conservation equation for the one-dimensional wormhole segment in the one-dimensional discrete wormhole topology is determined as follows: ； In the formula, For the node index in the directed acyclic graph; For the previous time step The control volume of the one-dimensional wormhole segment where each node is located; For the current time step Acid concentration at each node; For the previous time step Acid concentration at each node; outflow The set of flow rates for all pipe segments at each node; For outflow boundary flow; For the current time step Acid concentration at each node; Inflow is the first The set of flow rates for all pipe segments at each node; For inflow to the boundary flow; The acid concentration at the upstream node; For the previous time step The radius of the one-dimensional wormhole segment where each node is located; For the first The physical length of the one-dimensional wormhole segment where each node is located; This is the mass transfer coefficient of the previous time step; For the current time step Filtration at each node; Step 3.3: Substitute each node in the generated node sequence into the discrete mass conservation equation of the one-dimensional wormhole pipe segment for solution. During the solution process, the acid concentration of the upstream node corresponding to each node in the directed acyclic graph is pre-extracted. This is then transformed into a known constant term, causing the discrete mass conservation equation for the one-dimensional wormhole pipe segment to degenerate into one containing only a single unknown. The acid concentration at each node in the current directed acyclic graph is calculated using a linear algebraic equation in one variable.

5. The reservoir acidification simulation method at the reservoir scale according to claim 1, characterized in that, In step 5, during the reservoir acidification simulation, the radius of each one-dimensional wormhole segment and the cumulative extension length of each tip node in the one-dimensional discrete wormhole topology are updated at each time step. Then, a geometric topology evolution check is performed on the one-dimensional discrete wormhole topology to determine whether the cumulative extension length of all tip nodes in the one-dimensional discrete wormhole topology has reached a preset growth threshold. If any tip node has reached the growth threshold, new nodes and connected segments are added to the directed acyclic graph, and the cumulative extension distance corresponding to the tip node that has reached the growth threshold is deducted to generate a subgrid, thus completing the macroscopic cross-grid topology update and bifurcation evolution of the one-dimensional discrete wormhole topology.

6. The reservoir acidification simulation method at the reservoir scale according to claim 5, characterized in that, The formula for updating the radius of the one-dimensional wormhole segment in the one-dimensional discrete wormhole topology is as follows: ； In the formula, For the current time step The radius of a one-dimensional wormhole segment; For the previous time step The radius of a one-dimensional wormhole segment; The discrete time step; The volumetric solubility of the acid solution; The mass transfer coefficient of the previous time step; For the current time step The acid concentration within a one-dimensional wormhole segment; For matrix porosity; The formula for calculating the cumulative extension length update of the tip node is as follows: ； In the formula, The cumulative extension length of the current time step's tip node; This is the cumulative extension length of the tip node in the previous time step; The rate of spread of the unreacted acid solution.

Images