A simulation method for gas-water two-phase flow based on pore evolution

Abstract

This invention discloses a gas-water two-phase flow simulation method based on pore size evolution, belonging to the field of multiphase flow simulation technology. The method includes: discretizing the porous medium into a pore network model composed of pore nodes and connecting throats, and defining geometric parameters for each throat; updating the equivalent radius, cross-sectional area, perimeter, and shape factor of the throat at the next time step based on the pore size evolution equation, and further updating the single-phase fluid conductivity and capillary pressure of the throat at the next time step; solving for the fluid pressure field and flow rate field satisfying the law of mass conservation and Kirchhoff's equations in the pore network model according to the single-phase fluid conductivity and capillary pressure of the throat at the next time step, and advancing the gas-water two-phase interface until the termination condition is met, outputting the simulation results of the gas-water two-phase flow model. This method breaks through the static assumption of fixed pore throat geometry, achieving accurate characterization of the dynamic changes of pore throat structure in real underground environments, thus improving the accuracy of the simulation results.

Description

Technical Field

[0001] This invention relates to the field of multiphase flow simulation technology, specifically to a gas-water two-phase flow simulation method based on pore size evolution. Background Technology

[0002] Pore ​​network models, as an important tool for characterizing two-phase flow at the pore-throat scale in porous media, abstract complex porous media into a topological network composed of pores (nodes) and throats (links). They can efficiently reproduce typical phenomena such as fluid intrusion fingering, multiphase fluid residual distribution and trapping, and provide an interpretable pore-scale physical basis for intrusion threshold distribution, relative permeability curves, intrusive fluid breakthrough behavior and residual fluid saturation. They are widely used in oil and gas field development, geological storage, soil remediation and fuel cells.

[0003] However, existing pore network models are generally based on the assumption of "geostatic stability," meaning that the shape and size of pores and throats remain unchanged throughout the simulation. The model relies solely on fixed throat lengths, equivalent radii, and shape parameters to advance the fluid phase interface and solve the pressure field. While this assumption simplifies the calculation process, it ignores the objective fact that the pore-throat structure dynamically evolves with changing conditions in real underground environments, leading to systematic deviations between model predictions and physical reality. Summary of the Invention

[0004] The technical problem to be solved by the present invention is that existing models based on geometric static assumptions cannot reflect the changes in pore throat structure with working conditions in the real underground environment. The purpose is to provide a gas-water two-phase seepage simulation method based on pore size evolution, which solves the above-mentioned problem.

[0005] This invention is achieved through the following technical solution:

[0006] In a first aspect, the present invention provides a gas-water two-phase flow simulation method based on pore size evolution, comprising:

[0007] The porous medium is discretized into a pore network model consisting of pore nodes and connecting throats, and geometric parameters are defined for each throat; the geometric parameters include length, equivalent radius, cross-sectional area, perimeter, corner half angle and shape factor;

[0008] Based on the aperture evolution equation, the equivalent radius of the throat is updated at the next time step; the aperture evolution equation is used to describe the evolution of the equivalent radius of the throat driven by mechanical, hydration and chemical effects.

[0009] Update the cross-sectional area, perimeter, and shape factor of the throat at the next time step based on the equivalent radius of the throat at the next time step.

[0010] The single-phase fluid conduction capacity of the throat at the next time step is updated based on the length of the throat, the cross-sectional area of ​​the throat at the next time step, and the shape factor.

[0011] Update the capillary pressure of the larynx at the next time step based on the corner half angle of the larynx, the equivalent radius of the larynx at the next time step, and the shape factor.

[0012] Based on the single-phase fluid conductivity and capillary pressure of the throat at the next moment, the fluid pressure field and flow field that satisfy the law of conservation of mass and Kirchhoff's equations are solved in the pore network model, and the gas-water two-phase interface is advanced until the termination condition is met, and the simulation results of gas-water two-phase flow are output.

[0013] Optionally, the aperture evolution equation is a multiplicative equation, expressed as follows:

[0014]

[0015] in, The throat between adjacent pore nodes i and j The equivalent radius at any given moment; Let be the equivalent radius of the throat at time t; Indicates the time step; The effective stress of the throat at the current moment; The initial stress of the throat; The local water saturation of the throat at the current moment; The initial water saturation of the throat; The equivalent geometric rate of dissolution of the throat per unit time; The equivalent geometric rate of sedimentation in the throat per unit time; , , and The aperture evolution coefficient is determined experimentally.

[0016] Optionally, the aperture evolution equation is an additive equation, expressed as follows:

[0017]

[0018] in, This represents the rate of change of the equivalent radius of the throat between adjacent pore nodes i and j; Let be the equivalent radius of the throat at time t; The effective stress of the throat at the current moment; The initial stress of the throat; The local water saturation of the throat at the current moment; The initial water saturation of the throat; The equivalent geometric rate of dissolution of the throat per unit time; The equivalent geometric rate of sedimentation in the throat per unit time; , , and The aperture evolution coefficient is determined experimentally.

[0019] Optionally, the Calibration was performed through triaxial densification experiments, hydrostatic densification experiments, or pore pressure cyclic loading experiments; Calibration was performed using permeability or pore size imaging data under changes in moisture absorption / dehydration / water saturation; The effect of the thickness of the dissolved layer on the throat radius was experimentally calibrated; Experimental calibration of the effect of sediment thickness on throat radius variation.

[0020] Optionally, the formula for calculating the single-phase fluid conductivity is as follows:

[0021] in, The α-phase fluid conductivity of the throat between adjacent pore nodes i and j at time t; α is a phase characterization parameter of gas or water. , w represents gas, and w represents water; This represents the cross-sectional area of ​​the throat at time t; The shape factor of the larynx at time t; Represents a geometric function determined by the shape factor; Let be the phase viscosity of the α-phase fluid; The length of the larynx.

[0022] Optionally, the formula for calculating the capillary pressure is as follows:

