THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING WATERHAMMER ANALYSIS OF DISTRIBUTION SYSTEMS Victor L. Streeter October, 1966 IP-748

ACKNOWLEDGEMENTS This study has been supported by NSF Grant GP-340 to the University of Michigan. iii

TABLE OF CONTENTS Page ACKNOWLEDGEMENTS.......................................... iii LIST OF FIGURES......... oo.........O O...O.................... vii INTRODUCTION................ oo o o oooo....................... 1 HARDY CROSS METHOD FOR BALANCING HEADS........oo............... 3 ALGEBRAIC WATERHAMMER EQUATIONS.oo..o..oo......oo.oo........... 6 APPLICATION OF THE ALGEBRAIC METHOD TO A TYPICAL JUNCTION...... 8 DESCRIPTION OF COMPUTER PROGRAM AND EXAMPLE PROBLEMo........... 11 SUMARY AND CONCLUSIONS....... o.... o... o.................. 00 17 REFERENCES............ 00... o o O............................. 17 iv

LIST OF FIGURES Figure Page 1 Simple circuit to which head balancing is applied..... 18 2 Pipe of length L and wave speed a, with positive flow direction shown by arrow........................ 18 5 Junction 91 having four pipes and a valve outlet...... 18 4 Example system comprised of 45 pipes, 26 junctions, 20 simple circuits and 19 valved outlets. Junction numbers are encircled. Assumed positive flow direction shown by arrows.................................. 19 5 Computer program for balancing heads by Hardy Cross Method and for calculating transients by algebraic waterhammer equations................................. 20 6 Input data for first read statement.................. 21 7 Input data for seconds read statement................. 21 8 Transient solution for problem with data given in Figures 6 and 7................................e 22 9 Transients caused by rupture at junction 21 when reservoir is at 200 feet.......................... 23 vii

INTRODUCTION The great speed of a modern large digital computer permits very complex piping systems to be analyzed for transient flow provided a method for describing the network can be adapted to computer usage. For most systems the friction may be adequately treated either as proportional to the square of the velocity, or to some exponential value of the velocity. In the algebraic method, the end boundary condition equations for the C and C- characteristicsl developed from the method of characteristics equations are applied over the whole pipe length between two junctions. In effect this "lumps" the friction for the pipe as that given by the velocity at the opposite end of the pipe one wave-travel time previous to calculation of events at a junction. In case certain pipes are very long compared with others, any interior sections of those pipes may be considered as junctions. The combination of algebraic waterhammer equations plus indexing of the pipes leaving every junction permits a simple program to be prepared for complex systems. Each pipe is divided into a number Pi of equal reaches such that the wave travel time Dt for each reach is the same for every reach in the system. Hence Pi Dt = Li/ai in which Li is the length of pipe i and ai is the pulse wave speed in that pipe. This is a limitation on the method, but since Dt may be selected so that a pipe has up to 10 reaches, suitable values of Pi can usually be found for each pipe by minor adjustments of pipe lengths and of wave speeds. In general, wave speeds are not known with precision and small adjustments are permissible. Steady-state initial conditions in the piping system must first be established before solving the transient problem. For very large systems -1

-2the steady-state heads and flows should be found separately to conserve storage for the transient solution. For the example given both the steadystate determination and the transient solution are combined. The program will handle head-balancing and transient solution for up to about 300 pipes. Many works have been published on steady-state solution of networks; however, for the purposes of illustrating the transient methods, the Hardy Cross method- of balancing heads is used in this treatment. The method of balancing heads is first discussed, then the algebraic equations are developed. They are next applied to a typical junction of a network. The computer program is described and then an example problem is given consisting of 20 simple circuits, 26 junctions and 45 pipes.

HARDY CROSS METHOD FOR BALANCING HEADS The method for balancing heads, as developed by Hardy Cross, is discussed briefly with regard to its application to use by the computer. in this method, the elevation of hydraulic gradeline is known at one point in the system, and the inflows and outflows from the system are given. In addition the resistance properties of each pipe must be known. The head loss H. in pipe i is given by H = Ri Qi I n~l (1) in which R. is the head drop over the length of pipe for unit discharge; ni is the exponent in the head-loss equation, Qi is the discharge. Each pipe in the network is numbered and a positive flow direction indicated in an arbitrary manrner. If the flow changes direction, Equation (1) changes;:;i!n. To start the balancing procedure which yields the correct flow in each pipe, flow is assumed through each pipe such that continuity is satisfied at every junction. This can be accomplished rapidly by starting at a junction with known inflow. Many pipes may be assumed to have zero flow, since it is not important that values close to the final correct values be assumed. Continuity must be exactly satisfied at every junction, however, because it is not upset or improved by subsequent steps in the balancing procedure. The complex network is treated as a collection of simple circuits. In balancing a simple circuit, Equation (1) is applied to each pipe and the net head loss around the circuit is determined. In Figure 1 pipes are shown with their assumed flow directions; as Q's have been assumed for each pipe, the head drop around the system is given by 4 4 n.-l =1i = i-l Ri Qilil (2) -3

