A high-efficiency and non-oscillatory method for predicting two-phase flow patterns in marine oil and gas pipelines

Abstract

The present application relates to the technical field of computational fluid dynamics, in particular to a high-efficiency and non-oscillating two-phase form prediction method for ocean oil and gas pipes, which comprises setting oil and gas physical property parameters, setting time step and boundary conditions, initializing storage matrix of velocity and pressure; constructing coefficient matrix and performing preprocessing, judging whether the current time step is less than the set time step, if not, outputting the current oil and gas form distribution; otherwise, solving Cahn-Hilliard equation according to the applied boundary conditions and velocity field, reconstructing to obtain flow field physical property parameters; discretizing Navier-Stokes equation by finite volume method and constant coefficient projection algorithm and performing constantization, taking the flow field physical property parameters as the input of the constantized Navier-Stokes equation to calculate the flow field velocity distribution and pressure distribution; the present application effectively improves the solving efficiency, and realizes non-pressure oscillation calculation on the same position grid through the correction of velocity.

Description

Technical Field

[0001] This invention relates to the field of computational fluid dynamics, specifically to the field of numerical simulation in marine engineering, and particularly to a highly efficient and oscillatory method for predicting the two-phase morphology of marine oil and gas pipelines. Background Technology

[0002] Marine oil and gas resources are widely distributed, covering most of the Earth's oceans. It is estimated that over 90% of the world's undeveloped oil and gas resources are located in the ocean. Global marine oil and gas reserves account for more than two-thirds of global oil reserves, most of which remain undeveloped. Marine oil and gas pipelines, as vital channels for transporting these resources, play a crucial role in energy supply. The flow characteristics and interactions of the gas-liquid two-phase systems within these pipelines significantly impact their stability and transport efficiency.

[0003] Computational Fluid Dynamics (CFD) is a method for studying fluid motion based on computer numerical simulation. It discretizes a continuous fluid space into a grid system and uses numerical methods to discretize the governing equations of fluid dynamics, generating a set of algebraic equations. Solving these algebraic equations allows for the numerical simulation of fluid mechanics problems. In marine oil and gas pipeline flows, CFD can be used to predict the flow state and interactions of multiphase fluids, as well as quantify the magnitude and distribution of physical quantities such as velocity, pressure, density, and concentration.

[0004] In the field of computational fluid dynamics, decoupling pressure and velocity is also a key challenge in the simulation of incompressible flows. To date, algorithms for decoupling and solving incompressible Navier-Stokes flow can generally be divided into two categories: iterative methods (SIMPLE, PISO, etc.) and non-iterative methods (projection methods, etc.). However, all of these methods inevitably face the difficulty of suppressing pressure oscillations.

[0005] Marine oil and gas pipelines contain numerous gas-liquid two-phase interfaces, which are constantly in a delicate balance between surface tension, gravity, and pressure gradients. Therefore, even minute numerical oscillations can affect the accuracy of numerical simulations. Particularly in long-term simulations, small pressure disturbances can evolve into sources or sinks, driving non-physical movements of the fluid interface. This can introduce significant errors in the distribution of oil and gas morphology and flow rate measurements, leading to design deviations in the pipeline and even safety hazards.

[0006] Therefore, how to achieve pressureless oscillation calculation on the same grid is an urgent problem to be solved in this field. Summary of the Invention

[0007] In view of this, this invention proposes an efficient and oscillatory-free method for predicting the two-phase morphology of marine oil and gas pipelines. A two-dimensional, equidistant, orthogonal Cartesian grid is used to divide the computational flow domain. After division, the physical parameters of the flow field are discretely distributed at the center of each grid. By solving the partial differential equations of the flow and phase interface, the dynamic changes of the physical parameters on the grid are obtained. Based on these dynamic changes, the two-phase morphology is predicted. The solution for the dynamic changes of the physical parameters specifically includes the following steps:

[0008] Set the oil and gas physical property parameters, set the time step and boundary conditions, and initialize the storage matrix for velocity and pressure;

