THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING TWO AND THREE-DIMENSIONAL FLUID TRANSIENTS V. L. Streeter E. B. Wylie January, 1968 IP-806

ACKNOWLEDGEMENTS The studies reported in this paper have been supported by the National Science Foundation research grant GK721 to the University of Michigan, and by contracts with the Annapolis Division of the Naval Ship Research and Development Center. iii

TABLE OF CONTENTS Page ACKNOWLEDGEMENTS............................................ iii LIST OF FIGURES................................................ v NOMENCLATURE.................................................. vii 1. INTRODUCTION............................................... 1 2. ALGEBRAIC WATERHAMMER EQUATIONS............................ 2 3. CONFIRMATION OF LATTICEWORK METHODS........................ 4 One-dimensional Wave Transmission through Frictionless Grid s...................................................... 4 Line Source in Two Dimensions................................. 6 Point Source in Three Dimensions.......................... 8 Translation of Sphere Through an Infinite Fluid............ 9 4. EXAMPLES................................................. 10 5. SUMMARY AND CONCLUSIONS................................... 12 REFERENCESo.. 0................................... 13 iv

LIST OF FIGURES Figure Page 1. Prismatic Piping Element................................. 14 2. Rectangular Latticework of Pipes with Head at Left End a Function of t Only.................................. 14 3. Plane Pressure Wave Through Rectangular Latticework...... 15 4. Equilateral Triangular Latticework....................... 16 5. Pressure Head Resulting from Plane Pressure Wave at an End of Triangular Latticework............................ 17 6. Diamond-Shaped Latticework............................... 18 7. Pressure Head Resulting from Plane Pressure Wave at a Side of a Three-Dimensional Cubical Latticework......... 19 8. Pressure Head Resulting from Plane Pressure Wave at end of Axially-Symmetric Latticework......................... 20 9. Comparison of Heads Developed Due to Variable Source Flow at the Origin in Two-Dimensional Case. Plot of Head Versus Time at Point 32 Feet from Origin................. 21 10. Plot of Head Versus Distance from Source for Three Given Times in the Two-Dimensional Case........................ 22 11. Comparison of Heads Developed Due to Point Source in Three-Dimensional Case. Head Versus Time at Two Points.. 23 12. Head Versus Distance from Point Source for Three Times in Three-Dimensional Case.............................. 24 13. Plot of Head Versus Time at Two Points for Acceleration of Sphere............................................... 25 14. Plot of Head Versus Time at (8,9) and (12,3) for Impulsive Motion of Sphere along x-axis at 50 ft/sec............... 26 15. Head Contours for Steep-Fronted Pressure Wave Passing Circular Cylinder in a Channel (two-dimensional)......... 27 16. Two-dimensional Flow of Pressure Pulse Around an Arbitrarily-Shaped Body.......................................... 28 17. Pressure Head Contours for Motion of a Cylinder (twodimensional)........................................... 29 v