-4The clockwise direction may be assumed as positive; if ZHi = O, then the heads are balanced and no correction to the flows are needed. When ZHi does not equal zero, a correction to the discharge AQ is calculated and applied to every pipe in the simple circuit in either the clockwise or the counterclockwise direction depending upon the sign of Equation (2). To obtain the correction AQ, Equation (1) is differentiated AH. _~1, n.-l - = ni Ri IQ 1 i ( AQ. 1 (3) and AH./AQ is summed for every pipe in the simple circuit with all terms 1 i positive. Then the corrective flow AQ is given by Hi ZR i QiIQil ni-l AQ= ---- (4) Z(AHi/AQ) Z i KQi 1 If the ZHi in the clockwise direction is positive, then the flow in each pipe should be decreased by AQ in the clockwise direction, e.g. in Figure 1 the flow in pipe 1 would be decreased by AQ and the flow in pipes 2,3, and 4 would be increased by the same AQ. Each simple circuit is corrected in this manner, in any convenient sequence; then the procedure is repeated until the sum of the absolute values of AQ for all simple circuits becomes less than some predetermined tolerance. After the heads are balanced in all simple circuits, the elevation of hydraulic gradeline is calculated for each junction, starting with the junction with given head. For determination of steady-state flows for different boundary conditions, other balancing tedhniques may be needed. The transient equations, developed in the next section, may also be used to find steady-state initial conditions.' Whether this is an economical way depends upon the

-5relative cost of engineering time to program the steady-state solution compared with cost of computer time. To do this, for example, all flows could be taken as zero, all heads equal to reservoir head, then the assumed outflows started instantaneously at time zero. Steady-state conditions will be approached after several round-trip wave travel times across the network.

ALGEBRAIC WATERHAMMER EQUATIONS Solution of the partial differential equations of waterhammer by the method of characteristics2 leads to two finite difference equations, one valid along the Ca characteristic and the other valid along the Ccharacteristic. When expressed in terms of discharge Q and elevation of hydraulic gradeline H, applied to the pipe of Figure 2, they take the form (5) C+ H (t) = H (t L) _ a Q (t) - Q (t L) I - Q1 L ) aH Q - BA gA LB A a 2D gA2 Aa Aa a (6) nC-: H (t) = H (t-L) + a FQ(t) - Q(taL) + a Q a L)Q(tL) L A B a gA A B aj 2DgA a A a a a is the wave speed, A the pipe cross sectional area, f the friction factor, L the pipe length, t the time, D the pipe diameter, and g the acceleration due to gravity. These equations are valid for any time t, in steady or unsteady flow, as long as the pipe remains filled with liquid. In general, each equation has two unknowns, QB and HB at time t for Equation (5), and QA and H at time t for Equation (6). A For algebraic waterhammer these equations are simplified by collecting the constants. B = gAi Ri fiL (7) 1- 2 1i 2gDiAi The term Ri is the same head loss coefficient used in balancing head in the Hardy Cross method. Each pipe is divided into P. equal reaches such that Pi Dt = i (8) ai -6

-7with Dt the time increment for the wave to travel through any reach of the system. The discharge Q for a pipe is triple subscripted Q(J1,J2,J) in which Jl is 2 if the downstream end of the pipe is designated and is 1 if the upstream end of the pipe is designated, J2 is an integer representing the pipe number, and J is an integer representing the time T = (J + JC)Dt (9) with JC a constant integer. The junction heads are double subscripted, H(E,J) with E the junction number and J the time counter. In this notation, Equations (5) and (6) become C+: Q(2,J2,J) = Q(1,J2,J-P(J2)) + B(J2) [H(EA,J-P(J2))-R(J2)Q(l,J2,J-P(J2)) IQ(, J2 J-P(J2))In-1] - B(J2) H(EB,J) (10) ca: Q(1,J2,J) = Q(2,J2,J-P(J2)) - B(J2) [H(EB,J-P(J2))+R(J2)Q(2,J2,J-P(J2)) |Q(2,J2,J-P(J2))|n-l] + B(J2) H(EA,J) (11) These two equations, which are available for each pipe at any time, are used together with continuity at a junction, and a valve equation (if an outlet is present) to calculate the elevation of hydraulic gradtine at the junction and the flow through each pipe into the junction.

APPLICATION OF THE ALGEBRAIC METHOD TO A T:PICAL JUNCTION In Figure 3 junction 91 has four pipes connecting it with a valve outlet; two of the pipes, 7 and 13 have the sign convention for flow into the junction while the other two, 17 and 23 have the sign convention indicating flow out of the junction. The head H(91,J) and the flows Q(2,7,J), Q(2,13,J), Q(1,17,J) and Q(1,23,J) are unknowns to be found at the current time. There is one equation for each pipe, Equation (10) for pipes 7 and 13 and Equation (11) for pipes 17 and 23 as well as the equation of continuity for the junction, so for no valve outlet five equations in five unknowns are at hand. With a valve, as shown, an additional equation for flow through the valve QV(91,J) for given opening is available for determination of the additional unknown QV(E,J). In Equation (10), applied to pipe 7 Q(2,7,J) = Cl(7) - B(7)H(91,J) (12) in which Cl(7) is comprised of the first two terms on the right-hand side of Equation (10), all being known from previous calculations. For pipe 17, Equation (11) applies and yields for flow into the junction -Q(l,17 J) Cl(17') - B(17) H(91,J) (13) in which C1(17) is the negative of the first two terms on the right-hand side of Equation (11), a known quantity. For no valve at the junction the sum of the discharges into the junction must be zero. After summing up all the flows into the junction ZQi = 0 = ECli - H(91,J)jBi (14) Hence H(9l,J) -- Cli (15) ZBi Then from Equations (12) and (13) the discharges are calculated. -8

-9If junction 91 of Figure 3 has an outlet to atmosphere controlled by a valve, then the net inflow to the junction through the four pipes must just equal the flow through the valve. It is convenient to express the valve opening by a dimensionless coefficient T. First by writing the orifice equation for the valve for steady-state QW(E) = (CDAV)O g HH(E) (16) in which, for convenience, the network is assumed to be in a horizontal plane and the elevation datum for hydraulic gradelines is the system elevation. In general for non-horizontal systems QW(E) = (CDAV)O 2g(HH(E) -EL(E)) (16a) in which EL(E) is the junction elevation based on the same datum as the hydraulic gradelines. For any flow through the system QV(E,J) = C 2g H(E,J) (17) or. for the more general case QV(E,J) = CDAV 2g (H(E,J) - EL(E)) (17a) After dividing Equation (17) by Equation (16) QV(E,J) = T VCN(E) %H(EJ) = CV IH(E,J) (18) or dividing Equation (17a) by Equation (16a) QV(E,J) = T VC(E) JH(E,J) - EL(E) (18a) The dimensionless valve coefficient T is given by CDA V 7T (19) (CD AV)D in either case and the other constants are collected into VCN(E) = Q v(E) (20) %IHH(E) or VC(E) = - v(EEL) (20a)

-10T is unity for steady-state flow and goes to zero as the valve is closed, CV in equation (18) is T VCN(E). By writing the continuity equation including the valve - QV(9,J) - 91) - CV (91J) 0 = zc (91J), - CV JH(9,J (21) whih yields a quadratic equation in H(91,J) only and may easily be solved. The transient solution starts with steady-state at time t - 0, the time is incremented to Dt, then each junction is solved for heads and flows as indicated, in any order. Special boundary conditions are then solved, such as the flows from junctions of fixed hydraulic gradelines. The results may then be printed out, the time incremented by Dt and the procedure repeated, taking the valve openings as function of t, i.e. T(EJ), Special boundary conditions such as one or more pumping stations may be included.

DESCRIPTION OF COMPUTER PROGRAM AND EXAMPLE PROBLEM A network, Figure 4, is used to illustrate the methods. It contains 20 simple circuits, 26 junctions, and 45 pipes. All inflow is at junction 26 where the elevation of hydraulic gradeline is held at 500 feet by the presence of a reservoir. Ten outlets, numbered 1 to 10 discharge from the system to atmosphere through valves. These valves may be stroked in any manner, the T~time data being supplied as input data to the program. Since Hardy Cross balancing may be accomplished using relative values of flow and of pipe resistances, the inflow is taken as 100 and the desired steady-state outflow from the ten outlets is given as input data that must add to exactly 100. Each pipe is numbered in any manner with the pipe number shown at its midpoint. An arrow is drawn along each pipe to establish a positive flow direction. The relative pipe resistances R(l).. o R(45) are given in the input data. They may be the actual resistances, from Equation (1), calculated from Darcy Weisbach friction factors or from the various exponential formulas. The exponent may vary for each pipe as well as the resistance coefficient. The steady-state flows through each pipe and the junction head determinations are first consideredo By starting at junction 26 flows are assumed through each pipe such that continuity is satisfied at each junction. These are shown in the input data as QQ(1) - 60o, 40o, 10,, etc. listing the flow through each pipe from 1 through 45. The apportioning of flows is rapidly accomplished and does not need to be close to the final answer for the method to converge. One may, for example, select flows such that 28 of the lines are initially assumed to have no flow, and still satisfy continuity at every junctbno -11

-12The computer program, written in MAD compiler language (Michigan Algorithm Decoder), Figure 5, contains two read statements. The first read statement, No. 001 inputs circuit information for the Hardy Cross steadystate balancing. The input data for this statement is shown in Figure 60 The first part of the program, statements 001 through 007 solves for the steady-state initial conditions. Statement 002 is used in locating integers in the series X(l) = 1,2,2.. o which describes the simple circuits, For example G(4) is the number of integers neededto describe the 53sided circuits, G(5) equals G(4) plus the number of integers needed to describe the 4-sided circuits, etc. Statement 003 sets all exponents equal to the given exponent NN in the head loss term. Each NN(I) may be placed in storage directly with input data if they varyo Statement 004 is a compound iteration statement that solves Equations (2), (3), and (4) then corrects the flows, checks to determine if the tolerance has been met, and continues to balance all the circuits until it has gone through them Z times or until the desired tolerance is reached. Statement 005 converts the relative R's by multiplying by SS, and converts the relative discharges to actual discharges QQ(I) by multiplying by SFo Statement 006 calculates the heads at each junction of the system. The indexing system for specifying the circuits starts with X(1) = 1,2,2,12, 11l 3 These first 6 integers specify the 3-sided circuit comprised of pipes 12, 11, and 35 The first number 1 indicates flow in pipe 12 is clockwise, the seconds number 2 indicates the flow in pipe 11 is counterclockwise and the third number 2 indicates the flow in pipe 3 is counterclockwise. In the second line of statement 004 starting J = 0,1, etc., W determines the clockwise or counterclockwise term, and Y the pipe

-15number. QP collects certain terms that are used in the next two terms, DH collects EHi and HD is Z dHi/dQi. Then DQ the correction is calculated, and the absolute DQ's summed for all circuits. Finally the Q's are corrected for each circuit. This is continued until the tolerance is satisfied or Z is reached. In accomplishing the balancing, statement 004 takes up all 3-sided circuits in the order given by the data in X(1) =..., then continues to the 4-sided circuits, etc. Statement 006 uses the indexing given by the data for XX(1) = 1,1, 12,5,1,14... H = 500 is given as elevation of hydraulic gradeline at junction 26. With reference to Figure 4, the first 3 numbers of XX(1) indicate in order, the pipe number, the flow direction (1 if toward the next junction, 2 if toward the previous junction) and the junction at which the elevation of hydraulic gradeline is being computed. Similarly 5,1,14 gives indexing information for calculation of head at junction 14, etc. By judicious selection of order of junctions, all elevations of hydraulic gradelines are calculated and the final one is junction 26 which checks the program. In statement oo6, L is the junction number, Y is the sign on the head loss term and Z is the pipe number. The print results statement, 007, prints out the final steadystate discharges through each pipe and the elevation of hydraulic gradeline at each junction. The second read statement, 008, provides the additional data needed for the transient calculation. These data are shown in Figure 7. Statement 009 sets initial high values for minimum HGL elevation and low values for maximum HGL elevations, which are subsequently adjusted to their proper values. Statement 010 is similar to statement 002, now acting for the 2nd set of X(I) data read in. It computes the number of integers in

-14X(I), i.e. G(2) = 0 is given in the data, G(3) is the number of junctions with 2 pipes times the number of integers needed to describe the junction. For a junction having two pipes 8 numbers are needed. X(l) = 10,1,7,14,1, 8,15,2 states, in order, junction 10, 1 means a valve is present (0 for no valve), pipe 7, connected to juncticn14 has flow (1) into junction 10; pipe 8, connected to junction 15 has flow (2) out of junction 10. The indexing data for all other junctions having two pipes follow in any order; then all the junctions having three pipes are listed etc. Hence, G(4) represents the number of all the integers listed in X(I) up to the beginning of the listing for 4-pipe junctions. Statement 011 gives constants needed to handle the valves. DTA(I) is the spacing, in time, between data points in the tau-time data TA(I,0) = 1.,... listed in Figure 7. VCN(I) calculates Equation (20). Since the algebraic waterhammer method needs values of H and Q, Li/Ai seconds earlier than the current calculation, it is necessary to place in core storage the steady-state heads and Q's for the previous JJ increments of DT, where JJ is the largest Pi for any pipe,plus one. Statement 012 accomplishes this objective for discharges, and statement 013 does the same for heads. Statements 014, 015, and 016 sets the headings for the transient output. Statement 017 is the heart of the transient program, as it, together with the internal functions Pt. and AT., calculates transient heads and discharges at each junction. Jt is the number of time steps DT that may be taken without exceeding storage requirements. Statements 018 through 022 move the calculations back in the storage space and causes 017 to be executed again and again until the maximum time TM of the transient is obtained.

-15In reading statement 017, each opening parentheses starts an iterative through statement; each closing parentheses to one of these statements is the end of the iteration. J represents the time steps, M is used to denote the number of pipes at junction under consideration and K locates data in X(I). A is the ECli in Equation (14) and C is EBi in the same equation. E denotes junctionnumber. L is the counter in proceeding through all the pipes at a junction, from 1 to M, used in determining A and C. Cl(L) is the first two terms on the right hand side of Equation (10) or the negative of the corresponding terms of Equation (11). The itera tion through L for a single junction closes just before AT. in line number 4. AT.(T) is the internal function statements 026 through 041 that performs a parabolic interpolation for T and then calculates the head, and the discharge through the valves at the junction. The next L iteration calculates the discharges through the pipes at the junction. The M, K, and L loops close here, then the PT. internal function solves for special boundary conditions and prints out results. For the example problem, arbitrary valve closures have been selected for all ten of the valves, rapid enough so that a rather violent pressure fluctuation occurs. The computer output is shown in Figure 8 - - - it has been selected in an arbitrary manner to illustrate heads and flows at various junctions. To print out all the results would take excessive computer time. The maximum and minimum HGL elevations, however, are printed out for all junctions. Although calculations are made for each 0.1 sec the results are shown for 1.0 sec intervals. The core storage has been moved back after each ten time steps. This, however, is not necessary with the small number of pipes used. The program, with the dimension declarations shown (statement 025) is suitable for a 250-pipe

-16system with. up to 20 valve-controlled outletso An approximate formula for JT, the number of steps that may be retained in storage is JT 185=00 - (14JP + 5 JU + 17 JV) (22) 2 JP + JU + JV in which JP is the number of pipes, JU the number of junctions and JV the number of valves. This is based on the 32K - core storage of the IBM 7090o If more detailed information is desired on the transient conditions within a given pipe, any section between two reaches (or all sections) may be considered to be a valveless junction. In fact, column separation and subsequent rejoining could be calculated for the pipe. Computer execution time for larger systems would be roughly proportional to the number of pipes in the system. The example problem required 81 sec time, from which estimates may be made for larger systems and longer duration of the transient, Figure 9 shows the transient heads resulting from a pipe rupture at junction 21o For this calculation the head at junction 21 was set equal to 0. after steady-state conditions were established with the reservoir at 200. feet. Minimum heads were restricted to - 335 feeto Since these periods of vapor formation lasted only a short time, up to 1 2 seco, very little column separation would occur and high pressures would not be obtained.

-17SUMMARY AND CONCLUSIONS A procedure has been developed for calculating hydraulic transients in complex distribution systems containing up to two or three hundred pipes with multiple valved outlets~ The complete computer program, in compact MAD compiler language is presented. The basis of the waterhammer analysis is the use of algebraic waterhammer equations applied over each pipe of the system, thereby calculating flow at each end of the pipe and head at each junction at time increments of the order of one-tenth the wave travel time in the longest pipe. REFERENCES 1. "Water-hammer Analysis Including Fluid Friction", by V.Lo Streeter and Chintu Lai, Transactions ASCE, Vol. 128, 1963, Part I, pp, 1491 - 1552. 2. "Analysis of Flow in Networks of Conduits or Conductors," by Hardy Cross, University Illinois, Bullo 286, November, 1946.

-18Figure 1. Simple circuit to which head balancing is applied. ___Q - I.. 0 0 Figure 2. Pipe of length L and wave speed a, with positive flow direction shown by arrow. 1 ~~17 91 F. n 9 np aa2i Figure 3. Junction 91 having four pipes and a valve outlet.

4% 22 35~~~~~~~3 1~~~ 1rv3 \ T 2 12 39 25 33 37~~~~~~74 3'11 // 19 26/ \ \\ /42 L. 91'9 0: \\28 / ^ 55 8 9~~~~~~~~~~~~~~~~~3 4~~~~~~4 Figure 4, Example system comprised of 5 pipes, 26 junc~tjons, 20 simple circuits and 10 valved outlets. Junction numbers are encircled. Assumed positive flow direction shown by arrows.

-20HERE READ AND PRINT DATA *001 (1=4,l, I.G.SG( I GII-l)+2*( l-l)*N( l-l). __ *002 SPRAY.(NN,NN(L)...NN(JdP) *003 CALCULATION OF STEADY STATE HEADS AND FLOWS (M=Ol,,M.G.Z.OR.QD.L.TOLtQD~sO.~(L3,lL.G.N ~UZ2*L,~VN(L*U,( *004 2 I=1,U,I.G.V,DHO.,HD=O.,I(J-Ol,J.E.L~,W-XG(LIIl*J*L,Y-X(G(L) *004 3 +I+J),QP-R(W) *.ABS.QQ(W).P.(NN(W-l.),DH-DH+S(Y)*QQIWI*QPtHD *004 4 HD NN( ) *QP),DQDH/HO, QD=QD+.ABS.DsIJaLi1.j uA.UlaXlG(L I+J *Q004 5 -L ),Y=X(G(L )++J),QQ(Y)=QQI(Y-SIWI*DQ~ ))) *004 (I=1,,lI.G.JP,R(I) )R(I)*SSQQ(IIQQ(I)*SF) *0r05 (I=1,3, I.G.JXX,L=XX(1+2),Y=S(XX(I ),lZ=XX(I), *006 2 HH(L)=H-Y*R(Z)*QQIZ)*.ABS.QQ(Zh.P.(NN(Zi-l.,H-HH(L)I *006 PRINT RESULTSM,TMSS,SFTOL,QQ(1)......QQ(JP,R(...R(JPI,HH(l)...HH(JUI *007 READ AND PRINT DATA *008 (1=1,1, I.E.JU,HMIN(I)H,HMAX( 11=0.) *009 (I=3,1,I.G.6,GIl)=G(I-1*+N(I-ll*(3*I-1~I *010 (1=1,1, I.G.JV,DTA(I=TC(I)/DTV,VCN(I )=QVV( I)*SF/SQRT.(HH(I)) *011 (J2=1,1,J2.G.JP,(Jl=1,1JIl.G.2,(J=1,1,J.G.JJtQIJ1,J2,J =QQ(J2 *012 2 ) l) *('12 (J=1,1,J.G.JJ,(Jl=1,1,JI.G.JU,H(Jl.J)=HH(JI)1I *013 PRINT COMMENTSl HEADS AT SELECTED JUNCTIONS FOR INDICATED *014 2 TIMES. ALSO DISCHARGES AT DOWNSTREAM END OF SELECTED PIPESS *014 PRINT COMMENTSO J TIME HJI HJ2 HJ3 HJ4 HJ5 HJ6 *015 HJ7 HJ8 HJ9 HJIO HJIL HJ13 HJ14 HJ19 HJ21 Ql *015 3 Qd Q23 I31 Q38S *015 PRINT COMMENTSO$ *016 FOX (J=J7,1,J.G.JT,T=(J+JC) *DT,(M=2,,M.G.N,J22+3*MJ3:N(M)*J2, *017 ( iK=l,J2,K.G.J3,A=O. tC=0.,E=X(G(M)+K), (L= 1tL.G.M,V=K+3*L,W=X *017 2 (G(M)+V+l),Y=X(G(Mh+V-lI,ZX(G(M)+V),OD=OQ(WtYJ-P(Y I,C(L=D* *017 3 SW +Ba(Y)*(H( Z,J-Pl(Y I-S(W)*R(Y)*D*.ABS.D.P.(NNIY)-l., A=A+C *017 4 1(L) C=C-8(Y)),AT.(T).(L=I,I,L.G.M,VK+3*LtFK(XIG(NMI+V*1,Y *017 5 =X(G(M)+V-11,Q(FY,J)=S(FI*(BIYi* HIEJ-Cl(L1ll lPTI.(T I *017 (J=l,1,J.G.JJ,(V=1,1,V.G.JU,H(V,J)~H(V,J+ZZi)I *018 (I=ilt1 I.G.JP,(V=,l,tV.G.2,IJ=l,l,J.G.JJQ(V,I.JizQ(VI,J+ZZ) *019 2 ))) 019:J=,l,tJ.G.JJ,(V=Il,.V.G.JVQV(J)(V =QV(V,JZZ)I *020 JC=JC+ZZ *r21 J7=dJJ+l.______._ ___ 022 TRANSFER TO FOX *023 INTEGER ltJtJlJ2,J3,J7,JC,UVtWtYZt,N.P.XX.EXltJRtJPJUJV,JXX, *024 2 JJ,JTZZL,,KI,K2 ML,GF,X'*024 DIMENSION X(2000) O(QVVVCNDTATCC I, K( 20),Q(2*250*20),H(135*20, *025 2 QV(20*20),(P,B,R,NN)(2501,TA(20*(O...IIIItN(9,S(2IXX(420), *025 3 HH(135),QQ(2501 (HMAX,HMIN ( 250)G(10i _____.... 025 INTERNAL FUNCTION AT.(TT) *026 WHENEVER X(G(M?+K* ).E.I.AND.TT.L.TC(E) *027 -=TT/DTA( E +1 *028 TH=(TT-I*DTA( E)/DTA( E *029 CV=VCN(E)*(TA(E, I+.5*TH*(TA(E, l+l-TAIE, I-liTH*(TA(EE, +l+T *030 2 A(E,I-1)-2.*TA(EsI I)) *030 WHENEVER CV.LE.O.(TRANSFER TO FXI *031 HIE,J)=(.5*(CV/C).P.2-A./Ci*l.-SQRT.Il.-IA/(C*(.5*(CV/C.P.2- *032 2 A/C))).P.21) *032 QV(E,J)=CV*SQRT.(H(ErJII.033 OTHERWISE *034 FXI WHENEVER XlGtL3tKtLJi.l.Q ltULE.JIJ___.;10)5 H(E,Ji=-A/C *036 END OF CONDITIONAL *037 WHENEVER H(E,J).L.HMIN,(E)HMINIE(EH(E,J) *038 WHENEVER H(EJ)I.G.HtMAX(E),HMAX(E)=HIE*JI *039 FUNCTION RETURN *040 END OF FUNCTIOSN.___- ___ _.___. - *041 INTERJNAL FUNCTION PT.(TT) *042 H(JUJ)=HO *043 ( =, l,I.G.2,O=Q(2t I,J-P( )Q( l, I,J)D+BI )* (HO-H 13-I,J-P(I *044 2 ))-R(I)*D*.ABS.O.P.(NN(I)-l. I) _044 WHENEVER J/U*U.E.JPRINT FORMATSIH,I3,F6.61lSF6.0,5F6.2*$,J+JC,T,(V=1,1,V. *045 2 G. 1,H(WVeJl)Hi3JeH"4!,JjIHJig9.J eH(2 1J, 12_.l.,1Ql(2eR.*J.*045 3 ),Q(2,23,J),(2,31,J),Q(2,38,J) *045 WHENEVER T.G.TNM *046 PRINT RESULTS HMAXII)...HMAX(JU-1)tHMIN(1)...HMINHMI U-I) DTDTV,JTJXX *,47 TRANSFER TO HERE.*. r.48 END OF CONDITIONAL *049 FUNCTION RETURN - ____-_____ __ 050 END OF FUNCTION *051 END OF PROGRAM ___ _ __*__052 Figure 5. Computer program for balancing heads by Hardy Cross Method and for calculating transients by algebraic waterhammer equations.

-21INPUT DATA FOR STEADY STATE FLOW AND HEAD DETERMINATION J(1t)- 60.0 40.0 10.0 30.0 30.0 10.0 20.0 15.0 15.0 35.0 15.0 30.0 20.0 22.0 12.0 10.0 30.0 30.0 5.0 16.T 4.0 10.A 6.P 11.0 20.0 10.0 25.0 15.0 25.0.0 15.0 8.0 3.0.0 6.0 9.0.0 6.0 9.0 6.0 -5.0.C 10.r 25.0.0 R(lI) 1.0 1.4 15.0 7.0 6.0 8.0 5.0 6.0 6.5 3.0 9.0 4.0 6.0 5.5 10.0 8.0 4.5 4.0 20.0 5.0 25.0 7.0 In. 8.0 3.C 7.0 2.5 5.0 3.2 20.0 7.0 12.0 25.0 22.0 15.0 12.0 18.0 14.0 9.0 8.0 10.0 8.0 7.5 2.3 6.5 xlL)= 1 2 2 12 11 3 1 2 2 4 6 5 1 1 2 22 34 23 1 1 2 24 30 25 I 2 2 25 26 19 1 1 2?6 28 27 I 2 2 35 33 34 1 1 2 33 37 32 1 2 2 45 44 42 1 1 2 2 2 3 4 1 1 1 2 2 13 21 15 14 1 1 2 15 20 17 16 I 1 2 2 17 19 18 9 1 2 2 2 23 24 20 21 I 1 2 2 32 41 31 30 1 1 2 2 31 43 29 28 1 1 2 2 3ti 42 43 41 1 1 1 2 11 14 16 10 1 I 1 2 2 6 10 9 8 7 1 1 2 2 2 36 39 40 38 37 XX( L) 1 1 12 5 1 14 7 1 10 8 1 15 18 1 20 27 1 9 29 1 8 44 1 7 45 2 25 40 1 6 39 2 24 36 2 23 37 1 5 41 1 22 31 2 21 30 2 4 34 2 2 22 2 18 21 1 19 20 1 3 17 2 16 16 2 17 14 2 1 11 2 13 3 2 11 2 2 26 N(3)9,9,2,NN=2.00~Z16,N=5,TGL=.10,JP-45~S(I 1l=S(2~)=-1JXX=78,HH(26)=500,SS=.nl.SF=.60,JU=26,H-500G(3)=0,QOD2 Figure 6. Input data for first read statement. INPUT DATA FDOR ALGEBRAIC WATERHAMMER SOLUTION P(l) 7475o 66 56 55565 5574565 55696765679665 5 7675547699 (1l)-.04).04.040.C12.020.020.039.016.012.OC8.024.010.027.020.014.012.026.021.020.033.01.l01.011.O 10.00d.014.006.019.010.016.015.021.008.004.002.006.007.003.006.006.008.001.004.009.013.006 rC(I) 7.0 5.0 6.0 16.0 8.0 9.0 15.0 12.0 10.0 10.0 QVV(I)- 3.0 4.0 10.0 6.0 7.0 15.0 25.0 10.0 15.0 5.0 TA(1,0)l.00.90.80.70.60.50.40.30.20.10.00.00 TA(2,0)1-.00.93.85.76.66.55.43.32.20.07.00.00 TA(3,0-)1.00.91.82.73.60.48.37.27.15.04.00.00 TA14,0-l1.00.80.62.46.32.20.10.05.01.00.00.00 7A(5,0)-1.00.90.80.70.60.50.40.30.20.10.CO.00 TA(6,0)=1.00.93.85.76.66.55,43.32.20.07.00.00 TA(7l/O)1.00.91.82.73.60.48.37.27.15.04.00.00 TAI(8,O)-l.00.80.62.46.32.20.10.05.01.00.00.00 TA19,0)z1.00.93.85.76.66.55.43.32.20.15.10.00 TA(I0,0)-l.00.90.80.70.60.50.40.30.20.10.00.00 2 PIPES Xll1 10 1 7 14 1 8 15 2 24 0 36 23 1 39 6 2 6 1 39 24 1 40 25 1 7 1 44 8 1 45 25 1 3 PIPES 11 0 2 26 1 3 13 2 12 1 2 12 0 1 26 1 4 13 2 5 14 2 14 0 5 12 1 6 13 2 7 10 2 15 0 8 10 1 9 16 I 18 20 2 17 0 14 1 1 15 19 2 16 16 2 19 0 21 18 15 I17 1 20 3 2 9 1 27 20 1 28 21 1 29 8 2 2 1 22 18 1 34 4 2 35 23 2 22 0 41 5 1 31 21 1 43 8 2 4 PIPES I I 12 11 I 11 13 1 14 17 2 13 18 2 16 0 10 13 1 16 17 1 9 15 2 17 3 2 20 0 19 3 1 18 15 1 26 21 2 27 9 2 18 0 13 1 1 21 19 2 23 4 2 22 2 2 23 0 35 2 1 33 4 1 36 24 2 37 5 2 5 1 37 23 1 32 4 1 41 22 2 38 25 2 8 1 29 9 1 43 22 1 42 25 1 44 7 2 25 0 38 5 1 42 8 2 45 7 2 40 6 2 5 PIPES 13 0 4 12 1 3 11 1 6 14 1 11 1 2 10 16 2 3 1 17 16 1 20 19 1 19 20 2 24 4 2 25 21 2 21 0 25 3 1 26 20 1 30 4 1 31 22 2 28 9 2 6 PIPES 4 1 23 18.1 34 2 1 24 3 1 30 21 2 32 5 2 33 23 2 K( I)2,K(2)lU-1O0,TM-30,DT-.l,NI6,N(21~4,9,8,3,1,JVtlO,JT-20,ZZ=10,J7-10,HO'500,G(21=0,DTV10JJ-, 10,JC:-10 Figure 7. Input data for second read statement.

HEADS AT SELECTED JUNCTIONS FOR INDICATED TIMES, ALSO DISCHARGES AT DOWNSTREAM END OF SELECTED PIPES 4 TIME HJ1 HJ2 HJ3 HJ4 H1 5 HJ6 HJ7 HJ8 HJ9 HJ10 HJ11 HJi13 HJI14 HJ 19 HJ21 Ql Q8 Q23 Q31 Q38 0.0 471 451 454 451 447 441 441 443 448 468 488 475 476 457 451 30.39 9.84 7.47 7.67 5.46 10 1*0 474 469 463 464 472 472 468 469 464 478 489 476 478 461 455 30.39 9.92 7.34 7. 55 5.44 20 2.0 486 499 481 491 507 522 520 509 494 484 491 482 482 476 481 30.29 9.89 7.20 7.15 5.39 -30 3-. — 504 540 519 532 562 600 575 556 536 503 499 495 493 507 529 29.95 9,69 6.85 6.78 5.24 40 4.0 531 592 564 593 625 676 632 618 597 533 508 518 516 553 586 29.05 9.21 6.28 6.41 4.99 50 5.0 567 645 624 647 695 764 122 686 661 578 523 549 547 598 644 27.34 8.62 5.57 6.00 4 78 60 6.0 598 666 673 706 761 856 830 770 734 625 536 585 583 658 709 24.57 7.75 4.68 5.39 4.35 70 7.0 621 720 719 753 833 952 931 861 805 664 544 614 616 693 770 20.89 6.89 3,34 4.59 3. 79 80 8.0 642 768 757 813 891 1052 1025 926 864 691 554 629 630 719 828 16.58 5.86 1.94 3.70 3.05 so 9.0 65L7 815 794 853- -921 1012 -10L8 966 883 &99 5571 633 _630 755 864 12.16 4*60.21 2.77 2.18 100 10.0 675 829 803 861 938 980 1086 968 901 701 565 641 633 779 870 7.74 3.11 -1.42 1. 89 1.59 11O0 11.0 -671 774 792 828 845 872 1029 950 897 697 566 650 643 751 866 3.12 1.52 -2.77I 1.11 1.17 120 12.0 628 695 764 743 772 790 959 842 803 705 549 640 645 744 794 -1.58.19 -3.99.75 60 30 13.0 _S96 636 677 679 108 743 822 757 730 653 538 608 620 672 701 -5.85 -.37 -5.25.54.24 140 14.0 568 625 586 609 624 675 690 635 614 562 529 555 550 559 587 -8.53 -.94 -5.56 -.24 -. 16 150 1'0?12 S4S 8n 527 557 554 5?22 -500 478 __483 500 478 478 503 508 -8.90 -1.15 -5.48 -.37 -. 20 160 16.0 441 425 409 433 421 415 357 394 403 391 477 430 419 404 405 -7.18 -1.05 -5.35 -.29 -.01 ___10 17.a__ -395 329 337 303 279 251 288 290 330 376 458 399 399 338 322 -4.21 -1.23 -4.85 -.07 04 180 18.0 369 247 282 240 184 146 203 220 232 383 449 395 408 322 224 -.89 -.10 -3.98.08.16 — I90 19.0 372 247 254 206 173 119 i 128 146 201 329 457 405 404 274 210 2.03 1.15 -2.80.21.30 200 20.0 376 288 257 245 195 171 90 144 179 343 453 384 394 281 203 5.46 2.32 -1.51.77.62 210 21-0 381 33- 291 Z 27 2 101 2 5 192 - 233-3630_A52 374 373 320 252 9.31 3.64 -.35 1.11 1.07 220 22.0 409 341 346 340 320 323 326 319 331 336 469 404 390 345 331 12.96 3.94.38 1.39 1.23 -230 23.0- 437 407 411 403 415 427 471 462 430 439 474 443 449 423 422 15.36 4.01.74 1.20 1.07 240 24.0 491 498 503 506 517 524 567 541 527 521 494 499 509 492 517 15.99 4.32.58.97.98 250 25.0 565 600 613 622 632 652 599 616 622 567 528 560 556 591 616 14.82 4.07.27.75.74 260 26.0 624 704 683 698 720 715 678 682 681 637 547 603 597 697 694 11.97 3.29 -.17.85.92 27l0 27. 0 67 i727 730 756 752 769 -81 7 777 754 _674 557 - 633 632 696 759 8.16 2.63 -.65 1.06 1.05 280 28.0 622 710 735 751 792 803 854 842 824 684 549 635 641 710 788 3.59 1.76 -1.60.80.70 290 29.0 607 685 732 728 771 848 886 844 773 664 537 608 613 714 768 -.34.84 -3.04.32 48 300 30.0 609 688 683 717 750 795 819 778 756 611 541 582 570 660 749 -3.43 -.23 -4.17 03 -.12 -- -— _ — ___ _ H AX111...HIMAXL251 - -. 6.793659E 02 8.374125E 02 8.031355E 02 8.611248E 02 9.460591E 02 1.078578E 03 1.092752E 03 9.740511E 02 9.284932E 02 7.085588E 02 5.6905878E 02 5.860337E 02 6.504819E 02 6.472526E 02 7.667943E 02 7.247586E 02 7.241734E 02 7.648816E 02 7.793333E 02 8.328438E 02 8.785914E 02 9.280301E 02 9.269794E 02 1.032449E 03 1.049448E 03 HMlN(1)...HMINf25J 3.6795378 02 2.351517E 02 2.497660E 02 2.063129E 02 1.686469E 02 1.189392E 02 8.494810E 01 1.250823E 02 1.785323E 02 3.248854E 02 4.485489E 02 4.194394E 02 3.716610E 02 3.725363E 02 2.670132E 02 3.165466E 02 -3.183864E -02 2.92529BE 02 2.681525E 02 2.352643E 02 1.987157E D2 1.505149E 02 1.651205E 02 1.191523E 02 1.279425E 02 Figure 8. Transient solution for problem with data given in Figures 6 and 7.

HEADS AT EACH JUNCTION FOR RUPTURE AT JUNCTION NUMBER JI TIME Jl J2 J3 J4 J5 J 6 J7 J8 J9 JO1 Jll J12 J13 J14 J15 J16 J17 J18 J19 J20 J21 J22 J23 J24 J25.00 171 151 154 151 147 141 141 143 148 168 188 191 175 176 162 164 165 157 157 152 0 147 148 144 142.50 171 151 154 151 1471 141 141 143 109 168 188 191 175 176 162 164 165 157 157 152 0 147 148 144 142 1.00 171 151 115 77 147 141 141 143 109 168 188 191 175 176 162 164 165 157 157 129 9 -33 148 144 142 1.50 171 141 90 68 100 141 141 76 73 168 188 191 175 176 139 144 144 120 121 78 0 -33 119 144 142 2.00 145 98 67 92 80 141 117 71 57 153 188 191 166 176 81 127 128 110 99 82 0 20 93 114 97 2.50 131 87 76 60 96 118 116 57 41 116 173 186 155 160 92 93 100 107 71 44 0 67 80 102 83 3.00 118 90 40 61 73 111 104 85 42 118 158 173 136 139 73 93 69 91 64 47 0 41 79 85 89 3.50 102 81 40 46 70 100 98 63 38 106 163 158 108 Ill 71 73 89 77 42 44 0 -8 69 73 101 4.00 110 64 64 41 68 98 94 54 52 89 158 136 92 90 87 67 77 63 64 43 0 -24 68 91 105 4_.50 99 56 39 49 58 99 93 51 30 76 164 135 96 94 64 83 66 77 86 57 0 -18 60 98 73 5.00 110 67 50 43 55 90 80 55 39 75 153 147 107 92 51 62 108 76 43 29 0 51 65 79 56 5.50 119 70 38 45 54 81 74 57 27 65 165 170 114 115 47 89 74 78 64 37 0 44 62 68 61 6.00 108 68 39 50 52 76 71 51 33 90 168 175 136 129 63 92 72 85 67 41 0 22 54 71 71 1 6.50 128 69 51 43 54 75 72 43 31 100 161 173 149 149 81 83 93 81 50 33 0 -28 64 56 80 s 7.00 118 70 40 45 51 72 66 38 23 113 179 180 145 161 69 111 96 88 81 51 0 21 62 67 58 7.,50 124 74 65 46 46 68 61 37 29 126 172 185 144 148 97 98 124 86 75 53 0 31 67 69 38 8.00 134 75 65 50 51 62 57 39 39 130 166 178 147 143 112 110 120 94 83 62 0 45 67 60 53 8.50 131 80 61 57 50 62 54 46 34 129 182 176 146 143 92 121 103 93 99 74 0 26 63 70 58 9.00 134 78 73 53 54 62 55 34 41 115 172 168 135 139 93 99 119 94 82 57 0 -16 62 59 69 HMINl )...H IN (I25 _ 9.860570E 01 5.591483E 01 2.799770E 01 3.828903E 01 4.432028E 01 6.004538E 01 5.405500E 01 2.899357E 01 2.109700E 01 6.113739E 01 1.525240E 02 1.260319E 02 9.209563E 01 8.483140E 01 4.031638E 01 5.863247E 01 5.948784E 01 6.307012E 01 4.021127E 01 2.046601E 01.OOOOOOE 00 -3.300000E 01 5.025888E 01 4.737748E 01 3.773393E 01 Figure 9. Transients caused by rupture at junction 21 when reservoir is at 200 feet.