Numerical simulation method for multiphase fluid stratum seepage

Abstract

The invention relates to a numerical simulation method for multiphase fluid stratum seepage, which belongs to the technical field of numerical simulation and comprises the following steps: S1, generating an unstructured grid by adopting a mode of combining a leading edge propulsion method and Delaunay triangulation; s2, performing dynamic adjustment through an adaptive grid encryption technology according to the change of the seepage physical quantity; s3, establishing a dynamic phase change model, and obtaining a function of representing the phase change rate as temperature, pressure and concentration of each phase; s4, the mass and energy change of each phase in the phase change process is calculated in real time through a kinetic equation; s5, processing dynamic evolution of a phase interface by adopting a phase field method, and describing distribution conditions of different phases through phase field variables; and S6, coupling the grid system obtained by dividing the multi-scale unstructured grid with the dynamic phase change model through a numerical solution algorithm. According to the method, the complex geometrical shape of the stratum can be better fitted through multi-scale unstructured grid division, and then the change of the phase change process along with time and space can be tracked in real time through the dynamic phase change model.

Description

Technical Field

[0001] The present invention relates to the technical field of numerical simulation, and particularly relates to a numerical simulation method for multi-phase fluid formation seepage. Background Art

[0002] During the process of oil and gas exploitation, there is complex seepage of multi-phase fluids such as oil, gas, and water inside the reservoir. Traditional numerical simulation methods have many limitations when facing such multi-phase flow problems. For example, due to the regular geometric characteristics of conventional structured grids, it is difficult to accurately fit complex formation shapes. When encountering complex geological structures such as faults and folds in the formation, the grid division accuracy is severely insufficient, resulting in a large deviation between the simulation results and the actual situation. Relevant research shows that under some complex formation conditions, when using structured grids for seepage simulation, the pressure calculation error can be as high as more than 20%.

[0003] For the phase change phenomenon in multi-phase fluids, existing static phase change models cannot reflect the dynamic changes of phase change over time and space in real time. During the process of natural gas exploitation, as the formation pressure decreases, condensate oil in the gas will gradually precipitate, and the static phase change model cannot accurately capture this dynamic process, resulting in a large deviation between the prediction of gas well productivity and the actual situation. Summary of the Invention

[0004] The purpose of the present invention is to address the problems in the background art and propose a numerical simulation method for multi-phase fluid formation seepage that can better fit complex formation geometries through multi-scale unstructured grid division and can track the changes of the phase change process over time and space in real time through a dynamic phase change model.

[0005] The technical solution of the present invention: A numerical simulation method for multi-phase fluid formation seepage, the method comprising the following steps:

[0006] S1. Generate an unstructured grid by combining the advancing front method and Delaunay triangulation;

[0007] S2. Dynamically adjust according to the change of seepage physical quantities through adaptive grid refinement technology;

[0008] S3. Establish a dynamic phase change model to obtain the phase change rate expressed as a function of temperature, pressure, and the concentration of each phase;

[0009] S4. Calculate the changes in the mass and energy of each phase during the phase change process in real time through kinetic equations;

[0010] S5. Adopt the phase field method to process the dynamic evolution of the phase interface and describe the distribution of different phases through phase field variables;

[0011] S6. Couple the grid system obtained by multi-scale unstructured grid division with the dynamic phase change model through a numerical solution algorithm;

[0012] S7. Discretize the mass conservation equation, momentum conservation equation, and energy conservation equation using the finite volume method.

[0013] Preferably, the advancing front method starts from the formation boundary and gradually advances inward to generate triangular grid elements. During the advancement process, the Delaunay triangulation property is used to ensure good quality of the generated grid elements;

[0014] In the advancing front method, new nodes are generated based on the optimization function

[0015]

[0016] to determine the position,

[0017] where θi is the interior angle of the triangular element;

