THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING WATERHAMMER ANALYSIS WITH NONLINEAR FRICTIONAL RESISTANCE Vo Lo Streeter August, 1962 IP - 579

ACKNOWLEDGEMENT Use of the facilities of the Computing Center of the University of Michigan is gratefully acknowledgedo ii

TABLE OF CONTENTS Page ACKNOWLEDGEMENT....................... o.............. ii LIST OF FIGURES..........................o.................. iv INTRODUCTION. o o o.................................... 1 DEVELOPMENT OF EQUATIONS.................................. 2 Continuity Equation...................................... 2 Equation of Motion.................................. 4 Characteristic Equations............................... 6 Method of Specified Time Intervals.................... 9 Boundary Conditions............................... 10 APPLICATIONS.................................................. 13 Series Pipe System with Dead-End Branch................ 13 Parallel Pipe System............... 0...o............. 21 Pipe Network.................... 21 SUMMA RY................................................... 28 REFERENCES..................................... 29 iii

LIST OF FIGURES Figure Page 1 Segment of Pipe for Application of Continuity Equation............................................... 7 2 Segment of Pipe Showing Forces Exerted on Fluid in the x-Direction................................... 7 3 Characteristic Lines on the x -t Plane................ 7 4 Preselected Points on the xt - Plane for which Calculations of H and V are to be Made............... 12 5 Valve Closing Relationship........................... 12 6 Series Pipe System with Dead-End Branch. Steel Pipe... 14 7 Computer Program and Sample Calculations for Series Pipe with Dead-End Branch.......................... 15 8 Head Versus Time at Junction, Valve and Dead-End, for Series Pipe System of Figure 6. T 1 - t/tc, tc = 0.518 sec...................................... 19 9 Parallel Pipe System. Steel Pipe; Q1 = 20 cfs., 2 = 8.39 cfs o.,.o............o...................... 22 10 Plot of Head Versus Time for Sudden Valve Closure for the Pipe System of Figure 9....................... 23 11 Pipe Network for Waterhammer Calculation. Steel Pipe.. 24 12 Head Versus Time Plots for Junctions Jl, J2, and J3 for Network of Figure 11 for Uniform Valve Closure in 2 Seconds................................26 13 Head Versus Time Plots for Junctions J4, J5, and Valve for Network of Figure 11 for Uniform Valve Closure in 2 Seconds................................. 27 iv

INTRODUCTION Closure of a valve in a pipe network may cause a rupture a short time later at a point remote from the valve. Owing to the reflections from junctions and dead ends, the maximum pressure at a point in the system may not be attained at passage of the first pressure rise due to closure, and numerical methods of analysis formerly used1 23,4 are tedious to apply to complex systems, even with the neglect of fluid friction. The new methods presented here, worked out by Streeter and Lai,5 include friction and are based on the theory of characteristics with use of a high-speed digital computer. Courant and Friedrichs6 first applied these methods of computation to supersonic flow. The equation of continuity and the equation of motion are first derived in this paper. The resulting nonlinear partial differential equations are then transformed into total differential equations which are valid along certain "characteristic" lines. These equations are expressed as finite difference equations, and interpolations are used to apply the method of specified time intervals. Boundary conditions for a reservoir, a valve, and pipe junctions are then discussed, and the methods applied to three complex pipe situations: a) a series pipe system with a dead-end branch; b) a parallel pipe system; and c) a pipe network with three simple circuits. -I

DEVELOPMENT OF EQUATIONS The continuity equation and the equation of motion are developed for a segment of an elastic tube, assuming that deformations are small, as is the usual case with metal pipes. Continuity Equation The continuity equation is written for a short segment of pipe, Figure 1, of unstressed length Ax. When subjected to a change in pressure within the pipe, the walls are stressed both circumferentially and axially causing the segment to change in diameter and length. With i the Poisson ratio,`1 the axial stress and a2 the circumferential stress, the unit axial strain l1 due to change in stresses Aa1 and Ao2 is A1 -. A2 E and the unit circumferential strain 2 is Aa2 - 1. Aa 2 E in which E is Young's modulus of elasticity. The increase in volume of the pipe element is then A~1 Ax + DD 2 p2 Ax neglecting second order terms. A is the cross-sectional area of the pipe and D is the inside diameter. The extra volume that can be stored in the segment because of compressibility of the liquid is A Ax pg AH K -2

-3in which p is the mass density of liquid, g the acceleration of gravity, K the bulk modulus of elasticity of liquid, and AH is the increase in head in time At. The net inflow into the element in time At is then the sum of the last two expressions t2 pga,) A ax( +- + ) 2 K and is, from Figure lb equal to -A Vx Ax At in which V is the velocity at distance x along the pipe, and the subscript x denotes partial differentiation with respect to x. After equating the last two expressions and after substituting for ~1 and ~2 their values, 1 [Aa2(2 - 1) + Ac (1 - 2 t)] + pg l+ V At = 0 (1) E 1 K Three cases of pipe restraint2 are considered (Figure la): a) Pipe fixed at upstream end, free otherwise, A _ pg AH D AH = Pg gH D 1 4t? 2 2t' in which t' is the pipe wall thickness. After making these substitutions into Equation (1) and by writing dH/dt in place of AH/At, pg d [1 + +v (2) K dt Et - 4 x( b) Pipe anchored throughout its length, Aal A2 ACT2 - pg 2 t D 12 2 2tv

-4and pg dH [ + KD (1 - P2)] + Vx = (3) K dt Et c) Expansion joints throughout, Aa1 =0 A pgH = 2t4 and pg dH [ + KD (1 -)] + 0 (4) K dt Et? 2 x By defining c1 as the modifying expression involving the Poisson ratio (i.e., cl = - p for case a) the continuity equation may take the form dH a = dt Vx = 0 in which a2 K (6) c1 KD p (1 + 1 ) Et' a is the speed of the pressure pulse through the pipe. Now, by expanding dH/dt 2 H dx () x +t +t+ = ~ The first term may be shown to be small when compared with the second term,2 and is dropped from the equation, leaving Ht + a V 0 (8) g Equation of Motion The equation of motion, Figure 2, for the x-direction states that the resultant x-component of force equals the product of the mass

-5of the element by its acceleration, pA - [pA + (pA)x Ax] + pA Ax x - ToitDx p= A Ax dV dt T0 is the fluid shear stress at the wallo After simplifying 4Tro dV p + + = 0 X dtD dt By replacing px by pg Hx, To by pf V2/8 to introduce the Darcy-Weisbach friction factor8 f, and by expanding dV/dt fV2 gHx + - + x + Vt = O (9) The term VV may be shown2 to be small compared with Vt and is dropped from the expression, yielding gHx + Vt + -o (10) as the equation of motiono The absolute value sign is introduced into the frictional term so that it has the proper sign for reversal of flow direction in the pipe. If an exponential form of friction expression is desired in place of the quadratic term it may be substituted at any step in the derivationn If f is appropriate for the steady-state condition V = Vo, then an exponential form is f n-l _f vlvl 2D with losses varying as Vn in the exponential equation. Equations (8) and (10), labeled L1 and L2 respectively, are two hyperbolic partial differential equations, one nonlinear, that contains the continuity principle and Newtonns second law of motion, plus

-6elasticity and compressibility requirements. They are now solved by the method of characteristics. Characteristic Equations By combining Equations (8) and (10) linearly, using an unknown multiplier X, L = L2 + XL1 = X[L Hx + Ht] + [2 V ] + V f] = 0 (11) g 2D By selecting any two arbitrary (non-equal) values of X, L provides two new independent equations in V and H. However, by taking two particular values for X, L may be transformed into a pair of total differential equations. In the first set of brackets in Equation (11), by making the restriction g/X = dx/dt the quantity in brackets becomes the total derivative dH/dt. Similarly, by making the restriction dx/dt = Xa2/g, the quantity in the second set of brackets becomes dV/dt.'The two expressions for dx/dt must be identical, hence, g = Xa2 X g After solving for X, the unknown multiplier is X =+ a (12) Hence the restriction imposed on the equations is dx = + a (13) dt - The Equations (13) together with Equations (11) with X given by Equation (12) are

-7trD~' __' X t' ____ O1TDt ot x t' v v+v1Lx D P-'pPD AX I p Ox t (a) (b) Figure 1. Segment of Pipe for Application of Continuity Equation. pAxAX pA — i pA +(pA)XAX - T'1 DEX I _ AX | x —Ax -- Figure 2. Segment of Pipe Showing Forces Exerted on Fluid in the x-Direction. t P c+ C - S R Figure 3. Characteristic Lines on the x - t Plane.

-8g dH + dV + fVMVj o 1 (14) a dt dt 2D t a j (15) d-= g dH + dV + fVIVI =. s ^ + V + M o' (16) a dt dt 2D dx = -a (17) dt These four equations now replace Equations (8) and (10). Equation (14) is subject to the restriction of Equation (15)o Since a is generally constant for a given pipe, Equation (15), dx/dt = a, plots as a straight line on an x - t plane; and similarly Equation (16) is limited to lines dx/dt = -a on the x - t plane. These lines on the x - t plane are the "characteristic" lines along which Equations (14) and (16) are valido Equations (14) and (16) may be written in terms of finite differences, Figure 3, in which P is the unknown point and R and S are respectively points on Equations (15) and (17) that are known. t (Hp - H) + v. VR + (f ) (tp- t) - (18) x - a(tp - tR) = (8a) g (Hp Hs) + Vp- VS + (fV ) (tp t) = 0 (19) Xp - x + a(tp - t) = (9a) If the pipe thickness varies along the pipe, a becomes a function of Xo With conditions known at R and S (i.eo. with xR 9 tR, VR, HR, Xs, ts, VS, and HS known) the four equations permit solution for Xp, tp, Vp, and Hpo In this manner solutions may be built up from known initial conditions and from end conditions.

-9Method of Specified Time Intervals The procedure outlined for solution by the method of characteristics is not very suitable for systematic machine solution of transient 9 flow problems. By means of interpolation, one may find values of VR, VS, HR, and HS for points on the characteristic lines that go through specified points on the x - t plane, Figure 4. Consider that values of V and H are known for the evenly spaced points along the pipe, Ax apart, at time t = tc. From the values of V and H at A, C, and B, values of VR, Bi, VS and HS are computed by linear interpolation, as follows: XC xR, V VR xC XA VC VA from Figure 4. By use of Equation (18a), since xp = XC xC xR= xp =a(tp a tR) = a At Now, solving for VR V R VC + a(VA - VC) (20) in which 0 is At/Ax, the ratio of preselected time increment to distance increment. Similarly HR = HC + ~ a(HA - HC) (21) VS = V + 9 a (VB V) (22) H = H + ~ a(H - H) (23) S C B C For convergence of the solution At should be less than Ax/a. By use of the interpolated values, only two equations are needed, Equations (18) and (19) to solve for the two unknowns, Hp and Vpo

-104 The friction terms in Equations (18) and (19) may be evaluated for conditions at C without appreciable error, if increments Ax and At are kept sufficiently small. After solving Equations (18) and (19) for Hp and Vp HR + HS _ Hp = + (VR VS) (24) 2 2g VR + Vs - _ c VIV I Vp + (HR - HS) c ci t (25) 2 2a 2D These equations, together with Equations (20) to (23), permit computation of velocity and head at interior points in the pipe. Special relations for computation of end points are needed. known as boundary conditions, and are discussed in the following section, Boundary Conditions At the end of a pipe, only one of the two Equations (18) and (19) is applicableo At the upstream end (x = 0), the interpolation Equations (22) and (23) are appropriate, and with Equation (19), one equation is available in the two unknowns Hp and Vpo An auxiliary equation is needed that specifies Vp, Hp or some relation between them, so that two equations in two unknowns are available. The simplest case is probably that of a reservoir, where Hp is constant, or a function of time. At the downstream end of a pipe Equations (18), (20), and (21) are available and an auxiliary relation is needed to solve for Vp and Hp. The relationships between Vp and Hp for a valve are worked out to illustrate the procedures for determining boundary conditions.

-11N The valve is treated as an orifice. With the steady-state head loss Ho across the valve and the steady-state velocity Vo in the pipe, the orifice equation becomes VoA = (Cd AG)o 2g o (26) with A the pipe cross-sectional area, Cd the valve coefficient and AG the valve opening. In general vp = Cd AG r (27) After dividing Equation (27) by Equation (26) and substituting T = Cd AG/(Cd AG)o, the dimensionless valve opening VP - T H (28) Vo Ho In this dimensionless equation T varies from 1 to 0 for closure, in a manner depending upon the rate of closure., T is usually expressed as a function of time, as shown in Figure 5. Equation (28) provides the second relation between Vp and Hpo By eliminating Hp in Equations (18) and (28) VP - F;i 1 + 4g Ho (g R +, R c VcVcI t) i' (29) 2g L aT2 V a 2D In the solution of the quadratic in Vp thI positive sign was taken before the radical so that Vp = VO for T = 1 in the steady-state case. After Vp is determined from Equation (29) Hp is found from Equation (28)o For the special case T = 0 V = = 0 and Hp is found from Equation (18)o In general, boundary conditions are easy to formulate. Special cases are included in the applications to complex piping situationso

-12t At T- _ 7c - - _ - tc R C S At I x I LAX — ^1 Figure 4. Preselected Points on the xt-plane for which Calculations of H and V are to be made, t Figure 5. Valve Closing Relationship.

APPLICATI ONS Three applications are presented that show the general adaptability of equations and the boundary conditions. The Poisson s ratio effects have not been specifically taken into account, since its inclusion affects only the speed of pressure wave, and the pipe thicknesses have been arbitrarily assumed for the exampleso First a series pipe system with a deadend branch is taken up, followed by a 4 pipe system with two of the pipes in parallel, and concluding with a network having 9 pipes and 3 simple circuitso Series Pipe System with Dead-End Branch In Figure 6 a series pipe system with a dead-end branch, is taken for the first exampleo The valve closure is T l1 - tin which tc = 0o518 seco is the time of closure of the valve It is the round-trip travel time of a pressure pulse in pipe lo The computer program, written in MAD (Michigan Algorithm Decoder) language, together with a sample of the calculations, is shown. in Figure 7~ A plot of head versus time for the dead end of pipe 2, the junction of pipes 1 and 3, and the valve are shown in Figure 80 Constant friction factors have been assumed [F(l) ooo F(3)], which results in V losses. The steady-state solution is first obtained with. T = 1, and the head and velocity at 11 equally-spaced sections along each pipe are computed, stored and printed for time t = 0. The time is then incremented (DELT) by an amount such that the interpolation points R and S lie within A and C, Figure 4, for all three pipeso The interior points are next calculated, -13

-14605.68 1500'r30"D. 1000'- 36"D. ~ ~ f.0.03 t' 0.72 f0.025 VALVE, Ho-600' t' 0.72" Q =20cfs Figure 6. Series Pipe System with Dead-End Branch. Steel Pipe.

__ —-------— "'M —R ENVY 3t ^ -STAT,~^,E TO.S l ----— " *s —a'.^ ^^~ —--— H p l- 3 - -- ==VECTOR VALUESV(3 31- 00 DA:r.... —. — (9- 31- -, i DT t-(( 3 1.....-...... ATA DEHTI2t __ _......... 2- 3) _ VE RVLE I= GI —-- VEN DAT LAD I TS IN -. r A- - - ------ EEND OF C CO'NDIT ------ I1A I.. 1 ET9LA()/ -]-N ~NS..W............. -— RR-'"=RR:_____________2FO R_-10 R_ ---—........... = —L(I ( - LATHROUGH IFOR _Z_1________._- ~ — RcTj —-— "C-1-T2UE. AND AEPIT C_.OI-— WRE~EET...... INC, —It P E _ AU _ 2 VHR -,i F8 -. --- pRIT oMP ENU' _TA 1 um v THROUGH Z9 tF. R J=lti J.. a3

-16TH —---------- ROUGH+-Z9.-FOR i —, I ----- ----------- VR=V4J, I)+TH(J)*A(J)*(V(J I-1 )-V(J I)) HR<HJJ-T+THlT) *A CjJi* (H{JjTPIVHTtu W -RTJ_ __ _J _ HS=H(J I)+TH(J)*A(J)*(H(J +1)-H(J I ) -— ~ —-— 1 —— V'~; V$S-Vi-J +- H ( Jl:A (-J *-CV(' -JiI-I-) —V: i --- -- - - ---- - --- -- -- VP (J I)t =(VR+VS)/2 +16 1*( HR-HS)/A(J)-F(J)*V{(J )*ABS4(V-(Jt i) ---------------- 2)-*DELT/T2.*.D(J)') -. Z9 HP(J ) =A(J)*(VR-VS)/64.4+(HR+HS)/2. RBOUNDARY CONDITION AT RESERVOIR -V - - V(-3-0)-+TH-( 3')*A 3 )* (V-(-3 11-V3(O - - -. —--- HS=H(3tO)+TH(3)*A(3)*(H(31)-H (3.0)) - VP( 3,O)=VS+G*(H<(30)-HS)/A(3)-F(3) *V3t0)*.ABs.V(3,0)*DELT/(2 _ 2*D(3)) R RJUNCTION BOUNDARY CONDITION --- - VS=V(10,)+TH( 1)*A(l)(V(l1l)-V{,0I}) HS=H(1,O)+TH(1)*A(1)*(H(1{1)-H(1{0)) - VR=V(3,N)+TH(3)*A(3)*(V(3,N-1)-V(3,N)) HR=H(3,N)+TH(3)*A(3)*(H(3,N-1)-H(3 N)) C2=VR+G*HR/A(3)-F(3)*V(3,N)*'ABS.V(3,N)*DELT/(2,*D(3))..HP(3-,N)=(C2-Cl*D(l)*D(1)/(D(3)*D(3)))/(G/A(3)+G*D(l)*D(1)/(D( 23)*D(3)*A(1))) HP(-10t) =HP(3 N)'VP ( 3N)=C2-G*HP(3iN)/A(3)..............'-'- I -- -"C'i-'-+G*HP' ( 3 N-) /A(i R -RD-EAD- END —BOUNDAR Y COND IT ION VP(2#N)-Oe.....VR=V{2,N)+TH(2)*A(2)*(V(2,N-')-V (2-'N')r-} -- _HRH(2,N)+TH(2)*A(2)*(H(2,N-1)-H(2,N)) H — P(2-N)-=H R +A(2) *VR/G.. R --------— BO-UNDARY —CONDITI ON AT-VA-LVE WHENEVER T*LE.TCTAU=1.-SQRT.(T/TC) -W —" —V-W-HENEVETR- - TCG TTAU ------------------ VR=V lN)+TH(1)*A (1)*(V( 1N-1)-V(1 N)) HR=H ( N+THTl 1 N-1) -H ( N )Ti) VS$V(2,0)+TH(2)*A*2)*(V(2 1)-V\(20)) HSH( 2.0)+TH (2) *A( 2- *H( H(2 )-H (20)) C l=VS-G*HS/A(2)-F(2)*V(2 0)*. ABS V (2 0)*DELT/(2. *D(2)) -._ _ _ _ _ __ — -_ — _ _ -________________ _______ _............ _ _ _ _ _ __ _ _____ C2=G/A( 2) C3=VR+G'*HR/A (1) -F( 1).*V ( 1,N ) *.ABS V (1 N ) *DELT/( 2, * D (1) ) C4=-G/A(1) C6=(C3-Cl*D(2)*D(2)/(D(1)*D(l))l)/(C4-C2*D(2)*D(2)/(D(1)*D(l)) 2) WHENEVER TAU*E.Oe0.HP(1iN)=-C6 WHENEVER TAU.GeO.e C5TAU*V7( 1) ((C4-C2*D(2)*D( 2 )/(D( 1 )*D( 1) ) *SQRT. (HO)) C7=SQRT ( C5*C5/4.-C6) ___ _ _HP{ __N)={ C5/2_.+C7)(C5/2.C7) _ _ ______ __ END OF CONDITIONAL HP(2{:.O)=HP( l#N) ____VP<2^0)C1C2_HP(ltN)__l________. - ________ ___._ ~................................... VP( 1.N)z=C3+C4*HP( 1N) R~EP( OFL)=CT+C2NH P I R(lMAN) R RSTORA-6E OF PLOTTING INFORMATION

-17---------------------------------— _-artrrt -- - ---------------- ------------------------------- ---------- ---------------------------------- HJ(U)=HP( 3N) HV(U)=HP( 1,N) THROUGH Z0lFOR Jlolt1J.G*3 t —----------------------------— R ZT R —--------------------------------------------- -- ------------------- --- H(J.I)=HP(JI) ---------------— z —------— r. - ---- -- --- -------------------------------------------------------------------------- WHENEVER U/5*5.E.UTRANSFER TO Z3 TRANStER TO Zb PLOT PRINT RESULTS HJ* *HJ( W-1 )HV.*,HV(W-1 ) -----------------— _ —--- --— t — pT — cOm _T$T -. —-—. ----------— _-^A-A-A -- - -- --- - - - --— A —— A-VAV —----------- --- EXECUTE SETPLT (1lTTtHVWS$*S,24tORD) ----------------------------— VECTR - -UE-ORD-s ----- - ------ ------ PRINT COMMENTSO TIME IN SECONDS$ PRINT COMMENTS1 HEAD AT JUNCTION BETWEEN PIPES 1 AND 2 3 VS TIMES -----— X7 —---- ---------------------— EXCE-W $ - ------------------- ---- -------------- PRINT COMMENTSO TIME IN SECONDS$ ------------------------------- T — ------------------- -------------------- -- END OF PROGRAM SDATA L(1)s1000. 2000. 1500. D(1)=3. 2 *2.5F( 1)=l025*,03 t03 THI( 1 )06.05~ Q( 1)=20*0,,20.* GIVEN DATA Q(l) = 2C.000000, F(1) =.025000, F(2) =.030000, F(3) =.030000 L(1) = 100C.000000, L(2) = 2000.000000, L(3) = 1500.000000, D(1) = 3.000000 0(2) = 2.000000, D(3) = 2.500000, THI(1) =.060000, THI(2) =.050000 THI(3) =.060000, HO = 600.000000, RHO = 1.935000, E = 4.320000E 09 K = 4.32COOGE 07, N = 10, W = 300, G = 32.200000 CON = 1.000 CC HRR(O) = 605.675774, DELT =.025921, HF(1) = 1.035920, A(1) = 3857.942596 A(2) = 3993.349976, A(3) = 3969.790161, LA(1) =.259206, LA(2) =.500833 LA(3) =.377854, TC =.518411 _______ _______ ______ _______ _______ ______ _______ _______ ______ -_ _ _ _ _ _ _ _ _ - -- --— _-_____________ ______________,_______________ _________________ —- --- -- - ---- __

-18VELOCITIES AND HEADS AT TENTH POINTS FOR' EQUAL TIME INCREMENTS'IPE TIME TAU X/L=.0.1.2.3.4.5.6.7.8.9 1. 1.000 1.OOC H= 601.036 600.932 600.829 600.725 600.622 600.518 600.414 600.311 600.207 600.104 600.000 1 V= 2.829 2.829 2.89 T 2-.8292.89- -.829.2 —-- 2.9 2.82...829 2...-. 9 2..-2 — — 2 —--.00 100- o00 H= --—.)- -. —"-6007ooo60. 6- --- — 6 6C — 6-~-? —-.-6 -5-: —6~-.-" —D6 U~ —2 V=.000.000.000.000.000.000.000.000.000.000.000 3.000 1.000 H= 605.676 605.212 604.748 604.284 603.820 603.356 602.892 602.428 601.964 601.500 601.036.-^-3 V: 4.074 4. 074 - 4.074 4.074-?~ 4.~-074- 4.0..4 4.0774.-7- 4.0 74- 4.7OT4 -.4.0-74)-7,- 4.0-74 — T".130 " 500 H =601-0~6 — — 6. 93-2 —600.-829 -6-0-57-72 —6-6- — G.T-6-5 —--- 5-6 —1 V= 2.829 2,829 2.829 2.829 2.829 2.829 2.445 2.277 2.144 2.029 1.926 2.130.500 t= 708.384 680.888 647.781 618.438 603.307 600.000 600.000 600.000 600.000 600.000 600.000 - ----------------------------- 2 —---------- -------------—. —5 —----- —. — -.-^-_-__ —---—._U —------ O —% —- --— G —----— 0 —----— 0 —------- --- 3 ".130.~- 500 -H-= 65.67-6 605-.212 604.748 6.864-.8 3.V 4.074,4.074 4.074 4.074 4.074 4.074 4.074 4.074 4.074 4.074 4.074 1.259.293 H= 601.036 647.022 667.017 682.840 696.505 708.776 720.082 730.637 740.622 750.121 759.237 V- 2.829 2.445 2.277 2.144 2.030 1.927 1.832 1.743 1.659 1.57..-03 —-,^ —--------------------------- ------ --- -_^^^'2.259.-293 H = -75-9. 237- 740.537 ~7T9.3 —8 9-4. —T~ —-66. —-T 6-6- Z -r —5T.r — 6-U. -- 2 V= 1.283 1.133.962.763.538.318.147.049.010.001.000 3.259.293 H= 605.676 605.212 604.748 604.284 603.820 603.356 602.892 602.428 601.964 601.500 601.036 3 V= 4.074 4.074 4.074 —-- 4.0-45774.-7 — -...-T7-4 —- -.- 4..4; —-— 4.7_ —-!-7-.389-7.-i-4 —------- H — -7-2T6i —?73-9) —.-.zf55-75l8~ —75 7- -—. —T2-T —?~-.-6TE-s T — I V= 2.103 1.988 1.877 1.768 1.655 1.504 1.430 1.359 1.290 1.223 1.157 2.389.134 H= 800.841 785.188 768.527 750.440 730.295 707.540 682.453 656.968 634.625 619.339 613.918 2 V= 1.618 1.492 1.358 1.212 1.050.867.665.458.271.121.000 -^ —- ---------------------— " —- ----— " — -------------— ^^ —-------------— ^ —--- 3.389.134 H= 605.676 605.212 604.748 604.284 603.80 603.356 6 8 644.641 679.030 707.559 730.216 3, V= 4.074 4.074 4.074 4.074 4.074 4.074 3.975 3.732 3.450 3.215 3.029 1.518 -.000 H= 790.899 797.850 804.444 810.689 816.606 822.177 827.388 832.164 836.385 839.659 837.773 1 V= 1.764 1.675 1.588 1.503 1.418 1.334 1.250 1.165 1.078.985.851 2.518 -.000 H= 837.773 823.713 809.114 793.844 777.786 760.932 743.634 726.886 712.448 702.512 698.953 2 V= 1.914 1.801 1.684 1.561 1.429 1.284 1.114.907.649._ 341.000 3.518 -.000 H= 605.676 607.C63 615.390 634.504 661.826 690.746 716.613 738.587 757.677 774.910 790.899 3 V= 4.074 4.059 3.988 3.829 3.604 3.366 3.154 2.972 2.8'5 2.672 2.540 -- ---- ------ --------- - ^; — - ----- ~ ~ ~ ~ ~~ ~ ~ ~ ~.648.000 H= 840.615 846.677 852.534 858.209 863.697 869.013 870.768 870.441 869.701 868.694 867.457 i ______ __ V= 1.487 1.410 1.334 1.259 1.185 1.111 1.066 1.037 1.009.983.956 2.648.000 H= 867.457 859.696 848.976 836.197 824.739 817.013 813.305 813.113 815.225 817.706 818.774 2 V= 2.150 2.084 1.987 1.853 1.695 1.510 1.286 1.015.703.359.000.3648 0OC H= 605.676 653.004 693.180 723.680 746.667 765.51 782.437 798125 812932 827054 840.61 -J______3 _ V= 3.352 3 315 3 215 3,071 2.914 2.763 2.623 2.493 2.370 2.253 2.141 Figure 7. Computer Program and Sample Calculations for Series Pipe vith Dead-End Branch.

-191.200 DEAD END 1000 VALVE / JUNCTIONM \ \ \\\ I //V 400 200 0 2 3 4 5 6 7 8 TIME, Sec. Figure 8. Head Versus Time at Junction, V _alve and Dead End, for Series Pipe System of Figure 6. T = 1 -.4t/tc, t 0.518 sec.

-20then the various boundary conditionso At the reservoir the head [H(3,0)] is held constant, and the velocity [VP(350)] is calculated from Equation (19)o At the junction between pipes minor losses are neglected, and the conditions to be satisfied are continuity and equality of head (this neglects change in velocity head, as is customarily done). At the deadend the velocity is zero [VP(2,N) = 0], and the head is computed from Equation (18)o The friction term drops out of statement 067 (Figure 7) since V = 0. At the valve, the boundary conditions are continuity and equality of head. For continuity V1 A1 V2 A2 CD AG 2g H and for steady flow Vo A1 (CD AG)o 2 o After dividing the first equation by the second equation V V A'7 1 2 2 = H Vo V0 A1 Ho By substituting for V1 and V2 from Equations (18) and (19), a quadratic in fH is obtained. The proper root of the radical term in the quadratic solution must be selected so that whenf H is squared it equals Ho in the limit for steady flowo This completes the calculation of all points for this time. These values are stored temporarily, and after 5 time increments, the results are printed. Figure 8 illustrates the complex variation of head with time due to the reflections from the various boundaries.

-21Parallel Pipe System A parallel pipe system, Figure 9, is the second example of waterhammer in a complex circuito Boundary conditions must be worked out for the reservoir, the upstream and downstream junctions, and the valve At the reservoir the head remains constant, so Equation (19) provides the value of Vp after VS and HS are computed. The junction boundary conditions are continuity and a common head for the 3 pipeso Forthe downstream junction Equation (18) is applied twice, to pipes 2 and 35 and Equation (19) to pipe 4. With continuity this yields 4 equations in 4 unknowns, Vp, Vp, Vp, and Hpo A similar procedure applies to the upstream junction. For the value Equation (30)may be applied when T > 0. For T = 0 Vp = 0, and Hp is given by Equation (18). Owing to space limitations, the program is not presented, but head versus time for the two junctions and the valve are presented in Figure 10,9 with initial head at the valve Ho = 500 fto and for instantaneous valve closureo Pipe Network As the final example of flow through a complex system, a network, Figure 11, is selected having 9 pipes and 3 simple circuits. A valve at the downstream end is closed uniformly in twice the round-trip travel time of the pressure wave through pipe 9% The steady-state flow through the system is first determined by the Hardy Cross methodl1 of distributing flow in a network, and the results are shown in Figure 1

-22503.62 1200'-24"D f=0.03 1000'-36" D. J. At 70.025' /1200 -39'D. f 0.024 J2 00 f= 0.02'/(2) // t'0.03 2100'- 27.6"D f =0.018 t'0.02 Figure 9. Parallel Pipe System. Steel Pipe; Q1 20 cfs., Q2 = 8.39 cfs.

1000.. 800 1 VALVE VALVE DOWNSTREAM 0 / 1 2 3 4 JUNCTION 600 "'*: / \LL I UPSTREAM d- JUNCTION 400 \\ / / 200 0 2 3 4 5 6 7 8 9 10 TIME, Sec. Figure 10. Plot of Head Versus Time for Sudden Valve Closure for the Pipe System of Figure 9.

626.64' 0 o; 0 ^< t'= 0.035" ) Figure 11. Pipe Network for Wterh er Calculation. Steel Pipe. ~o,':oo \s00'-36". f =0.0 17.96 cfs -- = 30.cfs tf 0.055" t'=7. 075 ^ LH ~ Figure 11. Pipe Network for Waterhammer Calculation. Steel Pipe.

-25In writing the watehammer program double subscripting is used, i.eo, V(3,4) is the velocity in pipe 3 at the fourth section. The velocities and heads were computed for 9 equally-spaced sections in each pipe (ioe, 250 fto apart in pipe 1i and 375 fto apart in pipe 2), and the time increment At = 0.04695 seco was selected on the basis of pipe 8 calculations so that At does not exceed Ax/a in any pipe. The boundary conditions at the reservoir, the valve, and the 5 interior junctions were satisfied for each time increment. The head was held constant at the reservoir, and at each junction continuity was satisfied and the head made the same for end sections of each pipe entering or leaving the junction. The program is not included to conserve space, but plots of head versus time at the 5 interior junctions and the valve are shown in Figures 12 and 13o

1200.......... -26 Soo 200___ 0 2 4 6 8 10 12 TIME, SEC. Figure 12. Head Versus Time Plots for Junctions J1, J2, and J3 for Network of Figure 11 for Uniform Valve Closure in 2 Seconds. of Figure ll for Uniform Valve Closure in 2 Seconds.

-271200 ALVE 10ooo00 1 VALVE / / 200 200 0 Fi 2 4 6 8 10 l 12 TIME, SEC. Figure 13. Head Versus Time Plots for Junctions J4, J5, and Valve for Network of Figure 11 for Uniform Valve Closure in 2 Seconds.

SUMMARY Equations are derived for the numerical solution of transient liquid flow in a conduit with nonlinear friction by use of the theory of characteristics and the method of specified time intervals. Boundary conditions for various cases have been discussed to show their ease of formulation. The methods are applied to 3 complex pipe systems: series, parallels and network, using an IBM 709 high-speed digital computer to generate solutions which are presented as plots of head against time. -28

REFERENCES lo Allievi, L,, "Theory of Waterhammer," translated by Eo Eo Halmos, printed by Richardo Garoni, Rome, Italy, 19250 2o Parmakian, Jo, "Waterhammer Analysis," Prentice-Hall, Inc, New York, 1955o 35 Angus, RQ W., "Water-Hammer Pressures in Compound and Branched Pipes," Trans ASCE, Volo l04, ppo 340-401 (1959) 4. Rich, Go R., "Hydraulic Transients,', McGraw-Hill Book Co., Inco New York, 1951L 5. Streeter, V. Lo and Chintu, Lai, "Water-Hammer Analysis Including Fluid Friction," J. Hyd. Divo ASCE, Volo 88, No. HY3, May, 1962, Part Io 6. Courant, Ro, and Friedrichs, K. O., "Supersonic Flow and Shock Waves," Interscience Publishers, New York, 1948o 7. Timoshenko, S., Strength of Materials, 2nd Edo, Part 2, Do Van Nostrand, Inc., New York, 1941. 80 Streeter, VO Lo, "Fluid Mechanics," 3rd Ed., McGraw-Hill Book Company, Inc., New York, ppo 211. 214-222, (1962). 9. Lister, M., The Numerical Solution of Hyperbolic Partial Differential hEquations by the Method of Characteristics in "Mathematical Methods for Digital Computers," edited by Anthony Ralston and Herbert So Wilf, John Wiley and Sons, Inc., 1960. 10. Cross, Hardy, "Analysis of Flow in Networks of Conduits or Conductors, University of Illinois, Bull. 286, November 1946. -29