THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING A STUDY OF WATERHAMMER INCLUDING EFFECT OF HYDRAULIC LOSSES Chintu Lai A dissertation submitted in partial fulfillment of the requirement for the degree of Doctor of Philosophy in the University of Michigan Department of Civil Engineering 1961 December, 1961 IP-546

Doctoral Committee: Professor Victor L. Streeter, Chairman Professor Ernest Fo Brater Assistant Professor Bernard A. Galler Associate Professor Arthur G. Hansen Professor Chia-Shun Yih

ACKNOWLEDGMENTS The author wishes to express his sincere appreciation to Professor Victor L. Streeter, chairman of the commit-tee, for his suggestion of the topic, and for his unfailing encouragement, advice and guidance during the course of this study. The author also wishes to express his gratitude to each of the other committee members for their interest and cooperation in supervising his work. Thanks are extended to the many people with whom he has been associated while doing this work, and by whom his work has been greatly accelerated. The research was partly supported by the Civil Engineering Department of the University of Michigan. The author also received a fellowship and a scholarship from the University of Michigan, which enabled him to pursue his graduate study with no financial worries. He is deeply obliged to these sources of support. The work would have been impossible without the use of the facilities of the computing center of the University of Michigan. The author gratefully acknowledges this convenience rendered by the center. The aid of the Industry Program of the College of Engineering in the final preparation of the manuscript is also much appreciated. ii

TABLE OF CONTENTS Page ACKNOWLEDGMENTS oo.............................................. ii LIST OF TABLES......o........................................ v LIST OF FIGURESo.o...............................O..O.......O vi LIST OF PLATESo.......... ooo o.......... eo.................... viii NOMENCLATURE O O o o o o o o o o o o o o o o o o o o o o o ix Io INTRODUCTIONoooo. 0 o 0 e o................. oo o.......... 1 IIo THE WATER HAMMER EQUATIONS INCLUDING FRICTION EFFECTooooo 4 Ao Condition of Dynamic Equilibrium..o............. B. Condition of Continuity.O.......................... 8 CO SummaryO.... o o......... o...... o...o... o.... o..... o 9 IIIo A METHOD FOR THE SOLUTION OF WATER HAMMER EQUATIONS BY DIGITAL COMPUTERoo..oo.o...o.oooooo......... 11 Ao Characteristic Equations for Water Hammero o......... 11 Bo Finite Difference Approximations 0........0o.......... 13 Co Specified Time Intervalso....... o..................o 15 Do Boundary Conditionsoo....oo.......... o..... 18 Eo Extrapolation Procedures o.... ooo... o.... o....... 21 Fo The Procedure in Machine Computation and Flow DiagramO o 0..o... oo.. o o.....o... o.. 22 Go Further Remarks..o........ o..o......... 0.....oo....o. 24 IVo APPLICATIONS TO VARIOUS PIPE CONNECTIONS AND DIFFERENT GATE CLOSURES..ooo....o.... ooooo0ooooo0ooo0000 26 Ao Variable Friction Factoro............. o............ 27 Bo Rapid and Slow Gate Closures0....... o.....0.....o oo 28 CO Pipe Line with Stepwise Changes in Diameter.........o 29 Do Compound Pipes...........0000000000000000000000000. 35 Eo Nonuniform Pipeso0....o.oooo oooooo o 0000000........ 35 Vo EXPERIMENTAL SET-UP AND PROCEDURES.... o oooooooooo ooooo 38 VIo DISCUSSION OF THEORETICAL AND EXPERIMENTAL RESULTS.o...o 51 VIIo COMPARISON WITH EARLIER METHODSo....... O.O O.oo00000000000 63 VIIIo CONCLUSIONS0\ 0 0.000.000.0.ooo0o...OOOO 0 0 0o00oo 0o..00. 68 iii