[0018] Taking the pressure gradient ▽P as an example, for the adaptive grid refinement technique, when |▽P| > τ 1 , subdivision is carried out. The subdivision method is to divide the triangular element from the centroid to the midpoint connection of the three sides into four sub-triangular elements. When |▽P| < τ2, adjacent element merging is carried out;

[0019] τ1 is the refinement threshold;

[0020] τ2 is the coarsening threshold;

[0021] and τ2 < τ1.

[0022] Preferably, the adaptive grid refinement technique automatically refines the grid in the region where the physical quantities such as pressure gradient and flow velocity change violently, and uses a coarser grid in the region where the physical quantities change gently.

[0023] Preferably, the distribution of different phases is described by the phase field variable, and the phase field equation is solved to accurately capture the movement and deformation of the phase interface during the seepage process.

[0024] Preferably, the finite volume method integrates the control equations on the grid elements and discretizes the time and space derivatives using methods such as central difference.

[0025] Preferably, for the phase change between the gas phase g and the liquid phase l, the phase change rate

[0026] R gl = k1·f(P, T, C g , C l )

[0027] where k1 is the kinetic coefficient,

[0028] f(P, T, C g , C l )

[0029] is a function of pressure P, temperature T, gas-phase concentration C g and liquid-phase concentration Cl;

[0030] f(P, T, C g , C l ) = (P - P sat (T)) · C g · C l , P sat (T)

[0031] is the saturation pressure at the corresponding temperature.

[0032] Preferably, the dynamic phase transition model uses the phase field method to handle the dynamic evolution of the phase interface, describes the distribution of different phases through phase field variables, and solves the phase field equation

[0033]

[0034] to accurately capture the movement and deformation of the phase interface during the seepage process, where M is the mobility and F is the free energy functional

[0035]

[0036] is the local free energy density, is the interface thickness parameter.

[0037] Preferably, the solution process of the numerical solution algorithm includes initializing the calculation parameters, performing iterative solution within each time step according to the control equation and boundary conditions, and updating the grid system and variables related to the phase transition model during the solution process until the set simulation time or convergence condition is reached.

[0038] Preferably, the initialized calculation parameters include the initial pressure P0, initial temperature T0, fluid saturation S0 of the formation, and related parameters of the grid system and phase transition model, such as the kinetic coefficient k1 and mobility M in the dynamic phase transition model.

[0039] Preferably, in the data preparation stage, geological data such as the geometry of the formation, permeability distribution, porosity, etc., and physical property parameters such as the density, viscosity, and phase transition heat of the multiphase fluid are collected.

[0040] Compared with the prior art, the present invention has the following beneficial technical effects:

[0041] In the present invention, multi-scale unstructured mesh generation can better conform to complex formation geometries, and the simulation accuracy for irregular regions is improved by more than 50% compared to traditional structured meshes. The dynamic phase change model can track the changes of the phase change process over time and space in real time. Compared with the traditional static phase change model, the simulation error of the phase change process is reduced by about 40%, which more accurately reflects the actual situation of multiphase fluid seepage in the formation. The adaptive mesh refinement technology reduces unnecessary computational effort while ensuring simulation accuracy. By using coarser meshes in regions with gentle changes in physical quantities, the overall computational time is shortened by 30%-40% compared to full-region fine mesh simulation, improving the computational efficiency.

[0042] The above description is only an overview of the technical solution of the present invention. In order to understand the technical means of the present invention more clearly and to be able to implement it in accordance with the content of the specification, the following provides a detailed description of the preferred embodiments of the present invention in conjunction with the accompanying drawings. The specific implementation manners of the present invention are given in detail by the following embodiments and their accompanying drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] The drawings described herein are used to provide a further understanding of the present invention and form a part of this application. The schematic embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute an improper limitation of the present invention. In the drawings:

[0044] Figure 1 is a schematic flow chart of an embodiment in the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0045] The principles and features of the present invention are described below in conjunction with the accompanying drawings. The examples given are only used to explain the present invention and are not intended to limit the scope of the present invention. In the following paragraphs, the present invention is described more specifically by way of example with reference to the accompanying drawings. It should be noted that the drawings are all in a very simplified form and use non-precise scales, only for the purpose of facilitating and clearly assisting in explaining the objectives of the embodiments of the present invention.

