Gas drive effect evaluation method for gas drive oil exploitation multi-physical field coupling relationship

Abstract

This invention discloses a method for evaluating the multiphysics coupling relationship in gas-driven oil extraction. The method includes: acquiring reservoir parameters and constructing a discrete topological space based on a graph data structure in memory; inputting the state vector into a physical information neural network surrogate model to perform phase transition dimensionality reduction calculations, outputting solid phase precipitation and viscosity fields; calculating dynamic real porosity by combining mechanical strain and chemical dissolution; generating a dynamic topology trimming operator based on percolation theory, modifying conductivity properties in the underlying memory to sever the mathematical connectivity of physical blockage zones; assembling a fully implicit Jacobian matrix based on the updated conductivity and iteratively solving it; finally, eliminating isolated dead zone volumes and generating an evaluation index containing a tensor of real sweep efficiency and phase blockage damage index. This invention eliminates the matrix singularity divergence problem, significantly improving simulation computing power and gas-driven evaluation fidelity.

Description

Technical Field

[0001] This invention relates to the fields of computer-aided engineering and numerical simulation technology, and in particular to a method for evaluating the gas drive effect in gas-driven oil extraction based on the multiphysics coupling relationship. Background Technology

[0002] Existing computer numerical simulation and evaluation methods suffer from significant underlying algorithmic defects when dealing with extremely complex multi-field coupling (e.g., thermo-hydraulic-mechanical-chemical-phase). First, traditional phase equilibrium calculations heavily rely on iterative nested equations. When performing fluid-solid phase transition calculations on large-scale 3D meshes, the computational overhead 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 pore 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 divergence of 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 based on the multiphysics coupling relationship in gas-driven oil extraction, thereby addressing 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 gas drive effect in gas-driven oil extraction based on the multi-physics coupling relationship, 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 fluid-solid phase change and chemical dissipation, 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 true porosity of each node is calculated. The dynamic true 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 and generate a multidimensional dynamic gas drive performance evaluation index that includes dynamic real sweep efficiency and phase blockage damage index tensors.

[0006] As a preferred embodiment of the gas drive effect evaluation method for multi-physics coupling relationships in gas-driven oil extraction as described in this invention, wherein: the construction of a multi-physics discrete topological space based on a 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 fluid 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 gas drive effect evaluation method for multi-physics coupling relationships in gas-driven oil extraction as described in this invention, the multi-physics 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 fluid 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 gas drive effect evaluation method for multi-physics coupling relationships in gas-driven oil extraction as described in this invention, wherein: 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, and 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 gas drive effect evaluation method for gas-driven oil extraction based on multi-physics coupling relationships described in this invention, the calculation of the dynamic true porosity of each node based on the solid phase precipitation volume fraction field, rock mechanical volumetric strain, and chemical dissolution increment includes: Obtain the initial porosity of the nodes; Based on the initial porosity, rock compressibility coefficient, and current fluid 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 true porosity at the current time step is obtained by adding the porosity due to the effective stress compression effect in geomechanics to the chemical dissolution increment of carbonate rocks calculated based on hydrogen ion concentration, and subtracting the solid precipitation volume fraction field output by the physical information neural network surrogate model.

[0010] As a preferred embodiment of the gas drive effect evaluation method for multi-physics coupling relationships in gas-driven oil extraction as described in this invention, the method for determining the dynamic true porosity of each node based on the critical threshold of percolation theory and generating a dynamic topology trimming operator includes: Preset the percolation critical threshold that characterizes the complete rupture of the micropore-throat network. A permeability tensor over-permeability order reduction and dimensionality increase mapping rule is constructed. The difference between the dynamic real 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 true 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 gas drive effect evaluation method for multiphysics coupling relationships in gas-driven oil extraction according to the present invention, wherein: the conductivity matrix connectivity attribute of the discrete topology space is modified in memory by the dynamic topology pruning operator, including: When traversing any edge representing the interface in the undirected graph, determine the dynamic true porosity of the two adjacent vertices connected by that edge; If the dynamic true 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 of 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 gas drive effect evaluation method for multi-physics coupling relationships in gas-driven oil extraction as described in this invention, the step of assembling the fully implicit residual vector and Jacobian matrix of the nonlinear multi-physics control equations based on the conductivity matrix after modifying the connection 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 true 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 gas drive effect evaluation method for multiphysics coupling relationships in gas-driven oil extraction as described in this invention, wherein: 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 gas drive effect evaluation method for multi-physics coupling relationships in gas-driven oil extraction as described in this invention, the step of extracting the solution results data for all time steps and generating a multi-dimensional dynamic gas drive effect evaluation index containing dynamic true sweep efficiency and phase blockage damage index tensors includes: When calculating the dynamic true sweep efficiency, only the sweep volume of vertices whose dynamic true 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 phase 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 absolute stability 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 fluid-solid phase transformation nonlinear root-finding process into explicit tensor forward propagation on a 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 multi-dimensional evaluation mechanism based on graph theory-based edge truncation. When calculating sweep efficiency, this mechanism can accurately 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, accurately characterizing the degree of permeability damage in different reservoir directions, providing the most realistic 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 gas drive effect evaluation method for gas-driven oil extraction based on multi-physics field coupling relationships, 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. Example 1

