A Numerical Simulation Method and System for Multiphase and Multicomponent Fluid-Structure Coupling
By employing a multiphase, multicomponent fluid-structure interaction numerical simulation method, combining thermodynamic phase equilibrium and Biot's porosity elasticity theory, mechanical and mass transport equations are constructed and iteratively solved using the fully implicit Newton method. This approach solves the problem of efficient and stable numerical simulation of multiphase changes and rock deformation in unconventional oil and gas reservoirs, and achieves high-precision fluid-structure interaction processes.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- CENT SOUTH UNIV
- Filing Date
- 2026-07-02
- Publication Date
- 2026-07-31
AI Technical Summary
Existing numerical simulation techniques cannot accurately and stably reproduce the real-time coupling process of complex phase transitions of multiphase and multicomponent fluids and deformation of multiscale heterogeneous rocks while ensuring computational timeliness. In particular, in the development of unconventional and complex oil and gas reservoirs such as deep shale oil and gas, tight sandstone oil and gas, and condensate gas reservoirs, traditional methods suffer from problems such as large computational load, non-convergence of calculations, and memory overflow.
A multiphase, multicomponent fluid-structure interaction numerical simulation method is adopted. Based on the thermodynamic phase equilibrium theory and Biot's porosity elasticity theory, mechanical equilibrium equations and mass transport control equations are constructed. Nonlinear iterative solutions are performed using the fully implicit Newton method. Combined with the fully coupled solution format, the deformation of the rock mass skeleton, fluid flow and component content are solved iteratively at the same time. Jacobi matrix is constructed to perform cross-partial derivative calculations, so as to achieve efficient flash evaporation calculation and high-precision phase distribution.
It achieves high-precision and rapid convergence under strong nonlinear coupling conditions, accurately predicts the dynamic production characteristics of unconventional oil and gas reservoirs, improves computational stability and efficiency, and solves the problems of large computational load and non-convergence in traditional methods.
Smart Images
Figure CN122491160A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the technical field of oil and gas field development engineering and computational fluid dynamics, and in particular to a numerical simulation method and system for multi-component fluid-structure interaction. Background Technology
[0002] In the development of unconventional and complex oil and gas reservoirs such as deep shale oil and gas, tight sandstone oil and gas, and condensate gas reservoirs, there is a strong interaction between formation fluid flow and reservoir rock deformation. As wellhead fluids are extracted, the pore pressure within the porous medium continuously decreases, leading to an increase in effective stress, which in turn triggers formation rock deformation and fracture closure. Simultaneously, drastic changes in reservoir pressure can also cause complex phase changes in multi-component fluids, such as gas condensation, component fractionation, and supercritical fluid transformation in shale oil. This complex interplay of multiphase, multi-component phase transitions and nonlinear rock mechanics makes accurate numerical simulation extremely difficult.
[0003] Currently, although existing technologies have made significant progress in both fluid-structure interaction and multiphase multicomponent reservoir simulation, when deeply and organically integrating the two, the following technical challenges remain to be addressed in the face of the unique multi-scale and strongly nonlinear characteristics of unconventional reservoirs:
[0004] First, there is the "computational disaster" caused by the interplay between high-dimensional component flash evaporation calculations and nonlinear mechanics. Traditional black oil models cannot characterize the fractionation phenomena and gas condensation behavior of multiphase and multicomponent fluids under deep, high-temperature, and high-pressure conditions. Although introducing component models can solve phase transition problems, the computational load increases exponentially with the number of components and the grid size because high-dimensional phase equilibrium flash evaporation calculations are required for each time step and each grid of the multiphase and multicomponent fluid flow. When this complex phase equilibrium calculation is combined with the rock pore elastic stress equation, which also has strong nonlinearity, the Jacobian matrix required by the fully implicit method becomes extremely large and the condition number is extremely poor. This causes existing simulators to easily experience computational stagnation, memory overflow, or non-convergence when dealing with large-scale real oil reservoirs.
[0005] Second, there is a dilemma between stability and timeliness in multi-field strongly coupled solution architectures. While fully implicit coupling offers good numerical stability, it suffers from catastrophic increases in memory consumption and computation time when dealing with large mesh sizes and numerous components, failing to meet engineering timeliness requirements. Existing iterative coupling methods, while reducing the dimensionality of single-step calculations, exhibit extremely strong coupling between the flow and mechanics equations when facing highly nonlinear problems such as fracture closure and large, non-uniform deformation of reservoirs. In the early stages of mining under drastic pressure and stress changes, iterative solutions often require numerous internal loops to converge, frequently resulting in non-convergence and continuous data oscillations. There is a lack of semi-implicit or spatiotemporally split solution architectures that can guarantee absolute convergence while achieving efficient computation.
[0006] In summary, current numerical simulation techniques cannot accurately and stably reproduce the entire real-time coupling process of complex phase transitions in multiphase and multicomponent fluids and deformation of multiscale heterogeneous rocks while ensuring computational timeliness. Summary of the Invention
[0007] The purpose of this invention is to overcome the shortcomings of the prior art. This invention provides a multiphase and multicomponent fluid-structure interaction numerical simulation method and system. The method aims to accurately and stably reproduce the real-time coupling process of complex phase transformation of multiphase and multicomponent fluids and deformation of multiscale heterogeneous rocks while ensuring computational timeliness, thereby improving the solution efficiency of phase state calculation and fluid flow.
[0008] To achieve the above-mentioned technical objectives, the present invention provides the following technical solution: In a first aspect, the present invention provides a multiphase, multicomponent fluid-structure interaction numerical simulation method, comprising: S1: Construct a geometric model based on the geological data of the target research object and perform grid generation; S2: Construct the governing equations of the geometric model after mesh generation; wherein, the governing equations include mechanical equilibrium equations characterizing the deformation behavior of the porous media rock mass skeleton, and mass transport governing equations characterizing the flow and material transport behavior of multiphase and multicomponent fluids in porous media; S3: Based on the thermodynamic phase equilibrium theory, flash evaporation calculations are performed on each grid computing unit to determine the phase distribution and parameters of each component of the multiphase multicomponent fluid under the current thermodynamic state; S4: Construct a fully coupled solution scheme for the discretized control equations. Incorporate displacement variables representing rock mass skeleton deformation, pressure variables representing fluid flow, and component variables representing the content of each component as master variables to be solved into a unified solution framework. Simultaneously iterate and solve the master variables to be solved to obtain the numerical simulation results at the current time step.
[0009] Furthermore, the mechanical equilibrium equation in S2 is derived from Biot's porosity elasticity theory; The mechanical equilibrium equation is: ; in, ; For Cauchy stress tensor; For rock mass density, , = The total density of the fluid. for Phase saturation, for Phase density, This represents the current porosity of the rock mass. Density of the rock; It is the vector of gravitational acceleration; Cauchy stress tensor Effective stress of porous media rock mass skeleton The relationship is: ; The effective stress of the porous media rock mass skeleton was obtained by sorting out the data. for: ; in, It is a fourth-order elastic tensor; For small strain tensors, For solid displacement vectors; It is a second-order unit tensor; For Biot coefficient, The bulk elastic modulus of the rock mass skeleton. The bulk modulus of the solid particle matrix itself; Equivalent pore pressure for multiphase fluids; The formula for calculating the equivalent pore pressure of the multiphase fluid is as follows: ; In the formula: for Phase pressure.
[0010] Furthermore, the mass transport control equation in S2 is constructed by using the molar concentration of each component as a conserved variable; for any component The mass transport control equation is: ; in, This represents the current porosity of the rock mass. for The molar density of the phase; Components exist The mole fraction in the phase; for Phase saturation; For the body strain; for The relative Darcy velocity vector of the phase; Components exist Hydrodynamic dispersion-diffusion tensor in the phase; For source / sink; The number of components; The The relative Darcy velocity vector of the phase is: ; in, For the absolute permeability tensor of the rock mass; for The relative permeability of the phase; for Phase viscosity; for Phase pressure of the phase; for Average molar mass of the phase, , Components The molecular weight.
[0011] Furthermore, the flash evaporation calculation in S3 includes: For a fluid mixture containing multiple components, an equation of state is constructed; wherein the equation of state is: ; in, Equivalent pore pressure for multiphase fluids; This is the universal gas constant; Reservoir temperature; This represents the molar volume of the mixture. The energy parameters of the mixture; This refers to the volume parameters of the mixture; Based on the constructed equation of state, the fugacity coefficients of each component are calculated; exist Fugacity coefficient in phase satisfy: ; in, Components Volume parameters; for Phase compressibility factor for molar volume of the phase; and All are dimensionless state equation parameters; Components exist The mole fraction in the phase; Components and The interaction coefficient between them; Construct and solve the phase equilibrium equations The phase distribution and parameters of each component under the current thermodynamic state are obtained; among them, the phase equilibrium equation is... for: ; in, Components within a grid computing cell Total mole fraction; Components The phase equilibrium constant, , Components Mole fraction in the gas phase Components Mole fraction in the liquid phase; This refers to the number of components (i.e., the total number of components in the mixture).
[0012] Furthermore, the flash evaporation calculation also includes: Solving the phase equilibrium equation Previously, the composition was based on Wilson's empirical formula. Initialize the phase equilibrium constant: ; in, Components Critical pressure; Components The critical temperature; Components eccentricity factor; Solving the phase equilibrium equation After obtaining the phase distribution and parameters of each component under the current thermodynamic state, the fugacity coefficients are calculated, and the ratios of the calculated fugacity coefficients are used to assess the composition. The phase equilibrium constant is iteratively updated: ; in, Components Fugacity coefficient in the liquid phase; Components In the gas phase fugacity coefficient.
[0013] Repeatedly solve the phase equilibrium equations And update the components The phase equilibrium constant is calculated until the phase equilibrium convergence condition of the fugacity equality criterion is met; wherein, the fugacity equality criterion is: ; in, Components Fugacity in the gas phase ; Components Fugacity in the liquid phase ; This is the preset convergence tolerance.
[0014] Furthermore, it also includes updating the porosity parameters and permeability tensor of the porous medium based on the numerical simulation results at the current time step; The porosity parameter is based on volumetric strain. The updated result is: ; in, The initial porosity of the rock mass; For volumetric strain, For small strain tensors, For solid displacement vectors; This represents the current porosity of the rock mass. The bulk modulus of the solid particle matrix itself; for Phase saturation; for Phase pressure of the phase; The permeability tensor It is obtained by updating the Carman-Kozeny equation based on anisotropic correction: ; in, This represents the initial permeability tensor of the rock mass. For the current volumetric strain Porosity below; It is a second-order unit tensor.
[0015] Furthermore, it also includes: calculating the volumetric load of the rock mass skeleton based on the numerical simulation results of the current time step, according to the phase saturation distribution and pressure distribution of the multiphase fluid, and feeding the volumetric load as an external force term in the force balance equation to the rock mass deformation solution of the next time step; the volumetric load for: .
[0016] Furthermore, in S4, the discretized control equations include: The discretized mechanical equilibrium equations are: ; in, It is the volume integral domain; For effective stress tensor; For strain tensor; For the boundary; For the boundary The known surface force tensor on; Density of the rock mass; This is a test function for rock mass displacement. , For the displacement-dependent first-order homogeneous Sobolev space; Discretized mass transport control equations: ; in, for The molar density of the phase; Components exist The mole fraction in the phase; for Phase saturation; for Average mobility of the phase; For component k in The diffusion coefficient tensor in the phase; For source / sink; Scalar test function for component molar concentration. , Let be a first-order homogeneous Sobolev space related to the mole fraction.
[0017] Furthermore, in S4, the fully coupled solution scheme employs a fully implicit Newton method for nonlinear iterative solution: The residual vector after spatial discretization is constructed as follows: ; in, It is a vector transpose operator; It is the residual vector; The residual vector of the mechanical equilibrium equations: ; Let J be the residual vector of the mass conservation equation for the j-th component: ; The iterative solution involves: in each Newton iteration step, constructing the Jacobian matrix of the discretized control equation with respect to the main variables to be solved; and then updating the main variables to be solved synchronously based on the Jacobian matrix until the convergence criterion is met. The Jacobian matrix includes several cross-partial derivative terms between rock mass skeleton deformation, fluid flow, and component content variations: ; ; in, A vector composed of the molar numbers of different components; For displacement; This represents the current porosity of the rock mass. For rock mass permeability tensor; for Phase saturation; for The molar density of the phase; Components exist The mole fraction in the phase; The number of components, and These are the residual vectors of the mechanical equilibrium equations and the first... The residual vector of the mass conservation equation for each component can be expressed as: ; ; Secondly, the present invention provides a multiphase, multicomponent fluid-structure interaction numerical simulation system, the system being used to perform the steps of the method described above, including: Geometric model and mesh generation module: used to construct a geometric model based on the geological data of the target research object and to perform mesh generation; Equation Construction Module: Used to construct the governing equations of the geometric model after mesh generation; wherein, the governing equations include mechanical equilibrium equations characterizing the deformation behavior of the porous media rock mass skeleton, and mass transport control equations characterizing the flow and material transport behavior of multiphase and multicomponent fluids in porous media; Phase equilibrium calculation module: Based on thermodynamic phase equilibrium theory, it performs flash evaporation calculations on each grid calculation unit to determine the phase distribution and parameters of each component of the multiphase multicomponent fluid under the current thermodynamic state; Coupled solution module: Constructs a fully coupled solution format for the discretized control equations, incorporating displacement variables characterizing rock mass skeleton deformation, pressure variables characterizing fluid flow, and component variables characterizing the content of each component as master variables to be solved into a unified solution framework, and synchronously iterating to solve the master variables to be solved to obtain the numerical simulation results of the current time step.
[0018] This invention proposes a multiphase, multicomponent fluid-structure interaction numerical simulation method and system. The method first constructs the rock mass skeleton mechanical equilibrium equations and combines them with multi-component mass conservation and multiphase, multicomponent mass transport control equations to comprehensively characterize the multi-physics interaction behavior. Then, by introducing a thermodynamic phase equilibrium and efficient flash vaporization calculation model, fluid properties and phase distribution are dynamically updated, significantly reducing the dimensionality and frequency of high-dimensional phase equilibrium calculations and solving the problem of excessive computational load in multi-component models. For spatial discretization, the controlled volume finite element method is used to perform high-precision geometric discretization of complex fractures and unstructured meshes, maintaining local mass conservation. For computational coupling, it breaks through the limitations of traditional loose iteration and adopts a fully implicit Newton method nonlinear solution strategy, incorporating displacement, pressure, and component content into a unified equation set. Synchronous iteration is performed by accurately constructing a Jacobian matrix containing full cross derivatives, ensuring the numerical convergence of strong nonlinear coupling while achieving rapid alternating solutions for flow equations and nonlinear stress equations, thereby significantly reducing the time consumption for simulating large-scale complex fracture networks. By constructing a cross-scale fracture network fluid-structure interaction mechanical response model that considers micro-stress variations, the porosity and permeability tensors of the matrix and fractures can be updated in real time, accurately predicting the dynamic production characteristics of unconventional oil and gas reservoirs throughout their entire lifecycle. This invention achieves seamless, bidirectional, and strong coupling of energy and mass between physical fields, significantly improving convergence speed, numerical stability, and computational accuracy under strongly nonlinear conditions. Attached Figure Description
[0019] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.
[0020] Figure 1 This is a flowchart of the multiphase, multicomponent fluid-structure interaction numerical simulation method provided in this embodiment of the invention; Figure 2 This is a schematic diagram of the geometric model obtained by constructing the area where the three horizontal wells are located and the surrounding reservoirs provided in the embodiment of the present invention after meshing; Figure 3 This is a spatial distribution diagram of the well network, well trajectory, and artificial fractures provided in the embodiments of the present invention; Figure 4 This is a graph showing the change in bottom-hole flowing pressure of three wells provided in an embodiment of the present invention; Figure 5 This is a comparison chart of the cumulative gas production of three wells provided in an embodiment of the present invention; Figure 6 This is a comparison chart of the cumulative water production of three wells provided in an embodiment of the present invention; Figure 7 This is a numerical model diagram of CO2 injection for tight oil reservoirs provided in an embodiment of the present invention; Figure 8 This is a CO2 component distribution diagram at different times during the first cycle of CO2 huff and puff, provided by an embodiment of the present invention, without considering rock mass deformation; wherein, Figure 8 (a) The first cycle of gas injection has ended; Figure 8 (b) marks the end of the first cycle of well shut-in; Figure 8 (c) The first production cycle has ended; Figure 9 This is a CO2 component distribution diagram at different times during the first cycle of CO2 uptake and downtake, considering rock mass deformation, provided by an embodiment of the present invention; wherein, Figure 9 (a) The first cycle of gas injection has ended; Figure 9 (b) marks the end of the first cycle of well shut-in; Figure 9 (c) The first production cycle has ended; Figure 10 This is a comparison chart of production output considering the influence of rock mass deformation and not considering the influence of rock mass deformation, provided by an embodiment of the present invention; In the diagram: 1-overlying formation; 2-CO2 injection well; 3-large-scale artificial fracture; 4-tight oil reservoir; 5-underlying formation. Detailed Implementation
[0021] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be described in detail below. Obviously, the described embodiments are merely some embodiments of this invention, and not all embodiments. Based on the embodiments of this invention, all other implementation methods obtained by those skilled in the art without creative effort are within the scope of protection of this invention.
[0022] Example 1
[0023] like Figure 1 As shown, this embodiment provides a multiphase, multicomponent fluid-structure interaction numerical simulation method, including:
[0024] S1: Construct a geometric model based on the geological data of the target research object and perform mesh generation; the construction of the geometric model and mesh generation are existing technologies, and therefore will not be described in detail. In specific implementation, the appropriate method can be selected based on actual needs.
[0025] S2: Construct the governing equations of the geometric model after mesh generation; wherein, the governing equations include the mechanical equilibrium equations characterizing the deformation behavior of the porous media rock mass skeleton, and the mass transport governing equations characterizing the flow and material transport behavior of multiphase and multicomponent fluids in the porous media.
[0026] Specifically, the saturated fluid-laden rock is considered as a continuum composed of a rock skeleton and a multiphase fluid. In this continuum, the rock and fluid are spatially superimposed. It is assumed that the fluid is slightly compressible, the rock deformation is linearly elastic, and the rock deformation is a quasi-static process that satisfies the small strain assumption. Throughout the process, the system remains at a constant temperature and no chemical reaction occurs. Based on the conservation of momentum and the theory of porosity elasticity, the mechanical equilibrium equation for the deformation behavior of the porous media rock skeleton can be expressed in the form of the Cauchy stress tensor: ; in, ; For Cauchy stress tensor; For rock mass density, , = The total density of the fluid. for Phase saturation, for Phase density, This represents the current porosity of the rock mass. Density of the rock; Let be the gravitational acceleration vector. To connect macroscopic solid mechanics with microscopic porous fluid behavior, the effective stress principle from Biot's porosity elasticity theory is introduced.
[0027] Cauchy stress tensor Effective stress of porous media rock mass skeleton The relationship is: ; The effective stress of the porous media rock mass skeleton was obtained by sorting out the data. for: ; in, It is a fourth-order elastic tensor; For small strain tensors, For solid displacement vectors; It is a second-order unit tensor; For Biot coefficient, The bulk elastic modulus of the rock mass skeleton. The bulk modulus of the solid particle matrix itself; This refers to the equivalent pore pressure of a multiphase fluid. In practical implementation, the equivalent pore pressure of a multiphase fluid... The definition strictly follows the principle of energy conjugation: ; In the formula: for Phase pressure.
[0028] Preferably, in conventional industrial simulations, its first-order approximation, namely the saturation-weighted average form, is typically used: .
[0029] Specifically, for any component The conservation equations must be strictly maintained in a dynamic grid system where the fluid phase changes drastically and the solid skeleton deforms over time. The mass transport governing equations are constructed using the molar concentrations of each component as conservation variables. For any component... The complete continuity equation corresponding to the mass transport control equation is written as: ; in, This represents the current porosity of the rock mass. for The molar density of the phase; Components exist The mole fraction in the phase satisfies the normalization constraint. ; for Phase saturation; For the body strain; Components exist Hydrodynamic dispersion-diffusion tensor in the phase; For source / sink; The number of components; for The relative Darcy velocity vector of the phase, in its rigorous form, takes into account gravity and capillary pressure effects: ; in, For the absolute permeability tensor of the rock mass; for The relative permeability of the phase (its magnitude is mainly determined by...) (Determined by the saturation of the phase). for Phase viscosity; for Phase pressure of the phase; for Average molar mass of the phase, , Components The molecular weight.
[0030] S3: Based on the thermodynamic phase equilibrium theory, flash evaporation calculations are performed on each grid computing unit to determine the phase distribution and parameters of each component of the multiphase multicomponent fluid under the current thermodynamic state.
[0031] Specifically, phase saturation in macroscopic flow transport equations and phase component mole fraction They are not independent control variables; they are entirely determined by the local thermodynamic state. The closure of the continuum system consisting of the rock skeleton and multiphase fluid is achieved by solving a set of nonlinear algebraic equations (i.e., flash calculation) at each spatial Gaussian point (or control volume center).
[0032] This embodiment provides two methods for calculating fugacity. One method uses the Peng-Robinson (PR) equation of state to calculate the physical properties of the gas, liquid, and supercritical phases. For mixtures containing multiple components, the cubic equation of state is expressed as: ; in, Equivalent pore pressure for multiphase fluids; This is the universal gas constant; Reservoir temperature; This represents the molar volume of the mixture. The energy parameter of the mixture, Let be the volume parameter of the mixture. , It is through the critical parameters of each pure component ( , ) and eccentricity factor ( ) is obtained by pairing and combining binary interaction parameters, where and Components The critical pressure and critical temperature.
[0033] Based on the constructed equation of state, the fugacity coefficients of each component are calculated; exist Fugacity coefficient in phase satisfy: ; in, Components Volume parameters; for Phase compressibility factor for molar volume of the phase; and All are dimensionless state equation parameters; Components exist The mole fraction in the phase; Components and The interaction coefficients between them; the fugacity of the components is expressed as .
[0034] Because the Peng-Robinson equation of state cannot accurately calculate the molar ratio of the gas and liquid phases when solving for fugacity, nor can it determine whether it is a single-phase gas, a single-phase liquid, or a gas-liquid coexistence under the current pressure conditions—that is, it lacks "material balance constraints" and "phase stability analysis"—it cannot independently solve the phase determination and material distribution problems in multiphase coexistence. Therefore, the phase balance equation is transformed into the solution of the classic Rachford-Rice equation, assuming the composition within a certain grid computational cell... The total number of moles is The total volume is The corresponding total molar concentration is: ; in, Components within a grid computing cell Total mole fraction; denoted as the total number of moles of the j-th component.
[0035] Define the phase equilibrium constant (i.e., the equilibrium ratio) as follows: And introduce the mole fraction of the gas phase as Solve the Rachford-Rice equation: ; The phase distribution and parameters of each component under the current thermodynamic state are obtained, where, This represents the total number of moles in the gas phase. This represents the total number of moles in the liquid phase. Components The phase equilibrium constant, , Components Mole fraction in the gas phase Components Mole fraction in the liquid phase.
[0036] The detailed solution steps are as follows: Step 1: Initialize the balance ratio That is, assigning initial values to each component using Wilson's empirical formula. value: ; Step 2: Solve the Rachford-Rice equations, i.e., based on the initial conjecture obtained. Value, solve the Rachford-Rice equation .because about It is monotonically decreasing, and in engineering, the one-dimensional Newton-Raphson method is used to solve it. set up The initial guess value (in this embodiment, we take) =0.5), calculate and its related derivative : ; Update the mole fraction of the gas phase under the current thermodynamic state. : ; in, This represents the updated mole fraction of the gas phase; This represents the mole fraction of the gas phase before the update.
[0037] when At this point, the Rachford-Rice equation converges, yielding the mole fraction of the gas phase under the current thermodynamic state. .
[0038] Step 3: State update and fugacity coefficient calculation. Based on the mole fraction of the gas phase under the current thermodynamic state obtained in Step 2, calculate the instantaneous mole fraction of the components in both the gas and liquid phases: : : The current phase composition ( or Along with pressure Substituting temperature T into the Peng-Robinson equation of state, and then solving this cubic equation, the compressibility factors of the gas and liquid phases are obtained. and Then based on the fugacity coefficient The expression is used to calculate the fugacity coefficients of each component in the gas and liquid phases. and Furthermore, the true fugacity of each component can be calculated: ; ; Step 4: Physical equilibrium convergence determination. According to the requirement of ideal thermodynamic equilibrium, the fugacity of each component in the two phases must be completely equal. First, the fugacity residual is defined. : ; Repeat steps two and three until the fugacity residuals satisfy the convergence criterion. < Tol (suggested to use Tol= to If the fugacity calculation has reached convergence, then the flash calculation of the current grid cell is closed, and the saturation of each phase is output to the flow transport control equation. ,density and the mole fraction of each phase.
[0039] Step 5: Balance ratio replacement update. If step 4 does not converge, then explicitly update the balance ratio using the currently calculated fugacity coefficient ratio. : ; in, Components Fugacity coefficient in the liquid phase; Components In the gas phase fugacity coefficient; Components Fugacity in the gas phase; Components Fugacity in the liquid phase.
[0040] Repeat steps two through four until the fugacity residual convergence criterion is met in step four, at which point the calculation stops.
[0041] Fluid-structure interaction calculation: Solid deformation leads to changes in the geometry of micropores in the porous media framework. The porosity change caused by rock framework deformation can be calculated by the following formula: ; In the formula: The initial porosity of the rock mass. This refers to volumetric strain. Deformation of the rock skeleton causes changes in pore space and also alters the permeability tensor of the rock mass. In this embodiment, the anisotropically modified Carman-Kozeny equation is used to describe the change in permeability caused by rock mass deformation as follows: ; in, This represents the initial permeability tensor of the rock mass. For the current volumetric strain Porosity below; It is a second-order unit tensor.
[0042] Conversely, changes in spatial saturation and fluid pressure caused by fluid flow will be expressed as volumetric loads. The action of the rock mass reacts on the rock mass framework, causing corresponding rock mass deformation: .
[0043] S4: Construct a fully coupled solution scheme for the discretized control equations. Incorporate displacement variables representing rock mass skeleton deformation, pressure variables representing fluid flow, and component variables representing the content of each component as master variables to be solved into a unified solution framework. Simultaneously iterate and solve the master variables to be solved to obtain the numerical simulation results at the current time step.
[0044] To facilitate solving on a computer, we first construct a weak form of the fully coupled system of equations.
[0045] Selecting the test function for rock mass displacement and scalar test functions for component molar concentrations The virtual work principle equation of the rock mass skeleton mechanical equilibrium equation requires... This makes it possible for all All satisfy: ; in, It is the volume integral domain; For effective stress tensor; For strain tensor; This represents the current porosity of the rock mass. For the boundary; For the boundary The known surface force tensor on; Density of the rock mass; This is a test function for rock mass displacement. , For the first-order homogeneous Sobolev space of displacement: ; ; in, For the number of nodes, For example, the shape function of a hexahedral element can be expressed as: ; Let be the natural coordinates of any point within the cell, taking values between [-1, 1]. Let be the specific coordinates of the i-th node in the natural coordinate system (takes a value of +1 or -1), that is: ; ; ; ; Weak form requirements for the control equations of multiphase and multicomponent transport using the control volume finite element method. This makes it possible for all All satisfy: ; in, for The molar density of the phase; Components exist The mole fraction in the phase; for Phase saturation; for Average mobility of the phase; For component k in The diffusion coefficient tensor in the phase; For source / sink; Scalar test function for component molar concentration. , For the first-order homogeneous Sobolev space related to the mole fraction: ; ; in, Number of units It is a shape function. When i is the same as the cell number it belongs to, it takes the value 1; otherwise, it takes the value 0.
[0046] More specifically, the fully coupled solution scheme employs a fully implicit Newton method for nonlinear iterative solution: The residual vector after spatial discretization is constructed as follows: ; in, It is a vector transpose operator. It is the residual vector; The residual vector of the mechanical equilibrium equations: ; Let be the residual vector of the mass conservation equation for the k-th component; ; In each Newton iteration step, a Jacobian matrix of the discretized governing equations with respect to the principal variables is constructed. Then, the principal variables are synchronously updated based on the Jacobian matrix until the convergence criterion is met. The Jacobian matrix contains several cross-partial derivative terms between rock mass skeleton deformation, fluid flow, and component content changes. In the fully implicitly coupled solution, the convergence of the Newton iteration strongly depends on the accuracy of the Jacobian matrix, as flash evaporation calculations lead to phase saturation... Phase density and mole fraction of component k Become a key variable with respect to independent main variables (such as total component molar density). and fluid pressure For highly nonlinear implicit functions, the construction of their partial derivatives introduces the total differential chain rule: ; ; in, For the first The residual vector of the mass conservation equations for each component; A vector composed of the molar numbers of different components; For displacement; This represents the current porosity of the rock mass. For rock mass permeability tensor; for Phase saturation; for The molar density of the phase; Components exist The mole fraction in the phase; The number of components.
[0047] Example 2
[0048] This embodiment provides a multiphase, multicomponent fluid-structure interaction numerical simulation system, which performs the steps of the method described above, including: Geometric model and mesh generation module: used to construct a geometric model based on the geological data of the target research object and to perform mesh generation; Equation Construction Module: Used to construct the governing equations of the geometric model after mesh generation; wherein, the governing equations include mechanical equilibrium equations characterizing the deformation behavior of the porous media rock mass skeleton, and mass transport control equations characterizing the flow and material transport behavior of multiphase and multicomponent fluids in porous media; Phase equilibrium calculation module: Based on thermodynamic phase equilibrium theory, it performs flash evaporation calculations on each grid calculation unit to determine the phase distribution and parameters of each component of the multiphase multicomponent fluid under the current thermodynamic state; Coupled solution module: Constructs a fully coupled solution format for the discretized control equations, incorporating displacement variables characterizing rock mass skeleton deformation, pressure variables characterizing fluid flow, and component variables characterizing the content of each component as master variables to be solved into a unified solution framework, and synchronously iterating to solve the master variables to be solved to obtain the numerical simulation results of the current time step.
[0049] It is understood that the same or similar parts in the above embodiments can be referred to each other, and the contents not described in detail in some embodiments can be referred to the same or similar contents in other embodiments.
[0050] To further verify the technical solution of this application, the following simulation application examples are provided: Application Example 1: Taking the area containing three horizontal wells (YYH1-1, YYH1-2, and YYH1-3) and their surrounding reservoirs as the research object, a geological model was constructed based on the obtained physical geometric data of the research object. A three-dimensional corner grid system was adopted, with an overall grid size of 124×190×8, totaling 188,480 grid cells. According to the corner grid cylindrical coordinate statistics, the physical range of the model plane is approximately 2234 m×3094 m, effectively covering the area containing the three horizontal wells (YYH1-1, YYH1-2, and YYH1-3) and their surrounding reservoirs. This modeling method helps to realistically describe the lag process of pressure drop propagation to the outside of the well area and the fluid recharge mechanism at the edge of the stirred volume, and it balances the accuracy of the physical mechanism with the efficiency of numerical experiments. Vertically, the model is discretized into 8 calculation segments to characterize the horizontal wellbore cross-layer relationship, the non-uniform distribution of fractures in different segments, and the vertical gravity differentiation of gas and water phases, such as... Figure 2 As shown.
[0051] On the plane, the model mesh scale is generally controlled at the 15m level, with average adjacent column spacing of approximately 14.8m and 15.0m in the x and y directions, respectively. Minor fluctuations exist locally due to structural adaptation and mesh geometric orthogonalization adjustments. This resolution can sensitively capture the leading-edge changes in pressure propagation and gas-water two-phase flow, effectively avoiding the nonlinear surge in model solution cost caused by over-refinement. Vertically, the reservoir is divided into 8 layers. Key simulation parameters include: permeability ranging from 0.01 to 0.19 mD, with an average of 0.07 mD; matrix porosity characterized within a discrete range of 0.020 to 0.054; and horizontal permeability using a 1×10⁻⁶ m² / m² mesh. -5 Discrete values such as 0.0022, 0.007, 0.010, and 0.011mD are assigned to characterize the low porosity, low permeability, and heterogeneity of the target reservoir. Vertical permeability is taken as 0.10 of horizontal permeability. The pressure gradient is approximately 0.37-1.06 MPa / 100 m, generally increasing with depth. The pressure gradient of deep reservoirs generally exceeds 0.8 MPa / 100 m. The reference depth for this example is 2100 m, and the reference pressure is set to 180 bar. All three wells are long horizontal wells with a cumulative wellbore length of approximately 1.49-1.61 km. The wellbore radius is uniformly taken as 0.0762 m, and the wellbore skin factor is set to 0. The well network, well trajectory, and spatial distribution of artificial fractures are shown in the figure. Figure 3 As shown, the bottom-hole flowing pressure of the three wells is controlled using a stepped control method, such as... Figure 4 As shown in the figure, the initial gas production was approximately 170 bar, which gradually decreased to 160, 150, 140, 127-130, and 118-120 bar, eventually stabilizing at approximately 110 bar. A comparison of the cumulative gas production of the three wells under this extraction scheme is shown in the figure below. Figure 5 As shown in the figure, the cumulative water production comparison curves of the three wells are as follows: Figure 6As shown, the cumulative gas production of the three wells during the simulation period was approximately 1.10×10⁷, 7.8×10⁶, and 7.2×10⁶ m³, respectively. 3 The cumulative water production was approximately 120, 80, and 60 m³, respectively. 3 .
[0052] Application Example 2: Taking a tight oil reservoir undergoing CO2 injection as the target research object, a geological model is constructed based on the acquired physical geometric data of the target research object, and a mesh is generated. The composition of tight oil is much more complex than that of shale gas, and its flow process is also more complex. As extraction progresses, when the pressure decreases to a level close to its two-phase region in the formation conditions, phase transitions occur within the reservoir. These phase transitions lead to mass transfer between different phases, resulting in the formation of new phases and the disappearance of certain phases in local meshes, causing changes in the dimensions of the flow matrix. The mesh model of a tight oil reservoir is shown below. Figure 7 As shown, the formation includes overlying strata 1; CO2 injection well 2; large-scale artificial fractures 3; tight oil reservoir 4; and underlying strata 5. The tight oil reservoir area contains one CO2 injection well 2 and two large-scale artificial fractures 3. To better simulate the influence of the self-weight of the overlying rocks, this example establishes an overlying rock region, a tight oil reservoir region, an underlying rock region, and a confining pressure region.
[0053] Key simulation parameters: fracture porosity 0.1, matrix porosity 0.25, initial fracture permeability 400 mD, initial matrix permeability 0.1 mD, Young's modulus 10 GPa, Poisson's ratio 0.4, Biot coefficient 0.95, rock density 2025 kg / m³ 3 The fracture wall roughness is 8, the fracture wall protrusion strength is 250 bar, the initial reservoir temperature is set to 100 K, and the initial pressure is set to 250 bar; at a reference depth of 2000 m, the corresponding reference pressure is 250 bar, and the mole fraction of the components is shown in Table 1.
[0054] Table 1 Fluid Composition Parameters of Tight Oil Reservoirs
[0055] The tight oil reservoir was developed using a CO2 huff and puff method, employing a multi-cycle process of "gas injection-well shut-in-oil production" to achieve gas dissolution, viscosity reduction, expansion, and miscible / near-miscible displacement, utilizing the remaining oil in the matrix and fractures. For the first 1000 days, production wells were designed to operate at a constant bottomhole flowing pressure of 200 bar. After 1000 days, CO2 huff and puff was used to enhance production.
[0056] The production-increasing measures will be implemented in 5 cycles, each lasting 90 days. The specific work schedule for each cycle is as follows: 1) The first 30 days are the gas injection phase, with a daily gas injection volume of 800 m³. The purpose is to inject CO2 gas into the formation to replenish energy, build pressure, achieve dissolution or miscibility, and displace crude oil. 2) After that, the well is shut down for 20 days to allow the gas and crude oil to fully dissolve, diffuse, mix, reduce viscosity, and expand, and to establish a uniform pressure field in preparation for oil production; 3) Then, production is carried out for 40 days, mainly relying on the formation elasticity + dissolved gas expansion + gas displacement to produce the tight crude oil that is used.
[0057] This example simulates four scenarios: Scenario 1: No geomechanical influence and no increase in production; Scenario 2: No geomechanical influence, production increased by gas injection; Scenario 3: Geomechanical influence considered and no increase in production; Scenario 4: Geomechanical influence considered, production increased by gas injection.
[0058] The main mechanism of CO2 huff and puff for enhanced production is that it increases reservoir pressure while simultaneously improving fluid flowability by injecting lighter components. In the first 1000 days, the fluid composition is uniformly distributed. CO2 injection alters the fluid composition, thus its effect on reservoir fluid modification can be tracked by calibrating the molar concentration of CO2. Figure 8 The concentration distribution of carbon dioxide in the component model is shown, where, Figure 8 (a) The first cycle of gas injection has ended; Figure 8 (b) marks the end of the first cycle of well shut-in; Figure 8 (c) marks the end of the first production cycle.
[0059] Figure 9 The concentration distribution of carbon dioxide in the fluid-structure interaction model is shown below. Figure 9 (a) The first cycle of gas injection has ended; Figure 9 (b) marks the end of the first cycle of well shut-in; Figure 9 (c) The first production cycle has ended. The fluid-structure interaction model shows that after CO2 injection, due to its high lateral permeability, it initially improves the fluid in the perforated layer. After well shut-in, due to its low molar mass, it experiences buoyancy in the fluid, resulting in a funnel-shaped distribution in the later stages. In the rock mass deformation model, the early fracture closure leads to low permeability and high fluid pressure in the fractured reservoir. Therefore, the CO2-stimulated area predicted in the fluid-structure interaction model is smaller than that in the component model, and the CO2 concentration in the stimulated area is higher in the fluid-structure interaction model than in the component model.
[0060] Figure 10Comparing the cumulative oil production under different conditions, the figure shows that the effect of rock mass deformation caused the production predicted by the fluid-structure interaction model to be lower than that predicted by the component model. This is consistent with the observations in the previous two examples, because the fluid-structure interaction model takes into account the fracture closure effect. Nearly 1000 days ago, the production capacity of the well was already severely insufficient. After the implementation of the huff-and-puff measures, the production capacity was effectively improved, ultimately increasing the cumulative production by approximately 7.7%.
[0061] Although embodiments of the present invention have been shown and described above, it is understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those skilled in the art can make changes, modifications, substitutions and variations to the above embodiments within the scope of the present invention.
Claims
1. A multiphase multicomponent fluid-structure interaction numerical simulation method, characterized in that, include: S1: Construct a geometric model based on the geological data of the target research object and perform grid generation; S2: Construct the governing equations of the geometric model after mesh generation; wherein, the governing equations include the mechanical equilibrium equations characterizing the deformation behavior of the porous media rock mass skeleton and the mass transport governing equations characterizing the flow and material transport behavior of multiphase and multicomponent fluids in the porous media; S3: Based on the thermodynamic phase equilibrium theory, flash evaporation calculations are performed on each grid computing unit to determine the phase distribution and parameters of each component of the multiphase multicomponent fluid under the current thermodynamic state; S4: Construct a fully coupled solution scheme for the discretized control equations. Incorporate displacement variables representing rock mass skeleton deformation, pressure variables representing fluid flow, and component variables representing the content of each component as master variables to be solved into a unified solution framework. Simultaneously iterate and solve the master variables to be solved to obtain the numerical simulation results at the current time step.
2. The method of claim 1, wherein, The mechanical equilibrium equations in S2 are derived from Biot's porosity elasticity theory. The mechanical equilibrium equation is: ； wherein, ; is the Cauchy stress tensor; is the rock density, , = is the total fluid density, is the saturation of the phase, is the density of the phase, is the current rock porosity, is the rock density; is the gravity acceleration vector; Cauchy stress tensor The relationship between the effective stress of the porous medium rock matrix skeleton is: The relationship between the effective stress of the porous medium rock matrix skeleton is: ； ； The effective stress of the porous media rock mass skeleton was obtained by sorting out the data. for: ； in, It is a fourth-order elastic tensor; For small strain tensors, For solid displacement vectors; It is a second-order unit tensor; For Biot coefficient, The bulk elastic modulus of the rock mass skeleton. The bulk modulus of the solid particle matrix itself; Equivalent pore pressure for multiphase fluids; The formula for calculating the equivalent pore pressure of the multiphase fluid is as follows: ； In the formula: for Phase pressure of the phase.
3. The method according to claim 1, characterized in that, The mass transport control equation in S2 is constructed by using the molar concentration of each component as a conserved variable; for any component The mass transport control equation is: ； in, This represents the current porosity of the rock mass. for The molar density of the phase; Components exist The mole fraction in the phase; for Phase saturation; For the body strain; for The relative Darcy velocity vector of the phase; Components exist Hydrodynamic dispersion-diffusion tensor in the phase; For source / sink; The number of components; The The relative Darcy velocity vector of the phase is: ； in, For the absolute permeability tensor of the rock mass; for The relative permeability of the phase; for Phase viscosity; for Phase pressure of the phase; for Average molar mass of the phase, , Components The molecular weight.
4. The method according to claim 1, characterized in that, The flash evaporation calculation in S3 includes: For a fluid mixture containing multiple components, an equation of state is constructed; wherein the equation of state is: ； in, Equivalent pore pressure for multiphase fluids; This is the universal gas constant; Reservoir temperature; This represents the molar volume of the mixture. The energy parameters of the mixture; This refers to the volume parameters of the mixture; Based on the constructed equation of state, the fugacity coefficients of each component are calculated; exist Fugacity coefficient in phase satisfy: ； in, Components Volume parameters; for Phase compressibility factor for molar volume of the phase; and All are dimensionless state equation parameters; Components exist The mole fraction in the phase; Components and The interaction coefficient between them; Construct and solve the phase equilibrium equations The phase distribution and parameters of each component under the current thermodynamic state are obtained; among them, the phase equilibrium equation is... for: ； in, Components within a grid computing cell Total mole fraction; Components The phase equilibrium constant, , Components Mole fraction in the gas phase Components mole fraction in the liquid phase The number of components.
5. The method according to claim 4, characterized in that, The flash evaporation calculation also includes: Solving the phase equilibrium equation Previously, the composition was based on Wilson's empirical formula. Initialize the phase equilibrium constant: ； in, Components Critical pressure; Components The critical temperature; Components eccentricity factor; Solving the phase equilibrium equation After obtaining the phase distribution and parameters of each component under the current thermodynamic state, the fugacity coefficients are calculated, and the ratios of the calculated fugacity coefficients are used to assess the composition. The phase equilibrium constant is iteratively updated: ； in, Components Fugacity coefficient in the liquid phase; Components In the gas phase fugacity coefficient; Repeatedly solve the phase equilibrium equations And update the components The phase equilibrium constant is calculated until the phase equilibrium convergence condition of the fugacity equality criterion is met; wherein, the fugacity equality criterion is: ； in, Components Fugacity in the gas phase ; Components Fugacity in the liquid phase ; This is the preset convergence tolerance.
6. The method according to claim 1, characterized in that, Also includes: Based on the numerical simulation results at the current time step, update the porosity parameters and permeability tensor of the porous medium; The porosity parameter is based on volumetric strain. The updated result is: ； in, The initial porosity of the rock mass; For volumetric strain, For small strain tensors, For solid displacement vectors; For Biot coefficients; This represents the current porosity of the rock mass. The bulk modulus of the solid particle matrix itself; for Phase saturation; for Phase pressure of the phase; The permeability tensor It is obtained by updating the Carman-Kozeny equation based on anisotropic correction: ； in, This represents the initial permeability tensor of the rock mass. For the current volumetric strain Porosity below; It is a second-order unit tensor.
7. The method according to claim 6, characterized in that, Also includes: Based on the numerical simulation results of the current time step, the volumetric load of the rock mass skeleton is calculated according to the phase saturation distribution and pressure distribution of the multiphase fluid, and the volumetric load is fed back as an external force term in the force balance equation to the rock mass deformation solution of the next time step; the volumetric load for: 。 8. The method according to claim 6, characterized in that, In S4, the discretized control equations include: The discretized mechanical equilibrium equations are: ； in, It is the volume integral domain; For effective stress tensor; For strain tensor; For the boundary; For the boundary The known surface force tensor on; Density of the rock mass; This is a test function for rock mass displacement. , For the displacement-dependent first-order homogeneous Sobolev space; Discretized mass transport control equations: ； in, for The molar density of the phase; Components exist The mole fraction in the phase; for Phase saturation; for Average mobility of the phase; Components exist The diffusion coefficient tensor in the phase; For source / sink; Scalar test function for component molar concentration. , Let be a first-order homogeneous Sobolev space related to the mole fraction.
9. The method according to claim 8, characterized in that, In S4, the fully coupled solution scheme uses the fully implicit Newton method for nonlinear iterative solution: The residual vector after spatial discretization is constructed as follows: ； in, It is a vector transpose operator; The residual vector; The residual vector of the mechanical equilibrium equations: ； For the first The residual vectors of the mass conservation equations for each component: ； The iterative solution involves: in each Newton iteration step, constructing the Jacobian matrix of the discretized control equation with respect to the main variables to be solved; and then updating the main variables to be solved synchronously based on the Jacobian matrix until the convergence criterion is met. The Jacobian matrix includes several cross-partial derivative terms between rock mass skeleton deformation, fluid flow, and component content variations: ； ； in, A vector composed of the molar numbers of different components; The number of components.
10. A multiphase, multicomponent fluid-structure interaction numerical simulation system, said system being used to perform the steps of the method according to any one of claims 1-9, characterized in that, include: Geometric model and mesh generation module: used to construct a geometric model based on the geological data of the target research object and to perform mesh generation; Equation Construction Module: Used to construct the governing equations of the geometric model after mesh generation; wherein, the governing equations include mechanical equilibrium equations characterizing the deformation behavior of the porous media rock mass skeleton, and mass transport control equations characterizing the flow and material transport behavior of multiphase and multicomponent fluids in porous media; Phase equilibrium calculation module: Based on thermodynamic phase equilibrium theory, it performs flash evaporation calculations on each grid calculation unit to determine the phase distribution and parameters of each component of the multiphase multicomponent fluid under the current thermodynamic state; Coupled solution module: Constructs a fully coupled solution format for the discretized control equations, incorporating displacement variables characterizing rock mass skeleton deformation, pressure variables characterizing fluid flow, and component variables characterizing the content of each component as master variables to be solved into a unified solution framework, and synchronously iterating to solve the master variables to be solved to obtain the numerical simulation results of the current time step.