[0023]

[0024] in, The capillary pressure of the throat between adjacent pore nodes i and j at time t when gas and water coexist. The interfacial tension between air and water; Let be the equivalent radius of the throat at time t; The corner half-angle of the throat; It is the wetting contact angle between the gas-liquid two-phase interface and the wall of the throat. The shape factor of the larynx; A correction function to account for non-circular cross sections and corner film effects.

[0025] Optionally, the step of solving the fluid pressure field and flow rate field satisfying the law of conservation of mass and Kirchhoff's equations in the pore network model based on the single-phase fluid conductivity and capillary pressure of the throat at the next moment, and advancing the gas-water two-phase interface, includes:

[0026] Based on the single-phase fluid conductivity and capillary pressure of the throat at the next moment, the law of conservation of mass and Kirchhoff's equations are solved to obtain the fluid flow rate of the throat.

[0027] Based on the fluid flow rate and time step of the throat, calculate the fluid volume fraction of the throat and the axial advance distance of the gas-water interface in the throat.

[0028] The position of the gas-water two-phase interface in the pore network model is updated based on the axial advance distance.

[0029] Optionally, the law of conservation of mass and Kirchhoff's equations are as follows:

[0030]

[0031] in, The α-phase fluid flow rate is the throat between adjacent pore nodes i and j. The α-phase fluid saturation of the throat; The α-phase fluid conductivity of the throat between pore nodes i and j; The pressure of the α-phase fluid within pore node i; The pressure of the α-phase fluid within pore node j; This indicates the directionality of the capillary pressure of the α-phase fluid in the throat. The capillary pressure of the larynx.

[0032] Optionally, the α-phase fluid pressure within the pore nodes is obtained in the following manner:

[0033] Based on the fluid saturation of the throat at time t, the interconnected clusters of the gas phase and the water phase are identified respectively;

[0034] For each phase of the connected cluster, nonlinear equations concerning the pore node pressure are assembled separately;

[0035] The nonlinear equations were solved using a sparse linear algebra algorithm to obtain the α-phase fluid pressures within all pore nodes; the nonlinear equations are as follows:

[0036]

[0037]

[0038]

[0039] in, This is a vector consisting of the α-phase fluid pressures within all the pore nodes to be solved; It is a coefficient matrix; The right-hand term; Let α be the α-phase fluid elastic quantity of pore node i; For time step; Let be the pressure of the α-phase fluid inside pore node i at time t+1; Represents the cluster of neighboring nodes of pore node i; k is... Any pore node in the middle; The α-phase fluid saturation of the throat at time t; The α-phase fluid conductivity of the throat at time t; Let be the pressure of the α-phase fluid within pore node j at time t; The capillary pressure of the larynx at time t; The directionality of the capillary pressure of the α-phase fluid in the throat at time t.

[0040] Optionally, the time step satisfies: Where Comax is a preset time step limit index, with a value range of [0.1, 0.5]. The α-phase fluid flow rate of the throat; For time step; Let be the cross-sectional area of ​​the throat at time t; max() represents the maximum value function.

[0041] Compared with the prior art, the present invention has the following advantages and beneficial effects:

[0042] This invention provides a gas-water two-phase flow simulation method based on pore size evolution. By introducing a pore size evolution equation driven by mechanical, hydration, and chemical processes, it breaks through the static assumption of "fixed pore throat geometry" in traditional pore network models, achieving accurate characterization of the dynamic changes in pore throat structure in real underground environments. Through a linkage update mechanism, the evolution of the throat's equivalent radius is coupled in real time and bidirectionally with key geometric parameters (cross-sectional area, perimeter, shape factor) and core flow properties (single-phase conductivity, gas-water two-phase capillary pressure), thereby constructing a strongly coupled feedback closed loop of "pore size evolution - capillary force - fluid flow" within the pore network model framework. This closed-loop system fully reveals the dynamic interaction mechanism between porous media and fluid flow, significantly improving the model's prediction accuracy from a physical perspective. It narrows the gap between simulation results and actual engineering observations in key dynamics such as intrusion fingering, residual saturation, and breakthrough time, providing a reliable tool for accurate prediction and optimization decision-making in projects such as sequestration, energy storage, and oil and gas development. Attached Figure Description

[0043] To more clearly illustrate the technical solutions of the exemplary embodiments of the present invention, the accompanying drawings used in the embodiments will be briefly described below. It should be understood that the following drawings only show some embodiments of the present invention and should not be considered as a limitation of the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort. In the drawings:

[0044] Figure 1 A schematic flowchart of a gas-water two-phase flow simulation method based on pore size evolution is provided for an embodiment of this application.

[0045] Figure 2 A schematic diagram of a circular cross-section pore network model provided in an embodiment of this application;

[0046] Figure 3 A schematic diagram of the simulation results of aperture evolution provided in the embodiments of this application;

[0047] Figure 4 A schematic diagram of the simulation results of the static aperture provided in the embodiments of this application;

[0048] Figure 5 A schematic diagram of a triangular / square throat pore network model provided in an embodiment of this application;

[0049] Figure 6 A schematic diagram showing simulation results considering the corner membrane, provided for embodiments of this application;

[0050] Figure 7 This is a schematic diagram of the passivation effect of fluid intrusion provided in an embodiment of this application. Detailed Implementation

[0051] To make the objectives, technical solutions, and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the embodiments and accompanying drawings. The illustrative embodiments and descriptions of the present invention are only used to explain the present invention and are not intended to limit the present invention.

