Pipeline multiphase flow simulation method supporting slug flow

Abstract

The invention discloses a pipeline multiphase flow simulation method supporting slug flow, and relates to the field of pipeline multiphase flow simulation. The method comprises the steps of geometric information input and grid generation, boundary and initial condition setting, physical property parameter calculation and interpolation table construction, control equation discretization and solution, and time step size control and convergence judgment. The core improvement lies in that sectional staggered grid modeling is adopted to adapt to a complex geometric pipeline; a volume fraction self-adaptive correction mechanism is introduced to optimize an interphase drag force model, and the drag force underestimation problem under the low volume fraction is solved; the volume fraction is corrected through a remapping strategy, and the conservation is ensured. The method gives consideration to calculation efficiency and prediction precision, avoids numerical oscillation and divergence, is suitable for gas-liquid slug flow transient simulation of long-distance and complex geometric pipelines, and can provide support for pipeline system optimization and risk early warning.

Description

Technical Field

[0001] This invention relates to the field of pipeline multiphase flow simulation technology, and in particular to a pipeline multiphase flow simulation method that supports slug flow, applicable to transient simulation of gas-liquid slug flow in long-distance, complex geometric pipelines in industries such as oil extraction and chemical engineering. Background Technology

[0002] Pipeline transportation of multiphase fluids is widely used in oil extraction, chemical industry, and other industrial sectors. In long-distance or steeply inclined pipelines, slug flow is a typical unsteady flow pattern, characterized by the alternating formation of high-speed liquid plugs and bubble clusters between gas and liquid phases. This leads to drastic fluctuations in pressure and flow rate within the pipeline, reducing transportation efficiency and easily causing equipment vibration, accelerated corrosion, and even safety accidents. Therefore, accurate prediction and effective control of slug flow are crucial.

[0003] Existing simulation techniques for transient multiphase flow in pipelines are mainly divided into three-dimensional simulation and one-dimensional simulation. Three-dimensional simulation can capture the dynamics of the gas-liquid interface and local flow details with high precision, but it consumes large computational resources and is time-consuming, making it difficult to apply to full-scale simulation of long-distance industrial pipelines. One-dimensional simulation, based on simplified frameworks such as two-fluid models, reduces computational complexity through cross-sectional averaging, making it suitable for long-distance pipeline engineering applications, but it has significant drawbacks: First, the conservation of phase volume fractions is difficult to guarantee, and the sum of the volume fractions of the gas and liquid phases is prone to deviating from 1, causing numerical oscillations or even solution divergence. Second, the friction pressure drop model is highly dependent on empirical parameters, resulting in large errors under strongly unsteady and non-uniform flow patterns such as slug flow, affecting prediction accuracy and robustness. Third, it lacks adaptability to complex geometries; traditional methods assume the pipeline is a uniform straight pipe, making it difficult to handle complex topologies such as abrupt changes in inclination angle and cross-sectional area, and failing to accurately reflect the generation and propagation characteristics of slug flow under real-world conditions. Summary of the Invention

[0004] The technical problem to be solved by this invention is: addressing the issues of poor phase volume fraction conservation, reliance on empirical parameters for friction pressure drop models, and insufficient adaptability to complex geometries in existing one-dimensional simulation techniques for slug flow simulation, this invention provides a multiphase flow simulation method for pipelines that supports slug flow, thereby improving the accuracy and numerical stability of the simulation.

