The Colebrook-White formula for pipe networksG. Keady,Mathematics Department,University of Western Australia, Nedlands, 6907.January 1995SUMMARY Flow resistance laws, as used for example in water-supply pipe networks, areformulae relating the volume ow rate q along a pipe to the pressure-head di�erence t between itsends, q = (t). is monotonic. The simple Hazen-Williams power \law" is often used: it hasbeen claimed that the more complicated Colebrook-White law (CW) better represents aspects ofthe experimental data.Equilibrium ows in pipe networks are described by convex optimization problems, where theobjective function in a primal form involves the inde�nite integral of . This is analytically tractablefor both Hazen-Williams and (CW). There is also a dual of this optimization problem, the startingpoint for which is �, the inverse to . It turns out that for (CW) � can be expressed in terms ofa known transcendental function, the Lambert W -function. We explain how these facts might beuseful in computer codes.Some simple numerical work, and some phenomena of independent interest are brie y reported.See Appendices.NOMENCLATURE The following symbols are used in this article and in [18].A arc set for a nework; U dual objective function;bi consumption at node i; V primal objective function;c1, c2, c3 constants in CW ; W Lambert W -function;Da diameter of pipe a; � viscosity of waterg acceleration due to gravity; � inverse to ;Ka pipe roughness for a; � integral of �;La length of pipe a; ow function;N node set for a nework; integral of ;pi pressure head at node i;P power consumed;qa ow along pipe a; ka (separable) conductivity factor;ta pressure head di�erence along a; � (separable) ow function.These are not used in [18]. 1

E node-arc incidence matrix;f friction factor;Fr Froude number;Re Reynolds number;R integral of �; ra (separable) resistivity factor;S integral of �; � (separable) resistance function;1 INTRODUCTIONThere is a substantial, and growing, body of knowledge concerning general network ow problems.See, e.g. [27]. In network ow problems, there is a potential function de�ned over the nodes and a ow de�ned over the arcs. Furthermore the ow on any arc a is a given increasing function of thedi�erence in potential between the arc's initial and terminal node. In the case of equilibrium owin pipe networks, the potential is the head. Two of the commonly used empirical rules giving the ow as a function of head-di�erence are Colebrook-White and Hazen-Williams. Hazen-Williamsnetworks are special cases of what we call Type (H) networks. `Power-law nonlinearity' ow-lawsare where the ow qa on an arc a joining node i to node k satis�es, for some constant ka,qa = ka�(p(i) � p(k); sa) where �(t; s) = tjtjs�1; s > 0:A network is said to be Type (H) when there exists s > 0 such that for all arcs a, the ow-law isas above, with the same s on every arc. Hazen-Williams' means Type (H) with s = 27=50. (The`rough-pipe limit' of Colebrook-White (1) has c3 = 0 and is Type (H) with s = 1=2.)The �rst parts of this paper, in x1.1, are concerned with elementary mathematical facts about theColebrook-White rule which are of use, for example, in connection with the convex optimizationformulations of the network- ow problem. Because we feel that the convex optimization approachto the problem may yet further in uence methods of computing pipe network ows, we also collect,in x3, some of the results in order to publicise them to practitioners involved with pipe network ows.A later part of this paper, in x4, is concerned with an example of \Braess's paradox". In thepipe-network problem with Hazen-Williams' nonlinearities, that is the same s > 0 on each arc,given the consumptions and supplies, the power usage is a decreasing function of the diametersDa. (There is no Braess' paradox with Hazen-Williams.) In certain sorts of networks, this neednot be the case when the Colebrook-White rule is used. An example, involving smooth pipes andchosen for mathematical simplicity rather than engineering realism, is given in this paper. It maybe that here the Braess paradox is just a mathematical curiosity. In particular we expect that thephenomenon is likely to be rarely observed with water-supply networks in practice (except perhapsif there are a variety of di�erent sorts of nonlinear elements, di�erent from pipes, e.g. valves andpumps). However we feel that it may be worth recording if only because it might be noticed in thehigher precision obtainable with numerical computation, and it doesn't mean that the programsare necessarily wrong.1.1 A single pipeBefore treating networks, we consider the case of a single pipe. With the nomenclature as above,de�ne c1 = �ln(10)rgD52L ; c2 = 10K37D ; and c3 = �s L2gD3 :Let q denote the ow down the pipe, and t the tension, or head di�erence between the ends of thepipe. The Colebrook-White rule for the turbulent ow down a pipe states that an approximationto the ow rate q along the pipe is,q = CW (t) = �c1pt ln(c2 + c3pt ) = �c1r tgL ln(c2 + c3rgLt ); (1)2

when the tension t is su�ciently large positive that the argument of the ln is less than one, i.e.tgL � t0gL = c23(1� c2)2 :(ln denotes the natural logarithm. The coe�cients c1 = c1pgL, c2 and c3 = c3=pgL are functionsof D, K, g, � only and are independent of L.) The simplest way of extending this so that it isde�ned for all positive values of t is to de�ne CW (t) = 0 for 0 < t � t0.Later in this paper we will consider more general ow laws q = (t). As in the Colebrook-Whitecase, we consider only which are nondecreasing functions of t and which have (0) = 0. Theextension to negative t is always done so that is an odd function of t. However, for the remainderof this section we need only consider CW .It is sometimes useful to express the tension in terms of the ow, t = �(q). This is done as follows.THEOREM 1.1 CW is nondecreasing on (0;1), increasing on (t0;1), CW ! 1 as t ! 1. CW is concave on (t0;1) and log-concave on (0;1).For t > t0 the inverse function �CW to function CW de�ned by equation (1) is�CW (q) = � c3c1c3q W ( q exp( c2qc1c3 )c1c3 ) � c2 �2; for q > 0: (2)Here W denotes Lambert's W -function, de�ned by W (x) exp(W (x)) = x (and we need only x > 0,W (x) > 0).�CW is increasing on (0;1), �CW (q)!1 as q!1.This was discovered by running the single line of Maple:solve(q=-c1*sqrt(t)*log(c2+c3/sqrt(t)),t);and then checked by hand-calculation. �CW has, of course, been approximated numerically byengineers for some time. There are many other uses of the Lambert W -function, and it seemspossible that exchanges concerning the numerical methods used in engineering to calculate �CW ,and the numerical methods used in these other uses of W could be mutually bene�cial. (Byway of caution, we comment that di�erent applications have di�erent needs. In the pipe-networkapplication, speed is more likely to be important than great accuracy.)Proof. The items in the �rst sentence are elementary, and the proof is omitted. We verify (2). Itrearranges as follows:With � = q exp( c2qc1c3 )c1c3 ; W (�) = qc1c3 �c2 + c3p�CW � = �:The last equation in the above de�nes �. However, quite generally (not merely for the � and � wehave de�ned), if W (�) = � then � = � exp(�). Using this, and writingt = �CW ; q13 = qc1c3 ; �13 = �q13 = c2 + c3pt ;we �nd that q13 exp(c2q13) = � = � exp(�) = � exp(c2q13) exp(c3q13pt );so 1 = �13 exp(c3q13pt ) = �13 exp( qc1pt ):Finally, taking logs of the last equation, we �nd that q is expressed by the formula at the right of(1).For properties and uses of Lambert's W -function, see [9]. In practical use for machine-precision oat computation of the W -function, the routine in [13], appears to be adequate.3

1.2 The friction factorNo items in this subsection are needed in other parts of the paper. Useful nondimensional param-eters are the Froude number Fr and the Reynolds number Re whereFr2 = Q2gD5 ; Re = 4Q��D:Engineers sometimes are concerned with the friction factor f satisfyingtgL = f Fr22�2 :See [4]. The Colebrook-White rule gives the following equation:1pf = �2 ln10�10K37D + 251100Repf �:Either from the items earlier in this paper, or directly, it can be shown that this equation can besolved explicitly for f in terms of Lambert's W -function. This is likely to be an observation oftheoretical rather than practical signi�cance, as numerical approximations such as [8] are probablysatisfactory in those circumstances when the quantity f is needed. (The friction factor is not neededin the optimization formulations of the pipe network problem given in this paper.)2 NETWORKS2.1 NotationWe begin this section giving the notation of [27], which is now standard, used in network ow prob-lems. This is used in stating the `network equilibrium problem' and in some related optimizationproblems.Our notation follows that of Bertsekas et al. [1] and of Rockafellar [27]. A network G = (N;A), orthe directed graph associated with the network, consists of two �nite sets A and N and a functionthat assigns to each a in A an ordered pair of nodes (i; k) such that i 6= k. The elements of Acorrespond to the pipes and are called arcs (or edges or { in [4] { elements): the elements of N arecalled nodes (or vertices). We may label the arcs with the �rst jAj positive integers, i.e. startingfrom 1, and identify A with this set. We label the nodes with the jN j nonnegative integers from 0to jN j � 1. The correspondence of arc j with its nodes is written j � (i; k). We say that i is theinitial or start node of arc j and k is the terminal or end node of arc j. The de�nitions are fussy inorder to cover the cases which occur in engineering, e.g. parallel pipes such as in [4] p70, Example4.1.Let E be a node-arc incidence matrix for G, that is assign a column vector to each arc with 0severywhere except for a 1 for one node of the arc and �1 for the other. More precisely, the entrieseij of E are given by eij = (+1 if i is the initial node of arc j,�1 if i is the terminal node of arc j,0 otherwise.We always assume that G is connected, so that jAj � jN j � 1, and the rank of E is jN j � 1.For each node i, let p(i) denote the head (or, in other contexts, potential or voltage or time) atnode i. For arc a, let qa denote the ow on arc a, from start to end. For node i, let b(i) be theconsumption (to use the same word as in [4]) (or de�cit or input, or current supplied) at node i.For each arc a, suppose there is given a real-valued function a of a single real variable. The function a is the ow function and it will relate the head di�erences and ows by (4). For the purposes ofthis paper it su�ces to think of being given the diameters Da, lengths La and roughnesses Ka andusing the Colebrook-White formula (1). 4

DEFINITION. We say the network ow functions satisfy Assumption A( 0) if, for each a 2 A, ais Holder continuous on IR for some positive exponent, a(0) = 0, a(�t) = � a(t) for all t, a isnondecreasing on (�1;1), and either(i) a � 0 or(ii) a(t)!1 as t!1, t0 de�ned as the in�mum of ft > 0 j (t) > 0g is bounded and either,(iia) t0 > 0 and a 2 C1([t0;1)), and 0a(t) > 0 for all t > t0 or(iib) t0 = 0 and a 2 C1(0;1), with 0a(t) > 0 for all t > 0.DEFINITION. We say the network satis�es Assumption A+, if all functions satisfy assumptionA( 0) and if, with A+ = fa 2 A j a 6� 0g, then (N;A+) is a connected graph.In physical terms this corresponds to the network formed from pipes having nonzero diameter beingconnected. As all of our results depend on G being connected and assumption A+ being satis�ed,these are presumed to be satis�ed throughout the paper.In this paragraph we state a problem, the network equilibrium problem, { less general than in [4]Chapter 3 { but su�cient to illustrate the points we wish to make. The standard notation is [27,x8H]. Here is the de�nition. Given consumptions b satisfying the consistency condition,Xi2N b(i) = 0; (3)given a network satisfying Assumption A+, �nd a head vector p in IRjNj such that,qa = a(p(i) � p(k)); 8a � (i; k) 2 A; (4)and Eq = b: (5)We remark that the above only determines p up to a p + ce, c any real number, where e is thevector all of whose entries are one. When we want to make p unique we suppose its 0-th entry is 0.A function (t; c1; c2; : : : ; cn) is said to be separable if (t; c1; c2; : : : ; cn) = k(c1; c2; : : : ; cn)�(t):When � is (strictly) monotonic we describe the as separable and strict, and then always normaliseso that �(1) = 1. If there is just one signi�cant parameter c1, call this single parameter, and describek as the conductivity parameter. In an earlier paper [5] we treated networks where each arc had asingle-parameter separable and strict . The single-parameter separable and strict set-up includesHazen-Williams as a special case, where the ow laws on each arc are power laws. A `power-lawnonlinearity' is special single-parameter separable and strict ow-law where the ow qa on an arca � (i; k) satis�es, qa = ka�(p(i) � p(k); sa) where �(t; s) = tjtjs�1; s > 0: (6)A network is said to be Type (H) when there exists s > 0 such that for all a 2 A, �(t) = �(t; s).Hazen-Williams' means Type (H) with s = 27=50. The `rough-pipe limit' of Colebrook-White hasc3 = 0 and is Type (H) with s = 1=2.A function �(q; c1; c2; : : : ; cn) is said to be separable if�(q; c1; c2; : : : ; cn) = r(c1; c2; : : : ; cn)�(q):When � is (strictly) monotonic we describe the � as separable and strict, and then always normaliseso that �(1) = 1. If there is just one signi�cant parameter c1, call this single parameter, and describer as the resistivity parameter. In [6] we treated networks where each arc had a single-parameterseparable and strict �. As before, the single-parameter separable and strict � set-up includesHazen-Williams as a special case, where the ow laws on each arc are power laws. A similar5

de�nition of Type (H) can be given. The results of [6] apply when one takes pipe-lengths La of theColebrook-White formula �a as the signi�cant parameters, the others Ka, Da being held �xed.Some results in section 4 4.3 can be phrased in terms of `power-loss'. The power-loss P in thenetwork is de�ned by P = X(i;k)�a2A(p(i) � p(k))qa: (7)Where there are arcs in parallel joining i to k, the summation is over all of them. We always havei < k. Both of these conventions will be used elsewhere in this paper. Using our hypothesis thatall the a are odd functions, the power-loss P is nonnegative and, by [27, x1I],P = bTp: (8)In a two-terminal network, where b(i) = 0 except for the two nodes, one of which we always labeli = 0 and the other nt � jN j � 1, we have P = b(0)(p(0) � p(nt)). From this, in the language ofDC electric circuits, with b(0) > 0 the applied voltage di�erence p(0) � p(nt) is a nonincreasingfunction of any conductivity factor ka precisely when P is.[27, x8B] gives the following de�nition.DEFINITION. The set of pairs (b(0); p(0)� p(nt)) when p solves the equilibrium problem is calledthe characteristic curve of the two-terminal network with terminals 0 and nt.A small, but nontrivial, example of a two-terminal network is the Wheatstone graph, the twoterminals being 0 and 3.DEFINITION. The Wheatstone bridge graph is the graph G = (N;A) withN = f0; 1; 2; 3g A = f(0; 1); (0; 2); (1; 2); (1; 3); (2; 3)g:Four of the arcs form a quadrilateral, the other forms a diagonal. See Figure 1 on page 6.m0@@

@@

@@

@@

@@R

��

��

��

��

���

m3��

��

��

��

���

@@

@@

@@

@@

@@R

m1?m2Figure 1: The Wheatstone bridge network: in ow at node 0, out ow at node 3.Some curious properties of the behaviour or power-consumed P , when pipe parameters are variedwith b held �xed are discussed in section 4. 6

2.2 Nonuniqueness for Colebrook-White networksThe fact that �a is not strictly increasing gives us possible nonuniqueness of solutions. In particularconsider a network consisting of a single pipe, with prescribed input ow C = 0 at one end and thesame owing out at the other end. Then the solution set for t for this set-up is �t0 � t � t0.For a slightly less trivial case, consider the Wheatstone bridge network with in ow at node 0 andout ow at node 3. Suppose a small nonzero ow of C ows from node 0 to 1 to 3. Suppose thatthe other arcs (0; 2), (1; 2) and (2; 3) all have large t0;a and as a consequence do not have any ow.Then p(2) can take any value in a certain interval, yet be consistent with zero ows on the threearcs we have mentioned.3 THE CONVEX OPTIMIZATION PROBLEMFrequently the problem of �nding the ows in a pipe network is approached as a question of �ndingthe solution of the system of nonlinear equations f0(p) = b given by the network equilibriumproblem, equations (4), (5). Newton iteration is frequently used to solve this. See, for example, [4].Newton iteration, however, will only converge if the starting guess is close to the solution, and inpractice this is often not easy to satisfy. (In [4], it is also used that the jacobian of the nonlinearfunction f is positive de�nite.)There is a solution to this last practical di�culty. It turns out that the system of equationsf0(p) � b = 0 is such that f0(p)� b = gradp�V0(p)� pT b� = gradpVb(p);for some function V0. For reasonable ow laws �a and a are all nondecreasing: this has theconsequence that V0 and Vb are convex functions of p. Approaching the problem as an optimizationproblem, minimising Vb, has practical, numerical advantages in suggesting methods which convergefrom a wider set of initial guesses. Simply using the convexity of Vb we have that, while a full Newtonstep may give divergence, the Newton direction is a descent direction for Vb so an appropriate partialNewton step can be derived. (See, for example, [24].)Results concerning the pipe-network problem treated as an optimization problem are given in Hall[14], Collins et al. [11], Dembo et al. [12], [27], [28], [29].One observation about Colebrook-White facilitates the application of the convex optimizationmethod. The integral of is as follows. For t > t0,Z t0 CW (t)dt = CW (t)� CW (t0);where CW (t)c1 = �2c333c32 ln�c2pt+ c3� � 23 tpt ln�c2 + c3pt�+ 2c233c22pt � c33c2 t: (9)In connection with the dual form of the optimization problem, as described in subsection 3.2, wenote that the inde�nite integral of �CW can also be found. This is because an integration by partsgives, Z �CW (q)dq = Z t 0CW (t)dt = t CW (t) � Z CW (t)dt;or equivalently Z �CW (q)dq = tq � Z CW (t)dt:7

3.1 Primal form for general aIn Section 2.1 we de�ned the network equilibrium problem. [27, x8H] is the standard reference thatthis problem is equivalent to two others, namely the optimal di�erential problem, Problem (P) givenin equation (12) in the next paragraph, and a dual. The dual problem is the optimal distributionproblem, Problem (D) given in Section 3.2. The results in this section are not new, and for reasonsof space the proofs are suppressed. For detailed proofs and references see [16], [27].The optimal di�erential problem de�ned in [27, x8G] is essentially the following. De�nea(t) = Z t0 a(t)dt:De�ne V0 (summing over all arcs, pipes, a in A) then Vb by,V0 = X(i;k)�a2Aa(p(i) � p(k)); (10)Vb = V0 � jNj�1Xi=1 b(i)p(i): (11)De�ne X to be the set of vectors with coordinate indexing starting at 1 in IRjNj�1, augmentedwith a zero-th component which is zero. Problem (P) is to �nd p� minimizing Vb over all X, i.e.satisfying Vb(p�) = minp2X Vb(p): (12)THEOREM 3.1 With Assumptions A( ) on the a(ia) V0 is convex on X;(ib) V0(p)= k p k !1 as k p k!1 in any norm on X;(ii) solutions to Problem (P) exist and form a convex subset of X.(iii) p� solves the network equilibrium problem if and only if it solves Problem (P).Computations using Matlab's Optimization Toolbox, with formulation (P) and (9) are given in anAppendix. A favoured test example there is the Wheatstone bridge with ow in at node 0 and outat node 3, b(1) = 0 = b(2). (The results in this Wheatstone bridge case have been checked usingother approaches based on equation (2).)We will call items (iii)-(iv) Du�n's Existence Theorem.3.1.1 V (x ^ y) + V (x _ y) � V (x) + V (y)An additional property of V0, but which is not used in this paper, is the following. For a satisfyingAssumption A( ) in the pipe-network problem, we have the inequalityV0(p0 ^ p1) + V0(p0 _ p1) � V0(p0) + V0(p1) 8p0; p1:Here p0 ^ p1 denotes the minimum of the two vectors p0, p1 and p0 _ p1 denotes their maximum.3.2 Duality for general a, �aFor a satisfying Assumption A( 0), a has an interval-valued inverse. Call the inverse �a,�a( a(t)) = t; 8t 2 IR:Then �a is non-decreasing and 0 2 �a(0), and �a can be identi�ed with a strictly increasing C1function on (0;1). De�ne � to be the integral of � from 0, R(0) = 0. For q 2 IRjAj de�neU (q) = Xa2A�a(qa):8

The problem of minimising U (q) over q satisfyingEq = b;is Problem (D). Recall that eT b = 0, where e denotes the vector all of whose entries are 1. Problem(D) is a convex separable programming problem, which [27, x8D] calls the optimal distributionproblem.Both pairs of functions (Ra; Sa) and (�a;a) are Fenchel-conjugate convex function-pairs on IR.Collins et al. [11] establish the following Theorem. See also [27, x8L]. The notation involving the : IRjAj ! IRjAj, i.e. acting on vectors, is that (q) is the vector with components a(qa).THEOREM 3.2 Let G = (N;A) be a connected network.If p� solves the primal problem (P), de�ned above, then q = (ET p�) solves the dual problem (D).If q� solves the dual problem (D), and q� = (ET p) for some p with p(0) = 0, then p solves theprimal problem (P).At the solutions U (q�) = �Vb(p�).The constraints Eq = b restrict q to lie on a certain hyperplane.3.3 Mixed problemsThere are a greater variety of problems of interest to designers of networks than those describedearlier in this section. In [4] the code treats mixed problems, where, at the nodes other than node0 one is given one of p(i) or b(i), and is required to solve for the other. Presumably optimizationmethods remain equally appropriate. However, mixed problems are not considered in this report.4 BRAESS BEHAVIOURBraess [3] gave an example of a network, in the setting of tra�c ow, in which an extra arc - road -was added and the total travel time for any road user was increased. Braess's example has just oneorigin-destination pair so that it is, like the physical network ows in this paper, a single-commodity ow. (See Rockafellar [27, x1K] for de�nitions.)Nonlinear networks occur in many contexts: see Dembo et al. [12], Rockafellar [27].In the context of ows in physical networks, Braess's paradox can be restated informally as follows.BRAESS'S PARADOX. The power consumed in a nonlinear network can increase if an arc's con-ductivity is increased with consumptions held constant. In this section, by an arc's conductivityincreasing, we mean a pipe's diameter increasing.Theorem 4.1 states that Braess's paradox cannot occur in a two-terminal series-parallel network.One of the examples of physical networks is the ow of water in a pipe network. When the commonlyused Hazen-Williams' rule (with the same exponent on each arc, e.g. Brebbia and Ferrante [4]) isassumed, Theorem 4.3 shows that Braess's paradox cannot occur. (The result for the special cases = 1 has been known for several decades.)Whether a network is type (H), and whether a network can exhibit Braess's paradox are related:see the main theorems, Theorem 4.3 and Theorem 2 of [5], and also those of [6]. In a recent paperon Braess's paradox, Cohen and Horowitz [10] write: `The task remains of specifying the generalconditions under which such paradoxes can occur, for general network topologies and broad classesof components ...'. Our theorems in [5, 6] summarise some progress with this task: the nonseparablefunctions, such as diameter dependence in Colebrook-White considered in this paper, mostly exposethat more work remains to be done.The next two theorems from earlier papers, give some indications towards our proof of [5] Theorem2. 9

DEFINITION. Let n0 and nf be given nodes of a network G. The network G is said to be series-parallel with respect to n0 and nf if through each arc of G there is at leat one path from n0 to nfnot touching any node twice, and no two of these paths pass through any arc in opposite directions.An equivalent inductive de�nition is as follows. The one-arc graph G0 with A0 = (n0; nf ) is de�nedto be series-parallel. A network is series-parallel with respect to n0 and nf if it is either (i) aconnection of a series-parallel network with respect to n0 and ni in series with a second series-parallel network with respect to ni and nf , or (ii) a parallel connection of two networks bothseries-parallel with respect to n0 and nf .(The Wheatstone bridge graph not series-parallel with respect to f0; 3g but is series-parallel withrespect to f1; 2g.) A two-terminal network is said to be series-parallel if it is series-parallel withrespect to the two terminal nodes.THEOREM 4.1 For series-parallel two-terminal networks with each a single-parameter, sepa-rable and strict, the power-loss P decreases when any conductivity factor ka increases.This is Theorem 11 of [5]. See also Keady [16, Part I, Appendix A] for further discussion. Charac-terisations of series-parallel networks are also known.THEOREM 4.2 A two-terminal network is series-parallel if and only if there is no embeddednetwork having the Wheatstone bridge con�guration.This is Theorem 12 of [5].4.1 Main theorems for separable networksTHEOREM 4.3 Suppose (N;A) is a type (H) network and that b satis�es (3). Then, for any a,the power-loss P is a nonincreasing function of any arc's conductivity factor ka, or equivalently isa nondecreasing function of any arc's resistivity factor ra.The proof is given in [5] Section 5 together with a much stronger result with the same hypotheses.Before we state a converse, in Theorem 2 of [5], we give an additional example of Braess's paradox,this time involving only power-law nonlinearities on the arcs, which shows that it is necessary inTheorem 4.3 that the powers be the same on all the arcs.EXAMPLE of Braess's paradox. Consider the Wheatstone bridge graph, with the diagonal arc(1; 2) having variable-conductivity. The ow functions ki;j�i;j on the arcs (i; j) arek01�01(x) = k23�23(x) = x; k02�02(x) = k13�13(x) = xjxj�1=2;and the variable conductivity arc has k12�12(x) = kx:Choose b such that b(1) = 0 = b(2), so that the network is a two-terminal one.The result is that the power-loss P = b(3)(p(3)�p(0)) can either decrease or increase as k increases,depending on the value of b(3) = �b(0) and the value of k at the start.There is a lot of symmetry in the network and the ows, heads and power-loss can be calculatedexplicitly. The elementary calculations are in [16]. The three-line calculation when one specialisesto changing the k12 from 0 to 1 is given in the Appendix to this report.For networks with jAj � 5, the Wheatstone network is the only two-terminal instance where Braess'sparadox can occur. For jAj > 5, we have the results in x4.1.1 and in x4.1.2.We remark that the only instance when we have simultaneously a separable and a separable �is when both are power-laws. The converses of Theorem 4.3 which we have to date are �rstly forseparable , and secondly for separable �.An odd continuous function f : IR ! IR satis�es Assumption As if f is C1(0;1) and f 0 > 0 on(0;1). We remark that if � satis�es Assumption As, so does its inverse �.10

4.1.1 Separable THEOREM 4.4 Suppose that for all a 2 A, a is single-parameter separable and strict. Supposef�a j a 2 IN�g is a set of � functions, � � 6, all satisfying assumption As. Let N denote the setof all networks (N;A) with ow functions ka�'(a) on arcs a 2 A, where ' ranges over all the one-to-one mappings from A to IN� . Suppose that for any network (N;A) in N , and any b satisfying(3), the power-loss P is a nonincreasing function of each ka. Then there is an s > 0 such that forall a 2 IN� and t � 0, �a(t) = �a(1)ts.This is Theorem 2 of [5] and it follows immediately from the following result.THEOREM 4.5 Suppose that for all a 2 A, a is single-parameter separable and strict. Let N2denote the set of all two-terminal networks (N;A), with the terminal nodes being 0 and nt � jN j�1.Let � � 6 and let assumption As be satis�ed as in the preceding theorem. With, in the precedingtheorem, N replaced by N2 and b satisfying b(i) = 0 for i 6= 0; nt and b(nt) = �b(0), the conclusionof the preceding theorem holds. That is, if for any network (N;A) in N2 and any b satisfying thepreceding restriction the power-loss P is a nonincreasing function of each ka, then there is an s > 0such that for all a 2 IN� and t � 0, �a(t) = �a(1)ts.The preceding two theorems follow from their � = 6 versions. This is because when � > 6 we canconstruct networks which have e�ectively 6 arcs by (i) putting the surplus arcs in series and thenputting this series in parallel with one of the chosen 6, and (ii) setting the conductivity factors ofthe surplus arcs to zero.The proof of the � = 6 case of the last theorem is given in [5] depends in part on a detailed analysisof a general nonlinear Wheatstone bridge network.(If readers �nd the requirement � � 6 to be unpleasant, we remark that the theorems can bemodi�ed to allow its removal to � � 1. The modi�cation is to allow repetition of the conductivityfunctions in the elements in the test networks, i.e. to drop the requirement that ' be one-to-one.)4.1.2 Separable �Rather than, as in x4.1.1, �xing conductivity functions �a and varying conductivity factors ka, herewe �x resistance functions �a and vary resistivity factors ra.Here is a description of the results of subsection 4.1.1 in physical terms. It applies to Colebrook-White with La variations (but not, of course, to the case of Da variations which is the main concernin x4). Suppose that one is given � � 6 very large rolls of wire, with each roll being identi�able byits colour say. The rolls of wire are made of rather exotic electrical conductors. For each type of wirew, the form of the resistance function �w is constant along the wire in roll w. By cutting a lengthl of wire from roll w one obtains an element with resistance r(l)�w (with larger r(l) correspondingto longer l). Electrical networks can be built making their arcs from lengths lw of wire from roll wassembled in various kinds of networks. If, no matter which electrical network is built, the powerloss P increases when r is increased, the theorems of [6] state that the network would have to beof type (H).THEOREM 4.6 Suppose that for all a 2 A, �a is single-parameter separable and strict. Supposef�a j a 2 IN�g is a set of � functions, � � 6, all satisfying assumption As. Let N denote the setof all networks (N;A) with phia = ra�'(a) on arcs a 2 A, where ' ranges over all the one-to-onemappings from A to IN�. Suppose that for any network (N;A) in N , and any b satisfying (3),the power-loss P is a nondecreasing function of each ra. Then there is an s > 0 such that for alla 2 IN� and t � 0, �a(t) = �a(1)ts.This is Theorem 2 of [6] and it follows immediately from the following result.11

THEOREM 4.7 Suppose that for all a 2 A, �a is single-parameter separable and strict. Let N2denote the set of all two-terminal networks (N;A), with the terminal nodes being 0 and nt � jN j�1.Let � � 6 and let assumption As be satis�ed as in the preceding theorem. With, in the precedingtheorem, N replaced by N2 and b satisfying b(i) = 0 for i 6= 0; nt and b(nt) = �b(0), the conclusionof the preceding theorem holds. That is, if for any network (N;A) in N2 and any b satisfying thepreceding restriction the power-loss P is a nondecreasing function of each ra, then there is an s > 0such that for all a 2 IN� and t � 0, �a(t) = �a(1)ts.As in subsection 4.1.1, the preceding two theorems follow from their � = 6 versions. Again thisis because when � > 6 we can construct networks which have e�ectively 6 arcs by (i) putting thesurplus arcs in series and then putting this series in parallel with one of the chosen 6, and (ii)setting the conductivity factors of the surplus arcs to zero.The proof of the � = 6 case of the last theorem is given in [5] depends in part on a detailed analysisof a general nonlinear Wheatstone bridge network.4.2 Further results for networks with separable aThe following theorems (i) show, for the Wheatstone graph, that type (H) is not necessary to havepower-loss decreasing in ka, and (ii) show that � � 6 is necessary in Lemma 3.THEOREM 4.8 Consider the two-terminal Wheatstone bridge graph. Suppose that for all a 2A, a = ka� is single-parameter separable and strict, and � satis�es Assumption As. Here theconductivity functions on each arc are the same. If it is given that, at any k 2 IR5, k � 0 for whichthe network satis�es Assumption A+, the power-loss is a decreasing function of any ka then ln� isconcave on (0;1).THEOREM 4.9 Suppose that for all a 2 A, a is single-parameter separable and strict. Considerthe Wheatstone bridge graph, and two-terminal ows with b(1) = 0 = b(2). Suppose that the con-ductivity functions on each arc are the same, that is ka�a = ka� for some � satisfying AssumptionA(�). The power-loss doesn't increase if ka is increased(i) with no further restrictions on �, if a is any arc except (1; 2),(ii) if a is (1; 2) provided ln� is concave.4.3 A Colebrook-White exampleBraess's paradox cannot occur in the rough-pipe limit of (CW) when one neglects the c3 termcompared to the c2 term. (See [5, 6].) Thus, we think that good examples where the Braessphenomenon is noticeable will involve at least some pipes being smooth, K = 0, c2 = 0.Consider the Wheatstone bridge network, with ow in at node 0 and out at node 3, b(1) = 0 = b(2).Arcs (0; 1) and (2; 3) have parameters La and Da. Arcs (0; 2) and (1; 3) have parameters Lb andDb. In our numerical study (La; Da) corresponds to a long pipe with a large diameter, (Lb; Db)corresponds to a short pipe with a small diameter. We compare two particular extreme forms of theWheatstone network. In the �rst, the diameter of the pipe joining (1; 3) is zero (c1 = 0; t0 !1): inthe second, the diameter of the pipe joining (1; 3) is very, very large - e�ectively in�nite - (c1 !1,c3 = 0; t0 = 0).We illustrate the Braess phenomenon in the following form. We begin with sending a prescribed ow C0 through the �rst form of the network, so that, using the symmetry of the network, thepressure head at node 3 is (given p(0) = 0)p0(3) = �a(C02 ) + �b(C02 ):12

Here the subscripts a and b refer to the functions evaluated with (Da; La) and (Db; Lb) respectively.The subscript 0 reminds us that we are dealing with c1 = 0 on the central (1; 2) pipe. We nowapply the same pressure head across the second network, which has more installed capacity c1 !1.Then the ow becomes C1 = a(p0(3)2 ) + b(p0(3)2 ):We are interested in whether C1�C0 is positive or negative, and will see that both signs can occurwhen C0 > 0. It is easy to compute this as function of C0:C1 � C0 = a��a(C02 ) + �b(C02 )2 �+ b��a(C02 ) + �b(C02 )2 �� C0:We remind readers that, for smooth pipes, the functions �, are, for appropriate values of c1, c3determined from (D;L), (t) = � c1tpjtj ln� c3pjtj� for jtj > c23; �(q) = qjqj� jqjc1W ( jqjc1c3 )�2:We remark that the rhs of the �rst equation (de�ning ) can be written as a simple function oft3 = t=c23 : the rhs of the second equation (de�ning �) can be written as a simple function ofq1 = q=c1.See the appendix for numerical examples. These include the following.� In one set of computations we take K = 0 on all pipes and La = 100m, Da = 0:05m,Lb = 10m, Db = 0:1m.� In an other set of computations, arcs subscipted a are rough (c3 = 0) and arcs subscripted bare smooth (c2 = 0). This item is not yet completed.These �rst examples are not physically realistic as the ows are very small on some smooth arcs,and possibly t is just a little greater than c23. Work is in progress to attempt to �nd more realisticexamples: however, we do not expect the phenomenon to be more than a curiosity.5 IMPROVING THE COMPUTATIONS?The only network treated numerically in this report is the Wheatstone one. A special case usingthe formula from x1.1 was treated, in Maple, to give the results given in the preceding subsectionx4.3. The results for the optimization approach of x3 are described in an Appendix.The optimization formulation, problem (P), has advantages in not requiring a good starting guess.The author is defering any larger-scale computation until he has an established convex-optimizationsoftware exploiting the fact that the hessian is banded.For large-scale networks with hundreds, or even thousands, of nodes e�ciency may become an issue.This is particularly likely to be the case when repeated solution with the same network topologybut with di�erent diameters, say, is required. See x6. There are several ways that large-scalecomputation might be improved: our speculations are listed below.1. The sparsity might be exploited in ways which take into account the network origins of theproblem, e.g. the ASG algorithm. See [23].2. There may be some scope for Fortran/C code generation, e.g. from automatic di�erentiationor from Computer Algebra packages, in this general area. The production of Fortran/C codefor functions and derivatives could be made easier by code generators, an item which may bemore worthwhile as the number and variety of elements modelled - valves, pumps, etc. as wellas pipes - increases. However, of more use for improving e�ciency is that code generators cangenerate network-speci�c Fortran/C codes, so that, as one proceeds with Newton iteration no13

code testing which nodes are connected with which other nodes is necessary. This is discussedin [16, Part II] and in [17, 21, 22]. It is, however, conceded that the pipe network poroblemmay be too simple for the payo� to be important: signi�cant payo�s are more likely in thecontext of multicommodity ows.6 SENSITIVITY AND DIAMETER SELECTIONPart of the reason for recording the results concerning the explicit formulae related to the Colebrook-White ow law is that there are many calculations done for pipe networks. The equilibrium owcalculations are just an example. Investigations of sensitivity, and diameter selection are otherexamples, which, however, this author has only cursorily investigated.See [7], [23].

14

Appendix. NUMERICAL WORKAfter some preliminaries in this section, concerning a single pipe, the remainder of the numericalwork focuses on the Wheatstone bridge network.A single pipeThe correct place to begin the numerics is with physical quantities, and in this we use [4] as ourstandard reference. Viscosities of di�erent uids are given on p5 of [4]. The dynamic viscosity ofwater at 20 degrees C is nearly 1 centipoise, or 10�3 kg=(ms). The density of water is 103 kg=(m3).In MKS units the kinematic viscosity, at 20 degrees C is 10�6 (and at 10 degrees C it is 1:31�10�6).Test for function CWAn easy test problem for our Colebrook-White code is the following, based on problem 3.6 of p57of [4]. Suppose we have 103m of 0:5m concrete pipe, for which K=D = 0:006.(a) Find the ows, and Reynolds numbers, if the pressure head di�erence between the ends is given:tabulate this for pressure head di�erences 0.75m to 1.0m in steps of 0.05m. Do this (i) at 10 degreesC, (ii) at 20 degrees C.(b) Determine the head loss to pump water at 10 degrees C at a rate of 0.1 m3s�1 through thepipe.Here is the Matlab to solve part (a).%%g = 9.81%c =100% Use data similar to Brebbia p57 Problem 6L = 1000;d = 0.5;K = 0.006*d; % cement pipeglobal nunu = 0.00000131 % 10 degrees Cp = parameter1(L,d,K);global c1 c2 c3 t0c1 = p(1); c2 = p(2); c3 = p(3);t = 0.75:0.05:1.0q1 = CW1(t)r1 = reynol1(q1,d)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%nu = 0.000001 % 20 degrees Cp = parameter1(L,d,K);c1 = p(1); c2 = p(2); c3 = p(3);t = 0.75:0.05:1.0;q2 = CW1(t)r2 = reynol1(q2,d)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function p=parameter1(L,d,K);% calculate from geometric and physical quantities%constants 15

%nu=1.0e-6; % 20 degrees Cglobal nu;g=9.81;clg=pi*sqrt(g/2)/log(10);p(1)=clg*sqrt(d^5/L);p(2)=10*K/(37*d);p(3)=nu*sqrt(L/(2*g*d^3));%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function q = CW1(t)global c1 c2 c3 t0;st = sqrt(t);rt = 1./st;% 0.5.*log(10) = 1.1513if (abs(t)>t0)q = -c1.*st.*log(c2.*ones(1,length(t))+c3.*rt);elseq = 0;end;%% OUCH - Engineers used log_10%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function reyn = reynol1(q,d)global nu;cnst = 4.0/(nu*pi*d);reyn = cnst.*q;Here are the numerical results for part (a).t = 0.7500 0.8000 0.8500 0.9000 0.9500 1.0000q1 =0.0937 0.0968 0.0998 0.1027 0.1055 0.1082r1 =1.0e+05 *1.8216 1.8815 1.9395 1.9959 2.0508 2.1042nu = 16

1.0000e-06q2 =0.0938 0.0969 0.0998 0.1027 0.1056 0.1083r2 =1.0e+05 *2.3879 2.4664 2.5424 2.6163 2.6881 2.7581It can be observed that even though there is a 30% change in the viscosity (and hence also Re),there is only a tiny 0:1% change in the ow. The �gure, Moody diagram, indicating this behaviourin the fully rough zone is given in [4] p53, and later in this single-pipe part of the appendix.For part (b), which is the question asked in Brebbia, we �nd a slightly di�erent answer than Brebbia,namely that the pressure head di�erence is 0.853 compared to Brebbia's 0.845. However Brebbia'smethod is graphical and a 1% accurate read-o� from graphs is probably acceptable. Here is theMatlab.>> nu = 0.00000131 % 10 degrees Cnu =1.3100e-06>> p = parameter1(L,d,K);>> c1 = p(1); c2 = p(2); c3 = p(3); CW1(0.853)ans =0.1000The Re does agree with that given in Brebbia, 0:19� 106.CW and HWWe now check out how well CW and HW laws agree. Here we work with 100m of smooth (e.g.poly) pipe of 50mm diameter. The numbers are in approximate agreement (say within 30%) withSwan's table in `The Countryman Farm Handbook' (West Australian Newspapers).%%g = 9.81%c =100% Use data as in SwanL = 100;d = 0.05; % 50mm pipeK = 0.0; % Smooth pipes, e.g. polyglobal nunu = 0.000001 % 20 degrees Cp = parameter1(L,d,K); 17

global c1 c2 c3 t0c1 = p(1);c2 = p(2);c3 = p(3);t0 = (c3/(1-c2))^2;%t = 1:1:10;%tau = g.*t;%q2 = CW1(t)hold onfplot('CW1',[0.02 4])global chwchw = paramHW2(L,d,150);fplot('HW1',[0.02 4],'-+')xlabel(' Pressure-head difference')ylabel(' Flow')title(' Colebrook-White law, and Hazen Williams +')print CWHW1p -dpshold off%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%figure % Very large pressureshold onfplot('CW1',[1 100])fplot('HW1',[1 100],'-+')hold off%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Totally ridiculously large pressuresCW1(10^10)-HW1(10^10) % positive 7.7CW1(10^11)-HW1(10^11) % -14.8%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%figure % Very small pressureshold onfplot('CW1',[10^(-5) 0.04])fplot('HW1',[10^(-5) 0.04])xlabel(' Pressure-head difference')ylabel(' Flow')title(' Colebrook-White law, and Hazen Williams. Tiny flows')print CWHW2p -dpshold off% t0 = c3^2;figurehold onfplot('CW1',[t0/2 2*t0])fplot('HW1',[t0/2 2*t0],'-+')xlabel(' Pressure-head difference')ylabel(' Flow')title(' CW and HW +. Extremely tiny flows')print CWHW3p -dpshold off 18

function pHW = paramHW2(L,d,c)% Brebbia p54pHW = pi*0.3545*c*(d^2.63)/(4*(L^0.54));function qhw = HW1(t)global chw;qhw = chw.*t.^0.54;The numbers, in Figure 2 are in approximate agreement (say within 30%) with Swan's table in`The Countryman Farm Handbook' (West Australian Newspapers).0.5 1 1.5 2 2.5 3 3.5 4

0

0.5

1

1.5

2

2.5

3

3.5x 10

−3

Pressure−head difference

Flo

w

Colebrook−White law, and Hazen Williams +

Figure 2: Matlab plot, 100m of 50mm poly pipeThe plots in Figure 3 are for small ows and head-di�erences, one so small as to be unphysical,probably even being a low Re ow.19

0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.040

0.5

1

1.5

2

2.5x 10

−4

Pressure−head difference

Flo

w

Colebrook−White law, and Hazen Williams. Tiny flows

3 4 5 6 7 8

x 10−8

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2x 10

−7

Pressure−head difference

Flo

w

CW and HW +. Extremely tiny flows

Figure 3: Matlab plots, 100m or 50mm poly pipe20

The Moody diagramHere is Matlab giving the Moody diagram as in Figure 3.8 of [4].function [reynff] = reynff(ff)global epD;% CW for Re as function of friction factor, ff, and roughness ratio epD% epD = roughness/D; D is diameterrf = 1./sqrt(ff);% 0.5.*log(10) = 1.1513denom = 3.7.*exp(-1.1513.*rf)-epD.*ones(1,length(ff));reynff = 9.287.*rf./denom;%% OUCH - Engineers used log_10%lprint(solve(rf=-2*log((10*ep/(37*D))+251*rf/(100*R)),R));%9287*rf*D/(3700*exp(-1/2*rf)*D-1000*ep)%% FILE moodys.m%% Algorithm looks at equal spacing in friction factor f.% Consequently it doesn't cope well when f nearly constant in Re.%global epD;clfepD = 0.1;for i = 1:6%[X Y] = fplot('reynf01',[0.01 0.10 10^3 10^7], 'TOL=1.0e-4');epD = epD/10X = (0.01:0.0001:0.10);Y = reynff(X);loglog(Y,X)axis([10^3 10^7 0.01 0.10 ])hold onend;gridxlabel(' Reynolds Number Re')ylabel(' Friction factor f')title(' Moody diagram: Colebrook-White law')print moodyp -dpshold offSee Figure 4 for the result from Matlab's Postscript.21

103

104

105

106

107

10−2

10−1

Reynolds Number Re

Fric

tion

fact

or f

Moody diagram: Colebrook−White law

Figure 4: Matlab plot for KD = 10�1 down to 10�6, in geometric progression, ratio 10�1.22

NUMERICAL WORK FOR WHEATSTONE BRIDGEAll the numerics below need to be re-worked. The value of � seems to be wrong.Open questionsFor Hazen-Williams the power-consumed is a monotonic function of La. Is this still true forColebrook-White (not in general, but for Wheatstone network)? Note that �a is separable forvariations in La.The primal optimization problemHere is an example with a Wheatstone bridge solved by using Matlab's Optimization Toolbox,speci�cally the function fminu. See [30, 20].function c=CW0(x,C);% integrated Colebrook-White, Psi in Keady 95 paperxs=sqrt(abs(x));M=log(C(2)+C(3)/xs);R=C(3)/C(2);c=C(1)*( (-2/3)*M*(xs^3) - ...(2/3)*(R^3)*M - ...(1/3)*(R^3)*log(abs(x)) + ...(2/3)*(R^2)*xs - ...(1/3)*R*abs(x));function p=parameter(L,d);% calculate from geometric and physical quantities%constantsnu=1.5e-5; %At 20 degrees Cg=9.8;K=1.0e-5; %Roughness of PVC pipe% c0=0.1; %Total flow from earlier diagnosticsclg=pi*sqrt(g/2);p(1)=clg*sqrt(d^5/L);p(2)=10*K/(37*d);p(3)=nu*sqrt(L/(2*g*d^3));function f=fun(x);% fun.m - Objective function to be minimized, CW0 for Wheatstone bridgeglobal C01 C02 C12 C13 C23w=-1/10; % Consumption at node 3f=CW0(x(1),C01)+CW0(x(2),C02)+CW0((x(1)-x(2)),C12)+ ...CW0((x(1)-x(3)),C13)+CW0((x(2)-x(3)),C23)-w*x(3);23

clc;% script file for Wheatstone bridge with Colebrok-White flow law% Solved with fminu: primal optimization problem.global C01 C02 C12 C13 C23L01=100;L02=100;L12=141;L13=100;L23=100;d01=.05; % make same as d23 used to be 0.05d02=.1;d12=.1; %d13=d02; % with usual arrangementd23=d01; %C01=parameter(L01,d01);C02=parameter(L02,d02);C12=parameter(L12,d12);% C12=[0.0,2.0,2.0]; % if d12=0C13=parameter(L13,d13);C23=parameter(L23,d23);t1=clock;% Initial guess% x=[-38.6,-20.3,-59.0]; % failed, too close to soln?x=[-38.0,-20.0,-59.0];options(1)=1;type funx=fminu('fun',x,options)% Elapsed timet=etime(clock,t1)Results of doit.Optimization Terminated SuccessfullyGradient less than options(2)NO OF ITERATIONS=45x =-38.6833 -20.3922 -59.0755t = 8.4674Check using Maple fsolveHere we formulate the problem di�erently. We solve in exactly the same network and same param-eter values as with the Matlab which immediately precedes this.24

# Pipe network problem: Wheatstone bridge network# Colebrook-White## Basically, just an fsolve# Notation: 4 nodes, initial node 0, final node 3## This version with symmetry.## Global constantsnu:= 15/(10^6); # m^2/sec at 20degreesC 0.000015g:= 98/10; # m/sec^2K:= 1/(10^5): # m roughness PVC pipe 0.0001C0:= 1/10; # m^3/sec total flowc1g:= Pi*sqrt(g/2); # Check log_10 factor# Global functionsc1 := (L,D) -> c1g*sqrt(D^5/L);c2 := D -> 10*K/(37*D);c3 := (L,D) -> nu*sqrt(L/(2*g*D^3));psi:= proc(t,L,D)local at1,st1,tmpl,ans,c1l,c2l,c3l:at1:=abs(t): st1:=t/at1:c1l:= c1(L,D): c2l:= c2(D): c3l:=c3(L,D):tmpl:= min(1,c2l+ (c3l/sqrt(at1))):-st1*c1l*sqrt(at1)*log(tmpl);end:phi:= proc(u,L,D)local u1,au1,su1,c2d3:c1l:= c1(L,D): c2l:= c2(D): c3l:=c3(L,D):u1:=u/c1l: au1:=abs(u1): su1:=u1/au1:c2d3:=c2l/c3l:su1/(c2d3-W(au1*(exp(c2d3*au1))/c3l)/au1)^2;end:# Global expressionsh1:=-x:xmax:=phi(C0,L0,D01):q01:=psi(x,L0,D01):q02:=C0-q01:h2:= -phi(q02,L0,D02):q12:=psi(h1-h2,L12,D12):eqn:=q02+q12-q01:h3:=h1+h2:hfinal:= proc(Dv)local eq1,xmax1,xs: 25

eq1:=subs(L0=100,L12=141,D01=5/100,D12=Dv,D02=1/10,eqn):xmax1:=evalf(subs(L0=100,D01=5/100,xmax)):xs:=fsolve(eq1=0,x,1/100..xmax1):[xs,evalf(subs(L0=100,L1=141,D01=5/100,D12=Dv,D02=1/10,x=xs,h3))]:end:printlevel:=0:h1list:=[]: h3list:=[]: nvals:=2;for j from 1 to nvals doprintlevel:=0:hns:=hfinal(j/20):h1list:=[op(h1list),-hns[1]];printlevel:=1;h3list:=[op(h3list),hns[2]];od;h2list:=[seq(h3list[j]-h1list[j],j=1..nvals)];q01list:=[seq(evalf(psi(h1list[j],100,5/100)),j=1..nvals)];q02list:=[seq(evalf(psi(h2list[j],100,1/10)),j=1..nvals)];q12list:=[seq(evalf(psi(h1list[j]-h2list[j],141,j/20)),j=1..nvals)];chklist:=[seq(q01list[j]+q12list[j]-q02list[j],j=1..nvals)];## h1list; [-131.9568615, -38.68324151]# h3list := [-145.0782216, -59.07550128]# h2list := [-13.1213601, -20.39225977]# q01list := [-.03596618153, -.01818616458]# q02list := [-.06403381870, -.08181383546]# q12list := [-.02806763701, -.06362767084]#hlist:=[seq([j,h3list[j]],j=1 .. nvals)];interface(plotdevice=x11):# above for output to screen# below for output to file# interface(plotdevice=postscript,plotoutput=plotfile):## plot(hlist,x=1..nvals,style=POINT);# ghostview to look at postscript file#26

Numerical examples of Braess behaviourNeeds to be recomputed. However, perhaps it is impossibly hard to �nd a realisticexample with this network. Flows here are un-physically small.Power-law nonlinearitiesPower law nonlinearitiesThe method expounded for Braess behaviour for Wheatstone bridge arcs was applied earlier, [5],for power laws with di�erent powers of di�erent arcs. In [5] we considered the example a(t) = �(t; 1=2); �a(q) = �(q; 2); b(t) = t; �b(q) = q:This would be obtained by replacing s:=27/50 by s:=1 in the following.# C_infty - C_0 denoted CdCh:= C0/2:# Example in CK was s:=1s:= 27/50:Cd:= sqrt((Ch^2 + Ch^(1/s))/2) + ((Ch^2 + Ch^(1/s))/2)^s -C0:interface(plotdevice=postscript,plotoutput=plotfile):plot(Cd/C0, C0 = 1/2 .. 4);In Figure 5 the region below the axis, i.e. with C0 a little below 2 is the region where Braessparadox is exhibited.-0.02

0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

0.18

0.5 1 1.5 2 2.5 3 3.5 4C0Figure 5: Plot for s = 1=2 and s = 1The rough-pipe limit of CW is s = 1=2. The smooth-pipe regime has a di�erent behaviour, and,here we pretend that it might be like s = 27=50. In Figure 6 we give the result of running thecode we listed at the beginning of this subsection. The Braess phenomenon is still there, just notas large as before. 27

0

0.0005

0.001

0.0015

0.002

0.5 1 1.5 2 2.5 3 3.5 4C0Figure 6: Plot for s = 1=2 and s = 27=5028

Colebrook-White BraessThis is introductory to the actual numerical work in the following two subsections. Here is thegeneral approach.1. Fix Da, Db, La.2. Fix C�, say C� = 0:02.3. Find ta so that a(ta) = C�=2.4. Find Lb so that b(ta) = C�=2.5. With these values of Da, Db, La, Lb, plot C1 � C0 for C0 between C�=2 and 2C�.At present the numerical approach below does not implement this: it is just ad hoc. It will berevised to be more systematic, as this may be necessary to increase the chance of �nding a morerealistic example.In the two subsections immediately below we treat1. all pipes smooth;2. pipes a smooth and pipes b rough.Colebrook-White Braess, all pipes smoothThis is not a good example. Flows are ridiculously small.# Colebrook-White for smooth pipes, K=0# Illustration of Braess-style paradox:# "It is possible to provide extra capacity on a link -# a wider diameter pipe - yet get less flow through at# same pressure differences."# Remark: This does not happen with Hazen-Williams.## Here we consider Wheatstone-bridge graphs.# C0 flows in at node 0 and out at node 3.# Pipes (0,1) and (2,3) are the same length and diameter:# long and with a large diameter.# Pipes (0,2) and (1,3) are the same length and diameter.# short and with a small diameter.## We compare two particular extreme forms of Wheatstone-bridge graphs.# In the first the diameter of the pipe joining (1,3) is 0.# In the second the diameter of the pipe joining (1,3) is infinite.#Digits:= 20; # to preserve against possible subtractive cancellationg:=98/10:c1:= (L,D) -> Pi*sqrt(g*(D^5)/(2*L));nu:= 15/10^6:c3:= (L,D) -> nu*sqrt(L/(2*g*D^3));psi:= proc(t,c1LD,c3LD)local abt,st: 29

abt:=abs(t): st:=t/abt: s2t:=sqrt(abt):-st*s2t*c1LD*log(c3LD/s2t);end:phi:=proc(q,c1LD,c3LD)local cq,aq,sq:aq:=abs(q): sq:=q/aq:cq:=aq/c1LD:sq*(cq/W(cq/c3LD))^2;end:c11:=evalf(c1(100,1/10)); c31:=evalf(c3(100,1/10));c12:=evalf(c1(5,1/20)); c32:=evalf(c3(5,1/20));Ci0:= proc(C0,c11g,c31g,c12g,c32g)local h1,h2,h3,q1,q2:h1:= evalf(phi(C0/2,c11g,c31g)):h2:= evalf(phi(C0/2,c12g,c32g)):h3:= h1+h2:q1:= evalf(psi(h3/2,c11g,c31g)):q2:= evalf(psi(h3/2,c12g,c32g)):[h1,h2,h3,q1,q2];end:test:= C0 -> Ci0(C0,c11,c31,c12,c32):# When Digits was 20, results agree very well# Results ########################################################### test(4/10^5);# [.0003065948939, .0003020509501, .0006086458440, .00001988064212,# .00002011101244]# "[4] + "[5] - 4/10^5;# -8# -.834544*10###################### a5:=test(1/10^5);# a5 := [.00006673836989, .00005425473825, .0001209931081,# .4494882299*10^(-5) , .5512616171*10^(-5) ]# a5[4]+a5[5]-1/10^5;# -8# .749847*10####################### a4:=test(1/10^4);#a4 := [.001024876856, .001101602335, .002126479191, .00005134383823,# .00004880396209]# a4[4]+a4[5]-1/10^4;# -6# .1478003*10# 30

####################for j from 1 to 9 doa.j:=test(j/10^5);a.j[4]+a.j[5] -j/10^5;od;# a2, a3, and a4 give negative values, others positive.# This result continues when Digits is increased to 20# (done as a check).###################################################### The flow as a function of head difference has a slightly# peculiar shape near t = c3^2.# Below c3^2 the flow is ZERO.# Just above c3^2 the flow increases, and is a concave function# of t.# I think these numbers are down in this rather peculiar region.# Next exercise: seek more realistic example of the effect.# Caution. Might not be easy. Effect is likely to be small.# Caution is because of Hazen-Williams being used in real# pipe networks.Colebrook-White Braess, half the pipes smoothNot done yetPipes b are rough, q = b(t) = (c1 ln(c2))t=sqrtt:Lb small and Db large, say Db = 1. � = 10�6 helps make c3 << c2. Rough concrete hasK=D = 0:01.31

NUMERICAL WORK FOR 6-NODE LADDERNot done yet. Aim is to �nd a more realistic example of a Braess paradox. The ideais that series combinations of di�erent diameter CW pipes might have q(t) behavioursu�ciently di�erent from CW that it will make �nding Braess paradox in it easier.References[1] D.P. Bertsekas, P.A. Hosein and P. Tseng, \Relaxation methods for network ow problemswith convex arc costs", SIAM J. Control and Opt. 25 (1987), 1219-1243.[2] N. Bolland, C.J. Goh and A.I. Mees, \An algorithm for nonlinear network programming:implementation, results and comparisons", J. Operational Res. 43 (1992), 979-992.[3] D. Braess, \Uber ein paradoxen der verkehrsplannung", Unternehmenforschung 12 (1968),256-268.[4] C.A. Brebbia and A.J. Ferrante, Computational hydraulics (Butterworth, London: 1983).[5] B. Calvert and G. Keady, \Braess's paradox and power-law nonlinearities in networks", Jnlof the Australian Math. Soc., 35B, pp1-22.[6] B. Calvert and G. Keady, \Braess's paradox and power-law nonlinearities in networks, PartII", To appear in Proceedings of the First World Congress of Nonlinear Analysts, WC528,Aug. 92, editor W. de Gruyter, editorial secretary D. Harnish.[7] S. Carr and G.J. Savage, \Sensitivity analysis of nonlinear physical systems using Maple".Pp 118-127 in Mathematical computation with Maple V: Ideas and applications, 1993, editorT. Lee. (Birkhauser.)[8] H.-J. Chen, \An exact solution to the Colebrook equation", Chem. Eng. 94 (1987), 196-198.[9] R.M. Corless, G.H. Gonnet, D.E.G. Hare and D.J. Je�rey, \Lambert's W-function in Maple",Maple Technical Newsletter 9 (1993), 12-22.[10] J.E. Cohen and P. Horowitz, \Paradoxical behaviour of mechanical and electrical networks",Nature 352 (1991), 699-701.[11] M. Collins, R. Cooper, J. Helgason, J. Kennington and L. LeBlanc, \Solving the pipe networkanalysis problem using optimization techniques", Mgmt. Sci. 24 (1978), pp747-760.[12] R.S. Dembo, J.M. Mulvey and S.A. Zenios, \Large scale nonlinear network models and theirapplication", Operations Research 37 (1989), 353-372.[13] F.N. Fritsch, R.E. Schafer, and W.P. Crowley, \Algorithm 443: solution of the transcendentalequation wew = x", Comm. ACM 16 (1973), pp123-124.[14] M.A. Hall, \Hydraulic network analysis using generalised geometric programming",Networks,6 (1976), 105-131.[15] I.E. Idelchik, Handbook of hydraulic resistance. (Hemisphere publishing, 2nd ed.: 1986).[16] G. Keady, \SENAC and other symbolic front ends for subsequent numeric computation",Parts I to IV, University of Waikato, Mathematics Research Reports, (1991).32

[17] G. Keady, \FORTRAN subroutines produced from Computer Algebra systems: using GEN-TRANs from REDUCE and from Macsyma". Pp 265-272 in Proceedings CTAC-91. Computa-tional Techniques and Applications Conference, Adelaide, July 1991, ed. J. Noye, B. Benjaminand L. Colgan.[18] G. Keady, \The Colebrook-White formula for pipe networks". Submitted to Civil Eng. Trans-actions, Inst. Eng. Australia.[19] G. Keady and M. Laine, \Using genmex { examples of calls to NAG". Document availableby anonymous ftp, directory /pub/keady.[20] G. Keady and 3rd year students, \Matlab, genmex, etc. - Applications to equilibrium owsin pipe networks". Document available by anonymous ftp, directory /pub/keady.[21] G. Keady and G. Nolan, \Production of Argument SubPrograms in the Axiom-NAG link:examples involving nonlinear systems". Pp 13-32 in Proceedings of the Workshop on Symbolicand Numeric Computation, Helsinki, May 1993. (Ed. H. Apiola, M. Laine and E. Valkeila).Rolf Nevanlinna Research Report series, B10.[22] G. Keady and M.G. Richardson, \An application of IRENA to systems of nonlinear equationsarising in equilibrium ows in networks". Pp 311-320 in Proceedings of the Symposium onSymbolic and Algebraic Computing, ISSAC '93. (Ed. M. Bronstein). ACM, 1993.[23] P. Neame, \Honours Student project: Some design problems in water distribution networks".(Supervised by C.J. Goh and G. Keady. Oct. 1994.)[24] A.L. Peressini, F.E. Sullivan and J.J. Uhl, The mathematics of nonlinear programming.(Springer-Verlag, U.S.A.: 1988).[25] A.J. Reynolds, Turbulent ows in engineering. (Wiley: 1974)[26] R.T. Rockafellar, Conjugate duality and optimization (SIAM, Philadelphia: 1974).[27] R.T. Rockafellar, Network ows and monotropic optimization (Wiley, New York: 1984).[28] G. Strang, \A framework for equilibrium equations", SIAM Review 30 (1988), 283-297.[29] A.M. Watts, \Variational principles, duality, Legendre transformations and mine shaft venti-lation", J. Austral. Math. Soc. 34B (1993), pp274-281.[30] J. Yuryevich, \Student project: Pipe networks using Matlab's optimization toolbox". (Super-vised by G. Keady. May 1994.)[31] C.A. Zukowski and J.L. Wyatt, \Sensitivity of nonlinear one-port resistor networks", IEEETrans. Circuits and Systems 31 (1984), 1048-1051.33