[0052] Under practical engineering conditions, the geometry of pores and throats changes over time driven by multiple mechanisms. First, changes in effective stress can cause pore compression or expansion, altering the equivalent radius and conductivity of local channels, thus affecting the intrusion process. Second, changes in water content can drive hydration expansion or contraction of clay minerals, showing greater sensitivity to small throats, often leading to confinement of non-wetting phase channels and enhanced capillary trapping. Third, geochemical processes such as dissolution and precipitation can reconstruct the throat cross-section and cause changes in wettability, manifested as precipitation preferentially "passivating" small pore throats, and dissolution promoting channel expansion, exhibiting significant path dependence and selectivity. For non-circular cross-section throats (such as triangular or square ones), the formation and stability of corner films are also jointly regulated by wettability and corner geometry, affecting water connectivity and local trapping under low-velocity conditions.

[0053] Experimental observations and engineering data show that when considering geometric evolution and changes in wettability, the system behavior exhibits systematic deviations compared to the static geometric model: the overall difficulty of intrusion increases, the fractal complexity of fingering structures decreases, residual aqueous phase saturation increases, and breakthrough time is delayed; in scenarios with sedimentation, the early pressure drop is significantly increased, while the gas-connection pathway tends to stabilize in the later stages. These phenomena reflect a significant bidirectional feedback between "pore geometry-capillary pressure-connectivity-fluid flow," which cannot be accurately captured by static geometric assumptions. The main shortcomings are as follows:

[0054] First, static geometric assumptions fail to reflect the pore-throat evolution and wettability drift driven by stress, water content, and geochemical processes, leading to deviations in predictions of intrusion difficulty, relative permeability endpoints, and connectivity. Second, existing methods often employ weak coupling or split solutions, failing to fully describe the strong coupling feedback between geometric changes and capillary rules and flow behavior. Third, the coexistence of multiple timescales makes explicit or coarse-step updates prone to numerical instability and error accumulation, affecting the reliability of long-term simulations. Fourth, parameters such as pore size evolution, water content-induced expansion / contraction, and the geometric influence of precipitation largely rely on experience or table lookups, lacking a unified dimensional closure and traceable calibration process.

[0055] The shortcomings mentioned above directly affect the reliability of engineering predictions and decisions. Ignoring geometric evolution often underestimates the difficulty of intrusion and residual saturation, misjudges sweep efficiency and breakthrough time, and thus affects the optimization of injection and production strategies and risk assessment. Sedimentation and hydration expansion can lead to selective "passivation" of critical throats, triggering abrupt connectivity changes and channel migration, inducing early pressure drop anomalies, capacity reduction, or accidental leakage risks. Under long-term injection and production conditions, the cumulative effect of geometric evolution can also alter the system lifetime and maintenance cost curves, having a decisive impact on the economics of the scheme.

[0056] Therefore, this application provides a method for simulating gas-water two-phase flow based on pore size evolution. Please refer to... Figure 1 This is a schematic flowchart of a gas-water two-phase flow simulation method based on pore size evolution provided in an embodiment of this application. The following is a description of... Figure 1 The gas-water two-phase flow simulation method based on pore size evolution is introduced.

[0057] S1. Discretize the porous medium into a pore network model consisting of pore nodes and connecting throats, and define geometric parameters for each throat.

[0058] In the specific implementation process, a pore network model is used to discretize the porous medium. The pore network model consists of pores and throats. Pores are abstracted as nodes, representing the main spaces for fluid storage and pressure balance; throats are abstracted as channels connecting the pores, serving as key components controlling fluid flow and capillary pressure. Each node is connected to multiple neighboring nodes through different throats. The set of all neighboring nodes directly connected to any given node i is defined as the neighboring node cluster N(i) of that node.

[0059] In the porosity network model, for any throat between adjacent porosity nodes i and j It is necessary to define its complete geometric parameters to accurately characterize its morphology and transmission properties. These geometric parameters include, but are not limited to:

[0060] 1. Length: Characterizes the axial extension dimension of the throat.

[0061] 2. Equivalent radius: The characteristic radius used to characterize the flow capacity of the throat, usually defined as the radius of a circle with the same area as the cross-section of the throat.

[0062] 3. Cross-sectional area: The cross-sectional area of ​​the throat perpendicular to the flow direction.

[0063] 4. Perimeter: The total length of the boundary of the throat's cross-section.

[0064] 5. Corner half angle: Used to quantify the sharpness of the corner of the throat section.

[0065] 6. Shape factor: Used to comprehensively characterize the degree of deviation of the cross-sectional shape from a circle.

[0066]

[0067] in, The throat between adjacent pore nodes i and j The shape factor at time t; For the throat The cross-sectional area at time t; For the throat The circumference at time t.

[0068] S2. Based on the aperture evolution equation, update the equivalent radius of the throat at the next time step.

[0069] The pore size evolution equation describes the evolution of the throat's equivalent radius driven by mechanical, hydration, and chemical effects. In the pore network model, considering mechanical (effective stress), water film (water content), and chemical (dissolution / precipitation) effects, a multiplicative term suitable for gradual changes and guaranteeing positive properties, and an additive term suitable for small deformation linear approximations, are proposed to realize the evolution of the throat's equivalent radius at different time scales.

[0070] In one possible embodiment, the aperture evolution equation is a multiplicative equation, expressed as follows:

[0071]

[0072] in, The throat between adjacent pore nodes i and j exist The equivalent radius at any given moment; For the throat The equivalent radius at time t; Indicates the time step; For the throat Effective stress at the current moment; For the throat The initial stress; For the throat Local water saturation at the current moment; For the throat The initial water saturation; For the throat The equivalent geometric rate of dissolution per unit time; For the throat The equivalent geometric rate of sedimentation per unit time; , , and The aperture evolution coefficient is determined experimentally.

[0073] In one possible embodiment, the aperture evolution equation is an additive equation, expressed as follows:

[0074]

[0075] in, This represents the throat between adjacent pore nodes i and j. The rate of change of the equivalent radius; For the throat The equivalent radius at time t; The effective stress of the throat Ω at the current moment; The initial stress of the throat Ω; The local water saturation of the throat Ω at the current moment; For the throat The initial water saturation; For the throat The equivalent geometric rate of dissolution per unit time; For the throat The equivalent geometric rate of sedimentation per unit time; , , and The aperture evolution coefficient is determined experimentally.