[0005] The present invention solves the above-mentioned technical problems through the following technical solution, including the following steps: S1. Geometric Information Input and Mesh Generation: Input the pipe geometric parameters and divide the pipe into several geometric segments along the axial direction according to the pipe cross-sectional area and inclination angle characteristics. Each geometric segment is discretized using an interlaced mesh system containing velocity control volume (VCV) and pressure control volume (PCV). VCV is used to solve the momentum equation, and PCV is used to solve the mass / pressure / energy equation. Define the inlet and outlet boundaries. S2. Boundary and Initial Condition Setting: Apply pressure or mass flow boundary conditions at the inlet and outlet boundaries, and set the initial field of VCV velocity distribution and PCV pressure, temperature, and volume fraction at the initial moment. S3. Calculation of physical property parameters and construction of interpolation table: Determine the composition of the transport medium, perform phase equilibrium calculation, solve the equilibrium state parameters of each phase under different temperatures and pressures, and construct a fast interpolation table; S4. Discretization and solution of governing equations: Based on the two-fluid model and cross-sectional averaging method, a closed one-dimensional set of governing equations is constructed, and spatial and temporal discretization is performed. The wall friction force and interphase drag force are calculated and corrected. After solving, the corrected volume fraction satisfies the conservation constraints, and the convergence is judged. S5. Time step control and termination judgment: The time step is dynamically updated according to the stability condition to advance the solution time. The calculation ends when the preset termination time is reached; otherwise, it returns to the S4 iteration.

[0006] Furthermore, the geometric parameters mentioned in S1 include the total number of geometric segments, the inner diameter / outer diameter / length / inclination angle / material properties of each segment; when the inner diameter of adjacent geometric segments changes, the abrupt change in area is handled by the conservation of gas-liquid two-phase flux; the height of the center point of each VCV is calculated according to the pipe inclination angle; when the VCV is at the connection of geometric segments with different inclination angles, the gravity term in the momentum equation is calculated using the volume of the VCV in the two segments as the interpolation weight.

[0007] Furthermore, the closed one-dimensional governing equations described in S4 include the continuity equation, momentum equation, energy equation, and pressure equation, as detailed below: Continuity equation: ; Momentum equation: ; Energy equation: ; Pressure equation: ; in, In the set The middle represents a phase of gas or liquid; These are time, distance along the flow path, volume fraction, density, flow velocity, pressure, internal energy, and enthalpy, respectively. These are the gravitational constants, and the pipe is in The height and inclination angle of the pipe; These are the phase transition rates, The frictional force between the phase and the pipe wall, relatively Interphase pull and tube wall heat transfer rate; The phase transition rate is calculated using the Bendiksen model. The interphase drag and pipe wall friction can be constructed using a hydraulic model of pipe flow, and their expressions are as follows:

[0008]

[0009]

[0010] in , These represent the volume of the VCV and the contact areas between the gas phase and the wall, between the gas and liquid phases, and between the liquid phase and the wall within the VCV, respectively. , These are the friction coefficient and the Reynolds number, respectively.

[0011] further,.

[0012] Furthermore, the interphase drag force described in S4 adopts a volume fraction adaptive correction mechanism: first, correction parameters are established. Its expression is: ; Further correction is made by adjusting the interphase friction coefficient. ; c is an empirical constant calibrated by the Bayesian optimization algorithm, preferably c=0.01.

[0013] Furthermore, the volume fraction correction described in S4 is performed after each iteration step, using a remapping strategy to correct the volume fraction, the expression of which is: .

[0014] Furthermore, in S4, the spatial discretization adopts the finite volume method, the convective flux is constructed using the upwind scheme, the pressure gradient is constructed using the central difference scheme, and the linear equations are solved using the chasing method; the time discretization adopts the implicit Euler scheme, and the pressure and momentum decoupling is implemented using the SIMPLE algorithm.

[0015] Furthermore, the equilibrium state parameters mentioned in S3 include the density, viscosity, mass fraction, partial derivative of density with respect to temperature and pressure, phase fraction, specific heat capacity, and heat transfer coefficient of each phase; when updating the physical property parameters in S4, the parameters are updated by traversing the PCV and using the Flash Table, and the gas-liquid two-phase viscosity at the center point of each VCV is calculated by linear interpolation.

[0016] Furthermore, in S4, the convergence criterion is that the initial residuals of each linear equation system are less than a preset threshold.

[0017] Beneficial effects: 1. This invention uses a volume fraction conservation correction mechanism to force the physical constraint that the sum of the volume fractions of the two phases equals 1, thereby avoiding numerical oscillations and solution divergence caused by non-physical solutions and improving the stability and reliability of the simulation.

[0018] 2. This invention optimizes the interphase drag force model by introducing an adaptive correction parameter for volume fraction, which solves the problem of underestimation of drag force under low volume fraction and enhances the robustness and convergence of the algorithm under extreme working conditions.