[0009] Construct a coefficient matrix and perform preprocessing; determine whether the current time step is less than the set time step; if it is not less than the set time step, output the current oil and gas morphology distribution.

[0010] Otherwise, based on the applied boundary conditions and velocity field, the Cahn-Hilliard equations are used to capture the multiphase flow interface, that is, to solve the Cahn-Hilliard equations based on the set boundary conditions and reconstruct the flow field properties.

[0011] The Navier-Stokes equations are discretized using the finite volume method and constant coefficient projection algorithm. The coefficients of the discretized Navier-Stokes equations are then steadyened, and the flow field properties are substituted into the steadyened Navier-Stokes equations to calculate the flow field velocity and pressure distributions, thus completing the morphological prediction.

[0012] Furthermore, the process of reconstructing the flow field properties includes:

[0013]

[0014]

[0015] Where ρ is the density of the fluid particles; ρ1 is the density of the pure gas phase; ρ2 is the density of the pure liquid phase; φ is the phase fraction variable; μ is the viscosity of the fluid particles; μ1 is the viscosity of the pure gas phase; and μ2 is the viscosity of the pure liquid phase.

[0016] Furthermore, the phase fractional variables are obtained through the Cahn-Hilliard equations, which include:

[0017]

[0018]

[0019] Where φ is the phase fraction variable; t represents time, and the partial derivative of time t is... Or the derivative dt represents the time step; This indicates the calculation of divergence; This indicates calculating the gradient of the parameters on the right-hand side; denoted by ; M represents the velocity vector of the fluid particles; ψ represents the mobility; ε represents the chemical potential; and ε represents the thickness of the phase interface.

[0020] Furthermore, when solving the Cahn-Hilliard equation, the boundary conditions for both the phase variable and the chemical potential are Neumann boundaries.

[0021] Furthermore, the process of discretizing the Navier-Stokes equations using the finite volume method and the constant coefficient projection algorithm, and then stabilizing the coefficients of the discretized Navier-Stokes equations, includes:

[0022] In the prediction step, transient terms in the prediction step equation are discretized using a second-order backward difference scheme and explicitly processed using a second-order time extrapolation scheme. To obtain a steady coefficient matrix, a constant-coefficient dynamic viscosity term is added to both sides of the prediction step equation, and an approximation is introduced to decompose the unsteady pressure gradient into steady and unsteady parts. Numerical methods are used to discretize the partial differential equation of the prediction step, resulting in the prediction step equation AU = b with a steady coefficient matrix. u and AV=b v Where A is the first coefficient matrix, denoted as A = {a P ,a E ,a W ,a N ,a S}, a P The first coefficient value of the central grid is represented as h is the height of the Cartesian grid, v0 represents the structural viscosity of the fluid particles, denoted as... ρ1 is the density of the pure gas phase, ρ2 is the density of the pure liquid phase, μ1 is the viscosity of the pure gas phase, μ2 is the viscosity of the pure liquid phase, and dt is the time step of each cycle; a E The first coefficient value of the grid adjacent to the east side of the central grid is denoted as -v0dt; a W The first coefficient value of the grid adjacent to the west side of the central grid is denoted as -v0dt; a N The first coefficient value of the grid adjacent to the north side of the central grid, denoted as -v0dt; a S The first coefficient value of the grid adjacent to the south side of the central grid is denoted as -v0dt; U is the horizontal component of the flow field velocity to be solved; V is the vertical component of the flow field velocity to be solved; b u b is the right-hand side of the equation to be solved, representing the prediction step for a steady coefficient matrix when the velocity is in the horizontal direction. v Let be the right-hand side term of the equation to be solved, which is the prediction step for when the coefficient matrix is ​​constant when the velocity is in the horizontal direction.