[0076] Aperture evolution coefficients include mechanical evolution coefficients. Water film evolution coefficient Dissolution evolution coefficient and precipitation evolution coefficient These coefficients together ensure that the equivalent radius of the throat is accurate during numerical simulation. The evolution of these coefficients follows physical laws and remains consistently positive. The following sections detail how each evolution coefficient is calibrated.

[0077] (1) Mechanical evolution coefficient

[0078] It is used to characterize the sensitivity of effective stress changes to throat compression or expansion, and can be calibrated through triaxial compaction experiments, hydrostatic compaction experiments, or pore pressure cyclic loading experiments.

[0079] Specifically, triaxial densification experiments, hydraulic densification experiments, or pore pressure cyclic loading experiments can be used. In these experiments, the porosity or permeability changes of core samples under different effective stresses are measured, and the mechanical evolution coefficients can be determined through inversion analysis. .

[0080] (2) Water film evolution coefficient

[0081] The sensitivity of the throat wall to expansion or contraction caused by changes in water saturation can be characterized by permeability or pore size imaging data under hygroscopic-dehumidification / water saturation changes.

[0082] Specifically, by measuring the permeability changes of core samples in real time during the hygroscopic-dehumidification cycle, or by using high-resolution microscopic imaging technology to obtain dynamic data on pore structure, and by fitting the observed data with model predictions, the water film evolution coefficient can be calibrated. .

[0083] (3) Dissolution evolution coefficient

[0084] The efficiency of the throat radius expansion caused by the dissolution reaction can be characterized by experimental calibration of the effect of the dissolution layer thickness on the throat radius change.

[0085] Specifically, by conducting core displacement experiments and injecting fluids that can induce mineral dissolution, the thickness and spatial distribution of the dissolution layer within the throat can be measured using characterization techniques such as scanning electron microscopy. Combined with the reaction time, the dissolution evolution coefficient can be determined by inversion. .

[0086] (4) Precipitation evolution coefficient

[0087] The efficiency of the throat radius reduction caused by the precipitation reaction can be characterized by experimental calibration of the effect of the deposition layer thickness on the throat radius change.

[0088] Specifically, by conducting core displacement experiments and injecting supersaturated fluid to induce mineral precipitation, the thickness of the sedimentary layer within the throat is measured using characterization techniques such as scanning electron microscopy. Combined with the reaction time, the precipitation evolution coefficient can be determined by inversion. .

[0089] In this embodiment of the application, the above-described experimental calibration process provides reliable and traceable parameter inputs for the pore network model, thereby ensuring its accuracy and reliability in simulating the gas-water two-phase seepage process under real geological conditions.

[0090] S3. Based on the equivalent radius of the throat at the next time step, update the cross-sectional area, perimeter, and shape factor of the throat at the next time step.

[0091] In the specific implementation process, based on the equivalent radius of the throat Ω at time t... Update the cross-sectional area of ​​the throat at time t, based on the preset throat cross-sectional geometry. and perimeter Based on cross-sectional area and perimeter Update the shape factor of the larynx Ω at time t. .

[0092] S4. Update the single-phase fluid conduction capacity of the throat at the next time step based on the throat length, the cross-sectional area of ​​the throat at the next time step, and the shape factor.

[0093] The formula for calculating the conductivity of a single-phase fluid is as follows:

[0094]

[0095] in, The throat between adjacent pore nodes i and j The conductivity of the α-phase fluid at time t; α is a characteristic parameter of a certain phase of gas or water; g represents gas; w represents water; A ij (t) represents the throat. The cross-sectional area at time t; For the throat The shape factor at time t; Represents a geometric function determined by the shape factor; Let be the phase viscosity of the α-phase fluid; For the throat The length.

[0096] When the throat Ω is a circular tube:

[0097]

[0098] S5. Update the capillary pressure of the larynx at the next time step based on the corner half angle of the larynx, the equivalent radius of the larynx at the next time step, and the shape factor.

[0099] The formula for calculating capillary pressure is as follows:

[0100]

[0101] in, This refers to the throat between adjacent pore nodes i and j when gas and water coexist. The capillary pressure at time t; σ is the gas-water interfacial tension; For the throat The equivalent radius at time t; For the throat The corner half angle; The gas-water two-phase interface and throat The wetting contact angle between the walls; For the throat Shape factor; To account for the correction function for non-circular cross sections and corner film effects, when approximating a circular cross section, .

[0102] S6. Based on the single-phase fluid conductivity and capillary pressure of the throat at the next moment, solve the fluid pressure field and flow field that satisfy the law of conservation of mass and Kirchhoff's equations in the pore network model, and advance the gas-water two-phase interface until the termination condition is met, and output the simulation results of the gas-water two-phase flow model.

[0103] In one possible embodiment, based on the single-phase fluid conductivity and capillary pressure of the throat at the next moment, the fluid pressure field and flow field satisfying the law of conservation of mass and Kirchhoff's equations are solved in the pore network model, and the gas-water two-phase interface is advanced, including:

[0104] Based on the single-phase fluid conductivity and capillary pressure of the throat at the next moment, the mass conservation law and Kirchhoff's equations are solved to obtain the fluid flow rate of the throat. Based on the fluid flow rate and time step of the throat, the fluid volume fraction of the throat and the axial advance distance of the gas-water interface in the throat are calculated. Based on the axial advance distance, the position of the gas-water interface in the pore network model is updated.

[0105] In the specific implementation process, firstly, the fluid flow rate in the throat is calculated. The fluid flow rate in the throat obeys the law of conservation of mass and Kirchhoff's equations, as follows:

[0106]