[0046] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by those skilled in the technical field to which the present invention belongs. The terms used herein in the specification of the present invention are only for the purpose of describing specific embodiments and are not intended to limit the present invention. The term "and / or" used herein includes any and all combinations of one or more of the related listed items.

[0047] Embodiment 1

[0048] As Figure 1 shown in -9, a numerical simulation method for multiphase fluid formation seepage proposed by the present invention includes the following steps:

[0049] S1. Generate an unstructured mesh by combining the advancing front method and Delaunay triangulation;

[0050] S2. Dynamically adjust according to the changes in seepage physical quantities through adaptive mesh refinement technology;

[0051] S3. Establish a dynamic phase change model and obtain the phase change rate expressed as a function of temperature, pressure, and the concentration of each phase;

[0052] S4. Calculate the changes in the mass and energy of each phase during the phase change process in real time through the kinetic equation;

[0053] S5. Use the phase field method to handle the dynamic evolution of the phase interface and describe the distribution of different phases through phase field variables;

[0054] S6. Couple the grid system obtained by multi-scale unstructured grid division with the dynamic phase change model through a numerical solution algorithm;

[0055] S7. Discretize the mass conservation equation, momentum conservation equation, and energy conservation equation using the finite volume method.

[0056] The advancing front method starts from the formation boundary and gradually advances inward to generate triangular grid elements. During the advancement process, the Delaunay triangulation characteristics are used to ensure good quality of the generated grid elements;

[0057] In the advancing front method, new nodes are generated based on the optimization function

[0058]

[0059] to determine the position,

[0060] where θi is the interior angle of the triangular element;

[0061] Taking the pressure gradient ▽P as an example for the adaptive mesh refinement technology, when |▽P| > τ 1 subdivision is carried out. The subdivision method is to divide the triangular element from the centroid to the midpoint connection of the three sides into four sub-triangular elements. When ∣▽P∣ < τ2, adjacent element merging is carried out;

[0062] τ1 is the refinement threshold;

[0063] τ2 is the coarsening threshold;

[0064] And τ2 < τ1.

[0065] The adaptive mesh refinement technology automatically refines the mesh in the regions where the physical quantities such as pressure gradient and flow velocity change significantly, and uses coarser meshes in the regions where the physical quantities change gently. Through the adaptive mesh refinement technology, the computational effort is concentrated in the regions where the pressure gradient > τ1, such as near the wellbore, and the overall number of meshes is reduced by 40%. Combined with GPU parallel computing, efficient simulation of mesh elements is achieved, and the calculation time of a single node is shortened.

[0066] The distribution of different phases is described by the phase field variable, and the phase field equation is solved to accurately capture the movement and deformation of the phase interface during the seepage process.

[0067] The finite volume method integrates the control equation on the mesh elements and uses methods such as central difference to discretize the time and space derivatives.

[0068] The present invention generates unstructured meshes by combining the advancing front method and Delaunay triangulation. The advancing front method starts from the formation boundary and gradually advances inward to generate triangular mesh elements. During the advancing process, the characteristics of Delaunay triangulation are used to ensure that the generated mesh elements have good quality and avoid the appearance of distorted meshes. For example, in the simulation of a certain complex formation, the meshes generated by this method can closely fit the complex boundary of the formation, and the average angular deviation of the mesh elements is controlled within 5°. At the same time, according to the change of seepage physical quantities, an adaptive mesh refinement technology is introduced. This technology automatically refines the mesh in the regions where the physical quantities such as pressure gradient and flow velocity change significantly, such as near the wellbore and where the formation permeability changes abruptly, and uses coarser meshes in the regions where the physical quantities change gently. Taking the pressure gradient as an example, when the pressure gradient is greater than the set refinement threshold, the triangular element is divided from the centroid to the midpoint connection of the three sides into four sub-triangular elements for refinement; when the pressure gradient is less than the coarsening threshold, adjacent mesh elements are merged. Through this dynamic adjustment, while ensuring the overall computational efficiency, the simulation accuracy is significantly improved. In actual simulation, in the region near the wellbore, through adaptive refinement, the pressure calculation accuracy is improved by more than 30%