TABLE OF CONTENTS (CONT'D) Page APPENDIX I - A MAD LANGUAGE PROGRAM FOR CASE I................... 71 APPENDIX II - A MAD LANGUAGE PROGRAM FOR CASE IIo. 0o....0..... 78 APPENDIX III - A MAD LANGUAGE PROGRAM FOR CASE IIIo00 o0o0 00 0 0o o0 86 APPENDIX IV - A MAD LANGUAGE PROGRAM FOR THE MOODY DIAGRAMo....o o 92 APPENDIX V - THE TABLE OF SYMBOLS USED IN THE MAD STATEMENTSo.... 94 BIBLIOGRAPHYo............ 0 0...................... oo 95 iv

LIST OF TABLES Table Page I Experimental Data...................................... 52 II Symbols Used in the MAD Statements Through Appendices I -- IV and Their Equivalent Conventional Notations or Descriptions.......................... 94

LIST OF FIGURES Figure Page 1 Definition Sketch................................ 6 2 An Intersection of Two Characteristics................ 17 3 The Specified Time Interval........................... 17 4 The Procedure of Computation......................... 17 5 The Right-End Boundary................................ 20 6 The Left-End Boundary................................. 20 7 Abbreviated Flow Diagram for the Solution of Water Hammer Including Friction Effect in a Simple Pipe Line of Constant Cross Section and Constant Wall Thickness............................................. 23 8 Valve Closure Time Relations...................... 30 9 Variation of Water Hammer Pressure with Time, at Valve; End, Rapid Closure, Tc _.2L/a, for Different Valve Closure Time Relations...................... 31 10 Variation of Water Hammer Pressure with Time, at Valve End, Rapid Closure, 0 < Tc 2L/a, for Valve Closure Time Relation 1........................ 32 11 Variation of Water Hammer Pressure with Time, at Valve End, Instantaneous Closure, Tc = 0............. 32 12 Variation of Water Hammer Pressure with Time, at Valve End, Slow Closure, Tc> 2L/a, for Valve Closure Time Relation 1...O.........., o............o 33 13 Pipe Line with Stepwise Change in Diameter............ 34 14 A Junction Where Three Pipes Meet................ 34 15 Three Examples of Compound Complex Pipes.............. 36 16 Resistance Diagram for 1/2 in. O.D. Copper Tube....... 46 17 Resistance Diagram for 3/8 in. 0.D. Copper Tube....... 47 18 Pressure-Time Curve at P1; Case I; (No. 1),........... 56 vi

LIST OF FIGURES (CoNT'D) Figure Page 19 Pressure-Time Curve at Pi; Case I; (No. 2)...o......... 57 20 Pressure-Time Curve at P2; Case II; (No. 3)............ 58 21 Pressure-Time Curve at P3; Case III; (Noo 4)........... 59 22 Pressure-Time Curve at P1; Case I; (No. 5)............. 60 vii

LIST OF PLATES Plate Page I The Experimental Set-Up............................... 39 II The Close-Up of the Solenoid Valve and the Pressure Pick-Up Element.............................. 42 III The Calibration of Pressure Reading.................. 43 IV The Experimental Results (No. 1 and Noo 2)........ 53 V The Experimental Results (No. 3 and No. 4)............. 54 VI The Experimental Results (No. 5)....................... 55 viii

NOMENCLATURE Symbol Units Description A ft2 Cross sectional area of pipe a ft/sec Velocity of pressure wave b ft Thickness of pipe wall c1 Factor of pipe restriction D ft Inside diameter of pipe E lb/ft2 Modulus of elasticity of pipe wall material f Friction factor g ft/sec2 Acceleration of gravity H ft Piezometric head at a given time and place on the pipe; H = P + Z Ho ft Head at the reservoir hf ft Friction loss in feet K lb/ft2 Bulk modulus of elasticity of the liquid K,KX Minor loss coefficient L ft Total length of pipe Ln ft Length of any uniform section of pipe P lb/ft2 Pressure intensity at any point in a pipe R Reynolds number t sec Time Tc sec Time of gate closure V ft/sec Velocity in pipe for surge conditions Vo ft/sec Velocity in pipe for initial steady conditions ix

NOMENCLATURE*t (CONT'D) Symbol Units Description x ft Distance measured positive from upstream Z ft Elevation Angle of slope of pipe Y lb/ft3 Specific weight of the liquid ft Roughness height p. Poisson's ratio for the pipe wall material p.~ lb-sec/ft2 Viscosity V ft2/sec Kinematic viscosity p slug/ft3 Density of the liquid T The ratio of effective gate opening to the full gate opening See also Figure 1. t Somewhat different notations are used in the examples of "MAD" language programs because there are only twenty-six Roman capital letters available in the keypunch. They will not cause any confusion since they appear only in the appendices~

Io INTRODUCTION Not too many problems in water hammer including friction losses have been solved, partly owing to the comparatively small, insignificant quantity of friction loss in this phenomenon, where rapid changes in flow are most likely to occur, and partly owing to the fact that no entirely satisfactory method has been devised for their inclusion in water hammer analysis.(13,l5,16) However, this effect sometimes cannot be disregarded in cases of slow gate closure, a long pipe, or a high friction factor. The existing methods for estimating friction losses are inadequate, owing to some rough approximations. A few authors have developed graphical methods,(21l5) whereby the effect of losses can be approximated by placing one or more hypothetical obstructions at certain selected locations along the pipe line and by lumping the friction losses of each section at these points. These methods will produce an abrupt drop of pressure at each obstruction, which evidently is not the true case. Although the accuracy can be augmented by increasing the number of obstructions, this increases greatly the complexity of the graphical solution. Some analytical solutions for this class of study have been published also. All of them have employed a linear approximation for friction effect which is different from the actual situation, They generally involve difficult mathematical operations or tedious series.(l7,24) It is, viewed as necessary and useful to devise a better method which can be used to solve the problems in water hammer including friction losses more directly and more rapidly, with closer agreement to the -1

-2actual situation. An attempt has been made by the author to attack the problem directly from the basic partial differential equations, without resort to graphical methods, indirect mathematical transformations, or linear approximations. In recent years, the coming into existence of electronic digital computers has made accurate calculations describing many complex physical phenomena practical in cases formerly beyond reach. The numerical solutions by digital computer will play an important role in the analytical part of the present study. In following chapters, the basic partial differential equations for water hammer including friction losses are reviewed first. The direct solution of these nonlinear equations with aid of the computer is then presented with mathematical analysis, computer approach, flow diagram and some notes in programming, The method of characteristics has been employed in the present study. The next part is devoted to a study of the application of the theoretical work. This is necessary for better description of the computer method, for appraisal of the validity of the present method, for comparison with the earlier methods and for the guidance in obtaining experimental verification. The experimental equipment was built in the Fluid Engineering Laboratory at the University of Michigan for the investigation of the attenuation of pressure surges due to friction losses. The experiments were performed for flow in a long pipe line and their results were compared with the theoretical solutions. The discussion of these results is given in Chapter VI. This study ends with a comparison with earlier methods, conclusions, and brief suggestions for further investigations.

-3The analysis has been started from the assumptions of homogeneous and elastic liquid and pipe, uniform velocity and pressure over any cross section, and sufficient minimum pressure to exceed the vapor pressure. It has been assumed in this study that the effect of hysteresis is negligibly small in comparison with the wall friction. By applying the principle of energy to the problem, and comparing the original kinetic energy in the moving water column and the work done in compressing the water and stretching the pipe walls, it may be shown that during the action no energy is lost or converted into heat, if the pipe and water are perfectly elastico(14) When the stresses of pipe material are below the elastic limit, the assumption of perfect elasticity is very nearly true and the energy loss due to hysteresis of the pipe is very small~(6,12) This is at least quite true compared with the wall friction for the case of experiments made in this study. The elastic hysteresis of water, although no reliable data are available, is also considered to be very small, because the pressure changes take place in a very short time and comparatively little heat energy will be lost.

IIo THE WATER HAMMER EQUATIONS INCLUDING FRICTION EFFECT Fundamental water hammer equations are derived for a general case of variable flowo The elasticity of the pipe walls, the compressibility of the water9 and the hydraulic losses due to the pipe friction in both directions are taken into account. The method employed is generally the same as that described in Reference 15, except more terms are involved in this case. The following assumptions are made: (a) The pipe line remains full of water at all times, iLeo, the law of continuity holds. (b) The static pressure in the pipe is sufficiently high to sustain the minimum pressure above the vapor pressure. (c) The velocity of water in the axial direction of the pipe is uniform over any cross section of the pipe. (d) The pressure is uniform over any transverse cross section and is the same as that at the center line of the pipe0 (e) The water level at the reservoir remains constant during the period of investigation. For the derivation of equations, an element of water which is bounded by two parallel cross sections normal to the pipe axis is consideredo The general water hammer equations, then, can be obtained from the conditions of dynamic equilibrium and continuity~ -4

-5Ao Condition of Dynamic Equilibrium The condition of dynamic equilibrium can be set up for an element dx at a distance x along the pipe measured from the reservoir towards the valve, as illustrated by Figure lo If the crosssectional area at M is A, that at N is A + a dx, and the angle of slope of pipe is a, then the forces acting on these two faces can be expressed as 7A(H-Z) and 7(A + 6A dx)[H-Z + (H + sin a)dx] respectivelyo Here, y, H, and Z are the specific weight of water, the piezometric head and the elevation at point x. The gravitational force acting on the mass of element at its center of gravity is 1 aA y(A + 2 dx)dx 2 6x The wall resistance acting against the flow in the element is dhf feet of water. The unbalanced force acting on the element of water along the axis of the pipe is y(A + -A dx)[H-Z + (-H + sin a)dx] - 7A(H-Z) ax ax 1 A 1 A - 7(A + 7 a-d dx)dx sin ao + y(A + 2 ad dx)dhf, where the positive direction of the force is taken opposite to the direction of the normal flowo Neglecting the terms of higher order and simplifying, this force reduces to y[A aH dx + aB (H-Z) dx + Adhf]

-6HYDRAULIC GRADIENT 6 H dx FOR VALVE CLOSURE bx W.S. RESERVOIR H YA(H-Z) AI bA d+dh -(A + L d )dhx)dx 2 ubx 1 Fiur 1.Dfnto Skth DATUM PLANE THROUGH', OUTLET VALVE'Y'(A+ A dx )dx Figure 1. Definition Sketch.

-7The mass of the element of water to be moved by this unbalanced force is 1 6A 7(A + 2 ax dx)dx g dV and its deceleration is -t d in which g is the acceleration of gravity, V is the velocity of flow in pipe and t is the time0 Again, neglecting the term of higher order, this mass reduces to 7Adx Then from g Newton's second law of motion y[A6H dx + 6A (H-Z) dx + Adhf] = -Adx dV ax ax g dt Since V is a function of both x and t, a v v+ 6V dx =v + V av dt 8t ax dt 6t ax Then aH + 1 6A (H-Z) + dhf (_ _ aV), (1) ax A ax dx g t ax which is the equation of motion for the element of water. When the friction is negligible, this reduces to H + A (H-Z) =_1(a + V 6) (2) ax A ax g at ax If, A-1 A (H-Z) is negligible, as is always the case compared with other A 6x terms,(l5) Equation (1) reduces to H dhf ly V) 6H 1(6v + ( V + V (3) ax dx g at ax l6A dhf When both 1 _A and -- are negligible, then Equation (1) reduces to A ax dx the familiar form a= - (av + v g v(4) y~= g yf v

-8Since dhf is a very small quantity, the assumption dhf = f dx V2 D 2g for the element of water dx, is made. Here, f and D are the friction factor and the inside diameter of the pipe respectively. The friction factor is a function of Reynolds number R and relative roughness e/D, or f = f(R, E/D) = f(V,D,v,e/D), in which v is the kinematic viscosity and E is the roughness height of the pipe. The value of f can be obtained from the Moody diagram,(20,2l) or from experimental evaluations for a specific pipe of interest, Substituting dhf = f dx yV into (1) and (3) respectively, we D 2g obtain aH + 1 aA f v2 I(H-Z) + f V2 = 1 + V a, (5) ax A ax D 2g g at ax and aH + fV =2 1 (V V). (6) 6x D 2g g at ax Bo Condition of Continuity There will not be any change in the equation of continuity whether the wall friction is neglected or considered. Hence, the equations shown in Reference 15 will be used. aH + V aH a2 v ( at ax g ax where the velocity of pressure wave, a, can be expressed as 1 a =...... (8)

-9It is apparent from Equation (8) that the velocity of pressure wave depends on the diameter and the wall thickness of the pipe, the properties of the pipe material and of the fluid, and the ability of the pipe to move in the longitudinal directiono This indicates that a wave reflection occurs at every change in pipe thickness, area, or change in the pipe material, iLeo, a = a[b(x), D(x), E(x)]o The factor of pipe restriction, cl, takes the following values. (15) (a) c1 = 5 - for a pipe anchored at the upper end and without expansion joint, (b) c1 = 1 - 2 for a pipe anchored against longitudinal movement throughout its length, (c) Cl = 1 - for a pipe with expansion joints. 2 In these relationships, p represents the Poisson's ratio for the pipe wall material. Co Summary Summarizing the equations shown in the preceding pages, and denoting Vx - aV, Ht =, etc., we have the general differential dx = t

'-10equations for water hammer in a pipe, Hx + = 1 (Vt + W) (6) D 2g g Ht + VHX a2 Vx (7) g When the wall friction is negligible, Hx = - - (Vt + VVx) (4) g Ht + VHX -a Vx (7) It is often assumed that in Equations (4) or (6) the term VVx is small compared with Vt, and in Equation (7) the term VHx is small compared with Hto(15,25) Then Equations (6) and (7), and Equations (4) and (7) are again reduced to f V= 1 Vt Hx + _ - Vt (9) D 2g g Ht = - Vx (10) g and Hx = - g Vt (11) Ht -= g Vx (10) respectively. In these equations, a= 1 0 (8) ( K + Dcl) 9 R E-T

IIIo A METHOD FOR THE SOLUTION OF WATER HAMMER EQUATIONS BY DIGITAL COMPUTER The governing partial differential equations, presented in the previous chapter, can be solved by the,method of characteristics with the aid of a computer, using specified time intervals and an extrapolation procedure, as described in the following sections* The abbreviated flow diagram is shown at the end of this chapter, A. Characteristic Equations for Water Hammer The partial differential equations for water hammer including friction effect, (6) and (7), are rewritten in a modified form as follows: J1 = VVx + Vt + gHx 2 = (12) J2 = V + VHx + Ht = 0 (13) g t where av aH Vx' t These are two simultaneous quasi-linear partial differential equations of the first order with two independent and two dependent variables. Combining (12) and (13) linearly, J = J1 + XJ2 2 f V2= +)t = (V +X ) Vx +Vt + (g + (V)Hx + XHt + (14) * For the general and thorough mathematical treatment of hyperbolic partial differential equations by the method of characteristics, see Reference 10,o -11

-12If V = V(x,t) and H = H(x,t) are solutions to (12) and (13), then dV =LV dx +V dt, dH dx + dt, (15) or dV __ dH dx + Ht d V _ Vv, dHx + VHt (15) dt at at dt Now, by examination of Equation (14), with Equations (15') in mind, let 2 dx = dV (V + X )Vx + Vt = Vx dt + Vt dt (16) ( + V)Hx + Ht = Hx dt + Ht dH (17) then Equation (14) may be reduced to a form dtJ = dV + XdH + f V2dt = 0 (18) The conditions for Equation (14) to be in this form are, therefore, dx +xa2 dx g+ (19) dt g dt X from which V+Xa _ - + V, g X x =+, (20) -a and the characteristic equations for V and H become dtJ = dV + dH + f V2dt = O (21) a 2D dtJ = dV - dH + LV2dt = 0. (22) a 2D

-13Here, from (19) and (20), the two different characteristic directions at the point (x,t) are given by dt 1 dt 1 5+dx-V+a' dx V-a (23) Equations (21) and (22) are two separate total differential equations with t as an independent variable and V and H as two dependent variables. If V = V(x,t) and H = H(x,t) satisfy the Equations (12) and (13), then Equations (23) become two separate ordinary differential equations of the first order. These determine two families of characteristic curves, or shortly "characteristics", C+ and C_ in the (x,t) plane belonging to this solution V(x,t), H(x,t). Equations (21) to (23) can be rearranged to the following four characteristic equations, (24) dt dx (24) V+a along C+ dV + dH + V2dt = 0 (25) a 2D dt dx =0 (26) V-a along C_ dV- dH + V2dt =0 (27) a Equations (24) to (27) are of a particularly simple form and are satisfied, according to the derivation, by every solution of the original system (12) and (13). B. Finite Difference Approximations The characteristic equations for water hammer containing the friction terms, (24) to (27), may be solved by employing the first-order

-14or the second-order finite difference approximation. They are X1 f f(x)dx. f(Xo)(xl-xo), (28) xo and X1 J f(x)dx 2[f(xo) + f(xl)] (xl-xo) (29) xo respectively. The former is a linear approximation, whereas the latter uses the trapezoidal ruleo If the extrapolation procedures which is described later are employed with the linear approximation, the degree of accuracy will be increased to that obtainable from the second-order process. Inasmuch as the second-order approximation has the disadvantage of requiring some iterative method, the first-order approximation will be used here; Referring to Figure 2, C+ and C_ are two characteristics passing through R and S, and P is their intersection point. Denoting the value of x at R by xR, etc., for known values of xR, tR, VR, HR, xS, tS, VS and HS, the values of xp, tp, Vp and Hp can be found by applying the linear approximation (28) to Equations (24) to (27)o (tp - tR) (V-a) (Xp - XR) = O (3) (Vp - VR) + g (Hp - HR) + (Lf V2) (tp - tR) = 0 (31) aR t2D R (tp - ts) - (V ) (xp - xS) = 0 (32) (Vp - VS) a (Hps - HS) + ( V2S (to - ts) = 0 (33)

-15There are two most typical ways of using the set of Equations (30) to (33) to obtain an approximate numerical solution to the original set of partial differential equations; ioeo, (i) Grid of characteristics, (ii) Specified time intervals, The latter uses specified intervals in the t-direction and relates the values of V and H at the beginning of the interval to those at the end by means of (30) to (33). This method is preferred over the former in case of water hammer problems, because xp and tp are known exactly, or they will be assigned exactly, and only two values Vp and Hp are to be determined. Secondly, with method (ii) it is possible to apply extrapolation procedures* to increase the accuracy. Moreover, in the calculation of water hammer problems, method (ii) has an advantage of directly providing the velocity V and the pressure H along the pipe distance x at different times t, in tihe fiorm m-lost desirable for such studies.** Further, with method (ii) the quantities Ax, At are under the control of the individual user. This is a great benefit in treating many different pipe conditions during this study. The process of solving (30) to (33) by method (ii) is described in the following section. C. Specified Time Intervals Referring to Figure 3, ti and ti+l are the beginning and the end of the time interval At, and A, C, B are three adjacent points on the line t = ti, Ax apart from each other. Let the point P of * See po *-'See Appendices I to III.

Figure 2 fall on the intersection of t = ti+l and x = xc, and R, S on the line t = tio Two characteristics C+ and C_ pass through P, R and P, S as before. Here, the values of V and H at t = ti are assumed to be known and their values at P are to be found. The steps for the computation are as follows: (a) From Equations (30) and (32), remembering xp = Xc, tR = tC = tS, xR and xS can be readily evaluated. XR = XC - (V+a)c (tp - tc) (34) XS = XC - (V-a)c (tp - tc) (35) Here, it is assumed that (V+a)R = (V+a)c, (V-a)s = (V-a)c because of the sufficiently short distances AC and CB. (b) Using a linear interpolation, (b+)_ at =() vc-V (+)C VC-VR Ax VC-VR A~x VC-VA VC-VR = () ( +)C (VC-VA) VR-VC = - (V+a)C (X) (Vc-VA) VR = VC[1 - @(V+a)CI + VAG(V+a)C. (36) Similarly VS = VC[l + G(V-a)c] - VBg(V-a)C, (37) HR = HC[l - G(V+a)c] + HAG(V+a)c, (38) HS = HC[1 + 9(V-a)C] - HBO(V-a)C, (39) where At Ax

-17Pt d% Cx4.4- x R S ti A C B 0 X 0 X Figure 2. An Intersection of Figure 3. The Specified Time Two Characteristics. Interval. D3 D4 CI Cs C4 CXI,__ B B2 B3 B4 B$ Be, Ao AI A2 A3 A4 A5 As A7 t=to X:X0 Xo X2 X Figure 4. The Procedure of Computation.

As in (a), the assumption of sufficiently small At has been made here that C+ between P and R, and C_ between P and S are straight lines, and that (4+)R = (4+)C, (_-)S = (_-)C respectively. (c) Rewriting Equations (31) and (33) as (Vp-VR) + ac(HP-HR) + (2 V2)c(tp-tc) = 0 (40) (.Vpv) _ a(Hp-HS) + (-+ V2)C(tp-tc) = 0 (41) gc Solving (40) and (41) simultaneously, Vp and Hp are found as Vp = (VR+VS) + g(HR-HS) - (f V2)(tp-t), (42) VP = (VR+VS) + 2aC 2D C (42) 2 2aC 2D Hp = C(VR-vs) +( HR+HS) (43) With these equations, (34) to (43), the computation can be performed as follows: First V and H will be given at Ao, A1, oo., A7 in Figure 4,then using Equations (34) - (43), V and H at B1, o.,y B6 can be evaluated. In a similar manner V and H can be obtained at C2,, C4, and finally at D3, D4o The number of known points at t = to will determine how far the computation can proceed. However, if suitable boundary conditions on x = xo are given, V and H at those points marked by a black circle can be calculated. The same procedures apply to the right end. Those points are marked by a cross. It is then self-evident that once suitable boundary conditions are given at both ends, the computation can be carried out as far as desired. D. Boundary Conditions (a) The Right End: Let x = xC be the boundary line as shown in Figure 5. Then Equation (34), Equation (36) and Equation (38). can be used for calculating

-19the values xR, VR and HRo Finally, Hp can be obtained by rearranging Equation (40) as Hp = HR - aC(Vp-VR) - (_a v2)C(tp - tc), (44) when VP is obtainable from the given boundary conditions. Alternatively, Vp can be obtained from the following equation VP = VR -j( Hp HR) ) V2)C(tp-t), (45) when Hp is obtainableo If the boundary line is at the gate-end, x = xC = L, in which L is the length of the pipe. Then, Equation (34) can be written as xR = L - (V+a)c(tp-tc) (46) Other equations remain the same and Equation (44) should be used in this case, because Vp can be obtained, when t < Tc (time of closure), by solving o salV (47) g and VC -A HC+Tp (48) V Ho 0 0 simultaneously for AX and AV,(21) and b',-y ul;:i the relationship Vp = Vc -A (49) Here, Tp is the ratio of the effective gate opening at time P during the gate closure to the effective gate opening at time zero, and Ho is the head across the gate when V = Voo For t> Tc, vp o o (5so)

-20t rBOUNDARY ti+1 ti -A /R C B.&X % — X Figure 5. The Right-End Boundary. t rBOUNDARY ti+l P ti A C S\ B ax Figure 6. The Left-End Boundary.

-21(b) The Left End: Again let x = xC be the boundary line as shown in Figure 6. Equation (35), Equation (37) and Equation (39) are used to evaluate xS, VS and HS this timeo Now Vp can be found by rewriting Equation (41) as Vp = VS + g-(Hp-HS) - (f V2)c(tp-tC) (51) aC 2D if Hp is obtainable from the boundary conditions. Alternatively, Hp = HS + KVp-Vs) + (af V2)C(tp-tC), (52) if Vp is obtainable. Eo Extrapolation Procedures(l0,ll,18) Since the method of specified time intervals has been used, it is possible to employ extrapolation procedures to increase the accuracy of the computation. According to the study made by Mo Lister and L. Roberts,(ll) if a function f(x,t) is to be determined at t = 2nLt, in terms of its given value at t = 0, this may be achlieved by taking n steps of size 2At, repeating a linear process at a constant value of xo Let the value of f(x,t) thus obtained be fA(x,t)o Alternatively, we may use 2n steps of At and the value of f(x,t) thus obtained will be denoted as fB(x,t)o Now, we define f(2nAt) = 2fB(x,t) - fA(x,t) ~ (53) If T(2nAt) is computed, this value and the true value f(2nAt) agree when terms of the order (At)3 are neglected~ After n steps if nit = 0(1), then the error is of O(At)20

-22The linear procedure given by Equations (34) to (43) is of this form. Computing V and H for two different step sizes 2At and At, and combining them by using (53), the error in V and H after n steps is of O(At)2, This is the accuracy achieved when the second-order process of finite difference approximation is employed [see Equation (29)]. Fo The Procedure in Machine Computation and Flow Diagram The steps in the actual machine computation are stated briefly as follows: (a) Calculate V and H for a time interval 2At, using (34) to (43), (44) to (52). (b) Calculate V and H for At, using the same equations. (c) Using the results of step (b), calculate V and H for At againo (d) Combining the results of steps (a) and (c) by using (53), calculate extrapolated values of V and H. (e) Repeat steps (a) - (d) as far as desired. The output is the results obtained at (d). The abbreviated flow diagram for the solution of water hammer including friction effect in a simple pipe line of constant cross section and constant wall thickness (Figure 7) is shown on the following page. The emphasis is placed on a more logical presentation, omitting some machine programming details. One of the features in this basic structure of computer programming is that the same block of instructions, which is the core of the analytical solution in this study, can be used repeatedly for the

A \ |=x PRINT PRINT'L =- V H = 0, READn STAR L,Vo,HoD M e X i a Aft- t_ aa L,VoHo,D,b, 1 Dci)' b, T, v, a 9 Atl 1<- ) Q A Tc,n v c/D,, K E, c/D,,Ei / THROUGH A, x = x + L PRINT THROUGH, x A, =,L I FOR i = 01 1, Vli = Vo,Vli,Hli OR i = 0, H2i=Hoi > n Hli=Ho-igr hf io Wt 2 YES, > 3TC t -t+Atl J1M. >20L 5: NO PTHROUGH C, PRIN~~~~~~~t ~~Vi = Vli -R = 1 Hi = Hli C THROUGH D, FOR i = 1.9 11 COMPUTE D COMPUTE COMPUTE i = n VPi, HPi VPn, HPn VPontPo CrossNSectionOand Constant Wall Thickness. I~~ HOG, ~ /~oUoGo,'~ Iv. —---- ~, VU,, Y~ V oi~o ~_, L__Zi m: ~~nt=FRi,1 H FO R i = O, v = (vi )G(Vli ii > Hi = H (~i)_ili m~~~~l?=2 NO ROGHH V~i = VPi FOR i = 0 igr, f or = S ion of W IL i > n H ~ m-~~~~~~?THOUGH J, x =x + A~L PRINT FOR i = 0, i VCl B = 2(V22i)-(Vli) rx, VB, HB i > B=2H i -Hi 111Fi~gure 7. Abbreviated Flow Liagrain f or the Solution of Water Hammer Including Friction Effect in a Simple Pipe Line of Constant Cross Section and Constant Wall Thickness.

computation of V and H. By suitably storing the previous values of V and H and by using a number m as in Figure 7, the computation can be carried out ingeniously and economically to any desired length. With this flow diagram as a skelton, MAD language programs(3) can be written by isolating different ways of gate closure and the evaluation of friction factor, f = f(V,D,v,E/D), as subroutines. Some examples of MAD statements can be found in Appendices I to III. Go Further Remarks It is well known that the friction always acts against the direction of motion. That means for the reversal of flow the friction term must change its sign accordingly, Unfortunately the effect of friction is expressed in terms of V2 in the basic differential equations and is incapable of satisfying this requirement. This difficulty can be eliminated by writing V2 as VIVI which has the same magnitude as before and changes its sign automatically when the flow is reversed. This modification has been provided in the computer programming. As mentioned in the summary of the previous chapter, usually the nonlinear terms VVW and VHx may be neglected in engineering practice. These have been retained in this chapter, inasmuch as the omission of them will not mathematically reduce the difficulty of solution due to one more nonlinear term - V2 and by treating the solution more generally it may D 2g serve helpfully in some occasions when their effect comes to be appreciable. In the earlier analytical solutions, authors solved the differential equations linearly by omitting these two terms and applying a linear approximation to the friction term,

-25It is true that in the actual application of the present method to an engineering problem, some computer time can be saved by neglecting these two terms. This can be easily achieved by replacing all terms of (V+a) and (V-a) in the computer program by a and -a respectively. In case of gradually changing cross-section or wall-thickness a = a(x) must be computed at each point along x. In case of abrupt changes in the pipe line, the pipe can be divided into segments and the same water hammer computation can be performed in each segment with new boundary conditions at each junction. Some further details of computer solution and its application to a few compound systems are presented in the next chapter and Appendices.

IVo APPLICATIONS TO VARIOUS PIPE CONNECTIONS AND DIFFERENT GATE CLOSURES As briefly mentioned in the preceding chapter, with Equations (36) and (43) and the method of computation as a core, the solution of water hammer phenomena can be extended to a great many different problems, eog., to various pipe connections and different gate closures. These solutions can be found by adding modifications to the program, by providing appropriate boundary conditions, or by combining sets of computations. It is helpful and worthwhile here to introduce some applications, because: (a) The effectiveness and flexibility of this method can be best demonstrated through these applications. (b) It is useful for comparison with the results obtained by other methods. (c) Some details in the experimental verification necessitates this kind of study. Since the experiment must be conducted in conformity with the conditions studied in the theoretical analysis, when it is impossible to achieve the theoretical conditions, the theoretical solution should be modified accordingly to meet experimentally feasible conditions. -26

-27In the following parts of this chapter, several examples of solutions of water hammer problems are presented. Some of them are essential tools for the comparison of the analytical and experimental results. It is not the purpose of this chapter to attempt to cover the entire field of water hammer, but merely to indicate the usefulness of this toolo Ao Variable Friction Factor Since the friction factor f in the equation dhf = f dx V2 D 2g is a function of the Reynolds number and the relative roughness of the pipe, i.e., a function of D, V, v, and c/D, it is found more convenient to keep the evaluation of f values outside of the main program. Thus the main program can be kept neat and its structure be kept free from any influence due to change in evaluation of f. Unlike a constant friction factor used hitherto, the variable friction factor f, calculated from given D, V, v and E/D, should theoretically, be able to take care of ahy flow cases0 For the laminar flow, the value f in a pipe will be, 64 64 64v (54) R DV DV in which p and. are the density and the viscosity of the liquid. For other flow cases with different sizes and materials of pipe, the corresponding equations or values of f can be found from friction fac-tcor charts or tables.(20,21)

-28If for a specific pipe, the friction factor be evaluated in L V2 terms of hf/( V), then theoretically there should be no need for D 2g classifying the type of flow. In this study, a few subroutines for the evaluation of friction factor have been made. One is a subroutine for the Moody diagram, which employed the same equations used in producing that diagram. With this subroutine, the machine is able to provide any current value of f at any point along the pipe, from the given current data for that point. Needless to say, the usefulness of this subroutine is not limited to this study, but can be extended to any other pipe problems in connection with the computer. (See Appendix IVo ) Others are special subroutines for the particular pipes used in the experiments. They were obtained by careful determination of friction loss through those pipes, and were used in producing the theoretical results in Chapter VI, (See Appendix I.) Bo Rapid and Slow Gate Closures It was mentioned in the previous chapter that many different types of gate closures can be handled by imposing suitable boundary conditionso In this case, it is also found that the program can be kept orderly and the computation can be carried out more efficiently by providing boundary conditions for gate operation in a subroutine, because: (a) In case of an instantaneous gate closure, the main program can simply skip the step of calling this subroutine; thus the structure of the main program will remain unaffected. The same applies to the time after the gate closure is completed.

-29(b) By applying many different ways of gate closure, we can solve a number of water hammer problems produced by those closures. Since all these boundary conditions due to gate closure are supplied by subroutines and can be easily replaced without any change in the main program, great flexibility in computation, results especially in designing gate operations for keeping the pressure below a fixed maximum value. Some examples of water hammer produced by different gate closure time relations for a simple pipe line are shown in Figure 8 to Figure 10. The examples of instantaneous and slow gate closures for the same pipe are shown in Figure 11 and Figure 12. CO Pipe Line with Stepwise Changes in Diameter In the case of a pipe line with stepwise changes in diameter, the pipe line can be divided into segments at each step, and the same water hammer computation can be performed in each segment with suitable boundary conditions at each junctiono At junction A, the conditions, (see Figure 13) 2 2 (55) V1D1 = V2D2 (55) and HA= HA2' (56) in which V1, D1 and V2, D2 are velocity and pipe diameter before and after the junction respectively, may be used when the velocity head and the minor loss at the junction are negligible. An example of applying the relationships (55) and (56) can be found in Appendix II.

-301.0 0.8 0 T=(Tc-t)/Tc 0.6 \ 2 ~ T=0.5(1.5-t/Tc) -0.125 0.4 (rTc 0. 0 0.2 0.4 0.6 0.8 1.0 Tc Figure 8. Valve Closure Time Relations.

-31L: 200ft Vo =2.77 ft/sec Ho: 350 ft D =0.0365 ft b:0.00261 ft. Tc:O.09 sec v: 1.21xl1O ft/sec E:0.2448xlO lb/ft2 800 600... 400 - 200 ----------- - - - (a) VALVE CLOSURE TIME RELATION 800 _ _ _ _ _ _ _ _ _ 600 \ /I 1 i I /1 400 - 200_____________ __ - 0 ------- 1I (b) VALVE CLOSURE TIME RELATION 800 - 600 400- -. 200 c) VALVE CLOSURE TIME RELATION Figure 9 Variation of Water Hammer Pressure with Time, at Valve End, Rapid Closure, T 2L/a, for Different Valve Closure Time Relations (See Figure 8).

-32800 600 200 0 2 3 4 5 6 7 8 9 10 1 in L/a sec Figure 10. Variation of Water Hammer Pressure with Time, at Valve End, Rapid Closure, 0 < T 2L/a, for Valve Closure Time Relation iT (See Figure 8). L = 200 ft VO = 2.77 ft/sec Ho = 350 ft D = 0.0365 ft b = 0,00261 ft Tc = 0.022 sec y = 1.21 x 10-5 ft2/sec E = 0.2016 x 1010 lb/ft2 800 600 -'400 200 0 I 2 3 4 5 6 7 8 9 10 t in L/a sec Figure 11, Variation of Water Hammer Pressure with Time, at Valve End, Instantaneous Closure, Tc = 0. L = 300 ft vo== 2.77 ft/sec Ho = 350 ft D = 0.026 ft b = 0.00261 ft Tc = 0 Y = 1.21 x 10-5 ft2/sec E = 0.2448 x 1010 lb/ft2

L 200 ft. \4:2.77 ft/sec Ho=350 ft D=0.0365 ft b 0.00261ft Tc 0.2836 sec Y 1l.21x 105 ft2/sec E=0.2448 x Id' lb/ft' 500 400 300 200 0 2 4 6 8 10 12 14 16 18 20 t in L/a sec Figure 12. Variation of Water Huffner Pressure with lime, at Valve End, Slow Closure, T >2L/a, for Valve Closure Time Relation~ (See Figure 8)

-34Q A VI -DI DD2 V2 Figure 13. Pipe Line with Stepwise Change in Diameter. V2 VI DI A Di - V3 Figure 14. A Junction Where Three Pipes Meet.

-35Do Compound Pipes With the same principle mentioned in the previous section, but with more involved boundary conditions, the method of characteristics by the computer can be developed to almost any kind of pipe system. Three examples, which are the conditions employed in the following experiments, are shown here. They are: (a) a simple pipe line with a dead-end branch connection, (b) a pipe line with a stepwise change in cross section and a dead-end branch connected near th.ie valve end, (c) a pipe line with a stepwise change in cross section and a dead-end branch connected at the junction of two segments, as shown in Figure 15. The conditions to be satisfied at junction A, (see Figure 14) where three pipes meet, are VlD1 + V2D2 = V3D3 and HAl = HA2 = HA3 (58) if the velocity head and the minor loss are neglected. The detailed MAD Language Programs for cases (a), (b) and (c) are shown in the Appendices I to III. Minor losses and velocity heads have been neglected because they are very small compared with the loss due to pipe friction. Eo Nonuniform Pipes In case of gradually varying cross-section or wall-thickness, the computation introduces no major difficulty so far as the one-dimensional

-36(a) Case I (b) Case II I IIP3 ~ P4 (c) Case m Figure 15. Three Examples of Compound Complex Pipes.

-37method is applicable; i.e., the change of factors should not be so large that the assumption of uniform velocity and uniform pressure distribution at any cross section would not be valid. In those computations, instead of using constant a and D, a = a(x), D = D(x) should be used at each point along x.

Vo EXPERIMENTAL SET-UP AND PROCEDURES An experimental check was planned as the next phase in this study. In. order to make the friction effect easily observable, a pipe which would produce appreciable head loss between two ends had to be used. Although this can be achieved by increasing the length of pipe and reducing its diameter, there are some difficulties in actual installation in a laboratory. Owing to the limited space in a building, it is usually very difficult to place a sufficiently long straight pipe indoors. Moreover, a long pipe needs many supports, which will produce undesirable disturbances at each point of supporto It is also not advisable to bend the pipe, since it gives disturbances to the pressure wave alsoo Under these restrictions, two coiled copper tubes, as shown in Plate I, were finally used0 One was 300 ft. in length, 0032 in0 in thickness and 1/2" OoDo, the other one 300 fto* in length, 0032 in. in thickness and 3/8" O0Do A core, on, which the tubes were to be coiled, was built0 It consisted of two wooden plates at each end, about 3v-8" apart, and twelve 1/2" steel bars connecting the plates to frame a core. Each bar was wrapped with a soft rubber tube in order to protect the copper tubes from any damage due to direct contact with the steel bars and to prevent possible disturbance to the transmission of pressure waves at the point of direct contact0 The copper tubes were then coiled very carefully, uniformly and loosely around this core0 The diameters of the coils were about * At the end of the experiments, the former was reduced to 299o051 and the latter to 294 30, because by handling, in determination determination of friction factor, by connection and disconnection to various sources and fittings, the ruined tips had to be cut off each time~ -38

-39Plate I. The Experimental Set-Up.

two feet. It was considered that although the copper tubes were not straight but coiled, due to the uniform curvature along the entire length of the tubes, there would be no abrupt disturbance along the tubes, and they would act in a manner similar to straight tubes. The reservoir should have a capacity to provide a constant head during the experiment. A large pipe connected to a very small pipe in which the water hammer phenomena will be investigated, can serve as a reservoir. Two sources were used as intake in the experiments. One was from the 2" IoD. outlet pipe of a pressure pump which could produce about 200 psi. The other one was from a 8" pipe which was connected to a head tank, the elevation of which was about 44' above the outlet of the tube. The former has an advantage of being able to maintain a high pressure inside the tube and thus can provide a higher velocity without causing the minimum pressure to drop below the vapor pressure. However, there is a doubt of it possibly creating some transmission of pressure wave towards upstream from the junction of the tube and the outlet pipe, by which the reflection will be reduced by that amount. The latter, on the contrary, does not have such trouble, because the pipe is sufficiently large in comparison with the tube, but unfortunately the head is quite low in this case. A considerable time had been consumed for the selection of the valve. If the slow valve closure were also to be investigated, the area of valve opening had to be accurately controllable as a function of time, or at least the time and the way of valve closure should be known. Since there were no such valves in all manuals searched, and the primary interest in the present experiment was to check the analytical solution of

-41differential equations into which the nonlinear friction term was introduced, rather than to observe the difference caused by rapid and slow valve closures, it was decided to limit the test to the instantaneous closure, It was not too difficult to find a solenoid valve which could be closed practically instantaneously as compared to the time of the traveling of pressure wave; however, most of the valves were disc type and there was a possibility of bouncing of the disc during its operation. On that account, a slide-gate type solenoid valve was finally chosen, This type not only has the capability of tight instantaneous closure but also can possibly be used for the case of slow valve closure in the future if a good method is devised to slow down the speed of closure, or at least it can be used for the case of partial valve closure by providing some check. The sluice-gate type solenoid valve used in the test and its connection with the copper tube are shown in Plate II. The gate is closed when no voltage is applied, ioe,, it is a normally closed type. A hydrauliscope, as shown in Plates I and III, was used for the pressure measurements. This is a high speed electronic analyzer which can respond to pressure changes of high frequency and high rate, and is furnished with attachments to permit photographic recording of the curves traced on the screen of the cathode-ray tube, The pressure pick-up element used in connection with the hydrauliscope is of the resistance type, having a linear response with pressure variation thus permitting the direct static calibrationo It is also designed for temperature compensation. The friction factors of the experimental tubes were determined with great care~

-42N I A_ lI _lIk I Z _ c....'?i Plate TII The Clos e-U-p of the Solenoid Valve and the Pre ssure Pick-U-p Element *

i~~~~~~~~~~~~ Plate III. The Calibration::v~~~~~~~~~~~~~~~~~~~~~~~~~~ of Pressure Read~~~lin::: Plate III. The Calibration of Pressure Reading.

-44Four peizometer rings were made and attached at each end of the two tubes(8) Four radial holes were carefully drilled at a distance not too close to the end of the tube but within the reach of a small round file so that all burrs could be removed from the inside. The four orifices were placed symmetrically about the vertical bisecting plane with equal spacing, allowing no orifice to be located at the top of the tube, lest air bubbles should enter the hose connected with the manometer. Then the rings were soldered to the outside and at each set of measurements they were connected by rubber hoses to two ends of a differential manometer. Two kinds of liquid were used for the manometer. Acetylene tetrabromide, specific gravity of which is 2.94, was used for the sensitive measurements in the range of low velocity, and mercury for larger pressure drops, (22) Water was taken from a head tank located under the roof of the Fluid Engineering Laboratory. The water level in the tank was maintained constant during the measurement. The height of wreir in the tank above the tube was approximately 43 ft. Before each series of the measurements, the fastest flow of water, obtainable from the tank, was made to go through the tube to expel any possible bubbles trapped inside. The waste cocks on the top of two limbs of the differential manometer were also widely opened to bleed out the water in order to get free of any trapped air bubbles. The connecting hoses were made short and were continually sloping upward from the rings toward the waste cocks. The determination of friction factor was then made. The discharge, the temperature of water and the difference in height of the two

menisci of the manometer were measured for each different opening of the discharge cock. From these data, the velocity, the kinematic viscosity and the friction loss were evaluated and the curve of friction factor f vs. Reynolds number R was plotted. Each pair of these two values were obtained from f = hf (59) (L/D)(V2/2g) and R =DV (60) V The curves, and the MAD program made thereby, are shown in Figure 16, Figure 17, and Appendix I, respectively. Before proceeding with the experiments of water hammer, the hydrauliscope was checked and adjusted. Then it was calibrated very carefully with the aid of a dead weight gage tester, as shown in Plate IIIo The Y-sensitivity of the hydrauliscope, which is to adjust the height of the curve or to set the vertical scale at any desired value, say 200 psi for each vertical inch, was calibrated for the two pressure pick-ups used in the experiments. One of them had a 0-500 psi system operating range with 750 psi maximum continuous shock pressure and 1250 psi maximum static pressure permitted. The other one had a 0-200 psi pressure range with 300 psi and 500 psi corresponding permissible values. The pressure pick-up element was installed in the dead weight gage tester. The desired dead weights were applied on the tester and the Y-sensitivity was adjusted to give a desired vertical movement of the beam for the given weight0 This value of Y-sensitivity was noted for the experiments.

1 C 9 9 _ _ _ _ 8 7 _ 55 ~ -- t- -- __ -~ — + i-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 3.5. + ~ 2.5 2 CI. 6~~~~~~~~~~ _J _ I H 5 -8'~~~~~~ - D V — I/~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Figure 16. Resistance Diagram for 1/2 in. O.D. Copper Tube. 7 25 _ - _ -- - 1.5-_ -2 --- -....2 3 102 1.5 22.53 4 5 6 7 8910O 1.5 2 2.5 3 4 5 6 7 8 9I04 DV Ri/ Figure 16. Resistance Diagra~m for 1/2 in. O.D. Copper Tube.

-479 8 7 5 - 2.5 2 10 5 4 2.5 1.5 -2 10 102 1.5 2 2.5 3 4 5 6 7 8 910 1.5 2 2.5 3 4 5 6 7 8 9 104 R D V istance Diagram for 3/8 in. O. Figure 17. Resistance Diagram for 3/8 in. O.D. Copper Tube.

The experiments were planned for two different flow cases; one for a simple line of a constant cross-section, with the pressure pick-up attached at the valve end, the other one for a line with a stepwise change in cross-section, to which the pressure pick-up would be attached at the valve end and at the junction. However, the practical situation turned out to be the three cases shown in Figure 15 of Chapter IV, owing to the additional projection caused by attaching the pick-up element to the main line. The length of this element is so short in comparison with that of the copper tube that the friction loss caused by its wall shear and by the tee-joint (See Plate II) is quite negligible, yet the effect of storing the additional volume of water in this element can by no means be ignored. Because of the compressibility of water and the elasticity of the pick-up tube, the water stored in this part would be compressed and expanded alternatively, generating a pressure wave moving back and forth with v-ery high frequency as compared to that of the main lineo This buffering effect and disturbances caused by emitting successive tiny pressure waves in all directions from the junction soften what otherwise would be a more abrupt pressure versus time curve. This effect has been clearly demonstrated by comparing the photographs taken in the experiments and the two analytical solutions by omitting and includinglg tihe dead-end branch, (See Plate IV, Figure 11 and Figure 18o) Each flow case was investigated for both water sources, the pressure pump and the head tank. A micro-switch was connected to the "external" plug-in to facilitate the operation of the hydrauliscopeo The camera was installed oJer the screen of the hydrauliscope, as shown in Plate I, in order to take photographs of pressure versus time curves~ The camera consisted

-49only of a metal box with a fixed f 5~6 lens which could focus the screen image continuously on a standard 2-1/4 x 3-1/4 film pack, and with a black plate which could be pulled out for exposure, because there was no shutter or iris. The steps for the experiments are as follows: to start the flow, the solenoid valve was opened by pressing the switch; then the discharge cock attached at the downstream side of the valve was adjusted to obtain approximately the desired velocity in the tube. The discharge was then measured by using a stop watch and a measuring cylinder or weighing the water. At the same time, the temperature of water was also recorded. The t"X-selector" dial of the hydrauliscope was turned to the "external single" (see Plate III) so that only one sweep of the beam would appear on the screen for one push of the microswitcho To start taking a photograph the exposure plate was first pulled out of the camera, and five sweeps of beam were recorded successively to represent five distinct features. (See Plates IV to VIo) The first horizontal sweep is a pressure line at the pressure pick-up for the steady flowo Then the switch of the solenoid valve and the button of the micro-switch were pushed down almost at the same time so that the curve of the water hammer produced by valve closure could be traced out on the screen. After the pressure wave died out due to friction, the third horizontal line was recorded. This represents the static pressure in the reservoir. The zero-pressure line, or the reference line, was recorded fourth, then followed by the last line, the timing wave, which would indicate the time scale used, by providing sine waves of a definite frequency~ The screen illumination was increased before the exposure plate was inserted back into the camera, in order to

-50show the grid lines clearly on the film. When the water was taken from the pressure pump, tiny ripples indicating the pressure fluctuation caused by pumping action were observed. The experimental results thus obtained are shown in the next chapter for the comparison with the theoretical results.

VIo DISCUSSION OF THEORETICAL AND EXPERIMENTAL RESULTS A number of experiments for water hammer under different conditions were performed in the Fluid Engineering Laboratory at the University of Michigan~ Owing to the limitation of space and computer time, only several runs which best represent each flow case are listed in Table I. The photographs of pressure versus time curves taken in these runs are also shown in Plates IV, V, and VI. The frequency of timing wave on the hydrauliscope was set at 5 cps for those shown in the photographs. HoVever, a careful calibration of the timing wave, later, indicated that the actual frequency of the hydrauliscope used in the experiments was 50285 cps. From these photographs and the data given in Table I, the pressure head versus time curves were plotted by solid lines inFigure 18 to Figure 22. The same data used in the experiments were fed into the electronic computer. The factor of pipe restriction cl was taken as cl = 1 - p/2 in these computations. [See Chapter II, B., Case (c).] The velocity and the piezometric head (strictly speaking) at various points along the lines for every assigned time interval were computed and printed out in a tabular form as shown in a few examples given in Appendices I to III. Since in these experiments the datum plane was taken at the level of the pressure pick-up element and the velocity heads are so small, the total head, the piezometric head, and the pressure head at the pick-up are practically identical. From these answers, the pressure head versus time curves at the point of pressure pick-up were plotted by dotted lines in -51

TABLE I EXPERIMENTAL DATA Bulk Inside Initial Head at Kinematic Modulus of Frequency Theoretical Type Diameter Velocity Reservoir Temperature Viscosity Reynolds' Elasticity of Timing "Y" Scale Source and of Line D V Ho ~F v Numlber K Wave Factor of Experimental No. System ft. ft sec ft. ft2 /sec R lb/ft2 f psi/in Water Results Fig. 15 Fig.(a) 0.03633 15 2.94 451 106.0 0.692 15330 0.475 x 108 5.285 200 Pressure Fig. 18 Case I Pump Plate IV RD Fig. 15 1 2 (a) 0.02592 3.08 486 93.1 0.799 9990 0.466 5.285 200 Pressure Fig. 19 Case I Pump Plate IV Fig, 15 0 03633 1.112 5556 3 (b) 0.02592 2.183 4.8 0.72 0.471 5.285 00 Pressure Fig. Case TI PmIp Plate V Fig4 (c) 002365933 2.183 432 103.8 0.712 7740 0.473 5.285 200 Pressure Fig. 21 Case III Pump Plate V Fig. 15 5 (a) 0.03633 0.367 45.0 76.9 0.993 1340 o.462 5.285 40 Head Fig. 22 Case I Tank Plate VI

No. 1 No. 2 Plate IV. The Experimental Results. (See Table I and Figures 18 and 19)

-54No. 3 No. 4 Plate V. The Experimental Results. (See Table I and Figures 20 and 21)

No. 5 Plate VI. The Experimental Results. (See Table I and Figure 22.)

1000 - - NO I EXPERIMENTAL RESULTS.. THEORETICAL RESULTS -- - -- STATIC HEAD 800 - -- INITIAL HEAD AT PI 600 I /I~~~~/' I I - - - - -' - - - -~~I li I,.-A~~~I I / \''/ II / 200 I I III % I I I _ _ __\' II 20 r ~- r - -I.. -. - -. -. - -. -. 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 t, sec Figure 18. Pressure-Time Curve at P1; Case I (See Figure 15 a).

1000 -I — -I NO 2 EXPERIMENTAL RESULTS.. — THEORETICAL RESULTS -.- - STATIC HEAD 800 —80I - INITIAL HEAD AT P, I, I I I I III - I I/ I \ I I / i I 600i I 400 200 I I II 0 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 t. sec Figure 19. Pressure-Time Curve at Pl; Case I (See Figure 15a).

1000 - NO 3, EXPERIMENTAL RESULTS.....- THEORETICAL RESULTS STATIC HEAD ~~~~~~~~800 ~ — INITIAL HEAD AT P2 600, — I I~ I I. t sec Figure 20. Pressure-Time Curve at P2 Case II (See Figure b). 0 V 0 0.2 0.4 0.6 0.8 1.0 I.2 1.4 t, sec Figure 20. Pressure-Time Curve at P2; Case IL (See Figure 15b).

1000 l. -— INITIAL HEAD AT P4 600 t I I' A 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 t,SEC Figure 21. Pressure-Time Curve at P; Case III (See Figure 15c)....Fiur.2.'Pesr-ieCrea P;Cs I SeFgr 5).........

NO 5 EXPERIMENTAL RESULTS.. THEORETICAL RESULTS STATIC HEAD INITIAL HEAD AT PI 120 100 -" — I I,. I- I I ~,~~~~~~~~~~~~~~~~~~~~~~~~ / 80 4 0r~~ % 00, 60 000_ I I 0 I I I I I I i I I I 2I ( 40 F, I I I /,' 20 ~~~~~~~~~~~~~~~~~~~~~~~~~~~i~ I I I' I I I ~I \ 20 1- ___ 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 i, Sec Figure 22. Pressure-Time Curve at P1; Case I (See Figure 15a).

Figure 18 to Figure 22 for comparison. These are the theoretical solutions corresponding to the above-mentioned experiments. A few words worth mentioning here are the effect of attaching the pressure pick-up to the main line. The significance of this effect and what type of pipe system has to be employed have been described in the previous two chapters. However, a major difficulty arises in the actual computer application. Since the length of pick-up is so short compared with the length of main line, formidable computer time Will be required if the actual length of pick-up be used. This difficulty can be reduced considerably by using a much longer hypothetical dead-end branch so long as its length is sufficiently small relative to the main line, because the significant effect is caused not by the extra wall friction of the short branch, which is essentially negligible, but by its buffering effect and the emission of successive pressure waves therefrom, due to its existence, This effect should cause tiny ripples along the primary curve (this was observed in some plottings which are not shown here), but it did not appear in the theoretical curves shown in this chapter due to the scale used in plotting and the time interval used in printing. The comparison of the computed and experimental results shows good agreement for the flows with higher velocities. The computer programs, which were compiled according to the analytical procedures described in Chapters II, III, and IV, not only could depict the shape of the curve very closely to the experimental one - even for a very complicated curve -, but also gave numerical values of pressure heads very close to those obtained from the experiments.

-62The most important value for pipe design is usually the value of the first pressure peak after the valve closure. All of the computed values showed an excellent agreement with the corresponding experimental ones. They also indicate remarkable agreement in rate and way of attenuation of pressure surges. It is not uncommon in a complex or compound pipe system that the different concentrations of reflected pressure waves from different points cause the magnitude of a succeeding peak or trough to exceed that of the previous one.. The good agreements between computation and experiment in this respect will warrant the extensive application of this method to complicated pipe systems. The situation is somewhat different in flows with extremely low velocities. Observing No. 5 in Table I and Figure 22, we find a lack of resemblance in their shapes between the computation and experiment. The speed of pressure wave in the experiment is slower than that obtained by computation. The discrepancy in the case of laminar flow could be conjectured as the violation of the assumption of uniform velocity distribution over any cross section made in Chapter II. For the exact solution of such flow we may have to consider the three-dimensional case, starting from the general Navier-Stokes equations. It may also require some reexamination or refinement in the experimental technique. These are out of the scope of this study and will be left for a future investigation.

VII. COMPARISON WITH EARLIER METHODS As already mentioned in the introduction, the existing methods for estimating friction losses in water hammer are not only few but inadequate, The direct analytical solution of the basic partial differential equations is the only method capable of distributing the effect of friction losses along the pipe line in agreement with the actual situation. However, the effect of friction losses, which are assumed to vary with V2, makes the basic differential equations nonlinear even when the terms VV and Vm are neglected, as indicated in ax ax the equation of dynamic equilibriumn (6). The simultaneous solution of this equation with the equation of continuity(7), had been considered as impossible. (13,15,16,24) Thus some approximations, some of which appeared as graphical methods, some as operational mathematics and so forth, had been employed as substitutions because of no exact solution, Among those existing methods, the graphical methods are probably most widely used. The first approximation of this class is an estimation of loss by a hypothetical obstruction located either at the upstream end(l5) or at the valve end(l6) of the pipe line. This obstruction has the same total friction loss as the entire pipe line, The gist of the method is to proceed along the same fundamental network of solid lines as in the nonfrictional case of the graphical solution but adding the correction of head loss each time at the obstruction. The basic system of equations, from which the graphical method has been developed, remains the same as in the non-friction case, and the entire method is simply a combination -63

-64of nonfriction case and abrupt pressure deduction at the obstruction. Although this is a simple, good approximation for estimating pressure at a point at certain fixed times, it gives no information about other points along the pipe line. A better graphical approximation could be achieved by taking a number of hypothetical obstructions at various points along the pipe line lumping the friction losses of each section at these points. As in the first approximation, the same network must be applied, the correction of friction effect must be added at each point, and those lines would be interwoven and developed further as the time goes. In this case, it generally involves some trial-and-error method to locate two points at each side of an obstruction, except those defined by some boundary conditions. This method, as the first one, produces an abrupt drop of pressure at each obstruction which evidently is not be true case. In the graphical methods, although closer and closer approximations can be obtained by increasing the number of obstructions, this rapidly complicates the solution. This is quite evident if one considers the fact that each obstruction has two points of different heads at each side and each point of each pair emits two lines of different slope, and moreover, the distance between these two points, as well as their locations on the graph, has to be determined by trial-and-error method. These lines develop further and further into a great tangle as time progresses. Successive trial-and-error often results in the accumulation of errors and the values thus obtained may sometimes become doubtful.

-65Some authors have developed methods of estimating the effect of losses by operational mathematics or other analytical solutions (17,24,5) They differ in degree of accuracy and complexity, in method of approaching the problem, and in range of application, but they have one thing in common the linearization of the friction term, Since the effect of friction is expressed in terms of V2 in the basic differential equations, and since most cases in engineering practice are turbulent flows, the linearization of the friction term cannot exactly represent the actual situation, In other words, they are approximations and hence not wholly satisfactory. Wood(24) introduced Heaviside's operational calculus in the faield of water hammer and presented one example of a simple line with instantaneous gate closure at the lower end, allowing for friction~ He ass ured friction to be linear and for this purpose he replaced V2 by (KfVm)Vo Here, Vm is the initial velocity prevailing in the pipe and Kf is a constant for linear approximation, Besides the linear approximation, it gives only "'surge pressure"s and not the actual pressures and velocities along the pipe, The pioneer work of Wood was improved. by ich (17), superseding [!eavisidels operation by the Laplace-Mellin transformation~ It permitted working directly with total pressures and velocities instead of surge pressures and velocities, The solutions have rather involved series in Bessel functions, but the method required replacement of nonlinear tesrm as KV2, by a linear approximation0 The author also tried some Laplace transformations before the present study, and by taking advantage of the digital computer, he used several different values of Kf for corresponding v'lues of V, instead of

using a single constant Kf as was done by earlier authors. This is tantamount to dividing the parabolic curve of friction versus velocity into several segments, replacing each segment of curve by a straight line, This naturally should give better agreement with the actual situation. However, there are still some disadvantages in employing the operational mathematics, even with aid of the computer. The transformation and irnverse transformation involve many difficult mathematical manipulations and. often result in long series. It has certain limitations because of necessary approximations, and consequently is not quite flexible in application. It is sometimes very difficult or impossible to determine constants in the operational solutions for some given boundary conditions. The direct analytical solution of the basic partial differential equations, which has been presented in the previous chapters, naturally has greater accuracy over other methods because of its directness. It can handle the V2 term and consequently is able to distribute the effect of friction losses along the pipe line in accordance ^rith nature, without lumping, approximation or indirect trans formationo Besides its accuraey, this method. possesses some o ther advantages. Since the computation is carried out by the electronic digital computer, it is m:uch faster than old methodso Unlike the case in the present study, much computer time can be saved in actual applications, because there will not be a small pressure pick-up in an actual pipe, and even if the pipe system has a dead-end branch, it will be of the same order of length, Many other details, which were retained in this study- for greater accuracy, may be omitted in practical cases, For most purposes, the computer time will not be more than a few minutes to several minutes,

-67Great flexibility is one of the remarkable features of this method. With the same form of basic program, it can be developed to extensive problems of complicated pipe systems and gate operations, It may also serve in solving some problems formerly beyond scope. As can be seen in the examples given in the appendices, the anlswers can be printed out in tabular form so that it becomes extremely easy for anyone to understand the phenomena and the results are ready to plot in any form desired. It supplies all historical and geographical data for both velocity and pressure. These have all been included automatically in the course of computation, For any fixed point along the pipe, we can trace the velocity and pressure following the time variation, or for any particular instant, we can observe how velocity and pressure are distributed along the pipe, Indeed, by this method, we can gather extensive valuable information from answer sheets simultaneously, whereas the earlier methods could give answers for only limited points and limited time intervalsThe program is easily controllable both in time and in location according to accuracy and necessity, After we make up a computer program for one type of pipe system, we can solve as many as we desire for the same type of problems by simply feeding in the data, Those advantages mentioned above are also good for non-friction cases, If we prepare many programs for typical water hammer problems, a great number of solutions can be readily obtained in a form similar to data from mathematical tables.

VIIIo CONCLUSIONS From the foregoing chapters of this study the following conclusions may be drawn: 1. The basic partial differential equations of water hammer including effect of friction, which contain nonlinear terms and have not been solved, can be solved directly by the method of characteristics, with the aid of digital computer. 20 Solutions by the above mentioned method agree with solutions made by earlier methods when the friction term is neglected. Inclusion of this term clearly and sensitively depicts the damping effect. 3~ The variety of examples given in this study shows that the method is capable of handling many different gate operations and pipe connections easily and adequately. 4o The experimental studies made for the case of instantaneous valve closure check the validity of the theoretical solutions. The attenuation of pressure due to friction evaluated by this method9 agrees with experiment unless the velocity is too low. Hence, the method is able to provide an accurate evaluation of friction effect in most engineering cases. 5e The experiments also verify that the theoretically obtained curves of pressure versus time for different pipe systems and branch connections can correctly and sensitively describe the actual pressure variations. -68

6. From the analytical work and experimental verification, it is believed that the method presented in this study affords an adequate and rapid means of analyzing water hammer phenomena allowing friction in any pipe system and branch connection for any desired boundary conditions. 7o It has been observed that this method has remarkable advantages of exactness, directness, quickness, flexibility, and wide applicability over earlier methods, whether the friction is included or not. For a more thorough verification of the analytical work, further experimentation for slow valve closure and for non-uniform pipes are recommended, During the course of experimental investigations, it was observed that the minimum pressure inside the pipe dropped below the vapor pressure a few times. A better analytical study for this case, considering the energy dissipation due to both ordinary wall friction and resiurge effect, are also suggested. The theoretical solutions and experimental results showed unsatisfactory agreements in case of extremely low velocity. A further study for this case, both in theoretical analysis and experimental investigation, are highly desiredo It may require a three-dimensional analysis or re-examnination of assumptions. It may more likely require further modification or improvement in experimental equipment and technique as the first step of the study, Research on water hammer in a very elastic pipe will also be an interesting subject,

It may be worthwhile to prepare at a later date a library of problems in water harmner, in which complete workable computer programs for the rapid solution of many representative problems based upon the principle and method presented in this study will be collected, sa that a great number of solutions can be readily found by simply providing the type of problem and necessary data~ There will be certain cases in which the effect of hysteresis plays a significant role in the energy dissipation, Since there are little rigorous mathematical expressions and reliable data available for this effect, the'theoretical solution of these problems may becbme very difficult~ Nevertheless, experimental evaluations may be feasible by deducting the computed loss due to wall friction, minor losses, etco, from the total observed losso

APPENDIX I A MAD Language Program(3) for the solution of water hammer including effect of friction losses for the pipe system Case I and a part of its computed results, [See Figure 15a] Subroutines included * [See also APPENDSX VI] * F = 1 for R. Lo 64 is used in order to prevent F(:f) -, 00 In this case V 4 O. -71

-72CHINTU LAI Q042N 11 005 030 * EXECUTE, DUMP * COMPILE MAD, PUNCH OBJECT WH11 7 INTEGER I,JK,M,NTREC,TP DIMENSION V(40), H(40), V1(20)9 H1(20)9 V2(40), H2(40), VP(40 2), HP(40) PRINT FORMAT TITLE HERE READ FORMAT CARD9 LL2,L3,VOHOtNTCD9D2tBB2,NUS1S2,TPE, 2E29EW PRINT FORMAT GIVEN, L,L2,L3,VO,HO,NTCDD2,B,B2,NUS1,S29E9 2E2,EW Q=1./EW+0.81*D/(B*E) Q2=1./EW+0.81*D2/(B2*E2) A=SQRT.(1/(1.93*Q)) A2=SQRT.(1./(1.93*Q2)) J=2 K1l DELL=L/N HF=FR.(D,VO,NU,S1)*DELL*VO*VO/(D*64.332) HF3=HF*L3/DELL DELTI=L2/A2 DELV=O.01 DELH=O.O1 T:O HP(O)=HO PRINT FORMAT TEST, Q, A, Q2,A29 DELL, DELTI PRINT FORMAT TIME, T X=-DELL THROUGH ANN, FOR I=O, 1, I*GoN X=X+DELL V ( I ) =VO H1 I)=HO-I*HF ANN PRINT FORMAT RESULT,XV1(I)IH H1I) PRINT FORMAT BLANK THROUGH ARBOR9 FOR I=N+J+K,*1,I*GN+2*J Vl(I) =VO H1(I)=H1(N)-(I-N-J-K)*HF3 PRINT FORMAT RESULT,XVl(I),Hi(I) ARBOR X=X+L3 HE=Hi(N+2*J) PRINT FORMAT BLANK XO THROUGH MICH, FOR I=N+K,1,I*GeN+J V1(l)=O. HI(I)=H1(N) PRINT FORMAT RESULTX9V1(I)oH1(I) MICH X=L2 HF12=HF/2. HF32=HF3/2. THROUGH USA, FOR I=0, 1, I.G*2*N V2(I)VO0 USA. H2(I1)= HO-I*HF12 THROUGH CANADA, FOR I=2*(N+J+K),l,I.G2*(N+2*J) V2(I)=VO CANADA H2 ( I )H2(2*N)-(I-2*(N+J+K) )*HF32 THROUGH ENGLAD, FOR I=2*(N+K),1*I.G.2*(N+J) V2 ( I ) =0 4..........ENGLAD H2 ( I ) =H2(2*N) TH=DELTI/DELL