[0107] in, The throat between adjacent pore nodes i and j α-phase fluid flow rate; For the throat α-phase fluid saturation; For the throat The conductivity of the α-phase fluid at time t; For the throat Capillary pressure; Indicates the throat The directionality of capillary pressure in the α-phase fluid; when the capillary pressure is resistance. When capillary pressure is the driving force, =1, when the throat When the fluid inside is a single-phase fluid, there is no capillary pressure. . Let be the pressure of the α-phase fluid within pore node i.

[0108] Secondly, according to the throat Calculate the throat by considering the fluid flow rate and time step. The fluid volume fraction.

[0109]

[0110] in, For the throat α-phase fluid volume fraction; For the throat α-phase fluid flow rate; For time step;

[0111] Then, according to the throat Given the fluid volume fraction and the cross-sectional area of ​​the throat at time t, calculate the axial advance distance of the gas-water interface in the throat.

[0112]

[0113] in, Indicates the time step Inside, the gas-water interface is in the throat. The axial advance distance in the fluid; used to determine fluid breakthrough, invasion, and retreat.

[0114] Finally, based on the axial propulsion distance, the gas-water two-phase interface in the throat is updated. The position of the gas-water interface is determined in the time step, thus clarifying the starting position of the gas-water interface in the next time step and completing the interface advancement for the current time step. After the interface advancement, if the preset termination condition is met (such as the global mass conservation error being less than the preset error), the simulation results of the gas-water two-phase flow are output. If the termination condition is not met, the next time step Δt is adaptively adjusted, and S2-S6 are continued.

[0115] The flow stability and fingering intensity of gas-water two-phase fluids in a pore network model are jointly controlled by the number of local capillaries, wettability, and heterogeneity.

[0116] The formula for calculating the number of local capillaries is as follows:

[0117]

[0118] in, Indicates the number of capillaries in a given area; Let be the phase viscosity of the α-phase fluid; Let be the phase velocity of the α-phase fluid; For interface tension.

[0119] To characterize the effect of flow velocity on wettability, the number of local capillaries was used. against the larynx Static wetting contact angle Dynamic correction is performed to obtain the dynamic wetting contact angle. Its update formula is:

[0120]

[0121] in, For the updated dynamic wetting contact angle; This is the static wetting contact angle; It is a corner half angle; This indicates the number of capillaries in a given area.

[0122] In one possible implementation, the time step satisfies: ;in, The preset time step limit index has a value range of [0.1, 0.5]. For the throat α-phase fluid flow rate; For time step; For the throat The cross-sectional area at time t.

[0123] In one possible embodiment, the α-phase fluid pressure within the pore nodes is obtained as follows:

[0124] Based on the throat At time t, the fluid saturation is used to identify the connected clusters of the gas and water phases. For each connected cluster of phases, nonlinear equations about the pressure at pore nodes are assembled. The nonlinear equations are solved using a sparse linear algebra algorithm to obtain the α-phase fluid pressure within all pore nodes.

[0125] In the specific implementation process, firstly, based on the fluid saturation of the throat, the interconnected clusters of the gas phase and the aqueous phase are identified. The specific identification rules are as follows: if the fluid saturation in the throat is 100% gas phase, then the throat belongs to the gas interconnected cluster; if the saturation is 100% aqueous phase, then the throat belongs to the aqueous interconnected cluster; if both gas and aqueous phases coexist and neither has 100% saturation, then the throat is considered not to form a continuous single-phase pathway at the current moment.

[0126] Secondly, in the entire pore network model, the fluid flow in any pore node k satisfies the law of conservation of mass and Kirchhoff's laws. Therefore, for any pore node i, the assembled nonlinear equations are as follows:

[0127]

[0128] in, Represents the cluster of neighboring nodes of pore node i; k is... Any pore node in the middle; Δt is the time step; For the pressure at pore node i The change value within; It is the fluid elasticity quantity; ; The volume of the pore nodes; The combined compressibility coefficient of gas and water is calculated using the following formula:

[0129]

[0130] in, The saturation level of the aqueous fluid; This refers to the saturation level of the gas phase fluid. is the compressibility coefficient of water; The compressibility coefficient of gas; This is the ratio between the location of the gas-water interface and the length of the throat. .

[0131] Finite volume discretization on volume element i is performed using backward Euler time discretization (t is the current time, t+1 is the next time). This can be converted into the following formula:

[0132]

[0133] in, Let α be the α-phase fluid elastic quantity at pore node i; Δt is the time step. Let be the pressure of the α-phase fluid inside pore node i at time t+1; Represents the cluster of neighboring nodes of pore node i; k is... Any pore node in the middle; throat The saturation of the α-phase fluid at time t; For the throat The conductivity of the α-phase fluid at time t; Let be the pressure of the α-phase fluid within pore node j at time t; For the throat The capillary pressure at time t; For the throat The directionality of capillary pressure of the α-phase fluid at time t.

[0134] Then, combining the nonlinear equations, for each identified connected cluster of phases, nonlinear equations regarding pore nodal pressures are assembled, so the above equation can be written as:

[0135]

[0136] in, This is a vector consisting of the pressures of the α-phase fluid within all the pore nodes to be solved at the next moment; The coefficient matrix is ​​given by the following formula:

[0137]

[0138] For the right-hand side, the formula is as follows:

[0139]

[0140] Finally, a sparse linear algebra algorithm is used to solve the nonlinear equations, thereby obtaining... That is, the distribution of α-phase fluid pressure within all pore nodes, including and .

[0141] The following describes the overall process of the gas-water two-phase flow simulation method provided in this application.

[0142] I. Model initialization;

[0143] ① Initialize the model by generating the pore-throat topology and initial geometric parameters using computer speech or by importing experimental data: , , , , , ;

[0144] ②Setting fluid properties Contact angle Wettability and initial saturation;

[0145] ③ Set the aperture evolution coefficient and reference status;

[0146] ④ Specify the boundary conditions for gas-water two-phase flow (constant pressure boundary, constant flow boundary, inlet gas phase saturation, and outlet pressure).

