A method for evaluating gas drive effect of reservoirs by multi-physics coupling and topology tailoring

Abstract

The application discloses a kind of multi-physical field coupling and topological tailoring's oil reservoir gas displacement effect evaluation method, comprising: obtaining reservoir parameter and constructing the multi-physical field discrete topological space based on graph data structure in memory;State vector is input into physical information neural network agent model to carry out gas-liquid-solid phase equilibrium dimensionality reduction calculation, and output solid phase precipitation volume fraction field and fluid viscosity field;Dynamic effective porosity is calculated in combination with mechanical strain and chemical dissolution;Based on percolation theory, generate dynamic topological tailoring operator, modify flow conductivity attribute in underlying memory to cut off the mathematical connectivity of physical jamming area;Further, based on the updated flow conductivity, assemble full-implicit jacobi matrix and solve iteratively;Finally, exclude isolated dead zone volume, generate evaluation index including dynamic effective sweep efficiency and solid phase precipitation jamming damage index tensor.The application can reduce the risk of matrix singularity divergence, and improve the calculation efficiency and evaluation fidelity of gas displacement effect evaluation.

Description

Technical Field

[0001] This invention relates to the field of reservoir numerical simulation and computer-aided evaluation technology, and in particular to a method for evaluating the gas drive effect of reservoirs using multiphysics coupling and topology trimming. Background Technology

[0002] Existing computer numerical simulation and evaluation methods suffer from significant underlying algorithmic defects when dealing with extremely coupled multi-field systems (e.g., thermo-fluid-mechanical-chemical-phase systems). First, traditional phase equilibrium calculations heavily rely on iterative nested equations. When performing gas-liquid-solid phase equilibrium calculations on large-scale three-dimensional meshes, the computational cost increases exponentially, leading to system overload and extremely long evaluation cycles. Second, existing models often employ static mesh topology and weakly coupled empirical formulas, failing to accurately characterize the dynamic blockage mechanism of microscopic pores and throats caused by solid-phase precipitation at the computer graph theory level. When local porosity approaches its closed extremum due to mechanical compaction and solid-phase precipitation, traditional spatial discretization algorithms still force the calculation of extremely small fluid conductivity values. This results in severe ill-conditioned and singular divergences in the Jacobian matrix during the assembly of underlying algebraic equations, easily causing the collapse of large-scale linear equation systems. Finally, due to the lack of underlying operators that dynamically cut off and isolate physical "blockage dead zones" at the matrix topology level, existing methods are unable to accurately extract the high-dimensional solution space features under fully implicit computation, resulting in serious distortion of the output sweep efficiency and reservoir phase damage evaluation indicators, which cannot provide highly reliable digital decision support for the gas drive development of real heterogeneous reservoirs. Summary of the Invention

[0003] The purpose of this section is to outline some aspects of embodiments of the present invention and to briefly describe some preferred embodiments. Simplifications or omissions may be made in this section, as well as in the abstract and title of this application, to avoid obscuring the purpose of these documents; however, such simplifications or omissions should not be construed as limiting the scope of the invention.

[0004] In view of the aforementioned existing problems, this invention is proposed. Therefore, this invention provides a method for evaluating the gas drive effect of oil reservoirs based on multiphysics coupling and topology trimming, to solve the problems mentioned in the background art.

[0005] To solve the above-mentioned technical problems, the present invention provides the following technical solution: a method for evaluating the effect of reservoir gas drive based on multiphysics coupling and topology trimming, comprising: Obtain the heterogeneous parameters of the reservoir in three-dimensional space, and construct a multiphysics discrete topological space based on graph data structure in computer memory; The multiphysics state vectors of each node in the discrete topological space are input into a pre-configured physical information neural network proxy model to perform dimensionality reduction calculations of gas-liquid-solid phase equilibrium and chemical reaction-dissolution effect, and output the solid precipitation volume fraction field and fluid viscosity field of each node. Based on the solid precipitation volume fraction field, rock mechanical volume strain, and chemical dissolution increment of each node, the dynamic effective porosity of each node is calculated. The dynamic effective porosity of each node is determined based on the critical threshold of percolation theory, and a dynamic topology trimming operator is generated. The conductivity matrix connection attribute of the discrete topology space is modified in memory through the dynamic topology trimming operator. Based on the conductivity matrix after modifying the connection properties, the fully implicit residual vector and Jacobian matrix of the nonlinear multiphysics control equation are assembled, and iteratively solved by a preconditional combination solver to update the multiphysics state vector of each node. Extract the solution results data for all time steps to generate a multidimensional dynamic gas drive performance evaluation index that includes the dynamic effective sweep efficiency and solid phase precipitation blockage damage index tensors.

[0006] As a preferred embodiment of the reservoir gas drive effect evaluation method based on multiphysics coupling and topology trimming described in this invention, wherein: the construction of a multiphysics discrete topology space based on graph data structure in computer memory includes: The three-dimensional spatial discrete grid of the reservoir is transformed into an undirected graph structure. Each grid control volume is mapped to a vertex in the undirected graph, and the flow interface between adjacent grids is mapped to an edge in the undirected graph. Each vertex is assigned a multiphysics state vector and static property parameters, wherein the multiphysics state vector includes pore pressure, phase saturation, global component mole fraction, effective stress tensor, and temperature. Based on the permeability tensor of adjacent vertices and the geometric parameters of the flow interface, the initial fluid conductivity of each edge is calculated using the harmonic averaging method and stored in a sparse matrix.

[0007] As a preferred embodiment of the reservoir gas drive effect evaluation method based on multiphysics coupling and topology trimming described in this invention, the multiphysics state vectors of each node in the discrete topological space are input into a pre-configured physical information neural network surrogate model, including: The input tensor, composed of pore pressure, global component mole fraction, and temperature at each node at the current time step, is input into the physical information neural network proxy model. Through the forward propagation of the physical information neural network proxy model, the output includes the mole fractions of the gas phase, liquid phase, and solid phase, the component fractions within the gas and liquid phases, the phase change-corrected oil phase viscosity, and the solid phase precipitation volume fraction field. When the output of the physical information neural network proxy model of any node is found to be greater than the preset tolerance after being substituted into the thermodynamic residual equation, the phase calculation task of that node is rolled back to the Newton iterative solution process of the state equation.

[0008] As a preferred embodiment of the reservoir gas drive effect evaluation method based on multiphysics coupling and topology trimming described in this invention, the physical information neural network proxy model includes a loss function with embedded physical partial differential residual constraints during the training phase and online inference monitoring. The loss function includes: Data fitting residual term; The thermodynamic fugacity physical constraint residual term is used to calculate the L2 norm of the difference between the fugacity of each component in the liquid phase and the fugacity in the gas phase derived from the equation of state. Additionally, the mole normalization conservation operator term is used to constrain the summation of the predicted differences in mole fractions of gas and liquid phase components to approach a conservation constant.