LIST OF FIGURES (CONT'D) Figure Page 18. Heads Resulting from Two-Dimensional Explosion at Wall of a Liquid Filled Open Channel. Contours Given in 1000's of Feet and Time in Seconds................................ 30 vi

NOMENCLATURE a Pulse wave speed through continuum. B a/g g Acceleration of gravity. H Elevation of hydraulic gradeline (piezometric head). h Dimensionless head. Length of piping element of latticework. R Friction coefficient. r Distance from axis of symmetry, or distance from origin. t Time. V Average velocity at a section. T Constant in variable source strength equation. I Angle in Eq. (10). vii

1. INTRODUCTION A method for simulating various two- and three-dimensional transient fluid-flow cases is presented, based on the premise that the space may be represented by a latticework of piping elements, each of which obeys the one-dimensional transient flow equations. During the middle 1960's techniques have been developed to handle transient flow through distribution systems comprised of several (1)* hundred pipes ( By use of repeated configurations of piping elements, these piping systems may be extended to several thousands of pipes and may represent a plane (two-dimensional) or space (three-dimensional) continuum. The integrated (algebraic) waterhammer equations apply over the finite length of these pipes, and by satisfying boundary conditions at their junctions, transient pressure waves may be transmitted throughout the system. With the omission of certain pipes, interior boundaries are established and the translation of a pressure pulse around a body is simulatedo Similarly, by a suitable injection of source strength into the latticework, with appropriate pipes removed, the translation of any arbitrarily-shaped body through the fluid may be studied. There are certain restrictions and limitations in the methods which are discussed in connection with the various examples. The algebraic waterhammer equations are first presented, followed by examples that compare the latticework results to known hydrodynamic caseso Other examples illustrating the scope of the method are then presented. * Numbers in parentheses designate references at end of the paper. -1

2. ALGEBRAIC WATERHAMMER EQUATIONS Equations for flow through prismatic and tapering piping elements are presented, based on the use of finite difference solutions of the equations of motion and continuity by the method of characteris(2) tics ) The equations for tapered pipes are used to simulate some elements of an axially symmetric latticework. The following two equations are available to describe transmission of a pressure pulse through a pipe of length i with pulse wave speed of a, Fig. 1. C+; HPA- HB = - B(VpA - VB) - R VBIVBI (1) C-: HpB - HA B(VpB - VA) + R VAIVA (2) In these equations H is the elevation of hydraulic gradeline (piezometric head), and V is the average velocity at a section, both being functions of position and time. B = a/g with g the acceleration due to gravity; R is some suitable friction coefficient which when multiplied by the steady state velocity squared yields the head loss from B to A. The subscript P refers to some time t and the absence of subscript P on H or V denotes time t - j/ao An arbitrary flow direction is assumed, as shown by the arrow, Figo 1; the C+ equation then applies to the downstream end at time t relating Hp and VpA in terms of the previously determined quantities HB and VB for the upstream end at time t - W/a; and similarly for the C- equation for the upstream end. The absolute value signs on the resistance terms permit the equations to hold for flow in either directiono The two equations completely describe waterhammer transmission in a prismatic -2

-3piping element so long as column separation does not occur and when V is small compared with the pulse wave speed a. Latticeworks may be developed in which some of the piping elements are tapered. The differential equations must be derived for these conditions, which results in more complex forms for Eqs. (1) and (2). For a tube starting at zero area for r = 0 and increasing in area linearly with distance r, neglecting friction, the total differential equations, from the characteristics method are: - V + a V 0 (3) dt a dt r C+: dr a (4) dt /dV g dH a V dt a dt r C- d tdr- a (6) ~dt In finite difference form, by using a second order procedure to evaluate the term aV/r, Eqs. (3) and (5) become 1 VPA VB C+ VP VB H H + i PA B B A rA rB C- -HA)! Vp3 +VA) C-: Vp - VA - (HpB HA) - 2 () In these equations i — rA | - rBo1 The equations may be used in a similar manner to Eqs. (1) and (2).

3o CONFIRMATION OF LATTICEWORK METHODS Only a few hydrodynamic solutions of compressible flow cases are available for comparison with the latticework solutions. In this section one-dimensional wave transmission through various piping grids is compared with the classical solutions first, followed by a special variable source at the origin as given by Lamb(3) for the two-dimensional and the three-dimensional cases. Then the translation of a sphere, (4) given by Kirchhoff, is compared with the latticework method. Water is considered as the medium. One-dimensional Wave Transmission through Frictionless Grids In Fig. 2 a rectangular grid is shown, with the head Ho at the left end given as a function of time. If each pipe element has the length ~, with wave speed a the time for travel from one junction to an adjacent junction is At - i/a. If the system is at rest with head zero at time t = O then the head rises linearly at section 0 to 40 ft in 5 At, and is subsequently held constant, the head profile at t - 20 At is given by the light line of Fig. 3. The heavy line represents head profile for a continuum. It is to be noted that a considerable time delay occurs in the grid flow. The pulse wave travels through the grid at an apparent wave speed a' that is less than the wave speed for a single piping element or for the continuum. If there were no transverse elements in the grid this would no longer be true. The delay is caused by the fact that a wave of magnitude AH arriving at a junction of Fig. 2 is transmitted on as a wave of magnitude AH/2. This apparent change in wave speed is analogous to the change in resonant -4

