UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository   Help 
Material Information
Notes
Record Information

Full Text 
Publication No. 77 PIPE NETWORK ANALYSIS By MunFong Lee University of Florida Gainesville ''' ''* I:i " . ,i ', , ~% PIPE ;:','a'ORK ANALYSIS BY MUNFONG LEE NONTHESIS PROJECT REPORT IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE MASTER OF ENGINEERING DEGREE DEPARTMENT OF ENVIRONMENTAL ENGINEERING SCIENCES UNIVERSITY OF FLORIDA NOVEMBER 1983 ABSTRACT Analyses of large scale pipe networks are needed whenever significant changes in patterns or magnitudes of demands or supplies occur in municipal water or gas distribution systems. Changes of this nature occur whenever new industrial and residen tial areas are being developed or new sources of supply are tapped. In the absence of such analytical tools to determine the perform ance of an existing system under new demands, needlessly large investments are made for larger than necessary pipes, redundant lines or duplicate facilities. Another cause for concern is the ability of the numerous algorithms to provide reliable results without which deficient engineering judgments may be made in engineering applications deal ing with large scale pipe networks. Convergence and reliability problems of most of the algorithms are highlighted after the theoretical background has been presented. As an aid to more effective formulation of the loop and nodal equations, the essential concepts of network theory are also presented together with the fundamental hydraulic principles forming the backbone of the state of the art iterative procedures. This report concludes with a new approach which employs opti mization techniques to solve the pipe network problem as a viable and perhaps more versatile alternative to the widely used iterative methods. ACKNOWLEDGE EI l.i 1 ' I wish to express my appreciation to my graduate committee chairman, Dr. James P. Heaney, for his advice and enthusiasm in directing this report and my general course of study at the University of Florida, Special thanks are due to the members of my committee, Dr. Wayne C. Huber and Dr. Warren Viessman, Jr., for their assistance and review of this report. I wish also to record my grateful thanks to my employer, the Public Utilities Board, Republic of Singapore, for providing me this opportunity to further my studies. 97K F S5'r i K TABLE OF CONTENTS ABSTRACT ...................... ... ............. ........ ACKNOWLEDGEMENT ......... ................. ................. TABLE OF CONTENTS .......................................... PAGE i ii iii CHAPTER 1 1.1 1.2 1.3 CHAPTER 2 2.1 2.2 2.3 2. I 2,5 2,6 2.7 2,8 CHAPTER 3 3.1 3.2 CHAPTER 4 4.1 4,1.1 4.1.1.1 4.1.1.2 S INTRODUCTION ...... ... ................. ..... Problem Definition ......................... .... Significance ........ .... .. ................ Motivation .......... ... ...... ......... ........ NETWORK ANALYSIS & PIPE NETWORK HYDRAULICS ..... Fundamentals of Network Theory ................ Pipe Network Conservation Laws ............... Friction and Minor Losses ..................... Pumps .. .. ... ....... ....... .................... . Pressure Regulating Valves ................... Node Analysis ....... ........ ....... ............ Loop Analysis ....................... .......... Corrective Mesh Flow Analysis ................ SNEWTONRAPHSON METHOD ........ .............. Application to Node Equations ...... .......... Application to Corrective Mesh Flow Equations .. : LINEAR METHODS .... ...*........... .. .. .. ........ Gradient Method ......................... ....... Algorithms for the Solution of Loop Equations . Single Path Adjustment (P) Method ............ Simultaneous Path Adjustment (SP) Method ..... iii 4.1.1.3 Wood's Linear (L) Method ....................... 4,1.2 Algorithms for Solving Node Equations ........... 4.1.2.1 Single Node Adjustment (N) Method .............. 4.1.2,2 Simultaneous Node Adjustment (SN) Method ....... 4,2 Comments on Algorithms using Gradient Method .... 4.2.1 Accuracy of Solutions ........ .....,... .......... 4.2.2 Reliability of Algorithms ....... ............... 4,3 Linear Theory Method by Wood and Charles ........ 4.3.1 Inclusion of Pumps and Reservoirs ............... 4.3.2 Inclusion of Pressure Regulating Valves ......... CHAPTER 5 :ALTERNATIVE MATlHE.i ..TICAL APPROACHES ............. & COMPUTATIONAL EXPERIENCE 5.1 Mathematical Programming Techniques .............. 5.2 Computational Experience ..................... CHAPTER 6 : CONCLUSION ....................................... REFERENCES ... ... .... ............... APPENDIX .......................... .. PAGE 26 27 27 28 31 31 31 33 34 36 37 37 41 45 48 53 CHAPTER 1 1 INTRODUCTION Analyses and design of pipe networks create a relatively complex problem, particularly if the network consists of a range of pipes as frequently occurs in water distribution systems of large metropolitan areas, or natural gas pipe networks. In the absence of significant fluid acceleration, the behavior of a network can be determined by a sequence of steady state conditions, which form a small but vital component for assessing the adequacy of a network. Such an analysis is needed each time changing patterns of consumption or delivery are significant or addon features, such as supplying new subdivisions, addition of booster pumps, pressure regulating vales, or storage tanks change the system. The steady state flows of a network are governed by the laws of conservation of energy and mass and the classical pipe network analysis problem is to establish the steady state flows and pressures in a full flow closed conduit network of known physical characteris tics, Due to the complexity and the inherent nonlinearity of net works, solving the network analysis problem is not a trivial exercise. For over four decades, a number of algorithms havebeen developed since the pioneering work of Hardy Cross. All of these techniques are iterative in nature, differing only in the method in which an estimate of the true solution is obtained, A recent study (Collins, Cooper, Helgason and Kennington, 1978) uncovered a new approach to the pipe network analysis problem using optimization techniques which represent a radical departure from the traditional state of the art methods. This report attempts to provide a comprehensive write up of the theory behind some of the more commonly used algorithms and their efficiency and reliability. 1.1 PROBLEM DEFINITION A pipe network is physically a collection of interconnected elements such as pipes, pumps, reservoirs, valves, and similar appurtenances. Mathematically, the network is represented as an edge set consisting of pipes, pumps, valves and similar elements and a node set comprising reservoirs and element intersections. In most of the elements, a unique functional relationship between pressure and flow exists. Pressure, in incompressible flow net works, can be expressed in terms of an equivalent hydraulic head, a terminology which will be adopted throughout this report as is standard practice. The steady state condition of a network can be completely defined by the head at each node and the flow in each element. Having determined this unique set of flows/ heads for a given set of inputs and withdrawals, all other quantities of interest can be deduced therefrom. 1,2 SIGNIFICANCE Steady state network analysis is a basic tool in water dis tribution system management and design. It can also be used to develop operating policies and strategies to not only reduce operating costs but also increase reliability and reduce water wastage (Brock, 1970; Hudson, 1974; Rao et al. 1974, 1977; Shamir, 1974; Bree et al. 1975). Application of steady state network analysis in online system control is also receiving growing attention (Brock, 1963; Hudson, 1973; McPherson et al, 1974; Rao et al. 1974; Gerlt and Haddix, 1975; Eggener and Polkowski, 1976). 1.3 MOTIVATION Since Hardy Cross first provided a solution for the pipe network analysis problem, three general methods which are widely used today, have evolved: (i) Hardy Cross (Hoag and Weinberg, 1957; Graves and Branscome, 1958; Adams, 1961; Brock, 1963; Bellamy, 1965; Dillingham, 1967; Fietz, 1973; Williams, 1973; Chenoweth and Crawford, 1974; Eggener and Polkowski, 1976) (ii) NewtonRaphson (Martin and Peters, 1963; Shamir and Howard, 1968; Liu, 1969; Epp and Fowler, 1970; Zarghamee, 1971; Lam and Wolla, 1972; Lemieux, 1972; Donachie, 1973; Rao et al, 1974,1977) (iii) Linearization (McIlroy, 1949; Marlow et al. 1966; Wood, and Charles, 1972; Fietz, 1973; Collins and Johnson, 1975). These methods solve a set of nonlinear simultaneous equations iteratively beginning with an initial trial solution. The iteration is complete when a new solution differs from the trial solution by less than a specified amount; otherwise, the new solution becomes the trial solution and the procedure is repeated. Differences in the above methods arise because of the strategies used to determine a new solution. In view of the iterative nature of these methods, large scale networks with hundreds of nodes and elements require considerable computer efforts to solve. The choice of algorithm therefore, depends on the computational speed and reliability of a particular solution procedure. Matrices associated with water distribution networks, like most manmade systems, are sparse. One of the keys to faster conver gence and hence to greater computational efficiency and perhaps reliability for most, if not all, algorithms is the use of sparse matrix techniques in the solution procedures (Tewarson, 1973). In the following chapters, most of the essential tools.required for the analysis of incompressible flow in pipe networks are presented. Chapter 2 introduces graph theory which is useful in the formulation of pipe network simulators and also includes fundamental hydraulic principles governing pipe networks to provide the neces sary groundwork for the development of the loop and node system of equations. In Chapters 3 and 4, methods for solving these systems of nonlinear equations are described. Alternative mathematical approaches and the writer's computational experience are presented in Chapter 5. * :' *S * CHAPTER 2 NETWORK ANALYSIS & PIPE NETWORK HYDRAULICS There has been growing awareness that certain concepts and tools of network theory are very useful in the analysis of pipe networks especially in the formulation of computer simulators. The theory of network analysis is well established and several refer ences in this field are available (Gulliman, 1953; Belevitch, 1968; Karni, 1971; Clay, 1971; Shamir, 1973; Bazaraa,and Jarvis, 1977; Minieka, 1978), For consistency, the terminology used in this chapter has been adopted for pipe networks. Also presented in this chapter are some of the fundamental hydraulic principles which form the foundation of the three tradi tional methods described in Chapter 3. 2.1 FT'.iiI1;.MIENTALS OF NETWORK THEORY According to network terminology, a network is a graph consis ting of a set of junction points called nodes, with certain pairs of nodes being joined by line segments called edges (or arcs, branches or links). Edges joining the same two nodes are multiple edges and a node without an edge connected to it is an isolated node. If a node has only one edge connected to it, that edge is a pendant. An edge and a node at the end of the edge are said to be in_.dent. A subgraph is any collection of nodes and edges comprising only nodes and edges of a larger graph. The complement of a subgrap is the collection of nodes and edges remaining after the removal of the subgraph. A path between two nodes is a subgraph whose terminal nodes each have only one arc incident and all other nodes are incident to exactly two arcs. A graph is said to be connected if there is a path connecting every pair of nodes. A connected subgraph in which each node of the subgraph is incident to exactly two arcs of the subgraph is called a loon (or cycle). A tree is a connected graph containing no loops. The comp lement of a tree is a cotree. Edges of a cotree are links. A tree containing all nodes of a graph is a spanning tree. An edge of a graph is said to be directed (or oriented) if there is a sense of direction ascribed to the edge. If all edges of a graph are directed, it is called a directed graph. However, a network need not be directed because it may be feasible to have flow in either direction along an edge. The flow capacity of an edge in a specified direction is the upper limit to the feasible magnitude of the rate of flow in the edge in that direction. The flow capacity may be any nonnegative quantity, including infinity. An edge is directed if the flow capacity is zero in one direction. The topology of a directed graph of i nodes and E edges can be described by a ) X F node incidence matrix, A with typical element ( + 1, if edge j is directed away from node i ) ( ) ai ( 1, if edge j is directed towards node i ) ij( ) ( 0, if edge j is not connected to node i ) For a connected graph it is apparent each column of A will contain a 1 and a 1 and all remaining elements will be zero. As a check, addition of the rows of A should yield a zero row. Thus, the rank, r, of A is at most ) 1. If loops are formed, one by one, by adding links, one at a time, to a given spanning tree, it is apparent that each time a link is added a unique loop will be created. Such a loop is called a fundamental loop. A fundamental loop set for any connected graph, containing A loops, can be described by a A X C fundamental loop matrix, B, with typical element ( +1, if edge j is in loop i and the direction of edge j ) ( is clockwise, say ( ) bi ( 1, if edge j is in loop i and the direction of edge j ) 10 ( ) ( is counterclockwise ) ( ) ( 0, if edge j is not in loop i By performing elementary row operations on B, an identity submatrix of order, A can be obtained, implying the rank of B is A . Both the node incidence matrix and the fundamental loop matrix can be used to formulate the continuity and energy (or loop) sets of equations in a computer simulator. 2.2 PIPE NETWORK CONSERVATION LAWS Pipe network parameters are introduced to develop two conser vation laws utilizing graph theory. The following notation shall be adopted for convenience. A directed network will be described by a node set, N and an edge set, E of ordered pairs of nodes. Each node n C N is associated with a unique number called the head, Hn. For an edge directed from node i to node j, an edge head loss is defined as AHiji a Hi Hj in which AHHij = Hji. In each edge, a flow Qij exists, positive when the edge is directed from node i to node j. A basic law to be satisfied by the flows in a network is mass conservation, Q n Qin = r all nEN (2.1) (n,j)e E (i,j)$ E where rn is the requirement at node n, positive for inflows (supply) and negative for outflows (demand). Denoting the vector of Qij's by q, equation (2.1) can be rewritten as A q (2.2) where r = (rl, r2, .... r,). As noted in section 2.1, A has a rank of 9 1, implying one of the rows in A is redundant and can be arbitrarily omitted. The matrix Ar, obtained by deleting one row of A, say row '7, is defined as the reduced node incidence matrix, A corresponding element in T is also deleted and a demand vector, d, defined as d s (rl, r2,.....,r_ 1) is introduced. Then Ar d (2.3) It should be noted, in passing, that all the rows in A will be independent, that is, rank of A = if pumps and reservoirs are present in the network. However, a redundant row still exists if a junction is assumed at the reservoir or pumps. If the nodal head is unique, as assumed, then the summation of head losses around a loop is zero. This obvious proposition is used as the basis for the second fundamental network law. Thus, if Lk is the edge set for edges in fundamental loop k, k = 1, 2, ... A then SAH 0 all k (i,j)ELk Equation (2.4) can be written as BAh 0 (2.5) ifAh is defined as the vector of Hij's. If a mesh flow vector p = (Pl, P2, P3 *.... ) is defined, the following relationship can be written T BT q = B p (2.6) Thus if to each fundamental loop a unique mesh flow is associated, the flow in any edge is a linear combination of the mesh flows for fundamental loops containing the edge in question. If pumps and reservoirs are included in the pipe network, equation (2.5) is generalized as follows: AE = BAh h (2.7) where h = pump head vector, AE = vector representing the diff P erence in total grade between two reservoirs. In this generalized case, a junction node is assumed at a reservoir or pump and a pseudo loop is assumed to connect 2 reser voirs. 2.3 FRICTION & MINOR LOSSES The relation between head and discharge, that is, zh and q completes the number of equation sets required to define the net work problem. Total head loss in a pipe, H, is the sum of the line loss, hLp, and minor loss, hLM. The line loss expressed in terms of the discharge is given by: hLP = KpQn (2.8) where Kp is a pipe constant which is a function of line length, diameter and roughness and n is an exponent. Commonly used head loss equations include the DarcyWeisbach, HazenWilliams and Manning equations. Perhaps the most widely used of these equations is the HazenWilliams equation, which is, English System (ES) Q = 1.318 CHW AR"63 S054 ) 0.63 o.4 ) (2.9) S.I. Units Q = 0.849 CHW AR 63 S.5 ) in which CHW is the HazenWilliams roughness coefficient, S is the slope of the energy line and equals hLP/L, R is the hydraulic radius defined as the crosssectional area, A, divided by the wetted perimeter, P, and for full pipes equals D/4 (where D = diameter of pipe). Table 2.1 gives values for CHW for some common materials used for pressure conduits (Jeppson, 1977). Type of Pipe CHW PVC pipe 150 Very smooth pipe 140 New cast iron or welded steel 130 Wood, concrete 120 Clay, new riveted steel 110 Old cast iron, brick 100 Badly corroded cast iron or steel 80 Table 2.1 : Values of HazenWilliams Coefficient Equations (2.9) can be written in terms of hLP if Q is known. Thus, ES hLP 8.52 X 105L Q1852 1.852 I4,87 C D Hw with D in inches and L in feet 10.7 L Q1.852 SI hLP C1.852 DL 7 SI hyp 7W  HW Principles governing the flow of fluid as well as much exper im6ntal evidence indicates that the head loss due to added turbu lence or secondary flow in the presence of fittings, valves, meters and other components in a network, will be approximately proportional to the square of the velocity or the flow rate squared. Minor losses are commonly expressed in the form hL KM 2 (2.10) in which KM = M/(2gA2). Nominal values of M for various common al.pu,!i.enances are given in Table 2.2 (Jeppson, 1977). It is apparent from these loss coeffi cients that minor losses can be neglected if relatively long pipe lines are analyzed. However, in short pipelines, they may represent the major losses in the system, or if a valve is partly closed, its presence has profound influence on the flow rate. 2.4 PUMPS A number of alternative methods might be used to quantify the head, hp, produced by a pump. In some cases a constant power input is specified. In general, the relationship between pump head, hp and discharge, Q, can be expressed as hp P(Q) (2.11) For a constant power pump, P(Q) = Z/Q (2.12) DEVICE M Globe Valve (fully open) 10 Angle Valve (fully open) 5 Gate Valve (fully open) 0.19 Gate Valve (3/4 open) 1.0 Gate Valve (1/2 open) 5.6 Ball Check Valve (fully open) 70 Foot Valve (fully open) 15 Swing Check Valve (fully open) 2.3 Close Return Bend 2.2 Tee, Through Side Outlet 1.8 Standard Short Radius Elbow 0.9 Medium Sweep Elbow 0.8 Long Sweep Elbow 0.6 450 Elbow 0.4 Table 2.2 : Loss Coefficients for Valves and Other Pipe Fittings where the pump constant, Z = 550 HP/62.4 and HP = useful pump horsepower. Other functions have been suggested, and a common choice is a second order polynomial of the form P(Q) = AQ2 + BQ + H (2.13) in which A, B and Ho are constants for a given pump and might be determined by fitting Equation (2.13) to three points taken from a pump characteristic curve. 2.5 PRESSURE REGULATING VALVES A pressure regulating valve (abbreviated PRV) is designed to maintain a constant.pressure downstream from it regardless of how large the upstream pressure is. Therefore, it is apparent that the unique relationship that exists between head and discharge for line losses, minor losses and pumps, does not exist for a PRV. Solution of pipe networks which include control elements with non unique head discharge relationships using optimization techniques is still an active research area (Collins, Cooper, Helgason, Kennington, 1978). Exceptions to the above occurrence are; (1) If the upstream pressure becomes less than the valve setting, and (2) if the down stream pressure exceeds the pressure setting of the valve so that if the PRV were not present, the flow would be in the opposite direction to the downstream flow direction of the valve. If the first condition occurs, the valve has no effect on flow conditions. The PRV acts as a check valve, preventing reverse flow if the second condition occurs. By preventing reverse flow, the PRV allows the pressure immediately downstream from the valve to exceed its pressure setting. Thus, PRV's can perform both functions of reducing pressures in portions of a pipe distribution system if the pressures would otherwise be excessive, and may also be used to control from which sources of supply the flow comes under various demand levels. In the latter application, the PRV acts as a check valve until the pressure is reduced to critical levels by large demands at which time additional sources of supply are drawn upon. The analysis of a pipe network containing PRV's must be capable of determining which of these conditions exist. 2.6 NODE ANALYSIS To obtain the system of equations which contains the heads at the junctions/nodes of the network as unknowns, the 7 1 independ ent continuity equations are written as in Equation (2.3). The relationship between discharge and head loss is then substituted into the continuity equations yielding a set of r) 1 equations in 1 1 unknown nodal heads. Solving for Q from the exponential formula (Equartion 2.8), using double subscript notation, gives Qij (AHij/Kij)l/n (2.14) in which AZHij : (hLp)ij + (hLM)ij and Kij (Kp)j + (KM)ij Substituting Equation (2.14) into the junction continuity equations gives ( H H )1/n ( H H. ) /n L Kij ) o HK ) Sout iJ (2.15) 2.7 LOOP ANALYSIS If the discharge in each pipe is initially considered unknown instead of the head at each junction, the number of simultaneous equations to be solved is increased from (Y 1) to ( 1 +X) equations. However, this increase in the number of equations is somewhat compensated by a reduction in the number of nonlinear equations in the system. The analysis of flow in networks of pipes is based on the energy and mass conservation laws discussed in section 2.2. Math ematically, the continuity equations are concisely expressed as: Ar d where Ar is the reduced node incidence matrix. It is apparent that each of these continuity equations is linear. The remaining set of equations is formed by applying the energy conservation principle and expressed in terms of the funda mental loop matrix, B, as follows: BAh = 0 which has A independent nonlinear equations. Having solved the system of equations for the discharge in each pipe, the head losses in each pipe can be determined. From a known head or pressure at one junction, the heads and pressures at each junction throughout the network, or at any point along a pipe, can be determined by subtracting the head loss from the head at the upstream junction, and accounting for differences in elevations if this be the case. In some problems the external flows may not be known. Rather the supply of water may be from reservoirs and/or pumps. The amount of flow from these individual sources will not only depend on demands, but also will depend upon the head losses throughout the system. 2.8 CORRECTIVE !!i,: FLOW ANALYSIS This method of analysis yields the least number of equations. However, like the node analysis method, all the equations are non linear. These equations consider a corrective mesh flow as the unknowns and as discussed in section 2.2, the system of equations to solve is written as: T  q B p in which T is the mesh flow vector. Since there are A fundamental loops in a network, the corrective mesh flow system of equations consists of 2 equations. This method requires an initialization of the flow in each pipe which satisfies all junction continuity equations. Since these initial flow estimates generally will not simultaneously satisfy the A head loss equations, they must be corrected before they equal the true flow rates in the pipes. A flow rate adjustment can be added with due regard for sign, to the initially assumed flow in each pipe forming a loop of the network without violating continuity at the junctions. This fact suggests'establishing energy equations around the a loops of the network in which the initial flow plus the corrective mesh flow rate is used as the true flow rate in the energy equations. Upon satisfying these energy equations by finding the appropriate corrective mesh flow rates, 17 the 7 1 continuity equations would remain satisfied as they initially were. The corrective mesh flow rates may be arbitra rily taken positive in the clockwise or counterclockwise direction, but the sign convention must be consistent around any particular loop. ; Si 'r * CHAP "h 3 NEWTONRAPHSON METHOD The NewtonRaphson method is an iterative scheme which starts with an estimate to the solution and repeatedly computes better estimates. Unlike other methods which converge linearly, it has "quadratic convergence". Generally if quadratic convergence occurs, fewer iterations are needed to obtain the solution with a given precision than if linear convergence occurs. In addition to rapid convergence, the NewtonRaphson method is easily incorporated into a computer algorithm. Any of the three sets of equations defining the pipe network problem, that is equations considering (1) the flow rate in each pipe unknown, (2) the head at each junction unknown and (3) the corrective mesh flow rate around each loop unknown, may be solved by this method. An initial guess is required for the NewtonRaphson method. It is the best method to use for larger systems of equations because it requires'less computer storage for a given number of equations. 3.1 APPLICATION TO NODE EQUATIONS The iterative NewtonRaphson formula for a system of equations is, (me1) (m) 1 (m) x = x D .F(x (1) in which the superscripts.within parentheses are not exponents but denote number of iterations. The unknown vectors x and F replace the single variable x and function F and the inverse of the Jacob ian, D replaces l/t in the formula for solving a single equation. Adapting Equation (3.1) to solving the set of equations with the heads as unknowns, Equation (3.1) becomes _(m+l) (m) ,1 (m) H H D F(H ) (3.2) Making up the Jacobian matrix D, are individual rows consist ing of derivatives of that particular function with respect to the variables making up the column headings. For the system of head equations, the Jacobian is, a8F1 F 9PF aH1 ,H aH2 ......... aHj 9 OFaF F2 I HI H2 ......... 9Hj D > a aH1 j 9H2 *****..... 8H where J number of junction nodes. The Jacobian is a symmetric matrix and an algorithm for solving a linear system of equations with a symmetric matrix may be preferred for greater computational efficiency. 3.2 APPLICATION TO CORRECTIVE MESH FLOW EQUATIONS The NewtonRaphson method when applied to this set of equations becomes p(m+) p(m) D1 F(p(m)) (33) in which the Jacobian is F]2 9F2 9F2 9P1 1P2 ...**. aPL aPl 31P2 ...... aPL where L = number of loops and P = corrective mesh flow for each loop. The NewtonRaphson method suffers from a setback of requir ing a reasonably accurate initializations otherwise it may not converge. When PRV's are present in a pipe network, the procedure of using identical loops for the corrective flow rates and energy equations must be altered. The reasons are (1) the head drop across a PRV cannot be expressed as a. function of the P's circulating through that pipe, (2) continuity at some junctions will not be satisfied if the P's are assumed to circulate through pseudo loops from artificial reservoirs created by the PRV's to another reser voir in the network, The reason is that P in a pseudo loop would extract fluid from a junction, but not add an equal flow through another pipe joining at that junction. Consequently, some ofthe loops around which the energy equa tions are written cannot correspond to the loops around which the corrective flow rates, P, circulate. The individual P's will thus be assumed to circulate around the real loops satisfying continuity at all junctions. The energy equations will be written around loops containing pipes or other elements such as pumps or reservoirs whose head losses are functions of the discharge through them. r a *i *; CHAPTER 4 LINEAR METHODS Nonlinearity of the function relating head and discharge is the crux of the problem in solving a pipe network. Recall that in the loop analysis, there are n + 1 equations of which . number of equations describing energy conservation around loops, are nonlinear. The other two analyses, namely the node analysis and corrective mesh flow rate analysis, each of which having 2. energy equations written around each loop of the network; both involved nonlinear equations in each of its entire system of equations, Chapter 3 has dealt with the straightforward applica tion of the NewtonRaphson Method to linearizing the nonlinear equations associated with the latter two methods. This chapter will be devoted to other linearization techniques, some of which are variations of the NewtonRaphson Method. 4.1 GRADIENT METHOD The gradient method, which is given extensive coverage by Wood (Wood, 1981), is derived from the first two terms of the Taylor series expansion. Any function, f(x), which is continuous, that is, differentiable, can be approximated as follows: f(x) f(xo,) 4 f'(xo)(xxo) (4.1) It is apparent from the right hand side of Equation (4.1) that the approximation has reduced f(x) to a linear form. However, if f is a function of more than one variable, Equation (4.1) can be generalized as follows: f(x(l),x(2),...) = f(x(l)o,x(2),...) + (x( ) x( ) + *9a' T (x(2) x(2)o) .. (4.2) in which the partial derivatives are evaluated at some x(l) x(l)o, x(2) = x(2)o, etc. 4.1.1 AIGORITHMVIS FOR THE SOLUTION OF LOOP EQUATIONS To conform to the notation used in this chapter, Equation (2.3) which neatly describes the mass conservation (continuity) equation for each of the j* nodes in the network, is rewritten as follows: ZQout Qin e (j equations) (4.3) in which Qe denotes the external inflow or demand at the junction node, positive for inflows. The energy conservation equation in Equation (2.5) for fund amental loops without pumps, is now rewritten to include pumps as follows: ShL = Zhp (. equations) (4.4) where hL = energy loss in each pipe, including minor losses; hp en.:., input by pumps; t = number of fundamental loops. For any two fixed grade (or reservoir) nodes, the energy conservation equation written around this pseudo loop is written as: AE = jhL Zhp (f1 equations) (4.5) in whichAE difference in total grade between two fixed grade nodes; f = number of pseudo loops. If p equals the number of pipes in the network, then the mass and energy equations form a set of p simultaneous equations of which (A + f 1) equations constituting the set of energy equations are' nonlinear. *for networks with reservoirs Using Equations (2.8), (2.10), (2.11) and (4.5), the energy equations expressed in terms of the discharge, Q, are AE = E(KpQn 4 KQ2) P(Q) (4.6) It can be seen that Equation (4.4) is a special case of Equation (4.6) where AE is zero for a fundamental loop. Three algorithms are presently in significant use and gradient method is employed to handle.the nonlinear terms in Equation (4.6). For a single pipe section, Equation (4.6) can be written as f(Q) = KpQn + K2 P(Q) (4.7) which represents the grade difference across a pipe section carry ing flow Q. Substituting an estimate, Qi, f._r: Q, and denoting f(Qi) by Hi, Equation (4.7) becomes Hi a f(Qi) = KpQ + KM Q ' P(Qi) (4.8) Differentiating Equation (4.7) and setting Q = Qi, gives the gradient of the function at Q p Q1. Thus, n l f(Qi) = nKpQi I 2KM i P'(Q.) Denoting f'(Qi) by Gi, thus n1 G nKpQi + 2K i P'(Qi) (4.9) Both the function and its gradient, evaluated at Q = Qi, will be used in all three algorithms for solving loop equations. 4.1.1.1 Single Path Adjustment (P) Method This method was first described by Hardy Cross as the "Balanc ing Head Method" which was limited to closed loop systems and included only line losses. The procedure is generalized and summar ized as follows: (i) An initial set of flowrates which satisfies continuity at each junction node is determined, (ii) A flow adjustment factor is computed for each path (.+fl) to satisfy the energy equation for that path and continuity must be maintained when applying the correction factor. (iii) Step (ii) is repeated using improved solutions until the average correction factor is within a specified limit. Equation (4.6) is used to compute the adjustment factor for a path using gradient method to linearize the nonlinear energy equations. Thus, f(Q) = f(Qi) + f'(Qi) AQ (4.10) in which AQ =Q Q i, where Qi is the estimated discharge. Applying Equation (4.10) to Equation (4.6) and solving for AQ gives AE Hi AgQ ____ (4.11) EGi which is the flow adjustment factor to be applied to each pipe in the path. The numerator represents the imbalance in the energy relationship due to incorrect flowrates and this procedure reduces this to a negligible quantity. Flow adjustment is carried out for all i fundamental (closed) loops and (fl) pseudo loops in the network. 4.1.1.2 Simultaneous Path Adjustment (SP) Method This algorithm is similar to the corrective mesh flow method described in Section 2,8, the only difference is that gradient method is used here instead to linearize the energy equations. It is developed to improve convergence by simultaneously adjusting the flowrate in each loop representing an energy equation. The method is summarized as follows: (i) An initial set of flowrates which satisfy continuity at each junction node is determined. (ii) A flow adjustment factor is simultaneously computed for each loop to satisfy the energy equations without disturbing the continuity balance. (iii) Step (ii) is repeated using improved solutions until the flow adjustment factor is within a specified limit. The simultaneous solution of t + f 1 equations is required to determine the loop flow adjustment factors. Each equation in cludes the contribution for a particular loop as well ascontri butions from all other loops which have pipes common to both loops. For loop j, the head change required to balance the energy equation is expressed in terms of the flow change in loop j ( AQj) and the flow changes in adjacent loops (A Qk) as follows: af na, f(Q) = f(Qi) + aQ A Qj + aQ or, f(Q) = f(Qi) f(Q) AQj f'(Qi) A Q (4.12) Substituting f(Q) = E, fi) KHi f'(Qi) ZG, Equation (4.12) becomes AE li ( ci) A Qj + Z Z(GiAQk) (4.13) in which Z.Hi sum of the head changes for all pipes in loop j ( Gi) AQA j sum of all gradients for the same pipes times flow change for loop j. 4Z(G 1 Qk) u sum of gradients for pipes common to loops j and k multiply by the flow change for loop k. A set of simultaneous linear equations is. formed in terms of flow adjustment factors for each loop representing an energy equa tion. The solution of these linear equations provides an improved solution for another trial until a specified convergence criterion is met. 4.,1..3 Wood's Linear _(L Method This method developed by Wood (Wood, 1981) involves the solu tion of all the basic hydraulic equations for the pipe network. However, only the energy equations need to be linearized as the continuity equations are all linear. Using gradient approximation, the energy equations are linearized in terms of an approximate flowrate, Qi as follows: f(Q) f(Qi) 4 f'(Qi)(QQi) Introducing Hi and G. as before, the above equation becomes (E Gi)Q= ,(GiQji Hi) + AE (4.14) This relationship is employed to formulate (j + f 1) energy equations which together with the j continuity equations, form a set of p simultaneous linear equations in terms of the flowrate in each pipe. One significant advantage of this scheme is that an arbitrary set of initial flowrates, which need not satisfy continuity, can start the iteration. A flowrate based on a mean flow velocity of 4 ft/sec has been used by Wood (Wood, 1981). The solution is then used to linearized the equations and successive trials are carried out until the change in flowrates between successive trials become insignificant. 4.1.2 ALGORITHMS FOR SOLVING NODE EQUATIONS Two methods for solving the node equations are also widely used and are described here for completeness. 4.1.2.1 Single Node Adjustment (N) Method This method was also first described in the paper by Hardy Cross and is known as the "Balancing Flows Method". The procedure is outlined as follows: (i) A reasonable grade is assumed for each junction node in the system. The better the initial assumptions, the fewer the required trials. (ii) A grade adjustment factor for each junction node which tends to satisfy continuity is computed. (iii) Step (ii) is repeated using improved solutions until a specified convergence criterion is met. The grade adjustment factor is the change in grade at a particular node (AH) which will result in satisfying continuity and considering the grade at adjacent nodes as fixed. For conven ience, the required grade correction is expressed in terms of Qi which is the flowrate based on the values of the grades at adjacent nodes before adjustment. Thus, using gradient approximation, f(Q) = T(Qi) + f'(Qi).AQ with the usual substitution, AQ (l/Gi) AH (4.15) where AH = H Hi, the grade adjustment factor and ZAQ denotes the flow corrections required to satisfy continuity at nodes. From Equation (4.3), AQ Qi QG (4.16) Thus, from Equations (4.15) and (4.16), E Qi Qe a H (4.17) E (1/Gi) In Equation (4.17), inflow is assumed positive. The numerator re presents the unbalanced flowrate at the junction node. Q,, the flowrate in a pipe section prior to adjustment, is computed from iK/n Qi = ( Hi/K)* in which AHi grade change based on initial assumed values of grade. If pumps are included, the following expression is used to determine Qs SAHi KQi P(Qi) (4.18) Equation (4.18) is solved using an approximation procedure. Adjust ment of the grade for each junction node is made after each trial until a specified convergence criterion is satisfied, 4.1.2.2 Simultaneous Node Adjustment (SN) Method This method requires the linearization of the basic pipe net work node equations in terms of approximate values of the grade. If the discharge in Equation (4.3) is expressed in terms of the assumed heads, it can be written as: H H (4.19) KIab for any node, a, and b denotes an adjacent node. Equation (41.19) can be linearized with respect to grades if the flowrates are Vritten in terms of some initial values of the grades, Hai and Hbi, and the corrections in these grades. The gradient method is again used to calculate the flowrate in pipe section, ab. Thus, Q H Q + IQ A Ha + Q AnHb (4.20) Ha fHb in which Q ( (Ha Hb)/Kab)/n (4.21) A Ha = Ha Ha ( adjustment factor for head at node a) AHb = Hb Hbi ( adjustment factor for head at node b) Substituting the partial derivatives of the flowrate expression in Equation (4.21) in Equation (4.20) and simplifying, gives 1n Q = Qi(l 1/n) + (H Hb) (4.22) nKab The initial value of the flowrate, Qi is computed based on the initial values of the grades. Thus, Qi= ((Hai Hbi)/Kab)1 where Kab may include minor losses, if any. Using Equation (4,22), the continuity equation for each junction node can be expressed as a linear function of the variable and fixed grades of adjacent nodes and the variable grade of junction, a. Hence, NHence 1n N Q1n _Qi H Ha i b=1 nKab b; nKab a (4.23) N NV H NF 1ln Qe + Qi + ( ) Hb b=1 b1 n b=1 nKab where N refers to all adjacent nodes, Nv refers to adjacent variable grade nodes and NF refers to all adjacent fixed grade nodes. Qi is positive for outflow. Equation (4.23) is written for each junction node in the system resulting in a set of linear equations in terms of junction node grades. If pumps are included, two additional nodes may be assigned to a pump at the suction and discharge sides as shown schematically below: a c d Two additional equations can be written: ab Ha Hb aKcd (H d) (,2 ) He Hb P[((HHd)/Kcd)1/n ] (4.25) Equation (4.24) is just the continuity equation and Equation (4.25) relates the head change across the pump to the flow in either the discharge or suction line. Equation (4.25) can be linearized using gradient method as follows: Let Y = P(Qi) + Hb Hc 0 (4.25a) Using the gradient approximation, aY 9Y 3Y Y = Yi + i Hb C H4 AH + AH+ad H (4.26) Substituting the partial derivatives in Equation (4.26) and simpli fying, we have the following linearized equation: Hc(1 I+ ) Hb Hdp (4.27) where ac and p depend on the relationship used to describe the pump, P(Q), and are given by o = p(Q.) 1 P'(Qi) p = P (Qi)/(nKcd  A set of (j + 2N ) simultaneous linear equations (where Np number of pumps) is generated and solved starting with Qi's based on any assumed set of junction node grades. An improved set of junction node grades is then used to compute an improved set of Q. s and the procedure repeated until a specified convergence criterion is satisfied. 4.2 CC'. i.!l'i.'; ON ALGORITHMS USING GRADIENT METHOD Node equations are easier to formulate because the equations include only contributions from adjacent nodes. On the other hand, the loop equations require the identification of an appropriate set of energy equations which include terms for all pipes in funda mental loops and between fixed grade nodes. Computer formulation of this set of equations is considerably more difficult than form ulation of the node equations, Each of the procedures described is iterative in nature and computations terminate when a specified convergence criterion is met. The solutions are therefore only approximate although they can be very accurate, The ability of an algorithm to produce an acceptable solution is of prime concern and studies have demons trated that convergence problems exist and an accurate solution is not always possible. 4.2.1 ACCURACY OF SOLUTIONS A solution is considered accurate only when all the basic equations are satisfied to a high degree of accuracy. For the three methods based on loop equations, the continuity equations are exactly satisfied. Each of these methods then proceeds to satisfy the energy equations iteratively and the unbalanced heads for the energy equations is evidence of solution accuracy. For methods based on node equations, iterations are carried out to satisfy con tinuity at junction nodes and the unbalance in continuity is a signi'icant indication of solution accuracy. 4.2.2. RELIABILITY OF ALGORITHMS A study carried out by Wood, using an extensive data base, has shown that the P, N and SN methods exhibited significant convergence problems (Wood, 1981). Since these methods are widely used, great care must be exercised when using them. SN method failures are characterized by the inability to meet a reasonable convergence criterion and if this occurs in a limited number of trials, further trials are usually of no benefit. Failure rate was quite high and the use of results obtained employing this method is not recommended unless a good accuracy is obtained in a reasonable number of trials. It has been established that algorithms based on node equations (N and SN methods) failed to provide reliable results because of the inability of these methods to handle low resistance lines. This is attributable to the fact that solution algorithms for these equations do not incorporate an exact continuity balance. For each of the three methods singled out above, failure rates can be reduced if initial values closer to the correct values can be determined. However, this is no easy task and as evidenced in the study, even an excellent set of initial conditions does not guarantee convergence. Both the SP and L methods provide excellent convergence and the attainment of a reasonable convergence criterion is sufficient to assure great accuracy. Convergence failure is very rare. However, since a gradient method is used to handle nonlinear terms,, there is always the possibility of convergence problems. Illconditioned data such as poor pump descriptions are particularly prone. The L method has some advantages over the SP method. Assumed arbitrary flowrates need not satisfy continuity as the continuity conditions are already incorporated into the basic set of equations. This method also allows a more straightforward and reliable inclu sion of hydraulic components such as check valves, closed lines, and pressure regulating valves. Although the SP method has signi ficantly less equations to solve, the use of sparse matrix techni ques to handle the larger matrix generated by the L method has somewhat negated this advantage. 4.3 LINEAR THEORY METHOD BY WOOD AND CHARLES In this section, the linear theory method (Wood and Charles, 1972) will be described and used in solving the system of equations formulated by loop analysis which considers flowrates as unknowns (hereafter referred to as the Qequations). Like the other linear method described in Section 4.1.1.3, it has several distinct advantages over the NewtonRaphson or Hardy Cross methods. Firstly, it does not require an initialization, and secondly, according to Wood and Charles, it always converges in a relatively few itera tions. However, its use in solving the head oriented equations or the corrective loop oriented equations is not recommended. Linear theory transforms the nonlinear loop equations into linear equations by approximating the head in each pipe by n1 h = (KQ ) Q = K'Q (4.28) n1 in which Qi is an estimate of the flowrate, and KI = KQi . Combining these linearized loop equations with the j. junction continuity equations provides a system of p linear equations which can be solved by Gaussian elimination in conjunction with sparse matrix techniques (Tewarson, 1973). In applying the linear theory method it is not necessary to supply an initial estimate, as maybe implied. Instead, for the first iteration each K' is set'equal to K, which is equivalent to setting all flowrates Qi equal to unity. In developing the linear theory method, Wood observed that successive iterative solutions tend to oscillate about the final solution. Reasons for the oscil lation can be understood by observing that the linear theory method is a variation of the NewtonRaphson method described in Chapter 3 whereby K' in Equation (4.28) is simply the derivative of hL if multiplied by n. The oscillation could be prevented by multiplying each K by its n, which involves more computation than averaging consecutive solutions as proposed by Wood. Thus, the flowrate used in a trial is just the average flowrate for that pipe from the previous two solutions, or Qi(m) Qi(ml) + Qi(m2) / 2 in which m within parentheses denotes a trial number. 4.3.1 INCLUSION OF PUMPS AND RESERVOIRS When pumps (not booster pumps) and reservoirs are connected to a network, the flows in the two connected pipes become addi tional unknowns and therefore an additional equation is required beyond the j continuity equations and fundamental loop equations. The additional equation is obtained from a pseudo loop, which connects the two reservoirs (fixed grade nodes) by a "no flow" pipe. If f fixed grade nodes exist in a network, there would be fl independent equations. Energy conservation around a pseudo loop (of which a fundamental loop is a special case) is defined by Equation (4.6). Thus, AE = E(KpQn + KMQ2) P( SP(Q) If the expression for P(Q) in Equation (2.13) is adopted in Equation (4,6), the linear theory method does not give rapid conver gence as it does when pumps and/or reservoirs are not present. A modification will therefore have to be made to allow the linear theory method to converge rapidly. The reason for the modification is that the head produced by a typical centrifugal pump decreases nearly proportional to the reciprocal of the square root of the flowrate whereas the head loss in a typical pipe increases approxi mately proportional to the square of the flowrate. A consequence of using this typical pump relationship in Equation (4,6) is that if the equation is solved by the linear theory method, convergence may become very slow if at all. This situation can be improved by a transformation of variables so that the new unknown has an exponent close to n. Such a trans formation is G r Q + B/2A (4.29) in which G is the new variable and A and B are the same constants in Equation (2.13). The appropriateness of Equation (4.29) is demon strated by solving it for Q and substituting in Equation (2.13). After some simplification, h = AG2 + h (4.30) p o where h = H B2/4A o o Obviously, the exponent of G (that is, 2) is close to the typical n. Substituting Equation (4.30) in Equation (4.6) gives (K Qn + KQ2) CAG2 = AE h+ ho (4,31) Addition of Equations (4,29) and (4.31) produces a system with as many equations as unknowns. 4.3.2 INCLUSION OF PRESSURE REGULATING VALVES (PRV'S) Networks containing PRV's may be analyzed by the linear theory method by initially assuming that the pressure (or head) immediately downstream from a PRV is constant and equal to the valve setting. Junction continuity equations are then written as if no PRV's are present. To write the loop equations, pipes containing PRV's are disconnected from the upstream nodes and the PRV's are replaced by dummy reservoirs. After each iteration a check on the flowrate Q, in each pipe containing a PRV is made. If there is any negative Q, the pseudo loop equation which includes terms for that pipe is modified with Q replaced by an unknown grade (head) immediately downstream from that PRV. M* J !" A '4 CHAPTER 5 ALTERNATIVE MATHEMATICAL APPROACHES & COMPUTATIONAL EXPERIENCE 5.1 MA THEMATIC L PROGRAMMING TECHNIQUES As a prelude to introducing the alternative approaches to solving the governing equations for a pipe network using optimi zation techniques, it is convenient to define a network topology using notations which are consistent with those used in graph theory. Let the network topology be described by a node set N and an arc set (network element) Eo0 In each of the set Eo, let Qij denote the flowrate from node i to. j. Each node, n in the set N is associated with a hydraulic head, H n Let R, a subset of N, be the set of nodes corresponding to reservoirs (fixed grade nodes) and let HA for all nkR be the fixed head associated with a res crvoir. Also let r. for all nEN denote the flow requirements (that is, supply or demand) at node n, For an incompressible fluid, the governing network equations can be stated as: Z Qnj in rn all n&N (5.1) (n,j) E (i,n) &E Srn 0 (5.2) n &N Hi Hj F (Q ), all (ij)& E (5.3) Hn H4 all n&R (5.4) Equation (5.1) is just a statement of mass conservation at each node while Equation (5.2) stipulates mass conservation for the network as a whole. Equation (5.3) states that the head loss Hi Hj =AHij aL u( across an element is some function Pij of the discharge through the element while Equation (5.4) requires that at a reservoir node, the head is constant. The functional form of Fij or its inverse Eij (AHi. ) Qij is not specified and can represent any element in cluding simple pipes and minor loss devices as long as a unique relationship between head anddischarge exists. In general,, Fij for most or all (ij) in Eo is nonlinear, thus necessitating iterative techniques such as (i) Hardy Cross, (ii) NewtonRaphson, and (iii) linearization, to be used to solve the governing network equations. Most of these techniques are detailed in Chapters 3 and 4. Each of these methods is simply a technique for solving a set of nonlinear simultaneous equations which have been adapted to the network analysis problem. Each is iterative in nature and begins with an initial trial solution. A new solution is obtained by solving a set of linear equations using straight forward procedures, If the new solution differs from the trial solution by less than a specified amount then the iteration stops. Otherwise, the new solution becomes the trial solution and the procedure is repeated. In some of the algorithms, an initial trial solution sufficiently close to the true solution is required to en sure convergence. The differences in the methods result from the use of different strategies to determine the new solution. The new approach by Collins, Cooper, Helgason and Kennington (1978) represents a radical departure from the state of the art iterative methods as optimization models are employed to solve the network problem. Two alternatives models are formulated and these models play analogous roles to the node versus loop formulations for solution of the network equations used in the state of the art methods. The first of the two optimization models, called the Content Model assumes the form 5i7 G P 'iJ r nsn 1 Minimize G F (t)dt H*dt + (ioj)&E o (gn)& E1 n y7, rp ^Hg d (n,g)& E1 o Subject to , S c Qnj inQ. r all n &NU(g) (n,j)> EUE1 (i,n)&EUE n Qij 2 0, all (i,j) 8E U(E1) in which E is the arc set for a.network in which the arcs have been replaced by two equivalent oppositely directed oneway elements so that Qij can assume only positive values. This replacement is done as a mathematical convenience so that Qij can be treated as a con strained decision variable in the optimization model which has a nonnegativity condition imposed on Qij. It can be proven that the solution obtained by solving the E network will produce identical results as those which would be obtained by solving the original Eo network which permits Qij to be unconstrained. The arc set E1 is merely a set of arcs connecting all nodes in N to a ground node g and is introduced to satisfy mass conservation for the network as a whole (Equation (5.2)). Using the terminology of Cherry and Millar (1951), the above problem is to find a set of flows which satisfies flow conservation and minimizes system content, G, .hence the name Content Model. The second optimization model, the CoContent Model, is a complementary (but not dual) model which has the form Minimize J i f ng (i,j)&E E.ij(t)dt nSN r n dt o 1J o subject to S Hi + AH. AHig 0, all (i,j)SE AHng H i H all nE R ng n 9g In the terminology of Cherry and Millar (1951), the above problem is to find a set of head losses whichsums to zero around all loops and minimizes system CoContent, J, hence the name Co Content Model. Using KuhnTucker theory (Kuhn and Tucker, 1950), it can be proved that the solution to either of these models yields the solution to the pipe network problem, that is, the optimal solution satisfies the governing network equations. The proof is carried out by examining the derivatives of the objective function and show ing that the derivative conditions for a stationary point, along with other constraints, are identical to the network equations. In the proof, it is assumed that the Fij and Eij functions are monotonically increasing. This assumption insures the convexity of the objective function which in turn guarantees the existence of a unique solution to the optimization problem. The monotonicity of Pij and Eij merely implies the fact that energy losses in a net work element increase with increasing discharge. The Content Model has the special structure of a convex cost network flow problem for which efficient routines are available. i.i, nrous nonlinear algorithms such as (i) FrankWolfe method, (ii) piecewise linear approximation and (iii) convexsimplex method are available for solving such a problem. The use of mathematical programming techniques in pipe net work analysis has paved the way for potential research in the following areas: (i) Extension of mathematical programming techniques to solution of compressible flow pipe network analysis problem. (ii) Incorporation of time variable storage in network elements to solve transient network problems. (iii) Use of mathematical programming techniques, to solve complex open channel networks. (iv) Feasibility of using mathematical programming techniques to solve network parameter identification problems such as the headdischarge relationship in pipe network analysis. (v) Development of an economic model to minimize the operation al costs for a flow network with operational behavior given by one or more network problems. 5.2 COMPUTATIONAL EXPERIENCE A computer program was written based on the linear theory method (Wood and Charles, 1972) described in Section 4.3. The pro gram was designed to solve the system of loop and node equations using the iterative procedure described by the method. Two features this FORTRAN computer program may have for general application include (i) the capability of handling networks containing pumps and reservoirs, and (ii) an algorithm which analyzes networks containing pressure regulating valves. The use of the node incidence matrix and fundamental loop matrix described in Section 2.1 in the algorithm has provided an efficient means of translating information contained in any pipe network into a network simulator. Incidentally, the node equations and loop equations were formulated using the node incidence matrix and fundamental loop matrix respectively. In carrying out all computations, friction losses in pipes were assumed to be described by the HazenWilliams equation and pumps were described by the quadratic form (Equation 2.13). The convergence criterion employed was:  o0.0005 in which Qi is the flowrate obtained for a trial and Qi, is the flowrate obtained from the preceding trial. This appears to be a stringent requirement which may assure good accuracy if the condi tion is satisfied. However accuracy is achieved if and only if continuity at every node and the energy equations are exactly satisfied. A small scale network, taken from Jeppson (1977, p.109) and shown in Fig 5.1, was tested. This 8 pipe, 5 node network, with the properties given in Table 5.1, has 2 reservoirs, a pump and a pressure regulating valve. The solution to the test problem and the solution reported by Jeppson (1977, p.110), using the same theory are tabulated in Table 5.2. The results appear to be in good agreement. 180 ft 2.0 cfs (7) C1] (1) 2 e r Ofi 5450 ft PRV 1 200 ft ( )(3) (6) 2.0 cfs 48) (3) L] 1.0 cfs Fig 5.1 Test Problem Table 5.1 Network Parameters Pump Characteristics Discharge, cubic ft per sec Pipe Writer's Solution Jeppson's Solution 1 2.53 2.56 2 0.38 0.32 32.47 2,44 4 0.72 0.73 5 0.92 0.88 6 1.08 1.12 7 1.81 1.83 8 3.19 3.17 Table 5.2 Solution to Test Problem The detailed solution and program listing are contained in the Appendix. With this program, 'the test problem took 0.12 second of execution time oneanAmdahl 470 computer. The number of iterations required to meet the convergence criterion was 6. The subroutine used for solving the linearized set of loop equations and the linear node equations simultaneously was developed based on the Gaussian method of elimination improved by pivotal condensation (Tewarson, 1973). The capability of the program to handle a larger network has not been proven but it would have stretched the available storage of a computer to its limits if it has been tested. Storage space is primarily taken up by the final augmented matrix which comprises essentially the node and loop equations. The use of sparse matrix techniques instead of full matrix methods may extend the capability of the program to analyze larger networks of a few hundred pipes and nodes. *i* *}( ;* CHAPTER 6 CONCLUSION The Hardy Cross method which sparked off'the evolution of the numerous techniques of simulating pipe networks, is suitable only for relatively small networks. With the advent of the computer, and as larger and more complex networks were analyzed, the Hardy Cross method was found to frequently converge too slowly if at all. The classic method which is described in most hydraulics or fluid mechanics text books, is an adaptation of the NewtonRaphson method which solves one equation at a time before proceeding to the next equation during each iteration instead of solving all equations simultaneously. The single path and single node methods described in Sections 4.1.1.1 and 4.1.2.1 respectively, are basically the classic Hardy Cross methods. Procedures developed to improve the convergence of the single path method were described by Martin and Peters (1963) and later by Epp and Fowler (1970). The procedure involves the simultaneous computations of flow adjustments and was presented in Section 4.1.1.2. A similar approach has been developed for the node equations where all node equations are linearized and solved simultaneously. This method is described by Shamir and Howard (1968). All of the four methods mentioned so far require an initial guess as to the solution and the rate of convergence depends to a degree on how close this initialization is to the correct solution. For the system of equations which is flowrate oriented, two linearization techniques (Wood, 1981 and Wood and Charles, 1972) were described in Sections 4.1.1.3 and 4.3 respectively. Both of these procedures do not require an initializationn and have been reported to converge in a relatively few iterations. Significant convergence problems were reported for the Single Path, Single Node and Simultaneous Node Methods (Wood, 1981). It has been suggested that if a specified stringent convergence crit erion cannot be met using single adjustment methods, the solution is probably unreliable. For the simultaneous node adjustment method, it has been suggested that the best indication of an acceptable solution is that the average relative unbalanced flow at the junc tion nodes be less than 2%. Instances of failures have also been reported in cases where line losses vary greatly or pumps operate on steep curves even when good initial approximations are available. The simultaneous path methods and the linear method using gradient approximations, were reported to provide excellent conver gence and the attainment of a stringent convergence criterion is sufficient to assure great accuracy in most cases. In the study carried out by Wood (1981), in which a wide variety of situations was represented, some incorporating features which increase conver gence difficulties like low resistance lines; these methods were reported to attain accurate solutions in a relatively few iterations. However, if gradient approximations are used to handle nonlinearity, convergence problems are always a possibility, especially if ill conditioned data such as poor pump descriptions are employed. Of all methods, the linear methods developed by Wood and Charles (1972) and a later version by Wood (1981), who used gradient approximations, offer more advantages. A balanced initial set of flowrates is not required since the. continuity conditions are already incorporated into the basic set of equations. These algor ithms permit a more direct and reliable incorporation of hydraulic components such as check valves, closed lines and pressure regulat ing valves. For any pipe network simulators to be of general use, these components, which affect continuij.ty, and their effects on the hydraulics of the network must be incorporated into the basic set of equations. However, the set of equations solved by the linear methods involved significantly more equations which will be a setback if full matrix methods are used. The use of sparse matrix techniques has somewhat corrected this. disadvantage and has rendered it a more desirable algorithm to adopt for analysis of pipe networks. The use of mathematical programming techniques in pipe net work analysis holds a lot of promise for the future. One of the direct consequences of the theory described in Section 5.1 is the identification of a unitary measure by which the goodness of a solution can be gaged. Traditional methods described previously give no good insight into the goodness of an approximate solution, particularly for large scale problems. The optimization models remove the vagueness that inherently surrounds a definition of " close when an attempt is made to utilize a comparison of indi vidual flows, heads, or losses in individual elements. Optimization techniques also have their setbacks. One is that functions describing friction losses, minor losses in pipes and pump heads must necess arily be convex functions for a solution to be guaranteed. In addition, head loss must be a unique function of discharge. Such uniqueness may not exist for certain control elements such as check valves and pressure regulating valves. Until these problems are re solved, its application will be limited in scope. *i i i * REFERENCES (1) Adams RW., Distribution analysis by electronic computer, Institution of Water Engineers, 15, 415428, 1961. (2) Bazaraa M.S. and J.J. Jarvis, Linear programming and network flows, J. Wiley & Sons, New York, 1977. (3) Belevitch V., Classical network theory, HoldenDay, Inc., San Francisco, 122, 1968. (4) Bellemy C.J., The analysis of networks of pipes and pumps, Journal Institution of Engineers, Australia, 37(45), 111116, 1965. (5) Bree D.W.,Jr., AI. Cohen, L.C. Markel, & H.S. Rao, Studies on operations planning and control of water distribution systems, Final Technical Report to OWRR Project No. 5214, Systems Control, Inc., Palo Alto, Calif., 1975. (6) Brock D., Closed loop automatic control of water system operations, Journal American Water Works Association, 55(4), 467480, 1963. (7) Brock D.A., Metropolitan water system operation subsequent to nuclear attack or natural disaster, Report No. AD 711.956, Dallas Water Utilities, City of Dallas, Texas, 364 pp, 1970. (8) Chenoweth H., and C. Crawford, Pipe network analysis, Journal American Water Works Association, 66(1), 5558, 1974. (9) Cherry E.C. and W. Millar, Some general theorems for nonlinear systems possessing resistance, Phil. Mag., (Ser 7), 42(333) 11501177, 1951. (10) Clay R,, Nonlinear networks and systems, Wiley Interscience, New York, 139, 1971. (11) Collins A.G., and R.L. Johnson, Finiteelement method for water distribution networks, Journal American Water Works Association, 67(7), 385389, 1975. (12) Collins N.A., L. Cooper, V. Helgason, and J.L. Kennington, Solution of largescale pipe networks by improved math ematical approaches, IEOR 77016WR 77001, 1978. (13) Dillingham J.H., Computer analysis of water distribution systems, part 1, Water and Sewage Works, 114(1), 13; part 2, 114(2), 4345; part 4, Program application, 114(4), 141142, 1967. (14) Donachie RP., Digital program for water network analysis, Journal Hydraulics Division, Proc. Amer. Soc. Civil Engineers, 100 (HY3), 393403, 1973. (15) Eggener C.L., and L.B. Polkowski, Network models and the impact of modeling assumptions, Journal American Water Works Asso ciation. 68(4), 189196, 1976. (16) Epp R., and A.G. Fowler, Efficient code for steadystate flows in networks, J. Hydraulics Div., Proc. Amer. Soc. Civil Engineers, 96 (HYI), 4356, 1970. (17) Fietz T.R., Discussion of "Hydraulic network analysis using linear theory", J. Hydraulics Div., Proc. Amer. Soc. Civil Engineers, 99 (HT5), 855857, 1973. (18) Fietz T.R., Discussion of "Efficient algorithm for distribution networks", J. Hydraulics Div., Proc. Amer. Soc. Civil Engineers, 99 (HY1l), 21362138, 1973. (19) Gerlt J.L., and G.F. Haddix, Distribution system operation analysis model, J. Amer. Water Works Association, 67(7) 381384, 1975. (20) Graves O.B., and D. Branscome, Digital computers for pipeline network analyses, J. Sanitary Engineering Div., Proc. Amer. Soc. Civil Engineers, 84 (SA2), 118, 1958. (21) Guillemin E.A., Tntroductory circuit theory, John Wiley & Sons, Inc., New York,.510, 1953. (22) Hoag L.N., and G. Weinberg, Pipeline network analyses by elec tronic digital computer, J. American Water Works Association, 49(5), 517524, 1957. (23) Hudson W.D., Computerizing pipeline design, J. Transportation Engineering Div., Proc. Amer. Soc. Civil Engineers, 99(TE1) 7382, 1973. (24) Hudson W.D., A modern metroplex looks ahead, J. Transportation Engineering Div., Proc. Amer. Soc. Civil Engineers, 100(TE4) 801814, 1974. (25) Jeppson R.W., Analysis of flow in pipe networks, Ann Arbor Science, 1977. (26) Karni S., Intermediate network analysis, Allyn and Bacon, Inc., Boston, Mass., 84113, 1971. (27) Kuhn H.W., and A.W. Tucker, Nonlinear programming, in Proceedings of the Second Berkeley Symposium on Mathematical Statistics and Probability, Berkeley, Calif., 481492, 1950. (28) Lam C.F., and M.L. Wolla, Computer analyses of water distribution systems, part 1 formulation of equations, J. Hydraulics Div., Proc. Amer. Soc. Civil Engineers, 98 (HY2), 335344, 1972. (29) Lam C.F., and M.L. Wolla, Computer analysesof water distribution systems, part 11 numerical solution, J. Hydraulics Div., Proc. Amer. Soc. Civil Engineers, 98 (HY3), 447460, 1972. (30) Lemieux P.F., Efficient algorithm for distribution networks, J.iHydraulics Div., Proc. Amer. Soc. Civil Engineers, 98 (HY11), 19111919, 1972. (31) Liu K.T., The numerical analyses of water supply networks by digital computer, Thirteenth Congress of the International Association for Hydraulic Research, 1, 3643, 1969. (32) Marlow T.A., R.L. Hardison, H. Jacobson and G.E. Beggs, Improved design of fluid networks with computers, J. Hyd draulics Div., Proc. Amer. Soc. Civil Engineers, 92 (HY4), 4361, 1966. (33) Martin D.W. and G. Peters, The application of Newton's method to network analyses by digital computer, Institution of Water Engineers, 17(2), 115129, 1963. (34) McIlroy M.S., Pipeline network flow analyses, J. Amer. Water Works Association, 41, 422428, 1949. (35) McPherson M.B., E.C. Bolls, Jr., D.A. Brock, E.B. Cobb, H.A. Cornell, J.E. Flack, F. Holden, F.P. Linaweaver, Jr., R.C. McWhinnie, J.C. Neill, and R.V. Alson, Priorities in distri bution research and applied development needs, J. Amer. Water Works Association, 66(9), 507509, 1974. (36) Rao H.S., D.W. Bree, Jr., and R. Benzvi, Extended period simu lation of water distribution networks, Final technical report OWRR project no. C4164, Systems Control Inc., Palo Alto, Calif., 120 pp., 1974. (37) Rao H.S., D.W, Bree, Jr., and R. Benzvi, Extended period simu lation of water systems part A, J. Hydraulics Div., Proc. Amer. Soc. Civil Engineers, 103 (HY3), 97108, 1977. (38) Rao H.S., L.C. Markel and D.W. Bree, Jr., Extended period simu lation of water systems part B, J. Hydraulics Div., Proc. Amer. Soc. Civil Engineers, 103 (HY3), 281294, 1977. (39) Shanir U., Water distribution systems analysis, IBM Research Report RC 4389 (No. 19671), IBM Thomas J. Watson Research Center, Yorktown Heights, New York, 1973. (40) Shamir U., and C.D. Howard, Water distribution system analysis, J. Hydraulics Div., Proc. Amer. Soc. Civil Engineers, 94 (HY1),219234, 1968. (41) Shamir W., Optimal design and operation of water distribution systems, Water Resources Res., 10(1), 2736, 1974. (42) Tewarson R.P., Sparse Matrices, Academic Press, Inc., 1973. (43) Williams G.N., Enhancement ofconvergence of pipe network solutions, J. Hydraulics Div., Proc. Amer. Soc. Civil Engineers, 99 (HY7), 10571067, 1973. (44) Wood D.J., and C. Charles, Hydraulic network analysis using linear theory, J. Hydraulics Div., Proc. Amer. Soc. Civil Engineers, 98 (HY7), 11571170, 1972, (45) Wood DJ., An explicit friction factor relationship, Civil Engineering, 36(12), 6061, 1966. (46) Wood D.J., Algorithms for pipe network analysis and theirreli ability, University of Kentucky, Water Resources Research Institute Research report No. 127, 1981. (47) Zarghamee M.S., Mathematical model for water distribution systems, J. Hydraulics Div., Proc. Amer. Soc. Civil Engineers, 97 (HY1), 114, 1971. APP' 'IX Vf 44l; %4 3 '.1' C: 4~" 4~ 2 O i2 44r 4i 4 4` PIPE t iQ^c 44 1)4 B jSP st^'' 'I^ A ^ ~44 4)4 44l~ ~ :1~ 4P 4 A (. > ^ X ( a 444 X <. 4i 444 4 44 4 4 PEA4 t C.f 4  I~"" ANTS 1* A; f P NO. rA, 9 & Mr i A .y.c A~ 21 FH St \) 4. Xil ;:i d iiscr^  *' ;, 4~3 '" kiI%3 i:lj C,~M 4' :1 i. 8; iXI .L1 C ~1 i"'B '13~ p 1 F C(~ i~9 I TVIAONr; t~ii ?.4 I"~ ~ A;fiP~: 444 ~ih,  Ai / .4 2i"P , 4 F li~,l FtMi:~ari '.." ii I : ../* i f4 121 4 27 18 ID 14 4 T ( I 0 = < 'i:: 4 `4 1' 'i^ i 241 245 iS"S a 5 .4 ~444 % 4 5 'i~ 64. 5 1* K24 a 1'' S J 19r ~ *: IK 1 P ) HP 'A$ cl I 'I '1 i ','' 1 1 iTB._Y PR iWT TD! 4 C A,~ SOF Q A W Aj X 2 Mf) MA(1 . S ]' A A 9 ."  9 543 r r B i"~"i.i TE3 j7 rP:i HOL ( P!~ rI 0 R A U L. I c 1* 1 wvwwsmvysis 11 1M It: K! n! ] Pot Lv' / 4 C KI rA D Ec S G Hi~~~ 9 9%WNN @ S. 4 : t / J0yNS T ft E E" OWNSTRE""I F;r P R V Ei: "~ P f R~i i;X."~~ :s rsv n io e "%P"L`I i; g ~t: i~ ~ C 
Full Text 
PAGE 1 Publication No. 77 PIPE NETWORK ANALYSIS By MunFong Lee University of Florida Gainesville PAGE 2 PIPE NETWORK ANALYSIS BY MUNFONG LEE NONTHESIS PROJECT REPORT IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE MASTER OF ENGINEERING DEGREE DEPARTMENT OF ENVIRONMENTAL ENGINEERING SCIENCES UNIVERSITY OF FLORIDA NOVEMBER 1983 PAGE 3 ABSTRACT Analyses of large scale pipe networks are needed whenever significant changes in patterns or magnitudes of demands or supplies occur in municipal water or gas distribution systems. Changes of this nature occur whenever new industrial and residential areas are being developed or new sources of supply are tapped. In the absence of such analytical tools to determine the ance of an existing system under new demands, needlessly large investments are made for larger than necessary pipes, redundant lines or duplicate facilities. Another cause for concern is the ability of the numerous algorithms to provide reliable'results without which deficient engineering judgments may be made in engineering applications deal ing' with large scale pipe networks. Convergence and reliability problems of most of the algorithms are highlighted after the theoretical background has been presented .As an aid to more effective formulation of the loop and nodal equations, the essential concepts of network theoryare also presented together with the fundamental hydraulic principles forming the backbone of the state of the art iterative procedures. This report concludes with a new approach which employs optimization techniques to solve the pi'pe network problemas a viable and perhaps more versatile alternative to the widely used iterative methods. i PAGE 4 ACKNOWLEDGEMENT I wish to express mi appreciation to my graduate committee chairman, Dr. James P. Heaney, for his advice and enthusiasm in directing this report and my general course of study at the University of Flotida. Special thanks are due to the members of my committee, Dr. Wayne C. Huber and Dr. Warren Jre9 for their assistance and review of this report. I wish also to record my grateful thanks to my employer, the Public utilities Board, Republic of Singapore, for providing me this opportunity to further my studies. ****rr i1 PAGE 5 TABLE OF CONTENTS ABS TRAC Til., ...... ................... It 0 ...... ......... ... e ............. PAGE i ii iii ACKNOWLEDGEMENT II .. .. .. ft .. .. It .. .. .. .. 8 .. .. .. e TABLE OF CONTENTS .. .. ".' II .. III II ......... ... to .. .. lJ .. III ....... tI .. III III III CHAPTER 1 1.1 1.2 1.3 CHAPTER 2 2.1 2 2 2.3 2. LI2.5 2.6 2.7 2.p8 CHAPTER 3 3.1 3.2 CHAPTER 4 4.1 4.1.1 4.1.1.1 INTRODUCTION ... ., .... .......... .. It .............. Problem Definition ............. ............ Signi'ficance .................... ..... ........ A II ........ .. Moti va.tiOtl ............... ............... ........ .......... NETWORK ANALYSIS & PIPE NETWORK HYDRAULICS Fundamentals of Network Theory Pipe Network Conservation Laws II .. .. .. .. .. .. .. .. .. .. .'.111 ............ Friction and Minor Losses ..................... ............ 1 2 2 3 5 5 7 9 Pumps ... ....................... ,. ......... '. (I' ..... '" .. .. .. 11 Pressure Regula ting Valves .. .. .. .. .. .. .. .. Node Analysis Loop Analysis .. .. .. .. .. .. $ .. e O CI.'.& ., (II Corrective Mesh Flow Analysis ..... ,. rt NEWTONRAPHSON METHOD t ...................... 13 14 15 16 18 Application to Node Equations 18 Application to Corrective Mesh Flow Equations 19 LINEAR METHODS ., It II t Gradient Method flO .................... Algori thms' for the Solution of Loop Equations .Single Path Adjustment (p) Method ............. \I iii 21 21 22 23 24 PAGE 6 4.1.1.3 4.1.2 4.1.2.1 4.1.2.2 4.2 4.2.1 4.2.2 4.3 4.3.1 4.3.2 CHAPTER 5 5.1 5.2 Wood's Linear (L) Method lit .. .. ....... "lit III ....... fJ ....... Algorithms for Solving Node Equations .. ... Single Node Adjustment (N) I1Iethod .. .. ... .. .. ... It .. .. Node Adjustment (SN) Method Comments on Algorithms using Gradient Method a. Accuracy of Solutions ............................ Reliabili ty of Algori thms .; PAGE 26 27 27 28 31 31 31 Linear 'Theory Method by Wood and Charles .. 33 Inclusion of Pumps and Reservoirs ... 34 Inclusion of Pressure Regulating Valves 36 ALTERNATIVE MATHEMATICAL APPROACHES & COMPUTATIONAL EXPERIENCE Mathematical Programciing Techniques ., iI ..... II ....... II!! .. III ., ".8 Computational Experience ............................................... 37 37 41 CHAPTER 6 : CONCLUSION ........................................ q5 REFERENCES ." ................ II ........ \I ............ II It .............. .......... ,It ... " 'ft .. ". 48 APPENDIX ........ ". ill ............................... 0 ......... ........................... '" ... .. .. 53 iv PAGE 7 CHAPTER 1 INTRODUCTION and design of pipe networks create a relatively complex problem, particularly if the network consists of a range of pipes as frequently occurs in water distribution systems of large metropolitan areas, or natural gas pipe networks. In the absenc.e of significant fluid acceleration, the behavior of a network can be determined by a sequence of steady state conditions, which form a small but vital component for assessing the adequacy of a network. Such an analysis is needed each time changing patterns of consumption or delivery are significant or addon features, such as supplying new subcli visions, tion of booster pumps, pressure regulating or storage tanks change the system. The steady state flows of a network are governed by the laws I of conservation of energy and mass and the classical pipe network analysis problem is to establish the steady state flows and pressures in a full flow closed conduit network of known physical characteristics. Due to the complexity and the nonlinearity of works, solving the network analysis problem is not a trivial exercise. For over four decades, a number o.f algorithms have been developed since the pioneering work of Hardy C:::oss. All of these techniques are iterative in nature, differing only in the method in which an estimate of the true solution is obtained. A recent study (Collins, Cooper, Helgason and Kennington, 1978) uncovered a new approachio the pipe network analysis problem using optimization techniques which represent a: radical departure from the traditional state of the art PAGE 8 methods. This report attempts to provide a comprehensive writeup of the theory behind some of the more commonly used algorithms and their efficiency and reliability. 1.1 PROBLEM DEFINITION A pipe network is physically a collection of. interconnected elements such as pipes, pumps, reservoirs, valves, and similar appurtenances. Mathematically, the network is represented as an edge set consisting of pipes, pumps, valves and similar elements and a node set comprising reservoirs and element intersections In most of the elements, a unique functional relationship between pressure and flow exists. Pressure, in incompressible flow networks, can be expressed in terms of an equivalent hydraulic head, a terminology which will be adopted throughout this report as is standard practice. The steady state condition of a network can be completely defined by the head at each node and the flow in each element. Having determined this unique set of flows/ heads for a given set of inputs and withdrawals, all other quantities of interest can be deduced therefrom. 1.2 SIGNIFICANCE Steady state network analysis is a basic tool in water distribution system management and design. It can also be used to .... developopera ting policies and strategies to not" only reduce operating costs but also increase reliability and reduce water wastage (Brock, 1970; Hudson, 1974; Rao et ale 1974, 1977; Shamir, 1974; Bree et al. 1975). Application of steady state network 2 PAGE 9 analysis in online system control is also receiving growing attention (Brock, 1963; Hudson, 1973: McPherson et al. 1974; Rao et al. 1974; Gerlt and Haddix, 1975; Eggener and Polkowski, 1976) 1.3 MOTIVATION Since Hardy Cross first provided a solution for the pipe network analysis problem, three general methods which are widely used today, have evolved: 3 (i) Hardy Cross (Hoag and Weinberg, 1957; Graves and Branscome, 1958; Adams, 1961; Brock, 1963; Bellamy, 1965; Dillingham, 1967; Fietz, 1973; Williams, 1973; Chenoweth and Crawford, 1974; Eggener and Polkowski, 1976) (ii) NewtonRaphson (Martin and Peters, 1963; Shamir and Howard, 1968; Liu, 1969; Epp and Fowler, 1970; Zarghamee, 1971; Lam and Wolla, 1972;' Lemieux, 1972; Donachie, 1973; Rao et al. 1974,1977) (iii) Linearization (McIlroy, 1949; Marlow et ale 1966; Wood, and Charles, 1972; Fietz, 1973; Collins and ,Johnson, 1975). These methods solve a set of nonlinear simultaneous equations iteratively beginning with an initial trial solution. The iteration is complete when a new solution differs from the trial solution by less than a specified amount; otherwise, the new solution becomes the trial solution and the procedure is repeated. Differences in the above methods arise because of the strategies used to determine a new solution. PAGE 10 In view of the iteratiye nature of these methods, large scale networks with hundreds of nodes and elements require considerable computer efforts to solve. The choice of algorithm therefore, depends on the computational speed and reliability of a particular solution procedure. Matrices associated with water distribution networks, like most manmade systems, are sparse. One of the keys to faster convergence and hence to greater computational efficiency and perhaps reliability for most, if not all, algorithms is the use of sparse matrix techniques in the solution procedures (Tewarson, 1973). In the following chapters, most of the essential tools required for the analysis of incompressible flow in pipe networks are presented. Chapter 2 introduces graph theory which is useful in the formulation of pipe network simulators and also includes fundamental hydraulic principles governing pipe networks, to provide the necessary groundwork for the development of the loop and node system of equations. In Chapters 3 and 4, methods for solving these systems of nonlinear equations are described. Alternative mathematical approaches and the writerts computational are presented in Chapter 5. 4 PAGE 11 5 CHAPTER 2 NETWORK ANALYSIS & PIPE NETWORK HYDRAULICS There has been growing awareness that certain concepts and tools of network theory are very useful in the analysis of pipe networks especially in the formulation of computer simulators. The theory of network analysis is well established and several references in this field are available (Gulliman, 1953; Belevitch, 1968: Karni, 1971; 1971; Shamir, 1973; Bazaraa,and Jarvis, 1977; Minieka, 1978). For consistency, the terminology used in this chapter has been adopted for pipe networks. Also presented in this chapter are some of the fundamental hydraulic principles which form the foundation of the three traditional described in Chapter 3. 2.1 FUND.t\MENTALS OF NETWORK THEORY According to network a network is a granh consisting of a set of junction points called nodes, with certain pairs of nodes being joined by line segments dges (or arcs, branches or links). Edges joining the same two nodes are multiple edges and a node withou"t an edge connected to it is an isolated node. If a node has only one edge connected to that edge is a pendant. An edge and a node at the end of the adge are said to be im..iJlent. A 8ubgraph is any collection of nodes and edges comprising only nodes and edges of a larger graph. The complement of a subgraph is the collection of nodes and edges remaining after the removal of the subgraph. A path between two nodes is a subgraph whose terminal nodes PAGE 12 each have only one arc incident and all other nodes are incident to exactly two arcs. A graph is said to be connected if there is a path connecting every pair Of nodes. A connected in which each node of the subgraph is incident to exactly two arcs of the subgraph is called a (or cycle). A tree is a connected graph containing no loops. The comp lement of a tree is a cotree. Edges of a cotree are links. A tree containing all nodes of a graph is a spanning tree. An edge of a graph is said to be directed (or oriented) if there is a sense of direction ascribed to the edge'. If all edges, 6 of a graph are it is called a directed graph. However, a network need not be directed because it may be feasible to have flow in either direction along an edge. The flow capacity of an edge in a specified direction is the upper limit to the feasible magnitude of the rate of flow in the edge in that direction. The flow capacity may be any nonnegative quantity, including infinity. An edge is directed if the flow capacity is zero inane direction. The topology of a directed graph nodes and e edges can be described by a X e node incidence matrix, A with typical element ( + 1, if edge j is directed away from node i ) ( ) a .. :: ( 1, if edge j is directed towards 'node i ) lJ ( ) ( 0, if edge j is not connected to node i ) For a connected graph it is apparent each column of A will contain a 1 and a 1 and all remaining elements will be 2ero. As a ,check, addition of the rows of A should yield a zero row. Thus, rank, of A is at most? 1. PAGE 13 If loops are formed, one by one, by adding links, one at a time, to a given spanning .tree, it is apparent that each time a link is added a unique loop will be created. Such a loop is called a fundamental loop. A fundamental loop set for any connected graph, containing i\ ioops, can be described by a ?l Xc fundamental looI? matrix, B, with typical element ( +1, if edge j is in loop i and the direction of edge j ) ( ) ( is clockwise, say ) ( ) b .. ( 1, if edge j is in loop i and the direction edge j ) lJ ( ) ( is counterclockwise ) ( ) ( 0, if edge j is not in loop i ) By performing elementary row operations on B, an identity submatrix of order, A t can be obtained, implying the rank of B is I\. Both the node incidence matrix the fundamental loop matrix can be used to formulate the continuity and energy (or loop) sets of equations in a computer simulator. 2.2 PIPE NETWORK CONSERVATION LAWS Pipe network parameters are introduced to develop two conservation laws utilizing theory. The following notation shall be adopted for convenience. A directed network will be described by a node N and an edge set, E of ordered pairs of nodes. Each node n e. N is associated with a unique number called the head, Hn. For an edge directed from node i,to node j, an edge head loss is defined as ,,AHij 12 Hi Hj in which AHij :: LlHji' In each edge, a flow Qij exists, positive when the edge is directed from node i to node j. 7 PAGE 14 8 A basic law to be satisfied by the flows in a network is mass conservation, "'" Q.... Q rn nJ In (n,j)eE (i,j)6E all nN (2.1) where rn is the requirement node n, positive for inflows (supply) and negative for outflows (demand). Denoting the vector of Qij'S by q, equation (2.1) can be rewritten as A q iii .i. where r = (rl r 2 r n). As noted in section 2.1, A has a rank of ry 1, implying one of the in A is redundant and can be arbitrarily omitted. The matrix Ar obtained by deleting one row of A, say is defined as the node incidence matrix. A corresponding element in r is also deleted and a demand vector, d, defined as a = (rl' r2' ,ry)l) is introduced. Then ... Ar q :: d (2.3) It should be noted, in passing, that all the rows in A will be iridependent, that is, rank of A :: 7) if pumps and reservoirs are present in the network. However, a redundant row still exists if a junction is asswned at the reservoir or pumps. If the nodal head is unique, as assumed, then the summation of head losses around a loop is zero. This obvio.usproposition is used as the basis for the second fundamental network law Thus, if Lk is the edge set for edges in fundamental loop k, k ::: 1, 2, /\., then L ... 0 (i,j)(Lk lJ .... all k PAGE 15 Equation (2.4) can be written as B.6.h ::; 0 is defined as the vector of 6Hij s. If a mesh flow vector 6 (PI' P2, P3''1h ) the following relationship can be written BTp q ::: is defined, (2.6) Thus if to each fUndamental loop a unique' mesh flow is associated, the flow in any edge is a linear combination of the mesh flows for fundamental loops containing the edge in question. If pumps and reservoirs are included in the pipe network, equation (2.5) is generalized as follows: A E = BAh hp (2 7 ) where h :: pump head vector, l:::.. E :;: vector representing the, diffp erence in total grade between two reservoirs. In this generalized case, a junction node is assumed at a reservoir or pump and a pseudo loop is assumed to connect 2 reservoirs. 2.3 FRIC11ION & MINOR LOSSES The relation between head and discharge, that is, 6 hand q completes the number of equation sets required to define the network problem. Total head loss in a pipe, H, is the sum of the line loss, hLP' and minot loss, hLMo The line loss expressed, in terms of the discharge is given by: (2.8) where Kp is a pipe constant which is a function of line length,' diameter and roughness and n is an exponent.. Commonly used head 9 PAGE 16 loss equations include the DarcyWeisbach, HazenWilliams and Manning equations. Perhaps the most widely used of these equations is the equation, which is, C ARO.6J SO.54 ) English System (ES) Q := 1,.318 ) HW C ARO.63 SO.54 ) (2.9) S. I. Uni ts Q = 0..849 ) HW ) in which CHW is the HazenWilliams roughness coeffic,ient, S is the slope of the energy line and equals, hLP/L, R is the hydraulic radius defined as the crosssectional area, A, divided by the wetted perimeter, P, and for full pipes equals D/4 (where.D = diameter of pipe). Table 2.1 gives values for CHW for some commOn materials used for pressure conduits (Jeppson, 1977). Type of Pipe CHW PVC pipe 150 Very smootb pipe lLw New cast iron or welded steel 130 Wood, concrete 120 Clay, new riveted steel 110 Old cast iron, brick 100 Badly corroded cast iron or steel 80 Table 2.1 : Values of HazenWilliams Coefficient Equations (2.9) can be written in terms of hLP if Q is known. Thus, ES hLP 8.52 X lo5L Q1. 852 1.852 CHW D4 .87 with D in inches and L in feet. 10 PAGE 17 SI hLP IiII iO.7 L Cl.852 D4.87 HW Principles governing the flow of fluid as well as much imEmtal evidence indicates that the Flead loss due to ao.ded turbu11 lence or secondary flow in the presence of fittings, valves, meters and other components in a network, will be approximauruy to the square of the velocity or the flow rate squared. Minor losses are commonly expressed in the form hLM ... KM Q2 2 in which KM = M/(2gA ). (2.10 ) Nominal values of Mfor various common appurtenances are given in Table 2.2 (Jeppson, It is apparent from these loss coefficients that minor losses can be. neglected if relatively long pipelines are analyzed. However, in short pipelines, they may represent the major losses in the system, or if a valve is partly closed, its presence has profound influence on the flow rate. 2.4 PUMPS A number of alternative methods might be used to quantify the head, hpp produced by a ptunp. In some cases a constant power input is specified. In general, the relationship between pump head, hp and discharge. Q, can be expressed as (2.11) For a constant power pump, p(Q) = Z/Q (2.12 ) PAGE 18 r DEVICE : M Globe Valve (fully open) 10 Angle Valve (fully open) 5 Gate Valve (fully open) 0.19 Gate Valve (3/4 open) 1.0 Gate Valve (1/2 open) 5.6 Ball Check Valve (fully open) 70 Foot Valve (fully open) 15 Swing Check Valve (fully open) 2.3 Close Return Bend 2.2 Tee, Through Side Outlet 1.8 Standard Short Radius Elbow 0.9 Medium Sweep Elbow 0.8 Long Sweep Elbow 0.6 0 45 Elbow 0.4 Table 2.2 Loss Coefficients for Valves and Other Pipe Fittings 12 PAGE 19 where the pump constant, Z :: 550 HP/62.4 and HP :: useful pump horsepower. Other functions have been suggested, and a common choice is a second order polynomial of the form p(Q) = AQ2 T BQ + H o (2.13) in which A, Band Ho are constants for a given pump and might be determined by fitting Equation (2.13) to three points taken from a pump characteristic curve. 2.5 PRESSURE REGULATING VALVES A pressure regulating valve (abbreviated PRY) is designed to maintain' a constant, pressure downstream from it regardless of how large the upstream pressure is. Therefore, it is apparent that the unique relationship that exists between head and discharge, for line losses, minor losses and' pumps, does not exist for a PRV. Solution of pipe networks which include control with nonunique head discharge relationships using optimization techniques is still an active research area (Collins, Cooper, Helgason, Kennington, 1978). Exceptions to the above occurrence are; (1) If the upstream pressure,becomes less them the valve setting, and (2) if the downstream pressure exceeds the pressure setting of the valve so that if the PRY were not present, the flow would be in the opposite direction to the downstream flow direction of the valve. If the first condition occurs the valve has no effect on flow conditions. The PRV acts as a check valve, preventing reverse flow if the second condition occurs. By preventing reverse flow, the PRV 13 PAGE 20 14 allows the pressure immediately downstream from the valve to exceed its pressure setting. Thus, PRY's can perform both functions of reducing pressures in portions of a pipe distribution system if the pressures would otherwise be excessive, and may also be used to control from which sources of supply the flow comes under various demand.levels. In the latter applicationo the PRY acts as a check valve until the pressure is reduced to critical levels by large demands at which time additional sources of supply are drawn upon. The analysis of a pipe network containing PRY's must be capable of determining which of these conditions exist. 2.6 NODE ANALYSIS To obtain the system. of equations which contains the heads at the junctions/nodes of the network as 1 independent continuity equations are written as in Equation (2.). The relationship between discharge and head loss is then substituted into the continuity equations yielding a set of 1) 1 equations in Y) 1 unknown nodal heads. Solving for Q from the exponential formula (Eqtmtion 2.8). using double subscript notation, gives in which AHi j and K' lJ = (hLP)ij = (Kp)ij + (hLM)ij + (KM ) ij (2.14) Substituting into the junction continuity equations gives H. H. )l/n 1: 1 J ) d K ) lJ (2.15) PAGE 21 2.7 LOOP ANALYSIS If the discharge in each pipe is initially considered unknown instead of the head at each junqtion, the number of simul equations to be solved is .increased from 1) to 1 equations. However, this increase in the number of equations is somewhat compensated by a reduction in the number of nonlinear equations in the system. The analysis of flow in networks of pipes is based on the energy and mass conservation laws discussed in section 2.2. Mathematically, the continuity equations are concisely expressed as: Ar q ::: d where Ar is the reduced node incidence matrix. It is apparent that of these continuity equatiohs is linear. The remaining set of equations is formed by applying the energy conservation principle and expressed in terms of the fundamental loop matrix, B, as follows: BAh:: 0 which has ihdependent nonlinear equations. Having solved the system of equations for the discharge in each pipe, the head losses. in each pipe can be determined. From a known head or pressure at one junction, the heads and pressures at each imction throughout the network, or at any point along a pipe, can be determined by subtracting the head loss from the head at the upstream junction, and accouhting for differences in elevations if this be the case. In some problems the external flows may not be known. Rather the supply of water may be from reservoirs and/or pumps. The amount 15 PAGE 22 of flow from these indi vidu,al sources will not only depend on dem.ands but also will depend upon the head losses throughout the system. 2.8 CORRECTIVE MESH FLOW ANALYSIS This method of analysis yields the least number of equations. However, like the node analysis method, all the equations are nonlinear. These equations consider a corrective mesh flow as the unknowns and as discussed in section 2.2, the system of equations to solve is written as: 'I' q I:.li B P in which 15 is the mesh flow vector. Since there are I\. fundamental loops in a network, the corrective mesh flow system of equations consis ts of Il equations. This method requires an initialization of the flow in each pipe which satisfies all junction continuity equations. Since these initial flow estimates generally will not simultaneously satisfy the head loss equations. they must be corrected before 16 they equal the true flow rates in the pipes. A flow rate adjustment can be added with due regard for sign, to the .initially assumed flow in each pipe forming a loop of the network without violating continuity at the junctions. This fact suggests'establishing /\. energy equations around the of the network in which the initial flow plus the corrective mesh flow rate is used as the true flow rate in the energy equations. Upon satisfying these energy equations by finding the appropriate corrective mesh flow rates, PAGE 23 the 1 1 continuity equations would remain satisfied as they initially were. The corrective mesh flow rates may be arbitrarily taken positive in the clockwise or counterclockwise direction, but the sign convention must be consistent around any particular loop. * 17 PAGE 24 CHAPTER 3 NEWTONRAPHSON METHOD The Newton ... Raphson method is an iterative scheme which starts with an estimate to the solution and repeatedly computes better estimates. Unlike other methods which converge linearly, it has "quadratic convergence". Generally if quadratic convergence occurs, fewer iterations are needed to obtain the solution with a given precision than if linear convergence occurs. In addition to rapid convergence, the NewtonRaphson method is easily incorporated into a computer algorithm. Any of the three sets of equa tiOllS defining the pipe network problem, that is equations considering (1) the flow rate in each pipe lmlmovmp (2) the head at each jlmction unknown and (3) the corrective mesh flow rate around each loop unknown, may be solved 18 by this method. An initial guess is required for the Newton ... Raphson method. It is the method to use for larger systems of because it requiresless computer storage for a given number of equations. 3.1 APPLICATION TO NODE EQUATIONS The iterative Newton ... Raphson formula for a system of equations is, (mtl) ... (m) Dl F( (m) x x .... x (}.l) in which the superscripts.within parentheses are not exponents but denote number of iterations. The unknown vectors x and F replace the single variable x and function F and the inverse of the Jacob1 ;dF ian, D replaces 1 dx in the formula for solving a single equation. Adapting Equation (3.1) to solving the set of equations with PAGE 25 the heads as unknowns, Equation (3.1) becomes (m) D F(H ) (3.2) Making up' the Jacobian matrix D. are individual rows consisting of derivatives of that particular function with respect to the variables making up the column headings. For the system of head equations, the Jacobian is, 19 The Jacobian is a symmetric matrix and an algorithm for solving a linear system of equations with a symmetric matrix may be. preferred for greater computational efficiency. 3.2 APPLICATION CORRECTIVE MESH FLOW EQUATIONS The NewtqnRaphson method when applied to this set of equations becomeS in which the Jacobian is 'OFl aF1 "aF1 gp 1 8PZ ., .... 'OPL D 'dF2 .. 9F2 ()F 2 9Pl 9P2 aPL .. ., .. aFL 'aFL ,)FL aPl oP2 .. ., ... aPL PAGE 26 where L number of loops and P corrective mesh flovl for each loop. The NewtonRaphson method suffers from a setback of requiring a reasonably accurate initializationl otherwise it may not converge. When PRY's are present in a pipe network, the pro6edureof using identical loops for the corrective flow rates and energy equations must be altered. The reasons are (1) the head drop across a PRY cannot be expressed as a functioll of the piS circulating through that pipe, (2) continuity at some junctions will not be satis:fied if the are assumed to circulate through pseudo loops frolo artificial reservoirs created by the PRY's to reservoir in the network" The reason is that P in a pseudo loop would extract fluid from a junction, but not add an equal flo'N through another pipe joining at that junction. 20 Consequently, Some of the loops around which the energy equa tions are written cannot correspond to the loops around which the corrective flow rates, P, circulate. The individual piS will thus be assumed to circulate around the real loo}Js satisfying cantinui ty at all junctions. The energy equations will be written around loops containing pipes or other elements such as pumps or reservoirs whose head losses are functions of the discharge through them. PAGE 27 CHAPTER 4 LINEAR METHODS Nonlinearity of the ,function relating head and discharge is the crux of the,prob1em in solving a pipe network. Recall that in the loop analysis. there are .." + A. I equations of which A number of equations describirlg energy conservation around loops, are nonlinear. The other two analyses, namely the node analysis and corrective mesh flow rate analysis, each of which, having energy equations wri tten around each loop of the network; 'both involved nonlinear equations in each of its entire system of equations. Chapter 3 has dealt with the straightforward app1ica tion,of the NewtonRaphson Method to linearizing the nonlinear equations associated with the latter two methods. This chapter will be devoted to other linearization techniques, some of which are',variations of the NewtonRaphson Method. 4,.1 GRADIENT METHOD The gradient method, which is given extensive coverage by ,Wood (Wood, 1981), i$ derived from the first two terms of the Taylor series expansion. Any function, f(x). 'which is continuous, that is, differentiable, can be approximated as follows: (4.1) It is apparent from the 'right hand side of Equation (4.1) that the approximation has reduced f(x) to a linear form. However, if f is a function of more than one variable, Equation (4.1) can be generalized as follows: 21 PAGE 28 I ) = f(X(1)o,x(2)0"") + :!(l) (x(I) X(l)o) + (2)(x (2) x (2) 0) + (4.2) in which the partial derivatives are evaluated at some x(l) = 4.1.1 ALGORITHMS FOR THE SOLUTION OF LOOP EQUATIONS To conform to the notation: used in this chapter, Equation (2.3) which neatly describes the mass conservation (continuity) equation for each of the nodes in the network, is rewritten as follows: (j equations) (4.3) in which Qe denotes the external inflow or demand at the jUnction 22 I node, positive for inflows. The energy conservation equation in Equation (2.5) for fundamental loops without pumps, is now rewritten to include pumps as follows: (1 equations) (4.4) where hL := energy loss in each pipe, including minor losses; hp = energy input by pumps; t :: number of fundamental loops. For any t\VO fixed grade (or reservoir) nodes, the ,energy conservation equation written around this pseudoloop is written as I (f1 equations) (4.5) in ::: difference in total grade between two fixed grade nodes; f :: number of pseudo loops. If p equals the number of pipes in the netvIOrk. then the mass and energy equations form a set of p simultaneous equations of which (1 + f 1) equations constituting the set of energy equations nonlinear. *for networks with reservoirs PAGE 29 Using Equations (2.8), (2.10), (2.11) and (4.5), the energy equations expressed terms of the discharge, Q, are p(Q) (4.6) It can be seen that Equation (4.4) is a special case of Equation (4.6) where is zero for a fundamental loop. 23 Three algori thms are presently in sigriificant use and gradient method is employed to handle.the nonlinear terms in Equation (4.6). For a single pipe section, Equation (4.6) can be written a.s (4.7) which represents the grade difference across a pipe section ing flow Q. Substituting an estimate, Qi, Q, and denoting f(Qj) by Hi, Equation (4.7) becomes n :G Hi ;;: f(Qi) = KpQi + KM Qi P(Qi) (4.8) Differentiating Equation (4.7) and setting Q :=: Q., gives the 1 gradient of the function at Q ? Qi. Thus, nl f' (Q.)= nKpQ. t 1 1 Denoting f' (Q.) by G1 thus .1 Both the function and its gradient, evaluated at Q = Qi' will be used in all three algorithms for solving loop equations. 4.1.1.1 Single Path Adjustment (p) Method This method was first described Hardy Cross as the "Balancing Head Method 10 which was limited to .closed loopsysterns and included only line losses. The procedure is generalized and summarized as follows: PAGE 30 (i) An initial set of.flowrates which satisfies continuity at each junction is determined. (ii) A flow adjustment actor is computed for each path to satisfy the energy equation for that path and continuity must be maintained when applying the correction factor. (iii) Step (ii) is repeated using improved solutions until the average correction factor is within a specified lind t. Equation (4.6) is used to compute the adjustment factor for a path using gradient method to linearize the nonlinear energy equations. Thus, (4.10) in which .6. Q ::: Q Qi, where Qi is the estimated discharge. Applying Equation (it.lO) to Equation (4.6) and solving for 6Q gives 6 LHi LGi (4.11) which is the flow adjustment factor to be applied to each pipe in the path. The numerator represents the imbalance in the energy relationship due to incorrect flowrates and this procedure reduces this to a negligible quantity. Flow adjustment is carried out for all t fundamental (closed) loops and (fl) pseudo loops in the network. 4.1.1.2 Simultaneous Path Adjutment (SP) This algorithm is similar to the corrective mesh flow meJthod 24 PAGE 31 described in Section 2.8, the only difference is that gradient method is used here instead to linearize the energy equations. It is developed toimprove convergence by simultaneously adjusting the flowrate in each loop representing an energy equation. The method is summarized as follows: (i) An initial set of flowrates which satisfy at each junction node is (ii) A flow adjustment factor is simultaneously computed for each loop to satisfy the energy equations without disturbing the continuity balance. (iii) Step (i1) is repeated using improved solutions until the flow adjustment factor is within a specified limit. The simultaneous solution of 1 + f 1 equations is required to determine the loop flow adjustment factors. Each equation includes the contribution for a particular loop as well ascontributions from all other loops which have pipes common to both loops. For loop j, the head change required to balance the energy equation is expressed in terms of the flow change in loop j (6 Qj) and the flow changes in adjacent loops (.6 Qk) as follows: f(Q) :: f(Qi) + of of a Q L Q j + () Q Q k or, f(Q) = f(Qi) + f'(Qi) 1::::.Qj + f'(Qi) .6Qk (4.12) Substituting f(Q) :;;: f(Qi) OIl l":H i f'(Qi) :: LGi Equation (4.12) becomes 1::::. E L H" ::I (L G)." ) L Q. + Z (G.o 6 Qk) (4.13 ) ). J). "in which. LHi ::: SlIDl .of the head changes for all pipes in loop j ( E GOi) 6. Q j :: sum of all gradients for the same pipes times flow 25 PAGE 32 change for loop j. L (Gi 6Qk) :iii sum of gradients for pipes common to loops j and k multiply by the flow change for loop k. A set of simultaneous linear equa tionsis formed in terms of flow adjustment factors for each loop representing an energy equa The solution of these linear equations provides an improved solution for another trial until a specified convergence criterion is met. 4.1.1.3 Wgod's .. Linear (L) Method This method developed by Wood (Wood, 1981) involves the solution of all the basic hydraulic equations the pipe network. However, only the energy equations need to be linearized as the continuity equations are all linear. Using gradient approximation, the energy equations are linearized in terms of an approximate flowrate, Qi as follows: Introducing Hi and G. as before, the above equation becomes l Gi)Q= Z(GiQ i Hi) t 6E (4.14) This relationship is employed to formulate .... f 1) energy equations which together with the j continuity equations, form a set of p simultaneous linear equations in terms of the flowrate in 26 each pipe. One significant advantage of this scheme is that an arbitrary set of initial flowrates, which need not satisfy continuity, can start the iteration. A flowrate based on a mean flow velocity of h ft/sec has been used by Wood (Wood, 1981). The solution is then used to linearized the equations and successive trials are carried out until the change in flowrates between successive trials become insignificant. PAGE 33 4.1.2 ALGORITHMS FOR SOLVING NODE EQUATIONS TWo methods for solving the node equations are also widely used and are described here for completeness. 4.1.2.1 Single" Node Adjustment(N) Method This method was also first described in the paper by Hardy Cross and is known as the "Balancing Flows Method". The procedure is outlined as follows: (i) A reasonable" grade is assumed for each junction node in the$ystem. The better the initial assumptions, the the required trials. (ii) A grade adjustment factor for each junction node which tends to satisfy continuity is computed. (iii) Step (ii) is repeated using improved solutions until a specified convergence criterion is met. The grade adjustment factor is the change in grade at a particular node (.6. H) which will result in satisfying continuity and considering the grade at adjacent nodes as fixed. For conven"ience, the required grade correction is expressed in terms" ofQi which is the flowrate based on the values of the g!ades at adjacent nodes before adjustment. Thus, using gradient approximation, f(Q) ;;a f(Qi) + f'(Qi) .6.Q with the usual substitution, (4.15) where 6H H Hi' the grade adjustment factor and 6 Q denotes the flow corr"ee tions required to satisfy continuity at nodes. From Equation (4.3), 27 PAGE 34 .6 Q;lI4 E Qi Qe Thus,' from Equations (4.15) and (4.16), r:: Qi Qe L (l/Gi) 28 (4.16) (4.17 ) In Equation (4.17), ;'In,flow is assumed posi ti ve., The numerator represents the unbalanced flowrate at the junction node. Q. l the f10wrate in a pipe section prior to adjustment, is computed from ( .6 Hili<:) lin Qi :: in which 6. Hi grade change based on initial assumed values of grade. If pumps are included, the following expression is used t'o determine Q.: l n .6 H .. KQ P ( Q. ) l ].. l Equation (4.18) is solved using an approximation procedure. Adjustment of the grade for each junction node is made after each trial until a specified convergence criterion is satisfied, 11.1.2.2 Simultaneous Node 'Adjustment (SN) Method This method requires the linearization of the basic pipe network node' equations in terms of approximate values of the grade. If the discharge in Equation (lJ.. 3) is expressed in terms of the assumed heads, it can be written as: (4.19) for any node, a, and b denotes an adjacent node. Equation (4.19) can be linearized with respect td grades if the flowrates are Jri tten in terms of some initial values of the grades, Hai and Hbi' and the corrections in these grades. The gradient method is again used to calculate the flow.rate in pipe section, abo Thus, PAGE 35 Q '"" Q. L Ha oQ (Ll. 20) l + aHa t aH 6Hb b in which Q = ( (Ha Hb)/Kab)l/n (4.21) L:. Ha H Hai ( adjustment a factor for head at node a) 6Hb :: Hb Hbi ( adjustment factor for head at node b) Substituting the partial derivatives of the f10wrate expression in Equation (4.21) in Equation (4.20) and simplifying, gives In Qi nKab (4.22) The initial value of the flowrate, Q is computed based on the l ini tial values of the grades. r.rhus, where Kab may include minor losses, if any. Using Equation (4.22). the continuity equation for each junction 29 node can be expressed as a linearfunctioh of the variable and fixed grades of adjacent nodes and the variable grade of junction, a. Hence, N In Zv _.9L_ b=l nKab N L b;::l Hb H .a N = (4.23 ) l where N refers to all adjacent nodes, N v refers to adjacent variable grade nodes and NF refers to all adjacent fixed grade nodes. Qi is positive for outflow. Equation (4.23) is written for each junction node in the system resulting in a set of linear equations in terms of junction node grades. If pumps are included, two additional nodes may be PAGE 36 30 assigned to a pump at the suction and discharge sides as shovm schematically below: Qi a ... Qc ... d Two additional equations can be written: Kab (4,24) H Hb :a y(He H d ) a cd Hc Hb :::: [.. lin ] P HcHd)/Kcd) (4.25) Equation (4.24) is just the continuity equation and Equation (4.25) relates the head change across the pump to the flow in either the discharge or suction line. Equation (4.25) can be linearized using gradient method as follows: (4.25a) Using the gradient approximation, (4.26) Substituting the partial derivatives in Equation (4.26) and simplifying, we have the following linearized equation: (4.27) where and fo depend on the relationship used to describe the pump, p(Q), and are given by 0< :: P ( Q" ) Qi P' (Q" ) 1 n 1 = )/(nK 1 c 1 A set of (J" + 2N ) simultaneous linear equations (where N = p' p number of pumps) is generated.and solved starting'withQi's based on any assUmed set of jUl).ction node grades. An improved set of junction node grades is then used to compute an improved set of Qi's and the procedure repeated until a specified convergence criterion is PAGE 37 4.2 CmilMEN'rS ON ALGOHI'rHlvIS USING GHADIEN'r METHOD Node equations are easier to formulate because the include only contributions from adjacent nodes. On the other hand, the loop equations require the identification of an appropriate set of energy equations which include terms for all pipes in fundamental loops and between fixed grade nodes. Computer formulation of this set of equations is considerably more difficult than formulation of the node equations. Each of the procedures described is iterative in nature and computations terminate when a specified convergence criterion is met. The solutions are therefore only approximate although they can be very accurate. The ability of an algorithm to produce an acceptable solution is of prime concern and studies have demonstrated that convergence problems exist and an accurate solution is no t always possi lile. 4.2.1 ACCURACY OF SOLUTIONS A solution is considered accurate only when all the basic equations are satisfied to a high degree of accuracy_ For the three methods based on loop equations, the continuity equations are exactly satisfied. Each of these methods then proceeds to satisfy the energy equations iteratively and the unbalanced heads for the energy equations is evidence of solution accuracy. For methods based on node equations, iterations are carried out to satisfy continuity at junction nodes and the unbalance in continuity is a :;i n1t'icant indication of solution accuracy. 4.2.2. RELIABILITY OF ALGORITHMS A study carried out by Wood, using an extensive data base, 31 PAGE 38 has shown that the p. Nand SN methods exhibited significant convergence problems (Wood, 1981). Since these methods are widely used, great care must be exercised when using them. SN meth6dfailures are characterized by the inability tb meet a reasonable convergence criterioh and if this occurs in a limited nwnber of trials, further trials are usually of no benefit. Failure rate was quite high and iheuse of this method is not recommended unless a good accuracy is obtained in a reasonable number of trials. It has been established that algor.ithms based on node equations (N and SNmethods) failed to provide reliable results because of the inabili ty of these methods to handle low resistance lines. This is to the fact that solutionalgori thms for these equations do not incorporate an exact continuity balance. For each of the three methods singled out above, failure rates can be reduced if initial values closer to the correct values can be determined. However, this is rio easy task and as evidenced in the study, even an excellent set of initial conditions does not guarantee convergence. Both the SF and L methods provide excellent convergence and the attairlment of a reasonable convergence criterion is sufficient 32 to assure great accuracy. Convergence failure is 'very rare. However, since a gradient method is used to handle nonlinear terms., there is always the possibility of convergence problems. Illconditloned data such as poor pump descriptions are particularly prone. The L method has some advantages over the SP method. Assumed arbitrary flowrates need not continuity as the continuity conditions are already incorporated into the basic set of equations. PAGE 39 This method also allows a more straightforward and reliable inclusion of hydraulic components such as check valves, closed lines, and pressure regulating valves. Although the SP method has significantly less equations to the use of sparse matrix techniques to handle the larger matrix generated by the L method has somewhat negated this advantage. 4.3 LINEAR THEORY METHOD BY WOOD AND CHAHLES In this section, the linear theory method (Wood and Charles, 1972) will be described and used in the system of equations formulated by loop analysis which considers flowrates as unknowns (hereafter referred to as the Qequations). Like the other linear method described in Section 4.1.1.3. it has several distinct advantages over the NewtonRaphson or Hardy Cross methods. Firstly, it does not require an initialization, and secondly, according to Wood and Charles, it always converges in a relatively few iterations. However, its use in solving the head oriented or the corrective loop oriented equations is not recommended. Linear theory transforms the nonlinear loop equations into linear equations by approximating the head in each pipe by nl h ::: (KQ. ) Q :: K' Q L l (4.28) n1 in which Qi is an estimate of the flowrate, and"K' :; KQ. l Com1)ining these linearized loop equations wi th 'the jl junction continuity equations provides a system of p linear equations which can be solved by Gaussian elimination in conjunction with sparse matrix techniques (Tewarson, 1973). In applying the linear theory method it is not necessary to supply an initial estimate, as maybe implied. Instead, for the 33 PAGE 40 first iteration each K' is set equal to K, which is equivalent to setting all flowrates Q i equal to unity. In developing the linear theory method, Wood observed that successive iterative solutions tend to oscillate about the final solution. Reasons for the oscillation can be understood by observing that the linear theory method is a variation of the NewtonRaphson method described in Chapter 3 whereby K' in Equation (4.28) is simply the derivative of hL if multiplied by n. The oscillation could be prevented by multiplying each K by its n, which involves more computation than averaging consecutive solutions as proposed by Wood. Thus, the flowrate used in a trial is just the average flowrate fOT that pipe from the pl"evious D'IO solutions. or Qi(m) ::: [ Qi(ml) + Qi(m2) ] / 2 in which m within parentheses denotes a trial number. 4.3.1 INCLUSION OF PUMPS AND RESERVOIRS When pumps (not booster pumps) and reservoirs are connected to a network, the flows in the two connected, pipes become additional unknowns and therefore an additional equation is required beyond the j continui ty equations and .. /L fundamental loop equations. The additional equation is obtained from a pseudo loop, which connects the two reservoirs (fixed grade nodes) by a "no flow" pipe. If f fixed grade nodes exist in a network; there would be fl independent equations. Energy conservation around a pseudo loop (of which a fundamental loop is a special case) is defined by Equation (lj,. 6). Thus, .6. E ::: L (KpQn ... KIvlQ2) p(Q) PAGE 41 If the expression for p(Q) in Equation (2.13) is adopted in Equation (4.6), trte linear theory method does not give rapid convergence as it does when pumps and/or reservoirs are not present. A modification will therefore have to be made to allow the linear theory method to converge rapirlly. The reason for the modification is that the head produced by a typical centrifugal pump decreases u., nearly proportional to the reciprocal of the root of the flowrate whereas the head loss in a typical pipe increases approxi35 mately proportional to the square of the flowrate. A consequence of using this typical pump relationship in Equation (1+.6) is.. that if'the equation is solved by the linear methodt convergence may become very slow if at all. This situation can be improved by a transformation of variables so that the new unknown has an exponent close to n. Such a transformation is G :;:: Q + B/2A in whichGis the new variable and A and B are the $ame constants in Equation (2.13). The appropriateness of Equation (4.29) is demonstrated by solving it or Q and substituting in Equation (2.13). After some simplification, hp ;:; AG2 + ho (4.30) h :::l H B2/4A o 0 where Obviously, the exponent of G {that is, 2) is close to the typical n. Substituting Equation (4.30) in EquatiOl1 (4.6) gives 22' Z (KpQn t KMQ ) _L AG ::.6. E t Lho (4.31) Addition of Equations (4.29) and (4.31) produces a system with as many equations as unknowns. PAGE 42 4.3.2 INCLUSION OF PRESSURE REGULATING VALVES (PRV'S) Networks containing PRV's may be analyzed by the linear theory method by initially assuming that the pressure (or head) immediately downstream from a PRY is constant and equal to the valve setting. Junction continuity equations are then written as if no PRY's are present. To write the loop equations, pipes containing PRV's are disconnected fr6m the upstream nodes and the PRY's are replaced by dummy reservoirs. After each iteration a check on the flowrate Q, in each pipe containing a PRY is made. If there is any negative Q, the pseudo loop equation which includes terms for that pipe is modified with Q replaced by an unknOwn grade (head) immediately downstream from that PRY. 36 PAGE 43 & COMPUTATIONAL EXPERIENCE 5.1 WIA THEMATI CAl, PH OGRAMMING TECHNI QUES As a prelude to introducing the al ternative to solving the governing equations for a pipe network using optimizatior: techniques v. it is convenient to define a network topology using notations which are consistent with those used in graph theory. Let the network topology be described by a node set Nand an arc set EQ In each of the set Eo' let Qij denote the flowra te from node ito. j. Each node, n in the set N is associa ted wi th a hydraulic head, Hn. Let R, 2. subset of N, be the set of nodes corresponding to reservoirs (fixed grade nodes) and let for all n f.I R be the fixed head associated wi th a reservoir. Also let rn for all n S N denote the flow requirements (that is. supply or demand) at node n. For an fluid, the governing network equations can be stated as: Q;;;: "' in all n & N (i,n) uE H n 0 all (i,j)BEo (5.3) ::::: WC all noR n (5. h) Equation (5.1) is just a tement of maSi; conservation at each node while Equation (5.2) stipUlates mass conservation for the network as a v/hole. Equation (5.3) states that the head loss Hi Hj o::6Hij 37 PAGE 44 across an element is some function F .. of the discharge through the lJ element while Equatibn (S.4} requires that at a reservoir node, the head is constant. The functional form of Fij or its inverse Eij ( 6 H .. ) =: Q is not specified and can represent any element inlJ lJ cluding simple pipes minor loss devices as long as a unique relationship between head and discharge exists. In general" Fij most or all (i,j) in Eo is nonlinear, thus necessitating iterative techniques as (i) Hardy Cross, (ii) NewtonRaphson, and (iii) to be used to solve .the governing network equations. Most of these techniques are detailed in Chapters J and 4. Each of these methods is simply a technique for solving a set of nonlinear simultaneous equations which have been adapted to the network analysis problem. Each is iterative in nature and begins with an initial trial solutioll. A new solution is obtained by solving a set of linear equations using straightforward procedures. If the new solution differs from the trial solution by less than 8. specified amount then the iteration stops. Otherwise, the new solution becomes the trial solution and the procedure is repeated. In some of the algorithms, an initial trial salution sufficiently close to the true solution is required td ensure convergence. The differences in the methods result from. the use of different strategies to determine the neW solution. The new approach by Colljns, Cooper, Helgason and Kennington (1978) represents a radical departure from the state of the art iterative methods as optimization models are employed to the network problem. Two alternatives models are formulated and these models play analogous roles to the node versus loop formulations 38 PAGE 45 for solution of the network equations used in the state of the art methods. The first. of the two optimization models, .called the Content Model assumes the form Minimize G =: Subject to L Q.:::: rn all n&NU(g) (i n) E> EUEI 1n 39 Qij 0, all (i,j) /EoU(EI ) 1n which E is the arc set for a.network in which the arcs have been replaced by two equivalent oppositely directed.oneway elements so that Qij can assume only positive values. This replacement is done as a mathematical convenience so that Qij can be treated as a constrained decision variable in the optimization model which has a non.negativity condition imposed on Qij. It can be proven that the solution obtained by solving the E network will produce identical results as those which would be obtained by solving the original Eo network which permits Q i j to be unconstrained. The arc set El i:3 merely a set of arcs c'onnecting all nodes in N to a ground node g and is introduced to satisfy mass conservation fhr the network as a whole {Equation (5.2). Using the. terminology of Cherry and Millar (1951), the above is to find a set of flows which satisfies flow conservation and minimizes system content, G, .hence the name Content Model. PAGE 46 The second optimization model, the CoContent Model, is a complementary (but not dual) model which has the form Minimize J ) 0E ( t ) d tJ' 0 lJ subject to I:::. H. + 6 H 6 H. :: 0.. all (i, j ) S E 1 J J g 19 /\ H := IP all n S R ng n g' Ih the terminology of Cherry and Millar (1951), the above problem is to find a set of head losses which sums to zero around all loops and minimizes system CoContent, J, hence the name CoContent Model. Using KuhnTucker theory (Kuhn and Tucker, 1950), it can be proved that the solution to either of these models yields the solution to the pipe network problem, that is, the optimal solution satisfies the governing network equations. The proof is carried out by examining the derivatives of the objective function and showing that the derivative conditions for a stationary point, along with other constraintsp are identical to the network equations. In the proof it is assumed that the Fij and Eij functions are monotonically increasing. This assumption insures the convexity of the objective function which in turn guarantees the existence of a unique solution to the optimization problem. The mono tonicity of Fij and Eij merely implies the fact that energy losses in a network element increase with increasing discharge The Content Model has the special structure of a convex cost network flow problem for which efficient routines are available., Numerous nonlinear algorithms such as (i) FrankWolfe method. 40 PAGE 47 (ii) piecewise linear approximation and (iii) convexsimplex method are available for solving such a problem. The use of mathematical programming techniques in pipe network analysis has paved the way for potential research in the following areas = (i) Extension of mathematical programming techniques to solution of compressible flow pipe network analysis problem. 41 (ii) Incorporation of time variable storage in network elements to solve transient network problems. (iii) Use of mathematical programming techniques. to solve complex open channel networks. (iv) Feasibility of using mathematical programming techniques to solve netwox'k parameter identification problems such as the headdischarge relationship in pipe network analysis. (v) Development of an economTc model to minimize the operational costs for a flov" network with operational behavior given by one or more network problems. 5.2 COMPUTATIONAL EXPERIENCE A computer program was written based on the linear theory method (Wood and Charles, 1972) described in Section 4.3. The program was designed to solve the system of loop and node equations using the ,i tera ti ve procedure described by the method. Two features this FORTRAN computer program may have for application include: (i) the capability of handling networks containing pumps and reservoirs, and (ii) an algorithm which analyzes networks containing pressure regulating valves. PAGE 48 The use of the node incidence matrix and fundamental loop matrix described in Section 2.1 in the algorithm has provided an efficient means of translating information contained in any pipe network into a network simulator. Incidentally, the node equations and loop equations were formulated using the node incidence matrix and fundamental loop matrix respectively. In carrying out all computations, friction losses in pipes were assumed to be described by the HazenWilliams equation and pumps were described by the quadratic form (Equation 2.1]). The convergence criterion employed was: Z I Qi Qill t' t 'Qi r < 0.0005 in which Qi is theflowrate obtained for a trial and Qil is the flowrate obtained from the preceding trial. This appears to be a stringent requirement which may assure good accuracy if the condi tion_is satisfied. However accuracy is achieved if and only if continuity at every node and the energy equations are exactly satisfied. A small scale network, from Jeppson (1977, p.169) shown in Fig 5.1, was tested. This 8 pipe, 5 node netwol"'k, with the properties given in Table has 2 reservoirs, a pump and a pressure regulating valve. The solution to the test problem and the solution reported by Jeppson (1977, p.lIO), using the same theory are tabulated in Table 5.2. The results appear to be in good agreement. 42 PAGE 49 180 f t 11 2.0 cfs (1) ft 1 (it) (2) 200 ft 8) (3) (6 ) [4] [3 1.0 cfs Fig 5.1 Test Problem 1 Pipe Length (ft) Diameter (in) HazenWilliams Coefficient 1 2 3 4 5 6 7 8 1000 6 800 6 1000 6 800 6 1200 6 1000 6 500 8 500 8 Table 5.1 Network Parameters Discharge (cfs) Head (ft) 1.0 40.0 1.5 35.0 2.0 26.0 ., Pump Characteristics 110 120 110 120 120 1?:0 130 130 43 2.0 cfs PAGE 50 Table 5.2 Solution to Test Problem The detailed solution and program listing are contained in the Appendix. With this program, the test problem took 0.12 second of execution time on an Amdahl 470 computer The number of iterations required to meet the convergence criterion was 6. The subroutine used for solving the linearized set of loop equations and the linear node equations simultaneously was developed based on the Gaussian method of elimination improVed by pivotal condensation (Tewarson, 1973). 44 The capability of the program to handle a larger network has not been proven but it would have stretched the available storage of a computer to its limits if it has been tested. Storage space is primarily taken up by the final augmented matrix which comprises essentially the node and loop equations. The use of sparse matrix techniques instead of full matrix methods may extenrl the capability of the program to analyze larger networks of a few hundred pipes and nodes. PAGE 51 CHAP'I'ER 6 CONCLUSION The Hardy cross method which sparked off the evolution of the numerous techniques of simulating pipe networks, is suitable only for relatively small networks. With the advent of the computer. and as larger and more complex networks were the Hardy Gross method was found to frequently converge too slowly if at all. The classic method which is described in most hydraulics or fluid mechanics text is an adaptation of the NewtonRaphson method which solves one equation at a time before proceeding to the next equation during each iteration instead of solving all equations simultaneously. The single path and single node methods described in Sections 4.1.1.1 and 4.1.2.1 respectively, 2re basically the classic Hardy Cross methods. Procedures develolled to improve the convergence of the single path method were described by Martin and Peters (1963) and later by Epp and Fowler (1970). The procedure involves the simultaneous computations of flow adjustments and was presented j_n Section 4.1.1.2. A similar approach has been developed for the node equations where all node equations are linearized and solved simultaneously. This method is described by Shamir and Howard (1968). All of the four methods mentioned so far require an initial guess as to the solution and the rate of convergence depends to 45 a degree on how close this initialization is to the correct solution. For the system of aqua tions which is flowrate oriented, two linearization teChniques (Wood, 1981 and Wood and Charles, 1972) were describe4 in Sections 4.1.1.3 and 4.3 respectively. Both of these procedures do not require an initialization and have been reported to converge in a relatively few iterations. PAGE 52 Significant ,convergence problems were reported for the Single Path, Single Node and Simultaneous Node (Wood, 1981). It has been suggested that if a specified stringent convergence criterion cannot .be met using single adjustment methods. the solution is probably unreliable. For the simultaneous node adjustment method, it has been suggested that the best indication of an acceptable solution is that the average relative unbalanced flow at the junction nodes be less than 2%. Instances of failures have also been reported in cases where line losses vary greatly or pumps operate on steep curves even when 'good initial approximations are available. 46 The simultaneous path methods and the linear method using gradient approximations, were reported to provide excellent convergence and the attainment of a stringent convergence criterion is sufficient to assure great accuracy in most' cases. In the study carried out by Wood (1981), in which a wide variety of situations was represented, some incorporating features which increase gence difficulties like low resistance lines; these methods were reported to attain accurate solutions in a relatively few iterations. However, if gradient approximations are used to handle nonlinearity, convergence problems are always a possibility, especially if illcondi tiqned data such as poor pump description's are employed. Of all methods, the linear methods developed by Wood and Charles (1972) and a later version by Wood (1981), who used gradient approximations, offer more advantages. A balanced initial set of flowrates is not required since the, continuity conditions are already incorporated 'into the basic set of equations. These algorithms permit a more direct and reliable incorporation of hydraulic components such as check valves, closed lin,es and pressure regula ting valves. For any pipe network simulators to be of general use, PAGE 53 these components, which affsct contjnuity, and their effects on the hydraulics of the network must be incorporated into the basic set of equations. However, the set of equations solved by the linear methods involved significantly more equations which will be a setback if full matrix methods are used. The use of sparse matrix techrtiques has somewhat corrected this disadvantage and has rendered it a more desirable algoritrw to adopt for analysis of pipe networks. The use of mathematical programmihg techniques in pipe network analysis holds a lot of promise for the future. One of the direct consequences of the theorydesriribed in Section 5.1 is the identification of a unitary measure by which the goodness of a solution can be gaged. Traditional methods described previously give no good insight into the goodness of an approximate solution, for large scale problems. The remove the vagueness that inherently surrounds a definition of 47 It close II when an attempt is made to utilize a comparison of individual flows, heads, or in individual elements. Optimization techniques also have their setbacks. One is that functions describing friction losses, minor losses in pipes and pump heads must necess arily be convex functions for a.solution to be guaranteed. In addi tion, head loss must be a unique function of discharge. Such. uniqueness may not exist for certain control elements such as check valves and pressure regulating valves ,Until these problems are resolved, its application will be limited in scope. PAGE 54 REFERENCES (1) Adams R.W., Distribution ,analysis by electronic computer, Institution of Water Engineers, 15,' 415428, 1961. (2) Bazaraa M.S. and J,J. Jarvis, Linear programming and network flows, J. Wiley & Sons, New York, 1977. (J) Belevitch V. Classical netwQrk theory, HoldenDay, Inc., 'San Francisco, 122,1968. (4) Bellamy C.J., The analysis of networks of pipes and pumps, Journal Insti tutiori. of Engineers, Australia, 37(4 ... 5), 111116, 1965. 48 (5) Bree D.W.,Jr.g Cohen, ,L.C. Markelp & H.S. Rao, Studies on operations planning and control of water distribution systems, Final Technical Report to OWRR Project No. 5214, Systems Control, Inc., Palo Alto, Calif., 1975. (6) Brock D., Closed loop automatic control of water system operations, Journal American Water Works Association, 55(4), 467480, 1963. (7) Brock Metropolitan water system operation subsequent to nuclear attack or natural disaster, Report No. AD 711956, Dallas Water Utilities, City of Dallas, Texas, 364 pp, 1970. (8) Chenoweth H., and C. Crawford, Pipe network analysis, JOl.trnal American Water Works Association, 66(1), ,5558, 1974. (9) Cherry E.C. and W. Millar, Some general theorems for nonlinear systems possessing resistance, Phil. Mag., (Ser 7), 42(333) 1951. (10) Clay R., Nonlinear networks and systems, Wiley Interscience, New York, 139, 1971. PAGE 55 ,(II) Collins A.G.t and R.L. Johnson, Finiteelement method for water distribution networks, Journal American,Water Works Association, 67(7), 385389, 1975. (12) Collins N.A., L. Cooper, V. Helgason, and J.L. Kennington, Solution of largescale pipe 'networks by improved mathematical approaches, lEOR 77016o/R 77001, 1978. 49 (13) Di1lingham,J.H., Computer analysis of water distribution systems, part 1, Water and Sewage Works, 114(1), 13; part 2, 114(2), 4345; part 4, Program application, 114(4), 141142, 1967. (14) Donachie R.P., Digital program for water network analysis, Journal Hydtau1ics Division, Amer. Soc. Civil Engineers, 100 (HY3), 393403, 1973. (15) Eggener c. L., and L. B. Po1kowski, Netviork models and the impact of modeling assumptions, Journal American Water Works Association. 68(4), 189196, 1976. (16) Epp R., and A.G. Fowler, Efficient code for steadystate flows in networks, J. Hydra'ulics Di v., Proe. Amer. Soc Civil Engineers, 96 (HYl), 4356, 1970. (17) Fietz T.R., Discussion of "Hydraulic network analysis using ,linear theory", Hydraulics Div., Proc. Amer. Soc. Civil Engineers, 99 (HT5), 855857, 1973. (18) Fietz T.R., Discussion of "Efficient algorlthmfor distribution networks", J. Hydraulics Div., Proc. Amer. Soc. Civil Engineers, 99 (HYIl), 21362138, 1973. (19) Gerlt J.1J., and G.F. Haddix, Distribution system operation analysis model, J. Amer. Water Works Association, 67(7) 381384, 1975. PAGE 56 (20) Graves O.B and D. Branscome, Digital computers for pipeline network analyses, J. Sanitary Engineering Div., Proc. Amer. Soc. Civil Engineers, 84 (SA2), 118, r958. (21) Guillemin E.A. Introductory circuit theory, John Wiley & Sons, Inc.r New York, .510, 1953. 50 (22) Hoag L.N., and G. Weinberg, Pipeline network analyses by electronic digital computer, J. American Water Works Association, h9(5), 517524, 1957. (23) Hudson VI.D., Computerizing pipeline design, J. Transportation Engineering Div., Proc. Amer. Soc. Civil Engineers, 99(TEl) 7382, 1973. (24) Hudson W.D., A modern metroplex looks ahead, J. Transportation Engineering Div., Proc. Amer. Soc. Civil Engineers, 100(TE4) 801814, 1974. (25) Jeppson R.W., Analysis of flow in pipe networks, Ann Arbor Science, 1977. (26) Karni S., Intermediate network analysis, Allyn and Bacon, Inc., Boston, Mass., 84113, 1971. (27) Kuhn H.W., 2,nd A.W. Tucker, Nonlinear programming., in Proceedings of the Second Berkeley Symposium on Mathematical StatistIcs and Probahility, Berkeley, 481492, 1950. (28) Lam C.P., and M.L. Wolla, Computer of water distribution systems, part 1 formulation of equations, J. Hydraulics Div., Proc. Amer. Soc. Civil Engineers, 98 (HY2), 335344, 1972. (29) Lam G.P.; and M.L. Wolla, Computer analyses of water distribution systems, part 11 numerical solution, J. Hydraulics Div., Froc. Amer. Soc. Civil Engineers, 98 (HY3), 447460, 1972. PAGE 57 (30) Lemieux P. F., Efficient algorithm for distribution networks, J IHydraulics Div., Proc. Amer. Soc. Civil Engineers, 98 (HYll), 19111919,1972. (Jl)Liu K.T., The numerioal analyses of water supply networks by digital oomputer, Thirteenth Congress of the International Association for Hydraulic Research, 1, 3643, 1969. (32) Marlow T.A. R.L. Hardison, H. Jacobson and Beggs, Improved design of fluid networks with compl,1ters, J. Hyd draulids Div., Proc. Amer. Soc. Civil Engineers, 92 (HY4), 'n6l, 1966. (33 )l\1artin D. W. and G Peters. The applica t.ion of Newton is method 51 to network by digital computer,Institution of Water Engineers, 17(2), 115129, 1963. (34) McIlroy M.S.# Pipeline network flow analyses, J. Amer. Water Works Association, 41, 422428, 1949. (35) McPherson E.C. Bolls, Jr., D.A. Brock, E.B. Cobb, H.A. Cornell, J.E. Flack, F. Holden, F.P. Linaweaver, Jr., R.C. McWhinnie, J.C. Neill, and R.V. Alson, Priorities in distribution research and applied development needs, J.Amer. Water Works Association,. 66(9), 507509, 1974. (36) Rao H.S., D.W. Bree, Jr., and R. Extended period simulation of water distribution networks; Final technical report OWRR project no. C4l64, Systems Control Inc._. Palo Alto, Calif., 120 pp., 1974. (37) Rao H.S., D.W. Bree, Jr., and R. Benzvi, Extended period simulation of water systems part A, J. Hydraulics Div., Proc. Amer. Soc.divil Engineers, 103 (HY3), 97108, 1977. PAGE 58 (38) Ra6 H.S., L.C. Markel and D.W. Bree, Jr., Extended period simulation of water systems B, J. Hydraulics Div., Proc. Amer. Soc. Civil Engineers, 103 (HY3), 281294, 1977. (39) SharrlirU., Via ter distribution systems analysis, IBM Research Report RC 4389 (No. 19671), IBM Thomas J. Watson Research Yorktown Heights, New York, 1973. (40) Shamir U., and C.D. Howard, Water distribution system analysis, J. Hydraulics Dfv., Proc. Amer. Soc. Civil Engineers, 94 (HYl),219234, 1968. (41) Shamir W., Optimal design and operation of water distribution systems, Via tel' ResourcesRes., 10 (1). 2736, 1974. (1.12) ':r'ewarsonR.P., Sparse Matrices, Academic Press, Inc., 1973. (43) Williams G.N., Enhancement of'convergence of pipe network solutions, J. Hydraulics Div.,Proc. Amer. Soc. Civil Engineers, 99 (HY7) 10571067, 1973. (4lt) Wood D.J., and C. Charles, Hydraulic network analysis using linear theory, Hydraulics Div., Proc. Amer. Soc. Civil Engineers, 98 (HY7), 11571170, 1972. (45) Wood DolT. ,An explicit friction factor relationship, Civil Engineering, 36(12), 6061, 1966. (46) Wood D.J., Algorithms for pipe network analysis and their.reliability, University of Kentucky, Water Resources Research Institute Research report No. 127, 1981. 52 (47) Zarghamee M.S., Mathematical model for water distribution systems, J. Hydrattlics Di v., Proc. Amer. Soc Civil Engineers, 97 (HY1), 114, 1971. PAGE 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 $JOB A.PPENDIX 5' c c c c C c C C C C C C C C C C C C C C C DIMENSION HWC(100',XL(100),DIA(100), XKM(100),RFLOW(60', *IN ( 60, 10) I NN ( 60 ) I NVAL V ( 10) I KP ( 10) A ( 10) I B ( 10) I HO ( 10) I *LP ( 40, 20), NLP ( 40 ) I H2 ( 40 ) I HI (40) G ( 100 I 12), G ( 10, 12) J *FHH 110. 111 ), XKT (100) I X (100), KPRVS ( 10) I HGL (10} I *JAC(10), HP(10,3), GP(10.3), JNL(40,), NIK(40).NB..1(10), *XLPRV(10),HD(60),HGLU(10).VALST(10) NOTE: FOLLOWING DATA ARE FORMATTED (713) NP==NQ. OF. PIPES (MAX=100, COL 13); NJ=NO. OF NODES (MAX=60, COL 46);NL=TOTAL NO. OF LOOPS (MAX=40,COL 79); NPUMP=NO. OF PUMPS (MAX=10. COL 101.2) .. NRLP=NO .. OF REAL LOOPS(LOOPS CON TAINING NO RESERVOIRS/PRESSURE REGULATING VALVES(PRV'S},COL 1315); NPRV=NO. OF PRV'S(MAX=10, COL 1618)iMAX=NQ. OF ITERA TIONS (DEFAULT=8, MAX=12 .. COL 1921) READ I, NP I NJ, NL NPUMP I NRLP, NPRV, MAX 1 FORMAT (713) IF (MAX) 3.3,2 3 MAX=8 2 DO 5 I=l,NP X(L}=LENGTH OF PIPE 'I' IN FT (FOR PIPES WITH PRV'S,DOWNSTREAM LENGTH IS READ INSTEAD);DIA(I)=DIAMETER OF PIPE 'I' IN INCHES; HWC(I)=HAZENWILLIAMS COEFFICIENTS;XKMCI)=MINOR LOSS COEFFICIENTS FOR VALVES, BENDS & OTHER FITTINGS IN PIPE I'. 5 READ, XL(I), DIAn), HWC(!), XKM(I) DO 7 1=1. N..1 RFLOW(I)=DEMAND(VE}/SUPPLY(+VE) AT NODE 'I'iNNJ=NO. OF PIPES MEETING AT NODE 'I' (MAX=10)iJN(I,J)=I/D NUMBERS OF PIPES MEETING AT NODE 'I' (ANY ORDER) AT NODE 'I' (CAN BE ARRANGED IN ANY ORDER) READ. RFLOW(I',NNJ, (JNCI.J).J=I,NNJ) 7 NNCI)=NNJ DO 1000 I=l.NPRV NVALV(I}=I/D NUMBERS OF PIPES WITH PRV'S (FOR IDENTIFICATION,PRV'S TO BE NUMBERED CONSECUTIVELY FROM 1 & NUMBERS OF PIPES READ IN SAME ORDER AS PRV'S ARE CONSECUTIVELY NUMBERED) NBJ(I)=NODE UPSTREAM OF PRV; XLPRV(I}=PIPE LENGTH UPSTREAM 1000 OF PRV (FT), VALST(I)=PRV SETTING (FT) READ, NVALV ( 1) NBJ (I ) I XLPRV ( 1) VALST ( I) DO 6 I=l,NPUMP c c c C C C C C C C C C C C C C C C C NOTE: FOLLOWING DATA ARE FORMATTED (I3,3F8.3) KP(I}=NUMBER ASSIGNED TO PIPE 'I' WHICH CONTAINS A PUMP (COL 13); A(I)(COL 411>.B(I)(COL 1219),HO(I)(COL 2027}=PUMP CONSTANTS IN THE EGUATION, PUMP HEAD,HP=A*G**2+B*Q+HO (COLS ASSIGNED TO A,B & HO TO BE LEFT BLANK IF HP/S & GP'S ARE SPECIFIED LATER IN PROGRAM). READ 11, KP ( I ) A ( I ), B ( I ) I HO ( I ) 11 FORMAT (I3.3F8.3) IF (A{I 6,9,6 9 DO 15 J=1,3 HP(LJ)=PUMP HEAD IN FT. CORRESPONDING TO PUMP IN PIPE NO. KP(I); GP(I,J)=CORRESPONDING DISCHARGE IN CFS (3 SETS OF HP & GP TO BE SPECIFIED. EACH TO A LINE) (NEED NOT BE SPECIFIED IF A.B & HO ARE SPECIFIED EARLIER IN PROGRAM>' READ, HPCI,J),GP(I,J) FIM(J.1)=GP(I,J)**2. FIM( J. 2)=GP ( 1, J) FIM(..1,3)=1. 15 FIM(J.4)=HPCI,J) CALL GAUSSCFIM,X.3) An )=X (1) B(I)=X(2) HO(I)=X(3) 6 CONTINUE NAPC=NP+NPUMP+l NJLR=NJ+NL+NPUMP DO 8 I=l,NL IF (INRLP) 10,10,12 NNLP=NO. OF PIPES FORMING A REAL LOOP (MAX=20)iLP(I,.J)=NUMBERS OF PIPES IN REAL LOOP 'I', ARRANGED IN CLOCKWISE ORDER, +VE IF CLOCKWISE & VE, OTHERWISE. 10 READ, NNLP, (LP(I.J},J=l,NNLP) GO TO 14 NNLP=NO. OF PIPES FORMING A PSEUDO LOOP (LOOPS CONTAINING RESERVOIRS/PRV'S. MAX=20) i LP (I, ..1)=1/0 NOS. OF PIPES IN PSEUDO LOOP I' ARRANGED IN CLOCKWISE ORDER. +VE IF CLOCKWISE, & VE, OTHERWISEi H2,Hl=RESERVOIR/PRV HEADS IN LOOP 'I',ARRANGED IN CLOCKWISE ORDER ALONG PATH CONNECTING RESERVOIRS/PRV'S. 12 READ. H2(INRLP). Hi< INRLP) I NNLP, (LP (1, J) I J=1/NNLP) 14 NLP ( !) =NNLP 8 CONTINUE DO 1100 I=1,NL NNJN=NO. OF NODES IN LOOP (MAX=20)iJNL(I,J)=I/D NOS. OF NODES IN LOOP ARRANGED IN SAME ORDER AS PIPES IN THE LOOP IJITH EACH PAGE 60 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 C NODE I/D NO. READ IN AFTER A PRECEDING PIPE I/D NO. READ, NNJN. (JNL(I.J),J=1.NNJN) 1100 NIK(I)=NNJN LX=O 16 DO 18 I=l,NP 18 G( I,l }=1. DO 20 I=l.NPUMP 20 GU I 1>=1. KB=1 500 KB=KB+l DO 25 I=l,NJLR DO 25 J=l,NAPC 25 FIM( I I J )=0. C FORMULATION OF NODE INCIDENCE MATRIX DO 30 1=1. NJ NM=NN(I) DO 30 J=l,NM MN=IABS(JN(I,J IF (JNCI.J 26.28,28 26 FIM( I, !"IN) =1. GO TO 30 28 FIM( I, MN)=1. 30 CONTINUE C FORMULATION OF FUNDAMENTAL LOOP MATRIX DO 40 1=1, NL MN=NLP(I) DO 40 J=l.MN ML=IABS(LP(I,J IF (LP(I,J 33,37/37 33 FIMCNJ+I,ML)=l. GO TO 40 37 FIM(NJ+I. 1"1'.)=1. 40 CONTINUE DO 60 I=1.NL IF (INRLP) 44.44.42 42 FIM(NJ+I,NAPC)=H2(INRLP)Hl(INRLP) 44 LN=NLP(I) DO 62 J::::l/LN KPA=IABSCLPCI,J DO 61 K=l.NPUMP IF (KPAKP(K}) 61146,61 46 IF (LP(r,J 48,50,50 48 HPS=B(K)**2. /(4. *A(KHO(K) 49 FIM(NJ+r,NAPC)=FIM(NJ+I,NAPC)+HPS GO TO 62 50 HPS=HO(K)B(K)**2. /(4.*A(K GO TO 49 61 CONTINUE 62 CONTINUE 60 CONTINUE DO 65 I=l.NJ 65 FIM(I.NAPC)=RFLOW(I) C COMPUTATION OF LINE & MINOR LOSSES DO 70 1=1. NP 70 XKT(I)=(S. 52E5*XLCI)*CABS(GCI,KBl)**0.852) */(HWC(I)**l. 852*DIA(I)**4. 87)+8. *XKM(I>*ABS( *G( I, KB1/ (32.2*3. 141593**2. *DIA( I) **4. ) DO 75 I=l,NL DO 75 J=l,NP 75 FIMCNJ+I.J)=FIM(NJ+I,J)*XKT(J} NJL=NJ+NL DO 100 I=l.NL LM=NLP(I) DO 100 J=l.LM KPA=IABS(LPCI,J DO 100 K=l,NPUMP IF (KPAKP(K 100.80,100 80 IF (LP(I.J 82,84,84 82 FIM(NJ+I,NP+K)=A(K)*ABS(G(K,KBl GO TO 86 84 FIM(NJ+I, NP+K)=A(K)*ABS(GCK,KBl 86 FIMCNJL+K,KPA)=l. FIM(NJL+K.NP+K)=l. FIM(NJL+K,NAPC)=B(K)/(2. *A(K 100 CONTINUE IF (LX) 102,102,200 102 CALL GAUSS(FIM,X.NJLR} GTOT=O. GFLCH=O. DO 105 1=1. NP GO, KB)=X( I) GTOT=GTOT+ABS(X(I 105 GFLCH=ABS(GCI,KB)QCI,KBl+GFLCH ERR=GFLCH/GTOTO. 0005 PAGE 61 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 13'7 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 DO 108 I=l,NPUMP 108 G(I.K9>=X(NP+I) IF (ERR) 110, liD. 130 liO DO 115 I=i,NPRV KPRV=NVALV(I) IF (CHKPRV,KB 112,115,115 112 LX=1 KKB=KB GO TO 16 115 CONTINUE GO TO 900 130 IF (KBMAX) 140,140,150 140 IF (KB2) 500,500,141 141 DO 142 I=l,NP 142 GCI,KB)=(G(I,KB)+GCI,KBl/2. DO 145 I=l,NPUMP 145 GCI.KB)=(G(I,KB)+G(I,KBl/2. GO TO 500 150 PRINT 155 155 FORMAT (2X, 'DESIRED ACCURACY CANNOT BE ATTAINED") PRINT 160 160 FORMAT (2X. 'IN NO. OF ITERATIONS SPECIFIED 'I) PRINT 165, MAX 165 FORMAT (2X, 'NO. OF ITERATIONS SPECIFIED = ',12/) PRINT 170, ERR 170 FORMAT (2X, 'ERROR FLO. 511) GO TO 900 200 LX=O DO 250 I=l,NPRV KPRV=NVALV(I) IF (G(KPRV,KKB 210,250,250 210 LX=LX+l KPRVS(LX)=KPRV NRLP1=NRLP+i DO 240 J=NRLP1,NL LJ=NLP(J) DO 245 JJ=1,LJ IF (KPRVIABS(LP{J,JJ) 245,213,245 213 IF 215,217,217 215 FIMCNJ+J.NAPC)=FIM(NJ+J,NAPC)H2(JNRLP) FIM(NJ+J,KPRV)=l. GO TO 245 217 FIM(NJ+J, NAPC)=FIM(NJ+J, NAPC)+Hl(JNRLP) FIM(NJ+J.KPRV)=1. 245 CONTINUE 240 CONTINUE 250 CONTINUE CALL GAUSS(FIM,X,NJLR) DO 300 I=i.LX KNO=KPRVS(I) DO 300 J=l.NPRV IF (KNONVALV(.,J 300,270,300 270 HGL(J)=X(KNO) X(KNO)=O. JAC(I)=J 300 CONTINUE GTOT=O. GFLCH=O. DO 310 1=1, NP G ( I, KB ) =X ( I ) GTOT=GTOT+ABS(X(I)} 310 GFLCH=ABS(Q(I,KB)G(I,KBl+GFLCH ERR=GFLCH/GTOT0.0005 IF (ERR) 320.320,350 320 DO 327 I=l,LX PRINT 325 325 FORMAT (2X, 'HYDRAULIC GRADE IMMEDIATELY') PRINT 330, JACCI), HGL(I) 330 FORMAT (2X, 'DOWNSTREAM OF PRY', 13,' = ',F6.2/) 327 CONTINUE GO TO 900 350 IF (KBMAX) 360,360.400 360 IF (KB2) 500,500,361 361 DO 380 I::::I,NP DO 365 J=L LX IF (IKPRVS(J 365,370,365 365 CONTINUE Q(I,KB)=(G(I,KB)+GCI,KBl/2. GO TO 380 370 G PAGE 62 200 PRINT 420 201 420 FORMAT (2X. 'IN NO. OF ITERATIONS SPECIFIED') 202 PRINT 430. MAX 203 430 FORMAT (2X, 'NO. OF ITERATIONS SPECIFIED = ',12) 204 PRINT 440, ERR 205 440 FORMAT (2X, 'ERROR ',Fa.5/) 206 GO TO 320 207 900 PRINT 890 208 890 FORMAT ('1', 1.5X, Ip I P E D I 8 C H A R GEl) 209 PRINT 892 210 892 FORMAT (16X, '**************************'/) 211 DO 910 1=1. NP 212 PRINT 920, L G( I. KB) 213 920 FORMAT (16X, 'DISCHARGE IN PIPE NO. ',13. I == I, F5. 2, I CFS 'I) 214 910 CONTINUE 215 DO 2000 I=l,NP 216 DO 1500 J=1,NPRV 217 IF (INVALV(J 1500.1400,1500 218 1400 XL(I)=XL(I)+XLPRV(J) 219 GO TO 1450 220 1500 CONTINUE 221 1450 XKT(I)=(8. 52E5*XL(I}*(ABS(Q(I,KB)**1.852)/ * PAGE 63 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 3612 NB=NB+l GO TO 3600 3610 IF (Jl) 3600,3600,3615 3615 JP1=JNL(I,Jl) IF (HDeJP1 3600,4605,3600 4605 DO 4700 K=l, NPUMP IF (KIPKP(K 4700,4610,4700 4610 HP1=A(K)*(ABS(G(KIP.KB)**2.+B(K)* *ABS(G(KIP,KB+HO(K) GO TO 3625 4700 CONTINUE HP1=0. 3625 IF (LP(I,J 3630,3635.3635 3630 HD(JP1)=HD(JP)XKT(KIP)+HPl GO TO 3637 3635 HD(JP1)=HD(JP)+XKT(KIP)+HPI 3637 NB=NBl 3600 CONTINUE IF (NBNT+l) 3501,3505,3505 3501 IF (NB) 3500,3500,3550 3505 KZ=l 3500 CONTINUE IF (KZ) 3901,3901,3003 3901 DO 4000 I=l,NPRV NPV=NVALV(I) KPV=NBJ(I) IF (G(NPV,KB 4010,4010,4015 4010 HGLU(I)=HD(KPV) GO TO 4000 4015 HGLUCI)=HDCKPV)XKTCNPV)*XLPRVCI)/(XL(NPV)XLPRV(I') HGL(I)=VALST(I) 4000 CONTINUE PRINT 4350 .57 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 4350 FORMAT (1IIi6X, 'H Y D R A U L I C G R A DES 0 F NOD E 8 / ) PRINT 4360 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 4360 FORMAT (16X. '**********************************************'/) DO 4500 J=l,NJ PRINT 4400, J, HD(J) 4400 FORMAT (16X, 'HYDRAULIC GRADE OF NODE NO. I I 12. ::: ',F8. 2. *' FT'/) 4500 CONTINUE PRINT 4510 451 0 FORMAT (1/116 X, 'U PST REA MID 0 W N S T REA M') PRINT 4515 4515 FORMAT (16X. 'G R A D E 8 0 F P R V 8') PRINT 4520 4520 FORMAT (16X, '************************************') DO 4450 I=l,NPRV PRINT 4455 4455 FORMAT (/i6X, 'HYDRAULIC GRADE IMMEDIATELY') PRINT 4456, I, HGLU(I) 4456 FORMAT <16X, 'UPSTREAM OF PRV NO. ',14, 1== ',F8.2. *; FT 'I) PRINT 4460 4460 FORMAT (16X. 'HYDRAULIC GRADE IMMEDIATELY') PRINT 44,61. I, HGL( I) 4461 FORMAT (16X, 'DOWNSTREAM OF PRV NO. 12.' = ',F8.2, *' FT'/) 4450 CONTINUE KB=KBl PRINT 930,KB 930 FORMAT (f 1I16X, 'NO. OF ITERATIONS:::; ',12/) PRINT 990, ERR 990 FORMAT C16X, 'RELATIVE ERROR= ',FI0.6} PRINT 994 994 FORMAT ('I') STOP END 343 SUBROUTINE GAUSS(A,X,N) 344 DIMENSION A(llO, Ill), X(100), Y(IOO) 345 M=N+l 346 N2=Nl 347 DO 800 11=1. N2 348 111=11+1 349 DO 20 I=II.N 350 20 Y(I)=ABS(A(I, II 351 KK=O 352 TR=Y(II) 353 DO 11 I=III.N 354 IF (Y(I)TR) 11.11.12 355 12 TR=Y(I) 356 KK=I 357 11 CONTINUE PAGE 64 358 IF (KK) 13, 14. 13 359 13 DO 15 I=II,M 360 STORE=A(KK, I) 361 A (KK. I ) =A ( I r, I ) 362 15 ACII, I)=STORE 363 14 K=II+l 364 DO 800 I=K,N 365 DO 800 J=K,M 366 800 ACI.J)=ACI,J)A(I.Kl)*ACKi,J)/A(Kl,Kl) 367 X(N)=A(N.M)/A(N,N) 368 DO 86 K=2,N 369 J=MK 370 L=J+l 371 X(J)=O. 372 DO 87 I=L,N 373 87 X(J)='X(J)+A(J. I>*XCI) 374 86 X(J)=(A(J,M)X(J)/A(J,J) 375 RETURN 376 END $ENTRY PAGE 65 P I P E D I S C H A R G E ************************** DISCHARGE IN PIPE NO. 1 :::: 2. 53 CFS DISCHARGE IN PIPE NO. 2 0. 38 CFS DISCHARGE IN PIPE NO. 3 ::; 2.47 CFS DISCHARGE IN PIPE NO. 4 = O. 72 CFS DISCHARGE IN PIPE NO. S :: 0.92 CFS DISCHARGE IN PIPE NO. 6 == L 08 CFS DISCHARGE IN PIPE NO. 7 ::::: L 81 CFS DISCHARGE IN PIPE NO. 8 :::::: 3.19 CFS H E A D L o SSE S ******************** HEAD LOSS IN PIPE NO. t .. 128. 15 FT HEAD LOSS IN PIPE NO. 2 ::::: 2. 64 FT HEAD LOSS IN PIPE NO. 3 :::: 122.03 FT HEAD LOSS IN PIPE NO. 4 == 8. 50 FT HEAD LOSS IN PIPE NO. 5 19.90 FT HEAD LOSS IN PIPE NO. b == 22.65 FT HEAD LOSS IN PIPE NO. 7 :::; 6. FT HEAD LOSS IN PIPE NO. e ::;:; 17.73 FT H Y D R A U L. I C G R A DES a F NOD E S ********************************************** HYDRAULIC GRADE OF NODE NO. 1 = 173.77 HYDRAULIC GRADE OF NODE NO. 2 :::::: 57. 58 HYDRAULIC GRADE OF NODE NO. 3 ::::: 60.22 HYDRAULIC GRADE OF NODE NO. 4 :::: 182.27 HYDRAULIC GRADE OF NODE NO. :5 _. 37. 56 UPS T REA MID 0 W N S T REA M G R A DES 0 F P R V S ************************************ HYDRAULIC GRADE IMMEDIATELY UPSTREAM OF PRV NO.1::::: HYDRAULIC GRADE IMMEDIATELY DOWNSTREAM OF PRV NO.1:::: NO. OF ITERATIONS::::: 6 RELATIVE ERROR= 0.000268 50. 12 FT 50.00 FT FT FT FT FT FT 