[0009] As a preferred embodiment of the reservoir gas drive effect evaluation method based on multiphysics coupling and topology trimming described in this invention, the calculation of the dynamic effective porosity of each node based on the solid phase precipitation volume fraction field, rock mechanical volumetric strain, and chemical dissolution increment of each node includes: Obtain the initial porosity of the nodes; Based on the initial porosity, rock compressibility coefficient, and current pore pressure change, the rock mechanical volumetric strain calculated based on effective stress is subtracted to obtain the geomechanical effective stress compressibility porosity. The dynamic effective porosity at the current time step is obtained by adding the chemical dissolution increment of carbonate rocks calculated based on hydrogen ion concentration to the porosity of the effective stress compression effect in geomechanics, and subtracting the solid precipitation volume fraction field output by the physical information neural network surrogate model.

[0010] As a preferred embodiment of the reservoir gas drive effect evaluation method based on multiphysics coupling and topology trimming described in this invention, the method for determining the dynamic effective porosity of each node based on the critical threshold of permeation theory and generating a dynamic topology trimming operator includes: A pre-defined percolation threshold characterizing the loss of connectivity in the micropore-throat network is established. A permeability tensor over-permeability scaling mapping rule is constructed. The difference between the dynamic effective porosity and the over-permeability critical threshold is divided by the difference between the initial porosity and the over-permeability critical threshold. The resulting ratio is then subjected to a power operation of the topological fractal exponent, and the result is used as a tensor scaling operator to perform a tensor product operation with the initial permeability tensor to obtain the dynamic permeability tensor. Specifically, when the dynamic effective porosity of a node is less than or equal to the permeation critical threshold, all elements in the dynamic permeability tensor are forced to be zero.

[0011] As a preferred embodiment of the reservoir gas drive effect evaluation method based on multiphysics coupling and topology trimming described in this invention, wherein: the conductivity matrix connectivity attribute of the discrete topology space is modified in memory by the dynamic topology trimming operator, including: When traversing any edge representing the interface in the undirected graph, determine the dynamic effective porosity of the two adjacent vertices connected by that edge; If the dynamic effective porosity of any vertex is less than or equal to the percolation critical threshold, then the dynamic topology trimming operator containing the step penalty function is triggered. The dynamic topology trimming operator forces the corresponding value of the fluid conductivity containing the edge in the conductivity matrix to be directly assigned to zero, thereby severing the mathematical connectivity between the adjacent vertices in the discrete solution space of the partial differential equation.

[0012] As a preferred embodiment of the reservoir gas drive effect evaluation method based on multiphysics coupling and topology trimming described in this invention, the step of assembling the fully implicit residual vector and Jacobian matrix of the nonlinear multiphysics control equations based on the conductivity matrix after modifying the connectivity properties includes: For each vertex, a discrete residual equation is established that includes the conservation of mass and momentum. The cumulative term of the discrete residual equation is calculated based on the vertex volume, time step, and dynamic effective porosity. The flux term of the discrete residual equation is calculated by accumulating the phase mobility of adjacent vertices, the difference in fluid potential energy, and the fluid conductivity after the dynamic topology trimming operator has been zeroed. The partial derivatives of the discrete residual equations with respect to the multi-physics state vectors are calculated using automatic differentiation techniques, and then assembled to form the Jacobian matrix, which exhibits high-order asymmetric sparse block characteristics.

[0013] As a preferred embodiment of the reservoir gas drive effect evaluation method based on multiphysics coupling and topology trimming described in this invention, the step of iteratively solving using a pre-conditional combined solver to update the multiphysics state vector of each node includes: The Jacobian matrix is ​​mathematically decoupled into a pressure principal matrix exhibiting elliptic partial differential characteristics and a saturation and component principal matrix exhibiting hyperbolic partial differential characteristics using the constrained pressure residual technique. The pressure principal submatrix is ​​preprocessed by applying an algebraic multigrid preconditioner. The preprocessed matrix system is input into a generalized minimum residual solver containing an incomplete LU decomposition preconditioner to perform iterative calculations, solve for the state vector increment, update the multiphysics state vector, until the infinite norm of the residual vector satisfies the preset convergence tolerance threshold.

[0014] As a preferred embodiment of the reservoir gas drive effect evaluation method based on multiphysics coupling and topology trimming described in this invention, the step of extracting the solution results data of all time steps and generating a multidimensional dynamic gas drive effect evaluation index including dynamic effective sweep efficiency and solid phase precipitation blockage damage index tensors includes: When calculating the dynamic effective sweep efficiency, only the sweep volume of vertices whose dynamic effective porosity is greater than the percolation critical threshold is accumulated, excluding the volume of vertices whose connectivity is cut off by the dynamic topology pruning operator. Based on the edge data truncated by the dynamic topology pruning operator, the number of edges whose conductivity is set to zero in the orthogonal three-dimensional spatial directions is counted respectively, and then divided by the global initial total number of edges in the corresponding direction to construct the anisotropic solid-phase precipitation blockage damage index tensor.

[0015] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. By combining graph data structures with statistical physics percolation theory, a dynamic topology trimming operator is constructed. When local pore throats approach their closure extrema due to the precipitation of asphalt or other solid phases or mechanical compaction, this operator directly severs the mathematical connectivity between meshes from the computer's low-level memory. This eliminates the ill-conditioned Jacobian matrix and singularity divergence problems caused by retaining minimal conductivity in traditional numerical simulations, ensuring the stable operation of the fully implicit Newton iterative solver under extreme multi-field coupling conditions.

[0016] 2. To address the computational bottleneck of traditional equation-of-state flash evaporation iteration being time-consuming and prone to crashing, this invention introduces a Physical Information Neural Network (PINN) surrogate model with embedded thermodynamic partial differential residual constraints. This transforms the complex nonlinear root-finding process of gas-liquid-solid phase equilibrium into explicit tensor forward propagation by the computer. Under the premise of strictly adhering to the physical laws of matter conservation and isofagulability, dimensionality reduction and acceleration are achieved, significantly shortening the simulation cycle for gas-driven evaluation.

[0017] 3. Furthermore, this invention also constructs a multidimensional evaluation mechanism based on graph theory-based edge truncation. When calculating sweep efficiency, this mechanism can identify and forcibly remove physically isolated dead zones through topological filtering operators, avoiding artificially inflated sweep volumes. Simultaneously, it reduces the dimensionality of phase damage to an anisotropic second-order tensor, characterizing the degree of permeability damage in different reservoir directions, providing a verifiable and reliable engineering evaluation basis for adjusting water-gas alternation (WAG) or deploying horizontal wells in actual mining operations. Attached Figure Description

[0018] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the description of the embodiments 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. Wherein: Figure 1 This is a flowchart illustrating the overall process of a reservoir gas drive effect evaluation method based on multiphysics coupling and topology trimming, as described in one embodiment of the present invention. Detailed Implementation

[0019] To make the above-mentioned objects, features, and advantages of the present invention more apparent and understandable, specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of them. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the protection scope of the present invention.

[0020] Example 1

[0021] Reference Figure 1 This is the first embodiment of the present invention, which provides a method for evaluating the effect of reservoir gas drive based on multiphysics coupling and topology trimming, including: S1: Obtain the heterogeneous parameters of the reservoir in three-dimensional space and construct a multiphysics discrete topological space based on graph data structure in computer memory.