[0019] 3. This invention employs a segmented pipeline modeling method and a matched spatial discretization format to effectively handle complex geometric features such as pipeline undulations, bend connections, and cross-sectional changes, thereby improving the accuracy of slug flow simulation under complex operating conditions.

[0020] 4. The method adopted in this invention balances computational efficiency and prediction accuracy, and is suitable for full-scale simulation of slug flow in industrial-grade long-distance pipelines. It can provide strong support for pipeline system structure optimization, operating parameter control, and risk warning. Attached Figure Description

[0021] Figure 1 This is a schematic diagram of the overall process of the pipeline multiphase flow simulation method that supports slug flow.

[0022] Figure 2 This is a schematic diagram of the structure and parameters of the simulated object—a natural gas transmission pipeline—in Embodiment 2 of the present invention.

[0023] Figure 3 This is a graph showing the simulation results of the pipeline pressure at different time points in Embodiment 2 of the present invention.

[0024] Figure 4 This is a graph showing the simulation results of the gas-liquid two-phase volume fraction along the pipeline at different time points in Embodiment 2 of the present invention.

[0025] Figure 5 This is a graph showing the simulation results of the gas-liquid two-phase flow velocity along the pipeline at different time points in Embodiment 2 of the present invention. Detailed Implementation

[0026] The embodiments of the present invention are described in detail below. These embodiments are implemented based on the technical solution of the present invention, and provide detailed implementation methods and specific operation processes. However, the scope of protection of the present invention is not limited to the following embodiments.

[0027] Example 1 like Figure 1 As shown in the figure, this embodiment discloses a method for simulating multiphase flow in a pipeline that supports slug flow, including the following steps: Step 1: Input of geometric information and generation of mesh Input the geometric parameters of the pipeline, including the total number of geometric segments, the inner diameter, outer diameter, length, inclination angle, and material properties of each segment; based on the above geometric information, construct the pipeline computational mesh, specifically including the following sub-steps: (1) Mesh generation and control volume generation Each geometric segment is traversed sequentially, and a pressure control volume (PCV) is generated within each segment using an equally spaced partitioning method. Using the center point of the PCV as a reference, an interlaced mesh is constructed, and velocity control volumes (VCVs) are generated. The height of the center point of each VCV control volume is calculated based on the pipe inclination angle.

[0028] (2) Boundary condition definition The boundaries of the computational domain are defined by using the upstream boundary of the first PCV as the ingress boundary and the downstream boundary of the last PCV as the egress boundary.

[0029] Step 2: Setting Boundary and Initial Conditions Boundary conditions, including pressure boundary conditions or mass flow rate boundary conditions, are applied at the inlet and outlet boundaries of the computational domain. Simultaneously, the initial velocity distribution at each VCV center point and the initial fields such as pressure, temperature, and volume fraction at each PCV center point are set.

[0030] Step 3: Calculation of physical property parameters and construction of interpolation table The composition of the transported medium within the pipeline is determined, and phase equilibrium calculations (flash calculations) are performed to solve for the equilibrium parameters of each phase under different temperature and pressure conditions, including physical properties such as density, viscosity, and mass fraction. The calculation results are organized into a flash table for efficient lookup and updating of physical property parameters in subsequent calculations.

[0031] Step 4: Discretization and Solving of the Governing Equations The momentum, mass, energy, and pressure equations are spatially discretized and solved within a single time step, specifically including the following sub-steps: (1) Update of physical property parameters The process iterates through all PCVs within the pipeline, using the Flash Table constructed in step 3 to rapidly update physical properties such as the density of the gas and liquid phases, the partial derivative of density with respect to temperature and pressure, the phase fraction under equilibrium conditions, viscosity, specific heat capacity, and heat transfer coefficient. Simultaneously, a linear interpolation method is used to calculate the viscosity of the gas and liquid phases at the center point of each VCV.