[0023] In the correction step, the correction step equation is discretized using a second-order backward difference scheme. An approximation is introduced into the second-order time extrapolation scheme, and an incompressibility condition is included. Then, the divergence of both sides of the correction step equation is taken to obtain the unsolved equation for the pressure correction step, namely the pressure Poisson equation with a constant coefficient matrix. The pressure Poisson equation with a constant coefficient matrix is ​​discretized using numerical methods to obtain the discretized equation, expressed as LP = b. p Where L is the second coefficient matrix, expressed as L={L P ,L E ,L W ,L N ,L S},L P The second coefficient value of the central grid is expressed as L E The second coefficient value of the grid adjacent to the east side of the central grid is expressed as: L W The second coefficient value of the grid adjacent to the west of the central grid is expressed as: L N The second coefficient value of the grid adjacent to the north side of the central grid is expressed as: L S The second coefficient value of the grid adjacent to the south side of the central grid is expressed as: P represents the pressure distribution in the flow field to be solved; b p Represented as u * To solve for the horizontal component of the obtained flow field velocity, v * To solve for the components of the obtained flow field velocity in the vertical direction, This indicates taking the difference in the horizontal direction. This indicates that the difference is calculated with respect to the direction of the plumb bob.

[0024] Furthermore, after obtaining the velocity distribution U in the flow field, the velocity is corrected. The correction process includes:

[0025]

[0026] in, The horizontal velocity component of the central grid obtained from the prediction step equation; The horizontal velocity component of the first grid east of the central grid obtained from the prediction step equation; h is the Cartesian grid height; P W The pressure is the pressure of the first grid located west of the central grid; P P Pressure on the central grid; P E The pressure is the pressure of the first grid located east of the central grid; P EEdenoted as dt, where α is the pressure of the second grid located east of the central grid; α is a correction factor, and α = dt.

[0027] Furthermore, the equation to be solved for the prediction step of the coefficient matrix with time invariance is expressed as:

[0028]

[0029] in, ρ represents the fluid velocity vector at time step *, where * indicates the intermediate step between time step n and time step (n+1); n+1 P represents the density of the fluid particle at time step n+1; *,n+1 The structural pressure of the fluid particle n+1 is denoted as P. *,n+1 =2P n -P n-1 ; Represents the discrete computational form of the diffusion term; μ n+1 This represents the viscosity of a fluid particle at time step n+1; The construction velocity of a fluid particle is represented as... This represents the resultant force of gravity and surface tension. This represents the discrete calculation form of the convection term; ρ is the density of the fluid particles. Let be the velocity vector of the fluid particle at time step n; ρ0 is the structure density of the fluid particle, expressed as ρ0 = min(ρ1, ρ2); P n This represents the pressure of a fluid particle at time step n. This represents the Laplace operator.

[0030] Furthermore, Approximately The pressure Poisson equation with a constant coefficient matrix is ​​expressed as:

[0031]

[0032] P' = P n+1 -P n

[0033] Among them, U * Let be the velocity vector of the fluid particle at time *.

[0034] This invention proposes a constant coefficient matrix projection algorithm for solving two-phase flows on a co-located grid without oscillation based on the phase-field method. This effectively avoids the difficulty of reconstructing the coefficient matrix in each step of the semi-implicit solution of the Navier-Stokes equations. At the same time, the coefficient matrix is ​​preprocessed by combining LU decomposition, which can effectively improve the solution efficiency. In addition, by simply modifying the velocity, pressureless oscillation calculation on the co-located grid is realized, which can quickly and accurately predict the bubble morphology in oil and gas pipelines. Attached Figure Description

[0035] Figure 1 This is a flowchart of a highly efficient and oscillation-free method for predicting the two-phase morphology of marine oil and gas pipelines according to the present invention.

[0036] Figure 2 The image shows the rising morphology of bubbles in the oil and gas pipeline predicted using the algorithm of this invention when the oil-gas density ratio is 10.

[0037] Figure 3 This is a diagram showing the rising morphology of bubbles in an oil and gas pipeline, predicted using the algorithm of this invention when the oil-gas density ratio is 1000. Detailed Implementation