[0022] In this embodiment, the reservoir's three-dimensional spatial heterogeneity parameters include at least corner grids or unstructured grids, initial porosity, initial absolute permeability tensor, rock compressibility coefficient, relative permeability curve, capillary pressure curve, initial pressure and temperature fields, initial oil, gas and water saturation, fluid PVT component parameters, well location and injection-production regime parameters. These parameters can be derived from well logging interpretation results, core experiments, well test interpretation, PVT test reports, historical production dynamic fitting results, or existing geological modeling results, and are uniformly converted into data objects corresponding to grid nodes or grid edges before being imported into computer memory.

[0023] It should be noted that in traditional multiphase flow numerical simulations, computer memory typically uses three-dimensional multidimensional arrays to store grid data. However, this rigid data structure is prone to matrix singularities when dealing with dynamic discontinuities in micropore throats caused by solid-phase precipitation (such as bituminous precipitation) and ground stress compaction during gas-driven processes, and it is difficult to perform flexible spatial dimensionality reduction calculations. Therefore, this invention breaks with the traditional approach by reconstructing the continuous physical space into a highly flexible graph data structure during the underlying data loading stage. The specific steps are as follows: S101: Topological mapping from a 3D spatial mesh to an undirected graph structure.

[0024] Specifically, the three-dimensional spatial discrete grid data of the reservoir geological model (including unstructured grids or corner grids) is read and converted into an undirected graph structure in computer memory. Each grid control volume is mapped to a set of vertices in an undirected graph. For any vertex It occupies a defined control volume in physical space. Map the flow interfaces that allow fluid exchange between adjacent grids to a set of edges in an undirected graph. For connecting vertices and vertex edge It has a specific interface surface area. and the center distance between the two vertices .

[0025] It should be noted that, through this mapping method, the present invention can transform the solution domain of the three-dimensional partial differential equation into the connection relationship between nodes and edges from the mathematical topology level, thereby providing the underlying data foundation for the dynamic network pruning to be performed.

[0026] S102: Memory allocation and initialization of multiphysics parameters on discrete nodes.

[0027] Furthermore, regarding the multi-field coupling characteristics of "thermal-fluidic-mechanical-phase state", at the time step (i.e., the initial time), for each vertex Assigning dynamic multiphysics state vectors and static rock physical property parameters .

[0028] Furthermore, the dynamic multiphysics state vector is defined as:

[0029] in, Characterizes the pore pressure at this node; , , Characterizing the oil phase and gas phase (e.g., injected) respectively. and the saturation of the aqueous phase; For inclusion The global mole fraction vector of hydrocarbon and non-hydrocarbon components (including) (multi-component systems, from methane to heavy hydrocarbons, etc.), used for phase equilibrium calculations; Characterizes the effective stress tensor of the rock skeleton caused by changes in pore fluid bearing capacity; Characterizes the local reservoir temperature; It is represented as a transpose matrix.

[0030] Furthermore, the static physical property parameters of rocks are defined as follows:

[0031] in, The initial porosity of the node; For the initial absolute permeability tensor (a second-order tensor containing anisotropic information); and These are Young's modulus and Poisson's ratio in rock mechanics, respectively.

[0032] S103: Initial fluid conductivity calculation and sparse matrix pre-assembly based on tensor projection.

[0033] Furthermore, to quantify the initial fluid exchange capacity between adjacent physical nodes, based on the permeability tensor of adjacent vertices and the geometric parameters of the flow interface, for each edge in the undirected graph... Calculate its initial geometric fluid conductivity Because deep reservoirs exhibit strong anisotropy, scalar permeability cannot be directly used. Therefore, this invention employs a spatial normal vector projection combined with the harmonic averaging method for derivation:

[0034] in, The flow interface consists of nodes Pointing to node The unit normal vector. Operator Indicates that the node The absolute permeability tensor is projected onto the effective permeability in the actual flow direction. It is important to emphasize that this invention employs a series harmonic mean (i.e., the product in the numerator and the sum in the denominator) rather than an arithmetic mean because, in extremely heterogeneous reservoirs (such as fracture-matrix interfaces), fluid flux is typically dominated by the side with lower permeability. Therefore, the harmonic mean can guarantee the conservation of mass and continuity of fluid flux across the interface from the discrete scheme of the partial differential equation, thereby avoiding non-physical flow oscillations at points of abrupt property changes.

[0035] Furthermore, once the calculation is complete, all edges will be... Stored in a large sparse matrix in computer memory.

[0036] It should be noted that since the sparse matrix embodies the absolute static physical connectivity of the reservoir during initialization, it can be used as the basic framework for assembling the fully implicit Jacobian matrix in nonlinear iteration.

[0037] S2: Input the multiphysics state vector of each node in the discrete topological space into the pre-configured physical information neural network proxy model, perform dimensionality reduction calculation of gas-liquid-solid phase equilibrium and chemical reaction-dissolution effect, and output the solid precipitation volume fraction field and fluid viscosity field of each node.

