The calculation of steady state incompressible flow in large networks of pipes

Chemical Engineering

Science,

1973, Vol. 28, pp. 1857-1864.

Pergamon Press.

Printed in Great Britain

The calculation of steady state incompressible of pipes

flow in large networks

M. J. BENDING and H. P. HUTCHISON Department of Chemical Engineering, Pembroke Street, Cambridge, England (Received

5 October

1972)

Abstract--A method of linearisation is described which enables the calculation of the steady state flow in networks of pipes and pumps handling an incompressible fluid. In an Appendix a numerical method is outlined for the repeated solution of large sparse sets of linear equations.

INTRODUCTION

complex networks of pipes and pumps are a common feature of many chemical plants both for the distribution of process streams and for the provision of utilities such as cooling water. Examples of very large networks outside the chemical process industry are the petroleum, natural gas and water distribution systems. In any such system it is of interest to calculate the steady state ‘behaviour of the network under various conditions of output flowrates, pump pressure heads etc. This is a difficult problem to solve for several reasons. The equations describing the turbulent fluid flow in a pipe are non-linear. The network can be very large, perhaps comprising over 500 pipes, leading to storage difficulties. A final possible cause of difficulty is that the set of equations tend to be ill-conditioned with respect to the internal flow-rates. Slight changes in the network properties can lead to considerable changes in the flowrates. The original method for solving this type of problem was that developed by Hardy Cross [ 1, 21. This is based on Kirchoff’s Laws: (1) The algebraic sum of the flows at any pipe (or pump) junction is zero. (2) The algebraic sum of the pressure drops around any closed loop (mesh) of the network is zero. A set of flowrates which satisfy the first law is assumed and is then successively corrected for each mesh until the second law is obeyed. AlterLARGE

natively the pressures satisfying the second law can be assumed and corrected until the flows obey the first law. This algorithm suffers from the disadvantage that a suitable set of meshes must be obtained. An unsuitable choice can seriously impair convergence, which in any case cannot be guaranteed. The other published algorithm for solving pipe network problems is that using Kron’s method of diakoptics described by Gay and Middleton[3]. The method of diakoptics [4] was first developed in the 1930’s for the solution of large electrical circuits which are mathematically similar to pipe networks. The network is decomposed into a simple network usually consisting of disconnected loops and single branches. This new problem is solved and the results converted into the solution of the original network by the use of an appropriate transformation matrix. As pipe network problems are nonlinear, iteration is necessary. The major advantage of this method apart from its superior convergence properties is that user specified guesses for either pressure distribution or flowrates are unnecessary. Also it is relatively simple to make small adjustments to the network topology by the use of a different transformation matrix. Both the Hardy Cross method and the diakoptics method adopted by Gay and Middleton are essentially performance type calculations as the pressure characteristics of any pumps in the system and all input and output flowrates must be specified before the calculations can be under-

1857 CF!S Vol. 28, No. 10-F

M. J. BENDING

and H. P. HUTCHISON

taken. The method adopted here, which we call the linearisation method, is simpler in conception than either of these methods, is more general in its application, and appears to require smaller computation times. It is more general in that certain design type calculations can be undertaken. Examples of these are problems in which the performance of certain pumps is left unspecified in the input data and is calculated in the course of the solution, and problems in which input and output flowrates are determined so as to satisfy pressure specifications. THE

defined by the pressure at the nodes at each end. Similarly the pressure of an input/output is that of the node to which it is connected. From these definitions a set of equations can be derived for a network of X, pipes, X, nodes, X,,, pumps and X, input/outputs. These are given below. The subscripts k and 1 refer to the input and output nodal pressures respectively. The sum for each member i of a set Gj is indicated by &Gj.

(a) Mass balance over each nodej

METHOD

For the purposes of the linearisation method a pipe network can be considered as a set of units each of which is of one of four types: Node: This is the meeting point of at least two of the following: pipes, pumps, input/outputs. Its only specific property is the pressure. Pipe: This corresponds to a physical pipeline. It has a node at each end, one of which is called the input node, the other the output node. Flow is considered to be from the input to the output. A pipe’s properties are: (1) Diameter, (2) Equivalent length, (3) Absolute roughness, (4) Superficial fluid velocity, (5) Fluid friction factor. Pump: This may correspond to a physical pump or to any other internodal connection which cannot be characterised as a pipe. It has a node at each end and in this respect is identical to a pipe. Its only property is the volumetric throughput. Input/output: This corresponds to a physical input or output. An output is treated as an input with a negative flowrate. Inputs and outputs are considered as a single species of unit as their mathematical descriptions are identical and in some cases their particular nature may be unknown. It is attached to one node only and is defined by the volumetric flowrate. It should be noted that the pressure drop is not a property of the pumps and pipelines as it is

X, equations.

(1)

(b) Pressure drop for each pipe (i) If the flow is turbulent Pk -PI

= 4CipLi/D,(V*1V**

@a)

(ii) If the flow is laminar PI, - P1 = 32pLJD:Vi

X, equations.

(2b)

(c) Specified pressure drop for some (may be all) pumps Pk-P[=x

N,, equations

(3)

(d) Specified input (or output) flowrates Fi = x

NI equations.

(4)

(e) Additional nodal pressure specifications - sufficient to completely define the problem Pi =

X

X,, + XI-

N,, - Nr equations.

(5)

If the nodal pressures, pipe velocities, pump throughputs, and input/output flowrates (P, V, Q and F respectively) are variables the above set of equations is linear except for (2a). If an initial guess VJ”) for the velocity is available, Eq. (2a) can be rewritten as:

1858

P, - Pl = 4ClpL,/Di 1V$O)IV,‘“.