[0147] II. Update throat radius and physical properties;

[0148] ①From the previous step , , and renew ;

[0149] ②by renew , , , and .

[0150] III. Node cluster identification and nonlinear equation assembly;

[0151] ①Based on the current fluid saturation level in the larynx Identify the connected clusters of the gas and water phases respectively;

[0152] ② Assemble the connected clusters of each phase separately:

[0153]

[0154] IV. Solve for fluid pressure and flow rate;

[0155] ① Solve using sparse linear algebra ;

[0156] ② Calculate the number of larynxes. .

[0157] V. Interface advancement and flow state handling;

[0158] ① Calculate the location X of the gas-water interface. ij and fluid volume fraction Update the interface position and advance the interface.

[0159] ②If the fluid pressure difference Greater than capillary pressure This triggers the intrusion process, updating the fluid saturation of the throat. Corner half angle ;

[0160] ③ Update the local water saturation of the larynx for the next step. .

[0161] VI. Convergence and Time Progression;

[0162] ① Check the residuals against the global mass conservation error; if the termination condition is met, output the simulation results; otherwise, output the maximum value. Constraint adaptive adjustment Then return to step II.

[0163] The following specific embodiments illustrate the gas-water two-phase flow simulation method provided in this application.

[0164] Example 1:

[0165] Construct a pore network model with a circular cross-section of 100×100×200 pore nodes, as follows: Figure 2 As shown, the average length of the larynx is L = 200 μm, and the initial equivalent radius of the larynx follows... Fluid properties, gas phase viscosity, are taken as... The viscosity of the aqueous phase is taken as interfacial tension Take 0.072 N / m, and take the wetting contact angle as... The boundary conditions are constant gas volumetric flow rate at the inlet and constant pressure at the outlet, and the model is initially fully water-saturated.

[0166] Considering the multiplicative aperture evolution, the stress sensitivity coefficient is taken as... The rate coefficient of dissolution / swelling is And adopt the effective stress Biot coefficient The value is 0.7, and the simulation results are as follows: Figure 3 As shown, the pore network model with static pore size is as follows: Figure 4 As shown, numerical results considering aperture evolution show that the overall invasion threshold distribution shifts to the right (making invasion more difficult), the finger fractal dimension decreases, the residual saturation of the water phase increases by about 8-12%, and the breakthrough time is delayed by about 15-20%.

[0167] Example 2:

[0168] Add triangular / square throats to the network and introduce a shape factor. and corner half angle The established porous network model is as follows: Figure 5 As shown, the nodes are 30×30×30. Chemical precipitation is... Characterization, The sedimentation equivalent geometric rate per unit time; It is the surface chemical reaction rate constant; C represents the initial phase concentration. This represents the phase concentration at equilibrium. This is the geometric amplification factor for precipitation. , The "positive part" operator indicates that the operator is used only when there is oversaturation. Precipitation occurs; a value of 0 is taken when undersaturated. Simulation results are as follows: Figure 6 As shown, the corner membrane maintains water connectivity and enhances capillary effect; precipitation has a selective "passivation" effect on pore throats, such as... Figure 7 As shown, the macroscopic manifestations are an increase in capillary diameter in the early stage and a stabilization of the gas connection path in the later stage.

[0169] Based on the same inventive concept, this application also provides a gas-water two-phase flow simulation system based on pore size evolution, the system comprising:

[0170] The simulation module is used to discretize the porous medium into a pore network model consisting of pore nodes and connecting throats, and to define geometric parameters for each throat. The geometric parameters include length, equivalent radius, cross-sectional area, perimeter, corner half-angle, and shape factor.

[0171] The aperture evolution update module is used to update the equivalent radius of the throat at the next time step based on the aperture evolution equation. The aperture evolution equation describes the evolution of the equivalent radius of the throat driven by mechanical, hydration, and chemical processes. Based on the equivalent radius of the throat at the next time step, the module updates the cross-sectional area, perimeter, and shape factor of the throat at the next time step. Based on the length of the throat, the cross-sectional area, and shape factor of the throat at the next time step, the module updates the single-phase fluid conductivity of the throat at the next time step. Based on the corner half-angle of the throat, the equivalent radius, and shape factor of the throat at the next time step, the module updates the capillary pressure of the throat at the next time step.

[0172] The interface propulsion module is used to solve the fluid pressure field and flow field that satisfy the law of conservation of mass and Kirchhoff's equations in the pore network model based on the single-phase fluid conductivity and capillary pressure of the throat at the next moment, and propel the gas-water two-phase interface until the termination condition is met, and output the simulation results of the gas-water two-phase seepage model.

[0173] It should be noted that each module in the gas-water two-phase flow simulation system based on pore size evolution in this embodiment corresponds one-to-one with each step in the gas-water two-phase flow simulation method based on pore size evolution in the aforementioned embodiment. Therefore, the specific implementation of this embodiment can refer to the implementation of the aforementioned gas-water two-phase flow simulation method based on pore size evolution, and will not be repeated here.

[0174] Based on the same inventive concept, this application also provides a computer device, which includes a processor, a memory, and a computer program stored in the memory. The computer program is executed by the processor to implement the aforementioned gas-water two-phase flow simulation method based on pore size evolution.

[0175] Based on the same inventive concept, this application also provides a computer storage medium storing a computer program, which is executed by a processor to implement the aforementioned gas-water two-phase flow simulation method based on pore size evolution.

[0176] In summary, compared with the prior art, this application has the following beneficial effects:

[0177] (1) Explicit strong coupling, breaking the static assumption.

[0178] The model explicitly couples a two-way feedback loop of "pore size evolution - capillary force - fluid flow" in the pore network model, breaking through the limitations of the traditional "fixed pore throat geometry". It can dynamically characterize the pore throat reconstruction and wettability changes caused by effective stress, hydration effect and geochemical action, thus more realistically reflecting key behaviors such as fingering morphology, capillary trapping, breakthrough time and pressure drop.

