An unconventional oil and gas reservoir thermal fluid-structure coupling numerical simulation method and system

Abstract

The application provides an unconventional oil and gas reservoir thermal fluid-solid coupling numerical simulation method and system, and the method comprises the following steps: constructing a multiphase numerical model coupling a temperature field, a seepage field and a stress field; the multiphase numerical model integrates a Cahn-Hilliard equation, a heat conduction equation, a continuity equation, a Navier-Stokes equation and a mechanical equilibrium equation; determining a model calculation domain of the multiphase numerical model; performing grid discretization processing on the model calculation domain and performing grid refinement processing on the periphery of a mass flow inlet and a pressure outlet, thereby obtaining a plurality of grid units; inputting environmental parameters and initial condition parameters into the multiphase numerical model, controlling the multiphase numerical model to perform discrete calculation in each grid unit, judging whether a convergence condition is met, if the convergence condition is met, updating physical property parameters of a multiphase fluid according to a current temperature and pressure, and outputting a numerical simulation result of the unconventional oil and gas reservoir. The application can reduce the production capacity prediction deviation.

Description

Technical Field

[0001] This application relates to the field of unconventional oil and gas numerical simulation technology, and more specifically, to a method and system for thermal-fluid-structure interaction numerical simulation of unconventional oil and gas reservoirs. Background Technology

[0002] Unconventional oil and gas reservoirs are characterized by nanoscale pores (50-800 nm), strong heterogeneity, and multi-scale flow mechanisms, but their development faces challenges such as high temperatures. High pressure, such as High ground stress, such as The three plateau geological environments are characterized by strong thermal-fluid-solid multi-field coupling effects in fluid transport.

[0003] Existing methods primarily employ three numerical simulation approaches to guide the development of unconventional oil and gas reservoirs: the pore network model method, the embedded discrete model method, and the meshless method. However, these approaches all suffer from the following drawbacks: First, multi-field coupling is insufficient, often focusing on fluid-solid or thermal-fluid coupling and neglecting the influence of thermal expansion and high pressure on the PVT properties of fluids. Second, in-situ geological conditions are not adequately reproduced, failing to constrain parameters such as stress direction, bedding dip angle, and in-situ water saturation. Third, reservoir structure adaptability is poor, failing to achieve cross-scale flow connectivity and ignoring porosity loss caused by clay mineral expansion. These shortcomings make it difficult for current numerical simulation methods to accurately guide the development of unconventional oil and gas reservoirs, resulting in significant deviations in production capacity predictions. Summary of the Invention

[0004] In view of this, the purpose of this application is to provide a numerical simulation method and system for thermal-fluid-structure interaction of unconventional oil and gas reservoirs, which can solve the problem that current numerical simulation methods are difficult to accurately guide the development of unconventional oil and gas reservoirs and reduce the deviation of production capacity prediction.

[0005] In a first aspect, embodiments of this application provide a numerical simulation method for thermal-fluid-structure interaction in unconventional oil and gas reservoirs, the method comprising:

[0006] A multiphase numerical model is constructed that couples the temperature field, seepage field, and stress field; the multiphase numerical model integrates the Cahn-Hilliard equation, the heat conduction equation, the continuity equation, the Navier-Stokes equation, and the mechanical equilibrium equation;

[0007] The computational domain of the multiphase numerical model is determined. The computational domain includes multiple matrix particles and pores and throats between adjacent matrix particles. A mass flow inlet is provided at the first end of the computational domain, and a pressure outlet is provided at the second end of the computational domain. The pores and throats are used to store and transport fluid. The first end and the second end are located diagonally opposite to each other in the computational domain. All boundaries of the computational domain are non-flow boundaries and adiabatic boundaries.

[0008] The computational domain of the model is discretized into a grid, and the area around the mass flow inlet and the pressure outlet is refined into a grid to obtain multiple grid cells.

[0009] The environmental parameters and initial condition parameters are obtained. The environmental parameters include the physical properties of the rock matrix, the physical properties of the multiphase fluid, the initial saturation of the multiphase fluid, and the interfacial tension coefficient between different phase fluids. The initial condition parameters include the initial reservoir temperature and the initial reservoir pressure.

[0010] The environmental parameters and the initial condition parameters are input into the multiphase numerical model. The multiphase numerical model is controlled to perform discrete calculations in each grid cell and to determine whether the convergence condition is met. If the convergence condition is met, the physical property parameters of the multiphase fluid are updated according to the current temperature and pressure, and the numerical simulation results of the unconventional oil and gas reservoir are output. The numerical simulation results include phase saturation, temperature, pressure and flow rate.

[0011] Secondly, embodiments of this application also provide a numerical simulation system for unconventional oil and gas reservoir thermal-fluid-structure interaction, comprising:

[0012] The model building module is used to construct a multiphase numerical model that couples the temperature field, seepage field, and stress field; the multiphase numerical model integrates the Cahn-Hilliard equation, the heat conduction equation, the continuity equation, the Navier-Stokes equation, and the mechanical equilibrium equation.

[0013] The computational domain determination module is used to determine the computational domain of the multiphase numerical model. The computational domain includes multiple matrix particles and pores and throats between adjacent matrix particles. The first end of the computational domain is provided with a mass flow inlet, and the second end of the computational domain is provided with a pressure outlet. The pores and throats are used to store and transport fluid. The first end and the second end are located diagonally opposite to each other in the computational domain, and all boundaries of the computational domain are non-flow boundaries and adiabatic boundaries.

[0014] The mesh generation module is used to discretize the model computation domain and refine the mesh around the mass flow inlet and the pressure outlet to obtain multiple mesh cells.

[0015] The parameter acquisition module is used to acquire environmental parameters and initial condition parameters. The environmental parameters include the physical properties of the rock matrix, the physical properties of the multiphase fluid, the initial saturation of the multiphase fluid, and the interfacial tension coefficient between different phase fluids. The initial condition parameters include the initial reservoir temperature and the initial reservoir pressure.

[0016] The numerical simulation module is used to input the environmental parameters and the initial condition parameters into the multiphase numerical model, control the multiphase numerical model to perform discrete calculations in each grid cell, and determine whether the convergence condition is met. If the convergence condition is met, the physical property parameters of the multiphase fluid are updated according to the current temperature and pressure, and the numerical simulation results of the unconventional oil and gas reservoir are output. The numerical simulation results include phase saturation, temperature, pressure and flow rate.

[0017] Thirdly, embodiments of this application also provide an electronic device, including: a processor, a memory, and a bus. The memory stores machine-readable instructions executable by the processor. When the electronic device is running, the processor communicates with the memory via the bus. When the machine-readable instructions are executed by the processor, the steps of the unconventional oil and gas reservoir thermal-fluid-structure interaction numerical simulation method described above are performed.

[0018] Fourthly, embodiments of this application also provide a computer-readable storage medium storing a computer program, which, when executed by a processor, performs the steps of the unconventional oil and gas reservoir thermal-fluid-structure interaction numerical simulation method described above.

[0019] The unconventional oil and gas reservoir thermal-fluid-solid coupling numerical simulation method and system provided in this application, by constructing a multiphase numerical model that couples temperature field, seepage field, and stress field and integrates multiple key equations, can fully realize full dynamic coupling of thermal-fluid-solid processes, avoiding the shortcomings of insufficient multi-field coupling in existing technologies. At the same time, by combining the precise design of the model's computational domain, the differentiated discretization of the grid, and the comprehensive acquisition of environmental parameters and initial condition parameters, it can effectively restore in-situ geological conditions, adapt to the strong heterogeneity of unconventional reservoirs, and thus accurately simulate the micro-transport process of multiphase fluids, significantly improving the accuracy of production capacity prediction, reducing prediction bias, and accurately guiding the development of unconventional oil and gas reservoirs.

[0020] To make the above-mentioned objectives, features and advantages of this application more apparent and understandable, preferred embodiments are described below in detail with reference to the accompanying drawings. Attached Figure Description

[0021] To more clearly illustrate the technical solutions of the embodiments of this application, the accompanying drawings used in the embodiments will be briefly introduced below. It should be understood that the following drawings only show some embodiments of this application and should not be regarded 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.

[0022] Figure 1 A flowchart illustrating a numerical simulation method for unconventional oil and gas reservoir thermal-fluid-structure interaction provided in this application embodiment;

[0023] Figure 2 This is a schematic diagram of a model computation domain meshing provided in an embodiment of this application;

[0024] Figure 3 This is a schematic diagram of another model computation domain meshing provided in an embodiment of this application;

[0025] Figure 4 A schematic diagram illustrating the change in gas phase saturation versus calculation time, provided as an embodiment of this application;

[0026] Figure 5 A gas phase saturation distribution diagram at different times provided in an embodiment of this application;

[0027] Figure 6 A flow velocity distribution diagram at different times provided in an embodiment of this application;

[0028] Figure 7 An oil phase saturation distribution diagram at different times provided in an embodiment of this application;

[0029] Figure 8 A water phase saturation distribution diagram at different times provided in an embodiment of this application;

[0030] Figure 9(a) is a schematic diagram of reference points and reference lines in a model calculation domain provided in an embodiment of this application;

[0031] Figure 9(b) is a diagram showing the saturation variation of each phase at a reference point provided in an embodiment of this application;

[0032] Figure 9(c) is a graph showing the change of oil phase saturation at different reference lines at different times, provided in an embodiment of this application.

[0033] Figure 10(a) is a flow rate variation diagram provided in an embodiment of this application;

[0034] Figure 10(b) is a diagram showing the variation of a seepage region provided in an embodiment of this application;

[0035] Figure 11(a) is a diagram showing the gas phase saturation variation considering and not considering the temperature field-stress field according to an embodiment of this application.

[0036] Figure 11(b) is a diagram showing the change of gas phase area considering and not considering the temperature field-stress field according to an embodiment of this application.

[0037] Figure 12 A diagram showing the variation of gas phase saturation under different mass flow rates is provided in an embodiment of this application.

[0038] Figure 13This is a schematic diagram illustrating the change of gas phase region area with mass flow rate, provided in an embodiment of this application.

[0039] Figure 14 This is a schematic diagram illustrating the variation of gas phase saturation with particle diameter provided in an embodiment of this application.

[0040] Figure 15(a) is a schematic diagram of the change of gas phase area with particle diameter provided in an embodiment of this application;

[0041] Figure 15(b) is a schematic diagram showing the change of gas phase area ratio with particle diameter according to an embodiment of this application;

[0042] Figure 16 This is a schematic diagram of a horizontally layered reservoir computational domain and a vertically layered reservoir computational domain provided in an embodiment of this application;

[0043] Figure 17 A diagram showing the variation of gas phase saturation in a vertically layered structure in a layer with different particle diameters, provided as an embodiment of this application.

[0044] Figure 18 This application provides a fluid velocity distribution diagram in layers with different particle diameters at 60s, as shown in the embodiments of the present application.

[0045] Figure 19(a) is a graph showing the change in the percentage of a region in a vertical layer provided in an embodiment of this application;

[0046] Figure 19(b) is a graph showing the change in the percentage of a horizontally layered region provided in an embodiment of this application;

