Generalization of Methods for Calculating Steady-State Flow Distribution in Pipeline Networks for Non-Conventional Flow Models

This paper proposes generalized models and methods for calculating flow distribution in hydraulic circuits with lumped parameters. The main models of the isothermal steady-state flow of medium are classified by an element of the hydraulic circuit. These models include conventional, implicitly specified by flow rate, and pressure-dependent ones. The conditions for their applicability, which ensure the existence and uniqueness of a solution to the flow distribution problem, are considered. We propose generalized nodal pressure and loop flow rate methods, which can be applied regardless of the forms of specific element models. Final algorithms, which require lower computational costs versus the known approaches designed for non-conventional flow models, are substantiated. Proposed models, methods, algorithms, and their capabilities, are analytically and numerically illustrated by an example of a fragment of gas transmission network with compressor stations.


Introduction
The problems of flow distribution are the fundamental problems of analysis and justification of decisions on the management of operation modes of pipeline systems (PLSs) of various types and purposes (heat-, water-, oil-, gas supply, etc.) in their design, operation, and supervisory control. The requirements for the efficiency of the methods of solving these problems are constantly increasing. Such requirements include the following: (1) speed; (2) reliability, which manifests itself in the guaranteed solutions with a predetermined accuracy; (3) ability to solve high-dimensional problems; (4) universality with respect to an arbitrary structure of the object of the analysis, the composition of its elements, and changes in design conditions. The expediency of development and elaboration of unified methods for modeling operation modes of PLSs of various types is due to the commonality of the following: (1) their structural and topological properties; (2) physical laws of liquid (gas) flow in individual elements; (3) conservation laws of networks; (4) conceptual and mathematical statements of analysis problems. Therefore, in what follows, when considering models and methods, we will use the term "hydraulic circuit" (HC), irrespective of the specific purpose of the PLSs being modeled.
The conventional model of steady-state isothermal flow distribution in the HC with concentrated parameters includes Kirchhoff laws and closing relations. It represents a system of linear and nonlinear algebraic equations, and there are two basic forms of writing such models: the nodal and loop ones. The nodal model is as follows [1,2]: where A-a complete (m × n)-matrix of incidences of nodes and branches of the directed graph of the analytical circuit model with elements of a j i = 1(−1), if node j is initial (final) one for branch i and a j i = 0, if branch i is not incident to node j; x, y-n-dimensional vectors of flow rates and pressure drops at the branches of the analytical circuit model; A-a [(m − 1) × n]-matrix of incidences formed from A by crossing out one of the linearly dependent lines; f (x)-a n-dimensional vector-function with elements f i (x i ), i = 1, n, reflecting the laws of the pressure drop as a function of flow rate (flow model) at the branches of the HC; Q-a (m − 1)-dimensional vector of nodal flow rates with elements Q j > 0 (Q j < 0) if node j has an inflow (extraction) and Q j = 0, if node j is a simple branch connection point; P-a m-dimensional vector of nodal pressures. The loop model differs only in the way the second Kirchhoff law is written: instead of the second relation in (1) one uses By = 0. Here (B) is a (c × n)-matrix of loops, with elements b r i = 1(−1), if the orientation of branch i coincides (does not coincide) with the direction of traversal of loop r and b r i = 0, if branch i is not included in loop r; c is the number of main loops equal to the number of branches of the HC not included in the spanning tree of the HC, and called chords. Hence c = n − m + 1. The equivalence of both forms of notation follows from the fact that B(A T P = y) ⇒ 0 = By due to the well-known property of orthogonality of matrices A and B [3].
Let P = P P m , then the classical problem of flow distribution consists in determining vectors x, y and P given the known form of f i (x i ), i = 1, n, predetermined matrix A, vector Q, and pressure at one of the nodes (P m ), which for simplicity's sake is assumed to be zero. In the loop model, vector P is determined after solving the problem with respect to x and y. Numerous methods and algorithms for solving this problem are known, an overview of which can be found, for example, in [1,4], etc. However, the main ones are the classical methods: those of node and loop methods (NM and LM) [1,2,[5][6][7][8], etc., which have been widely adopted in modeling the operation modes of PLSs of heat- [1,9,10], etc., water- [1,9], etc., gas supply [11][12][13], etc.
Both methods are based on the Newton-Raphson method in conjunction with special methods of decreasing the order of the system of linearized equations to be solved. Let us introduce decomposition x = x C x T , where x C , x T are the flow rate vectors at the chords and branches of the spanning tree, respectively. Finding the direction of the Newton-Raphson decrement to the flow rates at the chords (∆x k C ) in the LM and to the nodal pressures (∆P k ) in the NM at the k-th iteration involves solving the following systems of equations: where f x = ∂ f ∂x is the diagonal matrix of partial derivatives at point x k . In both systems the coefficient matrices are symmetric. For two-dimensional circuits (which is typical of most real-life PLSs), matrix (2) of coefficients of the system has a strict diagonal predominance, and in (3) it is not strict. Therefore, the LM usually requires fewer iterations than the NM (although the situation is corrected if one applies a decrement length adjustment to the NM). Computational costs are also related to the order of systems to be solved, which are different in these methods (c in (2) and m − 1 in (3)). Figure 1 shows limiting examples of analytical circuit models illustrating the absolute dominance of the LM for the circuit in Figure 1a and that of the NM for the circuit in Figure 1b, since in this case each of systems (2) and (3) is reduced to a single equation. But these systems can be arbitrarily large if one applies the LM to the circuit in Figure 1b, and the NM to the circuit in Figure 1a. Therefore, in general, both methods have their place and the areas of applications where one of them proves more preferable than the other. Analysis of numerous dependencies and formulas for the description of the steady isothermal flow of the working fluid (liquid or gas) through pipelines and other elements of the PLS allows one to introduce the following classification as derived from the standpoint of their mathematical features [14]: (1) conventional (explicit); (2) implicit flow rate-based; (3) pressure-dependent flow models; (4) 2-and-3 hybrid cases.
In the conventional case where s i is the hydraulic loss of the branch, Y i > 0 is the pressure increment at active branches (e.g., with pumping stations) and Y i = 0 at passive branches (e.g., for pipeline sections). Here we have explicit expressions Implicit flow rate-based models of the flow are mainly applicable to pipelines. It follows from the basic Darcy-Weisbach formula [15] 2g that s = λ 8l ρπ 2 d 5 , since V = 4|x|/(πd 2 ρ) and y = ρgh. Here: h, V-head loss (m) and flow velocity; y, x-pressure loss and mass flow rate of the working fluid; l, d-pipe length and inside diameter; ρdensity of the transported fluid; g-gravitational acceleration. The coefficient of friction drag λ can be treated as a constant only in the region of fully-developed turbulence, while in the general case it depends on velocity (Reynolds number) Re(V) = Vd/ν, where ν is the kinematic viscosity coefficient of the fluid. Examples of such dependencies are the common Colebrook-White [16] and Altshul [17] formulas. Therefore, friction loss s i for the branch that models the pipeline will generally be a function of s i (x i ). Function ψ i (y i ) becomes implicit, and the derivatives required for the LM or NM to be applied must already be determined by the rules of differentiation of such functions, which produces [14]: Pressure-dependent flow models mainly arise from the working fluid (natural gas, water vapor, etc.) compressibility effect. Here the pressure loss depends not only on the flow rate but also on the value of the pressure itself. Let us give examples of such models, denoting by p H , p K the pressures of the working fluid at the beginning and at the end of the modeled PLS element.
The dependency of the pressure built up by a centrifugal gas compressor unit on its inlet volumetric flow rate (q H ) is represented as a graphical characteristic of the compression ratio p K p H = ε(q H ). Functions ε(q H ) are approximated with sufficient accuracy by algebraic polynomials of degree three or lower [18,19] ε(q) = a 0 + a 1 q H + a 2 q 2 H + a 3 q 3 H . Taking into account that q = x/ρ and the equation of state of the gas can be represented as ρ = p/k (where k = ZRT is a known coefficient for a given gas temperature T and its physical properties: gas constant R and compressibility coefficient Z) we arrive at the following flow model Numerous variations of this model are known [20,21], etc., they are introduced as a result of neglecting some of the terms of the polynomial or approximating ε 2 (q H ) instead of ε(q H ).
In [22], the authors present a gas pipeline section model that allows factoring in more adequately the difference in geodetic heights of its ends where θ-pipe angle; Z-compressibility factor; R-gas constant; T-average gas temperature. In [23] an alternative model is given where α is the coefficient that depends on the difference of geodetic heights of the ends of the section. When modeling the pipelines of water-steam circuits of steam-turbine plants, the following formula is used to calculate the pressure drop [24,25] where ρ( p) , ρ( p) -densities of liquid and vapor phases of the flow in the function of average pressure p; ω-cross-sectional area of the pipe; σ-coefficient that accounts for the effect of the structure of the steam-water flow on the friction loss; χ-mass vapor content of the flow. For turbine control valves the following formula applies [26] 1.09 The above examples (4)- (8) attest to the impossibility of reducing it to conventional form y = f (x), where y = p H − p K . Conceivably, the technique of double iteration cycles [1,2] is applicable to the analysis of the HC with such flow models. Here, in the inner iteration cycle, the classical flow-distribution problem is solved (by means of the LM or NM) at fixed values s * , Y * of pressure-dependent functions s(p H , p K , x), Y(p H , p K , x) of flow models y = s * x|x| − Y * , and in the outer cycle, the values of these functions are defined more precisely depending on the obtained values of flow rates and pressures. For example, in (5) one can assume The technique is sufficiently universal, but its application is associated with the ambiguity of arriving at functions s, Y, the need for special justification and ensuring convergence, a multiple increase in computational costs if compared to the conventional case.
The above analysis of the published research shows that many authors: (1) pay attention to the fundamental difference of the models of the flow of compressible fluid from that of the case of incompressible fluid [12,21,22,25,[27][28][29]; (2) note the impossibility or significant difficulties of applying conventional methods of analysis of the flow distribution in these cases [12,[21][22][23]25,[27][28][29][30][31]; (3) state the lack of general, universal methods [12,21,22,25,28,30]; (4) propose their own methods and algorithms for calculating the flow distribution in complex systems of conveying gas through trunk pipelines [12,21,22,[28][29][30][31], heat and steam supply systems [25], process pipelines systems [27] and other PLSs, transporting gases or gas-liquid mixtures. The methods proposed in the published research boil down either to the technique of double iteration cycles [25,27,28], or presuppose increasing the dimensionality of the systems of equations to be solved as compared to (2) or (3) [22,23,30]. So, for example, at each step of the solution of the flow-distribution problem by the Newton-Raphson method in [12,[21][22][23], etc., it is proposed to solve a system of linearized equations of order n + m − 1 to simultaneously find the step direction by flow rates (x) and pressures (P), and in [30] it is proposed to solve that of order n by flow rates (x).
A new method, called the modified node method (MNM), was proposed in the study [14]. It provides a solution to the problem with arbitrary flow models in the space of vector P. That is, at each step of the iterative process linear systems of equations of order m are solved. Below, for the first time, an attempt is made to generalize the problems and methods of calculating the steady-state isothermal flow distribution in hydraulic circuits for arbitrary models of the fluid flow. In particular, a new (modified) loop method (MLM) is proposed, which provides a solution to the problem in the space of loop flow rates (x C ) of the order of c.

Generalization of the Flow Distribution Model and the Node Method
All of the above types of liquid or gas flow models along the elements of the HC can be reduced to a general form of writing down pressure-dependent model ϕ(p H , p K , x) = 0. With this in mind, instead of (1), the following generalized model of flow distribution is proposed [14] Ax = Q, where p H , p K -n-dimensional vectors of pressures at the beginning and end of branches, respectively, A H , A K -(m × n)-matrices of incidences that capture separately the initial and final nodes of branches so that A H + A K = A, ϕ(p H , p K , x)-a vector-function with elements ϕ i (p H,i , p K,i , x i ), i = 1, n, reflecting arbitrary flow laws, including conventional ones, implicitly determined by flow and pressure-dependent.
Here are the main provisions of MNM for solving the problem of flow distribution based on model (9)- (12). Let at the k-th iteration there be value P k , to which we can map unequivocal values of p k H , p k K and x k from (10), (11), and (12), respectively. Linearization of Equations (9)- (12) yields where u k , ϕ x = ∂ϕ ∂x -diagonal matrices of partial derivative matrices of order n. Let us express ∆x k from (15) (14) we obtain Substituing (16) in (13) yields By the rules of differentiation of implicit functions Hence, ϕ PK and are the diagonal matrices of the corresponding partial derivatives of order n. Then system (17) can be represented as Thus, the proposed MNM is reduced to an iterative process of finding a solution in the space of nodal pressures P k+1 = P k + λ k ∆P k , where ∆P k is determined based on the solution of (18), λ k -step length as determined by the condition of one-dimensional minimization of some norm of vector u k 1 . We will also show that the considered NM modification covers the canonical case when In this case ϕ PH = E n , ϕ PK = −E n , where E n is a unit matrix of order n, and Given that ϕ x = − f x , instead of (17) we have (3). This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation as well as the experimental conclusions that can be drawn.

Conditions of Existence and Uniqueness of the Solution
Proofs of the existence and uniqueness of the solution to the conventional problem of flow distribution as based on model (1) are given in [1,32,33]. The existence and uniqueness conditions are reduced to the requirement of monotonicity of functions f i (x i ), i = 1, n. Namely, these functions must be smooth (continuously differentiable) and monotonically increasing, i.e., (Figure 2a). For the conventional model represented in the generalized form (19), these conditions are equivalent to simultaneously meeting requirements of monotone increasing ϕ i (p H,i , x i ) (for arbitrary fixed p * K,i ) and monotone decreasing ϕ i (p K,i , x i ) (for fixed p * H,i ) over the entire set of values p H,i , p K,i , x i (Figure 2b,c). Therefore, the requirements for monotonicity of ϕ i (p H,i , p K,i , x i ), i = 1, n, can be formulated as Let us consider the properties of coefficient matrix (18) where there is a strict inequality when node j (or t) is incidental to a branch that has node m with a given pressure at the other end. Thus (as in the classical case (3)) M has a non-strict diagonal predominance, is symmetric in structure, although is not symmetric in the values of its components. These circumstances, combined with step length adjustments, allow for confident convergence of the MNM [14].

Generalized Loop Flow Rates Method
To justify the method of solving the problem in the space of loop flow rates, let us introduce the following notations for the equations of model (9)-(12) Let there be a k-th approximation to the solution such that u k

The linearization of Equations (21)-(24) at the point of this approximation yields
A∆x k = 0, Decreasing the order of the linearized model (25)- (28) to the space of loop flow rates is reduced to the following steps.

1.
The exclusion of vectors ∆p k H , ∆p k K on the basis of (26) and (27) yields where

Decomposition of vectors and matrices of model (29) and (30) on the basis of belonging to chord branches (C and branches of the spanning tree (T) yields
3.
Vector ∆x k T exclusion. From (32) we have since B T T = −A −1 T A C [1,2]. Hence, instead of (34), we get

4.
Vector ∆P k exclusion. Let us express ∆P k from (36) Substituting this expression into (33) yields the desired linearized loop model Let us also show that the proposed method in the case of canonical closure relations (19) when ϕ PH = E n , ϕ PK = −E n and ϕ x = ∂ f ∂x = f x ,-coincides with the conventional LM. In this case Hence, the matrix of coefficients in (38) Let us show that in the case in question the right-hand side of (38) is −B f (x k ). Here (24) given (22) and (23) can be represented as Thus, in the case of relations (19), system (38) in the proposed method coincides with system (2) of the conventional LM.

Algorithmization of a Generalized Loop Method
A computational scheme of the proposed generalized method of loop flow rates is reduced to the following. Let a next approximation x k be given, then: (1) calculate value P k such that ϕ T (P k , x k T ) = 0; (2) calculate ϕ k C = ϕ C (P k , x k C ) and, if ||ϕ k C || < σ, where σ is a specified accuracy of residuals in loops, end the calculation; (3) build and solve a system of equations with respect to ∆x k C ; (4) calculate a new approximation x k+1 = x k + B T ∆x k C , k := k + 1, and go to step 1.
There are two main distinctions from the conventional scheme here.

1.
The conventional method of loop flow rates suggests that value P is calculated at the point of problem solution with respect to x. It is also necessary to determine P k in each iteration to find the derivatives at the point of running solution with respect to both x, and P. This operation does not pose any difficulties as P k can be calculated algorithmically by alternate application of relationships ϕ i (p k H,i , p k K,i , x k i ) = 0 along the branches of a spanning tree starting from its root with a given pressure P * m .

2.
The conventional method of loop flow rates contains finite relations for nonzero elements of the coefficient matrix K = B f x B T [1,2]. For example, K rr = ∑ i∈I r ( f x ) i r = 1, c, where I r is a set of branches that enter the loop r and K rt = ±( f x ) i , r, t = 1, c, r = t, where i is a branch that belongs simultaneously to loops r and t, and sign ± is a result of scalar product b ri b ti , i.e., it depends on the direction of i-branch with respect to the direction of a "bypass" of each loop r and t. Obviously, the construction of the coefficient matrix of the system (38) by formally performing operations of inversion, multiplication, and addition of high-dimensional matrices is associated with significant computational costs. Therefore, this issue requires special consideration.
. By virtue of (31), each i-th row of the n × (m − 1)-matrix Φ contains no more than two nonzero elements that correspond to the indices of the initial (in(i)) and finite (out(i)) nodes of the i-th branch of the hydraulic circuit and it has only one node if the branch is incident to node m. This means that and j = out(i). Given this, the elements of matrix (39) will be determined as It remains to clarify how to define the elements of Φ −1 T . Nodes and branches of the directed spanning tree of the hydraulic circuit can always be renumbered so that in(i) > out(i) = i. Then, the (m − 1) × (m − 1)-matrix Φ T will be an upper triangular one, and Accordingly, matrix Φ −1 T is also an upper triangular one, and its elements will be determined as The proof relies on the known relationship Φ T Φ −1 T = E, from which: 1.
at i = j, given (40), we have the equation at i < j, we obtain the equation, which, by virtue of (40), has only two terms Potentially, relationship (41) makes it possible to organize a recurrent scheme for calculating elements of Φ −1 T (guaranteeing that for each combination i < j, value (Φ −1 T ) H i,j in the right-hand part of (41) is already known) in two forms: column-by-column (from left to right) j = 1, . . . , m − 1, i = j, . . . , 1.
This guarantee follows from the ascending sequences for j. Nevertheless, here we suppose that (3) will be applied O = 0.5 ((m − 1) 2 + (m − 1)) times, i.e., for all elements of the upper triangle Φ −1 T , among which there can be many nonzero ones. Let us show (Φ −1 T ) ij = 0, if there is a path from node j to node i and (Φ −1 T ) ij = 0, if there is no such path. Assume that i 0 = out(i) = i. Then the recurrent formula for (Φ −1 T ) i j can be expanded as where i 1 = out(i 1 ) = in(i 0 ), i 2 = out(i 2 ) = in(i 1 ) and so on. Since in(i k ) > out(i k ) = i k the sequence of indices of branches is ascending (i 0 < i 1 < i 2 . . .), and the branches themselves belong to some path that finishes at node i 0 . At the k-th step of building such a sequence, there can be three cases: i k > j-as it follows from (41), (Φ −1 T ) ik,j = 0. Case 2 means that the path from node j = i k to node i = k 0 is found and here (Φ −1 T ) ij = 0, and in case 3, there is no such path, and (Φ −1 T ) ij = 0. From this, one can obtain a formula for calculating only nonzero elements of Φ −1 T : , where R (j,i) is a set of indices of branches that belong to the path from node j to node i. The main flaw of using this formula is a repetition of operations for included paths, since the sequence i 0 , i 1 , i 2 , . . . means that R (j,i0) ⊂ R (j,i1) ⊂ . . . ⊂ R (j,i) . In this regard, the following recurrent algorithm will be the most optimal: for each and as long as in(j) < m, perform j ← in(j) , j . This algorithm can be called "row-by-row bottom-up", since there can also be a column-by-column (left to right) calculation pattern. However, it works well with the formulas for calculating the matrix of coefficients K, if the order of summation is changed to a reverse one.

Analytical Case Study
Let us show what the main matrices for the scheme in Figure 3 look like. Figure 3. A scheme. Thick lines are spanning tree branches, and dashed arrows are chords.
As seen, the matrices A T T , Φ T , Φ −1 T are upper triangular. (Φ −1 T ) ij = 0, if there is a directed path from node j to node i, and (Φ −1 T ) ij = 0, otherwise. The structure of matrix T is identical to B T . Matrices M and K are diagonally dominant (by absolute value), and symmetrical in structure but not in values of their elements.

Numerical Case Study
Let us consider a numerical case study on the application of proposed methods to a fragment of the high-pressure gas transmission system [21], presented in Figure 4.
and values of their coefficients are indicated in Table 1.
The following technique is used to obtain the initial approximation P 0 in the MNM, which is equivalent to the initial approximation x 0 employed in MLM. Let an arbitrary value x 0 C be set, then x 0 , and P 0 is determined from ϕ T (P 0 , x 0 T ) = 0. The results of calculating the MLM by iterations are shown in Tables 2-4. In this study, x 0 C = x 0 1 , x 0 2 . As seen, the initial approximation (k = 0) is rather far from the solution. It has the flow values opposite to the orientations of branches, and negative pressure values that are beyond the solution. Nevertheless, both methods provide monotonic convergence, which is illustrated in Figure 5. Here, for MNM u 1 = log max   Table 3. Derivative values over iterations.   . As seen, the initial approximation (k=0) is rather far from the solution. It has the flow values opposite to the orientations of branches, and negative pressure values that are beyond the solution. Nevertheless, both methods provide monotonic convergence, which is illustrated in Figure 5. Here, for MNM ( )    Figure 6 shows the results of testing MNM and MLM on 100 vectors of initial approximations, which were generated by a generator of random uniformly distributed numbers: for MNM −100 ≤ P 0 j ≤ 100, j = 1, m − 1; for MLM −100 ≤ x 0 i ≤ 100, i = 1, c. Iteration stopping criteria are for MNM max 1≤j≤m−1 ∆P k j < 0.01, max 1≤j≤m−1 (u k 1 ) j < 0.01; for MLM max 1≤r≤c ϕ k C,r < 0.01. The study indicates that: (1) both methods demonstrate strong convergence; (2) the number of iterations weakly depends on the proximity of the initial approximation to the solution; (3) MLM requires on average fewer iterations than MNM. Note also that the number of iterations weakly depends on the scheme dimension; however, as indicated, the computational costs for each iteration depend on the topology of a particular scheme.

Conclusions
Analysis of the available variety of methods for calculating isothermal steady-state flow distribution in pipeline systems designed for various purposes shows their significant dependence on the specifics of models of the medium flow for individual elements. The analysis has also demonstrated the absence of methods that are equally applicable in the general case. The main types of such models are highlighted, including traditional, implicit by flow rate, and pressure-dependent ones. The research shows general conditions for their application to ensure the existence and uniqueness of a solution to the flow distribution problem.
Against the background of a brief description of the generalized MNM put forward in our previous studies, we propose a new loop method. Both methods are universal with regard to the specificity of flow models, coincide with conventional methods in the particular case of classical flow models, but each of them has its area of preferred application depending on the PLS topology.
We propose the methods for algorithmizing the generalized loop method, which provide the computational costs lower than those of the well-known technique of double iteration cycles and comparable to those of the classical MLM.