Modernization of Todini Global Gradient Algorithm for hydraulic analysis of networks with choked flow

The modification of well-known Global Gradient Algorithm for hydraulic network flow distribution problem is proposed. This modification is based on problem equations rewritten in “upstream” form and on modified form of linearization, and can be effectively used for piping networks with gas and multiphase gas-liquid flow with multiple choked flow.


Introduction
Different methods of solving hydraulic network flow distribution problem were proposed [1,2], and Global Gradient Algorithm (GGA) was one of the most effective. GGA was proposed in [3][4][5][6] for liquid steady state flow, and was later extended for gas networks, networks with control valves [7,8], and some types of non-steady state flow problems [9][10][11]. Application of GGA for process piping was investigated by the author and his colleagues in a number of publications [12][13][14], and PASS/HYDROSYTEM software [15], developed by PSRE Co, uses GGA successfully for thermal and hydraulic analysis of piping networks transporting both one-phase (liquid, gas) and multiphase phase gas-liquid fluids.
However, in its original form GGA cannot be applied to networks with choked or nearchoked flow -while analysis of such flow is an important practical problem for some types of process piping -for example pressure relieve discharge systems or multiphase gas-liquid transfer pipelines transporting oil from furnace to column. Such flow was investigated by the author in [16][17][18][19]. He also established conditions for existence and uniqueness of flow distribution problem solution for this type of flow [19].
This article proposes modified form of GGA (MGGA), which (along with decomposition method) hopefully allows to solve flow distribution problem for choked flows. Further in this article we will consider on hydraulic analysis, i.e. will consider 2 nd thermodynamic parameter (besides pressure) describing state of the fluid (for example temperature for isothermal flow or full enthalpy for adiabatic flow) to be fixed.
(2) Where -incidence matrix of network graph, -vector of mass flow rates, -node mass inflow vector, -vector of node potentials (pressures), ( ) -vector function, each element of which is function of losses on edges.
It is supposed that node pressures are set in part of nodes (at least one), and inflow are defined in remaining nodes. We can rewrite matrices and vectors in (1), (2) as Where index 0 corresponds to nodes with set potentials, and index 1 -nodes with defined inflows. From (1), (2) GGA searches solution of non-linear system of equations (4), (5) via and � , without any a-priory restrictions on vectors and in system (1) (this is the reason of the name "global"). On each iteration system (4), (5) is replaces by its linearization in the vicinity of current iteration, and linearized system is solved. Taking into account that (4),(5) and (6) GGA equations can be reduced is Maxwell matrix -symmetrical positive-definite sparse diagonally dominant М-matrix [20][21][22][23][24], which inversion in (7) can be done computationally effective.
GGA has some significant advantages: it converges very quickly, does not demand any preliminary topological analysis of the network, and also almost not sensitive to initial approximation quality.
However, in real process piping function ( ) can depends not only on flow rate, but also on node pressure values [16][17][18][19]. While this dependence is weak (| ⁄ | ≪ | ⁄ |), GGA works effectively. The more solution is close to choked flow, the worse is GGA convergency. And the worst difficulty is that iterations more and more often fall in the region where functions of losses on edges are not defined -some edges simply cannot transport current iteration flow rate value, as this value exceeds critical flow at current pressure in start edge node.

MGGA equations
To work around this problem, let's reformulate hydraulic network equations, switching from "downstream" calculation to "upstream" calculation: (17) Matrix * is not symmetrical, but (in cases, when solution is unique-see [19]) still is sparse nonsingular WCDD M-matrix and can be inversed using effective calculation methods.
To preliminary check proposed MGGA, a series of calculation experiments was performed in cooperation with D.D. Fedyunina on small hydraulic networks. Isothermic ideal gas choked flow was modeled (in this case function * and its derivatives can be calculated analytically). The results of experiments demonstrated that proposed method is working and keeps all its advantages -quick convergency and weak dependency on initial point. Moreover, with pressure as node potentials -not squares of pressures (which usually are used for gas networks) the convergence is even better. This can be explained by the fact that in case of choked flow, edge flow rate is proportional to start node pressure.
Currently the work is performing on implementation of general case of "upstream" calculation (with possible multiple choked flow and complete fluid vaporization or condensation) in PASS/HYSROSYSTEM program. In the same time additional computational experiments for adiabatic gas flow (Fanno flow) are scheduled, and on the base of their results MGGA implementation in PASS/HYDROSYSTEM is scheduled.