[0047] Figure 20 A graph showing the variation of gas phase saturation in a layer with different particle diameters under horizontal lamination, provided in an embodiment of this application.

[0048] Figure 21 A distribution diagram of fluid velocity in a layer of different particle diameters at 60 s in a horizontal layered structure, provided as an embodiment of this application;

[0049] Figure 22 This is a schematic diagram of the structure of an unconventional oil and gas reservoir thermal-fluid-structure interaction numerical simulation system provided in an embodiment of this application;

[0050] Figure 23 This is a schematic diagram of the structure of an electronic device provided in an embodiment of this application. Detailed Implementation

[0051] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. The components of the embodiments of this application described and shown in the accompanying drawings can generally be arranged and designed in various different configurations. Therefore, the following detailed description of the embodiments of this application provided in the accompanying drawings is not intended to limit the scope of the claimed application, but merely represents selected embodiments of this application. Based on the embodiments of this application, every other embodiment obtained by those skilled in the art without inventive effort falls within the scope of protection of this application.

[0052] Compared with conventional oil and gas reservoirs, unconventional oil and gas reservoirs have significant differences: First, the pore scale is mainly at the nanoscale, with the pore size concentrated in the range of 50-800 nm, and the pore structure is complex and poorly connected; second, the reservoir is highly heterogeneous, with extensive development of bedding and fractures, resulting in significant differences in the permeability distribution inside the reservoir; third, the fluid flow mechanism exhibits multi-scale coexistence characteristics, including diffusion flow in the nanomatrix, slip flow in microfractures, and Darcy flow in hydraulic fractures, with complex fluid transport patterns.

[0053] More importantly, the development of unconventional oil and gas reservoirs faces the challenge of high temperatures. High pressure, such as High ground stress, such as The high-altitude geological environment of the three plateaus results in strong thermo-fluid-solid multi-field coupling effects in fluid transport within the reservoir: changes in the temperature field directly alter fluid viscosity and phase state, such as crude oil cracking and a decrease in heavy oil viscosity; changes in effective stress induce significant deformation of nanopores, such as a pore compressibility as high as [missing information]. This, in turn, alters pore connectivity; the dynamic changes in the seepage field continuously adjust the pore pressure distribution, in turn affecting the effective stress and temperature transmission.

[0054] Currently, the core challenges facing the development of unconventional oil and gas reservoirs lie in the unclear microscale multiphase flow mechanisms and insufficient in-situ geological condition reconstruction, directly leading to production prediction deviations exceeding 30%. Specifically, the multiphase competitive transport mechanism is complex. Within the nanopore throat, oil, gas, and water competitively adsorb / desorb, with capillary forces reaching MPa levels, 10-100 times that of conventional reservoirs, significantly inhibiting the saturation of mobile fluids. The thermal-fluid-solid dynamic coupling effect is significant; during extraction, reservoir temperature changes trigger rock thermal expansion and stress redistribution, thereby altering pore connectivity, with permeability losses reaching up to 70%. The transformation of multiscale flow patterns is difficult to characterize. Fluid flow transitions from diffusion flow in the nanomatrix to slip flow in microfractures to Darcy flow in hydraulic fractures, spanning a scale difference of 10 orders of magnitude. Traditional single-scale models cannot adapt to the cross-scale changes from diffusion flow to Darcy flow. Therefore, there is an urgent need to establish a numerical simulation method capable of reconstructing in-situ geological conditions, coupling multiphysics fields, and connecting multiple scales to provide guidance for the efficient development of unconventional oil and gas reservoirs.

[0055] For the numerical simulation needs of unconventional oil and gas reservoirs, existing technologies mainly focus on three categories: pore network model methods, embedded discrete model methods, and meshless methods. The specific implementation methods are as follows:

[0056] Firstly, the porous network model method:

[0057] This method is based on actual observation data of core pore structure and constructs a pore network model using digital means to simulate the fluid transport process in the pores. For example, the quadrangular star-shaped cross-section model uses geometric parameterization to describe the throat cross-section and simulates two-phase flow through a piston / clamping filling mechanism. However, it only considers the quasi-static process dominated by capillary forces, does not couple the pore deformation caused by stress fields, such as the dynamic changes in pore throat radius, and cannot handle phase change behavior under thermal recovery conditions.

[0058] For example, the nuclear magnetic resonance reconstruction model uses nuclear magnetic resonance imaging (MRI) or T2 spectrum scanning data to obtain the three-dimensional distribution information of pores inside the core. By using the conversion coefficient α to correlate the imaging data with the measured permeability, the accurate reconstruction of the three-dimensional pore network is achieved. Although it realizes the core-scale flow simulation, it does not introduce temperature field parameters and cannot simulate key processes such as viscosity reduction of heavy oil in heat injection development.

[0059] Secondly, the embedded discrete model method:

[0060] This method targets unconventional reservoirs with fractured or vulnerable formations, employing unstructured mesh technology to couple flow domains at different scales, achieving co-simulation of matrix and fracture (or cavern) flow. For example, in the fractured-vulnerable fluid-structure interaction model, embedded discrete fracture facies (EDFM) technology couples the free flow within the cavern with the seepage field within the matrix. In model construction, free flow equations for the fluid within the cavern and seepage equations for the fluid within the matrix are established separately. However, this method does not consider the microscale effects of nanopores, such as boundary slip and diffusion transport, leading to distortion in the prediction of seepage in fine-grained reservoir matrix.

[0061] For example, the THMC fully coupled simulator focuses on the thermal injection exploitation scenario of hydrate reservoirs and realizes the simulation of temperature-pressure-phase equilibrium evolution in hydrate thermal injection exploitation. However, this method does not distinguish the wettability differences between the matrix and fractures, and does not introduce in-situ stress constraints, such as controlling the fracture aperture by the direction of the maximum horizontal principal stress.

[0062] Thirdly, the meshless method:

[0063] This method addresses the problem of grid distortion in highly heterogeneous reservoirs by employing meshless numerical methods such as the generalized finite difference method (GFDM) to avoid the insufficient adaptability of traditional grid generation to heterogeneous reservoirs.

[0064] For example, the gas-water two-phase fluid-structure coupling model divides shale reservoirs into primary modified zones, secondary modified zones, and unmodified zones according to the degree of modification. It constructs a unified transport model with multiple flow regimes based on the flow characteristics of different regions, such as Darcy flow in modified zones and diffusion flow in unmodified zones. It uses the generalized finite difference method for numerical solution to avoid the impact of grid distortion on calculation accuracy. However, this method does not integrate the temperature field and has insufficient sensitivity analysis to water saturation.

[0065] However, a comprehensive analysis of the three existing technical solutions reveals that their core shortcomings lie in three aspects: insufficient multi-field coupling, inadequate in-situ geological condition reconstruction, and poor adaptability to unconventional reservoir structures. These are detailed below:

[0066] First, the multi-field coupling is insufficient and cannot reflect the full dynamic interaction between heat, fluid, and solid structures:

[0067] Existing solutions mostly focus on fluid-structure interaction, such as the pore network model method, the cavity-type fluid-structure interaction model, or heat-fluid interaction, such as the heat-fluid-chemical coupling of the THMC fully coupled simulator. They generally lack a fully dynamic coupling mechanism among heat, fluid, and solid, and cannot fully characterize the mutual feedback between temperature, seepage, and stress, such as the influence of rock thermal expansion caused by temperature changes on the seepage field, and the reverse effect of pore deformation caused by stress changes on temperature transfer.

[0068] Ignoring the thermal expansion effect of rocks caused by temperature changes, the thermal stress term is missing from the mechanical equilibrium equations, leading to deviations in the calculation of effective stress and failing to accurately reflect the influence of temperature on reservoir skeleton deformation and pore structure. Furthermore, the dynamic effects of in-situ high-pressure conditions on fluid PVT properties, such as supercritical fluid dynamics, are not considered. Sudden density changes and the pressure sensitivity of crude oil viscosity cause fluid properties to become disconnected from actual reservoir conditions, further amplifying simulation errors.

[0069] Secondly, the in-situ geological conditions are not sufficiently reconstructed, resulting in a disconnect between the model and the actual reservoir:

[0070] Existing technologies do not impose mandatory constraints on key in-situ geological parameters during model construction, resulting in models that cannot accurately reproduce the actual geological environment of the reservoir. Specifically, this manifests as follows:

[0071] The simulation fails to incorporate the directional characteristics of in-situ stress, such as the maximum horizontal principal stress σH being greater than the minimum horizontal principal stress σh. This makes it impossible to simulate the control of stress direction on fracture propagation azimuth and pore throat deformation direction, which is inconsistent with the actual development of fractures along the direction of the maximum horizontal principal stress in reservoirs. Furthermore, the simulation does not consider the dominant role of reservoir bedding dip angle (30°-60°) on the anisotropic permeability tensor, failing to characterize the permeability difference between the bedding direction and the direction perpendicular to the bedding. Additionally, the simulation does not use in-situ water saturation, such as in-situ water saturation Sw being greater than 40%, as a constraint condition, neglecting its inhibitory effect on effective gas phase permeability, leading to overly optimistic gas phase migration simulation results.

[0072] Third, unconventional reservoir structures have poor adaptability and insufficient multi-scale and heterogeneous treatment:

[0073] Existing methods fail to achieve cross-scale flow continuity. They either focus on the nanomatrix and microfractures (e.g., pore network models) or on microfractures and hydraulic fractures (e.g., embedded discrete models). These methods cannot fully characterize the cross-scale transport path of fluid from the nanomatrix to the microfractures and then to the hydraulic fractures, leading to an overestimation of the fluid supply from the matrix to the fractures and significant biases in production capacity predictions. Furthermore, they do not incorporate the swelling effect of clay minerals such as montmorillonite and illite. Unconventional fine-grained reservoirs often have high clay mineral content, and their swelling upon contact with water can cause porosity losses of up to 30%, directly altering pore structure and permeability. Existing models ignore this effect, resulting in simulations of dynamic changes in porosity and permeability that are significantly inconsistent with actual reservoirs, further reducing the reliability of the simulation results.

[0074] The aforementioned shortcomings make it difficult for current numerical simulation technology to accurately guide the design of development schemes for unconventional oil and gas reservoirs. There is an urgent need to develop a fully coupled, in-situ constrained, multi-scale adapted numerical simulation method to address the deficiencies of existing technologies.

[0075] Based on this, the embodiments of this application provide a numerical simulation method for thermal-fluid-structure interaction of unconventional oil and gas reservoirs, which can solve the problem that current numerical simulation methods are difficult to accurately guide the development of unconventional oil and gas reservoirs and reduce the deviation of production capacity prediction.

[0076] Please see Figure 1 , Figure 1 This is a flowchart illustrating a numerical simulation method for unconventional oil and gas reservoir thermal-fluid-structure interaction provided in an embodiment of this application. Figure 1 As shown in the embodiments of this application, the method includes:

[0077] S101. Construct a multiphase numerical model that couples the temperature field, seepage field, and stress field; the multiphase numerical model integrates the Cahn-Hilliard equation, the heat conduction equation, the continuity equation, the Navier-Stokes equation, and the mechanical equilibrium equation.