[0032] (2) Calculation of drag force and friction force By iterating through the VCVs in the pipeline, the wall friction and interphase drag are calculated based on the updated viscosity values ​​of the formula (5-7) and physical property parameters. Then, according to the description in the physical model (3) in Section 4, the interphase drag is corrected using the formula (8-9) to ensure the accuracy and stability of the calculation when the phase fraction is small.

[0033] (3) Calculation of phase change rate and heat transfer The PCV within the pipe is traversed, and the phase change rate is calculated based on the partial derivatives of density with respect to temperature and pressure, the phase fraction under phase equilibrium state, and the Bendiksen model obtained in step 4(1). At the same time, the wall heat transfer Q is calculated using the wall heat transfer model.

[0034] (4) Pressure-velocity coupling solution A semi-implicit method for pressure-linked equations (SIMPLE) is employed to decouple pressure and momentum. The time term is discretized using the implicit Euler scheme, specifically including: [1] Momentum equation solution: Traverse all VCVs in the pipe and discretize the momentum equation (2) using the finite volume method, where the convective flux is constructed upwind. If the VCV is located at the junction of two geometric segments and the inclination angles of the geometric segments are different, the gravity phase in equation (2) is calculated using the volume of the VCV in the two geometric segments as the interpolation weight. After discretization, a tridiagonal matrix is ​​formed, and the linear equation system is solved using the pursuit method (Thomas algorithm) to obtain the estimated velocity of the center point of all VCVs.

[0035] [2] Solving the pressure correction equation: Traverse all PCVs in the pipeline, modify the pressure equation (4) using the SIMPLE algorithm, and establish the control equation for the pressure correction. Discretize the pressure correction equation using the finite volume method, where the pressure gradient is constructed using the central difference scheme. After discretization, a tridiagonal matrix is ​​formed, and the linear equation system is solved using the chasing method to obtain the corrected pressure at the center point of the PCV. Based on the corrected pressure, update the pressure value at the center point of the PCV and the velocity value at the center point of the VCV using the SIMPLE algorithm.

[0036] [3] Solving the volume fraction equation: Traverse all PCVs in the pipeline and use the finite volume method to discretize the mass equation (1), where the convective flux is constructed upwind. After discretization, a tridiagonal matrix is ​​formed, and the linear equation system is solved using the chasing method to obtain the gas and liquid volume fractions at the center points of all PCVs.

[0037] [4] Volume fraction constraint correction: Traverse all PCVs in the pipeline, and use formula (10) to remap and correct the volume fraction according to the method of physical model (4) in Section 4, to ensure that the physical constraint condition that the sum of the volume fractions of each phase is equal to 1 is met.

[0038] [5] Solving the energy equation: Traversing all PCVs in the pipeline, the energy equation (3) is discretized using the finite volume method, where the convective flux is constructed using an upwind scheme. After discretization, a tridiagonal matrix is ​​formed, and the linear equation system is solved using the chasing method to obtain the temperature values ​​of the center points of all PCVs.

[0039] [6] Update physical property parameters: Based on the updated temperature and pressure, repeat step 4 (1-3) to update physical property parameters, phase change rate, wall friction force and interphase drag force.

[0040] (5) Convergence judgment Repeat step 4(4) to monitor the initial residuals of each linear equation system. When all residuals are less than the preset convergence threshold, end the iterative solution of the current time step; otherwise, continue iterative calculation until the convergence condition is met.

[0041] Step 5: Time step control and convergence judgment The time step dt is dynamically updated based on the CFL (Courant-Friedrichs-Lewy) stability condition, and the current solution time T is advanced. It is then determined whether the current solution time has reached the preset termination time. If the termination condition is met, the calculation process ends; otherwise, it returns to step 4 to continue advancing the solution.

[0042] Formula (1) Continuity equation:

[0043] Equation (2) Momentum Equation:

[0044] Formula (3) Energy equation:

[0045] Formula (4) Pressure Equation:

[0046] subscript In the set The middle part represents a phase of gas or liquid. For time, distance along the flow path, volume fraction, density, flow velocity, pressure, internal energy, and enthalpy. As is the gravitational constant, the pipe is in The height of the location and the inclination angle of the pipe. The phase transition rate, The frictional force between the phase and the pipe wall, relatively The interphase pull and the heat transfer rate of the tube wall.