(6)

Steady state incompressible flow in large networks of pipes

This new set of equations is now linear and can be solved by any standard method e.g. Gaussian elimination. Thus new values of the pipe velocities can be obtained. The process is repeated until convergence is attained. This is, in principle, a very simple algorithm, but there are some difficulties. Firstly, the basic method outlined above converges very slowly. The nodal pressures converge fairly rapidly but the pipe velocities oscillate about their correct values. The reason for this can be seen from a simple example. The pipe pressure drop equations for turbulent flow are of the form: AP = KI I/(n)1I’(n+l)

(7)

If a value k’(/cn) is guessed: W+l)=

API(KII’(“‘l).

(8)

If the pressure drop remains constant the calculated velocity V(n+2)returns to Vcn)on the next iteration. Various relaxation methods have been tried to correct this- the most successful being the simplest, i.e. after each iteration the pipe velocities are taken to be the mean of the previous value and the calculated value. This leads to rapid convergence as can be seen from Fig. 3. Secondly, a major problem arises from the size of the set of linear equations that need to be solved each iteration. The example shown in Fig. 2 requires the repeated solution of a set of equations in 66 unknowns. Clearly this could not be readily achieved by conventional techniques. Fortunately, the set of equations is very sparse, all equations except those of type one have only two or three terms. Sparse matrix methods can therefore be used to great advantage. Various examples have been given in recent literature[%71. The method outlined in the Appendix has proved particularly successful for the repeated solution of sparse sets of linear equations in which some terms of the coefficient matrix change in value only between one solution and the next. THE FRICTION

FACTOR

EQUATION

The friction factor for turbulent flow is calculated from the Colebrook function[8]:

l/C = -2.5 Log, (0*27e/D +0*885/(RedC)). (9) This has the apparent disadvantage that the equation does not directly provide the factor from the Reynold’s Number and so iteration is necessary to achieve an accurate result. As Daniel notes, this iteration will converge under normal conditions [2]. However, it was found that it is not necessary to calculate exact friction factors each overall iteration. A pipe network problem can be solved satisfactorily with only single updating of the factors for each pipe each iteration. This saves a considerable amount of time as the Colebrook function requires both a log and a square root determination. Test examples of pipe network problems sometimes required fewer iterations with a single redetermination of the fraction factors than with accurate calculations each iteration. A problem with these types of fluid mechanics calculations is the arbitrary nature of the laminarl turbulent transition. Previous algorithms have assumed that an instantaneous change occurs at a Reynolds Number of 2000. This was found to be unsatisfactory for the linearisation method in some cases as the singularity in the pipe pressure drop equations gave rise to convergence difficulties. A more complex model was therefore adopted. Laminar flow is assumed for Reynold’s Numbers below Re = 2000, in the transition region a linear interpolation between the laminar factor at Re = 2000 (8/Re) and the turbulent friction factor at Re = 3000 is used. Above Re = 3000 turbulent flow is assumed. It must be recognised however, that whatever model is used, agreement with practice may be difficult to achieve if many pipes are in the transition range. THE PROGRAM

This algorithm has been implemented on the Cambridge University Titan (prototype Atlas II) Computer. The program consists of about 700 lines of Fortran and when compiled with enough store to calculate a network of 100 pipes, 100 nodes, 20 pumps, and 20 input-outputs uses

1859

M. J. BENDING

and H. P. HUTCHISON

12K 48 bit words. The maximum store usually available on Titan is 32K words, so considerably larger networks could be calculated without difficulty. The Data input system is very simple: (4 Integers specifying units of density, viscosity, length, pressure and flowrate W Density and viscosity (cl Number of nodes, pipes, pumps, and input-outputs (4 For each pipe -input node, output node, diameter, length and roughness For each pump - input node, output node For some pumps -pump number, pressure rise terminated by ‘0’ number 63)For each input-output-node (h) for some input-outputs-input-output number, input flowrate (-ve if output) terminated by ‘0’ (9 For each extra specification - node number, pressure. Erroneous data (e.g. specification of all input/ output flowrates) is detected if possible and an appropriate message is printed. The preparation of this data is straightforward as it can be read directly from the network flowchart. The initial guesses for the pipe velocities (1.0 mlsec) and friction factors (0.003) are generated within the program. The linear equations are set up and solved by the primary linear equation routine (see Appendix). If the matrix is singular information is printed on the console giving the undefined variables and redundant equations, the program then stops. These sets of undefined variables and redundant equations are to some extent arbitrary, depending on the route by which the set of equations are solved, but provide some useful diagnostic information. If the initial solution is successful, the results are substituted into the non-linear pipe pressure drop equations and the solution is found by the secondary linear equation routine. This process is repeated until every variable (i.e. nodal pressure, pipe velocity, pump throughput and input flowrate) differs by less than O-1 percent from its value at the previous iteration, a very tight convergence criterion. If the matrix becomes singu-

lar during this stage a message is printed on the console and the set of equations is resolved by the primary linear equation routine. After convergence has been attained the results are printed in the units specified by the user in his original data file. An example of this output is given in Fig. 1. The present program could be extended. Further data input options such as being able to specify pipe or pump throughputs could be included without any difficulty. A further extension would be to allow the specification of pump char-

[F,’

No of iterations = 10 Comp. time used = +2.9900 sec. Pressure at nodes (psi) 1 +2.3258E- 1 3 + 1.6699E1 5 +2.0497E- 1 7 +2.56668- 1 9 +2.1585E-1 11 + 1.2889E- 1 13 +1.2815E-1 15 +8.49138-2 17 +o 19 +8.7780E-2 21 +9.5189E-2 Stream velocities (ftlsec) 1 -3.7257 3 +4.3683 5 +3.9932E-1 7 +7.1063 9 +1.5616 11 -3.7257 13 +7.3089 15 +8.8070 17 - 1.2283E 1 19 -6.4628 21 -2.6524 23 +3.9729 25 -3.4937 27 +3.3622 29 -9.9013 31 + 1.7438E-2 33 +4.9245 35 -4.1768 37 +6.7889 Input flowrates (ftYsec) 1 +2*OftOO 3 +1.5000 5 -4.0000

TR TU LA TU LA TR TU TU TU TU LA TR TR TR TU LA TU TU TU

2 4 6 8 10 12 14 16 18 20 22

+1.83958-l +1.719lE-1 +2.4462E- 1 +2.2070E-1 + 1.9453E- 1 +l.l464E-1 +1.3830E-1 +7.57448-2 +6.88OOE-2 +8488lE-2 + 1.9505E- 1

2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36

+8.0943 +7.0207 -5.3133 +l.O469El -3.1057 +4.9966 -56639 -2.8015E-1 +5.8163 +3.7261 +7.5365 -3.8930 + 6.0729 +4.2374 -1.5441 -3.5134 +2.6121 +3.5303

2 4 6

TU TU TU TU TR TU TU LA TU TR TU TR TU TU LA TR LA TR

+3~5000 -1.0000 -2+lOOO

End of results

1860

Fig. 1. Typical results from the pipe network program.

Steady state incompressible

flow in large networks

acteristics. This should cause no mathematical problems but would complicate the data input system. RESULTS

The program has been tested with the network given by Gay and Middleton (see Fig. 2). It consists of 22 nodes, 38 pipes and 6 inputloutputs. All pipes are 100 ft long and 6 in. dia. and are smooth. The liquid is water with a density of 62.5 lb/ft.3 and a viscosity of 1 cp. The pressure at node 17 is set to zero. Various cases have been solved. These are given below, together with the number of iterations and the computation time required to achieve the criterion of all variables converged to 0.1 percent. Case (1) All input and output flowrates specified. This is the example given by Gay and Middleton. 6 iterations, 2.1 sec. Case (2) As (1) but with density set to 1.Olb/ ft.3 This is done so that an example of

of pipes

mixed turbulent, transition and laminar flows could be tested. 10 iterations, 2.8 sec. Case (3) As (1) but with density set to O-1 lb/ ft3. All flows are completely laminar. 2 iterations, 1.6 sec. Case (4) The density reset to 62.5, outputs 5 and 6 were left unspecified, node 11 was also set at O-0 psi. 6 iterations, 2.3 sec. Case (5) As (1) but with pipe 15 omitted. 8 iterations, 2.8 sec. Case (6) As (1) but with pipe 15 replaced by a pump of unspecified properties. The extra specification was obtained by setting the pressure at node 20 to 0.0 psi. 7 iterations, 2.1 sec. -_ . _. The indicated computation times depend to some extent on the number of jobs in the machine and so should only be regarded as approximate. Case (1) has also been solved with the initial guesses of

0 Node number Input 1 Input 2 Input 3

2

Pipe number -

2.0ftYsec 3.5 ftYsec 1.5 ftYsec

Input 4 Input 5 Input 6

Fig. 2. Gay and Middleton’s

1861

\ Input/output - 1.O ft3/sec - 4-O ftYsec - 2.0 ft+ec network.

number

M. J. BENDING

and H. P. HUTCHISON

the pipe velocity set to O-1 m/set and IO.0 mlsec. These calculations converged to give the same results in 6 and 5 iterations respectively showing that the values of these guesses have no significant effect on the rate of convergence. An interesting point to note is that the convergence does not seem to depend greatly on the network topology but only on the types of flow existing. Thus it can be stated that usually: Laminar flow problems converge in 2 iterations. Mixed flow problems converge in lo-13 iterations . Turbulent flow problems converge in 5-8 iterations. This has been confirmed by a solution of the network reported by Doland and subsequently investigated by Ingels and Powers[9]. The results obtained by the linearisation method were in good agreement with those given by Ingels and Powers taking into consideration the different friction factor relationships. Only 5 iterations and a computation time of 2.15 set were required for this problem which was completely turbulent except for three transitional flows. This inci-

dently reinforces Gay and Middleton’s suggestion that the flowrates used by Ingels and Powers were too small. The network was also solved with flowrates of l/60 their correct value. This created a mixed laminar, transitional and turbulent problem which converged after 13 iterations and 2.95 set computation time. It should be noted that the extra number of iterations required for mixed flow problems is not of great practical significance as all iterations after the first take comparatively little computation time. A COMPARISON

WITH

OTHER

METHODS

The only direct comparison of the linearisation method with others is that for the solution of Case 1 above. The results obtained agree very closely with those given by Gay and Middleton from their diakoptic method, assuming that their pressure drops are expressed in Psf as they do not specify the units. A comparison of the convergence of the two methods can be made from the table of values of selected pipe pressure drops after each iteration given in Fig. 3. As can

A. Linearisation method. (i) Pine velocities in ftlsec- initial guess of 3.28 ft/sec. Iteration Pipe2 - Pipe 5 Pipe 17 1 7,270 0.3665 13.61 2 7,474 04611 12.28 3 7.504 04682 12.11 4 7.505 04678 12.11

Pipe 28 3.051 4.152 4.276 4,273

Pipe 34 2.339 2.563 2.615 2.617

(ii) Pipe pressure drop in lbf/ft* Iteration Pine 2 1 70.71 2 146.8 3 149.1 4 148.9

Pipe 28 29.70 50.82 53.73 53.78

Pipe 34 22.78 21.98 22.23 22.27

Pipe 5 3.568 1.116 1.031 1.047

Pipe 17 13-25 345.8 355.4 355.0

B. Diakoptic method as given by Gay and Middleton Pipe pressure drop inibf/ft2 Pipe 2 Pipe 5 Pipe 17 Iteration 4.85 176 1 96 2.00 261 2 121 1.37 309 3 135 4 333 142 l-17 1.10 344 5 145 147 1.05 6 350 7 1.04 352 148 8 148 1.04 353

Pipe 28 40 47 52 53 53 54 54 54

Fig. 3. Results from the solution of problem case 1.

1862

Pipe 34 30 26 24 23 22.5 23 22 22

Steady state incompressible

flow in large networks of pipes

be seen the linearisation method converges faster. Gay and Middleton found that the convergence they obtained was better than that obtained from Daniel’s Hardy Cross method[2]. A comparison of the computation time required is impossible as they give no information about their computer. An advantage claimed for the diakoptics approach is that the effect of changes in the network can be determined relatively easily. This is clearly important for Gay and Middleton’s calculations because running times of the order of 20 min were required for a complete solution. This facility is not possible with the linearisation method as each change in the network topology or problem specification necessitates a complete recalculation. However as computation times of the order of only 3 set are required for this it does not seem to be a significant disadvantage. A major advantage of the linearisation method is that it is very versatile although its data input requirements are simple. The program merely requires the topology of the network and sufficient extra specifications to define the problem. This may be contrasted with both the Hardy Cross and diakoptics approaches. The Hardy Cross method requires the mesh specification whilst the Gay and Middleton method requires two sets of incidence matrices defining the decomposition of the network. Neither of these is immediately obvious from the network diagram. We conclude that the linearisation method, as we have developed it, is a satisfactory general method for the solution of the steady state

behaviour of pipe networks incompressible fluid.

filled with a single

Acknowledgement-The tenure of an S.R.C. Research Studentship by M. J. Bending during the course of this work is acknowledged. NOTATION

A

c D F

Gj Hj lj Li K

k L 1

NI N PU

L R

s V

XI XN XP X PU x P

P

pipe cross-sectional area friction factor pipe diameter input volumetric flowrate the set of pipes for which node j is the input the set of pumps for which node j is the input the set of input-outputs to node j pipe, node, pump or input-output number a constant input node number pipe equivalent length output node number number of input-output definition equations number of pump definition equations nodal pressure pump volumetric throughput the set of pipes for which node jis the output the set of pumps for which node j is the output superficial fluid velocity in pipe number of input-outputs in network number of nodes in network number of pipes in network number of pumps in network a specified constant fluid viscosity fluid density

REFERENCES [l] STREETER V ., Handbook ofFluidDynamics, pp. 3-27. McGraw-Hill, New York 1961. 121 DANIEL P. T.. Trans. Insm Chem. Engrs 1966 44 T77. ,3] GAY B. and MIDDLETON P., Chem.Engng Sci. 197126 109. [4] KRON G., Diukoptics. Macdonald, London 1963. [5] WILLOUGHBY R. A., Proc. Symp. on Sparse Matrices and heir Applications. IBM Thomas J. Watson Research Centre 1969. [6] REID J. K., Large Sparse Sets ofLinear Equations. Academic Press New York 197 1. [7] ROSE D. J. and WILLOUGHBY R. A., Sparse Matrices undrheir Applications. Plenum Press 1972. [8] COULSON J. M. and RICHARDSON J. F., ChemicalEngineering, Vol. 1, p. 49. Pergamon 1964. [9] INGELS D. M. AND POWERS J. E., Chem. Engng Prog. 1964 60 65.

1863

M. J. BENDING

and H. P. HUTCHISON

APPENDIX -1 A METHOD FOR THE REPEATED SOLUTION OF SPARSE SETS OF LINEAR EQUATIONS We describe the method in outline only and do not consider features that permit the recognition of singularity as this would complicate the description. In this description we presume that a rearrangement of rows and columns has been achieved so that an element of reasonable magnitude is in each diagonal position.

Mathematically, stages.

the method is Gaussian elimination in two

Triangularisation steps are: (1) Eliminate a,; create a new element a,, = -a5 . a,/~?, (2) Eliminate a,,; a, becomes a: = a, -a5 . a,,/~, (3) Eliminate a,; a: becomes a: * = a: - a, . as/us. This process can be characterised by a string of integers terminated by two zeros: -2

3

12 5 -11

3

The resulting equations are:

1 5 -9

8

1 4 0

0

Backsubstitution steps are: (1) x, =-1 (2) xi = -a, . x,/u, (3) .r, = (-a,~~ - a12.r1)laa (4) x3 = -a4xx,lu, (5) x, = -a,x,/u,. This can be characterised by a string of integers: 7 0 -1

1

10 0

12

1 -2

6

4

1 -3

8

5 -4

1 3

in which the negative sign is used to indicate the final calculation of each unknown. The right hand side of the set of equations is considered as the coefficients of an extra unknown with a value of -1 as this simplifies both the list structure and the linear equation routines. The combined lists of integers becomes a specification or ‘Operator List’ for the complete solution of any set of equations with non-zero elements arranged in the same matrix positions as in the original equations. Two routines have been written. The primary linear equation routine solves the set of equations for the first time and also stores the operator list for the solution process. To do this it requires three lists of numbers giving respectively the column number, row number and value of each non-zero element. There is no need to prearrange the equations so that the matrix has non-zero diagonal elements as the routine does this itself and also rearranges the matrix so that the number of created elements (e.g. a,2 above) is kept to a minimum. The secondary linear equation routine solves the equations using only the operator list and the (updated) list of element values. This operator list method can save a significant amount of time if the matrix is sparse as the secondary linear equation routine does not need to search for elements which probably do not exist. Large test matrices (say 300 equations, 1000 elements) have been solved by the secondary linear equation routine in l/40 the time required for the primary solution.

1864