[0020] Reference Figure 1 This is the first embodiment of the present invention, which provides a method for evaluating the gas drive effect of gas-driven oil extraction based on the multi-physics coupling relationship, 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.

[0021] 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.

[0022] 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 .

[0023] 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.

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

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

[0026] Furthermore, the dynamic multiphysics state vector is defined as: in, Characterizes the fluid 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.

[0027] Furthermore, the static physical property parameters of rocks are defined as follows: 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.

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

[0029] 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: in, The flow interface consists of nodes Pointing to node The unit normal vector. Operator Indicates 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.

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

[0031] 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 fully implicit Jacobian matrices in nonlinear iterations.

[0032] S2: Input the multiphysics state vector of each node in the discrete topological space into the pre-configured physical information neural network surrogate model, perform dimensionality reduction calculation of fluid-solid phase change and chemical dissipation, and output the solid precipitation volume fraction field and fluid viscosity field of each node.

[0033] 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.

[0034] 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. : 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.

[0035] Furthermore, the network's output layers are configured to output prediction tensors through forward propagation. : 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. The viscosity of the oil phase 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 pore throats.

[0036] 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 traditional real gas equation of state (PR-EOS), within a preset temperature range (e.g., 300K~450K) and pressure range (1MPa~80MPa) of the reservoir, combined with Latin hypercube sampling (LHS), a dataset covering multiple components (including...) was generated. methane to heavy hydrocarbons A million thermodynamic equilibrium sample pairs (e.g., [list of samples]) were used as the benchmark dataset. Next, normalization of the input data was performed. Due to the significant difference in magnitude between pore pressure and component mole fractions, a logarithmic range normalization strategy was adopted for the input tensor. All elements are strictly mapped to the interval [-1, 1] to avoid gradient explosion during neural network training. Next, a deep network topology is constructed. In this architecture, the Multilayer Perceptron (MLP) contains one input layer, 5-8 hidden layers (each with 128-256 neurons), and one output layer. Furthermore, to ensure the continuity and smoothness of the partial differential equation differentiation, the hidden layer uses the Swish or Tanh activation function with second-order continuous differentiability, and the output layer uses the Softmax activation function to ensure forced normalization of the summation of component mole fractions. Finally, a hybrid optimization strategy is adopted during offline training. That is, the Adam optimizer is used for global fast optimization in the early stage, with an initial learning rate set to... Once the loss function's decline plateaus, switch to a second-order L-BFGS optimizer for local exact optimization until the total loss function converges to a preset value. For orders below a certain size, complete the fixed configuration of model parameters.

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

[0038] 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: in, These are dynamic adaptive weighting coefficients, with values ​​ranging from [0,1]. .

[0039] 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.

[0040] 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: 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.

[0041] 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 absolutely conforms to the objective physical laws.

[0042] Specifically, for the mole normalization conservation operator in the above formula... This term is mainly used to constrain the conservation of matter, that is: 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.

[0043] S203: Perform online inference output and chemical dissipation feature extraction.

[0044] 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. The data is then processed in batches using the PINN model with embedded physical mechanisms to perform rapid forward propagation, directly extracting the solid precipitation volume fraction field of each node. With fluid viscosity field .

[0045] It should be noted that this step aims to skip the nonlinear root-finding process of the traditional Rachford-Rice equation, thereby achieving a generational leap in computing power at the level of dimensionality reduction of multiphysics fields.

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

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

