Equation-oriented multiphase flow pipe network steady-state modeling method

Abstract

The invention discloses an equation-oriented multiphase flow pipe network steady-state modeling method, and relates to the technical field of multiphase flow pipe network simulation. The method is executed by a computer system, an abstract pipe network is used as a node-branch mathematical model, node pressure and branch mass flow are defined as independent variables, and a nonlinear equation set comprising a branch pressure drop equation, a node mass conservation equation and a boundary condition equation is constructed; the method comprises the following steps: initializing based on a linearized Darcy equation, carrying out variable scaling, solving by adopting a Newton iteration method in combination with a line search technology, integrating a pipeline hydraulic model, a black oil model and a separator node model in the solving process, and processing pipeline pressure-temperature coupling calculation through a Runge-Kutta self-adaptive step length algorithm. The method solves the problems of poor convergence, difficulty in treating reverse flow and complex loop of the existing ordered through module method, has the advantages of high convergence efficiency, high robustness and wide adaptability, and provides technical support for pipe network design and optimization in the fields of oil and gas gathering and transportation, chemical process and the like.

Description

Technical Field

[0001] This invention relates to the field of steady-state simulation technology for multiphase flow pipeline networks, specifically to an equation-oriented steady-state modeling method for multiphase flow pipeline networks, applicable to the design, optimization, and safety assessment of complex multiphase flow pipeline networks such as oil and gas gathering and transportation networks and chemical process pipelines. Background Technology

[0002] Multiphase flow pipeline systems are core infrastructure in the energy and process industries. Steady-state simulation of these systems requires determining node pressure, branch flow rate, temperature, and phase distribution to provide crucial data for system design and operation. Due to the coupled effects of friction, gravity, and interphase mass and heat transfer within the pipeline network, the gas-liquid (oil, water) multiphase flow exhibits strong nonlinearity, making steady-state solutions a challenging engineering problem.

[0003] In existing technologies, the sequential modular method is widely used for steady-state solutions of multiphase flow pipe networks. This method divides the pipe network into independent computational units and calculates them one by one according to a preset topology order, coordinating the pressure balance and mass conservation of nodes through outer-layer iteration. However, this method has inherent defects: First, convergence depends on the selection of initial values, and it is prone to oscillation or divergence when dealing with strongly nonlinear and high gas phase rate conditions, resulting in poor robustness; second, it depends on the preset flow direction, and the solution is interrupted when the actual flow direction does not match the preset direction, making it unable to adapt to dynamic flow direction changes; third, local solutions indirectly handle unit coupling relationships, resulting in low efficiency for large and complex pipe networks; and fourth, it is difficult to coordinate the modeling when dealing with complex equipment with multiple outlets such as separators. Therefore, there is an urgent need for a globally synchronous, efficient, and robust steady-state modeling method for multiphase flow pipe networks. Summary of the Invention

[0004] The technical problem to be solved by this invention is to provide an equation-oriented steady-state modeling method for multiphase flow pipe networks, which addresses the problems of poor convergence, difficulty in handling reverse flow and complex loops, low computational efficiency and insufficient modeling flexibility of existing sequential module methods, so as to achieve efficient, accurate and robust steady-state solutions for complex multiphase flow pipe networks.