[0069] For the phase change between the gas phase g and the liquid phase l, the phase change rate

[0070] R gl = k1·f(P, T, C g , C l )

[0071] where k1 is the kinetic coefficient,

[0072] f(P, T, C g , C l )

[0073] is a function of pressure P, temperature T, and gas phase concentration Cg and the function of the liquid-phase concentration Cl;

[0074] f(P, T, C g , C l ) = (P - P sat (T))·C g ·C l , P sat (T)

[0075] is the saturation pressure at the corresponding temperature.

[0076] The dynamic phase transition model uses the phase field method to deal with the dynamic evolution of the phase interface, describes the distribution of different phases through phase field variables, and solves the phase field equation

[0077]

[0078] to accurately capture the movement and deformation of the phase interface during the seepage process, where M is the mobility and F is the free energy functional

[0079]

[0080] is the local free energy density, is the interface thickness parameter.

[0081] Based on the principles of thermodynamics and mass conservation, a dynamic phase transition model is established. The phase transition rate is expressed as a function of temperature, pressure, and the concentrations of each phase. For example, for the phase transition between the gas phase and the liquid phase, the phase transition rate is related to these physical quantities through a specific functional relationship. By calculating the changes in the mass and energy of each phase during the phase transition process in real time through the kinetic equation, the actual phase transition situation of multiphase fluids in formation seepage can be accurately reflected. At the same time, the phase field method is used to deal with the dynamic evolution of the phase interface. The distribution of different phases is described through phase field variables, and the phase field equation is solved to accurately capture the movement and deformation of the phase interface during the seepage process. The phase field method avoids the discontinuity problem in the treatment of the phase interface in traditional models, making the simulation of the phase interface more realistic and accurate. When simulating the mixing region of multiphase fluids, the phase field method can clearly show the dynamic changes of the phase interface, which is highly consistent with the actual observation results.

[0082] The solution process of the numerical solution algorithm includes initializing the calculation parameters, performing iterative solutions within each time step according to the governing equations and boundary conditions, discretizing the time derivative using the implicit finite volume method, and the time step Δt needs to satisfy Δt ≤ CFL·Δx / |u|, where CFL is the Courant number taking 0.5 - 1.0, Δx is the grid size, and u is the fluid velocity. For regions with a sharp pressure gradient such as near the wellbore, the grid is automatically refined and Δt is dynamically adjusted to ensure stability;

[0083] During the solution process, update the variables related to the grid system and the phase change model until the set simulation time or convergence condition is reached. The initialization calculation parameters include the initial pressure P0, initial temperature T0, fluid saturation S0 of the formation, as well as the relevant parameters of the grid system and the phase change model, such as the kinetic coefficient k1 and mobility M in the dynamic phase change model. In the data preparation stage, collect geological data such as the geometry, permeability distribution, and porosity of the formation, as well as physical property parameters of multiphase fluids, such as density, viscosity, and phase change heat.

[0084] Couple the grid system obtained by multi-scale unstructured grid division with the dynamic phase change model. Use the finite volume method to discretize the mass conservation equation, momentum conservation equation, and energy conservation equation. During the discretization process, integrate the control equations over the grid cells and use methods such as central difference to discretize the time and space derivatives. Through this coupling method, the seepage and phase change processes of multiphase fluids can be considered simultaneously, and the physical quantities of each phase at different times and positions can be obtained. The solution process of the numerical solution algorithm includes initializing the calculation parameters, which cover the initial pressure, temperature, fluid saturation of the formation, as well as the relevant parameters of the grid system and the phase change model. Then, perform iterative solutions within each time step according to the control equations and boundary conditions. During the solution process, update the grid system in real time to adapt to the changes in seepage physical quantities, and calculate the mass and energy exchange between each phase according to the dynamic phase change model. Finally, by continuously advancing the time step, obtain the dynamic simulation results of multiphase fluid formation seepage.