[0078] S102. Determine the computational domain of the multiphase numerical model. The computational domain includes multiple matrix particles and pores and throats between adjacent matrix particles. The first end of the computational domain is provided with a mass flow inlet, and the second end of the computational domain is provided with a pressure outlet. The pores and throats are used to store and transport fluid. The first end and the second end are located at opposite corners of the computational domain, and all boundaries of the computational domain are non-flow boundaries and adiabatic boundaries.

[0079] S103. The model computation domain is discretized into a grid, and the area around the mass flow inlet and pressure outlet is refined into a grid to obtain multiple grid cells.

[0080] S104. Obtain environmental parameters and initial condition parameters. Environmental parameters include the physical properties of the rock matrix, the physical properties of the multiphase fluid, the initial saturation of the multiphase fluid, and the interfacial tension coefficient between different phase fluids. Initial condition parameters include the initial reservoir temperature and the initial reservoir pressure.

[0081] S105. Input the environmental parameters and initial condition parameters into the multiphase numerical model, control the multiphase numerical model to perform discrete calculations in each grid cell, and determine whether the convergence condition is met. If the convergence condition is met, update the physical property parameters of the multiphase fluid according to the current temperature and pressure, and output the numerical simulation results of the unconventional oil and gas reservoir. The numerical simulation results include phase saturation, temperature, pressure and flow rate.

[0082] In the embodiments of the present application, by constructing a multiphase numerical model that couples the temperature field, seepage field, and stress field and integrates various key equations, the full dynamic coupling of heat-fluid-solid can be fully realized, avoiding the defect of insufficient multi-field coupling in the prior art; at the same time, combined with the precise design of the model calculation domain, the differential discretization of the grid, and the comprehensive acquisition of environmental parameters and initial condition parameters, the in-situ geological conditions can be effectively restored, adapting to the strong heterogeneity of unconventional reservoirs, and then accurately simulating the micro-migration process of multiphase fluids, significantly improving the accuracy of production capacity prediction, reducing the prediction deviation, and being able to accurately guide the development of unconventional oil and gas reservoirs.

[0083] The above steps are exemplarily illustrated through specific embodiments as follows:

[0084] In step S101: Construct a multiphase numerical model that couples the temperature field, seepage field, and stress field; the multiphase numerical model integrates the Cahn-Hilliard equation, heat conduction equation, continuity equation, Navier-Stokes equation, and mechanical equilibrium equation.

[0085] Among them, by coupling the temperature field, seepage field, and stress field, the interaction relationship between each field is established, so that the temperature field, seepage field, and stress field act synergistically to jointly reflect the physical processes in the reservoir. The change in the temperature field will affect the physical properties of the fluid and then change the state of the seepage field. The change in pore pressure in the seepage field will affect the effective stress distribution in the stress field, and the deformation of the rock matrix in the stress field will change the pore structure and thus affect the seepage field, forming a dynamic coupling relationship among the three.

[0086] In an optional implementation manner, coupling the temperature field, seepage field, and stress field is based on the actual physical mechanism of the reservoir, incorporating the influence of temperature on fluid viscosity and density, the influence of pore pressure on the effective stress of the rock, and the influence of rock deformation on pore throat size into the same model system to achieve the coupled calculation of multiple fields. For example, when the reservoir temperature rises, the viscosity of the aqueous phase will decrease, reducing the fluid flow resistance and increasing the flow velocity in the seepage field; while the increase in flow velocity will accelerate the heat exchange between the fluid and the rock matrix, affecting the distribution of the temperature field; at the same time, the increase in pore pressure caused by fluid injection will reduce the effective stress of the rock, causing the matrix particles to undergo small deformations and the pore throat size to increase, further changing the seepage characteristics and realizing the dynamic coupling of the three fields.

[0087] Here, the multiphase numerical model refers to a mathematical model used to describe and calculate the migration law of oil, gas, and water three-phase fluids in the reservoir under the coupled action of the temperature field, seepage field, and stress field, and is the core carrier integrating various physical equations and parameters.

[0088] In an optional implementation, the multiphase numerical model is a mathematical model system built based on the finite element method, which can integrate multiphysics equations to achieve quantitative calculation of multiphase fluid transport processes. For example, based on the physical properties of the Qingshankou Formation reservoir in the Songliao Basin, the model can set the structural parameters of blocky or layered reservoirs, and calculate the phase saturation, temperature, pressure, and flow velocity at various points in the reservoir at different times through integrated multiple equations, providing a quantitative tool for analyzing reservoir fluid transport patterns.

[0089] The Cahn-Hilliard equations are mathematical equations used to describe the evolution and phase distribution of immiscible multiphase fluid interfaces. They characterize the proportion of each phase in space through phase field variables, thereby simulating the migration process of multiphase interfaces.

[0090] Specifically, the Cahn-Hilliard equations, by introducing phase field variables and based on the principle of minimizing system free energy, simulate the smooth migration of multiphase fluid interfaces. For example, the Cahn-Hilliard equations are expressed by the following formula:

[0091]

[0092] in, Let t represent the phase field variable and t represent time. M0 represents the fluid velocity, and M0 represents the fluid flow regulation parameter. This represents capillary force, where i represents oil phase A, water phase B, and gas phase C, respectively. This represents the phase field variable of oil phase A. Represents the phase field variables of water phase B. Represents the phase field variables of gas phase C;

[0093] This can be expressed by the following formula:

[0094]

[0095] in, Denotes the coefficients of the Cahn–Hilliard equation. Parameters that control the thickness of the interface; This represents the partial derivative of the free energy with respect to the order parameter. To express summation;

[0096] ΣA, ΣB, and ΣC are represented by the following formulas:

[0097]

[0098] This represents the interfacial tension coefficient between the oil phase and the water phase. This represents the interfacial tension coefficient between the oil phase and the gas phase. This represents the interfacial tension coefficient between the aqueous phase and the gas phase.

[0099] F(ψ) represents the volume energy, expressed by the following formula:

[0100]

[0101] The additional free volume energy is represented by Λ;

[0102] Furthermore, when simulating the blocky reservoirs of the Songliao Basin, the Cahn-Hilliard equation can calculate the finger-like propagation process of the gas phase in the pore throat using the above formula. The gas phase field variable gradually expands from 0.1 at the inlet (i.e., an initial saturation of 10%) towards the outlet, intuitively reflecting the displacement range of the gas phase and avoiding the computational distortion of traditional sharp interface models at the microscale. The phase field variable, which characterizes the proportion of a certain phase (oil, water, or gas) at a specific point in space, ranges from 0 to 1 and is the core parameter of the Cahn-Hilliard equation describing phase distribution.

[0103] Optionally, the value of the phase field variable directly corresponds to the saturation of a certain phase. For example, in the simulation of layered reservoirs, the initial value of the oil phase field variable in the fine-grained layer is 0.65, which corresponds to an initial oil phase saturation of 65%. As the gas phase is injected, the gas phase field variable gradually increases, while the oil phase field variable gradually decreases. The distribution dynamics of the oil and gas phases can be tracked in real time through the change of the phase field variable.

[0104] The heat conduction equation is a mathematical equation used to describe the heat transfer process within a reservoir. It calculates the spatial and temporal distribution of temperature and reflects the evolution of the thermal field. For example, the heat conduction equation is expressed by the following formula:

[0105]

[0106] in, and These represent the effective volumetric heat capacity and effective thermal conductivity, respectively.

[0107]

[0108]

[0109] in, These represent the density of the solid and the density of the fluid, respectively. and These represent the specific heat capacity of a solid and the specific heat capacity of a fluid, respectively. and Let Q represent the thermal conductivity of the solid and the thermal conductivity of the fluid, respectively, and let Q represent the energy source. Indicates porosity;

[0110] The thermal conductivity of a fluid is expressed by the following formula:

[0111]

[0112] in, Indicates the density of the oil phase. This indicates the density of the aqueous phase. Indicates gas phase density, Indicates the thermal conductivity of the oil phase. Indicates the thermal conductivity of the aqueous phase. Indicates the thermal conductivity of the gas phase;

[0113] The specific heat capacity of a fluid is expressed by the following formula:

[0114]

[0115] in, This indicates the relative heat capacity of oil. This indicates the specific heat capacity of water. This indicates the relative heat capacity of the gas.

[0116] Here, the heat conduction equation calculates the reservoir temperature field distribution by considering the heat conduction characteristics of the rock matrix and fluids, as well as the heat exchange brought about by fluid convection. For example, in simulating a thermal injection development scenario, this equation can calculate the heat diffusion process of the injected hot fluid in the reservoir using the above formula. The effective volumetric heat capacity and effective thermal conductivity, combined with the influence of fluid flow velocity on thermal convection, can yield the temperature gradient change from the inlet to the outlet in the reservoir at different times, so as to analyze the influence of temperature on fluid properties.

[0117] The continuity equation, based on the law of conservation of mass, describes the constant mass of a fluid during flow and forms the basis for calculating fluid seepage velocity and mass conservation. The Navier-Stokes equations, on the other hand, describe the conservation of fluid momentum and can be used to calculate velocity distribution and flow resistance; they are the core equations for analyzing fluid motion.

[0118] In an alternative implementation, the continuity equation calculates the changes in pore pressure and fluid flow rate within the reservoir by controlling the balance between the mass inflow and outflow of the fluid. The Navier-Stokes equations calculate the fluid flow velocity within the pore throat by considering the fluid's viscous forces, inertial forces, and volume forces (such as interfacial tension). For example, in simulating a mass flow rate inlet (injection rate) Under these conditions, this equation can be combined with the Navier-Stokes equation to ensure that the mass of fluid flowing into a certain grid cell is equal to the sum of the outflow mass and the change in fluid mass within the cell, thus avoiding calculation results that do not conserve mass and providing a guarantee for accurate calculation of flow velocity and pressure distribution.

[0119] For example, the Navier-Stokes equations and the continuity equation are expressed by the following formulas:

[0120]

[0121] in, Let I represent the fluid density, I represent the identity matrix, and F represent the volume force vector. Represents volume force;

[0122] K represents the viscous stress tensor, expressed by the following formula:

[0123]

[0124] The density and viscosity of a fluid mixture are expressed by the following formula:

[0125]

[0126] The viscosity of a fluid mixture is expressed by the following formula:

[0127]

[0128] in, Indicates the viscosity of the oil phase. Indicates the viscosity of the aqueous phase. Indicates the viscosity of the gas phase;

[0129] Volume forces are expressed by the following formula:

[0130]

[0131] Here, when simulating the coarse-grained layer of a layered reservoir, this equation can be used to calculate the high-speed flow state of the fluid within the large pore throat using the above formula. The viscous stress tensor and interfacial tension volume forces are also calculated using the above formula, ultimately yielding the maximum flow velocity within the coarse-grained layer. The results (at 50 s) show a significant difference in flow velocity compared to the fine-grained layer.

[0132] Among them, the mechanical equilibrium equation refers to the mathematical equation used to describe the mechanical equilibrium of the rock matrix under stress. It can calculate the displacement and deformation of the rock and reflect the influence of the stress field on the reservoir structure.

[0133] In an alternative implementation, the mechanical equilibrium equations calculate the deformation of the rock matrix by taking into account the rock's elastic properties, such as Young's modulus, Poisson's ratio, and the effects of thermal expansion and pore pressure. For example, the mechanical equilibrium equations are expressed by the following formula:

[0134]

[0135] Where E represents Young's modulus and v represents Poisson's ratio. This represents the displacement of the boundary under ground stress in the initial state. This indicates the boundary displacement caused by production. This represents the Biot–Willi coefficient. This represents the fluid pressure inside the pores. The symbol for Kronecker, T represents the coefficient of thermal expansion, and T represents temperature. Indicates the initial temperature. This represents the force per unit volume.

[0136] Specifically, when simulating the Songliao Basin reservoir, the deformation of matrix particles under an initial reservoir pressure of 30 MPa can be calculated using the above formula. When the reservoir temperature increases by 5°C from 90°C, the matrix will undergo slight displacement due to thermal expansion, which will increase the pore throat size and thus improve the reservoir permeability, reflecting the influence of the stress field on the seepage field.

[0137] In step S102, the computational domain of the multiphase numerical model is determined. The computational domain includes multiple matrix particles and pores and throats between adjacent matrix particles. The first end of the computational domain is provided with a mass flow inlet, and the second end of the computational domain is provided with a pressure outlet. The pores and throats are used to store and transport fluid. The first end and the second end are located at opposite corners of the computational domain, and all boundaries of the computational domain are non-flow boundaries and adiabatic boundaries.

[0138] Here, the model computation domain refers to the physical space range in which the multiphase numerical model performs calculations. It is the spatial carrier for simulating reservoir structures including matrix particles, pore throats, and setting boundary conditions, and determines the scope and object of the calculation.

[0139] In an optional implementation, the model computation domain is a three-dimensional or two-dimensional spatial region defined according to the scale of the laboratory experiment and the actual characteristics of the reservoir, used to define the physical scope of the numerical calculation. For example... Figure 2 As shown, for unconventional fine-grained oil and gas reservoirs in the Songliao Basin, the computational domain is set as a square two-dimensional region with a side length of 1000 mm. This region contains matrix particles 201 with a diameter of 40 mm and a center-to-center distance of 50 mm between adjacent matrix particles 201. Pore throats 202 are formed between matrix particles 201, which accurately restores the structural characteristics of the reservoir after deformation near the wellbore and provides a clear spatial range for subsequent grid discretization and equation solving.

[0140] Matrix particles refer to the solid particles that constitute the reservoir framework and form the basis for pore-throat spaces. Their size and distribution directly determine the pore structure and permeability of the reservoir. Here, matrix particles are solid units that simulate the reservoir rock framework, and their diameter and arrangement can be determined according to the reservoir type, such as massive or layered. For example, in the computational domain of massive reservoirs, the diameter of matrix particles can be uniformly set to 40 mm, and the center-to-center distance between adjacent matrix particles is 50 mm, forming a uniform pore-throat network. In the computational domain of layered reservoirs, the particles of the coarse-grained layer (45 mm), medium-grained layer (40 mm), and fine-grained layer (30 mm) are arranged in horizontal or vertical layers to simulate the layering structure of the reservoir. Matrix particles of different sizes correspond to different pore-throat sizes, thus affecting fluid migration characteristics.

[0141] Specifically, the pore throat refers to the void space between adjacent matrix particles, serving as a channel for the storage and flow of oil, gas, and water three-phase fluids within the reservoir. Its size and connectivity directly affect the fluid permeability. Optionally, the pore throat is the void region between matrix particles, its size determined by the particle diameter and arrangement, and is the core channel for fluid storage and migration. For example, in a blocky reservoir with matrix particle diameters of 40 mm and a center-to-center distance of 50 mm between adjacent matrix particles, the narrowest point of the pore throat is approximately 10 mm. This space can store oil, gas, and water three-phase fluids, and fluids are only allowed to flow within the pore throat. After hydraulic fracturing of the reservoir, the pore throat size can expand from the nanometer scale to the millimeter scale. For example, when simulating a fracturing reservoir, the pore throat size can be set to 20 mm to significantly improve fluid flowability and increase gas-phase displacement efficiency by more than 30%.

[0142] Furthermore, the mass flow inlet is the boundary in the model's computational domain used for injecting fluid. Fluid is input into the reservoir at a fixed mass flow rate, serving as a key boundary condition for simulating the fluid injection process. The pressure outlet is the boundary in the model's computational domain used for discharging fluid. A fixed pressure value is set at the boundary pressure to create a pressure gradient for the fluid flow, forming a key boundary condition for simulating the fluid recovery process.

[0143] In an alternative implementation, the mass flow inlet is a boundary set at the boundary of the computational domain, where fluid is injected at a constant mass flow rate to simulate the fluid injection process in actual development. For example, in the Songliao Basin reservoir simulation, such as... Figure 2 As shown, the lower left corner of the computational domain is designated as the mass flow inlet 203, with a length of 80 mm and an injection mass flow rate set to 1.0 g / s. This inlet continuously injects gaseous fluid into the reservoir, simulating the gas-driven development process in the field. By controlling the injection flow rate, the influence of different injection intensities on fluid migration is studied. In another optional embodiment, the pressure outlet is set at the boundary of the computational domain, controlling the fluid outflow with a constant pressure value, to simulate the fluid production process in actual development. Figure 2As shown, the upper right corner of the computational domain is set as pressure outlet 204, the length of pressure outlet 204 is 80mm, and the outlet pressure of pressure outlet 204 is 2.0MPa lower than the reservoir pressure. For example, it can be set to 28MPa, which is 2.0MPa lower than the initial reservoir pressure of 30MPa, forming a pressure gradient from mass flow inlet 203 to pressure outlet 204, driving the fluid to flow in the pore throat 202. When the outlet pressure 204 decreases, the pressure gradient increases, the fluid velocity increases, and the gas phase sweep range expands, which can simulate the reservoir development effect under different production pressures.

[0144] Here, the first end of the model computational domain is provided with a mass flow inlet, and the second end of the model computational domain is provided with a pressure outlet. The first end and the second end are located at opposite corners of the model computational domain, forming a diagonal flow path, which maximizes the fluid transport distance in the computational domain and more comprehensively reflects the flow characteristics in the reservoir.

[0145] For example, the mass flow inlet is placed in the lower left corner of a square computational domain, and the pressure outlet is placed in the upper right corner. This forces the fluid, after being injected through the inlet, to traverse the entire computational domain diagonally before exiting through the outlet, thus extending the fluid's migration path. For instance, when simulating massive reservoirs, a diagonal flow path allows the fluid to fully contact the pore throats at different locations, more realistically reflecting the process of gas displacing the oil and water phases. This avoids insufficient displacement caused by excessively short flow paths, making the simulation results more closely resemble the actual fluid migration patterns in reservoirs.

[0146] In addition, the fact that all boundaries of the model's computational domain are non-flowable and adiabatic boundaries means that, except for the mass flow inlet and pressure outlet, all other boundaries of the computational domain are set as boundary conditions that prevent fluid flow and heat transfer, ensuring that fluid flow and heat transfer within the computational domain only occur within the set inlet-outlet path and avoiding boundary interference.

[0147] In an optional implementation, the upper boundary of the square computational domain (e.g., y=1000mm), the lower boundary (e.g., y=0, excluding the inlet section), the left boundary (e.g., x=0, excluding the inlet section), and the right boundary (e.g., x=1000mm, excluding the outlet section) are all set as non-flowing and adiabatic boundaries. This means that fluid cannot flow into or out through these boundaries, and heat cannot be transferred through them. For example, when simulating the reservoir temperature field, the adiabatic boundary ensures that the heat from the injected fluid diffuses only within the computational domain and does not dissipate to the external environment, making the temperature field calculation results more accurate. The non-flowing boundary prevents fluid from flowing out along undefined paths, ensuring that the fluid flows only along the inlet-outlet path, consistent with the closed characteristics of actual reservoirs.

[0148] In one optional embodiment, the model computation domain includes a blocky reservoir computation domain and a layered reservoir computation domain. The matrix particles in the blocky reservoir computation domain have a uniform diameter, ranging from 30mm to 45mm, and the center-to-center distance between adjacent matrix particles is 50mm. The layered reservoir computation domain includes a horizontally layered reservoir computation domain or a vertically layered reservoir computation domain. Both the horizontally layered reservoir computation domain and the vertically layered reservoir computation domain include a coarse-grained reservoir computation domain, a medium-grained reservoir computation domain, and a fine-grained reservoir computation domain. The matrix particles in the coarse-grained reservoir computation domain have a diameter of 45mm, the matrix particles in the medium-grained reservoir computation domain have a diameter of 40mm, and the matrix particles in the fine-grained reservoir computation domain have a diameter of 30mm.

[0149] Here, the model computational domain is divided into two categories: one is the blocky reservoir computational domain, where the matrix particle diameter is uniform and ranges from 30mm to 45mm, and the center-to-center distance between adjacent matrix particles is 50mm; the other is the layered reservoir computational domain, which includes horizontally layered reservoir computational domains or vertically layered reservoir computational domains, including coarse-grained reservoir computational domains with matrix particle diameters of 45mm, medium-grained reservoir computational domains with matrix particle diameters of 40mm, and fine-grained reservoir computational domains with matrix particle diameters of 30mm.

[0150] The above settings can specifically reproduce the real structure of blocky and layered reservoirs, avoiding the distortion of traditional single-model simulations, and ensuring that the fluid migration patterns of different types of reservoirs can be accurately characterized. The matrix particle diameter is directly related to the pore throat size and permeability. For example, coarse-grained layers have large pore throats and high permeability. The layered design can clearly simulate the differences in fluid migration in different particle layers. For example, coarse-grained layers are prone to forming high-speed channels. Unifying the particle spacing and clarifying the layered particle size allows the model parameters to correspond to the actual reservoir properties, reducing parameter ambiguity, helping to reduce the deviation in production capacity prediction and improve simulation accuracy.

[0151] In step S103, the model calculation domain is discretized by mesh and the periphery of the mass flow inlet and pressure outlet is refined by mesh to obtain multiple mesh elements.

[0152] The above steps divide the continuous model computational domain into multiple discrete mesh elements with specific shapes, such as triangles or quadrilaterals. For example, mesh discretization can be performed using the finite element method, dividing a square computational domain with a side length of 1000 mm into multiple triangular mesh elements, each with defined nodal coordinates and physical properties. For instance, for the computational domain of the blocky reservoir in the Songliao Basin, triangular meshing is used to generate 27,500 mesh elements, each corresponding to 255,200 degrees of freedom. These discrete elements can serve as the basic computational units for solving equations such as heat conduction and the Navier-Stokes equations, allowing for the approximate calculation of continuous physical field variables, such as temperature, pressure, and flow velocity, within each mesh element, thus achieving quantitative analysis of the entire computational domain.

[0153] Here, the mesh around the mass flow inlet and pressure outlet is refined. That is, the key areas of fluid inflow and outflow, such as the inlet and outlet, are discretized with smaller mesh cells to improve the calculation accuracy of the area and accurately capture drastic changes in parameters such as fluid velocity and pressure.

[0154] For example, such as Figure 3 As shown, mesh refinement is performed around the mass flow inlet 203 (e.g., the lower left 80mm segment) and the pressure outlet 204 (e.g., the upper right 80mm segment), reducing the mesh cell size from the conventional 5mm to 2mm to create a denser mesh region. For example, around the inlet, when fluid is injected from the outside into the throat, the flow velocity changes drastically due to the narrow channel, rapidly increasing from 0 to hundreds of mm / s. A refined mesh can more accurately calculate the velocity gradient, avoiding distortion caused by an overly coarse mesh. Around the outlet, the pressure drops from reservoir pressure to outlet pressure, resulting in a large pressure gradient. A refined mesh can accurately capture this sudden pressure drop, ensuring the accuracy of the continuity equation and the Navier-Stokes equations, making the simulation results more reliable. Furthermore, as shown... Figure 4 As shown, when the degrees of freedom reach 255,200 (corresponding to 27,500 grid elements), the gas phase saturation tends to stabilize, and the computation time remains within an acceptable range.

[0155] In this model, physical parameters within each grid cell, such as temperature, pressure, and phase saturation, are assumed to be uniform or distributed according to a specific pattern. The physical field distribution of the entire computational domain is obtained by solving the algebraic equations of each grid cell. For example, a grid cell can be a triangular discrete element, with each cell containing three nodes. The physical parameters between nodes are correlated through interpolation functions, which can be used to store and calculate parameters such as temperature, pressure, and phase saturation within the cell. For instance, in simulating gas-phase displacement, each grid cell records the current gas phase field variables, temperature T, pressure p, and flow velocity u. By solving the Cahn-Hilliard equation and Navier-Stokes equation for each grid cell, the changes in physical parameters within the cell are obtained. Then, through information transfer between grid cells, the fluid transport patterns of the entire model's computational domain are integrated to obtain the overall model's computational domain. For example, the flow velocity in grid cells around the inlet is generally higher than that in the central cell of the computational domain, reflecting the local high-speed flow characteristics during fluid injection.

[0156] In step S104, environmental parameters and initial condition parameters are obtained. The environmental parameters include the physical properties of the rock matrix, the physical properties of the multiphase fluid, the initial saturation of the multiphase fluid, and the interfacial tension coefficient between different phase fluids. The initial condition parameters include the initial reservoir temperature and the initial reservoir pressure.

[0157] Among them, environmental parameters refer to the set of parameters describing the reservoir environment and fluid characteristics. They are the basic inputs for multiphase numerical model calculations and determine the degree to which the model restores the actual reservoir conditions. These parameters include the physical properties of the rock matrix, the physical properties of the multiphase fluids, the initial saturation of the multiphase fluids, and the interfacial tension coefficients between different phase fluids.

[0158] Specifically, the physical properties of the rock matrix include its specific heat capacity, density, thermal conductivity, porosity, coefficient of thermal expansion, Poisson's ratio, and Young's modulus; the physical properties of the multiphase fluid include those of the aqueous phase, oil phase, and gas phase. The aqueous phase physical properties include water density, specific heat capacity, thermal conductivity, and hydrodynamic viscosity; the oil phase physical properties include specific heat capacity, density, thermal conductivity, and dynamic viscosity; and the gas phase physical properties include density, specific heat capacity, thermal conductivity, and dynamic viscosity. The interfacial tension coefficients between different phase fluids include the oil-water tension coefficient, the oil-gas tension coefficient, and the water-gas tension coefficient.

[0159] Optionally, the specific heat capacity of the rock matrix is Density is Thermal conductivity is Porosity is 0.65, and coefficient of thermal expansion is... Poisson's ratio is 0.25, and Young's modulus is... ; and / or, water density, specific heat capacity, thermal conductivity, and hydrodynamic viscosity change with temperature; and / or, oil specific heat capacity is The density of the oil phase is The thermal conductivity of the oil phase is The dynamic viscosity of the oil phase changes with temperature; and / or, the dynamic viscosity of the gas phase is... The relative heat capacity of gas is The density of the gas phase is The thermal conductivity of the gas phase is The oil-water tension coefficient is The oil-gas tension coefficient is The water-air tension coefficient is .

[0160] For example, regarding the Qingshankou Formation reservoir in the Songliao Basin, the specific heat capacity of the rock matrix was measured through core experiments to be... Density is The specific heat capacity of oil was measured through fluid experiments. The dynamic viscosity of the gas phase is The initial saturation of the multiphase fluid is determined by combining well logging data, such as 65% for the oil phase, 25% for the water phase, and 10% for the gas phase. These parameters together constitute the environmental parameters of the model, ensuring that the model is consistent with the actual reservoir conditions.

[0161] Optionally, the physical properties of the rock matrix refer to parameters describing the physical characteristics of the rock matrix. For example, the physical properties of the rock matrix can be quantitative indicators reflecting the physical properties of the rock, obtained through laboratory core experiments and geological tests. For instance, the physical properties of the rock matrix of unconventional oil and gas reservoirs in the Songliao Basin are set as follows: specific heat capacity is... Density is Thermal conductivity is Porosity is 0.65, and coefficient of thermal expansion is... Poisson's ratio is 0.25, and Young's modulus is... These parameters are used to calculate heat transfer using the heat conduction equation, rock deformation using the mechanical equilibrium equation, and fluid flow using the seepage equation, thereby accurately simulating multi-field coupled processes.

[0162] Optionally, the physical properties of multiphase fluids can be parameters describing the physical characteristics of three-phase fluids (oil, gas, and water), and their values ​​can directly affect the fluid's flow and heat exchange behavior. For example, the physical properties of multiphase fluids can be quantitative indicators reflecting the fluid's physical characteristics, obtained through laboratory fluid testing. Some parameters, such as the properties of the aqueous phase, change with temperature. For instance, among the properties of the aqueous phase, the dynamic viscosity of water... According to the following formula, it changes with temperature:

[0163]

[0164] density of water It changes according to the following formula:

[0165]

[0166] The specific heat capacity of water under constant pressure changes according to the following formula:

[0167]

[0168] The thermal conductivity of water varies according to the following formula:

[0169]

[0170] In the oil phase physical properties parameters, the oil phase heat capacity is: The density of the oil phase is The thermal conductivity of the oil phase is The dynamic viscosity of the oil phase decreases with increasing temperature. The gas phase physical properties are constants, such as the dynamic viscosity of the gas phase. The relative heat capacity of gas is The density of the gas phase is The thermal conductivity of the gas phase is These parameters are used to calculate the fluid's flow resistance and heat exchange efficiency, and directly affect the flow velocity and temperature field distribution.

[0171] Optionally, the initial saturation of the multiphase fluid refers to the volume percentage of oil, gas, and water phases within the reservoir pore throat at the start of the simulation (t=0), representing the initial state of the multiphase fluid migration simulation. Specifically, the initial saturation of the multiphase fluid is a quantitative indicator reflecting the initial distribution of each phase, set based on the original fluid-bearing state of the reservoir. For example, in the simulation of layered reservoirs in the Songliao Basin, the initial saturation of the oil phase within the fractures is set to 65%, the water phase to 25%, and the gas phase to 10%. That is, at the start of the simulation, 65% of the pore throat volume is occupied by the oil phase, 25% by the water phase, and 10% by the gas phase. This initial state provides a starting baseline for the subsequent gas phase injection process that displaces the oil and water phases. By comparing the saturation changes at different times, the displacement efficiency and migration patterns of the gas phase can be analyzed.

[0172] Optionally, the interfacial tension coefficient between different phase fluids is a parameter describing the magnitude of molecular forces at the interface of oil-water, oil-gas, and water-gas two-phase fluids. It directly affects the interface morphology and migration resistance of multiphase fluids and is a key parameter for calculating interface evolution using the Cahn-Hilliard equation. Specifically, the interfacial tension coefficient is a quantitative index reflecting the interfacial forces between two-phase fluids, obtained through laboratory interfacial tension tests. For example, in the Songliao Basin reservoir simulation, the oil-water interfacial tension coefficient was set to 35 mN / m, the oil-gas interfacial tension coefficient to 15 mN / m, and the water-gas interfacial tension coefficient to 45 mN / m. Thus, the larger the interfacial tension coefficient, the greater the migration resistance at the multiphase interface. If the oil-water interfacial tension coefficient is greater than the oil-gas interfacial tension coefficient, the resistance to the gas phase displacing the oil phase is less than the resistance to the gas phase displacing the water phase, making it easier for the gas phase to form finger-like propagation in the oil phase region.

[0173] Furthermore, the initial condition parameters are the physical state parameters of the reservoir at the start of the simulation, including the initial temperature and initial pressure, which serve as the starting reference for the temperature and pressure field calculations. For example, the initial reservoir temperature of unconventional oil and gas reservoirs in the Songliao Basin can be set to 90℃ and the initial reservoir pressure can be set to 30MPa. At the start of the simulation, the temperature is uniformly 90℃ and the pressure is uniformly 30MPa throughout the computational domain. As the simulation progresses, the temperature and pressure of the injected fluid at the inlet will change the initial state within the computational domain. By tracking the changes in temperature and pressure, the influence of thermal-fluid-solid coupling on the reservoir state can be analyzed.

[0174] Here, the initial reservoir temperature is the overall temperature within the reservoir at the start of the simulation. It is the initial value for the evolution of the temperature field and determines the initial physical properties of the fluid and the initial thermal state of the rock. For example, the Qingshankou Formation reservoir in the Songliao Basin has a large burial depth and a high geothermal gradient. The initial reservoir temperature can be set to 90℃. This temperature determines the initial viscosity of the water phase and the initial viscosity of the oil phase at the start of the simulation, and also determines the initial thermal expansion state of the rock matrix, providing a starting reference for the subsequent simulation of temperature field changes. The initial reservoir pressure is the overall pore pressure within the reservoir at the start of the simulation. It is the initial value for the evolution of the pressure field and determines the initial effective stress state of the rock and the initial flow dynamics of the fluid. For example, due to the large burial depth of unconventional oil and gas reservoirs in the Songliao Basin, the initial reservoir pressure is set to 30MPa. This pressure determines the initial effective stress of the rock at the start of the simulation, thereby determining the initial deformation state of the matrix particles and the pore throat size. At the same time, the initial pressure (30MPa) and the outlet pressure (28MPa) form an initial pressure gradient, providing initial dynamics for fluid flow.

[0175] In step S105, environmental parameters and initial condition parameters are input into the multiphase numerical model, and the multiphase numerical model is controlled to perform discrete calculations in each grid cell. It is then determined whether the convergence condition is met. If the convergence condition is met, the physical property parameters of the multiphase fluid are updated according to the current temperature and pressure, and the numerical simulation results of the unconventional oil and gas reservoir are output. The numerical simulation results include phase saturation, temperature, pressure, and flow rate.

[0176] Inputting environmental parameters and initial condition parameters into the multiphase numerical model is a prerequisite for the model to perform calculations. This involves inputting parameters such as the obtained rock matrix physical properties, multiphase fluid physical properties, initial saturation of the multiphase fluid, interfacial tension coefficients between different phase fluids, initial reservoir temperature, and initial reservoir pressure into the multiphase numerical model to ensure that the multiphase numerical model can perform simulations based on actual reservoir conditions.

[0177] For example, environmental parameters of the Songliao Basin reservoir, such as the specific heat capacity of rocks, can be input into numerical simulation software, such as COMSOL Multiphysics. Oil-gas interfacial tension Initial condition parameters, such as initial reservoir temperature of 90℃ and initial reservoir pressure of 30MPa, are entered one by one into the multiphase numerical model. The numerical simulation software automatically assigns these parameters to the corresponding grid cells and equations, such as the rock thermal conductivity. Substitute the initial saturation into the initial phase field variables of the Cahn-Hilliard equation by substituting it into the heat conduction equation.

[0178] In this method, the integrated physical equations are solved within each discrete grid cell using the finite element method to obtain physical parameters such as temperature, pressure, saturation, and flow velocity within each grid cell. In an optional implementation, COMSOL Multiphysics combined with MATLAB can be used for discrete calculations. Specifically, environmental parameters and the initial condition parameters are input into the multiphase numerical model. Numerical simulation software combined with MATLAB is used to perform discrete calculations on the multiphase numerical model within each grid cell, and the convergence condition is determined. Here, within each grid cell, COMSOL Multiphysics uses the finite element method to solve the heat conduction equation to obtain the temperature distribution, solves the Navier-Stokes equation to obtain the flow velocity distribution, and solves the Cahn-Hilliard equation to obtain the phase saturation distribution. MATLAB reads the calculation data in real time and performs statistical analysis on the parameters within the cell to ensure the efficiency and accuracy of the calculation process. For the layered reservoirs of the Songliao Basin, this calculation process can accurately capture the flow velocity differences between fine-grained and coarse-grained layers, with the flow velocity in the fine-grained layer cells being lower than that in the coarse-grained layer cells, reflecting the influence of the layered structure on fluid transport.

[0179] Specifically, the convergence criteria are set as follows: total computation time of 60s, time step of 0.05s, and relative tolerance of 0.001. This means the calculation process proceeds step-by-step in 0.05s time steps until the total time reaches 60s, and the difference between the calculated and iterated values ​​of the residuals within each time step, such as flow velocity and pressure, is less than 0.001. For example, in simulating gas-phase displacement, within each time step (0.05s), the model iteratively solves each equation. If the relative residual of the flow velocity within a certain step is 0.0008 (less than 0.001), the calculation is considered converged for that step, and the process proceeds to the next time step. If the residual is 0.0012 (greater than 0.001), the process returns to iterate the finite element method and stiffness matrix, adjusts the mesh accuracy or iterative algorithm, and recalculates until the convergence criteria are met.

[0180] Specifically, if the convergence condition is not met, the process returns to the previous step to iterate the finite element and stiffness matrix again until the multiphase numerical model meets the convergence condition. That is, when the calculated residual does not meet the convergence criterion, the finite element discretization parameters, such as the mesh size and stiffness matrix, such as the equation coefficient matrix, are adjusted, and the discretization calculation is performed again until the convergence condition is met, ensuring the reliability of the calculation results.

[0181] In an optional implementation, if the calculation at a certain time step fails to meet the convergence condition, the model automatically returns to the mesh discretization step, adjusts the mesh size around the inlet, reassembles the stiffness matrix (such as the coefficient matrix of the equation), and performs discretization calculations again. For example, when simulating fluid migration in layered reservoirs, the large difference in physical properties between fine-grained and coarse-grained layers can easily lead to non-convergence. In this case, by refining the mesh at the layer interface, the calculation accuracy of the stiffness matrix can be improved, reducing the residual to 0.0009, satisfying the convergence condition, and ensuring that the calculation results accurately reflect the abrupt changes in fluid velocity at the layer interface.

[0182] Furthermore, once the multiphase numerical model meets the convergence condition, the physical properties of the multiphase fluid are dynamically updated based on the current temperature and pressure, and the phase saturation, temperature, pressure, and flow rate at each time point are output as numerical simulation results for unconventional oil and gas reservoirs.

[0183] Here, once the convergence condition is met, the fluid properties are recalculated based on the temperature and pressure of each grid cell within the reservoir at the current moment, achieving dynamic updates of fluid properties and reflecting the influence of multi-field coupling on fluid characteristics. For example, after convergence at a certain time step, the multiphase numerical model updates the fluid properties according to a preset formula based on the temperature T and pressure p of the current grid cell. For instance, when the temperature of the current cell rises from an initial 90℃ to 95℃, the viscosity of the aqueous phase is recalculated; simultaneously, the gas phase density is finely adjusted based on the current pressure, and the updated properties are substituted into the calculation for the next time step, ensuring that the model can reflect the influence of temperature and pressure on fluid flow in real time.

[0184] Furthermore, the calculated physical parameters within the reservoir, such as phase saturation, temperature, pressure, and flow rate, are output in data or visualization form. For example, simulation results can be output in the form of data tables and contour maps. The data tables record the phase saturation, temperature, pressure, and flow rate at the center point of the computational domain at each time step (e.g., 1s, 5s, 15s, 35s, 60s), while the contour maps visually display the parameter distribution throughout the computational domain. For example, as... Figure 5 As shown, initially, the gas phase saturation is low, only 0.10. With the passage of time, the gas phase gradually expands from the inlet to the outlet. In the region far from the inlet-outlet diagonal, the gas phase saturation change is not significant. The output gas phase saturation distribution map clearly shows the finger-like propagation process of the gas phase from the inlet. At 60 seconds, the gas phase saturation near the outlet reaches 0.8. Figure 6 As shown, the velocity distribution map displays high-speed flows near the inlet and outlet. It can be observed that at 10 s, due to the oil and water phases occupying the pore throat, the high viscosity of the mixed fluid increases flow resistance, making it difficult to quickly establish fluid flow channels, thus creating a larger affected area. By 50 s, the gas phase occupies the dominant position in the pore throat, forming high-speed flow channels on both sides of the inlet-outlet diagonal, resulting in a corresponding reduction in the affected area. These results provide quantitative evidence for analyzing the gas phase displacement efficiency and dominant migration channels of the Songliao Basin reservoirs, and can guide on-site adjustments to injection flow rates or well location layouts to improve development efficiency.

[0185] For example, such as Figure 7 and Figure 8 As shown, the saturation of the oil and water phases in the pore throat gradually decreases, exhibiting an overall reverse finger-like phenomenon. Some oil and water phases remain trapped in the pore throat; for example, significant oil phase accumulation occurs near the upper and lower boundaries, with saturation even exceeding the initial value. Similarly, the water phase saturation in the region far from the inlet-outlet diagonal also exceeds the initial value. This indicates that with the injection of gas phase, some oil and water phases are displaced but cannot flow out, gradually accumulating in areas far from the main channel. This further illustrates that in actual production, key measures are to increase the horizontal drainage area and expand the vertical control volume.

[0186] To quantitatively describe the characteristics of phase saturation within the pore throat, a feature point and two feature lines were selected as reference points and reference lines within the model's computational domain: reference point c, reference line df, and reference line gh. Here, g, h, d, and f represent the endpoint markers of the reference lines set within the model's computational domain, used to quantitatively analyze the spatiotemporal variation of fluid phase saturation (e.g., oil, water, and gas phases). Figure 9(a) shows this. Figure 9(b) illustrates the changes in phase saturation at reference point C. As can be seen from the figure, the oil phase (phase A) and water phase (phase B) generally show a decreasing trend, while the gas phase (phase C) shows an increasing trend. The curves are not smooth but exhibit significant fluctuations and bends. This phenomenon is attributed to the aggregation and displacement of the mixed fluid at this point until one phase completely occupies that position. Figure 9(c) shows the evolution of oil phase (phase A) saturation along the reference line, with the curve exhibiting regular fluctuations. Fluid accumulates on the upstream side of the matrix particles, while it is squeezed out on the downstream side. This phenomenon is mainly caused by the disturbance caused by the fluid passing through the narrow throat.

[0187] Furthermore, as shown in Figure 10, the maximum flow velocities at 10s, 30s, and 50s are 492mm / s, 643mm / s, and 678mm / s, respectively, all occurring at the inlet or outlet due to the narrow channels at these locations. Because the initial affected area is relatively large, the average flow velocity at 10s is 32.5mm / s. The change in the seepage area is more pronounced in the early stages, from 10s... Reduced to 30 seconds The proportion of the seepage zone to the total pore throat area decreased from 57.7% to 43.4% (a reduction of 14.3%). At 50 seconds, the seepage zone was... The flow rate decreased by only 3.3% compared to 30s. Interestingly, at the geometric center, upper boundary, and right boundary, the lower flow rate caused the oil and water phases to accumulate at these locations.

[0188] To analyze the influence of stress and temperature fields on the seepage field, the saturation distribution of the gas phase (C phase) was compared with and without considering the temperature-pressure field. Figure 11(a) shows that without considering multiphysics coupling, the area covered by the gas phase expands, highlighted by a yellow circle, because the water and oil phases have lower viscosity and higher fluidity at the initial reservoir temperature. Figure 11(b) quantitatively shows the difference in area under the two conditions, with gas phase saturation at 0.5 as the dividing point. Without considering multiphysics coupling, the affected area is larger, and the difference peaks at 43 s. Under the conditions of this study, the contrast effect is not significant, mainly due to the limited size of the computational domain. However, this effect should not be ignored in practical engineering applications.

[0189] Furthermore, such as Figure 12 As shown, mass flow rate is a key factor affecting fluid transport. This study designed four different flow rates (0.5 g / s, 1.0 g / s, 1.5 g / s, and 2.0 g / s) and analyzed the results. At low flow rates, the sweep range of the gas phase is small, and it has not yet reached the outlet after 50 s. At high flow rates, the gas phase displacement effect is more significant; for example, at a flow rate of 1.0 g / s, the gas phase front reaches the outlet after 50 s.

[0190] Specifically, Figure 13 The changes in gas phase (C phase) saturation regions over time are shown, with 0.5 as the dividing point. At 10 s, the regions corresponding to the four flow rates are 0.58 × 10⁵ mm², 1.15 × 10⁵ mm², 1.72 × 10⁵ mm², and 2.30 × 10⁵ mm² (from smallest to largest). By 30 s, these values ​​increase to 1.66 × 10⁵ mm², 3.28 × 10⁵ mm², 4.13 × 10⁵ mm², and 4.29 × 10⁵ mm², with increases of 183.3%, 185.0%, 140.1%, and 86.3%, respectively. In the later stages, the increase in region area slows down, especially at high flow rates. For example, at flow rates of 1.5 g / s and 2.0 g / s, the corresponding areas at 50 s are 4.34 × 10⁵ mm² and 4.39 × 10⁵ mm², respectively, representing increases of only 5.0% and 2.5% compared to 30 s. When conducting unconventional oil and gas production, the performance of surface pumps must be considered, with a focus on designing and adjusting the mass flow rate.

[0191] In this application, the differences in rock particle diameter represent the physical characteristics of different unconventional fine-grained oil and gas reservoirs in terms of porosity and permeability. Therefore, under the condition of an initial velocity of 2.0 g / s, this study designed four different diameters of 30 mm, 35 mm, 40 mm, and 45 mm for evaluation and analysis. Figure 14 As shown, with smaller particle diameters, the swept area of ​​the gas phase is larger in both the horizontal and vertical directions. At 10 s, the 30 mm gas phase is about to reach the upper and right boundaries of the computational domain. This phenomenon occurs because smaller particle diameters correspond to greater porosity and permeability, resulting in lower resistance to the flow of the mixed fluid and smoother flow. With larger particle diameters (35 mm and 40 mm), although the final swept area is smaller than that corresponding to a 30 mm particle diameter, its proportion is higher.

[0192] It is important to note that for particles with a diameter of 45 mm, the pore throat is relatively small, with the narrowest point being only 5 mm. This results in low reservoir permeability, making fluid flow extremely difficult and prone to short-circuiting, ultimately reducing the sweep range. The initially larger sweep range is due to the gas phase occupying a larger area in a short time due to the smaller pore throat at the same injection mass flow rate.

[0193] Figure 15 shows the change in gas phase (C phase) saturation region over time for different particle diameters (with 0.5 as the cutoff point), which strongly supports the above analysis. At 10 s, the regions corresponding to the four particle sizes are 1.45 × 10⁵ mm², 1.37 × 10⁵ mm², 1.15 × 10⁵ mm², and 0.95 × 10⁵ mm² (from smallest to largest). By 50 s, these values ​​increase to 5.49 × 10⁵ mm², 4.97 × 10⁵ mm², 4.16 × 10⁵ mm², and 2.28 × 10⁵ mm², respectively. Although the influence region of larger particle sizes is smaller, the proportion of the region is higher. For example, at 50 s, the proportions of the region for 40 mm and 35 mm particles reach 81.1% and 79.5%, respectively, while the proportion for 30 mm particles is only 75.8%. Furthermore, the data corresponding to the 45 mm particle size indicate that there is a clear threshold for the influence of particle size on fluid transport. Therefore, necessary porosity and permeability tests should be conducted before a design scheme is developed.

[0194] The lithology of unconventional fine-grained oil and gas reservoirs mainly includes massive and layered mudstone. The first two sections of this report analyzed massive reservoirs, and this section will focus on analyzing the characteristics of different sedimentary layers in layered reservoirs, including bedding orientation (horizontal or vertical) and permeability (particle diameter) of each layer. Figure 16 A schematic diagram of the geometric model is provided, where M, C, and F represent medium-grained, coarse-grained, and fine-grained reservoirs, respectively, with corresponding grain diameters of 40 mm, 45 mm, and 30 mm. The geometric model is divided into two main types: horizontally layered and vertically layered reservoirs. Each type contains four subtypes with different grain diameters, for a total of eight subtypes.

[0195] Figure 17 The fluid migration characteristics of different types of vertical reservoirs are demonstrated. At 5 s, the expansion direction of the gas phase differs significantly across different bedding planes. When transitioning from meso-grained to coarse-grained layers, the gas phase primarily migrates vertically. Conversely, when transitioning from meso-grained to fine-grained layers, the gas phase primarily migrates horizontally. Furthermore, gas phase accumulation is observed in the transition zones from coarse-grained to meso-grained layers and from meso-grained to fine-grained layers. This indicates that within the transition zones, as the velocity of high-velocity fluid decreases (from low-permeability layers to high-permeability layers), the displacement effect of the gas phase on the oil and water phases is enhanced.

[0196] At 35 s, the gas phase almost completely occupied the pore throat region. It can be concluded that in the medium- and fine-grained layers, the gas phase effectively replaced the oil and water phases, with the gas phase almost completely occupying the fine-grained layer. However, in the coarse-grained layer, a large amount of non-gas phase remained. Significant short-circuiting of fluid transport was observed; for example, in the contour maps corresponding to the MFMCM and MMFFM scenarios, large areas of the medium-grained layer were short-circuited by fluid channels in the fine-grained layer. The gas phase contour map characteristics at 60 s were not significantly different from those at 35 s, with only a slight increase in gas phase saturation in certain regions.

[0197] To visually describe the fluid transport path, Figure 18 The velocity characteristics of the fluid flow field at 60 s are shown. It can be observed that the overall velocity distribution is consistent with the gas phase saturation distribution contour map, but some differences remain. For example, in the simulation case corresponding to the MMFMFM stratification, the vertical velocity distribution in the second column of the fine-grained stratification is relatively sparse, indicating that this region mainly experiences gas phase displacement in the early stages. Once high-velocity seepage channels are formed, the fluid is bypassed. The velocity contour map also shows that the gas phase sweeping range is wider in low-permeability reservoirs.

[0198] Figure 19(a) quantitatively characterizes the above analysis, showing the percentage of the overall pore-throat region with a gas phase saturation greater than 0.5 at 5 s, 35 s, and 60 s. The proportion of regions under each type increases over time. Taking the homogeneous stratification case as an example, the region percentages at each time point are 11.9%, 73.6%, and 81.4%, respectively. Among the comparisons of different types, only the medium- and fine-grained stratified reservoirs have the smallest region proportions at 60 s, both below 60%, while the corresponding regions in other cases are all above 75%. This is because high-permeability layers are more likely to form high-velocity fluid channels, thus reducing the effectiveness of gas phase displacement.

[0199] like Figure 20 As shown, the gas phase distribution characteristics of horizontal and vertical stratification are similar. When encountering low-permeability, large-particle stratification, the gas phase mainly migrates along the original stratification direction (horizontal); while when encountering high-permeability, small-particle reservoirs, the gas phase migration direction turns to the high-permeability layer (vertical direction). Figure 21 Figure 19(b) shows the fluid velocity distribution in layers with different particle diameters at 60 s (horizontal layered), while Figure 19(b) shows the percentage of the gas phase (saturation greater than 0.5) at different time points. Similar to vertically layered reservoirs, the MCMCM reservoirs had the highest proportions at 5 s and 60 s, at 24.4% and 86.7%, respectively.

[0200] The unconventional oil and gas reservoir thermal-fluid-solid coupling numerical simulation method provided in this application constructs a multiphase numerical model that couples temperature, seepage, and stress fields and integrates multiple key equations. This fully realizes dynamic thermal-fluid-solid coupling, avoiding the shortcomings of insufficient multi-field coupling in existing technologies. Furthermore, by combining precise design of the model's computational domain, differentiated mesh discretization, and comprehensive acquisition of environmental and initial condition parameters, it effectively recreates in-situ geological conditions, adapts to the strong heterogeneity of unconventional reservoirs, and accurately simulates the micro-transportation process of multiphase fluids. This significantly improves the accuracy of production capacity prediction, reduces prediction bias, and can accurately guide the development of unconventional oil and gas reservoirs.

[0201] Based on the same inventive concept, this application also provides an unconventional oil and gas reservoir thermal-fluid-structure coupling numerical simulation system corresponding to the unconventional oil and gas reservoir thermal-fluid-structure coupling numerical simulation method. Since the principle of the system in this application is similar to the unconventional oil and gas reservoir thermal-fluid-structure coupling numerical simulation method described above in this application, the implementation of the device can refer to the implementation of the method, and the repeated parts will not be described again.

[0202] Please see Figure 22 , Figure 22 This is a schematic diagram of the structure of an unconventional oil and gas reservoir thermal-fluid-structure interaction numerical simulation system provided in an embodiment of this application. Figure 22 As shown, the system 2200 includes:

[0203] The model building module 2201 is used to build a multiphase numerical model that couples the temperature field, seepage field and stress field; the multiphase numerical model integrates the Cahn-Hilliard equation, the heat conduction equation, the continuity equation, the Navier-Stokes equation and the mechanical equilibrium equation.

[0204] The computational domain determination module 2202 is used to determine the computational domain of the multiphase numerical model. The computational domain includes multiple matrix particles and pores and throats between adjacent matrix particles. The first end of the computational domain is provided with a mass flow inlet, and the second end of the computational domain is provided with a pressure outlet. The pores and throats are used to store and transport fluid. The first end and the second end are located diagonally opposite to each other in the computational domain, and all boundaries of the computational domain are non-flow boundaries and adiabatic boundaries.

[0205] The mesh generation module 2203 is used to perform mesh discretization processing on the model calculation domain and mesh refinement processing on the periphery of the mass flow inlet and the pressure outlet to obtain multiple mesh cells;

[0206] The parameter acquisition module 2204 is used to acquire environmental parameters and initial condition parameters. The environmental parameters include the physical properties of the rock matrix, the physical properties of the multiphase fluid, the initial saturation of the multiphase fluid, and the interfacial tension coefficient between different phase fluids. The initial condition parameters include the initial reservoir temperature and the initial reservoir pressure.

[0207] The numerical simulation module 2205 is used to input the environmental parameters and the initial condition parameters into the multiphase numerical model, control the multiphase numerical model to perform discrete calculations in each grid cell, and determine whether the convergence condition is met. If the convergence condition is met, the physical property parameters of the multiphase fluid are updated according to the current temperature and pressure, and the numerical simulation results of the unconventional oil and gas reservoir are output. The numerical simulation results include phase saturation, temperature, pressure and flow rate.

[0208] The system provided in this application embodiment, by constructing a multiphase numerical model that couples temperature field, seepage field, and stress field and integrates multiple key equations, can fully realize full dynamic coupling of heat, fluid, and solid, avoiding the shortcomings of insufficient multi-field coupling in existing technologies. At the same time, by combining the precise design of the model's computational domain, the differentiated discretization of the grid, and the comprehensive acquisition of environmental parameters and initial condition parameters, it can effectively restore in-situ geological conditions, adapt to the strong heterogeneity of unconventional reservoirs, and thus accurately simulate the micro-transport process of multiphase fluids, significantly improving the accuracy of production capacity prediction, reducing prediction bias, and accurately guiding the development of unconventional oil and gas reservoirs.

[0209] Please see Figure 23 , Figure 23 This is a schematic diagram of the structure of an electronic device provided in an embodiment of this application. Figure 23 As shown, the electronic device 2300 includes a processor 2310, a memory 2320, and a bus 2330.

[0210] The memory 2320 stores machine-readable instructions executable by the processor 2310. When the electronic device 2300 is running, the processor 2310 and the memory 2320 communicate via the bus 2330. When the machine-readable instructions are executed by the processor 2310, they can perform the operations described above. Figure 1 The steps of the unconventional oil and gas reservoir thermal-fluid-structure interaction numerical simulation method in the method embodiment shown are described in detail in the method embodiment, and will not be repeated here.

[0211] This application also provides a computer-readable storage medium storing a computer program, which, when executed by a processor, can perform the above-described actions. Figure 1 The steps of the unconventional oil and gas reservoir thermal-fluid-structure interaction numerical simulation method in the method embodiment shown are described in detail in the method embodiment, and will not be repeated here.

[0212] Those skilled in the art will understand that, for the sake of convenience and brevity, the specific working processes of the systems, devices, and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0213] In the several embodiments provided in this application, it should be understood that the disclosed systems, apparatuses, and methods can be implemented in other ways. The apparatus embodiments described above are merely illustrative. For example, the division of units is only a logical functional division, and in actual implementation, there may be other division methods. Furthermore, multiple units or components may be combined or integrated into another system, or some features may be ignored or not executed. Additionally, the shown or discussed mutual couplings, direct couplings, or communication connections may be through some communication interfaces; indirect couplings or communication connections between devices or units may be electrical, mechanical, or other forms.

[0214] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment according to actual needs.

[0215] In addition, the functional units in the various embodiments of this application can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit.

[0216] If the aforementioned functions are implemented as software functional units and sold or used as independent products, they can be stored in a processor-executable, non-volatile, computer-readable storage medium. Based on this understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, or a portion of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of this application. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0217] Finally, it should be noted that the above-described embodiments are merely specific implementations of this application, used to illustrate the technical solutions of this application, and not to limit them. The scope of protection of this application is not limited thereto. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that any person skilled in the art can still modify or easily conceive of changes to the technical solutions described in the foregoing embodiments, or make equivalent substitutions for some of the technical features, within the scope of the technology disclosed in this application. Such modifications, changes, or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of this application, and should all be covered within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A numerical simulation method for thermo-fluid-structure interaction in unconventional oil and gas reservoirs, characterized in that, include: A multiphase numerical model is constructed that couples the temperature field, seepage field, and stress field; the multiphase numerical model integrates the Cahn-Hilliard equation, the heat conduction equation, the continuity equation, the Navier-Stokes equation, and the mechanical equilibrium equation; The computational domain of the multiphase numerical model is determined. The computational domain includes multiple matrix particles and pores and throats between adjacent matrix particles. A mass flow inlet is provided at the first end of the computational domain, and a pressure outlet is provided at the second end of the computational domain. The pores and throats are used to store and transport fluid. The first end and the second end are located diagonally opposite to each other in the computational domain. All boundaries of the computational domain are non-flow boundaries and adiabatic boundaries. The computational domain of the model is discretized into a grid, and the area around the mass flow inlet and the pressure outlet is refined into a grid to obtain multiple grid cells. The environmental parameters and initial condition parameters are obtained. The environmental parameters include the physical properties of the rock matrix, the physical properties of the multiphase fluid, the initial saturation of the multiphase fluid, and the interfacial tension coefficient between different phase fluids. The initial condition parameters include the initial reservoir temperature and the initial reservoir pressure. The environmental parameters and the initial condition parameters are input into the multiphase numerical model. The multiphase numerical model is controlled to perform discrete calculations in each grid cell and to determine whether the convergence condition is met. If the convergence condition is met, the physical property parameters of the multiphase fluid are updated according to the current temperature and pressure, and the numerical simulation results of the unconventional oil and gas reservoir are output. The numerical simulation results include phase saturation, temperature, pressure and flow rate.

2. The method according to claim 1, characterized in that, The environmental parameters and initial condition parameters are input into the multiphase numerical model, which is then controlled to perform discrete calculations within each grid cell. The model is then checked for convergence conditions. If convergence conditions are met, the physical properties of the multiphase fluid are updated based on the current temperature and pressure, and the numerical simulation results of the unconventional oil and gas reservoir are output. These simulation results include phase saturation, temperature, pressure, and flow velocity. The environmental parameters and the initial condition parameters are input into the multiphase numerical model. Numerical simulation software and MATLAB are used to perform discrete calculations on the multiphase numerical model in each grid cell, and it is determined whether the convergence condition is met. If the convergence condition is not met, return to the previous step to iterate the finite element and stiffness matrix again until the multiphase numerical model meets the convergence condition. Once the multiphase numerical model meets the convergence condition, the physical properties of the multiphase fluid are dynamically updated based on the current temperature and pressure, and the phase saturation, temperature, pressure, and flow rate at each time point are output as numerical simulation results of unconventional oil and gas reservoirs.

3. The method according to claim 2, characterized in that, The convergence conditions include: a total computation time of 60s, a time step of 0.05s, and a relative tolerance of 0.

001.

4. The method according to claim 1, characterized in that, The Cahn-Hilliard equation is expressed by the following formula: in, Let t represent the phase field variable and t represent time. M0 represents the fluid velocity, and M0 represents the fluid flow regulation parameter. This represents capillary force, where i represents oil phase A, water phase B, and gas phase C, respectively. This represents the phase field variable of oil phase A. Represents the phase field variables of water phase B. Represents the phase field variables of gas phase C; This can be expressed by the following formula: in, Denotes the coefficients of the Cahn–Hilliard equation. Parameters that control the thickness of the interface; This represents the partial derivative of the free energy with respect to the order parameter. To express summation; ΣA, ΣB, and ΣC are represented by the following formulas: This represents the interfacial tension coefficient between the oil phase and the water phase. This represents the interfacial tension coefficient between the oil phase and the gas phase. This represents the interfacial tension coefficient between the aqueous phase and the gas phase. F(ψ) represents volume energy, expressed by the following formula: The additional free volume energy is represented by Λ; And / or, the heat conduction equation is expressed by the following formula: in, and These represent the effective volumetric heat capacity and effective thermal conductivity, respectively. in, These represent the density of the solid and the density of the fluid, respectively. and These represent the specific heat capacity of a solid and the specific heat capacity of a fluid, respectively. and Let Q represent the thermal conductivity of the solid and the thermal conductivity of the fluid, respectively, and let Q represent the energy source. Indicates porosity; The thermal conductivity of a fluid is expressed by the following formula: in, Indicates the density of the oil phase. This represents the density of the aqueous phase. This represents the density of the gas phase. Indicates the thermal conductivity of the oil phase. Indicates the thermal conductivity of the aqueous phase. Indicates the thermal conductivity of the gas phase; The specific heat capacity of a fluid is expressed by the following formula: in, Indicates the relative heat capacity of oil. This indicates the specific heat capacity of water. This indicates the specific heat capacity of the gas. And / or, the Navier-Stokes equations and the continuity equations are expressed by the following formulas: in, Let I represent the fluid density, I represent the identity matrix, and F represent the volume force vector. Represents volume force; K represents the viscous stress tensor, expressed by the following formula: The density and viscosity of a fluid mixture are expressed by the following formula: The viscosity of a fluid mixture is expressed by the following formula: in, Indicates the viscosity of the oil phase. Indicates the viscosity of the aqueous phase. Indicates the viscosity of the gas phase; Volume forces are expressed by the following formula: And / or, the mechanical equilibrium equations are expressed by the following formulas: Where E represents Young's modulus and v represents Poisson's ratio. This represents the displacement of the boundary under ground stress in the initial state. This indicates the boundary displacement caused by production. This represents the Biot–Willi coefficient. This represents the fluid pressure inside the pores. The symbol for Kronecker, T represents the coefficient of thermal expansion, and T represents temperature. Indicates the initial temperature. This represents the force per unit volume.

5. The method according to claim 1, characterized in that, The model calculation domain includes a blocky reservoir calculation domain and a layered reservoir calculation domain. The matrix particles in the blocky reservoir calculation domain have a single diameter, and the diameter of the matrix particles ranges from 30mm to 45mm. The center-to-center distance between adjacent matrix particles is 50mm. The layered reservoir computational domain includes a horizontally layered reservoir computational domain or a vertically layered reservoir computational domain. The horizontally layered reservoir computational domain or the vertically layered reservoir computational domain includes a coarse-grained reservoir computational domain, a medium-grained reservoir computational domain, and a fine-grained reservoir computational domain. The matrix particle diameter in the coarse-grained reservoir computational domain is 45 mm, the matrix particle diameter in the medium-grained reservoir computational domain is 40 mm, and the matrix particle diameter in the fine-grained reservoir computational domain is 30 mm.

6. The method according to claim 1, characterized in that, The mass flow inlet has a length of 80 mm, the pressure outlet has a length of 80 mm, and the outlet pressure is 2.0 MPa lower than the reservoir pressure.

7. The method according to claim 1, characterized in that, The physical properties of the rock matrix include its specific heat capacity, density, thermal conductivity, porosity, coefficient of thermal expansion, Poisson's ratio, and Young's modulus. The physical properties of the multiphase fluid include those of the aqueous phase, oil phase, and gas phase. The aqueous phase physical properties include water density, specific heat capacity, thermal conductivity, and hydrodynamic viscosity. The oil phase physical properties include specific heat capacity, density, thermal conductivity, and dynamic viscosity. The gas phase physical properties include gas density, specific heat capacity, thermal conductivity, and dynamic viscosity. The interfacial tension coefficients between the different phase fluids include oil-water tension coefficients, oil-gas tension coefficients, and water-gas tension coefficients.

8. The method according to claim 7, characterized in that, The specific heat capacity of the rock matrix is Density is Thermal conductivity is Porosity is 0.65, and coefficient of thermal expansion is... Poisson's ratio is 0.25, and Young's modulus is... ; And / or, water density, specific heat capacity, thermal conductivity, and hydrodynamic viscosity change with temperature; And / or, the relative heat capacity of oil is The density of the oil phase is The thermal conductivity of the oil phase is The dynamic viscosity of the oil phase changes with temperature; And / or, the gas phase dynamic viscosity is The relative heat capacity of gas is The density of the gas phase is The thermal conductivity of the gas phase is ; The oil-water tension coefficient is The oil-gas tension coefficient is The water-air tension coefficient is .

9. The method according to claim 1 or 7, characterized in that, The initial saturation of the multiphase fluid includes the initial saturation of the oil phase, the initial saturation of the water phase, and the initial saturation of the gas phase, wherein the initial saturation of the oil phase is 65%, the initial saturation of the water phase is 25%, and the initial saturation of the gas phase is 10%. The initial reservoir temperature is 90°C and the initial reservoir pressure is 30 MPa.

10. A numerical simulation system for thermo-fluid-structure interaction in unconventional oil and gas reservoirs, characterized in that, include: The model building module is used to construct a multiphase numerical model that couples the temperature field, seepage field, and stress field; the multiphase numerical model integrates the Cahn-Hilliard equation, the heat conduction equation, the continuity equation, the Navier-Stokes equation, and the mechanical equilibrium equation. The computational domain determination module is used to determine the computational domain of the multiphase numerical model. The computational domain includes multiple matrix particles and pores and throats between adjacent matrix particles. The first end of the computational domain is provided with a mass flow inlet, and the second end of the computational domain is provided with a pressure outlet. The pores and throats are used to store and transport fluid. The first end and the second end are located diagonally opposite to each other in the computational domain, and all boundaries of the computational domain are non-flow boundaries and adiabatic boundaries. The mesh generation module is used to discretize the model computation domain and refine the mesh around the mass flow inlet and the pressure outlet to obtain multiple mesh cells. The parameter acquisition module is used to acquire environmental parameters and initial condition parameters. The environmental parameters include the physical properties of the rock matrix, the physical properties of the multiphase fluid, the initial saturation of the multiphase fluid, and the interfacial tension coefficient between different phase fluids. The initial condition parameters include the initial reservoir temperature and the initial reservoir pressure. The numerical simulation module is used to input the environmental parameters and the initial condition parameters into the multiphase numerical model, control the multiphase numerical model to perform discrete calculations in each grid cell, and determine whether the convergence condition is met. If the convergence condition is met, the physical property parameters of the multiphase fluid are updated according to the current temperature and pressure, and the numerical simulation results of the unconventional oil and gas reservoir are output. The numerical simulation results include phase saturation, temperature, pressure and flow rate.

Images