-5frequency of a series piping system as compared with an equivalent simple pipe. In Fig. 4 an equilateral triangular grid is shown which has been used for several studies. If the head at the left edge is permitted to increase from 0 to 50 ft in 4 At, is then held constant at 50 ft for 10 At, then returns to zero in 4 At, the head profile at t = 18 At is shown in Fig. 5a and for t = 24 At in Fig. 5b. The heavy lines represent the head profile for the same conditions for a continuum. The points indicated by x are for the same conditions except that flow is from bottom to top of the grid of Fig. 4. This case, however, has been corrected by displacing each plotted point of Figo 5 14o to the right to compensate for the longer path. The diamond grid of Fig. 6 would transmit a one-dimensional wave of pressure from left to right through it at a speed of a 3/2, which is caused entirely by the additional length of path. A pulse AH is transmitted through a junction undiminished, as in the case of transmission through a continuum. If a one-dimensional wave were initiated along one of the diagonal lines of piping, then the waves would be retarded at junctions and an apparent smaller wave travel velocity would result, A plane pressure wave applied at a side of a three-dimensional cubical latticework produces the pressure diagram shown in Figo 7o In this case the head at the left face of the semi-infinite half space was made to vary in the same manner as for Fig. 5. The head variation at t = 18 At is shown by the light line in Fig. 7a, and for t = 25 At in Fig. 7bo The heavy lines indicate the head profile for the same conditions for a continuum,

-6An axially-symmetric latticework has been developed for several of the three-dimensional cases. It is a rectangular grid of pipes equally spaced and parallel and normal to the axis of symmetry. The transverse pipes are tapered, however, from zero area at the axis and with area varying linearly with distance r from the axis. Each pipe parallel to the axis is of constant cross-sectional area, having the area of the tapered pipe at the intersection. In effect the model is a wedge shaped half plane grid of pipes which accounts for the threedimensional flow by expansion of size of pipes with distance from the axis. Figure 8 shows the time delay due to propagation of a one-dimensional pressure pulse through the grid; again, the heavy line represents the pulse for flow through a continuum. Each of the foregoing examples represent cases of one-dimensional flow through two and three-dimensional latticeworks. Although time delays are evident, in each case the amplitude of the resulting wave is in very close agreement with the actual conditions in the continuum. The next examples treat cases of two and three-dimensional flows. Line Source in Two Dimensions The transient pressures that develop in a two-dimensional field as a result of a time-varying source can be calculated(3). The theory assumes potential flow with very small velocities. The strength of the source at the origin is described as a function of time by the equation f(t) t (9) t2 + T2 where T is a constant. The dimensionless pressure head as a function of time, t, and distance, r, is given by the relationship

-77 h = 1. + 1 a sin - - cos/2 (10) 4g H0 2rT3 4 2 where Tj is defined in the equation t = + T tan q. The gravitational acceleration is denoted by g, and Ho is the steady state pressure head on the system. Equation (10) is valid if the radius, r, to the point being considered is large compared with the quantity aT. The same problem can be handled by considering only a wedgeshaped portion of the infinite space. A numerical solution of the equations describing the unsteady flow, Eqs. (3) to (6), is obtained after the equations are placed in finite difference form using the method of characteristics. A semi-infinite container of ideal fluid is considered with Eqo (10) used as the inlet boundary condition. Calculated pressure heads at other radial distances in the model can be compared with those given by Eqo (10)o By use of a second order procedure in the numerical solution, the analysis yields the same results as Eqo (10)o The results of these computations and data from Eqo (10) are shown by the solid line in Fig. 9 at a distance of 32 feet from the origin. Figure 10 shows the pressure head as a function of distance at particular instantso For the wedge shaped container the input boundary condition was located at r = 16 feeto A latticework of equilateral triangles can be used to model the two-dimensional space and Eqo (10) serves as the input pressure head boundary condition. The grid covers a 60~ semi-infinite wedge with two-foot side lengths for the triangles. The input pressure head was utilized at 9 grid locations located along the chord 16 feet from the origin. Figure 9 shows the pressure-time pattern at a distance of 32 feet, and Figo 10 shows the pressure-distance graph at three different