[0085] In this embodiment, the porosity range of the geological data collected from the formation is 0.15 - 0.25, and the permeability varies between 10 -15 -10 - 12 m 2 . At the same time, determine the physical property parameters of the oil, gas, and water three-phase fluids. The density of the oil phase is 850 kg / m 3 , and the viscosity is 5 mPa·s; the density of the gas phase is 12 kg / m 3 , and the viscosity is 0.012 mPa·s; the density of the water phase is 1000 kg / m 3 , and the viscosity is 1 mPa·s; the formation has a complex heterogeneous structure. Using the multi-scale unstructured grid division algorithm of the present invention, generate the initial grid according to the formation geometry. After the initial grid is generated, through the analysis of the initial distribution of seepage physical quantities, in the area near the oil well, due to the large pressure gradient, the number of grid cells in this area is increased by 50% through the adaptive refinement technology, effectively improving the simulation accuracy.

[0086] Initialize the dynamic phase change model according to the collected physical properties of the multiphase fluid and the initial formation conditions. Set the oil phase saturation to 0.2, the gas phase saturation to 0.7, and the water phase saturation to 0.1 at the initial moment. At the same time, determine the relevant parameters of the phase change model, with the kinetic coefficient being 10 -3 s -1 , and the mobility being 10 -10 m 2 / (N·s); Inject dry methane gas with an injection pressure of 42 MPa and an injection rate of 5000 m 3 / d.

[0087] Use the coupled solution algorithm to discretely solve the governing equations within each time step. After 1000 time steps of calculation with a simulation time of 100 days, the dynamic changes in the multiphase fluid seepage in the oil production area are obtained.

[0088] Analyze the results obtained from the simulation. The pressure distribution, flow velocity distribution, and saturation distribution of the multiphase fluid are basically in line with the actual monitoring data; in the simulation, the dry gas front advances to 50 meters from the injection well, the pressure is maintained at 38 MPa above the dew point, the condensate oil saturation increases from 0.2 to 0.4, and the retrograde condensation phenomenon is effectively suppressed; the advancing path of the injected gas in the reservoir is consistent with the actual observation situation, verifying the effectiveness of this numerical simulation method in the simulation of multiphase fluid seepage in oil production.

[0089] Multi-scale unstructured grid division can better fit the complex formation geometry, and the simulation accuracy for irregular regions is improved by more than 50% compared to traditional structured grids. The dynamic phase change model can track the changes in the phase change process over time and space in real-time. Compared with the traditional static phase change model, the simulation error of the phase change process is reduced by about 40%, more accurately reflecting the actual situation of multiphase fluid seepage in the formation. The adaptive grid refinement technology reduces unnecessary computational effort while ensuring the simulation accuracy. By using coarser grids in areas with gentle changes in physical quantities, the overall calculation time is shortened by 30%-40% compared to the full-region fine grid simulation, improving the computational efficiency.

[0090] As described above, it is only the preferred embodiment of the present invention and does not impose any form of limitation on the present invention; any ordinary technician in this industry can smoothly implement the present invention according to what is shown in the accompanying drawings of the specification and what is described above; however, any equivalent changes such as slight modifications, decorations, and evolutions made by those skilled in this profession without departing from the technical solution scope of the present invention by using the technical content disclosed above are all equivalent embodiments of the present invention; at the same time, any equivalent changes, modifications, and evolutions made to the above embodiments based on the essential technology of the present invention still fall within the protection scope of the technical solution of the present invention.

Claims

1. A numerical simulation method for multiphase fluid formation seepage, characterized in that, The method comprises the following steps: S1, the unstructured grid is generated by combining the frontier advancing method with Delaunay triangulation; S2, dynamically adjust according to the changes of seepage physical quantities through adaptive grid encryption technology; S3, establish a dynamic phase change model to obtain the phase change rate as a function of temperature, pressure and concentration of each phase; S4, calculate the changes of mass and energy of each phase in the phase change process in real time through kinetic equations; S5. Use the phase field method to deal with the dynamic evolution of the phase interface and describe the distribution of different phases through phase field variables; S6, coupling the grid system obtained by multi-scale unstructured grid division with the dynamic phase change model through a numerical solution algorithm; S7. Use the finite volume method to discretize the mass conservation equation, momentum conservation equation, and energy conservation equation.

2. The numerical simulation method for multiphase fluid formation seepage according to claim 1, characterized in that The frontier advancing method starts from the stratum boundary and gradually advances to the inside to generate triangular grid units. In the advancement process, the Delaunay triangulation characteristics are used to ensure the quality of the generated grid units. In the frontier advancing method, new nodes are generated based on the optimization function Determine the location, Where θi is the internal angle of the triangle unit; The adaptive mesh encryption technology takes the pressure gradient as an example. When happens, subdivision is carried out. The subdivision method is to divide the triangular element from the centroid to the midpoint connection of the three sides into four sub-triangular elements. When happens, adjacent element merging is carried out; τ1 is the encryption threshold; τ2 is the coarsening threshold; And τ2<τ1.

3. A numerical simulation method for multiphase fluid formation seepage according to claim 2, characterized in that, Adaptive mesh encryption technology automatically encrypts the mesh in areas where the physical quantities change dramatically, and uses a coarser mesh in areas where the physical quantities change gently, according to the changes in the seepage physical quantities such as pressure gradient and flow velocity.

4. A numerical simulation method for multiphase fluid formation seepage according to claim 3, characterized in that The distribution of different phases is described by phase field variables, and the phase field equations are solved to accurately capture the movement and deformation of the phase interface during the seepage process.

5. The numerical simulation method for multiphase fluid formation seepage according to claim 4, wherein The finite volume method integrates the control equations on the grid cells and discretizes the time and space derivatives using methods such as central difference.

6. The numerical simulation method for multiphase fluid formation seepage according to claim 5, wherein, For the phase transition between gas phase g and liquid phase l, the phase transition rate R gl = k1·f(P, T, C g , C l ) Where k1 is the kinetic coefficient, f(P, T, C g , C l ) It is a function of pressure P, temperature T, gas-phase concentration C g and liquid-phase concentration C l ; f(P, T, C g , C l ) = (P - P sat (T)) · C g · C l , P sat (T) is the saturation pressure at the corresponding temperature.

7. A numerical simulation method for multiphase fluid formation seepage according to claim 6, characterized in that, The dynamic phase change model uses the phase field method to deal with the dynamic evolution of the phase interface, describes the distribution of different phases through phase field variables, and solves the phase field equation To accurately capture the movement and deformation of the phase interface during the percolation process, where M is the mobility and F is the free energy functional is the local free energy density, is the interface thickness parameter.

8. A numerical simulation method for multiphase fluid formation seepage according to claim 7, characterized in that The solution process of the numerical solution algorithm includes initializing the calculation parameters, performing iterative solutions in each time step according to the control equations and boundary conditions, and updating the grid system and phase change model related variables during the solution process until the set simulation time or convergence condition is reached.

9. A numerical simulation method for multiphase fluid formation seepage according to claim 8, characterized in that, The initialization calculation parameters include the initial pressure P0, initial temperature T0, fluid saturation S0 of the formation, and related parameters of the grid system and phase change model, such as the kinetic coefficient k1 and mobility M in the dynamic phase change model.

10. A numerical simulation method for multiphase fluid formation seepage according to claim 9, characterized in that, In the data preparation stage, geological data such as the geometry of the formation, permeability distribution, porosity, and physical property parameters such as density, viscosity, phase change heat, etc. of the multiphase fluid are collected.

Images