[0038] It should be noted that during the injection of carbon dioxide (… During the oil displacement process, supercritical Multiphase and multicomponent mass transfer between crude oil and other substances (such as light hydrocarbon extraction and heavy asphaltene precipitation) exhibits strong nonlinearity. Traditional numerical simulations use real gas equations of state (such as PR-EOS) combined with nonlinear equations for iterative flash vaporization calculations, which not only consumes nearly 70% of the computational resources in the entire simulation process but also easily leads to iterative non-convergence near the critical region. Therefore, this invention introduces a Physical Information Neural Network (PINN) with embedded thermodynamic mechanisms as a surrogate operator for dimensionality reduction calculations, transforming the originally complex phase-state Newton iteration process into an explicit tensor forward propagation at the underlying computer level. The process is as follows: S201: Constructing the topology and input-output mapping of a physical information neural network proxy model.

[0039] Furthermore, a deep neural network topology containing a multilayer perceptron (MLP) is constructed in computer memory. The input layer of the network is set to the current time step (…). The state data of each node are combined and normalized to form the input tensor. :

[0040] It should be noted that, from a physical perspective, based on the Gibbs phase law, the three state variables selected in this invention are sufficient to uniquely determine the thermodynamic equilibrium state of a multi-component system.

[0041] Furthermore, the network's output layers are configured to output prediction tensors via forward propagation. :

[0042] in, These represent the total mole fractions of the gas phase, liquid phase, and precipitated solid phase (such as asphaltene), respectively. and These represent the mole fraction vectors of each component within the gas and liquid phases, respectively. This is the oil phase viscosity after phase change mass transfer correction (reflecting the viscosity surge effect caused by the vaporization of light components). The solid-phase precipitation volume fraction field serves as the core input for evaluating the dynamic blockage of microscopic pores and throats.

[0043] Furthermore, in this embodiment, the specific topology and data preprocessing process of the physical information neural network proxy model are as follows: First, an offline training dataset was constructed. Using the Real Gas Equation of State (PR-EOS), thermodynamic equilibrium sample pairs covering CO2, methane, and heavy hydrocarbon components were generated within a preset temperature range of 300K to 450K, a pressure range of 1MPa to 80MPa, and the target oil sample composition range, combined with Latin hypercube sampling. Sample labels included phase fraction, mole fraction of each phase component, oil phase viscosity, and solid precipitation volume fraction of asphaltenes or heavy components. Second, input data normalization was performed. Logarithmic normalization was applied to pressure, and range normalization was applied to temperature and component mole fraction, mapping the input tensor to the [-1,1] interval. Next, a multilayer perceptron network was constructed. The multilayer perceptron includes one input layer, 5 to 8 hidden layers (each hidden layer containing 128 to 256 neurons), and one output layer. The hidden layers used Swish or Tanh activation functions, and the output layer used Softmax constraints on phase fraction and component fraction. Finally, during the offline training phase, Adam and L-BFGS were combined for optimization, with the initial learning rate of Adam set to 1×10⁻⁶. -3 Up to 1×10 -4 When the validation set loss no longer decreases after a preset number of consecutive rounds, switch to L-BFGS optimization until the total loss function is less than 1×10. -5 Alternatively, verify that the average error of the phase separation rate of the verification set is less than 1%, and complete the solidification configuration of the model parameters.

[0044] S202: Offline training and loss function configuration of models based on multi-objective partial differential residual constraints.

[0045] Furthermore, purely data-driven neural networks are prone to errors that violate physical laws in reservoir geological extrapolation scenarios. Therefore, this invention configures a total loss function with embedded physical partial differential residual constraints during the training phase and online inference monitoring. Joint optimization by automatic differential computation graph:

[0046] in, These are dynamic adaptive weighting coefficients, with values ​​ranging from [0,1]. .

[0047] Specifically, for the data fitting residual term in the above formula... This item is mainly used to measure the mean square error between the predicted value and the sample label set pre-generated based on the high-precision state equation.

[0048] Specifically, for the thermodynamic isofugacity physical constraint residual term in the above formula... This is mainly based on the thermodynamic phase equilibrium principle. When a multi-component system reaches equilibrium, any microscopic fluid component... The chemical potentials in the gas and liquid phases must be equal, i.e., their fugacity must be equal. In this embodiment, this can be constructed as a penalty term for the network layer:

[0049] in, and Components The fugacity in both the gas and liquid phases can be explicitly expressed by the partial derivatives of the equation of state.

[0050] It should be noted that this constraint aims to force the neural network to actively follow the underlying thermodynamic evolution boundary of the fluid phase during the weight optimization process, so as to ensure that the output phase distribution conforms to physical laws.

[0051] Specifically, for the mole normalization conservation operator in the above formula... This term is mainly used to constrain the conservation of matter, that is:

[0052] It is important to note that the sum of the mole fractions of each component within the gas-liquid phase predicted by this constraint needs to approach the conservation constant 1 in order to avoid mass loss during dimensionality reduction mapping.

[0053] S203: Perform online inference output and extract features of chemical reactions and dissolution effects.

[0054] Specifically, at each time step of the nonlinear numerical simulation, tensor computation cores deployed on GPUs / NPUs are used to compute the data from millions of discrete nodes across the entire region. Batch input into the PINN model with the physical mechanism already embedded is performed for fast forward propagation to directly extract the solid precipitation volume fraction field of each node. With fluid viscosity field .

[0055] It should be noted that this step aims to skip the nonlinear root-finding process of the traditional Rachford-Rice equations, thereby improving computational efficiency at the level of multiphysics dimensionality reduction.

[0056] S204: Thermodynamic fault tolerance and rollback mechanism based on phase equilibrium residual tolerance determination.

[0057] Furthermore, since the computational conditions of complex heterogeneous reservoirs are extremely variable, this invention also adds a double insurance at the underlying algorithm level to ensure the numerical robustness of the digital simulation.

[0058] Specifically, the prediction tensor output by the PINN surrogate model is substituted into the thermodynamic isofugacity residual equation to calculate the isofugacity residual value of each node in real time. When the isofugacity residual value of any node exceeds a preset physical tolerance, the neural network prediction of that node is deemed to have a risk of distortion; the preset physical tolerance can be set to 1×10 based on the PVT experimental fitting accuracy. -6 Up to 1×10 -4 At this point, the computer triggers the backoff operator, discards the PINN prediction result of the node, and backoffs the local phase state calculation task of the node to the Newton iterative solution process of the PR-EOS or SRK-EOS equation of state until the phase fraction, component mole fraction and fugacity residual all meet the convergence conditions.

[0059] It should be noted that by using a hybrid computing power scheduling strategy that combines dimensionality reduction as the primary approach with traditional iterative fallback, we can ensure the numerical convergence of the global solution of partial differential equations while significantly accelerating the overall simulation.

[0060] S3: Calculate the dynamic effective porosity of each node based on the solid precipitation volume fraction field, rock mechanical volumetric strain, and chemical dissolution increment.

[0061] It should be noted that in traditional commercial reservoir numerical simulators, porosity is often simplified to a single-variable linear function of pore pressure (i.e., considering only the rock compressibility coefficient). However, in real-world carbon dioxide (CO2) simulations... During oil displacement, the evolution of reservoir pore space is an extreme multi-field coupling (THMC) process characterized by intense competition among mechanical compaction closure, acidic chemical dissolution expansion, and heavy solid precipitation blockage. The absence of any single physical field mechanism will lead to severe distortion in the prediction of microscopic pore throat evolution trends. Therefore, this invention constructs a dynamic effective porosity evolution equation with multi-physics joint correction. The specific implementation steps are as follows: S301: Couple fluid pressure variation with effective stress in the skeleton to calculate porosity due to effective stress compression effect in geomechanics.

[0062] Specifically, first, obtain the nodes. initial porosity Extract the current time step Pore ​​pressure at the next node And the volumetric strain of rock mechanics obtained by solving the three-dimensional rock mechanics equilibrium equations. Subsequently, based on the theory of porous elastic media mechanics, the porosity of the effective stress compression effect in geomechanics was calculated. :

[0063] in, is the equivalent compression coefficient of the rock skeleton; The initial pore pressure at the node; is the Biot effective stress coefficient, which characterizes the contribution weight of fluid pressure to the overall deformation of the rock; This represents the volumetric strain of the rock at the current time step relative to the initial state. This value is positive when the rock volume shrinks under pressure.

[0064] It should be noted that, through the above formula, this invention can decouple the mechanical response of porosity into two parts. The first part is... The first part represents the expansion effect of pore pressure itself (where increases or decreases in fluid pressure directly cause pore expansion or contraction); the second part is... This represents the squeezing effect caused by the compaction and deformation of the rock skeleton due to changes in the macroscopic geostress field. Through a dual mechanical constraint method, it is possible to accurately capture the severe compaction and closure phenomenon of deep rocks caused by reservoir pressure depletion in the later stages of gas drive.

[0065] S302: Calculate the incremental chemical dissolution of carbonate rocks based on fluidization reaction kinetics.

[0066] Specifically, due to After being injected into the reservoir, it dissolves in the formation water to form carbonic acid, releasing a large amount of hydrogen ions. This process, in turn, causes intense dissolution of the reservoir matrix (especially carbonate reservoirs), leading to an irreversible increase in porosity. Therefore, this invention, based on non-equilibrium chemical reaction kinetics, calculates the node at the current time step. Chemical dissolution increment :

[0067] in, This represents the current simulation time step. For dissolved rock minerals (such as calcite) The molar volume of ). Specific surface area of ​​porous media, representing the effective contact area for chemical reactions; Let be the chemical reaction rate constant under these temperature and pressure conditions; The concentration of hydrogen ions in the aqueous phase at the current time step is calculated from phase equilibrium. Solubility is determined dynamically; denoted as the reaction order.