-8instants of timeo It can be observed that the peak amplitudes and general form agree quite wello A delay that increases with time is also clearly recognized. Point Source in Three Dimensions The dimensionless pressure head variation in three dimensional space resulting from a point source whose strength varies with time according to the relationship given in Eq. (9) is given by the following equation(3) _ _ (E - t) h = 1. + ( 2 T (1l) gH~ 2tr [(t - a + T ] The effect of this same source in three dimensions can be modeled with a semi-infinite tapered tube of small divergence angle. A numerical solution (similar to Eqso (7) and (8)) which describes the flow of an ideal fluid in the cone yields results that agree with Eq. (11)o Equation (11) is also used as the input boundary condition. The solid line in Figo lla shows the variation of pressure head given by Eqo (11)) at a distance of 10.63 feet from the origin for the flow case defined by T = 0002 and Ho = 100 feet. Figure llb displays conditions at 15o28 feeto The same variation of pressure is obtained at these points from the tapered tube model by placing the input boundary condition at r = 2 feet. The solid lines in Fig. 12 show the pressure head variation with distance at three different instants for the same caseo The light lines in Figs. 11 and 12 show the pressure head variation obtained from the application of the source to the axially

-9symmetrical rectangular grid described above. In this case Eq. (11) was used to define the input pressure head at a radius of 7 feet. The grid results plotted in Fig. 12 were taken at random from the latticework, not just along one particular radial line. Translation of Sphere Through an Infinite Fluid Kirchhoff(4) has developed the solution for pressure and velocity in an infinite, compressible fluid due to translation of a sphere through ito As this is an axially-symmetric flow case, the latticework of wedge-shaped pipes is used. By feeding in appropriate source flow at each of the nodal points of the grid near the surface of the sphere, the boundary condition of translation of a sphere is simulated. Figure 13a is a plot of head versus time at grid point (8,9) for a 20 foot diameter sphere with center of origin accelerating from rest at 10/000 ft/sec2 along the x-axis. Figure 13b compares conditions at point (12,3). Figure 14 is a plot of head against time at points (8,9) and (12,3) for a 20 foot diameter sphere with center at origin impulsively set in motion at 50 ft/seco As in the case of the three-dimensional point source applied to this latticework, the general form of the grid results is in close agreement with the hydrodynamic solution. If a time delay correction factor were applied to the grid results, the agreement would be very goodo Some scatter is observed in the results of the latter case. This is a result of reflections at the grid intersections when subjected to a loading that contains very high frequency componentso

4. EXAMPLES In this section several examples of flow cases are given to illustrate the versatility of the method. There are no known hydrodynamic solutions of the same problems for comparison. In Fig. 15, using a rectangular grid, the two-dimensional case of flow around a circular cylinder due to a steep-fronted pressure wave of maximum head 40 is portrayed for a particular time after passage of the front of the waveo Figure 16, based on the equilateral triangular grid, shows the impact of a pressure wave on an arbitrary cylinder. The trapezoidal wave approaching has a maximum head of 50. Indexing is used in the computer program n.so that any body shape may be specified by changing the input data. Figure 17 shows the dimensionless pressure head contours at the front of a circular cylinder as a result of a constant acceleration from rest to a constant velocity of 40 ft/sec in 4 Ato The pressure distribution shown is at 6 At after the beginning of the motion. A grid of equilateral triangles with 2 foot sides was used and the cylinder has a radius of 13.84 feet. Motion of the cylinder was parallel to the base of the triangleso Figure 18 illustrates a series of instantaneous head contours for the case of an explosion at the side of a liquid filled channel 36 feet wide and 28 feet deep with free surface and with a fixed suspended bodyo When the pressure in the triangular grid is reduced to vapor pressure (indicated by dots), an interior boundary condition is imposed holding the pressure at no lower value than vapor pressureo -10

-11Since the velocity in the system is not small with respect to wave speed, the results may be distorted.

5. SUMMARY AND CONCLUSIONS A method for simulating two- and three-dimensional flow cases by use of the one-dimensional waterhammer equations in a latticework is outlined with some comparisons with hydrodynamic solutions. Since the methods are numerical, great flexibility in handling boundary conditions is apparent. Translation of arbitrary bodies through fluids, or impact of pressure waves on bodies may be examined. Free surfaces can be handled, and effects of friction may be included in the one-dimensional equations. There is a time delay in a wave front passing through a latticework, however, as compared with a continuum, and latticeworks necessarily have directional properties. Also a high frequency wave (duration of order of magnitude of At for grid spacing) is not simulated properly. The procedures are more accurate for low Mach numbers, as no interpolations are utilized in handling the algebraic waterhammer equations. Although only single-grid systems were used in the cases shown, this is not a necessary restriction, the size and nature of the grid may be changed as desired, subject only to complexity of programming and to available computer storage. Much more research is needed into the most effective latticeworks, and in ways to simulate various boundary conditions. -12

REFERENCES 1. Streeter, Victor L., "Water-Hammer Analysis of Distribution Systems," Journal of the Hydraulics Division, ASCE, Vol. 93, No. HY5, Proc. Paper 5443, Sept., 1967, pp. 185-201. 2. Streeter, V. L. and Wylie, E. B., Hydraulic Transients, Chapter IV, McGraw-Hill Book Co., Inc., New York, New York, 1967. 3. Lamb, H., Hydrodynamics, Dover Publications, Inc., 1945. 4. Kirchhoff, G., Mathematische Physik-Mechanik, 4th ed., Chapter 23, 1897. -13

-1.4--- & B Q A Figure 1. Prismatic Piping Element. o 23 57I, 0 1 2 3 4 5 6 7 8 Figure 2. Rectangular Latticework of Pipes with Head at Left End a Function of t Only.

-15LATTICEWORK HYDRODYNAMIC SOLUTION 40 - u- 20.... 14 0 5 10 15 20 DISTANCE Figure 3. Plane Pressure Wave through Rectarlgulair Latt icework.

-16Figure 4. Equilateral Triangular Grid.

-1750 / -> ~KIX.I —X X 40 H 30 - / \X HYDRODYNAMIC \ SOLUTION H3-/ | \x \ I'/ 20 I 10 \I'.4 0 0 5 10 15 20 25 DISTANCE Figure 5a. t = At 50.-'-' — _' __. 40 I, / \ HYDRODYNAMIC SOLUTION H 30 20 10 ~ " 510 15 20 25 DISTANCE Figure 5b. t = 246At Figure 5. Pressure Head Resulting from Plane Pressure Wave at an End of Triangular Latticework.

-18Figure 6. Diamond Grid.

-i960 HYDRODYNAMIC SOLUTION 40 I _ H 20/ L ATTICEWORK RESULTS 0 1....__ 0 5 10 15 20 25 DISTANCE, ft Figure 7a. t = l8At HYDRODYNAMIC SOLUTION 40/ I I/ \ LATTICEWORK. ID \ 5^ 2 RESULTS 20 0 0 5 10 15 20 25 DISTANCE, ft Figure 7b. t = 25At Figure 7. Pressure Head Resulting from Plane Pressure Wave at a Side of a Three-Dimensional Cubical Latticework.

-2060 40........ HYDRODYNAMIC H / I GRID~ \\ SOLUTION 20... 0 5 10 15 20 25 DISTANCE, ft Figure 8c. t = 18At 60 - HYDRODYNAMIC i ~ _ V \ 4r-~ SOLUTION H'I / / GRID —0 \ 20 0 5 10 15 20 25 DISTANCE, ft Figure 8b. t = 22At Figure 8. Pressure Head Resulting from Plane Pressure Wave at End of Axially-Symmetric Latticework.

-21HYDRODYNAMIC SOLUTION 2I/\ GRID i h 4 —u-5 A t — - 0 i — 0.00222.00444,00666,00889 time, sec. Figure 9. Comparison of Heads Developed Due to Variable Source Flow at the Origin in Two-Dimensional Case. Plot of Head versus rime at Point 32 Feet from Origin.

-22H YDRODYNAMIC SOLUTION t 0.00444 sec. GRID 0 HYDRODYNAMIC Y DR SOLUTION h2 _ __ _ GRID'J t =0.0080 sec. 2 DISTANCE, ft. GRIDh tO.Dim010667 sec.as HYDRODYNAMIC SOLUTION o I_-I i I.. 20 28 36 44 52 60 DISTANCE, ft. Figure 10. Plot of Head Versus Distance from Source for Three Given Times in the TwoDimensional Case.

-231.4 1.2 h GRID 1.0 __ -— ____ HYDRODYNAMIC SOLUTION o \ 0.8 0.6. 0.002.004.006 time, sec. Figure lla. r = 10.63 ft. 1.4 _ _ _ _ _ _ _ _ _ 1. \__GRID1 HYDRODYNAMIC SOLUTION \ \ 0.8 0,002.004.006 time, sec. Figure lib. r = 15.28 ft. Figure 11. Comparison of Heads Developed Due to Point Source in Three-Dimensional Case. Head versus Time at Two Points.

-241.4 1 HYDRODYNAMIC SOLUTION GRID h 1.2 HYDRODYNAMIC SOLUTION 7I,1' GRID-:<< 0.6 10 12 14 16 18 20 22 DISTANCE, ft. Figure 12. Head versus DistaYDRODYNAMIC Source for Three Times OLUTION 0.8onal Case. t= 0.003778 sec. 10 12 14 16 18 20 22 DISTANCE, ft. Figure 12. Head versus Distance from Point Source for Three Times in ThreeDimensional Case.

1400 -GRID 1200 1200 -oo GRID 1000 1000 | anR // HYDRODYNAMIC SOLUTION 800 - 800 H HYDRODYNAMIC SOLUTION 600. 600, 400 I... 400 200. 200 0 0 0.002.004.006.008 0.002.004.006.008 time, sec. time, sec. Figure 13a. Point (8,9) Figure 13b. Point (12,3) Figure 13. Plot of Head versus Time at Two Points for Acceleration of Sphere.

6000. 5000 - 4000 4000, ________ H 3000 -.0 3000.H 2000 -- 2000. GRID 4- GRID 1000 1000. -—.SOLUTION -- L SOLUTION s K 0.002.004.006 0.002.004.006 time, sec time, sec Figure 14. Plot of Head versus Time at (8,9) and (12,3) for Impulsive Motion of Sphere Along x-axis at 50 ft/sec.

2 0 -27 - 4 0 o o Figure 15. Head Contours for Steep-Fronted Pressure Wave Passing Circular Cylinder in a Channel (two-dimensional).

-28HEAD CONTOURS 20 20... * -" —-30' - $0 40'0 40 - - 50 - 1 — 7T 50 -=0 - - 40 ~70 4 0....50.. ~ ----—'-.,20 --- -- 3, - -------------— ~ —-- 10 -20 Figure 16. Two-dimensional Flow of Pressure Pulse Around an Arbitrarily-Shaped Body.

-29UNDISTURBED 1.6 1.6 1.2 2.0 2.2 2.4 DIRECTION OF MOTION Figure 17. Pressure Head Contours for Motion of a Cylinder (Two-dimens ional).

-30I0 30 2oo00 72.4 50 465 3I0 10 tu.0004 t-.0021 t-.0037 ) K) 600 10 10 4050 30 I 02 2.5 2 1. 8 30 20 40 4 ~40 40\ 0, 50 40 o0 20 40 0 20 10 30 20 t.1 "^ //N/ - ^}H- I I I0 t.0054 t.0070 t.0087 _^-0~~ 30-,"-,, 20(. 1 1 ^ ^''*-, 20^^ 20 \2\ \ V K) 26 SO/SO ~r ^ 0 /0 Figure 8. Heads Resulting from Two-Dimensional Explosion at 20 20 20Tr 10 O so_~ 1 II 1, 2s Wall of a Liquid Filled Open Channel. Contours Given in 1000 of Feet and Time in Seconds?:~,0g20 uo,o~1 0 506 ) 3 tO Id I /,, Fiur 18 ed eutngfo w-iesinlEpoina Wall~~~~~3 ofaLqi ile pnCane.Cnor Givn n 000s f eetan Tme n ecnd.: