Salgado R. (1988) PH D Thesis Vol I
Salgado R. (1988) PH D Thesis Vol I
Salgado R. (1988) PH D Thesis Vol I
by
Ruben 0. Salgado-Castro
VOLUME ONE
**********
COB 22970 8
-4..12S-i s L.2k
November, 1988.
ACKNOWLEDGEMENTS
this thesis and contributed with many discussions and ideas; his
time.
created the gradient method and contributed with his advice for
the facilities for carrying out this work and I am grateful for
that too. Mr. R. Mackay, lecturer in Hydrology, kindly allowed me
true values. For the estimation of the unmeasured flows, the raw
roughnesses.
TABLE OF CONTENTS
VOLUME ONE
**********
CHAPTER ONE
INTRODUCTION 1
CHAPTER TWO
2.1. Introduction 18
2.2. Hardy Cross (1936) methods 20
2.2.1. Method of balancing heads 22
2.2.2. Method of balancing flows 26
2.3. Newton-Raphson-based methods 33
2.3.1. Introduction 33
2.3.2. The mathematical background of the Newton-
Raphson method 34
2.3.3. Martin and Peters formulation (1963) 39
2.3.4. Shamir and Howard nodal formulation (1968) 44
2.3.5. The Newton-Cross method of Liu (1969) 49
2.3.6. Epp and Fowler loop formulation (1970) 51
2.3.7. The mesh-nodal approach of Hamam and Brameller
(1971) 56
2.3.8. Other formulations 65
2.4. The linear theory method 68
2.4.1. Wood and Charles formulation (1972) 68
2.4.2. Nodal formulation of Isaacs and Mills (1980) 72
2.4.3. The gradient-based formulation of Wood (1981) 75
2.5. Other approaches 78
2.5.1. Optimisation-based methods 78
2.5.2. Methods based on unsteady state analysis 79
2.5.3. A method based on a finite element approach 80
2.5.4. The gradient method of Todini (1979) 80
2.6. Comparison of the performance of some of the
existing methods: need for a more reliable
algorithm 88
CHAPTER THREE
GENERALISATION OF THE GRADIENT METHOD
TO INCLUDE PUMPS 99
3.1. Introduction 99
3.2. Traditional formulation of the network analysis
problem 99
3.3. Minimum power dissipation formulation: necessary
conditions 100
3.4. Todini's gradient method 102
3.5. A general head loss/flow model to include pumps
and valves in the gradient method 104
3.6. Extension of the gradient method to include pumps:
derivation of the recursive scheme for the
extended method 107
3.7. On the existence and uniqueness of the network
analysis solution 113
3.8. Other gradient formulations 121
3.9. Comparison of the gradient method with some of
the existing methods of network analysis 130
3.9.1. Introduction 130
3.9.2. The computer programs used in the comparison 130
3.9.3. The set of test examples 132
3.9.4. Comparison of the algorithms for the
simultaneous path adjustment, linear theory
and gradient methods 137
3.10. Concluding remarks 141
CHAPTER FOUR
EXTENSION OF THE GRADIENT METHOD TO INCLUDE
REGULATING VALVES:
CHAPTER FIVE
CHAPTER SEVEN
CHAPTER EIGHT
REFERENCES
REFERENCES 375
VOLUME TWO
**********
APPENDIX A
vi
A.4.7. Standard conjugate gradient method for the
solution of linear systems of equations.
(Hestenes and Stiefel, 1952) A-127
A.4.8. Preconditioning A-133
A.4.9. Preconditioned (modified) conjugate gradient
method for the solution of linear systems A-142
APPENDIX B
APPENDIX C
APPENDIX D
FIGURE PAGE
viii
FIGURE PAGE
FIGURE PAGE
TABLE PAGE
Table 2.1. Results of the comparison of 60 small water
distribution networks (less than 100 pipes)
from Wood (1981a) 89
Table 2.2. Comparison of network analysis methods in
terms of global and relative CPU times, and
number of iterations, from Carpentier et
al. (1985) 92
Table 3.1. Summary of the performance of the programs
LOOP, LT and GRAD with test examples 138
Table 4.1. Summary of effects of PRV's in different
formulations of network analysis equations
[following proposition of Jeppson (1976)
and Jeppson and Davis (1976)] 165
Table 4.2. Comparison of results reported in the
literature and our results for networks
containing pressure control devices 180
Table 4.3. Comparison of results presented by Jeppson
and Davis (1976) and our results, for
example JEPPO 182
Table 4.4. Comparison of the results given by Jeppson
(1976) and our results for example JEPP1 184
Table 4 .5. Comparison of the results given by Jeppson
(1976) and our results, for example JEPP12. 186
Table 4.6. Comparison of the results given by Jeppson
(1976) and our results, for example JEPP13 188
Table 4 .7. Network solutions for example COLLINS, from
Collins (1980) 190
Table 4.8. Operational cases for example COLLINS,
based on all the possible combinations of
pump operating modes 191
Table 4.9. Results of running our program for the
operational modes possible in example
COLLINS 192
Table 5.1. Main data corresponding to examples used
in the comparison between different linear
solvers in the gradient method for pipe
network analysis 207
Table 5.2. Comparison of execution times and number of
iterations required by different linear
solvers in the gradient method for pipe
network analysis 208
Table 5.3. Comparison of storage requirements (Bytes)
for different linear solvers used with the
gradient method for pipe network analysis 210
Table 6.1. Network data for Examples A, B and C 279
Table 6.2. Measurement data for Examples A, B and C 285
Table 6.3. Measurement data for Examples D, E and F 285
TABLE PAGE
Table 6.4. Initial values of C (Hazen-Williams) used
for testing the calibration algorithm.For
Examples A and B and low and high
uncertainty 288
Table 6.5. Initial values of C (Hazen-Williams) used
for testing the calibration algorithm in
Example C 289
Table 6.6. Initial values of C (Hazen-Williams) used
for testing the calibration algorithm. For
Examples D, E and F, and for low and high
uncertainty 290
Table 6.7. Modified nodal demands (us), used to study
the impact of bad demand estimation on the
calibration algorithm. Examples A, B and C 292
Table 6.8. Modified nodal demands (us), used to study
the impact of bad demand estimation on the
calibration algorithm. Examples D, E and F. 293
Table 6.9. Comparison between estimated piezometric
heads and its residuals, for example A 303
Table 6.10. Summary of the comparison between different
piezometric head estimation procedures 304
Table 6.11. Estimated piezometric heads for Cases I, II,
III, IV and V, using the one-dimensional
interpolation procedure 305
Table 6.12. Performance indices for the different head
estimation procedures, for examples A, B,
C, D, E and F 308
Table 6.13 Summary of the comparison for the
calibrated piezometric heads: average
heads 313
Table 6.14 Summary of the comparison for the
calibrated piezometric heads: variance of
the heads 314
Table 6.15 Summary of the comparison for the
calibrated piezometric heads: maximum
residual 315
Table 6.16. Summary of the comparison for the
calibrated piezometric heads: variance of
the absolute value of the residuals 316
Table 6.17. Summary of the comparison for the
calibrated flows: average flows 322
Table 6.18. Summary of the comparison for the
calibrated flows: variance of the flows 323
Table 6.19. Summary of the comparison for the
calibrated flows: maximum residual 324
Table 6.20. Summary of the comparison for the
calibrated flows: variance of the absolute
value of the residuals 325
Table 6.21. Summary of the comparison for the
calibrated C's: average C's 332
Table 6.22. Summary of the comparison for the
calibrated C's: variance of the C's 333
Table 6.23. Summary of the comparison for the
calibrated C's: maximum residual 334
TABLE PAGE
VOLUME TWO
**********
TABLE PAGE
xiv
TABLE PAGE
NETWORK ANALYSIS
Symbol Description
valves, etc.).
vector Q.
component of vector H.
approach).
vector.
vector h.
ij = H i - H j I
Ih
component of vector g.
Sa Flow (mass) nodal imbalance, tNN-NS)*1) column vector.
Qij = R ij nijl/n
component of vector a.
source nodes).
A 11 * = A 11 diag.(13/Q)
B21 1
_, F 13 11_L B 12 1
F =
B 22 L A 21 : 0 J L B 21 : B 22 -I n•n••••n
xix
Symbol Description
I
A22 = diag.(-6q-/Op-).
3. 3.
g o = a s + A 22 -a + A 22 Ps
variable.
LINEAR SYSTEMS OF EQUATIONS
Symbol Description
n
<x,y> = xT Y = yT x = E xi yi
i=1
KRIGING
Symbol Description
Simplified notation:
Cov(Z(x1,H),Z(x2,H)) Covariance:
estimated Xo , i.e.:
n I
Yo* = Y0 * (X0 ) = E
%o i Y i I
i=1
SPLINES
Symbol Description
h+4
s(x) = E r i M(x)
i=1
Bi-cubic splines:
h+4 k+4
s(x,y) = E E r ij Mi(x)Nj(Y)
i=1 j=1
INTRODUCTION
on the one hand, and the pressing need for the efficient
can be evaluated with the help of a model, and only the most
efficient alternative can be implemented in the real system.
where increased consumptions are fed into the model to study the
networks, have been used for that purpose. The advent of digital
of a reality, this does not mean that the model itself is simple.
A computer-based model of a water supply distribution system
a) Physical component.
b) Demand component.
c) Operational component.
d) Measurement component.
key factor, since they enable us to determine how close our model
measurements are then the bridge between the model and reality.
Typical field network measurements are: p ressures, water
e) Hydraulic component.
etc.
f) Mathematical component.
small and simple problems, this leads to the need for a digital
computer.
g) Computational component.
example:
r-
I
Measurement
component
L_ _
Calibrated
Model
z
Usable
results/
Fig. 1.1. The water distribution network model and its main
components.
The work presented in this thesis, contributes to some key
aspects of the computational component of water distribution
followin g problems:
unmeasured variables.
algorithms.
10
the main algorithm, an unreliable network analysis routine is not
the initial aim of this work has been to expand the method to
the network analysis problem. A further aim of this work has been
11
that of its simulation and optimisation applications) relies
network model.
algorithm which can use the raw model (see Fig. 1.1.), the
12
1.4. Simulation of water distribution networks.
The gradient method for network analysis has been used in the
Reservoir levels
+ Consumptions
Network analysis Extended period
program simulation
Piezom. heads program
+ Flows
Results
14
definite linear systems of equations generated by the gradient
method.
the same purpose for the case of the estimation using bi-cubic
15
1.6. Summary of main achievements.
the inclusion of pumps and valves into the algorithm, and which
the widest possible range of sparse linear solvers, and the most
future.
16
d) A new algorithm for the automatic calibration of network
models has been proposed and tested with synthetic data. The
2.1. Introduction.
finding the flow in each pipe and the piezometric head in each
much smaller than the friction head losses, and they can be
smaller than the friction losses, and they can also be neglected;
network, and not along the pipes, as happens in the real world;
18
e) all the pipes are working under positive pressure;
the flow distribution does not depend on the ground level. Only
node into a lower one, but the same is not generally true when
19
simultaneous equations is fully exploited; most of these methods
introduced.
20
H. Cross stated that two conditions for the steady state
is zero.
There are two H. Cross methods, each one starting from the
21
within some accuracy limits. We shall review both methods in the
next paragraphs.
also known under such names as: the "head balance" method, loop
This method can only be used for looped networks with a single
following way:
Let the head loss/flow relationship for the pipe joining nodes
where:
Qij represents the flow from node "i" to "j", which initially
22
The second condition establishes that, under steady state
conditions, these flows must produce pipe head losses which are
balanced round any loop, i.e. the algebraic summation of the head
E h 1 J• = 0 V loop (2)
loop
correction SQk can be calculated for each loop "k", so that the
all the pipes in a loop does maintain flow continuity. The new
head loss in the pipe joining nodes "i" and "j" will be:
- E Kij Qij2
8Qk = (7)
2 E Kij Qij
i.e.
- E hii
eQk = (3)
2 E (hij/Qij)
23
programmable calculator.
In general, for a head loss/flow relationship of the form:
- E hij
8Qk = (10)
n E (hij/Qij)
method depended on the way the loops are chosen, suggesting that
the loops should be chosen in such a way that the common pipes
resistance".
6Qk [computed from the full binomial expansion (4)], reduced the
24
(10 nodes and 13 pipes).
- A E hij
8 Qk =
n E (hij/Qii)
nodes).
25
is not critical in the initial iterations. This was important in
hand calculations.
computer; these loops are the natural set of loops, i.e. those
As it will be seen later on, the Cross loop method can also be
method.
H. Cross proposed that at any stage when flow continuity has not
26
their resistances".
satisfactory".
by the user), let us assume that the node (flow) balance is not
8Q i = E Q ij qi (12)
where:
flow through the pipe joining nodes "i" and "j".
27
from the following relationship:
SQi j
h-1J- + SH-1 = K-1 •J Q..
1J [ 1 + --- -I (15)
or,
h•1J
• + SH-1 = 1J (16)
6 j I 1 +n SQij
1+ die
(17)
Hence,
then,
SQij
1J 1 1J n 1J -- (19)
SQij
SHi = n hij ---- (20)
28
which is reordered to give:
We do not know the value of SQij for each pipe, but we do know
their summation (EQij) for any node "i", which is exactly the
1
E SH- = E 8Qij (23)
j n h-lj
• j
E Q ij - qi
i
SH i = (25)
E
j n h•1J-
29
estimates can be as simple as linearly interpolated values
ground level, produces the initial head. Let this initial head be
number.
in use.
4. Calculate the ratio Qij ( k )/n h 1 (k) for each pipe and add up
"i", as:
R i (k+1) = H i (k) + s1j1(k)
(29)
7. The process can be stopped either if the maximum 6'H(k) is
30
imbalance SQ i is less than a maximum limit. Otherwise, increment
the advantage that the initial solution for the heads is much
convergence can be very slow, and this is the main reason why
this method was not used initially, when computers were not
available.
approach.
given.
31
possible solutions. According to them, a faster conver g ence can
be achieved, if the head corrections (29) are not made in the
considered.
32
are carried out one node at a time, without allowing proper
different loops.
does not take full advantage of them, because it does not solve
2.3.1. Introduction.
33
2.3,2. The mathematical_ background of the Newton-Raphson
method.
Let
f(x) = 0 (30)
be a non-linear function, whose solution is sou ght within a
34
In fact, if x( k) is close to the true solution x( 0 ), the term
(8x( k )) 2 and its higher order terms will be very small, and they
textbooks; the user must be aware of them, but they are not
f i (x l , x2 , x3 , . . , xn ) = 0
f 2 (x 1 , x2 . x3 , .. • , xn ) = °
(35)
fn (x l , x2 , x3 , • , xn ) = °
35
which can be represented in vector notation as:
f(x) = Q (36)
ef l ell
- 1 Sxn ( k ) = 0
fi(xl,...,xn)+ --- ex i ( k )+ --- 8x2(k)+...+ !f
exl Ox2 erxn
6f 2 ef2
f2(xl,...,xn)+ --- ex 1 (k).1. ___ ex2(1).4....4. n 6x(k) = 0
Ox i153(2 Oxn
(37)
Of: Of
40t xn)
On. Ox 2 Oxn
ef 2 af 2 6f2
6x2(k)
On. Ox 2 exn
=-1* (38)
Of n Ofn afn
ex ]. Ox2emn
where, again, all the expressions are evaluated at x=x(k)
36
The s y stem (38) is a linear system of equations in the unknown
vector Sx(k) = ( sxl ( k) , 6x 2 (k) , ... ,sxn(k) ) T , the column vector
of corrections.
af l 6f1 Of1
Ox i Ox2 OXn
6f 2 6f2 Of2
j(k) = Ox i tpx2 dxn (39)
:
Of n Ofn 6fn
1 Oxn
- 6)( 15x2
37
As in the one-dimensional Newton method, the iterative
p rocedure is stopped when two successive solutions are close
I 8X(k) I < 6
39
m
1Hi-Hil
E Qij - q i = E (Hi - Hj) - qi = 0 (46)
J j K.13.1/n
m
Illi-Hil
f i (H 1 , H 2, ..., H NN )= E (H i - Hj) - qi = 0 (47)
j 1J
IH1-Hil
f 1 (H 1 , H 2 , ..., HNN )= E (H 1 - Hi) cil = 0
i Kiil/n
m
IH2-HI
f 2 (H 1' H 2 , ..., H NN )= E. (H 2 - Hi) - q2 = 0
• i K 2 , '/fl (48)
IHNN-HjI
fNN (H 1 , H 2' ...' H NN )= E (H NN - Hj) q-NN = 0
i K NNj l/n
Note that for all source nodes (reservoirs for example), the
40
equations, where NS is the number of sources. This condition is
not introduced into the formulae, in order to keep the notation
6f 1 (H) 1 IHi-Hilm
= for iXj (49)
6 Hj n Kij 1/n
and
Ofi(l) Of(H)
= E for i=j (50)
6 Hi j 6 Hj
0f1(fl) Of(H)
= (51)
6 H -J 6 Hi
Of i ( H)1 1 IHi(k)-iii(k)Im
_ for i7j (52)
6 H J- Klj
11=11(k T
41
Ofi(H)] 1 IHi(k)-Hi(k)im
_ E for i=j (53)
6
Hi n j 1J
H=H(k)
42
e.g. the maximum norm of the residuals, the maximum number
of iterations, etc.
when two adjacent nodes have almost the same piezometric head,
the Newton-Rap hson solution of the network, and the real network.
They recognised the fact that, normally, agreement will not be
43
2 3,4. Shamir and Howard nodal formulation (1968).
done at MIT. Also, some work carried out by Pitchai (1966) was
referenced.
One of the main reasons why this paper became classical, and a
computational cost.
fi = EQ ij - q i = 0 (56)
nodal equations, and taking just the first order terms for the
extended set of unknowns:
44
Z. : set of unknown nodal consumptions (subset of g).
the following linear system of equations is obtained:
Ofi Ofi
fi(X,Y,Z) + E --- 8Xj + E --- 8Yj + E 8ij 8Zj = 0 (57)
j OXj j elYj j i=1,... NU
where:
NO : total number of unknowns.
Kronecker delta (8ij=1 if i=j and zero otherwise).
8Xj : correction to the unknown piezometric head in node "j".
Let us assume that there are, as before, NH nodes with
unknown heads.
8Yj : correction to the unknown resistance in branch "j". Let
us assume that there are NR branches with unknown
resistances.
8Z j : correction to the unknown nodal consumption in node "j".
Note that, apart from the third and fourth terms of the left-
hand-side of equation (57), the problem can be dealt with in
exactly the same way as Martin and Peters (1963).
NN = NH + NR + NC (58)
eYN?(k)
6Z 1 k)
af NN MNN MNN 6fNN M NN MNN
---
6X 16XNH 617 1OYNH 6Z 1OZNc
ezNc(k) -fNij
( X1(°),..,XNH(0),Y1(0),-.,YNR(°),Z1(°),--,ZNC(°) )T
Peters' (1963) method, and both the Jacobian coefficients and the
(Hi-Hi)
fi = E qi = 0 (60)
j Kii 0.54 ifii_Hil0.46
46
and
6fi(H) 6f1(H)
= E for i=j (62)
6 Xi j 6 Hj
and
6f1(H) - 0.54 (Hi= Hj)
= ------------------ (63)
Kij1.54 ificHil0.46
a Ilk
where pipe "k" is connecting nodes "i" and "j".
and, finally:
6fi(H)
= - 1 if i=k, and 0 otherwise. (64)
6 Zk
Z i ( k+ 1) = Z 1 (k) + szi(k)
Z 2 ( k+1 ) = Z 2 ( k ) + sz2(k) z(k+1) = z(k) + 6'Z(k)
z Nc (k+1) =
Z NC (k) + SZNC(k)
47
Because of the presence of unknown resistances and nodal
consumptions, the condition that the number of nodes must be
the existence of the solution will depend on the way the unknowns
(1973):
unknown:
should not have more than one unknown in addition to the unknown
resistance.
Shamir and Howard (1968) stressed the point that any device can
48
2 3.5- The Newton-Cross method of Liu (1969)
sHi(k) - fl(H)
OH2 °NH
M2 M2 M2
6H 2 (k) - f2 (H)
6H 1 151-12 aHNH
• (66)
afl Of Of1
0 0 0
OH]. 6E12 aHNH
49
of i (H) mi(fi)
=-E
0 111 j 6 Hi
the diagonal terms are much bigger (in absolute value) than the
app l y ing the same concept behind the method of Jacobi (or
simultaneous dis p lacements method) for the iterative solution of
6f1
0 •. 0 - flap
6f2
0 ... 0 - f2(H)
OH2
* = (68)
:
6f NH
0 0 ••• 6HNH(k) - fNH(11)
6H NH
_
with all the expressions evaluated at 11=11(k)
- fi(H(k))
6H i ( k ) = V i=1,2,...,NH (69)
6f i 1
6H-
[ 1 H=H(k)
50
It is interesting to note that Liu's method is, in essence,
(69) and (25), both are completely equivalent. The point here is
nodes.
the fact that, normally, the number of loops is smaller than the
51
standard (minimum) connectivity data. Also, they presented a loop
starts in node "i" and ends in node "j"; thus a flow or head loss
opposite direction.
E hij = 0 (70)
loop
where the summation is carried out throughout all the pipes i-j
52
belonging to the loop.
done in the Cross' loop method, but now is done for all the NL
written as:
NL
E K ij (S ji Q ij + E S ijh 8Q 1 (k))n = (72)
Vi,jEl h=1
1=1,2,..., NL
where:
0 otherwise
Note that the second summation in (72) accounts not only for
the loop flow correction in the current loop, but also for the
53
corrections in the neighbouring loops (for these pipes belonging
to more than one loop).
NL
f 1 (Q,8Q) = E K ij (8 ij1 Q ij .4. E rs ijil eQh (k))n = 0
Vi,jEl h=1
NL
f2 (Q,6Q) = E Ki j (e ij2 Q ij + E 6 ijh 8Q h (k))n = 0
Vi,jE2 h=1
(73)
NL
f NL (§2,6Q) = E K-V Og ijNL Qij + E S
ijh 8Qh ( k )) n = 0
Vi, jEN h=1
solved successively:
-
of i af l afl
... 6'Q1 (k) - f1(42,8Q)
084:4 1 '58 % 68QNL
6f 2 af 2 6f2
... 8Q2(k) - f2(Q,8Q)
06Q 1 68Q2 • 68QNL
* = (74)
:
: :
af NL tlfNL afNL
• • 8QNL(k) - fNL(Q,6Q)
66Q 1 158Q 2 68QNL
54
the initial solution is balanced.
afi
= + E nK-
ij Q-
ij-n-1 for 1=1,2, ..., NL (75)
6 8Q 1 Vi, jEl
afl
= - E n Kij Qi3n-1 for 1X1( and 1=1,2,... NL (76)
OeQk Vi, jEk
Vi, jEl
If Qij (k+1) and Q1 (k) are close enough, say their difference
is less than E, a pre-specified accuracy, or if I5( k )1<E, then
more than one fixed head source can be performed via the
two nodes. They also gave a procedure for including pumps within
this formulation.
evidence that the loop approach performs better than the nodal
55
problems (less than 200 nodes); the nodal approach was apparently
better for larger problems.
Epp and Fowler (1971) recognised the possibility that the loop
56
The branch-to-node incidence matrix for the source nodes A10,
as loop j.
A 13 (i,j) = -1 if pipe i flows in the opposite direction
as loop j. (80)
0 if pipe i is not in loop j.
A21 = A l2 T(81)
A01 = A lO T(82)
A31 = A13 T(83)
Nodal balance:
A21 Q = a (84)
Pressure drop in a pipe:
h = A l2 H (85)
57
Loop balance:
A31 h = (86)
The head loss flow relationship in pipes, can also be
simplified with matrix notation; thus:
h = A ll Q (88)
where:
A11 = diag(K i lQ i I n-l ) i=1,2,.., NP (89)
with Ki and Qi denoting the resistance and flowrate at pipe
A31 Al 2 H = Q (90)
A31 A 10 = 0
A21 A 13 = ° (92)
A01 A 13 = °
58
introduced and used , e.g. the concepts of node (vertex), branch
the subgraph.
a) It is a connected subgraph.
defines a tree with only (NN-1) edges on it; this can be proved
by induction.
59
to as "branches", and the corresponding edges of the co-tree are
pump, valve, etc. Fig. 2.1. helps to visualise some of the graph-
NP n 16 edges
NN 11 nodes
65 m
edge
o) Network graph.
(NN — 1) m 10 branches
branch
NP — (NN — 1) .- 6 chords
60
The concepts of tree and co-tree are relevant to the mesh-nodal
partition the flow vector Q into two sub-vectors, say Q2 and Q3,
the tree) and the latter to the independent flows (i.e. flows in
the co-tree). Hence:
flow vector (93) into the nodal balance equation (84), we get:
Q2
A21 Q = A21 HI = a (94)
Q3
61
A 21 = [ A 22 : A23 ] (95)
A22 Q2 + A23 Q3 = g
Q = Qa + M 13 Q3 (99)
where:
A22 -I 1
Qa = (100)
03
[
and
[ - A 22 -1 A231
M 13 = (101)
13
62
A22 1
A l2 = [ ---- (102)
A32
A22
A l2 H = [ 7, ---I H = h (103)
A 32
h = -h2- (104)
h3
{
A22 H = h2 (105)
and
A32 H =h3 (106)
H = A22 -1 h2 (107)
which, when introduced into equation (106) gives:
63
[-- A 22 -1 A 23 1
A 13 =
13
Q = Qa + A 13 Q3 (112)
Now, introducing (112) and (88) into the loop equation (86), we
get:
[ A 31 A ll A 13 ] Q 3 + A 31 A ll Q a =0 (113)
[ A22 -1 g ]
Qa = (114)
03
J 33 = [ A 31 (N A ll ) A13 ] (116)
with:
64
or
J33 6' Q 3( k) = -f( Q3(k)) (119)
solved, since the head losses can be computed with (88) and the
have been made through the years, some of them are reviewed in
Lam and Wolla (1972a and b), based on previous work by Broyden
where the the Jacobian matrix (and its inverse) is not explicitly
65
time spent in the assembly and inversion of the matrix. Another
attributed the failure to the fact that the initial values of the
For the cases where the algorithm may present a low convergence
achieved.
66
independent of the numerical method used to solve the non-linear
Six.
system for the nodes (or loops, depending on the approach in use)
67
(1974), Chin et al. (1978), Gay et al. (1978), and others,
concentrated on this subject.
68
In the loop approach, it has been already noted that the
potential continuity around the loops (or energy paths). This can
be expressed as
E h 1„1
. - .-
= E K 1„1 Q ij
. - n = 0 loop=1,2,.., NL (121)
loop loop
The main idea behind the linear theory method of Wood and
69
respectively, it can be assumed that both are close to the
equilibrium flow vector, and that both are very similar. Hence,
the non-linear term in equation (121) can be approximated by:
Q ij n z [Qij(k)r-lQij(k+1) (123)
which is linear in Qij(k+1).
E K ij [Q ij (k)](n-1)Q ij ( k+ 1 ) = 0 (124)
loop loop=1,2,... NL
which can be simplified to:
(128)
E . -(k) Q 1. J.(k+1) = 0
K 1J loop=1,2, ..., NL
loop
70
At the beginning of the iterative process (k=0), any imbalanced
and compute:
Q(k) = 1.i, ( Q(k-1)+Q(k-2)) (129)
It has been argued and recognised [Fietz (1973), and Wood and
71
important with the advent of bigger and cheaper computer memory,
theory method of Wood and Charles (1972), Isaacs and Mills (1980)
p resented the corresponding nodal-based version of the method.
is highly desirable.
72
The rationale behind Isaacs and Mills' linear theory method is
expressed as:
and, if the node continuity equation for node "i" is written as:
conditions:
Q(0), Q(1) , . . , Q(k), Q(k+1) (133)
approximated as:
Q ij (k+1) lQij(k)in-1 = (1/K 1j ) [H 1 (k+1) _ Hj(k+1)] (134)
or,
Q ij (k+1) = 0 1j (k) [H 1 ( k+1) _ Hj(k+1)] (135)
with
c.1 J.(k) = 1/ (Kij 1Qij(k)In-1) (136)
E Q ij (k+1) _ q i = 0
1=1,2,.. NN-NS (137)
73
symmetric, and appropriate storage handling can reduce
For updating the flows at each stage, Isaacs and Mills (1980)
proposed the following method:
a) Compute:
[Hi(k+1) _ Hj(k+1)]
Qij* = (139)
1Qii(k)1n-1
updated flow:
Q ij (k+1) = (Q ij * Q ij (k)) (140)
continuity condition .
74
presence of nearly zero flows will produce unbounded values, and
performance.
original linear theory method (Wood and Charles, 1972), the main
75
initial and final nodes: 8E 1 . A loop is simply a path where
The fact that the loop balance is not exactly met, means that:
E h(Q(k)) = E K ij [Q ii (k)]n 6E1 (142)
1 1 1=1,2,..., NL
where the summation is carried out over all the NL loops (or
"energy paths") in the network.
order terms:
Oh
E h(Q(k+1)) 2., Eh(Q (k) ) 4- E ---- (Q(k+1) _ Q (k) ) = 6E, (147)
1 1 1 OQ
Q=Q(k)
1=1,2,... NL
76
ordering (147) leads to:
Oh Oh
E [____1 Q (k+1) = E 1[____I Q(k) _ h(Q(k))1 + 8E1 (148)
1 OQ 1 OIQ
Q=Q(k) Q=Q(k)
1=1,2,... NL
which represents a linear system of NL equations in NP unknowns.
E . Q ii (k+1) qi = 0
i=1,2,..., NN-NS
[ Oh Oh
Q(k+1) = E 1[____I Q(k) -
h(Q(k))1 + 6E1 (150)
1 (5Q 1 6,Q4
Q=Q(k) Q=Q(k)
1=1,2,... NL
for averaging the new flows with those of the previous stage.
Wood (1981a) pointed out that the new gradient-based method has
77
By this time (1981), the linear theory method had a reputation
view.
research carried out in this field during the last decade, but
they have not been widely used. Only for the sake of
78
fundamentals of the network analysis problem. From the practical
We believe that this is not the best way to approach the steady
79
technique. should be used only in situations where slow or fast
(2.3.7), as:
80
subject to A 21Q = a (152)
where:
allQ1In-1
*21Q2In-1
A 11 = (153)
aNPIQNpin-1
(1979) found that the necessary conditions for the steady state
form:
[ A li ' A l2 1[Q 1 [Q I
(154)
A21: ° H a
81
system of NN-1 independent linear equations in the flows Q.
E(Q11) . Q
(155)
a(Q,H) . Q
E(Q,H) = A 11 Q + A l2 H = D (156)
and
imbalance function.
and (157) will not be exactly zero, this means that, if Q (i) and
11.• II
a (i) are the flow and piezometric head vectors at iteration 1,
then:
82
8E (i) : head loss imbalance per pipe at iteration "i".
head loss and nodal flow imbalances via their respective total
da = A21 dQ (161)
and
sg(i) z da (165)
83
which is presented below. An alternative derivation is proposed
n A 11 : Al2 d Q d E
(166)
A2 1 : o
[ I * [ -c-1-i I = [ " i
where the right-hand-side vector represents the head loss
1 -1 [dEl
n A 11 1 A 12
= 1
1 (169)
: ° dg
[ A21
84
-1
n A ll : Al2 [ B 11 1 B12]
,, = : (170)
[ A2 . 0 B 21 : B22
following is obtained:
d H = H( i) + ( A 21 A 11 -1 A l2) -1 a) (176)
which, when considering (162) and (163), allow us to compute
85
and
network.
Q ( ° )= - A 11 -1 A l2 H(°) (179)
Introducin g equations (162) and (163) into (160) and (161), the
head loss and nodal flow imbalances can be computed as:
and
86
and
[-- [ A-
21 A 11 -1 A l2 ] H(i+1) = - g (188)
different methods: the loop and nodal Cross' methods (single path
that the computer codes used were similar from the programming
point of view. All the programs were designed to work with the
88
techniques for matrix manipulation were used in all the programs.
100 pipes and 21 systems with more than 100 pipes. Changes in the
for networks with less and more than 100 pipes, respectively.
methods was carried out for all the systems with less than 100
Table 2.1.
89
characteristic curve was considered. No failure whatsoever was
systematic comparison.
(1969).
distribution networks.
90
c) The computational costs of the simultaneous path adjustment
5.4 and 6.9 sec. were found in the case of the larger networks,
initial solution.
91
order to establish some theoretical differences.
NETWORK CHARACTERISTICS
92
The main conclusions that can be drawn from the comparison by
also the best in terms of CPU time for initialisation (not shown
in Table 2.2).
did with his linear theory method, the mesh-nodal method requires
38% of the CPU time needed for the looped Newton-Raphson for
networks with less than 100 pipes (see Table 2.2), whereas the
Although the programs used were not exactly the same, this
app ears to suggest that the mesh-nodal performs faster than the
linear theory method. The results are even more favourable for
reaching 13% for the 443 pipes example (see Table 2.2.).
93
clearly requires less time to converge than the looped Newton-
6 hij
Looped methods : = n Kij Q1j n-1(192)
6 Qij
6 Qij 1 (1-n)/n
Nodal methods : - (Hi - Hj) (193)
6 Hi n
6 Qij Rij
= (194)
6 Hi n (Hi _
would explain why the nodal methods are found to be more trouble-
94
1
The results in terms of global CPU time in Table 2.2 can also
be exp lained in terms of equations (194) and (192). In fact, if
0 [i.e. equation
(194) is pot "blowing up"], the derivatives for the nodal methods
[equation (194)] are, in general, "stronger" than those for the
case, with a pipe of low head loss (i.e. Hi Hi, but not close
will tend to zero, since the flow is "small", thus spoiling the
al. (1978)].
95
methods [equation (194)] is the worst we can expect, because it
is subjected to cancellation errors and also to unbounded values
when the heads Hi and Hi are close enough. Expressions of this
The linear theory method of Wood and Charles (1972) should not
flows and so on, may work for some cases, but there is no
given.
96
In the limited comparison carried out by Carpentier et al.
ones.
pipes and pumps (at most), without including other devices like
97
chapters.
3.1. Introduction.
the seventies.
establishes mass balance at the node level, while the latter has
99
Indeed, as shown in Chapter Two (Section 2.3.7., equations 90
achieved for any set of piezometric heads, and is due to the fact
matrix A31 and the branch to node incidence matrix Al2 is zero,
i.e.
A31 A l2 Li = 0 , V H
due to:
A31 Al2 = 0
100
power) which controls the steady state flow in a water
for the steady state flow condition are simply the simultaneous
101
obtained via the application of the gradient operator both to the
gradient method.
Yang and Song (1979), Song and Yang (1980) and Song and Yang
102
the nodal balance and the non-linear head loss/flow relationship
1 A 11 : A l2 11 Q 1 1 -A 10 HO 1
(3)
I_ A91 : 0 i 1 H
[ g
and Alo is the equivalent branch to node matrix for the source
transpose of Al2.
either with the heads or the flows. This leads to a linear system
103
The gradient method fully recognises the fact that all the
where
: head loss in pipe joining nodes "i" and "j"; a component
network.
nodal piezometric head; a component of the (NNx1)
104
is operating (laminar or turbulent).
passes through the origin in the flow/head loss plane. The model
particular device, is the value of"K-1 J- " and "n", which are
tests.
been proposed in the literature, relating the head gain with the
105
rather than concentrating on a more mathematicall y sophisticated
model.
h•ij •n +
• = a- • Q• li 15-7-1
li (8)
distribution systems.
checks that the pump is always working with positive flows and
106
3 6 Extension of the gradient method to include PUMPS:
A ll Q A l2 H = - A 10 HO (9)
where now
and where, upon using the g eneral head loss (gain)/flow model
El l in1-1+ 131/Q1
*11
* 21 Q21 n2-1+ 132/Q2
A 11 = (10)
5NP/NP
A ll : Al2 -A10 HO
A21:0 lu a
107
Because the system represented by (3) is a non-linear one, a
direct solution is not possible, and an iterative procedure must
A21 0 d H d g
where:
c'11Q1In1-1 13 1/ Q1
Q21n2-1
c<2I 132/Q2
A ll * = = A n - (12)
aNpAipInNP-1
13NP/QNP
where Q( i ) and HU) represent the pipe flow and the nodal heads
108
at the iteration "i". A11") stands for the matrix A ll evaluated
at Q=Q (i) , but from now on we drop the iteration superscript in
the matrices All and A11* for the sake of notational simplicity.
-1
[1:1Q1 N All* : Al2
(15)
d H A2 0 dg
where A1 1 * is evaluated at
I -1
A 11 * : Al2
N [ B 11 : B 2 1
1
I = : (16)
A21 ' 0 B21 I B22
[
d =B11 dE -I- B 12 dg
(19)
d H = B 21 dE + B22 dg
109
On introducing e quations (13) and (14), and the results of
(18) into equation (19), we get:
d Q = [I- G-1 A l2 (A21 G-1A l2) -1 A 21 ]G-1 (A llEi (1)+A l2H " )+A 10 HO ) +
+ ( G-1 Al2 (A21 G-1A l2 )-1) ( A21 ii(i)-g) (20)
Similarly:
The system of equations (20) and (21) can be solved for dQ and
dH via an iterative procedure, on defining:
d Q = Q( i ) - Q(i+1) (22)
and
d H = H(i) - H(i+1) (23)
And upon introducing (22) and (23) into (20) and (21), we
obtain:
from (23):
110
Starting with any flow distribution, not necessarily
gradient operator being applied over both the heads and the
111
( start (i=0) )
4
Read data:
- General data: title, no. of links, nodes and
sources, convergence criterion, formula,
printout options.
- Link data: initial and final nodes, length,
diameter, roughness, type of link, a, 13.
- Nodal data: type of node, consumption, ground
level, fixed heads.
4
Assemble matrix of coefficients: A
= [A 21 G-1 Al2]
Assemble right-hand-side vector: h
f A 21 G-1(A 11 Q " )+A 10 11 0 ) - (A 2 Q(i)-a) 1
Set
Q(i)=1
or other
Print results
no Stop
= i+1
112
c) The system of e quations (11) can be solved explicitly by
inverting the partitioned matrix, producing the system (19)
which, after some algebra, leads to the coupled equations (26)
and (25).
d) The most computationally expensive step of the algorithm now
becomes the solution of the symmetric positive definite linear
system (26). Chapter Five deals with the efficient computer
solution of the system (26). Equation (25) is simply the product
of matrices and vectors and no linear system needs to be solved.
The inversion of matrix G is straightforward, since it is a
diagonal matrix.
e) Iterations: starting from any initial flow distribution,
equations (26) and (25), in that order, are applied recursively
until some convergence criterion is met; see Figure 3.1.
solution.
113
From a mathematical point of view, since the solution to the
think that the network analysis problem might have more than one
argument).
Most of the work done in this area has been adapted from work
114
origin.
unlimited flow).
branches.
written as:
Content: G = f
0
h dQ (28)
hij
115
Both the content, and its dual, the co-content, have units of
power and they are actually proportional to the dissipated power.
From a geometrical point of view, the content and co-content
functions can be interpreted as the area underneath and above the
h-Q curve. Millar (1951) demonstrated that both quantities
(content and co-content) are constants when steady state flow
occurs, i.e.:
dG dJ
— = — = 0 (30)
dt dt
116
does have a unique solution.
117
In the 1960's and 1970's, the question of the existence and
uni queness of the solution remained without an explicit
valves (P. R. V.'s), "the experience has shown" that the network
link real reservoirs with the reference (ground) node. Then, the
Qij
118
subject to the constraints:
where all the flows are assumed to be positive (i.e. in the pre-
119
represented by equations (34), (35) and (36), is that G(Q), the
content function given by equation (34), must be a strictly
convex function.
Using a small example network with four pumps and one check
120
variational approach, provided the proof of the uniqueness when
the functions Oii(x) in (34) are twice differentiable, strictly
not been reviewed so far, and which also fit into the class of
from the gradient method [Todini (1979), and Pilati and Todini
121
linear theory method by Wood (1981a), can also be classified as
Stimson and Brameller (1981) recognised the fact that the loop
122
Head loss per pipe: h = A l2 H A lo Ho (38)
and
A 21 Q(i+1) = a (42)
Upon introducing (40) and (41) into (42), we get:
A 21 8Q(i) = (43)
introduced as:
(44) .
where:
123
Similarl y , the flow correction vector SQ can be partitioned as:
8Q2 I
SQ = (45)
SQ3
dependent flows.
6Q 2 I
[ A 22 : A 23 ]*[ --- =0 (47)
SQ3
The importance of equation (48) lies in the fact that the set
124
issue, since (48) ensures that, starting with a balanced flow
Taylor's expansion.
Ohl
J (i) is the tangent matrix J( i ) = [-- (50)
aQ
Q=C1( i)
J = N Ali* = G (51)
Because
h(i) = A ll Q(i) (52)
equation (49) can be introduced into the expression for the
head loss (38), giving:
keeping the Taylor's expansion with all its terms, including the
125
SQ (i) = [ J(i) ] -1 f A l2 H 4. A 10 HO - A 11 Q(i) - E(i) 1 (54)
(3- 4-1) --> Q and ji( i+1 ) --> H if E (i) --> 0 (56)
and equation (55) becomes:
The flows can be updated from (54) which, when introducing the
assumptions (56), gives:
it(i+1) = {I-EJ(]- 1A 11
L 1Q(i) - [J(i)]-1{Al211(i+1)+A10H0} (61)
1
126
ii (i+1) = II 3 -N3( i) ]-1 A 3 31 Q 3( i ) _[J 3 (i) ] - 1 {A 32 1:1( i +1) + AnHol (62)
where:
off errors in the solution of the linear system (58), where they
algorithm is as follows:
127
b) Compute the tangent matrix J( i ) and the head loss vector:
main structure of the linear system for the head update H (1+1) in
use (61) for the flow updating; they use (62) which comes out
128
The integrated mesh-nodal approach requires the generation of a
129
3 9 Comparison of the gradient method with some of the
3.9.1. Introduction.
comparison carried out by Wood (1981a) [see also Wood and Rayes
theory method.
to the fact that their structures are similar; this is also due
purposes.
parameters as well.
130
linear theory method and the gradient method . The computer
characteristics:
(1981a); the code requires both the initial flow distribution and
Press et al. (1986). The program does not take advantage of the
(LOOP) was found to have the best p erformance, in the sense that
131
Because the main objective of this comparison has been to
fair comparison has been carried out. Indeed, the codes could
Section 3.6 of this chapter. The code does not need the
132
i) Example A
ii) Example B
the branched part connected at node 25 (see Fig. 3.2) has been
The steep characteristic curves have been taken from data of real
iii) Example C
133
Fi g . 3.2. Basic network for examples A, B, C and D, from Chin et
al. (1978).
134
iv) Example D
The previous example C has been modified to include 2 valves
The valves are located in such a way that, when fully closed, the
v) Example E
vi) Example F
This is a network of 13 pipes, 10 nodes and 1 source, reported
distribution used.
12
•
5 (1/s)<=0
13
V
5 (us)
P
A16 15
.6
10 (Vs)
5 (us) 5(VO
14
V
5 (Vs)
1;t1)
100m
136
3.9.4. Comparison of the algorithms for the simultaneous path
flow solution as part of the input data. For the sake of a fair
distribution networks:
(program LT) converged for all the examples, even under cases
valves (examples D and E), where the high head loss (eventually
137
Table 3.1. Summary of the performance of the programs LOOP,
LT and GRAD with test examples.
PROGRAM
EXAMPLE
LOOP LT GRAD
,
converged converged converged
A iter=16 iter=17 iter=16
t=69.53 sec. [1] t=788.51 sec. t=29.77 sec.
Notes:
138
of "infinite" value) produced unreliable piezometric heads in the
section of the network which does not remain connected to the
reference node. This is in fact a problem inherent to all loop or
path-based methods, and is produced because the network is
broken-up into two (or more) disconnected sub-networks and
consequently the path to the reference node is lost. A solution
to this problem can be produced, by re-defining the reference
node, but it is not clear if this makes the program problem-
dependent, because, in more complex networks, the disconnection
of the system may happen anywhere (i.e. due to pressure
regulating valves, maintenance work, etc.).
139
These are approximate values computed for a regular square
network.
The global gradient method did not fail to conver g e for any
140
3 10. Concluding remarks.
141
In summary, it has been shown that the gradient method for
determining the steady state flow in water supply distribution
networks, offers some advantages when compared with other
gradient-based algorithms, such as the "simultaneous path
adjustment" method and the linear theory method. The reliability
of the method, in terms of its ability to converge for both the
flows and the piezometric heads, under extreme cases based on
ill-conditioned problems, makes it desirable for optimum design
and simulation applications, where manual intervention to avoid
convergence or disconnection problems is not possible.
142
CHAPTER FOUR
4.1. Introduction.
pipes, pumps and valves. The kind of valves we have been dealing
fully closed within the network, and the state of each valve is
the operator can alter the setting of any of these valves (by
the valve type, its design and its operational status (fully
(1981).
143
= a ij Q ij n + (1)
where:
The model (1) is then able to describe any type of device with
a fixed head loss/flow characteristic curve (valves, meters,
fittings, etc.).
144
to reduce leakages).
145
one direction), but that automatically closes when the flow
exclusive modes:
outlet pressure.
where:
= constant, dependent on the valve design.
146
ii) Fully closed (active mode):
plus fittings. The main valve and the relay system are basically
the same for most of the pressure control valves; they differ
Ratcliffe (1986).
a) Main valve.
Although there are different design concepts for the valve body
The main valve body is a cast iron globe valve design, modified
element that fits into the main valve seat. Either Figures
147
1 Valve Element
2 Upper Cylinder
3 Isolating Cocks
4 Strainer
5 Orifice
6 Needle Cock
7 Indicator
8 Relay Valve
9 Diaphragm
10 Spring
11 Adjusting Screw
12 Pressure Gauge
13 Alternative Constant
Pressure Connection
14 Electric Motor
VALVE FULL OPEN AND DOWNSTREAM STATIC PRESSURE WITH VALVE SHUT
IN ACTIVE WHEN INLET w
R
SET OUTLET
PIRLEVIEJ IR?
RANGE OF GRADIENTS WITH CORRECT
UT LET PRESSURE
SET OUTLET
PRESSURE RANGE OF GRADIEN'T---
RANGE OF GRADIENTS WITH
VALVE FULL OPEN UNDER EXCESS DRAW OFF
148
1 Valve Element
2 Up p er Cylinder
3 Isolating Cocks
4 Strainer
5 Orifice
6 Needle Cock
7 Indicator
8 Relay Valve
9 Diaphragm
10 Spring
11 Adjusting Screw
12 Pressur• Gauge
a) Valve No. 1310 Standard Valve with manually operated relay valves
b) Operation
•
LOSS
THAT.
POSSIBLE GRADIENTS VALVE
WITH coRRECT
INLET PRESSURE TERMINAL UPSTREAM PRESSURE
C)
VALVE SHUT WHEN TERMINAL
UPSTREAM PRESSURE IS LESS
THAN SET INLET PRESSURE
1
SET INLET
PRESSURE
INLET
RANGE OF GRADIENTS
149
on the relative values of the pressures at the inlet, outlet and
This component (the valve element) allows the valve to move from
b) Relay system.
The inlet, outlet and upper cylinder of the valve, and the
150
valves. A change in their connectivity implies a change in the
valve body into a PRV, a PSV, altitude valve, constant flow valve
the connecting piping system, added to the fact that the pressure
setting s can be adjusted on-site or by remote control, makes this
151
PRV. This operational mode is achieved in two situations:
ii) The outlet pressure (P2 in Fig. 4.1.a.) is lower than the
setting pressure, but greater than the inlet pressure (P1 in Fig.
4.1.a.). The relay valve keeps open and P1 < P3 < P2, then the
pressure in the upper cylinder is greater than the inlet one (P1)
this case the PRV operates fully open in its "inactive mode",
pressure, the relay valve remains open and, then, P1 > P3 > P2
which means that the valve element tends to stay fully open.
local head loss, while, when in the active mode and check valve
152
Another way of understanding the PRV's behaviour, is by
In the operating point "1" of Fig. 4.4 (Hi < H setting) the
valve is fully open, and the difference between the inlet and
outlet pressures depends only on the flow (or velocity) and the
heads.
ill is equal to the local head loss produced through the valve in
153
P. R. V. in the ACTIVE or
h ii = H• —H
I i CHECK VALVE operating mode
Qij
H setting
Hi
inlet head
H setting H setting + hl
..
1( >
Valve fully open Valve partly closed
Fi g . 4.4. Relationship between inlet and outlet heads in a
pressure reducing valve. (For a constant flow).
154
Figure 4.5 is valid for positive constant flows only. In the
When Hi < Hsetting + hl, the PRV is in the active mode (fully
= Hi - Hj = ao Qijn
where
For Hi > Hsetting + hl, the valve increases its aij , i.e. we
are in the active mode, and the relationship between hij and aij
the flow Qij are not independent, and the flow is not constant,
in the case of Figure 4.5, for showing the non linear effects,
155
but its practical usefulness is negligible. Figure 4.6 has been
obtained from a system where the reservoir is not connected
directly to the valve inlet, but through a network; then, the
variable inlet pressure is due to the head losses throughout the
network. Starting from zero flow, we move from the point "3" (in
Fig. 4.6) to the left.
The main difference between Figures 4.6 and 4.4 is that the
ascending limb of the curve is no longer a straight line, since
now it has to include the effect of the variation of flow.
156
•
Hsettin9.,
Hi
Hi
outlet head
•
hii= — H setting
el 3
H setting as
2
.1= °CoOrilj
Q>0
Hi
inlet head
•
H setting H setting + hi Hi max.
157
In the case when the outlet pressure is dependent on the head
loss through the valve, Fig. 4.7 clearly suggests how can we find
Hsetting.
subsequent sections.
following duties:
158
H-
J
downstream head
4
1
/ outlet pressure
independent of
valve resistance
H--+
H setting n ..
a . outlet pressure
is dependent on
valve resistance
4.2.a.), until the head losses within the valve produce the
i) The inlet pressure (P1, see Fig. 4.2.a.) is lower than the
pressure setting; then, the difference between spring tension and
ii) The outlet pressure (P2 in Fig. 4.2.a.) is greater than the
inlet pressure (P1 in Fig. 4.2.a.) and lower than the inlet
pressure setting. The relay valve keeps open and P1 < P3 < P2,
then the pressure in the upper cylinder is greater than the inlet
than the pressure setting, and greater than the outlet pressure.
In this case the PSV operates fully open in its "inactive mode",
pressure the relay valve remains open and, then, P1 > P3 > P2
which means that the valve element tends to stay fully open.
160
represent graphically. Figure 4.8 is the equivalent of Figure 4.4
for the case of a PSV. We shall O m it the rest of the figures,
the response of the system at the inlet of the PSV versus the
4.9.
actually modelled the set PRY plus the pipe in which the device
was inserted and, for the case when the PRV is next to the
where:
outlet head
H setting
H
inlet head
upstream head
inlet pressure
b
independent of
valve resistance
inlet pressure
is dependent on
valve resistance
H sew"
(Kips
162
When Hj < Hi < H setting , then H setting must be re p laced by Hi
in equation (1), and the same model can be used to represent the
PRV's behaviour. It is implicit that the valve closes when Hi<Hj
or Hj > Hsetting (i.e. Qij = 0.0).
163
constant outlet pressure equal to ps, p rovided that pi > ps.
Jeppson (1976) and Jeppson and Davis (1976), reviewed all the
Fig. 4.10.
Hsetting
164
Table 4.1. Summary of effects of PRV's in different
formulations of network analysis equations
[following prop osition of Jeppson (1976) and Jeppson
and Davis (1976)).
165
As it is clear from Table 4.1 that every PRV implies that the
original pipe (from nodes "i"- "j") is split up into two parts,
we need to restore one loop equation for every PRV introduced.
For this purpose, it is necessary to link every PRV with some
existing (real) reservoir, create a path between them, and write
down the equation establishing that the head loss in this path
must be equal to the difference between the pseudo-reservoir head
(Hsetting corresponding to the PRV) and the head of the real
reservoir.
Q ix = ((Hi-Hx)/aix)1/n (7)
Qyi = ((Hsetting41p/ayj)1/n (8)
then, from equation (6):
[(Hi-Hx)/a ix] lin = rtu
L, -setting-H Wayil l /n(9)
or
(Hi - Hx)/ a ix = (H setting - (10)
H)/ c(Yj
166
Lekane (1979) followed a similar approach to that proposed by
Zarghamee (1971) in modelling PRV's within a nodal formulation.
He also reported the results of the application of the method to
a small network ( 18 nodes and 24 branches with 2 PRV's).
Another application to a bigger network (314 nodes and 424
branches) is also mentioned.
167
carried out at the beginning of the iterations or after having
reached a certain degree of convergence ( if so, what is the
"ri ght" degree of convergence ?); clearly, because in the first
iterations, the network can be working in conditions that are
very far from the final ones, the checking procedure can lead to
the models for introducing NW's, CV's and booster pumps into the
the PRV setting. This author also gave some warnings about the
168
depends on the heads at the nodes, something which is determined
at the end of the process in the loop approach. Consequently,
Gessler recommended a nodal formulation. He also gave a couple of
detailed examples on how to include PRV's; we agree with his
results, as discussed in a subsequent section. Finally he
addressed the problem of uniqueness, recommending a reassessment
of the physical problem whenever mathematical and numerical
problems are found.
All the referenced authors deal only with PRV's, without even
mentioning PSV's or other pressure controlling devices, like
altitude valves, constant flow valves, etc. As it can be seen
later on, in the section corresponding to the examples found in
the literature, not a single example dealing with PSV's was
reported and, more importantly, there is no example dealing with
PRV's and PSV's simultaneously. The reason for this may lie in
the fact that, roughly speaking, PSV's are dealt within a similar way
to PRV's, but instead of the maximum downstream pressure we are
concerned with the minimum uPstream pressure. No reason
whatsoever can explain the lack of evidence on the behaviour of
real networks having more sophisticated controlling devices.
169
All the authors reviewed, with the exceptions of Jeppson (1976)
and Jeppson and Davis (1976), referred to a nodal formulation of
the set of equations. The attempt of Jeppson and Davis is
commendable, in the sense that they focused on all the possible
formulations.
170
because its downstream head is found to be greater than the
necessary in the case of the PRV's, when sometimes (in the nodal
clear that the fact that all these valves are variable resistance
checked.
that we can fix the PRV setting does not necessarily mean that
172
solution can be obtained. We believe that because pressure
controlling valves are in essence variable resistance devices,
consequently dissipative devices, their steady state solution
must be unique. What we shall try to do in our model is to
follow, as closely as possible, what we think should be the
response of the physical water distribution network.
173
to STATUS = 2 in the program.
on/off pump switches, with minor changes in the data, and without
X = IH j - Hil
higher than the inlet head. Also, we determine the two maximum
next most imbalanced one (if any) are known. When there is only
174
one valve out of balance, X* is set equal to the accuracy
between the outlet head and the setting for a PRV (or inlet head
with the setting. For CV's the correction stage simply changes
on the lumped curve presented in section 4.2 (Figure 4.7 for the
PRV's and Figure 4.9 for the PSV's). This explanation will
relate to the case of the PRV's only, because the case of PSV's
is completely similar.
suggests that, for finding a*, the a which matches outlet head
175
represented by the lumped curve is not known, the derivative of
corresponding head. Let us call the new outlet head H2, of course
needed, which means that we have to run the main program again to
176
but this was not found worthwhile during the development of the
program. If we find that that we have to stop the Newton
iterations, we go back to the detection stage, in order to check
if any other device has gone out of balance. Checking the valve
status is not as computationally expensive as re-running the
program.
H outlet
H1
H2
X -1 X2
H3
H4
Hsetting
177
In summary, the method for determining the resistance parameter
in Figure 4.12.
I
I.( fixed ai Program modelling
Standard network pressure regulating
analysis algorithm devices: performs
(gradient method detection and correction
or ) stages, determining ai
Heads: by Newton's method.
178
To establish the robustness of our proposed model, we have rtAn
results. At the end of Table 4.2. some of our own test examples
are also included.
including one PRV, was analysed by Jeppson and Davis (1976) under
two possible settings for the PRV:
a) H setting = 149 m
b) H setting = 140 m
179
Table 4.2. Comparison of results reported in the literature and
our results for networks containing pressure control
devices.
180
Fig. 4.13. Example JEPPO, from Jeppson and Davis (1976).
According to our findings, the PRY closes for any setting under
151.34 m. This is the limit where the downstream head of the PRV
is determined by the system, rather than by the PRV
characteristics.
182
4.5.2. Example JEPPl. from JePPson (1976).
4.14 and the comparison of the results obtained here with those
of Jeppson is shown in Table 4.4.
30.48 m
183
Table 4.4. Comparison of the results given by Jeppson (1976)
and our results for example JEPPl.
184
node 9) above is correct. Had the solution given a
negative flow rate in pipe 6, this assumption would
be incorrect, since the PRV would then have acted as
of this network is shown in Fig. 4.15, and Table 4.5 presents the
Figure 4.16 shows the schematic of the network and Table 4.6
presents the comparison between Jeppson's results and ours.
300 m
300 m
187
Table 4.6. Comparison of the results given by Jeppson (1976)
and our results, for example JEPP13.
188
4 5 5. Example COLLINS. from Collins (1980)
the PRV's are always acting as. non-return (CV) valves; thus, the
189
Table 4.7. Network solutions for example COLLINS, from Collins
(1980).
-
Operating modes Network
M of control elements solution (a)
o
D PRV's Pumps Head (m) Flowrates (l/s)
E status P 1P 2 P 3 H6 Q26 Q46 Q56 Q67
1 Passive on on off 238.0476 3.5813 0.8927 0.0000 4.4739
2 Passive on on on No solution
3 Passive on off off 190.0000 2.0000 0.0000 0.0000 2.0000
4 Passive off on off No solution
5 Passive off off off No solution
6 Passive off on on No solution
7 Passive off off on No solution
8 Passive on off on No solution
9 Active on on off No solution
10 Active on off off No solution
11 Active off on off No solution
12 Active off off off No solution
13 Active on off on No solution
14 Active off on on No solution
15 Active off off on No solution
16 Active on on on 245.0000 --(b) --(b) 1.0000 4.7258
_
Notes: (;,6"- 1 A
#.7,12 = Q26 ; Q34 = Q46
t ° 1 '426 + Q46 = 3.7526
215 rn
190
possible operating modes suggested by Collins are reduced to
eight cases (see Table 4.8), this number being the possible
flow can be possible, unless we accept that the pumps can operate
first seven modes defined above, assuming that the PRV's are set
to a head of 245 m in their outlet (node 6).
with the conditions imposed to the problem. Only in his mode 16,
although he finds that the flows between nodes 2,6 and between
4,6 admit more than one solution. In our case we found that this
1 on on on
2 on on off
3 on off on
4 on off off
5 off on on
6 off on off
7 off off on
8 off off off
,
191
solution in also unique, the flows being controlled by the
difference in the levels (available potential energy of the
reservoirs connected to nodes 3 and 1, see Fig.4.17). In his mode
1, which is equivalent to our case 2, the head in node 6 was
238.0476 m, while we found that this head is the PRV's setting.
In his mode 3, equivalent to our case 4, Collins found that the
head in node 6 is 190.0000 while we obtained 205.000. We
reproduced exactly Collins results in modes 3 and 6, provided
that we change the settings of the PRV's to the heads found by
Collins, but this is clearly a change in the original settings
adopted for the PRV and, therefore, is another problem.
192
In our opinion, what is happening in this example is a
193
there is no point in having a fast program leading to the wrong
solution.
We also have to emphasise the fact that the algorithm does not
different settings, one for the day and the other for the night),
194
1 Valve Element
2 Upper Cylinder
3 Isolating Cocks
4 Strainer
5 Orifice
6 Needle Cock
7 Indicator
8 Relay Valve AS FOR 1301
9 Oiaphragm
10 Spring
11 Adjusting Screw
12 Pressure Gauge
13 Alternative Constant
Pressure Connection
14 Electric Motor
15 Solenoid operated valve
1 Valve Element
2 Upper Cylinder
3 Isolating Cocks
4 Strainer
5 Orifice
6 Needle Cocks
7 Indicator
8 Riley Valve
9 Diaphragm
10 Spring
11 Adjusting Screw
12 Bill Cock
13 Non-Return Valve
b)
Valve No. 1330 Standard Altitude Valve — Flow into tank only
Valve No. 1331 Altitude Valve allowing flow in and out of tank
(Illustrated by dotted lines)
1 Valve Element
2 Upper Cylinder
3 Isolating Cocks
4 Strainer
5 Pressure Orifice
6 Needle Cock
7 Indicator
8 Relay Valve
flow 9 Diaphragm
=-4.-
10 Spring
11 Adjusting Screw
12 Flow Orifice
c)
Valve No. 1340 Standard Constant Flow Valve
195
CHAPTER FIVE
flows through the data set. This is possible because the coupled
196
equations which are the basis of the gradient method
equations are:
Heads update:
Flows update:
that, with the heads computed via (1), the next flows Q( i+1) will
197
Thus, after the first iteration, equation (2) produces a
balanced flow distribution [Q( 1+1 )] which, when used during the
imbalance to vanish.
devices (see Chapter Four for details), where the previous state
198
As a result, the problems of generating an initial flow
199
possible, leadin g to more or less efficient factorizations,
the data structures required are set up; this refers to all the
comprising:
200
- Band methods.
them are based on the observation that fill-in occurs within the
the objective is to
reduce this "area".
A=
4m.
201
Although, for general matrices, the envelope methods are better
than the band methods, the main shortcoming is concerned with the
in the Gauss elimination process, from the row (or column) with
been quite bad, particularly for networks larger than 900 nodes,
where it actually failed.
202
original matrix are explicitly needed, the rest being handled
implicitly.
graph through the shorter path (i.e. in only one direction). This
the algorithm.
systems of equations.
definite systems.
( L -1 A LT ) ( L T x ) = L -1 12 (6)
where the new transformed matrix ( L -1 A L -T ) is still
symmetric, and is better conditioned than the original one. The
204
same original conjugate gradient algorithm, when applied to the
A z L LT (6)
where L has exactly the same sparsity pattern as A. The non-zeros
outside the original pattern are simply thrown away, thus solving
the fill-in problem. Note that with (6), the transformed matrix
( L -1 A L -T ) in (5) becomes nearly an identity matrix, thus
In sections 5.2.1. and 5.2.2. the two main alternatives for the
205
- Nested dissection method.
- Multifrontal method.
- Preconditioned conjugate gradient method.
approaches.
206
Table 5.1. presents the main characteristics of the examples
used in the numerical comparison, while Table 5.2 presents the
results, including the execution times of the tested programs in
the Amdahl 5860 mainframe computer at Newcastle University; these
NETWORK DATA
Data file Branches Nodes Sources Equations CV PRV PSV
net5.dat 74 48 2 46 0 0 0
[1]
prv5.dat 74 48 4 44 1 2 2
[2]
bofn.dat 298 266 26 240 0 0 0
net50.dat 420 225 2 223 0 0 0
[4]
net51.dat 760 400 2 398 0 0 0
[4]
net52.dat 1200 625 4 621 0 0 0
[4]
net53.dat 1740 900 4 896 0 0 0
[4]
net54.dat 4900 2500 4 2496 0 0 0
[4]
net56.dat 9660 4900 4 4896 0 0 0
[4]
Notes:
[1] : Data from real network, published by Chin et al. (1978),
"Solution of water networks by sparse matrix methods ",
Int. J. Num. Meth. Eng., 12, 1261-1277, (1978).
[2] : Based on the data of previous example, with pressure
regulating devices added.
[3] : Real main network of Bogota, Colombia. Source: G. Gonzalez
"The gradient method", M. Sc. Dissertation, Civil Eng.
Department, University of Newcastle-upon-Tyne, 1987.
[4] : Synthetic network. Network generated as a square mesh with
SQRT(Nodes) vertices per side.
207
Table 5.2. Comparison of execution times and number of
iterations required by different linear solvers
in the gradient method for pipe network analysis.
208
execution times are, in general, average times of at least a
couple of runs and correspond to the time needed for the solution
accounts for all the extra auxiliary vectors and pointers needed
to implement the algorithm.
storage shown in Table 5.3 and Fig. 5.3 refers to that required
by the linear solver itself; it does not include that of the non-
for use with regular finite element meshes. This is a point that
209
Table 5.3. Comparison of storage requirements (Bytes) for
different linear solvers used with the gradient
method for pipe network analysis.
2,576 [3] 2,464 13,440 12,488 22,288 34,776 50,176 139,776 274,176
ICCG 1,748 [4] 1,672 9,120 8,474 15,124 23,598 34,048 94,848 186,048
------
4,328 [5] 4,136 22,560 20,962 37,412 53,374 84,224 234,624 460,224
Notes:
[1] : Refers to the linear solver algorithm:
ICES : modified Conjugate Gradient method, with Incomplete Choleski factorization ( Kershaw factorization).
RCM : envelope ("skyline') method with Reverse Cuthill-McKee ordering.
(16) : Datient Minimum Degree algorithm.
ROT : Refined Ouotient Tree algorithm.
1WD : one-Way Dissection method.
ND : Nested Dissection algorithm.
MA27 : Harwell subroutine MA27.
[2] : Corresponds with data sets in Table 5.1.
[3] : Upper figure corresponds to the primary storage (Bytes).
[4] : Corresponds to the overhead storage (Bytes).
[5] : Corresponds to the total storage: primary 4. overhead storage (Bytes).
[6] : No solution for the first iteration after 970.0 secs., program was interrupted and cancelled. No attempt was made with
the larger examples NET54.DAT and NET56.DAT.
.c
4-)
••••n
•••••n
• n••••••
211
i
0)
C
0
i
I- (N 0 0 47-;
0 < 0 0
c .1-- 0 0
0
D
0-
0 CD
*
lt
11)
% Lf-
1
%
%
0
1
. 1
a)
%. 1 ..0
0
%
.. 1
E
n
Ct
.
. 0z n
.
.
1
1
0
0
.
. 1 d- as
.
. 1 .4
4-)
. 1
% 1 e4.4
.
0
r.t4
o o o o o
o o o o o
IL) o in o to
N C\1 ,-- .—
212
There is an additional point that is worth noting, which is
convergence.
of the methods, i.e. "good" methods are still good when a more
213
the algorithms, while the "true" performance of the algorithms in
etc.).
respect. The total storage (see Fig. 5.3) behaves linearly for
214
execution time viewpoint, followed by the dissection algorithms
of George and Liu (1981). The nested dissection is about 100%
already taking place, makes the case for direct solvers even
stronger.
215
5.3. Concluding remarks.
gradient).
216
CHAPTER SIX
6.1. Introduction.
since the early sixties, when these machines became available for
217
calibration being responsible for guaranteeing that the results
is going to be applied.
is due to the fact that demands are always changing, pumps are
where hij and Qij are the head loss and flow between nodes
dependent.
219
The non-linear model represented by equation (1) is valid for
6,000). For laminar flow (i.e. Reynolds number less than 2,000),
In the case of valves, Bij z 0 and aij has to be obtained from the
220
to the uncertainty in the determination of the characteristic
a) Pipe resistance:
The roughness of a pipe, represented by Cii in the case of the
The older the pipe, the greater the uncertainty associated with
its resistance.
b) Diameter:
materials of the pipes are known (which is not always the case in
221
c) Length:
This variable is usually determined from a map of the city
under study (scales 1:5000 or 1:10000 are usual), where the pipes
d) Valve settings:
The valve characteristic parameters are also subject to errors
uncertain parameters.
due to the normal wear during the operation of the pumps or due
f) Minor losses:
222
fittings are present. For these cases, the introduction of an
223
there is also a daily variation, due to the fact that the
blurred.
model, thus reducing the effort needed to solve the problem. The
224
reduction was obtained via the simple elimination from the model
the network in order to lower peak demands and also to reduce the
modelling practice.
225
* Misinterpretation of recorded data: wrong dimensions and units
e.g. open in the real network, but closed in the model, or vice
versa.
* Pipes by-passing pressure regulating valves and meters.
procedure.
6 3. Measurements available in a real network and their
accuracy.
network:
field measurements) .
* effective diameters.
227
To simplify the terminology, pressure and water level
measurements will - from now on - be referred to as "piezometric
network measurements.
not expect to have more than 10% of the nodes of the network with
228
important consumers (or sectors of the network), etc.
229
be calibrated when the results obtained with the model match
flow and head difference should be less than 5%) or, we could ask
5%), whereas the rest (say 20%) should be within a bigger error
meaningful analysis
230
combinations of operating modes of devices like pumps, regulating
previous context.
231
chapter 1)]. These more severe prerequisites can be expressed as:
implementation on microcomputers.
calibration.
distribution network, these being the most widely used models for
in the geometrical sense; this means that both the shape and the
232
at pumping stations, etc.) are replicated by the model. Reduced
cases, the way to deal with the calibration problem in the early
achieved.
233
non-linear equations is produced, now containing heads,
found.
234
linear equations in the resistance corrections, which is
scenarios.
Donachie (1974) also shows how the same procedure can be used
that look like a good calibration, when this may not be the
235
case.
* The method considers that the pressure measurements are exact
calibration.
now the time varying nodal consumptions and reservoir levels are
adjusted to match model results with the real network. The new
that the SPT algorithm has been previously used to calibrate the
resistance p arameters. Thus, the remaining uncertainties in the
236
dynamic case are the reservoir levels and the nodal consumptions.
because the SPT algorithm does not reproduce exactly the level
squares of the errors. This needs to be carried out for one nodal
valid for the SPT and DPT algorithms. In the case of the DPT, we
sensitive demand.
Donachie (1974) and Rahal et al (1980) " do not make full use of
237
the available information and can lead to results which disagree
238
location of the "head generators"; additionally, they found that
the problem is solvable only under certain conditions, otherwise
the new problem either has no solution or has many solutions.
Almost the same arguments as in the previous methods are
applicable here, i.e. lack of capability for handling the error
of the pressure measurements, no use of flow measurements, etc.
239
arbitrary and many solutions are possible, or none at all, as
shown by Gofman and Rodeh (1981).
recognised that both the pipe resistance parameters and the nodal
and measured heads and flows. His approach is pot based on re-
240
pressures in the hydrant before and during the test, plus the
levels in the upstream reservoir and the flows measured at the
variables.
241
which renders the method better suited for pie-shaped networks
with a central reservoir feeding a nearly circular area. Another
can link this reference node with any of the pressure measured
242
losses and the measured total head loss between the reference and
zone and the corrections in the demands per zone are distributed
are basically the set of non linear equations relating head loss
and flows per branch, the mass balance per node equations and
243
viewpoint, Lansey does not include the set of non linear head
constraint set.
energy equations) and also by minimum and maximum bounds for the
244
I
pipes and 197 nodes) and the other corresponds to the Federally
nodes); both systems link together. The authors did not present
All the methods reviewed in the previous section for full scale
that most of them are not able to handle explicitly the errors
estimated parameters.
245
problem, as presented in Section 6.4.
246
techniques as well and we will compare them in order to evaluate
heads and pipe flows). In this case the reservoir levels and
247
First, let us assume that we know exactly the nodal demands and
their spatial distribution at a certain time instant; this
[see Brandon (1984)] for nodal demand allocation and/or using the
the model.
single pipe joining node "i" and "j" then, in order to estimate
248
the flow in the Pipe. With these measurements, say Ht - HI and
Ideally, having measurements for all the head losses and flows
per pipe in a larger network, the application of equation (2)
all the nodal piezometric heads and pipe flows, we have to find
the best estimates of them and use equation (2) to produce our
measurements.
Note that if, in our single pipe system with nodes "i" and "j",
Once the a- 1J
-* have been estimated, using equation (2), the
249
e is explicitly related to the estimated characteristic
resistance parameter aij* via the followin g relationship:
allow C-1 J-* to move within an "allowable band" around the initial
where:
C-
1J -(°) : initial roughness estimate for pipe "i-j".
has been set to 0.10, thus allowing Cij* to move within 90 and
In theory, when having head and flow estimates near the true
250
larg e as possible).
holds for every pipe joining nodes "i" and "j", we can write:
a. .t = h . .t /(Q . .tln
ij lj lj • (5)
where hii t is the true head loss, i.e.: hii t = Hi t - Hit.
Similarly, for the estimated variables, a*, Ii* and Q*, the
where hii* is the estimated head loss, i.e.: hii* = Hi* - Hi*
a.1Jj-* h. .* /(Q.lj* * ) n
---- =
a. .t hiit /(Qii t ) n
2.3
* h li
- -*/h--t
1J
(7)
a..ttal, . 4C /on. .tln
13 L'-'13 ''"1J J
[
a* ___> at
if h* ___> ht - (5)
and if Q* ___> iit
251
It is not difficult to see that equation (7) also holds for the
case when:
E[h* ] = K ht (9)
and
E[Q* ] = K Qt(10)
where K is a constant.
can be equal to the true parameters when both heads and flow
algorithm.
calibration Problem.
252
technique known as "Kriging", which is able to explicitly handle
noisy data and to produce an estimate of the error of its own
assumption that the modelled flows are close to the true flows is
estimates.
253
Start
i=1
assign:
C--
1J 1J
no i.e.:
c(i) <-- c*
yes
( Stop )
254
provided that we input optimal estimates of the nodal piezometric
heads and pipe flows, as indicated in (8). Additionall y , because
the estimation processes (of heads and flows) could be carried
future.
6.9.1. Introduction.
255
looped networks, in general). Another situation where Kalman
filters have been successfully used in water distribution systems
256
associated with the measurements cant be handled and, as a
6.9.2.1. Introduction.
257
Hamberg and Shamir (1988 a,b) have used the same concept for
level near the borders of the network. This may be valid for a
contains all the pipes and nodes fed from this particular source.
The piezometric plane actually exists only over the pipes, and
does not exist over the rest of the space but, since we are only
258
particular points of the space. On the other hand, current
models assume that the head loss between two connected nodes is
measurements taken not only over the nodes, but also measurements
unmeasured nodes.
259
network.
We are aware that Kriging might not be the best solution when
are located at key points of the network and that they are
(i.e. a hill, a hole, etc.), we shall miss them when the time
words, this means that any estimation method (for the piezometric
260
Following Appendix B, the main steps in estimating the
unmeasured piezometric heads using Kriging can be summarised as
follows:
261
set of unknowns weights X's and P's in equation (62) or (63) of
E
j
ko j k
k
K(Xi-Xj) + E Pk P (Xi) = K(Xi - X0)J
1=1,...,n
and j=1,...,n
E'-'
i
k'i k
(Xi) = P (Xo)] k=1, . .. , m
K 11 K 12 K 13 • 1 D xol
• • Ir.ln -12— P1k K10
xo2
K 21 K 22 K 23 . . • K in P 22 ••• P2k K20
1: • : . • • :
• • . . :
• •••
'Ikon
Kn1 Kn2 Kn3 • • Knn 1Pn2 — Pnk KnO
.1 0 0... 0 * =
1 (12)
1 1 1. -411
P21
. P 22
. P 2
.3 • • P 2n 9 0 ...
• • •
0
:
-P2 P20
: • • •
• • Pmn 0 0 .. 0 -P m Pm0
Pm1 Pm2 Pm3 • -
where
= K(xi - xj)
and
Pij = Pi(Xj)
262
adequate.
Z0* = E X0 i Zi (13)
where:
data set.
weights [Vs and P's in equations (11) and (12)] , the Kriging
that shown in Fig. 6.1, in fact, with the exception of the input
are there
more points to b
estimated yes
no
( Stop )
264
6.9.3. Estimating the piezometric heads using a deterministic
6.9.3.1. Introduction.
The information contained in the raw network model has not been
265
a) Read the data corresponding to the piezometric head
measurements and the initial values for the pipe roughnesses.
pipes.
each reservoir (if more than one). This gives us the paths
between nodes has been replaced by the head loss. The reason why
the pipe head loss has been chosen as the criterion for
following paragraphs.
of those nodes which are included in more than one spanning tree,
266
necessaril y a reservoir) or a previously interpolated node. When
information on the accuracy of the head measurements is
first, so that the lower accuracy nodes (and the unmeasured ones)
The minimum head loss criterion has been used to generate the
f) Having the paths which cover all the nodes of the network,
next section.
that shown in Fig. 6.1, in fact, with the exception of the input
are there
more paths to be
interpolated yes
no
( Stop)
as in Fig. 6.4.
estimated node), we search along the path for the next measured
node, defining a reach between two known head nodes [see Fig.
the reach (DIF1 in Fig. 6.4) and also that corresponding to the
269
Piezometric
Hood
(1.4)
Fig. 6.4. Minimum head loss path of nodes, before the head
interpolation has been carried out.
PiezcmetrIc
Hood
(H)
Length (14)
500 1000
1 3 5 9 16 20 31 19 Node Number
Fig. 6.5. Minimum head loss path of nodes after the head
interpolation has been carried out.
270
The result of this one-dimensional interpolation procedure, for
The algorithm relies on the fact that the piezometric head has
considered as indispensable.
characteristics.
6.9.4.1. Introduction.
random).
heads has been through a set of broken straight lines, with the
272
representing the piezometric heads. It can be argued that, as far
as the estimation of the piezometric heads is concerned, Kriging
In the search for such an approach, and having in mind the need
interpolation method.
the following:
273
a) Read the field data points representing the measured nodal
piezometric heads:
directions
[AT A] r = AT f (15)
where:
h+4 k+4
E E rij Mi(xr) Nj(yr) = f(xr, yr) = fr (16)
1=1 j=1 r = 1,2, ..., m
274
r : (h+4)(k+4)xl column vector of unknowns, with components rij
(x l , y l ) 1 = 1, 2, ..., n
h+4 k+4
s(xl, yl) = E E rij Mi(xl) Nj(yi) (17)
i=1 j=1
275
Read field data points:
location and measured
piezometric heads.
1 ,
Compute the estimate of the
unmeasured piezometric head:
s(x,y)= E E rij
M(x)N(y) 4
i j
and the variance, if needed
are there
more points to be
estimated yes
?
no
CStop )
276
6.10. Examples of applications of the proposed calibration
method.
techniques.
i) Example A:
the following two examples. Table 6.1. b) shows the nodal data.
Note that the nodal consumptions for this example (example A) are
the same as that for the third (example C), and that they
ii) Example B:
iv) Example D:
characteristics of this network are such that the length and the
is 1.0 (1/s) for all non-reservoir nodes and the ground level is 10
m. for all the nodes. This basic network, with minor changes, is
and F.
65 m
.1.
6 1
9 e:f.‘
n w 14
. 0
Fig. 6.7. Network for Examples A, B and C.
278
Table 6.1. Network data for Examples A, B and C.
a) Pipe data:
279
In Example D only one reservoir feeds the whole network and is
v) Example E:
that four reservoirs are feeding the network and they are
the reservoir levels are: 100.0, 90.0, 90.0 and 110.0 m.,
vi) Example F:
valves, separating the network into high and low pressure zones.
50.0 and 40.0 m., respectively). The node and link labelling has
zones. The P.R.V.'s are located at links 175, 176, 177 and 178
1J=
1
0.111 I I I 1). I I I I lig I lard I I 1111 I I I
0:0 30 2:031=1=0
0•
. D,SIMITIN411
1
crairaterglimill
ess o e
57
75
Rind
1795
G
N
'211
IS. 114
e
133
0
on ma
CD
171
281
Fig. 6.9. Network for Example E.
282
110.0 M. 50.0 M.
283
Due to the uniformity of the data for the last 3 network
Williams roughnesses were set to 80, 85, 90, ..., 160 for the 16
Examples D, E and F.
284
Table 6.2. Measurement data for Examples A, B and C.
a) Piezometric head (m) measurement data:
Example A and C Example B
Node Piez. Head Node Piez. Head
1 65.000000 1 65.000000
4 62.818108 4 56.614697
9 60.988101 9 47.368164
11 60.988960 11 47.364705
Note: flowrates are considered positive when they flow from the
initial into the final node of the pipe (see Table 6.1.a)
285
Table 6.3. Measurement data for Examples D, E and F.
(Continuation).
b) Flow (l/s) measurement data:
286
6 10.4. Defining the initially assumed roughnesses for the
network examples.
proceeded to generate a new set of C's, with the same mean as the
true values, but with some added Gaussian noise. This was done
for example C, which have been obtained for the initially assumed
true values, but they have been over-estimated by 20%. Table 6.6.
presents the initial C's assumed for examples D, E and F, for low
Uncertainty in C
Low High
Pipe
Values of C
1 81.4482 73.7684
2 89.6747 88.9714
3 95.4802 96.5187
4 98.6072 95.5956
5 103.7588 101.0750
6 116.2208 129.6719
7 116.7340 120.4833
8 121.4841 124.6930
9 126.4144 129.4729
10 130.4846 131.5323
11 133.3617 129.8191
12 137.0585 130.6980
13 148.0070 154.5090
14 152.2471 157.1059
15 153.0364 148.7907
16 162.2371 167.0742
Mean 122.8909 123.7362
Variance 596.1455 702.9258
For the true C's:
Mean 122.5000 122.5000
Variance 566.6667 566.6667
288
Table 6.5. Initial values of C (Hazen-Williams) used
for testing the calibration algorithm in Example C.
.
Pipe Values of C
1 97.738
2 107.610
3 114.576
4 118.328
5 124.511
6 139.465
7 140.081
8 141.781
9 151.697
10 156.582
11 160.034
12 164.471
13 177.608
14 182.696
15 183.643
16 194.684
.nn
Mean 147.219
Variance 860.341
289
Table 6.6. Initial values of C (Hazen-Williams) used for testing
the calibration algorithm. For Examples D, E and F,
and for low and high uncertainty.
290
Table 6.6. Initial values of C (Hazen-Williams) used for testing
the calibration algorithm. For Examples C, D and E, and
for low and high uncertainty. (Continuation).
Summary:
*******
For Examples D and E:
Low High
Mean 139.959 139.9418
Variance 7.499 14.9984
For example F:
Low High
Mean : 139.9697 139.9572
Variance : 7.4917 14.9832
For the true C's:
Mean 11X000
Variance 0.000 11410!18000
0.0000
____ \ n
291
6.10.5. Defining a perturbed set of nodal demands.
and F, where the nodal consumption is 1.0 (l/s) for all the
Table 6.7. Modified nodal demands (us), used to study the impact
of bad demand estimation on the calibration algorithm.
Examples A B,and C.
Demands (1/s)
the estimates of the pipe roughnesses, but also the flows per
i) Case I.
294
case establishing the lowest possible level in terms of flow
measurement data.
and C) due to their small size (11 nodes), the proportion has
by the fact that the networks are fairly regular both in their
295
ii) Case II.
To the same amount and distribution of piezometric head
The flow measurements are located in exactly the same pipes for
are concerned.
296
v) Case V.
In order to evaluate the behaviour of the calibration algorithm
characteristics.
of each network and for each one of the five Cases studied, has
nodal consumption data, these two sets of head estimates are the
same for all the five Cases studied, and they depend only on the
297
assumed to be available, the head estimates may be different from
case to case. Because Cases I, II and III are based on the same
have the same head estimates, but those estimates are different
run the network analysis program for the raw model corresponding
0 Calibration stage.
298
calibration results, corresponding to each one of the piezometric
numerical results are presented, both for the head estimation and
variable.
299
This performance index Rx2 , which is represented in Fig.
6.11, has the following properties:
ii) Rx2 =1, if the true value has been determined (i.e. when
X=Xt).
the variable, its variance, the maximum residual and the variance
improves all those parameters, over all the tested cases. Hence,
300
Rx2
1.0
0.0 X
and III are based on the same initial roughness estimate, their
301
are exactly the same, but they differ from the head estimates for
splines are used as head estimators, where the head estimates are
the same for all the five study cases.
and III, and for the three different head estimation procedures,
interpolation method. Both Table 6.9. and Fig. 6.12 show how the
Note that, from now on, when dealing with the residuals, we are
302
Table 6.9. Comparison between estimated piezometric heads and
its residuals, for exam p le A.
303
Table 6.10. Summary of the comparison between different
piezometric head estimation procedures.
EXAMPLE PROCEDURE AVERAGE VARIANCE STAN.DEV. MAXIMUM AVERAGE VARIANCE RATIO AVERAGES
RES.[1] RES.[2] RES. [3] EST/TRLE (4]
NOTES:
===
(1] : MAXIMUM ABSOLUTE VALUE RESIDUAL = MAX.(ABS(TRUE-ESTIMATED)1
[2] : AVERAGE OF THE RESIDUALS = (SUMMATION OF ABSOLUTE VALUES OF THE RESIDUALS)IN
[3] : VARIANCE OF THE ABSOLUTE VALLE OF THE RESIDUALS = ( SSAVR - N I AVRtt2 1 / (N-2)
SSAVR : SUMMATION SQUARES ABSOLUTE VALUES RESIDUALS
AVR : AVERAGE ABSOLUTE VALUE RESIDUALS
N : AUMBER OF DATA POINTS
(4] : RATIO AVERAGE ESTIMATED VALUES/AVERAGE TRUE VALLES
(5] : ALL THE VALUES SHOW IN THIS TABLE FOR THE ONE-DIMENSION-1 INTERPOLATION PROCEDURE
CORRESPOND TO THOSE FOR CASE STUDIES I, II AND III.
304
Table 6.11. Estimated piezometric heads for Cases I, II, III, IV
and V, using the one-dimensional interpolation
procedure.
70 and 91 to 100 (see Figure 6.9). Figure 6.13 shows how close
305
the one-dimensional interpolation method "follows" the true
piezometric plane. Also evident from Fig. 6.13 is the oscillatory
parameters: the average of the piezometric head for all the nodes
Tables 6.10. and 6.12. show clearly that, on the whole, the
306
•
04
0
rn
.p
-1-)
z
-P o
Ca
rn
1:1
• 0
P
-
C.)
Ca
•ri (1)
• in
fat
cci
o
cr)
fr-4
sq
g
B MW *Pau10Mci lowligr3
sq
3 2 •
p0914 VILZW pa4Du.ing3
-
z
1,1111/1111111111111TIII1IIIIITITIIIIII11
• q sq q
0
;
- 8 • 8 2 --2
Ivlisiv“Dilivvi.vvvi.mivroviv,11,11111 J ill •
?
v. S
Y; ; SY ; ;
307
Table 6.12. Performance indices for the different head estimation
procedures, for examples A, B, C, D, E and F.
a) Performance indices.
PERFORMANCE INDEX
EXAMPLE PROCEDURE AVERAGE VARIANCE MAXIMUM VARIANCE
RESID. RESID.
INTERPOL. 0.981 0.751 0.825 0.901
A KRIGING -3.593 -35.666 -23.955 -999.999
SPLINES -999.999 -999.999 -999.999 -999.999
INTERPOL. 1.000 -5.504 0.833 0.966
B KRIGING -5.060 -268.385 -29.088 -999.999
SPLINES -844.286 -999.999 -999.999 -999.999
INTERPOL. 0.992 0.972 0.948 0.996
C KRIGING 0.922 0.955 0.628 0.645
SPLINES -295.459 -999.999 -999.999 -999.999
INTERPOL. 0.990 -1.321 0.421 -14.800
D KRIGING -0.340 -999.999 -408.497 -999.999
SPLINES 0.711 -691.198 -133.725 -999.999
INTERPOL. 0.944 0.912 0.155 0.000
E KRIGING -956.712 -999.999 -999.999--999.999
SPLINES -158.527 -999.999 -999.999 -999.999
INTERPOL. 0.986 0.939 0.408 0.945
F KRIGING -88.996 -7.101 -450.198 -999.999
SPLINES -999.999 -999.999 -999.999 -999.999
NOTE: A VALUE OF -999.999 INDICATES THAT THE INDEX
IS ACTUALLY LESS THAN OR EQUAL TO -999.999
b) Average indices.
AVERAGE INDEX
PROCEDURE AVERAGE VARIANCE MAX.RES VAR.RES
INTERPOL. 0.982 -0.542 0.598 -1.832
KRIGING -175.630 -385.032 -318.518 -833.225
SPLINES -549.593 -948.532 -855.620 -999.999
c) Frequency of improvement.
FREQUENCY OF IMPROVEMENT (%)
PROCEDURE AVERAGE VARIANCE MAX.RES VAR.RES
INTERPOL. 100.000 66.667 100.000 66.667
KRIGING 16.667 16.667 16.667 16.667
SPLINES 16.667 0.000 0.000 0.000
308
procedure; the worst case corresponds to example D. Table 6.12.
of 1.0.
Table 6.10. shows that when compared with Kriging and splines,
309
The results for the estimated heads using Kriging and the bi-
For the six networks specified in Section 6.10 and for the five
cases detailed in Section 6.11.2, the calibration exercise has
computed before the calibration program has been run, i.e. their
values have been determined using the initial estimates for the
flows determined with the calibrated model, rather than with the
raw model.
For each one of the six examples, and for each study case
(Tables D.1, D.2, D.3, D.4, D.5 and D.6). In order to quantify
310
the performance of the calibration algorithm, for each one of the
D.11 and D.12). Each one of the last 6 tables (with the
across the five study cases (mid part of the tables) and a
of the "performance index" for the five cases considered, and its
all the study cases (last rows of the Tables). Also, a "global
average" has been computed, across all the examples and study
311
0 el-I
• rl 0
O $4
O +5 0.1
•,-i 0 0
0 4-) E 0
.0 be 0
a
1 -4
os
r-I
(n
0
N
a)
•,-4
4-)
c.)
> c C4 a) •r 0
•ffi $4 c f:14
(a
- 0 to 0 ...4
O ..-, 4-) ,-i ICJ a)
c.,) 14
0-1 E- W I-I
0
I f,4 0 OA
cn 4-)
a)
(1:i
0
a)
, r $.4 >
cC) , t
r * .0 0
1 g
›-• •1 I •,-4
cr) :1 ,-4
(J $4 I-I
n I
(•) 43
cHEn
o -4.
(i)
cri
U
,.
0
o
P 0
4-,
Ca
•r
U)
Cd
01 Ti
a
S
0 a)
-• 0 .0
/
di
n ,-1
(D
./
vt
2
I,,''' 111111111111TM
et
; ; i it 6
so c•H 31-q•uxizaid PewigiD3
- g
- z
-S
; a a
lIlIllIlIlil
a
II i
LI, a 2
0 n.
sPo pH *t l i owoze tal Palem el.P0
_.
.
,..N
0
111.1110 11111101011110111/101.10110.11III-
3 2 ; ;^
- - -
sP
i 2
o• 3 1UX49
H
-
in
; ;
14 PelcuclUD3
ort
6
312
Table 6.13. Summary of the comparison for the calibrated
piezometric heads: average heads.
NOTE
tilt
313
Table 6.14. Summary of the comparison for the calibrated
piezometric heads: variance of the heads.
NOTE
Mt
A VALUE OF -999.999 INDICATES THAT THE INDEX IS ACTUALLY
LESS THAN OR EQUAL TO -999.999
314
Table 6.15. Summar y of the comparison for the calibrated
piezometric heads: maximum residual.
NOTE
Mt
A VALUE OF -999.999 INDICATES THAT THE INDEX IS ACTUALLY
LESS THAN OR EQUAL TO -999.999
315
Table 6.16. Summary of the comparison for the calibrated
piezometric heads: variance of the absolute value of
the residuals.
-16.390 1
I
INTERPOL. -2.033 -5.326 0.997 -91.730 -1.023 0.773
AVERAGE KRIGING -999.999 -999.999 0.563 -999.999 -999.999 -999.999 -833.239
SPLINES -999.999 -999.999 -999.999 -999.999 -999.999 -999.999 -999.949
NOTE
heads (see Table 6.12, for example), in the sense that the
method are the best, both in the average and in the maximum
example E (case V), the calibrated heads computed using the one-
dimensional interpolation scheme always improve the initial (raw
6.14., 6.15. and 6.16. show that the same holds for the head
well). We shall discuss this point later on. Table 6.13 also
in Table 6.13.).
well even when the average of the initial roughnesses has been
317
Table 6.14 shows that, on average, the calibration algorithm
318
does not produce any improvement at all.
319
I, II and III). The corresponding "global success index" when
Kriging was used is 16.67 % (5 cases out of 30) and 0.0 % for
splines.
Fig. 6.15. represents the pipe flows for example E, Case I, for
- 31
0Z
41.
In
.-1
cd
r1111111IIIIIIII.1/7111.1111WII
4
q 0
r. $1.
3 3 ; 3
I
(sil) 2.01.4 rancuquD9
1
I
E
_g
o q
4
I
I
tr
to
—6
o
4
—:
i
(s/I) "%o ll PaIDAIPO
321
Table 6.17. Summary of the comparison for the calibrated
flows: average flows.
NOTE
tits
A VALUE OF -999.999 INDICATES THAT THE INDEX IS ACTUALLY
LESS THAN OR EOUAL TO -999.999
322
Table 6.18. Summary of the comparison for the calibrated
flows: variance of the flows.
NOTE
323
Table 6.19. Summary of the comparison for the calibrated
flows: maximum residual.
PERFORMANCE INDEX:MAXIMUMRESIDUAL
NOTE
slat
A VALUE OF -999.999 INDICATES THAT THE INDEX IS ACTUALLY
LESS THAN OR HUAI TO -999.999
324
Table 6.20. Summary of the comparison for the calibrated
flows: variance of the absolute value of the
residuals.
NOTE
tttt
A VALUE OF -999.999 INDICATES THAT THE INDEX IS ACTUALLY
LESS THAN OR EQUAL TO -999.999
325
Contrary to what we have seen in the case of the estimated and
calibrated piezometric heads, the behaviour of the calibrated
6.17 shows that, despite the fact that the global performance
flow.
all the examples (last column) of Tables 6.17 to 6.20, the values
find out that in examples A (Cases II, III and V), B (Cases II,
III and IV) and C (Cases IV and V), the average flow, its
both cases.
326
corresponding Table D.16, from Appendix D. Basically, what is
happening in this particular example is that the initial flows
are quite close to the true values (this is a consequence of
assuming that we know the mean pipe roughnesses which, by the
way, is one of the reasons why we use the raw model flows as
estimates of the true flows); this can be seen from Table D.16,
where the true and average calibrated flows differ only by 5x10-6
(1/s). This implies that a calibrated flow which on average is
onl y 3.9x10 -5 (1/s) different with respect to the true value (as
in Case I) is worse than the initial estimate, producing a large
negative performance index; in practice, a difference of this
magnitude is irrelevant.
As can be seen from Table D.16 (for example D), all the
calibrated flows computed either with the one-dimensional
327
interpolation, Kriging or splines give flows which, on average,
and III, with respect to case I (when only head measurements were
performance index near the true values, Examples C and E did not
examples (except D), for the average flows, and some smaller
328
Thus, in summary, as far as the estimation of the true flows is
work.
6.22, 6.23 and 6.24 condense the results of the performance index
for the average C's, their variance, the maximum residual and the
(shown by the asterisks in Fig. 6-16) and with the true values
329
that we are dealing with a discrete variable. Fi g ures 6-16 b) and
c) correspond to the calibrated C's obtained using Kriging and
splines, respectively.
Figure 6-16 a) shows clearly that the C's produced using the
one-dimensional interpolation tend to "follow" the initial
estimates in some pipes. It is also clear from that Figure that
the average calibrated roughnesses are close to the true values,
though their dispersion (variance) is larger than that of the
initial estimates.
For the C's produced using Kriging and splines Fig. 6-16 b) and
c) and Table D-29 show clearly that they are, on average,
underestimating the true roughnesses. What seems to be even
worse with the roughnesses produced with Kriging and splines is
the fact that they are "jumping" from one border to the other of
the "allowable band" [determined by the maximum variation factor
x=0.10, see equation (4)], but they never get close to the true
values. This explains why, in general, the variances of the
absolute value of the residuals are always smaller (when using
Kriging and splines) than those produced using the one-
dimensional interpolation; roughl y speaking, the residuals
(absolute value) produced by using Kriging and splines are
consistently wrong, while those C's produced using the one-
330
co
to
(I)
0)
r--1
0 1-1
ed
4.)
0 C/1 *el
1-3 C
.•••••n
CD 1-4 E4
CO
0 t
0 n r *
›"‘ 1
E
.• oa
• 0
•
0
0 In 0 In 0 o In o
ID in in a PI PI
0 0
-co -co
• ••• •
•
...nn••••
.••nnn•
•0
n111
o 0
0
.s-i
0 t
-o .0
- E E
3
:co
- _o
o co b:
0
4-)
•el •
03
-0 0 _o
a
• el
(NI 0
• r•I
•
filTfl .1,8 0 ;III 1•11111“111111•IIIII 0
o 0 o o In 0 In 0 In
••••
to In * PI PI CNI C-4 ID in in •P •t• Pn
41,
0
S.0 swomiJA -uazoil pa l oicth oo puo lomui 5,0 SU.101. 111A- u o il pa l oig n i oo puo
331
Table 6.21. Summary of the comparison for the calibrated
C's: average C's.
1
V KRIGING -230.244 -19.121 0.585 -999.999 -999.999 -999.999 -541.463
SPLINES -773.546 -395.178 0.819 -999.999 -999.999 -999.999 -694.650
NOTE
Itti
A VALUE OF -999.999 INDICATES THAT THE INDEX IS ACTUALLY
LESS THAN OR EQUAL TO -999.999
332
Table 6.22. Summary of the comparison for the calibrated
C's: variance of the C's.
I
INTERPOL. -37.185 -13.043 -1.196 -185.404 -134.856 -63.950 -72.606
AVERAGE KRIGING -7.865 -12.819 0.668 -512.155 -374.333 -259.790 -194.382
SPLINES -19.279 -7.398 0.892 -512.637 -424.031 -422.013 -230.745
/
NOTE
Mt
A VALUE OF -999.999 INDICATES THAT THE INDEX IS ACTUALLY
LESS THAN OR EQUAL TO -999.999
333
Table 6.23. Summary of the comparison for the calibrated
C's: maximum residual.
I
INTERPOL. -7.806 -7.227 -1.438 -5.466 -4.673 -4.737 -5.224
V KRIGING -6.704 -7.806 -0.457 -4.737 -5.466 -5.466 -5.106
SPLINES -6.704 -6.704 0.800 -4.737 -5.466 -5.466 -4.713
NOTE
MI
A VALUE OF -999.999 INDICATES THAT THE INDEX IS ACTUALLY
LESS THAN OR EQUAL TO -999.999 .
334
Table 6.24. Summary of the comparison for the calibrated
C's: variance of the absolute value of the residuals.
I
INTERPOL. -268.537 -190.796 -14.355 -150.180 -200.832 -92.442 -152.857
III KRIGING -47.169 -42.248 -6.104 -17.321 -9.159 -22.785 -24.131
SPLINES -31.983 -34.296 0.931 -11.248 -17.220 -10.806 -17.437
NOTE
stiR
and 16.7 % when the bi-cubic splines have been used for that
one-dimensional interpolation.
336
overestimated by 20% with respect to the true values. This can be
seen from Tables D.25 to D.30, where the ratio between estimated
and true averages of the C's are shown in the last column of
those Tables; on averaging those values across the different
study cases (I to V), we obtain Table 6.25.
337
Table 6.25. Summary of the ratios between average calibrated
roughnesses and average true roughnesses.
338
To study the effect of the "maximum variation factor" on the
results of the calibration exercise for example C, we ran the
variation for C (i.e. x=0.30 in equation (4)] and the results can
be summarised as follows:
worse.
- The new calibrated flows are much worse than when x=0.10.
On the other hand, the factor "x" can be used to spot pipes or
339
program displays a message in that case), it means, unless there
calibration procedure.
the way in which the calibration exercise was planned, means that
340
6.11.5. Summary of the results.
estimating the unmeasured heads (see Tables 6.10 and 6.12) the
341
The main reasons for the success of the one-dimensional
interpolation scheme, seems to be attributable to the fact that
it is the only method which explicitly considers the initially
assumed shape of the piezometric plane. Kriging and splines-based
estimation algorithms seem to be smoothing the estimated
piezometric plane, whereas the interpolation method has the
ability to keep the main structure of the initial piezometric
plane (raw model) throughout the whole estimation procedure.
342
determine accurately at least the relative value of the
first run of the calibration program has been made; the new head
calibrated model, not with those of the raw model. The calibrated
of the network. Perhaps this can help to reduce both the variance
As expected from what has been noted in section 6.7, and from
(compare Figures 6.13 and 6.14). Indeed, the use of the one-
343
0) Calibrated pipe flows.
the heads, and this is not good enough for our purposes,
approach has emerged, which avoids the use of the raw model by
The basic idea is to relate the flow measurements with the co-
344
as defined in Chapter Two, section 2.3.7, Fig. 2.1), and to
comp ute the rest of the flows (i.e. the dependent ones) directly
from the flow measurements. However, this has yet to be tested
piezometric heads.
and flows provides all the true nodal heads and pipe flows, the
Additionally, and due to the fact that the head and flow
not reasonable for some pipes. For example, if for pipe joining
0 •
nodes 1 and "j", the head estimator Hi > Hj and the flow
345
(when available). Undoubtedly, this should improve the
performance of the calibration algorithm.
For the flow estimation, the raw network model, obtained with
initial estimates of the roughnesses, provides the estimates of
the flows for the unmeasured pipes. Because those estimates, when
merged with the flow measurements, are neither balanced nor
compatible with head estimates, an iterative algorithm
approximates those flow estimates to the true ones ( as shown in
Fig. 6.1.), also giving the estimates of the C's.
346
ii) Following on from i), the calibration algorithm produces
calibrated piezometric heads which follow closely the estimated
and splines, but not good enough, and an alternative new approach
for unmeasured pipe flow estimation has been suggested for future
implementation.
347
To summarise, we believe that the proposed calibration
algorithm has shown to be the adequate framework for solving the
parameters.
networks.
346
Some way of testing the quality of head and flow estimates,
measurement system has not been studied here, but the need for a
7.1. Introduction.
discuss the implications that this approach can have for extended
Coulbeck and Orr (1983), Rao and Bree (1977)]. The main objective
350
has been the study of the performance of the gradient method in
such applications, particularly with some of the implementation
features discussed in Chapter Five.
not been covered in the present work. The program computes the
flow) only one flow is needed, while for pumped sources an on/off
diagram with the maximum flow is enough.
351
and the outflow ( Qi ), i.e.:
d hi QINi - Qi
(2)
dt (d Voli/d hi)
Voli = ai hi 2 + bi hi + ci (3)
where ai, bi and ci : are constants, dependent on the geometrical
d Volt
= 2 ai hi + bi (4)
d hi
d hi QINi - Qi
= fi(t, hi) (5)
dt 2ai hi + bi
Voli = bi h i + c i(6)
and
d Voli = bi
(7)
d hi
352
d hi QINi - Qi
- fi(t, hi) (8)
dt bi
be approximated as:
h 1 (k + 1) _ hi( k)
fi(t, h i ( k )) (9)
St
functions to be obtained.
353
the error accumulated during a 24 hours simulation has been
computed as 5.8 m 3 , in a total inflow of 12,960 m 3 , for a
simulation with a time step of one hour. The error was defined as
the difference between the total inflow, on the one hand, and the
the other hand. This shows that the error due to the simple
\ Read data
Initialise
time step
k = 0
yes
L'="1
Print results
( Stop)
355
7.3. Extending the gradient method for Pressure-dependent nodal
demands.
well known fact that the actual consumptions in the network are
leakages.
where the pressures at the nodes change over a wide range. The
match the performance of the network model with that of the real
very low flows (i.e.: very low and very high pressures), to cover
356
Lam and Walla (1972 b) proposed a general model to include the
influencehe pressures on the nodal demands:
qi = ai + bi H i c(11)
qi = ki Pi x(12)
where
emitter operating pressure.
(emitters).
357
recurrent algorithm is used, based on the repetitive solution of
a standard network analysis algorithm, verifying if the head at
proposed by Lam and Walla (1972 b), allowing the system itself to
find the equilibrium point between the actual nodal outflow and
pressure.
358
produces an increase in the size of the leaking cracks or holes,
reason deals with a more conceptual topic, connected with the way
where
pressure.
pressures:
359
Theoretical
w relationship
0
<
<
w Actual
_J
/ / relationship
/
/
/
/
I i r 1 1
0 20 40 60 80 100
PRESSURE (m)
q A
Demand
(I/s)
Actual demand
Linear model
—4 \----
Constant demand
nn•
p : Pressure (M)
P s P ma.
At
Service pressure
360
P = H - Z (15)
where:
H : piezometric head (NN-NS)xl column vector.
g = % - A22 H (17)
with
state flow in the network [see Chapter Three, equation (3)], for
[ A ll : A l2 I [ I
=
[ -A10 11 0 I
A 21 : 0 H a
and introducing the pressure-dependent nodal demands from (17),
we get:
361
{ A 11 A l2 Q -A10 HO
(19)
A 21 : A 22 I I Hgo
N A 11 * 1
Al2 dE
: * dQ = (20)
A 21 : A 22 I d HI d g
[
u ll Ql1n1-1 51/Q1
Q21n2-1
*21 132/Q2
A 11 * = =A11
aNFIQNprINP-1
1-6NP/QNP
d Q(i) (21)
E = All") Al2 H(i) +A 10 HO
and
d g = A 21 Q (i) 4- A22 H (i) - go (22)
362
The solution of (20), for the flow and head increments, can be
obtained as:
-1
[dCil N All * 1 A 12
= : * (23)
d H A21 1 A22
I [ : : 1
N A 11 * 1 A 12 [ B ll : B 12 1
1
1 = : (24)
[ A21 : A22 B 21 1 B22
On defining:
G = N A11* (25)
d Q =B11 dE + B 12 dg
(27)
d H = B21 d+ B 22 dg
1
which, when introducing equations (21) and (22), becomes:
d Q = B il( A llO i4A l211 (i) + A lollo] + 13 12[ A 210 i)+A22H(i)-go] (28)
and
363
d H = B 21 [A 11 Q( i)+A l2H" )+A 100] B 22 (A21Q (i) +A20 (i) -20] (29)
gi(=
i +1) (I- G-
A ll ) Q(i) - G-1 ( Al2 H(1+1) + A 10 Ho) 1 (33)
hand side the previous nodal demands g are now replaced by go, as
defined in equation (18) of the present chapter. Of course, when
364
At the present time the computer implementation of this
365
CHAPTER EIGHT
8.1. Summary.
The present work has been mainly concerned with the development
366
Kriging, another based on bi-cubic splines and a third based on
8.2. Conclusions.
367
comparison with most of the existing methods. A comparison with
the near future. Other advantages of the algorithm are such that
368
maximum of 1,000 nodes, though is four times faster than the
formulation.
369
deterministic interpolation method. The results of such a
pipe flows does not give results which are as good as the head
obtained are far from the true ones, indicating that further work
solution.
370
well-known fact that many engineers are not very keen on highly
sophisticated mathematical developments and, unfortunately, this
may be the case for the gradient method and its related
users.
371
6.2.2 Calibration.
has not been developed in the present work, is to use the graph-
head estimates after the calibration algorithm has been rwn was
the future.
372
d) One of the assumptions made in the proposed calibration
method is that the nodal consumptions are known. This hypothesis
the estimates. In this context, only Kriging and splines are able
devised. The same should be required for the flow estimates, thus
resistance parameters.
estimates.
and flows) has not been explicitly tackled here. This subject is
373
Thus, it seems that something should be done in the future in
this respect. Kriging provides a neat way of determining the
374
REFERENCES
375
BHAVE P. R.: (1981), "Node flow analysis of water distribution
systems", Transportation Engineering Journal, ASCE, Vol. 107,
No. TE4, pp. 457-467.
CHIN K. K., GAY R. K. L., CHINA S. H., CHAN C.H. AND YO S.Y.:
(1978), "Solution of water networks by sparse matrix method",
International Journal for Numerical Methods in Engineering,
Vol. 12, pp. 1261-1277.
376
CLEMENTS K.A., KRUMPHOLZ G.R. AND DAVIS P.W.: (1983), "Power
system state estimation with measurement deficiency: an
observability/measurement placement algorithm", IEEE
Transactions on Power Apparatus and Systems, Vol. PAS-102,No.
7, July, 1983.
377
March, 1986, Institute of Physics Conference Series Number 81,
Bristol, IOP Publishing Ltd., pp. 103-125. Also NFL Report
DITC 81/86.
COX M. G. : (1987a), Data approximation by splines in one and
two independent variables", in Iserles A. and Powell M.J.D.
(Editors): The state of the art in Numerical Analysis",
Oxford University Press, 1987. Also NFL Report DITC 77/86,
October 1986.
378
DONACHIE R. P.: (1974), "Digital program for water network
analysis", Journal of the Hydraulics Division, ASCE, Vol. 100,
No. HY3.
379
pp. 1145-1151.
380
GEORGE A. AND LIU J. W-H.: (1981), "Computer solution of large
sparse positive definite systems", Prentice Hall.
381
the geostatistical approach to the inverse problem in two-
dimensional groundwater modeling ", Water Resources Research,
Vol. 20, No. 7, pp. 1003-1020, July 1984.
HOEKSEMA R. J. AND KITANIDIS P. K. : (1985), "Analysis of the
spatial structure of properties of selected aquifers ", Water
Resources Research, Vol. 21, No. 4, pp. 563-572, April 1985.
382
KARMELI D., PERI G. AND TODES M.: (1985), "Irrigation systems
design and operation", Oxford University Press, Cape Town.
383
LEWIS J. G. : (1977), "Algorithms for sparse matrix eigenvalue
problems", Computer Science Department, Stanford University
Report STAN-CS-77-595.
LO K.L., ONG P.S., McCOLL R.D., MOFFATT A.M. AND SULLEY J.L.:
(1983 a), "Development of a static state estimator. Part I:
Estimation and bad data suppression", IEEE Transactions on
Power Apparatus and Systems, Vol. PAS-102, No. 8,pp. 2486-
2491, August, 1983.
LO K.L., ONG P.S., McCOLL R.D., MOFFATT A.M. AND SULLEY J.L.:
(1983 b), "Development of a static state estimator. Part II:
Bad data replacement and generation of pseudomeasurements",
IEEE Transactions on Power Apparatus and Systems, Vol. PAS-
102, No. 8, pp. 2492-2500, August, 1983.
384
MONTICELLI A. AND WU F.F.: (1985), "Network observability:
identification of observable islands and measurement
placement" IEEE transactions on Power Apparatus and Systems,
Vol. PAS-104, No. 5, May 1985.
385
ORMSBEE L.E. AND WOOD D.J.: (1986), "Explicit pipe network
calibration", Journal of Water Resources Planning and
MAnagement, ASCE, Vol.112, No. 2, April 1986. pp 166-182.
386
RATCLIFFE B.: (1986), The performance and selection of pressure
reducing valves", Technical Report TR 238, WRC Engineering
Centre.
387
SHERMAN A. H.: (1975), On the efficient solution of sparse
linear and non-linear equations", Report 46, Ph. D. Thesis,
Department of Computer Science, Yale University.
SONG C.C.S. AND YANG C.T. (1980), "Minimum stream power: theory",
Journal of the Hydraulics DiViSi011, ASCE, Vol. 106, No. HY9,
pp. 1477-1487.
SONG C.C.S. AND YANG C.T. (1982), "Minimum energy and energy
dissipation rate", Journal of the Hydraulics Division, ASCE,
Vol. 108, No. HY5, pp. 690-706.
388
reti idrauliche", Bolletino degli Inaegneri della Toscana, No.
11 , pp. 11-14 (in italian).
389
their reliability", Research Re p ort No. 127, Water Resources
Research Institute, University of Kentucky.
WOOD D. J.: (1981b), discussion to Isaacs and Mills (1980):
"Linear theory methods for pipe network analysis", Journal of
iy drauLia r
th el r l'aiom ASCE, March 1981, pp. 384-385.
WOOD D. J. AND CHARLES C. O. A.: (1972), "Hydraulic network
analysis using linear theory", Journal of the Hydraulics
Division, ASCE, Vol. 98, No. HY7, pp. 1157-1170.