Pipe network analysis

MISSING IMAGE

Material Information

Title:
Pipe network analysis
Series Title:
Florida Water Resources Research Center Publication Number 77
Physical Description:
Book
Creator:
Lee, Mun-Fong
Publisher:
University of Florida
Place of Publication:
Gainesville, Fla.
Publication Date:

Notes

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 performance 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 dealing 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 optimization techniques to solve the pipe network problem as a viable and perhaps more versatile alternative to the widely used iterative methods.

Record Information

Source Institution:
University of Florida Institutional Repository
Holding Location:
University of Florida
Rights Management:
All rights reserved by the source institution and holding location.
System ID:
AA00001575:00001


This item is only available as the following downloads:


Full Text

















Publication No. 77


PIPE NETWORK ANALYSIS


By


Mun-Fong Lee


University of Florida
Gainesville


'''
''*
I:i


" .


,i ', ,


-~%














PIPE ;:','a'ORK ANALYSIS


BY





MUN-FONG LEE





NON-THESIS 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 ................

SNEWTON-RAPHSON 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 add-on 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 non-linearity 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 on-line 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) Newton-Raphson (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 non-linear 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 man-made 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 non-linear 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_.d-ent. 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 sub-matrix

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 Darcy-Weisbach, Hazen-Williams and

Manning equations. Perhaps the most widely used of these equations

is the Hazen-Williams 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 Hazen-Williams roughness coefficient, S is the

slope of the energy line and equals hLP/L, R is the hydraulic

radius defined as the cross-sectional 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 Hazen-Williams 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 a-l.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 )
S-out 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 non-linear

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 non-linear 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 counter-clockwise

direction, but the sign convention must be consistent around

any particular loop.


-;- Si -'r *








CHAP "h 3

NEWTON-RAPHSON 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 Newton-Raphson 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 Newton-Raphson

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 Newton-Raphson 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 Newton-Raphson method when applied to this set of

equations becomes

p(m+-) p(m) D-1 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 Newton-Raphson 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 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 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


Non-linearity 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 non-linear. 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 non-linear equations in each of its entire system of

equations, Chapter 3 has dealt with the straight-forward applica-

tion of the Newton-Raphson Method to linearizing the non-linear

equations associated with the latter two methods. This chapter

will be devoted to other linearization techniques, some of which

are variations of the Newton-Raphson 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)(x-xo) (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 (f-1 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'

non-linear.
*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 non-linear 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

n-1
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

(.+f-l) 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 non-linear 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 flow-rates and this procedure reduces

this to a negligible quantity. Flow adjustment is carried out for

all i fundamental (closed) loops and (f-l) 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 as-contri-
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)(Q-Qi)

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)
K--Iab


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
1-n
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 1-n N Q1-n
-_Qi- H- Ha i-
b=1 nKab b;- nKab
a (4.23)

N NV H NF 1l-n
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[((H-Hd)/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 non-linear terms,, there is

always the possibility of convergence problems. Ill-conditioned

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 straight-forward 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 Q-equations). Like the other linear

method described in Section 4.1.1.3, it has several distinct

advantages over the Newton-Raphson 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 non-linear loop equations into

linear equations by approximating the head in each pipe by

n-1
h = (KQ ) Q = K'Q (4.28)
n-1
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 Newton-Raphson 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(m-l) + Qi(m-2) / 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

f-l 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 non-linear, thus

necessitating iterative techniques such as (i) Hardy Cross, (ii)

Newton-Raphson, 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 non-linear 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 one-way 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

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 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 Co-Content 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 Co-Content, J, hence the name Co-

Content Model.

Using Kuhn-Tucker 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, n-rous non-linear algorithms such as (i) Frank-Wolfe method,









(ii) piece-wise linear approximation and (iii) convex-simplex

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 head-discharge 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 Hazen-Williams 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 Newton-Raphson 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 non-linearity,

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, 415-428, 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, Holden-Day, Inc.,
San Francisco, 1-22, 1968.

(4) Bellemy C.J., The analysis of networks of pipes and pumps,

Journal Institution of Engineers, Australia, 37(4-5),

111-116, 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), 467-480, 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), 55-58, 1974.

(9) Cherry E.C. and W. Millar, Some general theorems for non-linear

systems possessing resistance, Phil. Mag., (Ser 7), 42(333)

1150-1177, 1951.

(10) Clay R,, Non-linear networks and systems, Wiley Interscience,

New York, 1-39, 1971.









(11) Collins A.G., and R.L. Johnson, Finite-element method for

water distribution networks, Journal American Water Works

Association, 67(7), 385-389, 1975.

(12) Collins N.A., L. Cooper, V. Helgason, and J.L. Kennington,

Solution of large-scale pipe networks by improved math-

ematical approaches, IEOR 77016-WR 77001, 1978.

(13) Dillingham J.H., Computer analysis of water distribution systems,

part 1, Water and Sewage Works, 114(1), 1-3; part 2, 114(2),

43-45; part 4, Program application, 114(4), 141-142, 1967.

(14) Donachie RP., Digital program for water network analysis,

Journal Hydraulics Division, Proc. Amer. Soc. Civil Engineers,

100 (HY3), 393-403, 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), 189-196, 1976.

(16) Epp R., and A.G. Fowler, Efficient code for steady-state flows

in networks, J. Hydraulics Div., Proc. Amer. Soc. Civil

Engineers, 96 (HYI), 43-56, 1970.

(17) Fietz T.R., Discussion of "Hydraulic network analysis using
linear theory", J. Hydraulics Div., Proc. Amer. Soc. Civil

Engineers, 99 (HT5), 855-857, 1973.

(18) Fietz T.R., Discussion of "Efficient algorithm for distribution

networks", J. Hydraulics Div., Proc. Amer. Soc. Civil

Engineers, 99 (HY1l), 2136-2138, 1973.

(19) Gerlt J.L., and G.F. Haddix, Distribution system operation

analysis model, J. Amer. Water Works Association, 67(7)

381-384, 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), 1-18, 1958.

(21) Guillemin E.A., Tntroductory circuit theory, John Wiley & Sons,

Inc., New York,.5-10, 1953.

(22) Hoag L.N., and G. Weinberg, Pipeline network analyses by elec-

tronic digital computer, J. American Water Works Association,

49(5), 517-524, 1957.

(23) Hudson W.D., Computerizing pipeline design, J. Transportation

Engineering Div., Proc. Amer. Soc. Civil Engineers, 99(TE1)

73-82, 1973.

(24) Hudson W.D., A modern metroplex looks ahead, J. Transportation

Engineering Div., Proc. Amer. Soc. Civil Engineers, 100(TE4)

801-814, 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., 84-113, 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., 481-492, 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), 335-344, 1972.

(29) Lam C.F., and M.L. Wolla, Computer analyses-of water distribution

systems, part 11 numerical solution, J. Hydraulics Div.,

Proc. Amer. Soc. Civil Engineers, 98 (HY3), 447-460, 1972.









(30) Lemieux P.F., Efficient algorithm for distribution networks,
J.iHydraulics Div., Proc. Amer. Soc. Civil Engineers,

98 (HY11), 1911-1919, 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, 36-43, 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),

43-61, 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), 115-129, 1963.

(34) McIlroy M.S., Pipeline network flow analyses, J. Amer. Water

Works Association, 41, 422-428, 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), 507-509, 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. C-4164, 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), 97-108, 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), 281-294, 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),219-234, 1968.

(41) Shamir W., Optimal design and operation of water distribution

systems, Water Resources Res., 10(1), 27-36, 1974.

(42) Tewarson R.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), 1057-1067, 1973.

(44) Wood D.J., and C. Charles, Hydraulic network analysis using

linear theory, J. Hydraulics Div., Proc. Amer. Soc. Civil

Engineers, 98 (HY7), 1157-1170, 1972,

(45) Wood DJ., An explicit friction factor relationship, Civil
Engineering, 36(12), 60-61, 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), 1-14, 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 <.
4-i- 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


ii-scr^


- *'


;,
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!~
r-I -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 -r-sv n io e


"%P"L`I


i; g-- ~t:
i~ ~ C




Full Text

PAGE 1

Publication No. 77 PIPE NETWORK ANALYSIS By Mun-Fong Lee University of Florida Gainesville

PAGE 2

PIPE NETWORK ANALYSIS BY MUN-FONG LEE NON-THESIS 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 NEWTON-RAPHSON 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 ........................................ q-5 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 add-on 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 non-linearity 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 .... d-evelopo-pera 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 on-line 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) Newton-Raphson (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 non-linear 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 man-made 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 non-linear 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 sub-matrix 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 BT-p 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, 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 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 Darcy-Weisbach, Hazen-Williams 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 Hazen-Williams roughness coeffic,ient, S is the slope of the energy line and equals, hLP/L, R is the hydraulic radius defined as the cross-sectional 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 Hazen-Williams 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 turbu-11 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 non-linear 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 non-linear 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 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: '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 arbitra-rily taken positive in the clockwise or counter-clockwise direction, but the sign convention must be consistent around any particular loop. * 17

PAGE 24

CHAPTER 3 NEWTON-RAPHSON 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 Newton-Raphson 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 Jacob--1 ;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 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, 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 Newtqn-Raphson 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 Newton-Raphson 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 Non-linearity 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 non-linear. 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 non-linear equations in each of its entire system of equations. Chapter 3 has dealt with the straight-forward app1ica tion,of the Newton-Raphson Method to linearizing the non-linear equations associated with the latter two methods. This chapter will be devoted to other linearization techniques, some of which are',variations of the Newton-Raphson 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 fund-amental 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 (f-1 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 netv-IOrk. 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 non-linear. *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 non-linear 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, n-l 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 "Balanc-ing Head Method 10 which was limited to .closed loopsysterns and included only line losses. The procedure is generalized and summar-ized 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 non-linear 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 flow-rates and this procedure reduces this to a negligible quantity. Flow adjustment is carried out for all t fundamental (closed) loops and (f-l) 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 solu-tion 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 re-presents 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 I-n 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 I-n 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 Hc-Hd)/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 non-linear terms., there is always the possibility of convergence problems. Ill-conditloned 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 straight-forward 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 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 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 Q-equations). Like the other linear method described in Section 4.1.1.3. it has several distinct advantages over the Newton-Raphson 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 or the corrective loop oriented equations is not recommended. Linear theory transforms the non-linear loop equations into linear equations by approximating the head in each pipe by n-l h ::: (KQ. ) Q :: K' Q L l (4.28) n-1 in which Qi is an estimate of the flowrate, and"K' :; KQ. l Com1)ining these linearized loop equations wi th 'the j-l 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 oscil-lation can be understood by observing that the linear theory method is a variation of the Newton-Raphson 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(m-l) + Qi(m-2) ] / 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 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 f-l 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 approxi-35 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 trans-formation 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 in-lJ 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 non-linear, thus necessitating iterative techniques as (i) Hardy Cross, (ii) Newton-Raphson, 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 non-linear 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 en-sure 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.one-way 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 Co-Content 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 Co-Content, J, hence the name Co-Content Model. Using Kuhn-Tucker 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 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 non-linear algorithms such as (i) Frank-Wolfe method. 40

PAGE 47

(ii) piece-wise linear approximation and (iii) convex-simplex 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 head-discharge 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 pro-gram 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 Hazen-Williams equation and pumps were described by the quadratic form (Equation 2.1]). The convergence criterion employed was: Z I Qi -Qi-ll t' t 'Qi r < 0.0005 in which Qi is theflowrate obtained for a trial and Qi-l 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 1------1 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) Hazen-Williams 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 Newton-Raphson 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 non-linearity, 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,' 415-428, 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, Holden-Day, Inc., 'San Francisco, 1-22,1968. (4) Bellamy C.J., The analysis of networks of pipes and pumps, Journal Insti tutiori. of Engineers, Australia, 37(4 ... 5), 111-116, 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), 467-480, 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), ,55-58, 1974. (9) Cherry E.C. and W. Millar, Some general theorems for non-linear systems possessing resistance, Phil. Mag., (Ser 7), 42(333) 1951. (10) Clay R., Non-linear networks and systems, Wiley -Interscience, New York, 1-39, 1971.

PAGE 55

,(II) Collins A.G.t and R.L. Johnson, Finite-element method for water distribution networks, Journal American,Water Works Association, 67(7), 385-389, 1975. (12) Collins N.A., L. Cooper, V. Helgason, and J.L. Kennington, Solution of large-scale pipe 'networks by improved mathematical approaches, lEOR 77016-o/R 77001, 1978. 49 (13) Di1lingham,J.H., Computer analysis of water distribution systems, part 1, Water and Sewage Works, 114(1), 1-3; part 2, 114(2), 43-45; part 4, Program application, 114(4), 141-142, 1967. (14) Donachie R.P., Digital program for water network analysis, Journal Hydtau1ics Division, Amer. Soc. Civil Engineers, 100 (HY3), 393-403, 1973. (15) Eggener c. L., and L. B. Po1kowski, Netviork models and the impact of modeling assumptions, Journal American Water Works Association. 68(4), 189-196, 1976. (16) Epp R., and A.G. Fowler, Efficient code for steady-state flows in networks, J. Hydra'ulics Di v., Proe. Amer. Soc Civil Engineers, 96 (HYl), 43-56, 1970. (17) Fietz T.R., Discussion of "Hydraulic network analysis using ,linear theory", Hydraulics Div., Proc. Amer. Soc. Civil Engineers, 99 (HT5), 855-857, 1973. (18) Fietz T.R., Discussion of "Efficient algorlthmfor distribution networks", J. Hydraulics Div., Proc. Amer. Soc. Civil Engineers, 99 (HYIl), 2136-2138, 1973. (19) Gerlt J.1J., and G.F. Haddix, Distribution system operation analysis model, J. Amer. Water Works Association, 67(7) 381-384, 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), 1-18, r958. (21) Guillemin E.A. Introductory circuit theory, John Wiley & Sons, Inc.r New York, .5-10, 1953. 50 (22) Hoag L.N., and G. Weinberg, Pipeline network analyses by electronic digital computer, J. American Water Works Association, h9(5), 517-524, 1957. (23) Hudson VI.D., Computerizing pipeline design, J. Transportation Engineering Div., Proc. Amer. Soc. Civil Engineers, 99(TEl) 73-82, 1973. (24) Hudson W.D., A modern metroplex looks ahead, J. Transportation Engineering Div., Proc. Amer. Soc. Civil Engineers, 100(TE4) 801-814, 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., 84-113, 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, 481-492, 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), 335-344, 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), 447-460, 1972.

PAGE 57

(30) Lemieux P. F., Efficient algorithm for distribution networks, J IHydraulics Div., Proc. Amer. Soc. Civil Engineers, 98 (HYll), 1911-1919,-1972. (Jl)Liu K.T., The numerioal analyses of wat-er supply networks by digital oomputer, Thirteenth Congress of the International Association for Hydraulic Research, 1, 36-43, 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), 'n-6l, 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), 115-129, 1963. (34) McIlroy M.S.# Pipeline network flow analyses, J. Amer. Water Works Association, 41, 422-428, 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), 507-509, 1974. (36) Rao H.S., D.W. Bree, Jr., and R. Extended period simulation of water distribution networks; Final technical report OWRR project no. C--4l64, 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), 97-108, 1977.

PAGE 58

(38) Ra6 H.S., L.C. Markel and D.W. Bree, Jr., Extended period simu-lation of water systems B, J. Hydraulics Div., Proc. Amer. Soc. Civil Engineers, 103 (HY3), 281-294, 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),219-234, 1968. (41) Shamir W., Optimal design and operation of water distribution systems, Via tel' ResourcesRes., 10 (1). 27-36, 1974. (1.1--2) ':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) 1057-1067, 1973. (4lt) Wood D.J., and C. Charles, Hydraulic network analysis using linear theory, Hydraulics Div., Proc. Amer. Soc. Civil Engineers, 98 (HY7), 1157-1170, 1972. (45) Wood DolT. ,An explicit friction factor relationship, Civil Engineering, 36(12), 60-61, 1966. (46) Wood D.J., Algorithms for pipe network analysis and their.reli-ability, 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), 1-14, 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 1--3); NJ=NO. OF NODES (MAX=60, COL 4-6);NL=TOTAL NO. OF LOOPS (MAX=40,COL 7-9); NPUMP=NO. OF PUMPS (MAX=10. COL 10-1.2) .. NRLP=NO .. OF REAL LOOPS(LOOPS CON TAINING NO RESERVOIRS/PRESSURE REGULATING VALVES(PRV'S},COL 13-15); NPRV=NO. OF PRV'S(MAX=10, COL 16-18)iMAX=NQ. OF ITERA TIONS (DEFAULT=8, MAX=12 .. COL 19-21) 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)=HAZEN-WILLIAMS 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 1-3); A(I)(COL 4-11>.B(I)(COL 12-19),HO(I)(COL 20-27}=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 (I-NRLP) 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(I-NRLP). Hi< I-NRLP) 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 I-JITH 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 (I-NRLP) 44.44.42 42 FIM(NJ+I,NAPC)=H2(I-NRLP)-Hl(I-NRLP) 44 LN=NLP(I) DO 62 J::::l/LN KPA=IABSCLPCI,J DO 61 K=l.NPUMP IF (KPA-KP(K}) 61146,61 46 IF (LP(r,J 48,50,50 48 HPS=B(K)**2. /(4. *A(K-HO(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,KB-l)**0.852) */(HWC(I)**l. 852*DIA(I)**4. 87)+8. *XKM(I>*ABS( *G( I, KB-1/ (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 (KPA-KP(K 100.80,100 80 IF (LP(I.J 82,84,84 82 FIM(NJ+I,NP+K)=A(K)*ABS(G(K,KB-l GO TO 86 84 FIM(NJ+I, NP+K)=-A(K)*ABS(GCK,KB-l 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,KB-l+GFLCH ERR=GFLCH/GTOT-O. 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 (KB-MAX) 140,140,150 140 IF (KB-2) 500,500,141 141 DO 142 I=l,NP 142 GCI,KB)=(G(I,KB)+GCI,KB-l/2. DO 145 I=l,NPUMP 145 GCI.KB)=(G(I,KB)+G(I,KB-l/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 (KPRV-IABS(LP{J,JJ) 245,213,245 213 IF 215,217,217 215 FIMCNJ+J.NAPC)=FIM(NJ+J,NAPC)-H2(J-NRLP) FIM(NJ+J,KPRV)=-l. GO TO 245 217 FIM(NJ+J, NAPC)=FIM(NJ+J, NAPC)+Hl(J-NRLP) 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 (KNO-NVALV(.,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,KB-l+GFLCH ERR=GFLCH/GTOT-0.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 (KB-MAX) 360,360.400 360 IF (KB-2) 500,500,361 361 DO 380 I::::I,NP DO 365 J=L LX IF (I-KPRVS(J 365,370,365 365 CONTINUE Q(I,KB)=(G(I,KB)+GCI,KB-l/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 (I-NVALV(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)/ ***1. 852*DIA(I)**4. 87)+8. *XKM(I>*ABSCQ(I,KB/ *(32.2*3.141593**2. *DIA(I>**4.) 222 2000 CONTINUE 223 PRINT 4100 224 4100 FORMAT (/1116X, 'H E A D LOS S E 8') 225 PRINT 4102 226 4102 FORMAT (16X, '********************'/) 227 DO 4300 I=1,NP 228 PRINT 4320, I, XKT(I) 229 4320 FORMAT (16X. 'HEAD LOSS IN PIPE NO. ',13. I = ',F8.2 *," FT 'I) 230 4300 CONTINUE 231 DO 2500 I=I.NL 232 MZ=NLP(I) 233 DO 2400 J=1.MZ 234 MZP=IABS(LP(I,J 235 IF (Q(MZP,KB 2410,2400,2400 236 2410 LP(I,J)=-LP(I,J) 237 2400 CONTINUE 238 2500 CONTINUE 239 DO 3002 J=I,NJ 240 3002 HD(J)=O. 241 NRLP1=NRLP+1 242 DO 3000 I=NRLPl,NL 243 NT=NIK(I) 244 DO 3100 J=1.NT 245 JP=JNL(I,J) 246 KIP=IABS(LP(I.J 247 DO 3200 K=I,NPUMP 248 IF (KIP-KP(K 3200.3150.3200 249 3150 HP1=A(K)*(ABSCGCKIP,KB)})**2.+B(K)* *ABS(Q(KIP,KB+HOCK) 250 GO TO 3210 251 3200 CONTINUE 252 HP1=0. 253 3210 IF (J-1) 3212,3212,3240 254 3212 DO 3218 L=1,NPRV 255 IF (KIP-NVALV(L 3218,3214.3218 256 3214 XKT(KIP)=XKT(KIP)*(XL(KIP)-XLPRV(L/XL(KIP) 257 GO TO 3220 258 3218 CONTINUE 259 3220 IF (LP(I,J 3225,3230.3230 260 3225 HDeJP)=H2(I-NRLP)+XKT(KIP)+HP1 261 GO TO 3100 262 3230 HD(JP)=H2(I-NRLP)-XKT(KIP)+HPl 263 GO TO 3100 264 3240 JP1=JNL(I,J-l) 265 IF (LP(I,J 3245.3247,3247 266 3245 HDeJP)=HD(JPl)+XKT(KIP)+HPl 267 GO TO 3100 268 3247 HD(JP)=HD(JPl)-XKT(KIP)+HP1 269 3100 CONTINUE 270 3000 CONTINUE 271 3003 KZ=O 272 DO 3500 I=1,NRLP 273 NT=NIKCI) 274 3550 NB=O 275 DO 3600 ..1= 1. NT 276 JP=JNL (I, \oJ) 277 KIP=IABS(LP(I,J 278 IF (HD(JP 3610,3612.3610 56

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 (J-l) 3600,3600,3615 3615 JP1=JNL(I,J-l) IF (HDeJP1 3600,4605,3600 4605 DO 4700 K=l, NPUMP IF (KIP-KP(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=NB-l 3600 CONTINUE IF (NB-NT+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=KB-l 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=N-l 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.K-l)*ACK-i,J)/A(K-l,K-l) 367 X(N)=A(N.M)/A(N,N) 368 DO 86 K=2,N 369 J=M-K 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