[0068] It should be noted that by calculating the chemical dissolution increment of carbonate rocks, the expansion effect of the acidic fluid front propagation in porous media on the micropore throat can be quantified from the perspective of chemical reaction kinetics. That is, in the near-wellbore zone of the injection well, due to continuous contact with fresh acidic fluid... fluid, Extremely high, therefore Rapid accumulation will significantly improve local fluid conduction capacity.

[0069] S303: Integrates multiphysics variables to solve for the dynamic effective porosity at the current time step.

[0070] Furthermore, the porosity obtained from the aforementioned S301 due to effective geomechanical stress compression effect is... Chemical dissolution increment obtained from S302 And the solid precipitation volume fraction field directly output by the forward propagation of the Physical Information Neural Network (PINN) surrogate model in S2. (The volume percentage representing asphaltene or heavy hydrocarbons precipitated in solid form and occupying the pores) can be algebraically superimposed to obtain the node at the current time step. The final dynamic effective porosity :

[0071] Furthermore, by fully expanding it, we can obtain the multiphysics coupled porosity master equation in the present invention:

[0072] It should be noted that the master equation reflects the complex physical competition mechanism of the reservoir during the gas drive process. Mechanical compaction and solid precipitation together act as a negative feedback mechanism, reducing porosity and potentially blocking microscopic pore throats; pressure support and chemical dissolution act as a positive feedback mechanism, expanding porosity and improving local conductivity. When the combined effect of these multiple physical effects causes the local dynamic effective porosity to approach the permeability critical threshold, it will trigger the topology trimming operation of the bottom layer operator.

[0073] S4: Based on the critical threshold of percolation theory, determine the dynamic effective porosity of each node and generate a dynamic topology trimming operator. Modify the conductivity matrix connectivity property of the discrete topology space in memory through the dynamic topology trimming operator.

[0074] It should be noted that in traditional reservoir numerical simulators (such as commercial software based on the finite difference method), the change in permeability with porosity is usually fitted using empirical continuous functions (such as the classic Kozeny-Carman formula). However, during gas injection, due to severe solid phase precipitation (asphaltite blockage) and in-situ stress compaction, the porosity of local nodes approaches its minimum value, and traditional formulas calculate extremely small but non-zero permeability (such as millidarcy). This not only physically violates the objective fact that the pore throat is completely broken and disconnected, but also has adverse effects on the underlying computer algorithm. That is, the extremely small permeability will cause the condition number of the large sparse Jacobian matrix to increase sharply during the subsequent discretization and assembly of partial differential equations, leading to serious matrix ill-conditioning and singularity divergence, ultimately causing the Newton iteration solution to collapse directly. To address this, this invention introduces percolation theory from statistical physics and combines it with computer graph theory algorithms to design a dynamic topology pruning operator that can achieve physical disconnection and mathematical decoupling. The implementation process is as follows: S401: Constructing a dynamic permeability tensor evolution model based on percolation theory.

[0075] Furthermore, a permeability critical threshold is preset to characterize the loss of connectivity in the microscopic pore-throat network. This permeability critical threshold can be determined based on core mercury intrusion porosimetry, nuclear magnetic resonance pore size distribution, porosity-permeability curves, or digital core pore-throat network simulation results. In the absence of measured data, an initial value can be selected within a porosity range of 0.02 to 0.05, and then corrected using core displacement experiments or historical production dynamics. Physically speaking, when the dynamic effective porosity decreases to near the permeability critical threshold, the fluid struggles to form continuous macroscopic flow channels in three-dimensional space, resulting in a loss of permeability connectivity.

[0076] Furthermore, obtain the node at the current time step calculated in S3. Dynamic effective porosity Construct the percolation scaling mapping rule for the percolation tensor and compute the node... dynamic permeability tensor :

[0077] in, This is the initial absolute permeability tensor (second-order tensor) of this node. This represents the initial porosity of the node. This is the preset percolation threshold; is the topological fractal index, used to characterize the geometric fractal features of the microscopic pore throat connectivity of porous media, and its value is generally between 2.0 and 4.0; 0 is a second-order tensor with all zeros.

[0078] It should be noted that, through the above formula, the present invention is able to... This normalized effective flow porosity ratio, after being increased in dimensionality through exponentiation, is used as a tensor scaling operator for tensor multiplication with the initial tensor. In other words, when... In this case, the algorithm can force all elements in the dynamic permeability tensor to be set to zero. Based on this, the design not only macroscopically recreates the nonlinear seepage barrier caused by pore throat closure, but also provides mathematical triggering conditions for generating digital topology trimming operators.

[0079] To avoid affecting implementation due to unclear threshold selection, the topological fractal index can be determined by fitting the core pore-permeability experimental curve, simulating the digital core pore-throat network, or fitting historical production dynamics. If experimental data is lacking for the target block, an initial value can be selected in the range of 2.0 to 4.0, and the fitting error of the gas content breakthrough time, oil production, injection-production pressure difference, or pressure recovery curve can be used as the correction target.

[0080] S402: Generate a dynamic topology pruning operator that includes a step penalty function.

[0081] Specifically, porosity is determined by mechanical compaction (reducing porosity), solid precipitation (reducing porosity), and chemical dissolution (increasing porosity). If a traditional step function is used, assuming that in the first time step, solid precipitation causes porosity to fall below a threshold, the edge is cut off; however, in the second time step, strong acidic dissolution causes the porosity of that node to exceed the threshold again, and the edge recovers. But in practical applications, after asphaltene and other solid phases block the pore throat, even with slight expansion, the original seepage channels are difficult to recover immediately (permeability hysteresis effect); and in fully implicit solutions, frequent edge disconnections and reconnections cause the Jacobian matrix structure to oscillate violently, directly leading to non-convergence. To achieve dead zone cutting at the underlying data structure, this invention constructs a dynamic topology pruning operator for undirected graph edge connections in memory. Furthermore, this operator incorporates a step penalty function based on the physical states of neighboring nodes:

[0082] in, and These represent the undirected graph formed by any edge. The two adjacent vertices connected; This indicates taking the extreme value with smaller porosity between two adjacent vertices, reflecting the flow bottleneck effect; This is represented as a step penalty function, which depends not only on the porosity of the current time step but also on the topological state of the previous time step. The specific decision logic is as follows: If the edge was pruned in the previous time step (i.e., its state was 0), and the expansion increment in the current time step fails to exceed the preset hysteresis recovery threshold, the current value is forcibly kept at 0; a binary mutation from 1 to 0 only occurs when a blockage or disconnection occurs for the first time. This avoids frequent disconnection and reconnection of edges due to small fluctuations in porosity near the threshold during nonlinear iteration, ensuring the stability of the non-zero element structure of the Jacobian matrix.

[0083] It should be noted that if the dynamic effective porosity of any adjacent vertex or any vertex is less than or equal to the critical threshold, the value will undergo a binary abrupt change, switching from 1 to 0.

[0084] S403: Perform dynamic reconstruction and disconnection of the conductance matrix in a discrete topological space in memory.