[0038] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0039] This invention proposes an efficient and oscillatory-free method for predicting the two-phase morphology of marine oil and gas pipelines. It employs a two-dimensional, equidistant, orthogonal Cartesian grid to divide the computational flow domain. After division, the physical parameters of the flow field are discretely distributed at the center of each grid. By solving the partial differential equations of the flow and phase interface, the dynamic changes of the physical parameters on the grid are obtained. Based on these dynamic changes, the two-phase morphology is predicted. The specific steps for solving the dynamic changes of the physical parameters include:

[0040] Set the oil and gas physical property parameters, set the time step and boundary conditions, and initialize the storage matrix for velocity and pressure;

[0041] Construct a coefficient matrix and perform preprocessing; determine whether the current time step is less than the set time step; if it is not less than the set time step, output the current oil and gas morphology distribution.

[0042] Otherwise, based on the applied boundary conditions and velocity field, the Cahn-Hilliard equation is used to capture the multiphase flow interface, and the Cahn-Hilliard equation is solved according to the set boundary conditions to reconstruct the flow field properties.

[0043] The Navier-Stokes equations are discretized using the finite volume method and constant coefficient projection algorithm. The coefficients of the discretized Navier-Stokes equations are then steadyened, and the flow field properties are substituted into the steadyened Navier-Stokes equations to calculate the flow field velocity and pressure distributions, thus completing the morphological prediction.

[0044] In this embodiment, the marine oil and gas pipeline flow is physically modeled, meshed, numerically discretized, and solved sequentially to predict oil and gas distribution. Specifically, this includes initializing the two-phase interface, setting the time step, meshing, and initializing the phase interface. Then, with each time step, the Cahn-Hilliard equation and the Navier-Stokes equation are solved iteratively. The algorithm flow is as follows: Figure 1 As shown, this embodiment describes the present invention from three aspects: multiphase flow interface capture algorithm, incompressible multiphase flow constant coefficient matrix projection algorithm, and pressureless oscillation correction technology.

[0045] (I) Multiphase Flow Interface Capture Algorithm

[0046] The phase-field method is chosen as the interface capture algorithm, and the motion of the two-phase interface is described by the Cahn-Hilliard equations:

[0047]

[0048]

[0049] Where φ is a phase fraction variable, ranging from [-1, 1], where 1 represents one phase and -1 represents another phase. Fractions greater than -1 and less than 1 represent the interface between two phases. In this embodiment, the two phases are the gas phase and the liquid phase, respectively; t represents the time step. This indicates the calculation of divergence. This indicates calculating the gradient of the parameters on the right side. This indicates the discretization of the right-hand Laplace operator, i.e. denoted by ; M represents the velocity vector of the fluid particles; M is the mobility that defines the diffusivity of the phase interface region, and in this embodiment, a fixed value M = 0.01 is used; ψ is the chemical potential, derived from the free energy functional; and ε represents the thickness of the phase interface.

[0050] (II) Algorithm for Projecting Constant Coefficient Matrix of Incompressible Multiphase Flow

[0051] In this embodiment, based on the incremental pressure correction algorithm, the momentum equation is solved in two parts sequentially:

[0052]

[0053] in, The resultant force of gravity and surface tension is represented by ρ, where ρ is the density of the fluid particle, and is expressed as... In this embodiment, ρ1 is the density of the pure gas phase and ρ2 is the density of the pure liquid phase; This indicates taking the partial derivative with respect to the parameters on the right. Let momentum flux be expressed as The Kronecker product is represented by μ, where μ is the viscosity of the fluid particle, expressed as... In this embodiment, μ1 is the pure gas phase viscosity and μ2 is the pure liquid phase viscosity; * indicates a virtual time step located between time step n and time step n+1.

[0054] The process of discretizing the Navier-Stokes equations using the finite volume method and constant coefficient projection algorithm, and then stabilizing the coefficients of the discretized Navier-Stokes equations, includes:

[0055] In the prediction step, a second-order backward difference scheme is used, i.e. The format is: Taking the transient term in the discrete prediction step equation and shifting the other terms to the right, we get:

[0056]

[0057] Let N and D represent the convection and viscous terms respectively, that is, let Explicit processing is performed using a second-order time extrapolation format, and the explicit processing is represented as follows:

[0058] χ *,n+1 =2χ n -χ n-1

[0059] get:

[0060]

[0061] To ensure the stability of the numerical solution, a semi-implicit scheme is used to handle the viscosity term. However, to obtain a steady coefficient matrix, this embodiment adds a constant-coefficient dynamic viscosity term to both sides of the equation. Here, the dynamic viscosity v0 is taken as:

[0062]

[0063] Then, the equation with solution becomes:

[0064]

[0065] At this point, an approximation is introduced to decompose the unsteady pressure gradient into two parts: steady and unsteady.

[0066]

[0067] The structural density ρ0 of the fluid particles is expressed as ρ0 = min(ρ1, ρ2).

[0068] The final equation to be solved for the prediction step with a constant coefficient matrix is:

[0069]

[0070] The partial differential equations of the prediction step are discretized using numerical methods. In this embodiment, discretization is performed using an equidistant orthogonal structured grid, i.e., dx = dy = h, to obtain the prediction step equations with a constant coefficient matrix:

[0071] AU = b u

[0072] AV=b v

[0073] Where A is the first coefficient matrix, denoted as A = {a P ,a E ,a W ,a N ,a S}, a P The first coefficient value of the central grid is represented as h is the height of the Cartesian grid, v0 represents the structural viscosity of the fluid particles, denoted as... ρ1 is the density of the pure gas phase, ρ2 is the density of the pure liquid phase, μ1 is the viscosity of the pure gas phase, μ2 is the viscosity of the pure liquid phase, and dt is the time step of each cycle; a E The first coefficient value of the grid adjacent to the east side of the central grid is denoted as -v0dt; a W The first coefficient value of the grid adjacent to the west side of the central grid is denoted as -v0dt; a N The first coefficient value of the grid adjacent to the north side of the central grid, denoted as -v0dt; a S The first coefficient value of the grid adjacent to the south side of the central grid is denoted as -v0dt; U is the horizontal component of the flow field velocity to be solved; V is the vertical component of the flow field velocity to be solved; b u b v This is a vector consisting of transient, convection, viscous, and source terms. These two parameters are constructed using parameters obtained from solving the Cahn-Hilliard equations. u b is the right-hand side of the equation to be solved, representing the prediction step for a steady coefficient matrix when the velocity is in the horizontal direction. v Let be the right-hand side of the equation to be solved for the prediction step when the velocity is horizontal. That is, under horizontal and vertical velocities, the right-hand side of the equation to be solved for the prediction step with a constant coefficient matrix is ​​finally obtained as a vector consisting of transient, convection, viscous, and source terms.

[0074] Discretizing the correction step equation using a second-order backward difference scheme yields:

[0075]

[0076] Right now:

[0077]

[0078] From the above formula, we can obtain:

[0079]

[0080] At this point, we introduce and approximate, namely:

[0081]

[0082] Then we can get:

[0083]

[0084] P' = P n+1 -P n

[0085] Introducing the incompressibility condition (the irrotational velocity formula), i.e. Taking the divergence at both ends simultaneously, we obtain the equation to be solved for the pressure correction step, which is the pressure Poisson equation with a constant coefficient matrix:

[0086]

[0087] Among them, U * Let be the velocity vector of the fluid particle at time *.

[0088] The pressure Poisson equation with a constant coefficient matrix is ​​discretized using numerical methods to obtain the discretized equation, which is expressed as:

[0089] LP = b p

[0090] Where L is the second coefficient matrix, denoted as L={L P ,L E ,L W ,L N ,L S},L P The second coefficient value of the central grid is expressed as L E The second coefficient value of the grid adjacent to the east side of the central grid is expressed as: L W The second coefficient value of the grid adjacent to the west of the central grid is expressed as: L N The second coefficient value of the grid adjacent to the north side of the central grid is expressed as: L S The second coefficient value of the grid adjacent to the south side of the central grid is expressed as: P represents the pressure distribution in the flow field to be solved; b p The right-hand side of the pressure Poisson equation with a constant coefficient matrix is ​​expressed as follows: u * To solve for the horizontal component of the obtained flow field velocity, v * To solve for the components of the obtained flow field velocity in the vertical direction, This indicates taking the difference in the horizontal direction. This indicates that the difference is calculated with respect to the direction of the plumb bob.

[0091] Finally, the fluid velocity vector at time step n+1 can be obtained, expressed as:

[0092]

[0093] In summary, the coefficient matrices of the equations to be solved in both the prediction and pressure correction steps of the constant coefficient matrix projection algorithm are no longer related to the dynamically changing density and viscosity of the flow field. The coefficient matrix can be constructed after the program has completed mesh generation and material parameter settings, but before the time-step solution begins. This avoids the difficulty of reconstructing the coefficient matrix in each iteration of general projection algorithms. Furthermore, combining this with LU decomposition for preprocessing significantly accelerates computational efficiency.

[0094] (III) Pressure-free oscillation correction technique

[0095] This invention employs a simpler data storage structure called a co-located grid, but directly using a co-located grid inevitably leads to pressure oscillations. Based on the idea of ​​momentum interpolation, this invention achieves tight coupling between pressure and velocity through modification, realizing oscillation-free numerical calculation of two-phase flow in the phase field on a co-located grid.

[0096] Within the framework of the constant coefficient matrix projection algorithm, a simpler and more effective pressure oscillation suppression scheme is achieved. This is achieved after completing the prediction step solution. Then, the speed was corrected as follows:

[0097]

[0098] in, The horizontal velocity component of the central grid obtained from the prediction step equation; The horizontal velocity component of the first grid east of the central grid obtained from the prediction step equation; h is the Cartesian grid height; P W The pressure is the pressure of the first grid located west of the central grid; P P Pressure on the central grid; P E The pressure is the pressure of the first grid located east of the central grid; P EE The pressure is the pressure of the second grid located east of the central grid; α is a correction factor. In the variable coefficient matrix algorithm, α is related to various factors such as density, viscosity, grid, and time step. However, in the constant coefficient matrix of this invention, α = dt. In this invention, the effect of pressure-free oscillation can be achieved through simple correction.

[0099] The constant coefficient matrix non-oscillation algorithm proposed in this invention also has the ability to simulate the rise of bubbles with high density ratios, showing good agreement. To compare the efficiency improvement of this algorithm, this embodiment compares it with a general coefficient matrix projection algorithm based on three different grid resolutions, using the same hardware conditions (AMD Ryzen 7 4800H, Radeon Graphics 2.90GHz) and settings. Table 1 lists the CPU time under the three grid resolutions. Because it completely avoids matrix construction for each time step and pre-decomposes LU, the constant coefficient matrix algorithm proposed in this paper is far more efficient than the variable coefficient matrix algorithm while ensuring computational accuracy. Due to the small density ratio in case 1 and the simple computation, the acceleration effect is more obvious; while case 2 has a large density ratio and complex topology changes, and the acceleration effect is not as good as case 1, but it is still more than twice as fast, which also greatly saves computational resources and improves computational efficiency.

[0100] Table 1 Comparison of computational efficiency between the algorithm presented in this paper and the variable coefficient matrix algorithm.

[0101]

[0102] The efficiency and accuracy of the algorithm were tested by simulating the rising of bubbles inside an offshore oil and gas pipeline. Figure 2 and Figure 3 The initial interface state is as follows: the computational domain is a fully enclosed cavity with a length of [1m × 2m]. The initial position of the bubble is (x, y) = (0.5m, 0.5m), and the initial diameter is D = 0.5m. The top and bottom are non-slip solid walls, while the left and right sides are free-slip solid walls. Gravity points towards the bottom of the defined domain. Since the density of fluid 1 (bubble) is less than that of the surrounding fluid 2 (liquid), the bubble will slowly rise under the action of buoyancy. In order to satisfy the CFL condition and an appropriate interface width, the computational settings are slightly different at different resolutions. When the mesh width is h = 1 / 64, h = 1 / 128, and h = 1 / 256, the time steps are 1e-4, 5e-5, and 2.5e-5, respectively, and the interface widths are 0.02m, 0.01m, and 0.005m, respectively. The physical property parameters of Example 1 and Example 2 are shown in Table 2.

[0103] Table 2. Property parameter settings for the bubble rise example.

[0104]

[0105] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A highly efficient and oscillation-free method for predicting the two-phase morphology of marine oil and gas pipelines, characterized in that, The computational flow domain is meshed using a two-dimensional equidistant orthogonal Cartesian grid. After meshing, the physical parameters of the flow field are discretely distributed at the center of each grid. By solving the partial differential equations of the flow and phase interface, the dynamic changes of the physical parameters on the grid are obtained. Based on the dynamic changes of the physical parameters, two-phase morphology predictions are obtained. The solution for the dynamic changes of physical parameters includes the following steps: Set the oil and gas physical property parameters, set the time step and boundary conditions, and initialize the storage matrix for velocity and pressure; Construct a coefficient matrix and perform preprocessing; determine whether the current time step is less than the set time step; if it is not less than the set time step, output the current oil and gas morphology distribution. Otherwise, based on the applied boundary conditions and velocity field, the Cahn-Hilliard equations are used to capture the multiphase flow interface, that is, to solve the Cahn-Hilliard equations based on the set boundary conditions and reconstruct the flow field properties. The Navier-Stokes equations are discretized using the finite volume method and the constant coefficient projection algorithm, and the coefficients of the discretized Navier-Stokes equations are then steadied. The specific steps include: In the prediction step, transient terms in the prediction step equation are discretized using a second-order backward difference scheme and explicitly processed using a second-order time extrapolation scheme. To obtain a steady coefficient matrix, a constant-coefficient dynamic viscosity term is added to both sides of the prediction step equation, and an approximation is introduced to decompose the unsteady pressure gradient into steady and unsteady parts. Numerical methods are then used to discretize the partial differential equation of the prediction step, yielding the prediction step equation AU=b with a steady coefficient matrix. u and AV=b v Where A is the first coefficient matrix, denoted as A={a P , a E , a W , a N , a S }, a P The first coefficient value of the central grid is represented as h is the height of the Cartesian grid. The structural viscosity of a fluid particle is represented as . , The density is that of the pure gas phase. The density is that of the pure liquid phase. For pure gas phase viscosity, Let be the viscosity of the pure liquid phase, and dt be the time step of each cycle; a E The first coefficient value of the grid adjacent to the east side of the central grid is represented as: ;a W The first coefficient value of the grid adjacent to the west side of the central grid is represented as: ;a N The first coefficient value of the grid adjacent to the north side of the central grid is represented as: ;a S The first coefficient value of the grid adjacent to the south side of the central grid is represented as: U represents the horizontal component of the flow field velocity to be solved; V represents the vertical component of the flow field velocity to be solved; b u b is the right-hand side of the equation to be solved, representing the prediction step for a steady coefficient matrix when the velocity is in the horizontal direction. v Let be the right-hand side term of the equation to be solved, which is the prediction step for when the coefficient matrix is ​​constant when the velocity is in the horizontal direction. In the correction step, the correction step equation is discretized using a second-order backward difference scheme. An approximation is introduced into the second-order time extrapolation scheme, and an incompressibility condition is included. Then, the divergence of both sides of the correction step equation is taken to obtain the unsolved equation for the pressure correction step, namely the pressure Poisson equation with a constant coefficient matrix. The pressure Poisson equation with a constant coefficient matrix is ​​discretized using numerical methods to obtain the discretized equation, expressed as LP=b. p Where L is the second coefficient matrix, expressed as L={L P , L E ,L W , L N , L S }, L P The second coefficient value of the central grid is expressed as L E The second coefficient value of the grid adjacent to the east side of the central grid is expressed as: L W The second coefficient value of the grid adjacent to the west of the central grid is expressed as: L N The second coefficient value of the grid adjacent to the north side of the central grid is expressed as: L S The second coefficient value of the grid adjacent to the south side of the central grid is expressed as: P represents the pressure distribution in the flow field to be solved; b p Represented as u * To solve for the horizontal component of the obtained flow field velocity, v * To solve for the components of the obtained flow field velocity in the vertical direction, This indicates taking the difference in the horizontal direction. This indicates taking the difference with respect to the direction of the plumb bob; The flow field properties were then substituted into the steady-state Navier-Stokes equations to calculate the velocity and pressure distributions, thus completing the morphological prediction.