[0179] (2) The prediction accuracy and verifiability are significantly improved.

[0180] Through modular evolution laws and convenient calibration of field / experimental data, the influence of aperture evolution on the invasion process, fingering mode and trapping mechanism can be effectively captured, significantly improving the prediction accuracy of indicators such as residual gas / water saturation, breakthrough time and pressure drop evolution; at the same time, it provides interpretable outputs such as invasion threshold distribution, connectivity path and trapping statistics, which facilitates consistency verification with microfluidics, core testing and in-situ imaging results.

[0181] (3) Numerical robustness and engineering reliability of multi-physics coupling.

[0182] While ensuring computational efficiency, a modular solution process and robust coupling strategy are adopted to support stable calculations under multiple time scales, strong nonlinearity, and complex wetting conditions. Through mechanisms such as event-triggered threshold capture, time step adaptation, and local incremental updates, numerical oscillations and divergence risks are reduced, and long-term simulations are stable and reliable.

[0183] (4) It is scalable, integrable and easy to implement in engineering.

[0184] The modular design of the solution process and evolution law in this invention facilitates the integration of field monitoring and experimental calibration data, supports complex geometries and wetting features such as non-circular cross sections and corner films, and is easily extended to different reservoir types and operating conditions. At the same time, it natively supports uncertainty quantification and sensitivity analysis, outputs confidence intervals and key parameter rankings, and improves the credibility of scheme design and risk assessment.

[0185] (5) High efficiency in computing and low maintenance costs.

[0186] By employing pre-computation and lookup table mechanisms, parallel pressure field solving, and on-demand hot-loading of physical / geometric databases, the system significantly reduces single-step overhead and improves computational efficiency for large-scale networks. Through the integration of methods, systems, and media, it enables cross-platform deployment and rapid migration, reduces secondary development and maintenance costs, and possesses good industrial applicability and sustainable upgrade capabilities.

[0187] Based on the same inventive concept, the present invention also provides a gas-water two-phase flow simulation system based on pore size evolution, comprising:

[0188] The initialization module is used to discretize the porous medium into a pore network model consisting of pore nodes and connecting throats, and to define geometric parameters for each throat. The geometric parameters include length, equivalent radius, cross-sectional area, perimeter, corner half-angle, and shape factor.

[0189] The first update module is used to update the equivalent radius of the throat at the next time step based on the aperture evolution equation; the aperture evolution equation is used to describe the evolution of the equivalent radius of the throat driven by mechanical, hydration and chemical effects.

[0190] The second update module is used to update the cross-sectional area, perimeter, and shape factor of the throat at the next time step based on the equivalent radius of the throat at the next time step.

[0191] The third update module is used to update the single-phase fluid conduction capacity of the throat at the next moment based on the throat length, the cross-sectional area of ​​the throat at the next moment, and the shape factor.

[0192] The fourth update module is used to update the capillary pressure of the larynx at the next time step based on the corner half angle of the larynx, the equivalent radius of the larynx at the next time step, and the shape factor.

[0193] The propulsion module is used to solve the fluid pressure field and flow field that satisfy the law of conservation of mass and Kirchhoff's equations in the pore network model based on the single-phase fluid conduction capacity and capillary pressure of the throat at the next moment, and to propel the gas-water two-phase interface until the termination condition is met, and output the simulation results of gas-water two-phase seepage.

[0194] It should be noted that each module in the gas-water two-phase flow simulation system based on pore size evolution in this embodiment corresponds one-to-one with each step in the gas-water two-phase flow simulation method based on pore size evolution in the aforementioned embodiment. Therefore, the specific implementation of this embodiment can refer to the implementation of the aforementioned gas-water two-phase flow simulation method based on pore size evolution, and will not be repeated here.

[0195] Based on the same inventive concept, this application also provides a computer device, which includes a processor, a memory, and a computer program stored in the memory. The computer program is executed by the processor to implement the aforementioned gas-water two-phase flow simulation method based on pore size evolution.

[0196] Based on the same inventive concept, this application also provides a computer storage medium storing a computer program, which is executed by a processor to implement the aforementioned gas-water two-phase flow simulation method based on pore size evolution.

[0197] In some embodiments, the computer-readable storage medium may be a memory such as FRAM, ROM, PROM, EPROM, EEPROM, flash memory, magnetic surface memory, optical disk, or CD-ROM; or it may be a device including one or any combination of the above-mentioned memories. The computer may be a variety of computing devices, including smart terminals and servers.

[0198] In some embodiments, executable instructions may take the form of a program, software, software module, script, or code, written in any form of programming language (including compiled or interpreted languages, or declarative or procedural languages), and may be deployed in any form, including as a standalone program or as a module, component, subroutine, or other unit suitable for use in a computing environment.

[0199] As an example, executable instructions may, but do not necessarily, correspond to files in the file system. They may be stored as part of a file that holds other programs or data, for example, in one or more scripts in a HyperText Markup Language (HTML) document, in a single file dedicated to the program in question, or in multiple co-located files (e.g., a file that stores one or more modules, subroutines, or code sections).

[0200] As an example, executable instructions can be deployed to execute on a single computing device, or on multiple computing devices located in one location, or on multiple computing devices distributed across multiple locations and interconnected via a communication network.

[0201] It should be noted that, in this document, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or system that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or system. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or system that includes that element.

[0202] The sequence numbers of the embodiments in this application are for descriptive purposes only and do not represent the superiority or inferiority of the embodiments.