[0005] This invention solves the above-mentioned technical problems through the following technical solution. This invention is executed by a computer system integrating numerical calculation, physical property database, and pipeline topology processing modules, and includes the following steps: S1. Pipeline System Abstraction and Variable Definition: Receive the pipeline topology, equipment parameters, and boundary conditions, abstract them into a mathematical graph model composed of nodes and branches, and define the pressure of all nodes. (j=1,2,..., , (Total number of nodes) and quality flow of all branches (k=1,2,..., , (Total number of branches) are independent solution variables. Mass flow rate, positive and negative indicate the direction of flow; S2. Construction of Nonlinear Equations: Construction includes... A set of independent equations, the set of equations including the pressure drop equations of each branch, the mass conservation equations of each internal node, and the constraint equations of each boundary node. S3. System Initialization and Variable Scaling: Based on the linearized Darcy equation, the mass balance and boundary condition equations are solved simultaneously, and the reference physical property parameters provided by the black oil model are used to obtain the initial estimates of pressure and flow rate; the scaling factor is calculated. and All variables are dimensionless. S4. Newton Iteration Solution: Based on the Newton iteration method combined with line search technology, solve the nonlinear equation system. During the iteration process, the residual vector and Jacobian matrix are calculated simultaneously. Solve the linear equation system to obtain the variable correction step size and update the variables until the residual norm is less than the preset tolerance, and output the steady-state parameters of the pipeline network.

[0006] Furthermore, the nodes mentioned in S1 include internal nodes and external nodes. Internal nodes are multi-branch intersection points, and external nodes are boundary nodes. Branches include pipes or equipment, and equipment includes at least separators and control valves.

[0007] Furthermore, the residual form of the pressure drop equation for the branch described in S2 is as follows: Based on a pipeline hydraulic model, which includes coupled pressure gradient equations and temperature gradient equations: The expression for the pressure gradient equation is: , Where P is the pressure along the pipeline (Pa), x is the length along the pipeline (m), f encapsulates the frictional pressure drop, gravity pressure drop, and acceleration pressure drop, and the mixing density is calculated based on the black oil model. Physical properties such as viscosity (μ); The expression for the temperature gradient equation is: , Where T is the fluid temperature (K), d is the pipe diameter (m), and u is the overall heat transfer coefficient (W / m). 2 •K), The ambient temperature (K) is the ambient temperature. For gravity-related heat transfer terms, Specific heat capacity at constant pressure of the fluid (J / kg·K).

[0008] Furthermore, the solution to the pipeline hydraulic model employs the Runge-Kutta adaptive step-size control algorithm, which includes two calculation modes: CPOUT mode: Given inlet pressure and flow rate, the positive integral pressure P and the pressure-to-flow derivative are calculated. Derivative of pressure with respect to inlet pressure And specific enthalpy H, simultaneously calling the black oil model to perform isenthalpic flash evaporation to obtain temperature T and ; CPIN mode: Given the outlet pressure and flow rate, first based on the inlet temperature... Estimate the initial temperature profile, and then solve iteratively through "reverse pressure integration - forward temperature profile calculation" until the inlet pressure or temperature profile change is less than 10⁻⁵ K.

[0009] Furthermore, the initial step size of the Runge-Kutta adaptive step size control algorithm is 5% of the pipe segment length. The error is estimated by comparing the results of the fourth-order and fifth-order Runge-Kutta formulas. When the error is small, the step size is increased by 5 times, and when the error is large, the step size is reduced. Sub-step points are set at 20%, 30%, 60%, 100%, and 87.5% for each calculation segment.

[0010] Furthermore, the pipeline hydraulic model also includes special operating condition handling: Reverse flow processing: When Q<0 is detected, the pipe segment sequence, elevation and ambient temperature distribution are reversed to convert the reverse flow into forward flow calculation; Critical flow treatment: Monitor flow velocity and Mach number, and limit the velocity to not exceed the speed of sound using a decay function; Negative pressure treatment: Internal negative pressure is maintained to ensure numerical continuity, and a limiting pressure slightly higher than zero is used when calling the physical property model.

[0011] Furthermore, the residual form of the mass conservation equation for the internal nodes described in S2 is: ; If it is a separator node, the equation expands to the total mass balance equation and Phase separation equations ,in The outlet flow rates of gas, oil, and water are obtained through isenthalpic flash evaporation calculations using the black oil model.

[0012] Furthermore, the boundary condition equations in S2 include a constant pressure. Fixed flow Or pressure-flow curve Three forms.

[0013] Furthermore, the linearized Darcy equation in S3 is: , in , , The coefficients are estimated based on the initial density ρ and viscosity μ provided by the black oil model under preset reference conditions; the dimensionless processing formula is as follows: , .

[0014] Furthermore, the Jacobian matrix described in S4 During the calculation process: Partial derivatives of the branched pressure drop equation , Partial derivatives of physical properties provided by combining the pipeline hydraulic model with the black oil model. , Solve this problem; The partial derivatives of the separator node equations are solved by using the partial derivatives of the outlet flow rate with respect to the feed pressure, composition, and enthalpy provided by the separator node model and the black oil model. Variable updates use , where λ is the step size factor determined by the line search technique.

[0015] Beneficial effects: 1. This invention transforms the pipeline system into a unified nonlinear equation set and solves it using Newton's method with quadratic convergence characteristics. The Jacobian matrix contains global variable coupling information, making the iteration direction more accurate and overcoming the shortcomings of linear convergence and error accumulation in the sequential module method. For pipelines with strong nonlinearity, high gas phase rate, or complex loops, the convergence success rate is increased from below 60% in the prior art to over 95%, and it does not rely on the accuracy of the initial value, greatly reducing the operational threshold. 2. This invention defines branch mass flow rate as a signed algebraic quantity, and reverse flow is only represented by negative flow rate. With the automatic reversal mechanism of pipe profile, there is no need to preset flow direction or complex logic judgment, and the reverse flow problem is uniformly transformed into forward flow calculation. Compared with the existing technology, which requires interruption and adjustment of the model, this invention can automatically adapt to any flow direction change condition, and improve the continuity of solution by 100%. It is especially suitable for simulation of dynamic conditions such as start-up and shutdown and load fluctuation. 3. This invention reduces the number of iterations by more than 70% compared to the outer loop count of the sequential module method by leveraging the quadratic convergence characteristic of Newton's method; the Runge-Kutta adaptive step-size algorithm maximizes the integration step size while ensuring accuracy, reducing the number of nodes in the pipeline profile calculation and improving the efficiency of single residual evaluation by 40%-60%; global synchronous solution avoids the repeated module traversal of the sequential method, and for large pipeline networks with more than 100 nodes and more than 200 branches, the total calculation time is reduced by 60%-80%.

[0016] 4. This invention deeply integrates the black oil model and the separator node model, supports isenthalpic flash evaporation, phase separation calculation, and virtual connector modeling of control valves, and is adaptable to various scenarios such as oil and gas mixed transportation and chemical medium transportation; linearized Darcy equation initialization provides a reasonable iteration starting point, variable scaling improves the condition number of the Jacobian matrix, and line search technology avoids oscillations caused by excessive step size. The three factors work together to reduce the probability of solution failure by more than 80%, and can stably handle extreme conditions such as critical flow and negative pressure. Attached Figure Description

[0017] Figure 1 This is a schematic diagram of the process of this invention.

[0018] Figure 2 This is a schematic diagram of the variables, equations, and Jacobian matrix.

[0019] Figure 3 This is a schematic diagram of a pipe heat transfer model. Detailed Implementation

[0020] 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.

[0021] like Figure 1 As shown in the figure, this embodiment discloses an equation-oriented steady-state modeling method for multiphase flow pipeline networks, integrating numerical calculation, physical property database, and pipeline network topology processing modules. The core of the system lies in adopting an equation-oriented architecture, transforming the physical laws of the entire pipeline network into a large-scale nonlinear equation system and solving it simultaneously.

[0022] The system is built upon and relies on three key fundamental physical models: the black oil model, the pipeline hydraulic model, and the separator node model. These models are dynamically invoked throughout the solution process, providing accurate physical relationships and data support for the construction and solution of the equations.

[0023] Step 1: Pipeline system abstraction and variable definition.

[0024] The system receives user input of pipeline topology, equipment parameters, and boundary conditions, and abstracts them into a mathematical graph model consisting of nodes and branches; where nodes represent pipeline junctions or boundary points, and each node is assigned a pressure variable. (j=1,2,..., , (Total number of nodes), where the number of internal nodes is The number of external nodes is , A branch represents a pipe or device connecting two nodes, and each branch is assigned a mass flow rate variable. (k=1,2,..., , (Total number of branches); mass flow rate Mass flow rate, with its positive and negative signs indicating the direction of flow; the purpose of this step is to discretize the physical system into solvable variables with well-defined mathematical characteristics.

[0025] Step 2: Constructing the nonlinear equation system.

[0026] Based on the variables defined in the first step, the system automatically constructs a list containing... A system of three independent equations. This system of equations consists of the following three types of equations: (1) Branch pressure drop equation: The system constructs a pressure drop equation for each branch, and its general residual form is as follows: The solution to this equation is heavily dependent on the pipeline hydraulic model; when the system needs to calculate the residuals and partial derivatives of this equation, it will call the pipeline model.

[0027] The pipe model is a comprehensive model that couples hydraulics (pressure) and thermals (temperature). Its core lies in solving a set of one-dimensional steady-state differential equations describing fluid flow and heat transfer within the pipe, and providing accurate residuals and their partial derivatives for the Newton iteration of the network layer.

[0028] ① Model Establishment. The system establishes two core differential equations for each pipe segment, forming a coupled system: Hydraulic model, i.e., pressure gradient equation

[0029] in, P is the pressure along the pipe (Pa), which is the main dependent variable to be solved. x is the length (m) along the geometric direction of the pipe, and is the independent variable. Q is the mass flow rate (kg / s), determined by the network solver at the branch level, and is used as a known quantity in the pipeline calculation. The function f represents the total pressure drop gradient, which encapsulates the frictional pressure drop, gravitational pressure drop, and acceleration pressure drop. Its specific form is determined by the chosen multiphase flow correlation (such as Darcy-Weisbach, Beggs & Brill, etc.) and depends on physical properties (such as mixing density) calculated in real-time by the black oil model. (e.g., viscosity μ).

[0030] Thermodynamic model, i.e., temperature gradient equation

[0031] in, T is the fluid temperature (K), which is the dependent variable to be solved. d is the pipe diameter (m). u is the overall heat transfer coefficient (W / m²) 2 ·K) cp is the specific heat capacity at constant pressure of the fluid (J / kg·K), provided by the black oil model. T0 is the ambient temperature (K). g is the acceleration due to gravity, g / s². θ is the angle of inclination of the pipe relative to the horizontal plane. ② Model solution.

[0032] The solution strategy for the pipeline model depends on the hydraulic calculation mode (CPIN or CPOUT) used, which determines how the two coupled variables of pressure and temperature are solved in a coordinated manner.

[0033] Numerical integration algorithm. Regardless of the mode, the system employs the Runge-Kutta (RK_ASLC) algorithm with adaptive step size control for integrating the pressure gradient equation. This algorithm automatically divides the user-defined pipe segment into a variable number of computational segments of varying lengths to ensure integration accuracy and efficiency. The initial step size is set to 5% of the current pipe segment length, and the error is estimated by comparing the results of the fourth and fifth order Runge-Kutta formulas. If the error is small, the next step size is increased to 5 times; if the error is large, the step size is decreased and the process is retried. The system uses a smoothing mixture function to ensure continuous step size changes and avoid discrete jumps. Within each computational segment, the Cash-Karp method performs evaluations at five fixed sub-step points (located at 20%, 30%, 60%, 100%, and 87.5% of the segment length) to provide high-precision results.

[0034] In CPOUT mode, given the inlet pressure and flow rate, the outlet pressure is calculated. The system performs forward integration starting from the pipe inlet and proceeding along the flow direction.

[0035] The system integrates the following four variables simultaneously: y1 = P (pressure) y2 = ∂P / ∂Q (Pressure derivative with respect to flow rate) y3 = ∂P / ∂Pin (Derivative of pressure with respect to inlet pressure) y4 = H (specific enthalpy, used to derive temperature) In this mode, temperature and pressure are integrated synchronously in the positive direction. In each integration step, the system calls the black oil model, uses the current pressure P and enthalpy H to perform isenthalpic flash evaporation, and obtains parameters such as temperature T and cp required to update physical properties.

[0036] CPIN mode calculates the inlet pressure given the outlet pressure and flow rate. This mode is the system default and more stable. Because this is a boundary value problem, and the temperature transfer direction is opposite to the inverse integration direction, iterative solutions are required.

[0037] The system first uses the analytical solution of the thermodynamic equation to estimate an initial temperature profile along the flow direction, starting from the known inlet temperature Tin. The system then enters a pressure-temperature iterative cycle: The system uses the currently estimated temperature profile to perform pressure integration against the flow direction, starting from the known outlet pressure Pout. This process calculates the pressure P and its derivatives y2 and y3, thereby obtaining a new inlet pressure estimate and pressure profile.

[0038] Using the new pressure profile obtained in the previous step, the system recalculates the accurate temperature profile along the flow direction using the thermodynamic equation (or its analytical solution), starting from the inlet.

[0039] Repeat the above steps until the change in inlet pressure (or the entire temperature profile) obtained from two consecutive iterations is less than the specified tolerance (e.g., 10⁻⁵ K).

[0040] ③ Reverse flow, critical flow, and negative pressure treatment.

[0041] When the system detects a negative mass flow rate Q, it automatically performs a pipe profile reversal. The original pipe segment sequence, elevation, ambient temperature distribution, etc., are completely inverted. Subsequently, all calculations (whether CPIN or CPOUT) are based on this reversed profile, thus transforming the reverse flow problem into an equivalent forward flow problem. During integration, the system monitors the flow velocity and Mach number. To prevent non-physical supersonic flow, the system uses a decay function to smooth and limit the velocity, keeping it close to but not exceeding the calculated speed of sound. Similarly, if negative pressure occurs during iteration, the system internally continues to use this negative pressure to ensure numerical continuity, but when calling the property model, it uses a slightly higher than zero limiting pressure processed by the decay function to ensure the rationality of the property calculations. (2) Node mass balance equation: The system constructs a mass conservation equation for each internal node, and its residual form is as follows: For conventional internal nodes, the equation only involves summing the flow rates. However, for separator nodes, the system invokes the separator node model. This model, supported by the black oil model, performs isenthalpic flash evaporation calculations on the mixed fluid flowing into the separator to determine the flow rates Qg, Qo, and Qw of each outlet phase (e.g., gas, oil, water). The mass balance equation for the separator node is therefore extended to include the total mass balance. and Phase separation equations The purpose of this step is to enforce the conservation of mass. For the separator node, the separation effect is determined by both the separator node model and the black oil model.

[0042] (3) Boundary condition equations: The system constructs constraint equations for each boundary node, in the form of fixed pressure. Fixed flow or pressure-flow curve The purpose of this step is to provide external drives or constraints for the system.

[0043] Step 3: System Initialization and Variable Scaling. The system approximates each branch using a linearized Darcy equation: The estimation of coefficients K1,k,K2,k,Bk requires initial fluid density ρ and viscosity μ, which are provided by the black oil model under preset reference conditions. The system simultaneously solves all mass balance and boundary condition equations to form a linear system of equations, obtaining initial estimates of pressure and flow rate. Subsequently, the system calculates scaling factors SP and SQ based on the initial values ​​and performs dimensionless processing on all variables. , The purpose of this step is to provide reasonable initial values ​​for iteration and improve numerical stability; its physical property basis comes from the black oil model.

[0044] Step 4: Newton's Iterative Solution Core Loop: The system enters a nonlinear iterative loop centered on Newton's method, with the goal of solving for R(x) = 0, where x is the vector of all scaled variables. Step 4 is divided into the following sub-modules: Residual and Jacobian Matrix Calculation: In each iteration, the system calculates the residual vector R(x) and the Jacobian matrix J for the current variable x. ).

[0045] This process deeply integrates three basic models: Calculate the branched pressure drop equation Rp,k and its partial derivatives. , At that time, the system invokes the pipeline hydraulic model. This model, while performing pressure traversal, utilizes the physical properties and partial derivatives provided by the black oil model. , By using analytical or numerical differentiation techniques, the required partial derivatives can be calculated simultaneously.

[0046] When calculating the separator node mass balance equation Rq,sep and its partial derivatives, the system calls the separator node model. Supported by the black oil model, this model not only performs flash evaporation calculations but also provides the partial derivatives of the outlet flow rate with respect to the feed pressure, composition, and enthalpy.

[0047] The black oil model is frequently invoked as a common sub-model throughout the process, providing real-time, accurate, and derivative-continuous fluid properties for the pipeline and separator models.

[0048] This module is the core of the three basic models working together, and its role is to construct the search direction required for Newton's iteration with high accuracy.

[0049] Solving Linear Systems and Updating Variables: Solving Systems of Linear Equations The corrected step size Δx is obtained, and the step size factor λ is determined using a line search technique. Update variables; this module's function is to intelligently update variables to ensure convergence.

[0050] Convergence check. The system calculates the norm of the updated residuals. If it is less than the preset tolerance, the iteration terminates and the results are output (including the complete pressure and temperature profile calculated by the pipeline hydraulic model, the flow rates of each phase calculated by the separator node model, and the physical properties of each point given by the black oil model); otherwise, the loop is repeated. The function of this module is to determine whether the solution is complete.

[0051] Although embodiments of the present invention have been shown and described above, it is understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those skilled in the art can make changes, modifications, substitutions and variations to the above embodiments within the scope of the present invention.

Claims

1. A method for steady-state modeling of multiphase flow pipe networks oriented towards equations, characterized in that, Performed by a computer system integrating numerical calculation, physical property database, and pipeline topology processing modules, the process includes the following steps: S1. Pipeline System Abstraction and Variable Definition: Receive the pipeline topology, equipment parameters, and boundary conditions, abstract them into a mathematical graph model composed of nodes and branches, and define the pressure of all nodes. (j=1,2,..., , (Total number of nodes) and quality flow of all branches (k=1,2,..., , (Total number of branches) are independent solution variables. Mass flow rate, positive and negative indicate the direction of flow; S2. Construction of Nonlinear Equations: Construction includes... A set of independent equations, the set of equations including the pressure drop equations of each branch, the mass conservation equations of each internal node, and the constraint equations of each boundary node. S3. System Initialization and Variable Scaling: Based on the linearized Darcy equation, the mass balance and boundary condition equations are solved simultaneously, and the reference physical property parameters provided by the black oil model are used to obtain the initial estimates of pressure and flow rate; the scaling factor is calculated. and All variables are dimensionless. S4. Newton Iteration Solution: Based on the Newton iteration method combined with line search technology, solve the nonlinear equation system. During the iteration process, the residual vector and Jacobian matrix are calculated simultaneously. Solve the linear equation system to obtain the variable correction step size and update the variables until the residual norm is less than the preset tolerance, and output the steady-state parameters of the pipeline network.

2. The method according to claim 1, characterized in that, The nodes described in S1 include internal nodes and external nodes. Internal nodes are multi-branch intersections, and external nodes are boundary nodes. Branches include pipes or equipment, and equipment includes at least separators and control valves.

3. The method according to claim 1, characterized in that, The residual form of the pressure drop equation for the branch described in S2 is: Based on a pipeline hydraulic model, which includes coupled pressure gradient equations and temperature gradient equations: The expression for the pressure gradient equation is: , Where P is the pressure along the pipeline (Pa), x is the length along the pipeline (m), f encapsulates the frictional pressure drop, gravity pressure drop, and acceleration pressure drop, and the mixing density is calculated based on the black oil model. Physical properties such as viscosity (μ); The expression for the temperature gradient equation is: , Where T is the fluid temperature (K), d is the pipe diameter (m), and u is the overall heat transfer coefficient (W / m). 2 •K), The ambient temperature (K) is the ambient temperature. For gravity-related heat transfer terms, Specific heat capacity at constant pressure of the fluid (J / kg·K).

4. The method according to claim 3, characterized in that, The solution to the pipeline hydraulic model employs the Runge-Kutta adaptive step-size control algorithm, which includes two calculation modes: CPOUT mode: Given inlet pressure and flow rate, the positive integral pressure P and the pressure-to-flow derivative are calculated. Derivative of pressure with respect to inlet pressure And specific enthalpy H, simultaneously calling the black oil model to perform isenthalpic flash evaporation to obtain temperature T and ; CPIN mode: Given the outlet pressure and flow rate, first based on the inlet temperature... Estimate the initial temperature profile, and then solve iteratively through "reverse pressure integration - forward temperature profile calculation" until the inlet pressure or temperature profile change is less than 10⁻⁵ K.

5. The method according to claim 4, characterized in that, The initial step size of the Runge-Kutta adaptive step size control algorithm is 5% of the pipe segment length. The error is estimated by comparing the results of the fourth-order and fifth-order Runge-Kutta formulas. When the error is small, the step size is increased to 5 times. When the error is large, the step size is reduced. Sub-step points are set at 20%, 30%, 60%, 100%, and 87.5% for each calculation segment.

6. The method according to claim 3, characterized in that, The pipeline hydraulic model also includes special operating condition handling: Reverse flow processing: When Q<0 is detected, the pipe segment sequence, elevation and ambient temperature distribution are reversed to convert the reverse flow into forward flow calculation; Critical flow treatment: Monitor flow velocity and Mach number, and limit the velocity to not exceed the speed of sound using a decay function; Negative pressure treatment: Internal negative pressure is maintained to ensure numerical continuity, and a limiting pressure slightly higher than zero is used when calling the physical property model.

7. The method according to claim 4, characterized in that, The residual form of the mass conservation equation for the internal nodes described in S2 is: ; If it is a separator node, the equation expands to the total mass balance equation and Phase separation equations ,in The outlet flow rates of gas, oil, and water are obtained through isenthalpic flash evaporation calculations using the black oil model.

8. The method according to claim 1, characterized in that, The boundary condition equations in S2 include constant pressure. Fixed flow Or pressure-flow curve Three forms.

9. The method according to claim 1, characterized in that, The linearized Darcy equation in S3 is , in , , The coefficients are estimated based on the initial density ρ and viscosity μ provided by the black oil model under preset reference conditions; the dimensionless processing formula is as follows: , .

10. The method according to claim 1, characterized in that, Jacobian matrix as described in S4 During the calculation process: Partial derivatives of the branched pressure drop equation , Partial derivatives of physical properties provided by combining the pipeline hydraulic model with the black oil model. , Solve this problem; The partial derivatives of the separator node equations are solved by using the partial derivatives of the outlet flow rate with respect to the feed pressure, composition, and enthalpy provided by the separator node model and the black oil model. Variable updates use Δx , where λ is the step size factor determined by the line search technique.

Images