[0085] Furthermore, the large sparse conductivity matrix initialized in S1 is retrieved at the computer's lower level. A graph-based traversal algorithm is then used to traverse any edge representing the interface in the undirected graph. .

[0086] Furthermore, during the traversal, the initial geometric fluid conductivity at the current time step is... (or based on) The conductivity calculated by re-harmonizing the average (and the dynamic topology pruning operator mentioned above) Perform scalar multiplication to modify the connectivity properties of the conductivity matrix in the discrete topological space:

[0087] It should be noted that, through the above steps, if a percolation disconnection is determined to have occurred, the topology pruning operator is set to 0, and the computer sets the corresponding value of the fluid conductivity containing that edge in the conductivity matrix to zero. From a mathematical topology perspective, this is equivalent to pruning the edge between adjacent vertices in the memory graph data structure. Through this algebraic operation, the mathematical connectivity between adjacent vertices can be severed in the discrete solution space of partial differential equations, avoiding the generation of minimal conductivity when assembling flux terms later. This reduces the risk of singularities or oscillatory non-convergence in the nonlinear Newton iterative equation system, and improves the fault tolerance and computational efficiency of multi-field coupled simulations.

[0088] S5: Based on the conductivity matrix after modifying the connection properties, assemble the fully implicit residual vector and Jacobian matrix of the nonlinear multiphysics control equation, and iteratively solve it through a preconditional combination solver to update the multiphysics state vector of each node.

[0089] It should be noted that the multiphysics coupling of heat, fluid, mechanical, chemical, and phase states induced by carbon dioxide gas drive is a highly nonlinear process. Traditional explicit or implicit pressure-explicit saturation solution schemes are easily subject to Courant-Friedrich-Lévy (CFL) conditions, leading to computational divergence. Therefore, this invention employs a fully implicit method (FIM) for time and space discretization. Because this invention uses a dynamic topology pruning operator in S4 to sever the mathematical connections corresponding to the blocked region, it can reduce the risk of ill-conditioned or singular Jacobian matrices when the orifice throat is completely blocked, thereby ensuring the assembly and solution of large-scale nonlinear equation systems. The specific processing procedure is as follows: S501: Construct the cumulative and flux terms of the discrete residual equations for multiphysics.

[0090] Furthermore, for each vertex (node) in the undirected graph space... and every fluid component in the system Based on the laws of conservation of mass and momentum, a fully implicit discrete residual equation based on the control volume is established. :

[0091] in, Source and sink terms (representing injection or production from a well). For the cumulative terms of the discrete residual equation, This is the flux term in the discrete residual equation.

[0092] Furthermore, the cumulative terms of this discrete residual equation Based on vertex volume, time step, and dynamic effective porosity obtained in S3. The calculation yielded:

[0093] in, For nodes The mesh volume; For time step; Represents the fluid phase (oil phase, gas phase, water phase); Phase density; Phase saturation; Components exist The mole fraction in the phase. Superscript and These represent the latest time step and the previous time step in the current nonlinear iteration, respectively.

[0094] Furthermore, the flux term of this discrete residual equation Based on the phase mobility of adjacent vertices, the difference in fluid potential energy, and the fluid conductivity after being zeroed out by the dynamic topology trimming operator in S4. Cumulative calculation:

[0095] in, For nodes The set of all adjacent nodes; The phase mobility between adjacent nodes is (relative permeability divided by phase viscosity, where phase viscosity is directly output by forward propagation of the PINN surrogate model of S2). It is the fluid phase potential energy (including the difference between pore pressure and gravitational potential).

[0096] It should be noted that, due to After the step penalty function operator filtering, when a local pore throat becomes completely blocked due to asphalt precipitation or mechanical compaction, It is strictly equal to 0. Therefore, when the computer performs flux accumulation, it will directly skip the disconnected edge (i.e., no fluid flux exchange), thus avoiding the accumulation of spurious flow and rounding errors caused by tiny conductivity from a physical source.

[0097] S502: Calculate partial derivatives and assemble higher-order Jacobian matrices using automatic differentiation techniques.

[0098] Furthermore, within the Newton-Raphson nonlinear iterative solution framework, it is necessary to calculate the partial derivatives of the residual vector with respect to all unknown state variables. Given that this invention involves extreme coupling between phase transition, mechanics, and chemistry, traditional analytical differentiation is almost impossible. Therefore, this invention embeds forward-mode automatic differentiation (AD) technology into the in-memory solver.

[0099] Specifically, using a dual algebraic system, the input multiphysics state vector is... This is converted to a dual form with perturbation terms. Next, the discrete residual equation is calculated. Simultaneously, it automatically carries and outputs precise partial derivatives on the underlying computational graph, assembling them to form a Jacobian matrix exhibiting high-order asymmetric sparse block characteristics. :

[0100] It should be noted that due to the strong convection term in the seepage process (requiring the use of an upwind scheme), the assembled Jacobian matrix has a high degree of asymmetry; and due to the presence of multiple field variables, the dimension of each grid submatrix block is extremely high.

[0101] S503: Jacobian matrix mathematical decoupling and preprocessing based on constrained pressure residual (CPR) technology.

[0102] Furthermore, due to the extremely rapid propagation speed of pore pressure (exhibiting elliptic partial differential equation characteristics) and the relatively slow phase saturation and component migration speeds (exhibiting hyperbolic partial differential equation characteristics), directly solving the overall Jacobian matrix of the above assembly is extremely difficult. Therefore, we utilize the Constrained Pressure Residual (CPR) technique to mathematically decouple the Jacobian matrix in two stages.

[0103] Specifically, firstly, a pressure-dependent principal submatrix is ​​separated from the full Jacobian matrix using an algebraic transformation matrix. Due to the globally strong coupling of the pressure field, this pressure principal submatrix is ​​preprocessed using an Algebraic Multigrid (AMG) preconditioner. It should be noted that AMG, through constraints and interpolation operations between different grid coarseness levels, can significantly attenuate high-frequency and low-frequency error components in the pressure field.

[0104] S504: Perform generalized minimum residual method iteration and state update through combined solvers.

[0105] Furthermore, the matrix system, after the two-stage decoupling preprocessing (at which point the condition number of the matrix has been significantly reduced), is input into the Generalized Minimum Residual Method (GMRES) solver, which includes an incomplete LU Factorization (ILU(k)) preconditioner, to perform Krylov subspace iterative operations. The state vector increment for the current Newton iteration step is then calculated. And update the multiphysics state vector:

[0106] Furthermore, the updated version Substitute into S501 to calculate the new residual vector Repeat steps S501 to S504 until the residual vector reaches its infinite norm. (That is, the maximum absolute value of the residuals of all nodes and all equations in the system) satisfies the preset convergence tolerance threshold (e.g.) If the nonlinear multi-field coupled equations of the current time step converge, the state variables are saved, and the simulation proceeds to the next time step.

[0107] S6: Extract the solution results data for all time steps and generate a multidimensional dynamic gas drive performance evaluation index that includes the dynamic effective sweep efficiency and solid phase precipitation blockage damage index tensors.