[0048] Specifically, the prediction tensor output by the PINN proxy model Substituting into the thermodynamic isofugacity residual equation (i.e., the above) The difference defined in the code is used to calculate the tolerance response value of each node in real time. When any node is detected ( In this embodiment, the preset minimum physical tolerance can be set to... When the system determines that the neural network prediction for that node is at risk of distortion, the computer will trigger a rollback operator, discard the prediction result of that node, and seamlessly roll back the local phase calculation task of that node to the traditional Newton-Raphson iterative solution process for state equations.

[0049] 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 absolute numerical convergence of the global solution of partial differential equations while significantly accelerating the overall simulation.

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

[0051] It should be noted that in traditional commercial reservoir numerical simulators, porosity is often simplified to a single-variable linear function of fluid 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 three factors: 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 real 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.

[0052] Specifically, first, obtain the nodes. initial porosity Extract the current time step Fluid 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. : in, is the equivalent compression coefficient of the rock skeleton; The initial fluid 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.

[0053] 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 the fluid pore pressure itself (where an increase or decrease in fluid pressure directly causes the pores to expand or contract); 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.

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

[0055] 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 : 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.

[0056] It should be noted that by calculating the chemical dissolution increment of carbonate rocks, the expansion effect of the propagation of the acidic fluid front in porous media on the micropore throat can be quantified from a molecular dynamics perspective. 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.

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

[0058] 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 true porosity : Furthermore, by fully expanding it, we can obtain the multiphysics coupled porosity master equation in the present invention: It should be noted that the master equation perfectly reflects the reservoir's state during... The complex physical competition mechanism in the gas drive process. Among them, mechanical compaction ( ) and solid precipitation ( ) together act as a negative feedback mechanism, attempting to shrink pores or even block microscopic pore throats; while pressure support and chemical dissolution ( This, acting as a positive feedback mechanism, attempts to expand the porosity. When the interplay of these multiple physical effects leads to local... When the descent approaches the minimum value, the topology pruning operation of the underlying operator will be triggered directly.

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

[0060] 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 an empirical continuous function (such as the classic Kozeny-Carman formula). However, during gas injection, due to intense solid phase precipitation (asphaltite blockage) and in-situ compaction, the porosity at local nodes can approach its minimum value. In such cases, traditional formulas will calculate an extremely small but non-zero permeability (e.g., ... (millidarcy). This not only physically violates the objective fact that the orifice throat is completely broken and disconnected, but also causes disastrous consequences at the underlying computer algorithm level. Specifically, an extremely small conductivity will lead to a sharp increase in the condition number of a large sparse Jacobian matrix during the subsequent discretization and assembly of partial differential equations, causing severe 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 topological pruning operator capable of achieving physical disconnection and mathematical decoupling. The implementation process is as follows: S401: Constructing a dynamic permeability tensor evolution model based on percolation theory.

[0061] Furthermore, a pre-defined percolation critical threshold is set to characterize the complete rupture of the micropore-throat network. From a physical perspective, due to the presence of bound water and dead pores in porous media, when the actual porosity decreases to... When the value is not absolute zero, but usually between 0.02 and 0.05, the fluid can no longer find any continuous macroscopic flow channel in three-dimensional space, that is, percolation phase transition has occurred.

[0062] Furthermore, obtain the node at the current time step calculated in S3. Dynamic true porosity Construct a permeability tensor over-permeability reduction and dimensionality-upgrading mapping rule, and compute the nodes. dynamic permeability tensor : 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.

[0063] 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 absolute 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.

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

[0065] 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: 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.

[0066] It should be noted that as long as adjacent vertices or The dynamic true porosity of any vertex is less than or equal to the critical threshold. ,but The value will undergo a binary mutation, collapsing directly from an absolute 1 to an absolute 0.

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

[0068] 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. .

[0069] 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: It should be noted that, through the above steps, if it is determined that a percolation disconnection has occurred, then... If the value is 0, the computer will force the fluid conductivity of the edge to be zero. The corresponding values ​​in the conductivity matrix are directly assigned absolute zero (rather than extremely small values). From a mathematical topological perspective, this is equivalent to directly severing vertices in the in-memory graph data structure. With vertex The edges between them. Through this simple and efficient algebraic operation, the mathematical connectivity between adjacent vertices in the discrete solution space of partial differential equations can be completely severed. This directly avoids the generation of minimal values ​​when assembling flux terms, thereby completely eliminating the singularity and oscillatory non-convergence problems in the nonlinear Newton iterative equation system, enabling a leap in fault tolerance and computational efficiency of multi-field coupled simulations.

[0070] 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.

[0071] It should be noted that carbon dioxide ( The multi-physics coupling of "thermal-hydraulic-mechanical-chemical-phase" induced by oil displacement is a highly nonlinear process. Traditional explicit or implicit pressure-explicit saturation (IMPES) solution schemes are easily constrained by the Courant-Friedrich-Lyuvi (CFL) condition, leading to computational divergence. Therefore, this invention employs a fully implicit method (FIM) for time and space discretization. Furthermore, it is important to emphasize that traditional fully implicit methods collapse due to the singularity of the Jacobian matrix when the pore throat is completely blocked. However, the solution in this invention, by using a dynamic topology pruning operator in S4 to sever the mathematical connections that cause singularities at the underlying data structure level, perfectly ensures the smooth 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.

[0072] 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. : 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.

[0073] Furthermore, the cumulative terms of this discrete residual equation Based on vertex volume, time step, and the dynamic true porosity obtained in S3. The calculation yielded: 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.

[0074] 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: 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 fluid pore pressure and gravitational potential).

[0075] 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.

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

[0077] 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.

[0078] 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. : 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.

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

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

[0081] 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.

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

[0083] 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: 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.

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

[0085] 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 gas saturation and pore volume across the entire grid. However, in real-world injection... During oil displacement, if severe solid-phase precipitation (such as asphaltene blockage) or extreme mechanical closure occurs in local pore throats, the fluid in these areas will completely lose its physical flow capacity, becoming isolated dead zones. Traditional statistical algorithms, unable to identify network disconnections, still include the residual fluid in these dead zones in the swept volume, resulting in a significantly inflated final gas displacement efficiency assessment. Therefore, this invention introduces a fully fidelity multidimensional dynamic evaluation method based on the aforementioned dynamically clipped records of the graph data structure. 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.

[0086] 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 ) (such as dynamic real 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).

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

[0088] Furthermore, a Boolean conditional decision mechanism based on graph theory vertices is introduced when calculating the dynamic true sweep efficiency. That is, only the sweep volume of vertices whose dynamic true 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: in, For the first The dynamic real sweep efficiency at each time step, a dimensionless quantity; The total number of vertices in an undirected graph. For nodes The volume of rock is controlled by its volume; and Representing nodes respectively In the Dynamic real 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 by 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.

[0089] Furthermore, when the node of ( When the percolation threshold is 1, and it has at least one edge whose conductivity is not set to zero, 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 .

[0090] It should be noted that the introduction The issue lies in the fact that when the injected gas front triggers severe asphaltene deposition, causing some mesh nodes that originally contained gas to be completely surrounded and severed by the surrounding clogging zone, the operator... The volume contribution of that node in the numerator is reduced to zero. Algorithmically, this eliminates visible but unminable spurious swept volumes, providing the most conservative and reliable numerical basis for predicting actual mine production.

[0091] S603: Construct an anisotropic phase blockage damage exponential tensor based on graph theory-based edge truncation and dimensionality reduction mapping.

[0092] Furthermore, because reservoir damage often does not occur uniformly due to differences in sedimentary facies, the initial permeability of the formation varies greatly in the horizontal and vertical directions. Therefore, due to supercritical... The blockage behavior induced by phase transitions exhibits strong anisotropy in three-dimensional space. Therefore, this step, based on the edge data truncated by the dynamic topology pruning operator at the bottom layer, counts the number of edges whose conductivity is set to zero in orthogonal three-dimensional directions. These counts are then divided by the total initial global edge count in the corresponding direction to construct the anisotropic phase blockage damage exponent tensor. : Furthermore, for any directional component on the main diagonal of the matrix (in... Direction For example, its calculation formula is: 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.

[0093] 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).

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

[0095] 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.

[0096] 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.

[0097] 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 gas drive effect in gas-driven oil extraction based on multi-physics coupling relationships, 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 fluid-solid phase change and chemical dissipation, 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 true porosity of each node is calculated. The dynamic true 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 and generate a multidimensional dynamic gas drive performance evaluation index that includes dynamic real sweep efficiency and phase blockage damage index tensors.

2. The gas drive effect evaluation method for multi-physics coupling relationship in gas-driven oil extraction 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 fluid 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 gas drive effect evaluation method for multi-physics coupling relationship in gas-driven oil extraction 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 fluid 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 gas drive effect evaluation method for multi-physics coupling relationship in gas-driven oil extraction 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. 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.

5. The gas drive effect evaluation method for multi-physics coupling relationship in gas-driven oil extraction as described in claim 3, characterized in that, The dynamic true 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 fluid 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 true porosity at the current time step is obtained by adding the porosity due to the effective stress compression effect in geomechanics to the chemical dissolution increment of carbonate rocks calculated based on hydrogen ion concentration, and subtracting the solid precipitation volume fraction field output by the physical information neural network surrogate model.

6. The gas drive effect evaluation method for multi-physics coupling relationship in gas-driven oil extraction as described in claim 5, characterized in that, The determination of the dynamic true porosity of each node based on the critical threshold of percolation theory, and the generation of a dynamic topology trimming operator, includes: Preset the percolation critical threshold that characterizes the complete rupture of the micropore-throat network. A permeability tensor over-permeability order reduction and dimensionality increase mapping rule is constructed. The difference between the dynamic real 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 true 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 gas drive effect evaluation method for multi-physics coupling relationship in gas-driven oil extraction 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 true porosity of the two adjacent vertices connected by that edge; If the dynamic true 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 of 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 gas drive effect evaluation method for multi-physics coupling relationship in gas-driven oil extraction 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 true 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.

9. The gas drive effect evaluation method for multi-physics coupling relationship in gas-driven oil extraction 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 gas drive effect evaluation method for multi-physics coupling relationship in gas-driven oil extraction 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 dynamic true sweep efficiency and phase blockage damage index tensors, including: When calculating the dynamic true sweep efficiency, only the sweep volume of vertices whose dynamic true 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 phase blockage damage index tensor.

Images