2. The efficient and oscillation-free method for predicting the two-phase morphology of marine oil and gas pipelines according to claim 1, characterized in that, The process of reconstructing the flow field properties includes: in, The density of the fluid particles; This is the density of the pure gas phase; The density is that of the pure liquid phase. For phase fractional variables; The viscosity of the fluid particles; Viscosity of pure gas phase; This is the viscosity of the pure liquid phase.

3. The efficient and oscillation-free method for predicting the two-phase morphology of marine oil and gas pipelines according to claim 2, characterized in that, The phase fractional variable is obtained by calculating using the Cahn-Hilliard equation, which includes: in, The phase fraction variable; t represents time; This indicates the calculation of divergence; This indicates calculating the gradient of the parameters on the right-hand side; The velocity vector represents the velocity of a fluid particle; M is the mobility. Chemical potential; This indicates the thickness of the phase interface.

4. The efficient and oscillation-free method for predicting the two-phase morphology of marine oil and gas pipelines according to claim 3, characterized in that, When solving the Cahn-Hilliard equation, the boundary conditions for both the phase variable and the chemical potential are Neumann boundaries.

5. The efficient and oscillation-free method for predicting the two-phase morphology of marine oil and gas pipelines according to claim 1, characterized in that, After obtaining the velocity distribution U in the flow field, the velocity is corrected. The correction process includes: ； in, The horizontal velocity component of the central grid obtained from the prediction step equation; The horizontal velocity component of the first grid east of the central grid obtained from the prediction step equation; h is the Cartesian grid height; P W The pressure is the pressure of the first grid located west of the central grid; P P Pressure on the central grid; P E The pressure is the pressure of the first grid located east of the central grid; P EE The pressure is the pressure of the second grid located to the east of the central grid; As a correction factor, and .

6. The efficient and oscillation-free method for predicting the two-phase morphology of marine oil and gas pipelines according to claim 1, characterized in that, The equation to be solved for the prediction step of the coefficient matrix with invariance is expressed as: in, This represents the fluid velocity vector at time step *, where * indicates the intermediate step between time step n and time step (n+1). This represents the density of the fluid particles at time step n+1; The structural pressure of the fluid particle n+1 is expressed as: ; Represents the discrete computational form of the diffusion term; This represents the viscosity of a fluid particle at time step n+1; The construction velocity of a fluid particle is represented as... ; This represents the resultant force of gravity and surface tension. This represents the discrete calculation form of the convection term; ρ is the density of the fluid particles. Let n be the velocity vector of the fluid particle at time step n; The structure density of the fluid particles is expressed as... ; This represents the pressure of a fluid particle at time step n. This represents the Laplace operator.

7. The efficient and oscillation-free method for predicting the two-phase morphology of marine oil and gas pipelines according to claim 6, characterized in that, Will Approximately The pressure Poisson equation with a constant coefficient matrix is ​​expressed as: in, Let be the velocity vector of the fluid particle at time *.

Images