-73TH2=DELTI/L2 TH3=DELTI/L3 PHILA TREC=O NEW T=T+DELTI WHENEVER ToG.3*TC*AND*T.G*20**(L+L3)/A, TRANSFER TO HERE TREC=TREC+1 WHENEVER TREC*NE.TP, TRANSFER TO YORK PRINT FORMAT TIME, T YORK M=1 DELT =DELTI TI=T THROUGH WASH, FOR I-0,1,I*.GN+2*J V(I):V1(I) WASH H(I)=H1(I) BOSTON THROUGH LONDON, FOR I=1, 1, I*E.N VR=V( I )*( -TH*(V(I)+A))+V(I-)*TH*(V( I )+A) VS=V(I)*(1.+TH*(V(I)-A))-V(I+)*TH*(V(I)-A) HR=H(I)*(1.-TH*(V(I)+A))+H(I-1)*TH*(V(I)+A) HS=H(I)*( 1.+TH*(V(I)-A))-H(I+1)*TH*(V(I)-A) VP(I)=0O5*(VR+VS)+(16.083/A)*(HR-HS)-0.5*FR.(D,V(I)NUS1)*A 2BS.V(I)*V(I)*DELT/D LONDON HP(I)=(A/64.332)*(VR-VS)+0O5*(HR+HS) VS=V(O)*(1+TH*(V(O)-A))-V(1)*TH*(V(0)-A) HS=HO*(1.+TH*(V(O)-A))-H(1)*TH*(V(O)-A) VP(O)=VS+(32.166/A)*(HO-HS)-0.5*FR.(D,V(O),NUS1)*V(0)*.ABS.V 2(O)*DELT/D WHENEVER M*.E.1 TRANSFER TO SPAIN THROUGH BELG, FOR I=N+K+1, 1, I.E.N+J VR=V(I)*(l.-TH2*(V(I)+A2))+V(T-1)*TH2*(V(I)+A2) VS=V(I)*(l.+TH2*(V(I)-A2))-V(I+1)*TH2*(V(I)-A2) HR=H(I)*(l*-TH2*(V(I)+A2))+H(I-1)*TH2*(V(I)+A2) HS=H(I)*(1+TH2*(V(I)-A2))-H(I+1)*TH2*(V(I)-A2) VP(I)=OS5*(VR+VS)+(16.083/A2)*(HR-HS)-0.5*FR.(D2,V(I),NUS2)* 2V(I)*.ABS.V(I)*DELT/D2........ BELG HP(I)=(A2/64.332)*(VR-VS)+O.5*(HR+HS) THROUGH CZECHO, FOR I=N+J+K+1,lI.E.N+2*J VR=V(I)*(1.-TH3*(V(I)+A ))+V(I-1)*TH3*(V(I)+A) VS=V(I)*(1.+TH3*(V(I)-A ))-V(I+1)*TH3*(V(I)-A) HR=H(I)*(1.-TH3*(V(I)+A ))+H(I-1)*TH3*(V(I)+A) HS=H(I)*(1.+TH3*(V(I)-A ))-H(I+1)*TH3*(V(I)-A) VP(I)=O.5*(VR+VS)+(16.083/A)*(HR-HS)-0o5*FR.(DV(I)iNUgSl)*A 2BS.V(I)*V(I)*DELT/D CZECHO HP(I)=(A /64.332)*(VR-VS)+0.5*(HR+HS) SPAIN VS=V(N+K+1)*TH2*A2 HS=H(N+K)*(1-TH2*A2)+H(N+K+1)*TH2*A2 HP(N+K)=HS-(A2/32.166)*VS VP(N+K)=O WHENEVER TI.G.TC, TRANSFER TO ITALY EXECUTE BC.(TCTI,A,V0,HE, V(N+2*J), H(N+2*J),DELVDELH) VP(N+2*J)= V(N+2*J)-DELV HP(N+2*J)= H(N+2*J)+DELH TRANSFER TO OHIO 6.. —-.ITALY VP(N+2*J)=O OHIO VR=V(N+2*J)*(1l-TH3*(V(N+2*J)+A))+V(N+2*J-1)*TH3*(V(N+2*J)+A) HR=H(N+2*J)*(1TH3*(V(N+2*J)+A))+H(N+2*J1-)*TH3*+(VN+2*J)+A) VENUS HP(N+2*J)=HR-(A/32.166)*(VP(N+2*J)-VR)-A*FR.(DV(N+2*J),NU.S1........ --. 2)*V(N+2*J)*.ABS.V(N+2*J)*DELT/(64.332*D) WHENEVER TI.G.TC. TRANSFER TO STAR

EXECUTE BC*(TC.TIA,VOHEVP(N+2*J),HP(N+2*J),DELVDELH) VP(N+2*J)=VP(N+2*J)-DELV HP(N+2*J)=HP(N+2*J)+DELH WHENEVER (*ABS.DELV.GEOQ 1) *AND. (oABS.DELH.GE#OOO1). TRA 2NSFER TO VENUS STAR VP(N+J+K)=V(N+J+K). VP(N+J+K+1)-V(N+J+K+l) VR1=V(N)*(1*-TH*(V(N)+A))+V(N-1)*TH*(V(N)+A) HR1=H(N)*(1.-TH*(V(N)+A))+H(N-1)*TH*(V(N)+A) VR2=V(N+J)*(1l-TH2*(V(N+J)+A2))+V(N+J-1)*TH2*(V(N+J)+A2) HR2=H(N+J)*(1l-TH2*(V(N+J)+A2))+H(N+J-1)*TH2*(V(N+J)+A2) VS3=V(N+J+K)*(1l+TH3*(V(N+J+K)-A))-V(N+J+K+1)*TH3*(V(N+J+K)-A 2) HS3=H(N+J+K)*(1.+TH3*(V(N+J+K)-A))-H(N+J+K+1)*TH3*(V(N+J+K)-A 2) GREECE HP(N+J+K)=HS3+(A/32.166)*(VP(N+J+K)-VS3)+A*FR,(DV(N+J+K),NU, 2S1)*V(N+J+K)*.ABS.V(N+J+K)*DELT/(64.332*D) HP(N)=HP(N+J+K) HP(N+J)=HP(N+J+K) VP(N)=VR1-(32e166/A)*(HP(N)-HR1)-HR.5*FR(DtV(N),NU,S1)*V(N)*o 2ABS.V(N)*DELT/D VPCN+J)=VR2-(32.166/A2)*(HP(N+J)-HR2)-0.5*FR.(D2,VCN+J) NU,S2 2)*V(N+J)*.ABS.V(N+J)*DELT/D2 VEL=VP(N)+VP(N+J)*D2*D2/(D*D) DIF=VEL-VP(N+J+K) WHENEVER *ABS.DIF *LE0O0011 TRANSFER TO POLAND VP(N+J+K)=VP(N+J+K)+DIF/3. TRANSFER TO GREECE POLAND WHENEVER M.NE1,. TRANSFER TO BERLIN THROUGH PARIS. FOR I=O,1,I.G.N+2*J V1( I )VP(I) PARIS H1(I)=HP(I) M=2 N=2*N J=2*J K=2*K DELT=DELT/2. TI=T-DELT THROUGH ROME. FOR I=O.1,I*GN+2*J V(I)=V2 (I) ROME H(I):H2(I) 12- TRANSFER TO BOSTON BERLIN WHENEVER M*NE*2, TRANSFER TO UNIV THROUGH TOKYO, FOR I=Ol,*I.G.N+2*J V( I)=VP(I)........ TOKYO -- H(I)=HP(I) M=4 9 —-~ ~TI=T TRANSFER TO BOSTON UNIV THROUGH SUN9 FOR I=O,1,ING*N+2*J V2( )=VP( I) SUN H2(I)=HP(I) N=N/2 J=J/2 K=K/2 WHENEVER TREC.NE.TP, TRANSFER TO NEW X=-DELL THROUGH MOON.t FOR I=0,1,I.G.N X=X+DELL

-75VB=2*V2(2*I )-Vl( ) HB=2*H2(2*I )-H1( I ) MOON PRINT FORMAT RESULT, X,VB,HB PRINT FORMAT BLANK THROUGH TURKEY. FOR I=N+J+Ktl,I.G.N+2*J VB:2*V2(2* )-V1l( ) HB=2*H2(2*I )-H1 ( I ) PRINT FORMAT RESULT, XVB,HB TURKEY X:X+L3 X=O PRINT FORMAT BLANK THROUGH CONGO, FOR I=N+K,1,I.G.N+J VB2*V2(2*I )-V1 ( I ) HB=2*H2(2*I )-Hl I ) PRINT FORMAT RESULT, XVB.HB CONGO X —L2 TRANSFER TO PHILA VECTOR VALUES TITLE=$64H1WATERHAMMER PRODUCED BY GATE CLOSURE 2 INCLUDING FRICTION EFFECT /1HO510,50HPIPE LINE WITH A SHOR 3T DEAD-END BRANCH CONNECTION *$ VECTOR VALUES CARD = $5F10.3,I4/5F8.5,3F10.8,I3/3El2.5*S VECTOR VALUES GIVEN 2 $5H4L = FlO.3,S8#5HL2 = F10.3,S895HL3' 2 F10.3,.S85HV0 = F10.3,S8,5HHO = F10.3/5HON = I4,585HTC = FS 3.5,5S84HD = F8*.5S8,5HD2 = F8.5/5HOB = F8.5,S8,5HB2 = F8.5,S8 4,5HNU = FlO.8,S8,5HS1 = F1O.8,58,5HS2 = F10.8/5HOE = E12.5,S4 5,5HE2 = E12.5ES4+5HEW = E12,5*$ VECTOR VALUES TEST = $lHOt4HQ = E12.5,S5,4HA = E125,S55,5HQ2 2= E12.5,S5,5HA2 = E12*5/8HODELL = E12,5S55t8HDELTI = E12*5*$ VECTOR VALUES TIME = $8H4TIME = F1O.5/lHO,512,1HX,515,8HVELOC 2ITY S12,8HPRESSURE*$ VECTOR VALUES RESULT = $1H,S7,F1l03,S8,El24,S8+,El24*$ VECTOR VALUES BLANK = $2H0 *$ END OF PROGRAM * COMPILE MAD, PUNCH OBJECT WHS8 3 EXTERNAL FUNCTION (D,V,NU,K) ENTRY TO FR. R-D*.ABS.V/NU WHENEVER D *G. 0.027, TRANSFER TO LARGE WHENEVER R L 64. F=1*0 OR WHENEVER R~L~200. F=64 /R OR WHENEVER RL.1050, F=17.77*R.P. (-0.7581) OR WHENEVER R~L.3250. F=6.948*R.P. (-0.6232) OR WHENEVER R*L.3600. F=498.31*Re.P (-1.1516) OR WHENEVER R*L.4200. F=8.0349*R.P. (-0.6476) OR WHENEVER R.L.5000. F=0.51216*R.P, (-0.3176) OR WHENEVER R.L.6000. F=0.13672*R.P. (-0.1625) OTHERWISE F=0.076236*ReP. (-0.09538) END OF CONDITIONAL TRANSFER TO BACK

LARGE WHENEVER R*.L64. F=1.O OR WHENEVER R.L.200. F=64./R OR WHENEVER R.L.1050. F=14.49*RP. (-0.7197) OR WHENEVER R*L.3650. F=7.525*RP. (-0.6258) OR WHENEVER R*L.4000. F=624.98*R.P.(-101643) OR WHENEVER R.L.4700. F=6.3212*R.P.(-0.6104) OR WHENEVER R*.L5700. F=0.6013*R.P.(-0.3322) OR WHENEVER R.L.7000, F=0.2424*R.P.(-0.2271) OTHERWISE F=O.07819*R.P. (-009933) END OF CONDITIONAL BACK FUNCTION RETURN F END OF FUNCTION * COMPILE MAD, PUNCH OBJECT WHS7 EXTERNAL FUNCTION (TC,TIA*VOHOHO,*VH.DELVDE LDELH) ENTRY TO BC. TAU=O.5*(1.5-TI/TC)*(1.5-TI/TC)-. 125 C=A*VO/(32.2*HO) W=V/VO Y=TAU*TAU X=2 *W+C*Y DLV=O.5*(X-SQRT.(X*X-4.*(W*W-Y*H/HO))) DLH=C*DLV DELV=VO*DLV DELH=HO*DLH FUNCTION RETURN END OF FUNCTION 12

!.-' ii 3 i i} 0 I,' i: J f'{ i_- C.' I_-!' J-J L-1 i-t ~:.i E::'.,_.Z-C!~~~~~~~~~~~~~~C ij. 2 i3.9 E' % L?' 0 7 9 3 1:3 iS!'' 0' S L~ U L L L'!~~~~~~~~~~~~~~~ E:,:i i.........._.~...._... ~!....;~ 5 0'~'1.i~ CI ~ 0 =.~L!, J_.. T 717~~~~~i'78'3 7'!i 3E'.b. -'C _i 3."" E' ~ [' i_-! I-I i i-i'? ~' n1.._ 3 U'-L.. _' C(~................~-..','i!~~17.,Er-pj&O _. __."2.' n..C'"Yr:~._;;j'!Ti 3 (!,;~ 2''Pi' Ci I. (l 3 F:,L!. E' ~ 0 0 l!3',-I O, L 3 i. Cl 2'' i isI S I l 4( 3 {L i- i 3 4,'1 i- iI- P',,~' ~-I [-1,r'-~~~~~~~~~C-I- ---- - - Fj.i'%'.t I i i' Ii 3 N I'.; C- 3! 0'! n 1- r' L I';'o~ ~ ~ ~ ~ ~ ~ ~ ~ ~~~~~~~~~~~~~~~- -;S~F-!:.... -.- -a:',,' -c.;',-:-,,':, ---- 31~~~~~~~~~~~~~~~~~~~~ L:-,'.3~,...;."!39''':!.':'' 0 T:-: 3 2.'i'0~~: e'q2' F!~~~~~~~; J. 0 9U,::-' i*} ~' [:-:.?:3.? ~ - 71: i[! [ ~.]! 3!]i~ ~ ~ ~~~~~~~C,:I C''D 1z 0E -/ I...................~~~~~~~~~~~~~~~~~~~~~~~~~:-?..-~ -, -~ -:' --- -,-, -: -,- -,.-~-,... F! F~~~~~~~~~~-!? -. -iC i.....- ~h LL~ ~ ~ ~ ~ ~ ~ ~~~~~~~~7. - - -:'. -. -~:. - - - - -i - - - 0 O!~~~~~~~~~~~~~~~~~~~~ ~~~ ~ ~~~~~~ ~~~~-i I- - - -5 -..... - - - - -- - - i-i - -!.-: -' -! --- -- ----- LI =,3 C8,"-'.QFi0=: Q ~ ~~~~ i. 7 - N- I - F'i H 0:=', n ~:;:- --,,',, 5i.:i 9 ";."T.:= -' 9 n 0n ~,."Z. -,IT0 -----------------------------------------------------------!..l~~~ ~~~~~ 70 i..i 9 -..~ ".[U:.:. 8 0l ~] J.;a-. H'.-.,.. H LI I T'' 7.;/1 3,-1,_

APPENDIX II A MAD Language Program(3) for the solution of water hammer including effect of friction losses for the pipe system Case II and a part of its computed results. [See Figure 15b] Subroutines omitted. [See also APPENDIX V.] -78

-79CHINTU LAI Q042N 12 005 030 * EXECUTE. DUMP * COMPILE MAD, PUNCH OBJECT WH12 5 INTEGER G.I, J,K.M,NR,TREC,TP DIMENSION V(60), H(60), V1(30), H1(30), V2(60), H2(60)t VP(60 2), HP(60). L(4), D(4), B(4), S(4)9 E(4), Q(4), A(4), DELL(4). 3VI(2), HF(4). TH(4) PRINT FORMAT TITLE START READ FORMAT CARD, VO, HO0 TC. NU, N. L(1)...L(4)t D(1)...D( 24), B(1)...8 (4), S(1)o*,S(4), E(1)...E(4), TP, EW PRINT FORMAT GIVEN, VO, HO0 TC, NY. N. L(1)***L(4). D(1)..*D( 24), B(1)o**B(4), S(1)*..S(4). E(1)...E(4), EW Z=1./EW THROUGH WORLD, FOR G=l, 1. GG.*4 Q(G)=Z+0.81*D(G)/(B(G)*E(G)) WORLD A(G)=SQRT(1,/(1.93*Q(G)) ) J:2 K=1 THROUGH TOUR, FOR G=1: 1, GG. 2 TOUR DELL(G)=L(G)/N DELL(3)L (3) DELL(4)L (4) VI ( 1)=VO C1D( 1)*D( 1 )/(D(2)*D(2) ) C2:D(3)*D(3)/(D(2)*D(2) ) VI(2)=VO*C1 THROUGH FROM. FOR G=l, 1. G.G#2 FROM HF(G)=FR. (D(G),VI(G),NU,S(G))*DELL (G)*VI (G)*VI (G) / (D(G)*6433 22) HF(4)=HF(2)*L(4)/DELL(2) DELTI=L(3 )/A(3) DELV=O O01 DELH=O.01 Hl(O)=HO..... H2(0)=HO HP (O)=HO PRINT FORMAT TEST. Q(1)e..Q(4). A(1)...A(4). DELL(1). DELL(2) 2. DELTI PRINT FORMAT TIME. T R=O 12. -- - -... X=O THROUGH ARBOR. FOR G=1, 1, G.G.2 X=X-DELL (G THROUGH ANN, FOR I=R, 1. I*G.N+R. XX+DELL ( G ) V1( I )VI (G) H1( I )-H1(R)-( I-R)*HF(G) ANN PRINT FORMAT RESULT, X. Vl(I)s H1(I) 8............... PRINT FORMAT BLANK R=R+N+J ARBOR Hl(R)=H1(R-2) THROUGH MICH, FOR I=R+K, 1. I.GeR+J Vl ( I ) =VI 2) Hl(I)=Hl(R)-(I-R-K)*HF(4) PRINT FORMAT RESULT. Xt V ( I) H1 I) MICH X=X+L(4) 4. II=HEH1 (RI+J ). PRINT FORMAT BLANK 9 ___ _._ ~.

-80X=O THROUGH USA, FOR I=R-K. 1, I.G.R V1(I)=0, H( I )=H1(R) PRINT FORMAT RESULT, Xt V1(I), H1(I) USA X=L(3) THROUGH OHIO, FOR G=1, 1, G.G.4 OHIO HF(G)=HF(G)/2. R=0 THROUGH MD* FOR G=i, 1, G.G.2 THROUGH PA. FOR I=R, 1, I.G.R+2*N V2( I )=VI (G) PA H2( I )H2(R)-( I-R)*HF(G) R=R+2*(N+J) MD H2(R)=H2(R-4) THROUGH WASH, FOR I=R+2*K, 1, I*GR+2*J V2( I )=VI(2) WASH H2 ( I )=H2(R)-(I-R-2*K)*HF(4) THROUGH NEW, FOR I=R-2*K,1,IeG.R V2( I ) =0. NEW H2 ( I ) =H2(R) THROUGH YORK, FOR G=1, 1, G.G.4 YORK TH(G)=DELTI/DELL(G) BALTMR TREC=O BOSTON T:T+DELTI WHENEVER T.G.3*TC *AND. T.G.16.*(L(1)/A(1)+(L(2)+L(4))/A(2)), 2TRANSFER TO START TREC=TREC+I WHENEVER TREC.NE.TP, TRANSFER TO PHILA PRINT FORMAT TIME, T PHILA M=l DELT=DELT I TI=T THROUGH CANADA, FOR 1=0, 1, I.G2*(NN+K+J) V( I )=V ( I ) CANADA H(I)=H1( I ) SCOT R=0 THROUGH ENGLAD, FOR G=1, 1, G.G.2 THROUGH LONDON, FOR I=R+1, 1, I*E.N+R VR=V( I )*( 1.-TH(G)*(V( I )+A(G)))+V( I-1 )*TH(G)*(V(I )+A(G)) VS=V( I )*( I+TH(G)*(V( I )-A(G)) )-V( I+)*TH(G)*(V(I )-A(G) ) HR=H( I )*( 1-TH(G)*(V( I )+A(G)))+H( I-1 )*TH(G)*(V( I )+A(G) ) HS=H(I )*( 1.+TH(G)*(V( I )-A(G) ) )-H( I+1)*TH(G)*(V( I )-A(G) VP ( I )=0.5*(VR+VS)+(16.083/A(G))*(HR-HS)-0.5*FR.(D(G)V( I ),NU 2S(G))*V(I)*ABSV( I )*DELT/D(G) LONDON HP( I )=(A(G)/64332 )* VR-VS)+. 5*(HR+HS) ENGLAD R=R+N+J VS=V(O)*(.+TH( )*(V(O)-A( l)) )-V(1)*tTH( 1)*(V(O )-A( ) ) HSHO*(1.+TH( 1)*(V(O)-A(1) ) )- H(1)*TH()*(V(O)-A(1) ) VP(O)=VS+(32.166/A(1) )*(HO-HS)-0.5*FR. (D 1) V(O) NUS( 1)*V(O - 2)-.*ABS.V(0 )*DELT/D( 1 ) WHENEVER M.E#1i TRANSFER TO SWEDEN THROUGH NORWAY, FOR G=3, 1, G.G.4 THROUGH OSLO. FOR I=R-K+1. 1, I.E.R VRV( I )*( 1.-TH(G)*(V( I +A (G) )+V( I-1 )*TH(G)*(V( I )+A(G)) VS=V( I )*(.+TH(G)*(V( I )A(G)))-G) )V( I+1)*TH(G)*(VI )-A(G)) HR=H( I )*(1.-TH(G)*(V( I )+A(GI) )+H( I-1 )*TH(G)*(V( I )+A(G ) HS=Ht I *( 1.+TH(G)*(V( I ) -A(G) W )-H( I+1 )*TH(G)*(V I )-A(G)) 3- ~ -- -...

-81VP(I)=.05*(VR+VS)+(16.0 83/A(G) )*(HR-H S)- 5*FR(D(G)V(I)NU 2S(G) )*V( I )*.ABSV( I )*DELT/D(G) OSLO HP(I)=(A(G)/64.332)*( VR-VS)+O.5*(HR+HS) NORWAY R=R+J SWEDEN R=2*N+3*K VS=V(R+1)*TH(3)*A(3) HS=H(R)*(1l-TH(3)*A(3))+H(R+1)*TH(3)*A(3) HP(R)=HS-(A(3)/32.166)*VS VP(R)=O R=2*(N+K+J) WHENEVER TI.G.TC, TRANSFER TO FINLAD EXECUTE BC.(TCTI,A(4)VI(2),HEoV(R), H(R)o DELVDELH) VP(R)=V(R)-DELV HP(R):H(R)+DELH TRANSFER TO MOSKVA FINLAD VP(R)=0 MOSKVA VR=V(R)*( 1.-TH(4)*(V(R)+A(4)))+V(R-1)*TH(4)*(V(R)+A(4)) HR:H(R)*(I.-TH(4)*(V(R)+A(4) ))+H(R-1)*TH(4)*(V(R)+A(4)) RUSSIA HP(R)=HR-(A(4)/32.166)*(VP(R)-VR)-A(4)*FR,(D(4)tV(R)gNU,$(4)) 2*V(R)*oABS.V(R)*DELT/(64.332*D(4)) WHENEVER TI.G.TC9 TRANSFER TO POLAND EXECUTE BCe(TC,TI,A(4),VI(2).HE,VP(R)HP(R),DELVDELH) VP(R) =VP(R)-DELV HP(R)=HP(R)+DELH WHENEVER (~ABSeDELV.GE.O.OO1)eAND.(.ABS.DELHeGE.O.001)9 TRANS 2FER TO RUSSIA POLAND R=2*(N+K) VP(R+J+K) =V(R+J+K)+VP(R+J+K+1 )-V(R+J+K+1) VR2=V(R)*( 1.-TH(2)*(V(R)+A(2)F) )+V(R-1)*TH(2)*(V(R)+A(2) ) HR2=H(R)*(1.-TH(2)*(V(R)+A(2)))+H(R-1)*TH(2)*(V(R)+A(2)) VR3=V(R+J)*C(1.-TH(3)*(V(R+J)+A(3)) )+V(R+J-1)*TH(3)*(V(R+J)+A( 23)) HR3=H(R+J)*(1q'TH(3)*(V(R+J)+AC3)))+H(R+J-1)*TH(3)*CV(R+J)+AC 23)) - - ~ "~ ~~VS4=V(R+J+K)*(1.+TH(4)*(V(R+J+K)-A(4)))-V(R+J+K+1)*TH(4)*(V(R 2+J+K)-A(4)) HS4=H(R+J+K)*(1.+TH(4)*(V(R+J+K)-A(4)))-H(R+J+K+1)*TH(4)*(V(R 2+J+K)-A (4)) ERLIN...H FCR+J+K ).=HS4+.(A (4) / 32. 166)* (VP (R+J+K)-V54)+A( 4)*FR. (D (4), V (R 2+J+K),NUS(4))*V(R+J+K)*.ABS.V(R+J+K)*DELT/(64.332*D(4)) 2, - -..................................P (R).=HPI(R+Ji+K) HP(R+J)=HP(R+J+K),....VP(R)-VR2( 32. 166/A(2 ) )*(HP(R)-HR2 )-O.5*FR (D(2),V(R),NUS(2) 2)*V(R)*oABS*V(R)*DELT/D(2) 10o Rl.3LP..i V R3- ( 32.166/A ( 3) )* ( HP ( R+J) -HR3 ) -O ~ 5*FR. ( D ( 3 ) V ( R+J),N 2US(3))*V(R+J)*.ABSeV(R+J)*DELT/D(3) ~~~9 V~~~~..%....EL=VP(R )+VP ( R+J) *C2 DIF:VEL-VP(R+J+K)...WHENEVER *ABS.DIFL.O.OO1 TRANSFER TO DEUTCH VP (R+J+K)=VP(R+J+K)+DIF/3. 7..........TRANSFE.R..TD. BERLIN.. DEUTCH VP(N+J)=V(N+J)+VP(N+J+1)-V(N+J+1)._____6 A/............R15.VN.(N)*1r —TH(1)*(V(N)+A(1) ))+V(N-1)*TH( 1)*(V(N)+A( 1)) HR1:H(N)*(1.-TH(1)*(V(N)+A(1)))+H(N-1)*TH(1)*{V(N)+A(1)) - -.;_ -_._._........2=.(N+J ).2,J( 1+TH(2)*(V(N+J)-A (2)))-V(N+J+1)*TH( 2)*(V(N+J)-A( 22)) 4,._HJS2-H ( N+._Ll,_1 IH( 2 ) *( V( N+J )-A ( 2 ) ) )-H( N+J+1 ) * TH( 2 ) * C VC N+J )-A( 22)) a- ~......... —.. —~ - ~. — --.- -. — -~ -- ~ ~ — - ~- ~ - ~- —. —- -.

-82NETH HP(N+J)=HS2+(A(2)/32.166)*(VP(N+J)-VS2)+A(2)*FR.(D(2) V(N+J), 2NU,S(2))*V(N+J)*.ABS.V(N+J)*DELT/(64332*D( 2)) HP(N):HP(N+J) VP(N):VRl-(32166/A(1) )*(HP(N)-HRl)-O.5*FR(D(1),V(N) NUS(1) 2)*V(N)*.ABS.V(N)*DELT/D(l) VEL=VP(N) *C DIF=VEL-VP(N+J) WHENEVER ABS.DIF.L.0*001, TRANSFER TO BELG VP(N+J)=VP(N+J)+DIF/2. TRANSFER TO NETH BELG R=R+2*J WHENEVER M.NEIl, TRANSFER TO SWISS THROUGH PARIS, FOR I=O0 1, I.G.R V1(I)=VP(I) PARIS H1(I)=HP(I) M=2 N=2*N J=2*J K=2*K DELT=DELT/2. TI=T-DELT THROUGH FRANCE* FOR I-0, 1. I.G.2*R VtI)=V2(I) FRANCE H(I)=H2(I) TRANSFER TO SCOT SWISS WHENEVER M.NE*.2 TRANSFER TO ITALY THROUGH ROME. FOR I=O,0II.G.R V( I )=VP(I) ROME H(I)=HP(I) M=4 TI=T TRANSFER TO SCOT _ ITALY THROUGH SPAIN, FOR 1=0, 1, I*GeR V2(I)=VP(I) SPAIN H2(I)=HP(I) N=N/2 J=J/2 K=K/2 WHENEVER TREC.NE.TP, TRANSFER TO BOSTON R=O 12................... X=O THROUGH AUST. FOR Ga1, 1, G*G*2,,........X=X-DELL(G) THROUGH WIEN, FOR I=R, 19 I.G.N+R 10 X-=X+DELL (G) VB=2*V2(2*I)-V( I) HB=2*H2(2*I)-H1(I) WIEN PRINT FORMAT RESULT. X,VB,HB PRINT FORMAT BLANK AUST R:R+N+J 7 I T. —-) TIHROUGH GREECE. FOR I=R+K. 1 IG, R+J VB=2*V2(2*I )-V( I) HB:2*H2(2*.I)-H ( I ) PRINT FORMAT RESULT, X.VB,HB GREECE...X=X+L (4.)........... X=O PRINT FORMAT BLANK THROUGH TURKEY, FOR I=R-K, 1. I.G.R

-83VB=2*V2(2*I)-V1(I) HB:2*H2(2*I)-H1( I) PRINT FORMAT RESULT. X.VB.HB TURKEY X=X+L(3) TRANSFER TO BALTMR VECTOR VALUES TITLE = $75H1WATER HAMMER SOLUTIONS FOR COMPOUN 2D SYSTEMS INCLUDING FRICTION EFFECT (A) *$ VECTOR VALUES CARD = $2F10.32Fl0.8,14/4F103,4F8.6/4F8.6,4F1 20#8/4E12.4I3/E12,24*$ VECTOR VALUES GIVEN = $6H4VO = F10.3,S8,5HHO = F10O.3S8#5HTC 2= F1O.8tS895HNU = F10.8,S894HN = 14/14HOL(1,2,3,4) = 4(S55F10 3.3)/14HOD(1,29394) = 4(S7,F8.6)/14HOrB(129394) = 4(S7#F8.6)/1 44HOS(1,2,3,4) = 4(55,F10.8)/14HOE (1,2,3q4) = 4(S3,E124)/6HOE E 5W = E12.4*$ VECTOR VALUES TEST = $14HOQ(1,2,3,4) = 4(S59E12.5)/14HOA(1t,2 23,4) = 4(S5,E12*5)/11HODELL(1) = E12.5,S5,1OHDELL(2) = E12.5, 3S5,8HDELTI = E12.5*$ VECTOR VALUES TIME = $8H4TIME = F10.5/1HOS1291HXS15o8HVELOC 2ITY S12 8HPRESSURE*$ VECTOR VALUES RESULT = $1H,S79F10.3,58,El2.4,S8,E12.4*$ VECTOR VALUES BLANK = $2HO *$ END OF PROGRAM --- -.. ~ -- - - - - --- 3-. -—. — -.

-84Wt P, T E F. H A M t,1E R' - ~ r]: qHF -tUt]E_,.- ST O t'i E; ]: t'-!:: L Ij D1 I "IJ. FFR I F: T.I O t'-4,EF F El-: T,_!:: J,~llR-TER -!7:, 4 J:, -. 0. -.:f' _'z~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~O 252r':i u102' —..~ 2!.':i~~~~_~~_ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~. - 0 r: ~'T.: (_- ]:2.3: 4.:) - i_'-!j?'i 0~ I"- i-'i-' i-1i~i Ci-i I)~(' ~ "1i0t'ti i".0i` " ]" i")': i-`i,C L-i?: 0~ I'-i -li~i 0_ i'E.,::1 ~2,3, 2...3:'4: E 1 O E i4 i E!..I -. 4.71 LE 0', I~.!,:21:, 2..,.3., 4 ":, - Fi, 2 5 7.':': o E- C 7 FI, 4 4 7 E - Ji'7'!, 2,4.4' E- E F!, 2 4 4 4;' E.- 0 T,q~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~',:] -', 2.., -_;. 4:, - z 04. 4 4 _: 6 EC" 0 4-.,r3 7E.4 --------------------------------------------------------------------------------- E! E L I O:.' 1 Z:' -!':~. 2':)9 0 5 E 0 2 D E L L (. 21:, -" O. 2 9 4.2; 0 E 0 2 D LI T?.. - O. I 0 (:., 5 4 E —0 2 T T HIE - 0 O OQ FJ0-i FJ TIME =~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ E' X~ I,., E L i3t I.: T T ","; F:' F;.' ES %'; Ijr R E i..- 9. o.5 o.111'2 E 0_1 O.4.4. ~~i E (:.!. 0.1112E;_if~~~~~~~~~~~~~~~~~UEI ---------- -- f- T ----------— i-di ~ --- ----------------- ~p~- 7 -------------'' ----— r —------- ----------- R~ _, T!i. 1 1 C i ~ t"i 1 O. 4. 4. 09E Ei.3~ 71 1:, 05'O.111 2E. 4, E. 14-9.55,i11' 0.'4 E C.3'2'.1 9:':Z47L~ O. 1112' E 01 O.437 5E J _-. — ~ —---- -------------— ~'1 —- --- ---------- -—.j —-~~ — -— ~3iiE -- I 7 I I 11-3 E 1 O.,.. 2:.SE 1!.Z,.E _...357 J.I1II. 21::E I014 7.471E u:.3-.453, 7..' 0 a. 21E: 5 E.' 4 2' - - - - -`! - - - - - - - - - O. 425- - - -'' d. 0 0~ 18:5E L-!i I. 4 229E 0.3 -- ------------— l~i~-~;. —-75i L —------- 5 E 1i 4::..!.E0 7.0.2:.', E'1 i.4175E r!3 5.3 4.'4'n 077. 1:-' 5 E 1i. 4!7E 0. 3 $:,.39: 1:- 016 =7! 41f- E!i......................................... 56 -i3.'.5. 2 i - 5E 1 -,- _3 5? l., 2,_..5"-i ~ 2'~ 8 5 E (.? 1 Ci. 4. ii8 E E I.3 U1.., E F'. 0 -.F": 3 E 4.. 9!35 O. 0~~~~0 i0 IF -- - iF. C-l.:.. E n.?. - - - - - - - - - - - - - - - -- - - - - - - - - - 0 0. 000! 1 rZE I' - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - - - - - -

C~", ~i! C) -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~... -:.. -'-.!..-.....).~~~~~~~~~~~ ~ ~ ~ ~~~~~............... ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~: " " " ~......I: f~ i ~-.....:.....:.....:'.... L.B~~~~~~~~~~~~~~~~~~~~~~~~D E:'II (J..!EZ;. ~~~~~~~~~~~~~~-:-... (....'.I%:"..).....;;)'.: -: C.) C:- (~~~~~~~~~~~~~~~~~~~~~~~~~~:D................... q;q;..'.)'.!)."" O.;T... 4:...O (. i. —I:.!-::.'-''.,.!.;t'' j ";Z)C:";..J'._: I; ~;.. _)-. ~~~~~~~~~~~~~~~~~~~~~~~~~~l.j C~ ~ ~ ~ ~~~~~~~~~~D...)........................... ~. —— I C;I~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

APPENDIX III A MAD Language Program(3) for the solution of water hammer including effect of friction losses for the pipe system Case III and a part of its computed results. [See Figure 15c] Subroutines omitted. [See also APPENDJX V.] -86

-87CHINTU L\ I 7$42 N 1 25 Cr. FXECUTE, l,>! CO P I L fIl'-\" DUJCH ^-'-JFC!T CH13 2 1TEGE R G I, J,K, M, N, f' R TREC, T DI?';E.!I \ V('-C) V(6(i) Vl(3) i ( 3 ) V2( 6.) H2 (6) V VP(6 2) H P(G6.))-, (3) C(3),)' (3) 5H(.3P)33) (3)) S3(). DELL(3) 2VI( 3 ) F (3) Ti-( 3) 3,DTNT'-R.VTT TITLE ST A.RT R E, F. V' TC" il! A,, L (1) ~, L 3 ) 2) L (1)..)( 23), (1). [(3.(!).. (3 ) (1. ).( 3) T P F-',.f P I 1 T.F T.T'.' /,...., L(1)oo.L(3). ( 1)....l( 2.... i _., ), Z=1 / /E' TH iROUGH'LLr:,, 9 1 = 1'C (G)=Z +') *18! r) ( tG ) / (.' (X)'' — ( G ) ) J= 2 \O L D A (C) R T ( 1 / (1' a "> ( K=1 THRClGH TCI!R F'R =.. 1P, G.G. 2 TOU.JR D E L ( L ) t (J ) /!' L) f () = (3) V I ( V )VF CI=D( 1) ( ) / ( D(2 ) (2 ) ) C2= ( 3) () / ( (2 )*D ( ) ) VI (2) =V C1 THrOJGH - FRI.'in cO 0-=, 1~ G+G F..2 FR()O F))F ( L....., )'() )'-" —LL (2)':(',.'? (r:) -) V(T ( -)/(2(i)-,.64.33. 22) DELT I=L( ) /A( 3) DELV='. * 0 1 D LH= 0 o 1 T=O H ] (()= H 9 H2 ( 0) =H HP( 0 ) =HO PR I NT FOR',AT TE.ST 1 ) o. ( 3),9 ( ).. " ( 3 ) 3 -'DFL ( ), DELL(2) 2, DELTI PRINT F Oi' AT TI -E T 0=\ X=O THROUGH ARROR, FOR =-1, 1, GI.22 X=X-DELL (G) Ti ROUGH AN f'l, FOR I R, 1, I.GI.'+R X=X+DELL (() V1(I)=VI (C) H( I )I) == HI(R)-(I-R)-) HF (G) ANN PR I IN!T FOR;MAT RESU LT ), X V( )I,!1 ( I PRIN T FOR A T BL A.!NK R = R + N + J ARBOR H1(R)=H ( R-2 ) HE=H1 ( R) X=O THROUGH M I CH, FOR I =R+k, 1, I.. V1(I )=C* H1( I ) =H1 (N) PRINT FOR'"AT RESULT, N, V( ) I I ( I ) ti I CH X=L ( 3 ) THROUtGH USA FOPR I a = o (1 0 USA HF(G)=HF (G)/2

-88R= = THROUG-H,rD, FC)R G=1! i, GG.r2 T HP'R 0 U(-IGH DA~ D r F _.uR GS 2'.THROJUGH A, FOR I=R, 1, I.(.] 2* V2( I )=V! (0) PA H2 (I) =H2 ( R )-( -R ) HF (G) R=R+2-,(N+J) At.'O H(H2 () =H2 (R-4) THROUGH 2NEW!, FORP 1I=r+2'K 1, I.GR+2*J V2( I ) = b~rlvNEW H2 ( I )=H2 ( 2" N ) THROUIGH YORK, FOR (=1, 11 GG..3 YORK TH ( G ) =- L T I / L L ( G ) 3ALT 1R TREC=O BOSTON T=T+DELTI IHE.N.VER T. 3. 3*TC ~ANI'* T., 16 (L (1 ) /A( 1 )+L(2) //(2)) T2T TRANSFER TO STA TRE= TREC+1!.,'HENE\VER Ti —TEC.i-. TP, TF-:ANSFER TCO HI LA. PR I Ni!T FORCAT T IE, T PHILA V1=l DELT = E LT I T I=T THROUGH CANAA, F OR I=, 1,.G.2- (N+K+J ) V( I )=V I ) CANADA H(I)=H1(I) SCOT R=U THROUGH EN',!GLAD, FOrR G=1, 1, G.G.2 THROUGH LONi')ON, F I R+, 1, I E +R VR=V (I ) (1.-TH ( G) (V(I)+A ((' ) ) +V(I -I ) -TiH (G)"- (V (I)+( ) ) VS=V( T )' (1.+TH ()'- V ( I )-A () ))-V (I +1)'TH- (G) (\ (I )-A(G) HR = I) ( 1.-T ( G ) ( V ( I +A ) +A ) +H (I -1) TH G ) ( \/ (I) +A ( G ) ) HSH( I )-' ( 1.+TH(G )'- ( I )-A(G) ) )-H( I +1 ) -TH( G ) " (V ( I )-A (G) ) VP(I)=0.5 (VR+VS )+( 16. 83 / () )-( R-H )-.5"FR. ( (G), V( I ),NUJ, 2S ( G ) ) -V( I ) -.ABS.*V( I )'-DFLT/D ( G ) LONDON HP ( I ) ( A (G) /6432 ) V(VR-V S) + 5 R+HS) ENGLAD R=R+N+J VS=V (0) "(1.+TH (1) * (V ( )-A( 1 ) )-V(1 ) -TH( 1 ) ( V ( )-(1) ) HS=H C-(. +TH (1) - (V ( 0) -A (1) )-H ( 1 ) TI(1) * C ()-At1) ) VP ( ) =vS+ ( 32. 166/A ( 1 ) ( H-HS ) -;.5 *FR. ( ( 1), V( (: ), ( 1 ) ) V ( 2 ) "-.ABS.V( ) ELT/D'( 1.'HENEVER,. 1, E TRi!SFER S T SWED:EN THROUGH NORWA.Y 9 FOR I=R+iK1,+ 1 I.E.R+J VR=V(I) (1.-TH(3) * (V( I ) +A (3) ) )+V (-I *TH( 3) ( (I) +A (3) ) VS=V( I ) - C( 1 +TH ( 3 ) (V( I )-A ( 3 ) ) ) -V( C + ) TII ( 3 ) " (V ( I A ( 3) Hr=H I ) ( 1 -TH(3) * (V (I) +A ( ) ) +H (I-1)'Tl ( 3) - (V ( I) +' 3 HS=H (I) -( 1.+TH ( 3) I ) - A ( ) 3) -H (I +"1 ) T ( ) V (I) ) ) V P () = 5 ( V R + VS ) 5- (. 6 83 / ( ) ) *, ( i-i R -H S ) - * ~ 5 -' FR ( D ( 3 ), V ( I ), i". U 2S(3) ) *V (I )*.ABS.V( I DE LT / ( 3) ORWAY HP ( I ) =( ( 3 ) /64.332 ) (VR-VS ) + 5- ( R+HS ) Sv,!EDCIN - R=R+IK VS=V ( R+1 ) TH ( 3 )' (3) HS=H(R)"(1.-TH(3)-A( 3)) H(R+)1 TH ( 3)A(3) HP ( R) =HS- ( A ( 3) /32 166) VS VP(R)=O R=2-' ( N+K) WIHENEVEER TI.G.TC, TRANSFER TO FINILAD - EXECUTE RC.(TC,TI,A(2),VI(2).,HE,V(R), H(R), DELV,DELH) VP ( R) =V ( R ) -DELV