[0108] It should be noted that traditional reservoir numerical simulators, when outputting evaluation indicators such as sweep efficiency, typically only perform a simple scalar summation of the gas saturation and pore volume across the entire grid. However, in actual oil injection and displacement processes, once severe phase solid precipitation (e.g., bituminous blockage) or extreme mechanical closure occurs in local pore throats, the fluid in these areas will essentially lose its physical flow capacity, becoming isolated dead zones. Traditional statistical algorithms, unable to identify network disconnections, will still include the residual fluid in these dead zones in the sweep volume, resulting in a significantly inflated final gas drive efficiency assessment. Therefore, this invention, based on the aforementioned dynamic trimming records of graph data structures, introduces a multidimensional dynamic evaluation method based on topological connectivity. The specific steps of this method are as follows: S601: Extracting multidimensional physical field solution sets and topological map features with time series data from computer memory.

[0109] Specifically, after completing the numerical solution for the entire lifecycle (from the initial moment to the set end of gas-driven extraction), each time step is extracted according to the set output frequency. The state data of each vertex (mesh node) under the substructure (such as dynamic effective porosity) Gas phase saturation ) and the connectivity identifier for each edge in the undirected graph (i.e., the dynamic topology pruning operator generated in S4). (Current value).

[0110] S602: Combine topological filtering operator and saturation field to calculate dynamic effective sweep efficiency.

[0111] Furthermore, a Boolean conditional decision mechanism based on graph theory vertices is introduced when calculating the dynamic effective sweep efficiency. That is, only the sweep volume of vertices whose dynamic effective porosity is greater than the percolation critical threshold is accumulated, and the volume of vertices whose connectivity is severed by the dynamic topology pruning operator is forcibly excluded in the underlying algorithm. The underlying discrete calculation formula is as follows:

[0112] in, For the first The dynamic effective sweep efficiency at each time step, a dimensionless quantity; The total number of vertices in an undirected graph. For nodes The volume of the rock is controlled by its volume; and Representing nodes respectively In the Dynamic effective porosity and initial porosity at each time step; and These represent the gas phase saturation at the current time step and the initial gas phase saturation (characterizing the incremental displacement due to gas injection), respectively. The initial bound water saturation is represented by the denominator, which represents the total effective pore volume available for hydrocarbons to occupy in the initial state of the reservoir. A topological filtering operator to characterize isolated dead zones.

[0113] Furthermore, when the node of ( When the percolation threshold is 1, and the percolation threshold is 0 ... then the percolation threshold is 0. If all edges of a node are forced to zero by the dynamic pruning operator in S4 (i.e., physically completely blocked and isolated), then .

[0114] It should be noted that when the injected gas front triggers severe bituminous precipitation, causing some grid nodes that originally contained gas to be completely surrounded and cut off by the surrounding blockage zone, the operator will collapse to 0. At this time, the volume contribution of that node in the numerator is cleared to zero. Algorithmically, this means eliminating those visible but unminable spurious swept volumes, providing the most conservative numerical basis for predicting the actual production of the mine.

[0115] S603: Construct an anisotropic solid-phase precipitation blockage damage index tensor based on graph theory-based edge-truncated dimension reduction mapping.

[0116] Furthermore, reservoir damage is often not uniform. Due to differences in sedimentary facies, the initial permeability of the formation differs significantly in the horizontal and vertical directions. Therefore, the plugging behavior caused by supercritical CO2 extraction and solid-phase precipitation exhibits anisotropy in three-dimensional space. To address this, this step, based on the edge data truncated by the dynamic topology pruning operator in the bottom layer, counts the number of edges with zero conductivity in orthogonal three-dimensional directions, and divides them by the total initial global edge count in the corresponding direction to construct the anisotropic solid-phase precipitation plugging damage index tensor:

[0117] Furthermore, for any directional component on the main diagonal of the matrix (in... Direction For example, its calculation formula is:

[0118] in, Indicates along The set of all connecting edges arranged along the dominant axis direction; For the first At time step, applied to the edge Dynamic topology pruning operator state values ​​(when connected) When leakage and disconnection occur due to blockage ); For simulation initialization (in S1) along The total number of edges whose initial fluid conductivity was successfully assigned in the direction.

[0119] It should be noted that traditional evaluations often provide a scalar penetration rate decrease percentage for the entire region, which is a vague assessment and lacks any guiding significance. Therefore, this invention utilizes the directional properties of edges in a graph data structure to abstract the damage into a second-order diagonal tensor. For example, if the output displays… This represents supercritical. The precipitation of heavy components caused by extraction mainly seals off the fluid interfaces between vertical layers. Therefore, this tensor index can guide engineers to adjust subsequent production plans in a timely manner (e.g., stop implementing water-gas alternating injection (WAG) and switch to horizontal wells to avoid the risk of loss of vertical interlayer connectivity).

[0120] S604: The integrated data structure outputs a comprehensive gas-driven decision map.

[0121] After outputting the multidimensional dynamic gas drive effect evaluation index, the consistency of the dynamic effective sweep efficiency and solid phase precipitation blockage damage index tensor with the field production dynamics can be verified. The field production dynamics include injection pressure, gas cut at the production end, gas channeling breakthrough time, daily oil production, daily water production, and inter-well tracer response data. When the deviation between the evaluation index and the production dynamics exceeds the preset verification threshold, the permeability critical threshold, topological fractal index, relative permeability curve, or the input sample range of the PINN surrogate model are recalibrated to ensure that the evaluation results can be repeatedly implemented and verified by those skilled in the art.

[0122] Furthermore, the computer extracts the aforementioned dimensionality reduction data. (Scalar evolution curve) and The tensor evolution matrix, along with the component mobility field data output by the Physical Information Neural Network (PINN), are serialized into a structured data file such as JSON or HDF5. Finally, this file is rendered and output as a multidimensional dynamic gas drive performance evaluation index panel.

[0123] It should be noted that by opening up the underlying data flow of the algorithm, this invention enables technicians (reservoir engineers) to intuitively observe how micro-phase phase transitions are disconnected through graph theory topology, and ultimately quantify the entire evolution process that affects macro-sweeping efficiency.

[0124] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A method for evaluating the effect of gas drive in oil reservoirs using multiphysics coupling and topology trimming, characterized in that, include: Obtain the heterogeneous parameters of the reservoir in three-dimensional space, and construct a multiphysics discrete topological space based on graph data structure in computer memory; The multiphysics state vectors of each node in the discrete topological space are input into a pre-configured physical information neural network proxy model to perform dimensionality reduction calculations of gas-liquid-solid phase equilibrium and chemical reaction-dissolution effect, and output the solid precipitation volume fraction field and fluid viscosity field of each node. Based on the solid precipitation volume fraction field, rock mechanical volume strain, and chemical dissolution increment of each node, the dynamic effective porosity of each node is calculated. The dynamic effective porosity of each node is determined based on the critical threshold of percolation theory, and a dynamic topology trimming operator is generated. The conductivity matrix connection attribute of the discrete topology space is modified in memory through the dynamic topology trimming operator. Based on the conductivity matrix after modifying the connection properties, the fully implicit residual vector and Jacobian matrix of the nonlinear multiphysics control equation are assembled, and iteratively solved by a preconditional combination solver to update the multiphysics state vector of each node. Extract the solution results data for all time steps to generate a multidimensional dynamic gas drive performance evaluation index that includes the dynamic effective sweep efficiency and solid phase precipitation blockage damage index tensors.

2. The method for evaluating the reservoir gas drive effect based on multiphysics coupling and topology trimming as described in claim 1, characterized in that, The construction of a multiphysics discrete topological space based on graph data structures in computer memory includes: The three-dimensional spatial discrete grid of the reservoir is transformed into an undirected graph structure. Each grid control volume is mapped to a vertex in the undirected graph, and the flow interface between adjacent grids is mapped to an edge in the undirected graph. Each vertex is assigned a multiphysics state vector and static property parameters, wherein the multiphysics state vector includes pore pressure, phase saturation, global component mole fraction, effective stress tensor, and temperature. Based on the permeability tensor of adjacent vertices and the geometric parameters of the flow interface, the initial fluid conductivity of each edge is calculated using the harmonic averaging method and stored in a sparse matrix.

3. The method for evaluating the reservoir gas drive effect based on multiphysics coupling and topology trimming as described in claim 2, characterized in that, The multiphysics state vectors of each node in the discrete topological space are input into a pre-configured physical information neural network surrogate model, including: The input tensor, composed of pore pressure, global component mole fraction, and temperature at each node at the current time step, is input into the physical information neural network proxy model. Through the forward propagation of the physical information neural network proxy model, the output includes the mole fractions of the gas phase, liquid phase, and solid phase, the component fractions within the gas and liquid phases, the phase change-corrected oil phase viscosity, and the solid phase precipitation volume fraction field. When the residual value obtained by substituting the output of the physical information neural network proxy model of any node into the thermodynamic residual equation is greater than the preset tolerance, the phase calculation task of that node is rolled back to the Newton iterative solution process of the state equation.

4. The method for evaluating the reservoir gas drive effect based on multiphysics coupling and topology trimming as described in claim 3, characterized in that, The physical information neural network proxy model includes a loss function with embedded physical partial differential residual constraints during the training phase and online inference monitoring. This loss function includes: Data fitting residual term; The thermodynamic fugacity physical constraint residual term is used to calculate the L2 norm of the difference between the fugacity of each component in the liquid phase and the fugacity in the gas phase derived from the equation of state. Additionally, the mole normalization conservation operator term is used to constrain the summation of the predicted differences in mole fractions of gas and liquid phase components to approach a conservation constant.

5. The method for evaluating the reservoir gas drive effect based on multiphysics coupling and topology trimming as described in claim 3, characterized in that, The dynamic effective porosity of each node is calculated based on the solid precipitation volume fraction field, rock mechanical volumetric strain, and chemical dissolution increment, including: Obtain the initial porosity of the nodes; Based on the initial porosity, rock compressibility coefficient, and current pore pressure change, the rock mechanical volumetric strain calculated based on effective stress is subtracted to obtain the geomechanical effective stress compressibility porosity. The dynamic effective porosity at the current time step is obtained by adding the chemical dissolution increment of carbonate rocks calculated based on hydrogen ion concentration to the porosity of the effective stress compression effect in geomechanics, and subtracting the solid precipitation volume fraction field output by the physical information neural network surrogate model.

6. The method for evaluating the reservoir gas drive effect based on multiphysics coupling and topology trimming as described in claim 5, characterized in that, The determination of the dynamic effective porosity of each node based on the critical threshold of percolation theory, and the generation of a dynamic topology trimming operator, includes: A pre-defined percolation threshold characterizing the loss of connectivity in the micropore-throat network is established. A permeability tensor over-permeability scaling mapping rule is constructed. The difference between the dynamic effective porosity and the over-permeability critical threshold is divided by the difference between the initial porosity and the over-permeability critical threshold. The resulting ratio is then subjected to a power operation of the topological fractal exponent, and the result is used as a tensor scaling operator to perform a tensor product operation with the initial permeability tensor to obtain the dynamic permeability tensor. Specifically, when the dynamic effective porosity of a node is less than or equal to the permeation critical threshold, all elements in the dynamic permeability tensor are forced to be zero.

7. The method for evaluating the reservoir gas drive effect based on multiphysics coupling and topology trimming as described in claim 6, characterized in that, Modifying the connectivity properties of the conductance matrix of the discrete topological space in memory using the dynamic topology pruning operator includes: When traversing any edge representing the interface in the undirected graph, determine the dynamic effective porosity of the two adjacent vertices connected by that edge; If the dynamic effective porosity of any vertex is less than or equal to the percolation critical threshold, then the dynamic topology pruning operator containing the step penalty function is triggered. The dynamic topology trimming operator forces the corresponding value of the fluid conductivity containing the edge in the conductivity matrix to be directly assigned to zero, thereby severing the mathematical connectivity between the adjacent vertices in the discrete solution space of the partial differential equation.

8. The method for evaluating the reservoir gas drive effect based on multiphysics coupling and topology trimming as described in claim 7, characterized in that, The assembly of the fully implicit residual vector and Jacobian matrix of the nonlinear multiphysics control equations based on the conductance matrix after modifying the connectivity properties includes: For each vertex, a discrete residual equation is established that includes the conservation of mass and momentum. The cumulative term of the discrete residual equation is calculated based on the vertex volume, time step, and dynamic effective porosity. The flux term of the discrete residual equation is calculated by accumulating the phase mobility of adjacent vertices, the difference in fluid potential energy, and the fluid conductivity after the dynamic topology trimming operator has been set to zero. The partial derivatives of the discrete residual equations with respect to the multi-physics state vectors are calculated using automatic differentiation techniques, and then assembled to form the Jacobian matrix, which exhibits high-order asymmetric sparse block characteristics.

9. The method for evaluating the reservoir gas drive effect based on multiphysics coupling and topology trimming as described in claim 8, characterized in that, The iterative solution obtained by using a pre-conditional combined solver to update the multiphysics state vector of each node includes: The Jacobian matrix is ​​mathematically decoupled into a pressure principal matrix exhibiting elliptic partial differential characteristics and a saturation and component principal matrix exhibiting hyperbolic partial differential characteristics using the constrained pressure residual technique. The pressure principal submatrix is ​​preprocessed by applying an algebraic multigrid preconditioner. The preprocessed matrix system is input into a generalized minimum residual solver containing an incomplete LU decomposition preconditioner to perform iterative calculations, solve for the state vector increment, update the multiphysics state vector, until the infinite norm of the residual vector satisfies the preset convergence tolerance threshold.

10. The method for evaluating the reservoir gas drive effect based on multiphysics coupling and topology trimming as described in claim 7, characterized in that, The solution results from all time steps are extracted to generate a multidimensional dynamic gas drive performance evaluation index that includes tensors of dynamic effective sweep efficiency and solid-phase precipitation blockage damage index, including: When calculating the dynamic effective sweep efficiency, only the sweep volume of vertices whose dynamic effective porosity is greater than the percolation critical threshold is accumulated, excluding the volume of vertices whose connectivity is cut off by the dynamic topology pruning operator. Based on the edge data truncated by the dynamic topology pruning operator, the number of edges whose conductivity is set to zero in the orthogonal three-dimensional spatial directions is counted respectively, and then divided by the global initial total number of edges in the corresponding direction to construct the anisotropic solid-phase precipitation blockage damage index tensor.

Images