[0047] The phase transition rate is calculated using the Bendiksen model. The interphase drag and pipe wall friction can be constructed using a hydraulic model of pipe flow, and their expressions are as follows: Formula (5)

[0048] Formula (6)

[0049] Formula (7)

[0050] in , These represent the volume of the VCV and the contact areas between the gas phase and the wall, between the gas and liquid phases, and between the liquid phase and the wall within the VCV, respectively. , These are the friction coefficient and the Reynolds number, respectively.

[0051] Formula (8)

[0052] Formula (9)

[0053] Formula (10) .

[0054] Example 2 The implementation process of this invention will be described in detail below, taking the natural gas transmission pipeline shown in Figure 2 as an example: Pipeline parameter settings: The pipeline consists of four geometric sections. Sections 1 and 4 are horizontal, with an inner diameter of 0.186m, an outer diameter of 0.21m, and a length of 200m, made of stainless steel. Section 2 has an inclination angle of -5.8°, and section 3 has an inclination angle of 5.8°. Both sections have an inner diameter of 0.186m, an outer diameter of 0.21m, and a length of 25m, also made of stainless steel. The inlet boundary is for natural gas inflow at a temperature of 30℃ and a flow rate of 5.2kg / s; the outlet boundary is for natural gas outflow at a pressure of 45bar. Initial gas phase volume fraction: Section 1 is 1, Sections 2 and 3 are 0, and Section 4 is 0.5.

[0055] Simulated execution steps: Following step 1, input the above geometric parameters to generate a computational mesh. Generate equally spaced PCVs within each of the four geometric segments. Construct staggered meshes and VCVs based on the PCVs, calculate the height of the center point of each VCV, and define the boundaries of the entrance (the beginning of the first segment) and the exit (the end of the fourth segment).

[0056] Step 2: Set the inlet boundary as the mass flow rate boundary (5.2 kg / s) and temperature boundary (30℃), and the outlet boundary as the pressure boundary (45 bar). Set the initial velocity field to 0 and the initial pressure field to be uniformly distributed according to the outlet pressure of 45 bar.

[0057] Step 3: Determine the medium as natural gas (the main components are set according to the actual working conditions), perform phase balance calculations, and construct Flash Tables for different temperature (20-100℃) and pressure (10-100bar) ranges.

[0058] Step 4: Follow the procedure to update physical property parameters, calculate drag and friction (using formula 8-9 to correct interphase drag), calculate phase change rate and heat transfer, solve pressure-velocity coupled solution (SIMPLE algorithm), correct volume fraction and solve energy equation, set the convergence threshold to 10⁻⁶, until the residuals meet the requirements.

[0059] Step 5: Based on the CFL condition, set the initial time step dt=0.01s, dynamically adjust the step size, and advance the solution time to 42.2s.

[0060] Simulation results, such as Figure 3-5 As shown: Pressure simulation results: The pressure along the pipeline at different time points shows periodic fluctuations, which is consistent with the pressure change characteristics of slug flow. The fluctuation amplitude is in good agreement with the actual operating conditions, and there are no non-physical abrupt changes.

[0061] Phase fraction simulation results: The gas and liquid phase volume fractions are alternately distributed along the pipeline, forming a distinct slug structure. The sum of the volume fractions is always maintained at around 1, which verifies the effectiveness of the correction mechanism.

[0062] Simulation results of flow velocity: The gas phase velocity is generally higher than the liquid phase velocity. The velocity fluctuation in the slug region is significant, which matches the flow characteristics of slug flow, and there is no numerical divergence.

[0063] This embodiment verifies the accuracy and stability of the method of the present invention in slug flow simulation in complex geometric pipes, and can effectively reproduce the dynamic evolution process of slug flow.

Claims

1. A method for simulating multiphase flow in a pipeline that supports slug flow, characterized in that, Includes the following steps: S1. Geometric Information Input and Mesh Generation: Input the pipe geometric parameters and divide the pipe into several geometric segments along the axial direction according to the pipe cross-sectional area and inclination angle characteristics. Each geometric segment is discretized using an interlaced mesh system containing velocity control volume (VCV) and pressure control volume (PCV). VCV is used to solve the momentum equation, and PCV is used to solve the mass / pressure / energy equation. Define the inlet and outlet boundaries. S2. Boundary and Initial Condition Setting: Apply pressure or mass flow boundary conditions at the inlet and outlet boundaries, and set the initial field of VCV velocity distribution and PCV pressure, temperature, and volume fraction at the initial moment. S3. Calculation of physical property parameters and construction of interpolation table: Determine the composition of the transport medium, perform phase equilibrium calculation, solve the equilibrium state parameters of each phase under different temperatures and pressures, and construct a fast interpolation table; S4. Discretization and solution of governing equations: Based on the two-fluid model and cross-sectional averaging method, a closed one-dimensional set of governing equations is constructed, and spatial and temporal discretization is performed. The wall friction force and interphase drag force are calculated and corrected. After solving, the corrected volume fraction satisfies the conservation constraints, and the convergence is judged. S5. Time step control and termination judgment: The time step is dynamically updated according to the stability condition to advance the solution time. The calculation ends when the preset termination time is reached; otherwise, it returns to the S4 iteration.

2. The method according to claim 1, characterized in that, The geometric parameters described in S1 include the total number of geometric segments, the inner / outer diameter / length / inclination angle / material properties of each segment; when the inner diameter of adjacent geometric segments changes, the abrupt change in area is handled by the conservation of gas-liquid two-phase flux; the height of the center point of each VCV is calculated according to the pipe inclination angle; when the VCV is located at the connection of geometric segments with different inclination angles, the gravity term in the momentum equation is calculated using the volume of the VCV in the two segments as the interpolation weight.

3. The method according to claim 1, characterized in that, The closed one-dimensional governing equations described in S4 include the continuity equation, momentum equation, energy equation, and pressure equation, as detailed below: Continuity equation: ; Momentum equation: ; Energy equation: ; Pressure equation: ; in, In the set The middle represents a phase of gas or liquid; These are time, distance along the flow path, volume fraction, density, flow velocity, pressure, internal energy, and enthalpy, respectively. These are the gravitational constants, and the pipe is in The height and inclination angle of the pipe; These are the phase transition rates, The frictional force between the phase and the pipe wall, relatively Interphase pull and tube wall heat transfer rate; The phase transition rate is calculated using the Bendiksen model. The interphase drag and pipe wall friction can be constructed using a hydraulic model of pipe flow, and their expressions are as follows: ； ； in , These represent the volume of the VCV and the contact areas between the gas phase and the wall, between the gas and liquid phases, and between the liquid phase and the wall within the VCV, respectively. , These are the friction coefficient and the Reynolds number, respectively.

4. The method according to claim 3, characterized in that, The interphase drag force described in S4 adopts a volume fraction adaptive correction mechanism: first, the correction parameters are established. Its expression is: ； Further correction is needed for the interphase friction coefficient. ; c is an empirical constant calibrated by the Bayesian optimization algorithm, preferably c=0.

01.

5. The method according to claim 4, characterized in that, The volume fraction correction described in S4 is performed after each iteration step, using a remapping strategy. The expression for this correction is: 。 6. The method according to claim 1, characterized in that, The spatial discretization described in S4 uses the finite volume method, the convective flux is constructed using the upwind scheme, the pressure gradient is constructed using the central difference scheme, and the linear equations are solved using the chasing method; the time discretization uses the implicit Euler scheme, and the pressure and momentum decoupling is implemented using the SIMPLE algorithm.

7. The method according to claim 1, characterized in that, The equilibrium parameters mentioned in S3 include the density, viscosity, mass fraction, partial derivative of density with respect to temperature and pressure, phase fraction, specific heat capacity, and heat transfer coefficient of each phase; when updating the physical property parameters in S4, the parameters are updated by traversing the PCV and using the Flash Table, and the gas-liquid two-phase viscosity at the center point of each VCV is calculated by linear interpolation.

8. The method according to claim 1, characterized in that, In S4, the convergence criterion is that the initial residuals of each linear equation system are less than a preset threshold.

Images