-89WP (D) -H('R ) +OELH t: TC SKVA 1F I.LA D VP) )....,'.... V,3,, T..... 5KV/t "I 2 V ( ) - 1-. AR ( (2) V ( il )'L (:,4 )3 T (2" ( ) + 2 ) )) ( tJ1955 I/ P 1 ( ) -H C ( 22 ),l( \- -A ( )N ( ), ( F) * <S ( 2 2-'$V( -A' iRs. (rV [L' Gi3" (2! vi - e VR TI.E - TC, TR' -" " IN S F' R T' P OLA''XF'CUTE 7-C. TC,TI,/ (2) /,I ('2)!F',V (: ),VP( ~D),-L VP,,!l., i) VP () =V ( R )-DLV -P ( =() HDP ( [- L ) + LELF'.....'.NE,,:: — ( oA'-S.DFLV.SE.".'..:1 ),,AN (."',B.'" 7.1):: - ":' - ",' TlF/:,'FS L2:-; T' RUS I A PC L Ai,!D' — ( 1+J )- V ( I4-J ) +VD ( + J L 1 )- V\ ( i J+1 ) R 1 V= (', )' ( i T,' ( ) V ( 1 ) ) " 1 V ( 1) V' ) I V. —:.- +..) - \,,, 1....:,(-'])'"'T -( ]_ J'". \ T +1) ) VS2=V(,+J ) (1 *+T (2), (V (+,'2)-) / )( ) ) _,,' 2.. " —J 1,) )T I( ) -?V( \2J)-' A( 22)) R=R+2- J VR3V =(X)"-(1.-T () -+T (' + (V (3) ) )+'( ( - )()T)-) (3 ) (Y' ( C ) ++ " (3) -) 7-,.33. T ( i ),3'- ( V ( R ) +( j! ( 3 ) -(T) F -i ) ( 2 ) -' ( \V (!R ) +, ( 3 ) ) E r L I N ( N + J ) = H S 2 + A A ( 2 ) (D)' ( ) ( -- ) -_ iy ) + A ( F ) " F \. ( /i ( 2 ) 3 ) ( )\.2... S ( 2 ) ) -V ( J ) -.AS.V ( I: +J ) -D T/ 64.3:'.2 " -" ( ) ) D- r" ( D1 ) - th ( 1LoJ ) PI-P-!or(?,) =HD ( PO+J ) 2 ),V ( N ) " o A ( S V ( (H )'' ) >. L T /; ( 1 ) V P ( R ) = V R 3 - ( 3 2. 6 / A( 3 ) ) ( H P ( ) - i " *:) )',. i:- "r ( () V ( R )' U 9 S ( 3 ) 2)7(... )' -V:. )) DLT/ " (L 3) /E L= 1 ] VP ( "V ) + C2,-/D) P ( R ) D F=VEL-V ( N+J ).~ I':E[VEcR b I' a L.. 1, i I TRANSFE, TO BEL'VP (J +) =VP ( + J ) + D F/ I B-FG' I r~,A-. l r F E o T,. E R L I r,.3..LG!._ _ —,i T Fr r,,,t1; F TC S'!.- I 1. I TH-ROI. JG'-H. /tiR:;IS. FOR I=,, 1, I~!.G.R'V' ( I)=VP I), T. I S ( ) )' = 2` J=2"-J <-2 K' L T = LT /2 ~ T I =T-DELT THROUGH FRA,'!C, FOR I=, 1, I.. 2V( I )=V2( I ) FRANCE H ( I )H2 ( I ) TRANSFER TO SCOT IS 7I SS!,,HENEVER N. NE.2 TRANSFER TO ITAL v THIROUGH RO'.1E, FOR 1=, 1, I.CR V( I )=VP ) RO?'E H( I)=HP( I) TI =T TRA NSFER TO SCOT ITALY THROUGH SPAIN,~ FOR i=. 1 I.G.R

-90V2(I)-=V(I) = I () ( I) I J/2 f7i-, / 2': \/E'VE!? TDRE:".x!F..TR, Tc'A.",t221F1332 TOq i:.OSTC?.),,,,.j,, X = L0 X=X-I),LL (G( X X=X+.LL (i) \/ 2-'V 2 ( 2-.- I ) -V ( I ) HT,=2- H, 7(2-'-I )-H ( ) P R J T fIiT L':5113 tC, T -D,;"-V2.J.,- 1. I * F 1 4/.'it- 2' x/, (.7-" I )-, ( 2 ) h.,-, = 2' ~ ( 2 -"- I' )- H' ( I ) PRQI;X T F'.Ri':..T R -? LT, XN,\' ^-,F: 001EECE X=_ X+L ( 3! TI4VVFF~I) I., )A,....i1.. T:.' S. I.. V/' — T \ C /'. L cTI TITL - = 7. — r7 -,... ", - OL - I-,> -5 F.) SY T'5." i?, LUDIC!:i-iCTIO F-T (1 ) -I.T..........'. (: -:ii2 r,).} 8a TZ.,./2. F1C:,,,''RF's 6/... 6,r,.&........... -2., &..2'. 4 I = i'. oHT-'" ) 1 2 l r ( *,- ~i X i1 3nLTI 12 V{/F C~nri) JAL!_ "t4 = )1 /' I'T r T Fl r'"/ 1. /':: 1(X 11.1,' HVLC V/+.,2,3 - CT6F,'. F)L / ].T E (1H,117 2 - 3S%,.? / E ~-%-E:,' =,L12.,-': VECT:) VALUES TI'L' 2 -T I". VrE?D) F P R

-91-.......................................................................... Ilij~~~~~~~~~~~~~~~~" F-... q.. -':;: Ci 0 L-.'T ('': (}, ("i (.. E,~;3 i"' (", ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~i I...!::i"i'0 i' i'2..............-~ ~.~. t.....................:r i ~~~~~~~~~~~~~~~.'.i C i ii"'" i'i'} i () i~'] "' - i ~ ~ ~ ]'],{'i?0i}i; ] ~ E: i' i;t ~~~~~~~~~~~~~~............... ~~1] 4-8 F-' i i-i 8, i;- 5.4 8 E! 3 O, 2.4 -; -': - - - - - - - - - -- --- - -- - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - - - - - - - - - - -- - - - --- - - - -- - - - - - - - - - - - - --- - - - - - - - - - - - - - - - - - - - 3~~~~~~~~~~~~D:T 0-;:".i'Q!335-! — ~7~ -=~'-:r: -II —-~~~~~~~~~~ —--- - - -- - -~~~~~~~~~ -:?-? -' -......... ~ -i ~]-3!- - -77-~...........-;-' -,.?:;- -J E"-57 T.............. -- - r-:.:-.-:~ T 0 o.'f~:~;~3:~.'.' —T i............. S: i o i::D I. - - I -i -.-E - - -, - - - - - -. - 9..?- --. i 7.2;E (,...:'i, -. S E i;~~~~~~~~~~~~~~~~~~~~~~~~i 9..~.:.:.:; 5 O.'i iE~?:'- [ —:! 2 -iii: - - - E -)'3 9. I,':i- 5 i-i.! Fi 8.3 i7 1 4:: 4' E 0''? U5ui. 1...7.: ]'i i: _: -''!. ------ 1'-~,'.-:). 5t () O. i - ---- - - - -. - -.' -: -.-------------— F" —_~~~~~~~~~~~~~~~~~~~~- -- --- -- --- --- -~~~ CL~~~~~~~~~~- ----

APPENDIJX IV The Moody Diagram written in MAD Language (3) for the determination of friction factor in pipe flow problems, in the form of a subroutine,* [See APPEND]IX V] * F M 1 for R.L. 64 is used in order to prevent F(=f). 00. In this case V I 0. -92

-93* COMPILE MAD, PUNCH OBJECT EXTERNAL FUNCTION (D#V,NUK) ENTRY TO FR. R=D*.ABS.V/NU WHENEVER R.L.64. F=1,O OR WHENEVER R*L*2000# F=64./R OR WHENEVER KE.*O. TRANSFER TO SMOOTH OTHERWISE TRANSFER TO ROUGH END OF CONDITIONAL tRANSFER TO BACK SMOOTH WHENEVER R.L.3250. TRANSFER TO CRITIC F=0.043 STURB KP=O 86858896*ELOG ( R*SQRT (F) )-08 FO=1./(KP*KP) WHENEVER oABS (F-FO) *LeO.OO1, TRANSFER TO BACK F=FO TRANSFER TO STURB ROUGH RS=D/K WHENEVER RS.G.400..AND.R.G*3250. *OR. RSG. 100..AND.R.G.3350. 2 *OR. RS.G.25..AND.R.G.3550..OR. RS.LE,25*.AND.R*G.3900.t TR 3ANSFER TO RTURB CRITIC W=ELOG.(R/2313.) X3.05396*W*W F=O.03*EXP. (X) TRANSFER TO BACK RTURB W=O086858896*ELOG (RS) +114 F=1./(W*W) WHENEVER R*SQRT.(F)/RS.G.200o, TRANSFER TO BACK TRASN X=1.74-0.86858896*ELOG.(2./RS+18.7/(R*SQRT.(F) ) FO=1./(X*X) WHENEVER.ABS.(F-FO).L.O.0001, TRANSFER TO BACK F=FO TRANSFER TO TRASN BACK FUNCTION RETURN F END OF FUNCTION

APPENDIX V TABLE II SYMBOLS USED IN THE MAD STATEMENTS THROUGH APPENDICES I -' IV AND THEIR EQUIVALENT CONVENTIONAL NOTATIONS OR DESCRIPTIONS Symbols Used In Equivalent Con- Symbols Used In Equivalent ConMAD Statements ventional Notations MAD Statements ventional Notations A,Al or A(1),... a,al,... VP Vp BB1 or B(1),... byb,... VR,VS VR,VS C1 (D1iD2)2 = A1/A2 VR1,VS2,.. (VR)l,(Vs)2Y,C2 (Dj/D2) = A3/A2 X x D,D1 or D(1),... D,D1,... DELH AlH DELLDELL ( 1 ),.... TaLL1.. DELT,DELTI At Symbols Used In DELV AV MAD Statements Descriptions E,E1 or E(1),... E,El,... EW K DIF Difference, A F,FR. f G Integer H,HO,... H,Ho,... HB Extrapolated value HF hf of H HP Hp HE Head across the gate HR,HS HR,HS when V = Vo HR1,HS3,... (HR)l, (HS)3,... H1 H for n steps of 2At I i H2 H for 2n steps of At J j R Integer K k R Reynolds number K c RS Relative smoothness LLl or L(l),... L, L1,... D/e M m TREC Integer (for time N n recording or counting) NU v TP Integer (for the Q 1/K + Dcl/Eb control of time inSS1 or S(1),... E, e l'. terval in printing) TAU T VB Extrapolated value TC Tc of V TH,TH1 or TH(1),.. Q, Q1,.. V1 V for n steps of 2At T,TI t V2 V for 2n steps of At V,VO,... V,Vo,.. VEL V -94

BIBLIOGRAPHY 1. Allievi, L. "Theory of Waterhammer," translated by Eo E. Halmos, printed by Richardo Garoni, Rome, Italy, 1925. 2. Angus, R. W. "Water-Hammer Pressures in Compound and Branched Pipes," Proceedings ASCE, 64, January 1938, p. 133-169. 3. Arden, B,, Galler, B,, and Graham, R. "Michigan Algorithm Decoder," The University of Michigan, Ann Arbor, Michigan, June 1960. 4, Barbarosa, Nicholas L. "Hydraulic Analysis of Surge Tanks by Digital Computer," Journal of the Hydraulics Division, Pro.c of ASCE, 85, No, HY4, April 1959. 5. Binnie, A. M. "The Effect of Friction on Surges in Long Pipe-Lines," Quarterly Journal of Mechanics and A plied Mathematics, Vol. IV, Pt, 3, 1951 )5330-343. 60 Cazard, R, "Fatigue of Materials," translated by Fenner, A. J,, Philosophical Library, Inc., 19535. 7, Durand, W. F, "Hydraulics of Pipe Lines,"'D. Van Nostrand Company, Inco, New York, 1921. 8. Freeman, John R. "Flow of Water in Pipes and Pipe Fittings," Published by ASME, 1941. 9. Joukowsky, N. "Waterhammer," translated by Miss O. Simin, Proc. Amero Water Works Assoc., 24, (1904) 341-424. 10, Lister, M. The Numerical Solutions of Hyperbolic Partial Differential Equations by the Method of Characteristics in "Mathematical Method for Digital Computers" edited by Anthony Ralston and Herbert S. Wilf, John Wiley and Sons, Inc., 1960. 11. Lister, IM, and Roberts, L, "On the Numerical Solutions of Spherical Waves of Finite Amplitude," Technical Report on the Project of Machine Method of Computation, MIT-Proj. DIC 6915, June 1956. 12. Love, A. Eo H. "A Treatise on the Mathematical Theory of Elasticity," Dover Publications, 1944. 13. McNown, J.S. Surges and Water Hammer in "Engineering Hydraulics," edited by H. Rouse, John Wiley and Sons, Inc., New York, 1950, 14. Moody, Lewis F. "Simplified Derivation of Water-Hammer Formula," Symposium on Water Hammer, ASME-ASCE, (1933) 25-28. -95

-9615. Parmakian, J. "Waterhammer Analysis," Prentice-Hall, Inc., New York, 1955. 16. Rich, George, R. "Hydraulic Transients," McGraw-Hill Book Co., Inc., 1951. 17. Rich, George, R. "Water-Hammer Analysis by the Laplace-Mellin Transformation," Trans. ASME, 1945. 18. Roberts, L. "On the Numerical Solution of the Equations for Spherical Waves of Finite Amplitude, II," J. Math. and Phys., Vol. XXXVI, No, 4, January, 1958. 19. Rouleau, W. T. "Pressure Surges in Pipelines Carrying Viscous Liquids," Journal of Basic Engineering ASME, December, 1960. 20. Rouse, Hunter. "Elementary Mechanics of Fluids," John Wiley and Sons, Inco, 1953. 21. Streeter, Victor L. "Fluid Mechanics," McGraw-Hill Book Company, Inc., 1958. 22. Streeter, Victor L. "Friction Resistance in Artificially Roughened Pipes," Trans. ASCE, 101 (1936) 681-713. 23. "Symposium on Water Hammer," ASME-ASCE, 1933. 24. Wood, F. M. "The Application of Heavisides Operational Calculus to the Solutions of Problems in Water Hammer," Transactions ASME, 59, Paper Hyd-59-15 (Nov., 1937) 707-7135