[0203] The above specific embodiments further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above are merely specific embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for simulating gas-water two-phase flow based on pore size evolution, characterized in that, include: The porous medium is discretized into a pore network model consisting of pore nodes and connecting throats, and geometric parameters are defined for each throat. The geometric parameters include length, equivalent radius, cross-sectional area, perimeter, corner half angle, and shape factor. Based on the aperture evolution equation, the equivalent radius of the throat is updated at the next time step; the aperture evolution equation is used to describe the evolution of the equivalent radius of the throat driven by mechanical, hydration and chemical effects. Update the cross-sectional area, perimeter, and shape factor of the throat at the next time step based on the equivalent radius of the throat at the next time step. The single-phase fluid conduction capacity of the throat at the next time step is updated based on the length of the throat, the cross-sectional area of ​​the throat at the next time step, and the shape factor. Update the capillary pressure of the larynx at the next time step based on the corner half angle of the larynx, the equivalent radius of the larynx at the next time step, and the shape factor. Based on the single-phase fluid conductivity and capillary pressure of the throat at the next moment, the fluid pressure field and flow field that satisfy the law of conservation of mass and Kirchhoff's equations are solved in the pore network model, and the gas-water two-phase interface is advanced until the termination condition is met, and the simulation results of gas-water two-phase flow are output. The aperture evolution equation is a multiplicative equation, and its expression is as follows: ； Alternatively, the aperture evolution equation can be an additive equation, expressed as follows: ； in, The throat between adjacent pore nodes i and j The equivalent radius at any given moment; Indicates the time step; This represents the rate of change of the equivalent radius of the throat between adjacent pore nodes i and j; Let be the equivalent radius of the throat at time t; The effective stress of the throat at the current moment; The initial stress of the throat; The local water saturation of the throat at the current moment; The initial water saturation of the throat; The equivalent geometric rate of dissolution of the throat per unit time; The equivalent geometric rate of sedimentation in the throat per unit time; and The aperture evolution coefficient is determined experimentally.

2. The gas-water two-phase flow simulation method based on pore size evolution according to claim 1, characterized in that, The Calibration was performed through triaxial densification experiments, hydrostatic densification experiments, or pore pressure cyclic loading experiments; Calibration was performed using permeability or pore size imaging data under changes in moisture absorption / dehydration / water saturation; The effect of the thickness of the dissolved layer on the throat radius was experimentally calibrated; Experimental calibration of the effect of sediment thickness on throat radius variation.

3. The gas-water two-phase flow simulation method based on pore size evolution according to claim 1, characterized in that, The formula for calculating the single-phase fluid conductivity is as follows: ; in, The α-phase fluid conductivity of the throat between adjacent pore nodes i and j at time t; α is a phase characterization parameter of gas or water. , w represents gas, and w represents water; This represents the cross-sectional area of ​​the throat at time t; The shape factor of the larynx at time t; Represents a geometric function determined by the shape factor; Let be the phase viscosity of the α-phase fluid; The length of the larynx.

4. The gas-water two-phase flow simulation method based on pore size evolution according to claim 1, characterized in that, The formula for calculating the capillary pressure is as follows: ； in, The capillary pressure of the throat between adjacent pore nodes i and j at time t when gas and water coexist. The interfacial tension between air and water; Let be the equivalent radius of the throat at time t; The corner half-angle of the throat; It is the wetting contact angle between the gas-liquid two-phase interface and the wall of the throat. The shape factor of the larynx; A correction function to account for non-circular cross sections and corner film effects.

5. The gas-water two-phase flow simulation method based on pore size evolution according to claim 1, characterized in that, The process involves solving for the fluid pressure and flow fields satisfying the law of conservation of mass and Kirchhoff's equations in the pore network model based on the single-phase fluid conductivity and capillary pressure of the throat at the next moment, and advancing the gas-water two-phase interface, including: Based on the single-phase fluid conductivity and capillary pressure of the throat at the next moment, the law of conservation of mass and Kirchhoff's equations are solved to obtain the fluid flow rate of the throat. Based on the fluid flow rate and time step of the throat, calculate the fluid volume fraction of the throat and the axial advance distance of the gas-water interface in the throat. The position of the gas-water two-phase interface in the pore network model is updated based on the axial advance distance.

6. The gas-water two-phase flow simulation method based on pore size evolution according to claim 5, characterized in that, The law of conservation of mass and Kirchhoff's equations are as follows: ； in, The α-phase fluid flow rate is the throat between adjacent pore nodes i and j. The α-phase fluid saturation of the throat; The α-phase fluid conduction capacity of the throat; The pressure of the α-phase fluid within pore node i; The pressure of the α-phase fluid within pore node j; This indicates the directionality of the capillary pressure of the α-phase fluid in the throat. The capillary pressure of the larynx.

7. The gas-water two-phase flow simulation method based on pore size evolution according to claim 6, characterized in that, The α-phase fluid pressure within the pore node is obtained in the following manner: Based on the fluid saturation of the throat at time t, the interconnected clusters of the gas phase and the water phase are identified respectively; For each phase of the connected cluster, nonlinear equations concerning the pore node pressure are assembled separately; The nonlinear equations were solved using a sparse linear algebra algorithm to obtain the α-phase fluid pressures within all pore nodes; the nonlinear equations are as follows: ； ； ； in, This is a vector consisting of the α-phase fluid pressures within all the pore nodes to be solved; It is a coefficient matrix; The right-hand term; Let α be the α-phase fluid elasticity of pore node i; Δt is the time step. Let be the pressure of the α-phase fluid inside pore node i at time t+1; Represents the cluster of neighboring nodes of pore node i; k is... Any pore node in the middle; The α-phase fluid saturation of the throat at time t; The α-phase fluid conductivity of the throat at time t; Let be the pressure of the α-phase fluid within pore node j at time t; The capillary pressure of the larynx at time t; The directionality of the capillary pressure of the α-phase fluid in the throat at time t.

8. The gas-water two-phase flow simulation method based on pore size evolution according to claim 5, characterized in that, The time step satisfies: Where Comax is a preset time step limit index, with a value range of [0.1, 0.5]. The α-phase fluid flow rate of the throat; For time step; Let be the cross-sectional area of ​​the throat at time t; max() represents the maximum value function.

Images