THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING FREQUENCY DEPENDENT FRICTION IN TRANSIENT PIPE FLOW Werner Zielke A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in the University of Michigan Department of Civil Engineering 1966 December, 1966 IP-762

Doctoral Committee: Professor Victor L. Streeter, Chairman Professor Glen V. Berg Professor Ernest F. Brater Associate Professor Walter R. Debler Associate Professor E. Benjamin Wylie ii

ACKNOWLEDGEMENTS The author wishes to express his gratitude to the members of his doctoral committee for their assistance, and in particular to Professor V. L. Streeter, chairman of the committee, under whose guidance and supervision the work was carried out, and to Professor E. B. Wylie for the interest provided in many discussions. The author is indebted to the University of Michigan for granting teaching and tuition fellowships and to the National Institutes of Health, and the National Science Foundation for their financial support of the work. Lastly, the author wishes to thank the staff of the Industry Program of the College of Engineering for the preparation of the manuscript. iii

TABLE OF CONTENTS Page;. AC KNOWLEDGEME N IS................................................. iii LIST OF TABLES........................................ v LIST OF FIGURES...o........o........,...........................o vi NOMENCLATURE.......................................................viii I INTRODUCTION.................................................. 1 II LITERATURE............................................... 8 III FUNDAMENTAL EQUATIONS OF WATERHAMMEER......................... 16 IV WALL SHEAR STRESS IN UNSTEADY LAMINAR FLOW................. 23 A. Fundamental Equations......,..o......................,. 23 B. Solution for Periodic Flow........................... 27 C. Solution for Arbitrary Flow........................... 33 V APPLICATION OF FREQUENCY DEPENDENT FRICTION TO STEADY OSCILLATING FLOW.........,................................. A. Method of Solution...................................... B. Experimental Verification................................ C. Numerical Example............................... VI APPLICATION OF FREQUENCY DEPENDENT FRICTION TO TRANSIENT FLOW CASES.....o................................... A. Method of Solution.............................o... B. Experimental Verification............................... VII APPLICATION OF FREQUENCY DEPENDENT FRICTION TO THE FLOW OF BLOOD IN ARTERIES................................o.......o A. Basic Assumptions....................................... B. Computed Results..................................... VIII SUMMARY AND CONCLUSIONS................................. APPENDIX DISCUSSION OF THE FUNCTION (C T-i), a REAL........... BIBLIOGRAPHY..........................,.............. o o o 44 44 49 51 56 56 63 67 67 71 77 81 84 iv

LIST OF TABLES Table Page I T and 5 for a = 0.1, 0.2,... 10.0................... 30 J0(x) II The First 100 Roots of the Equation x - 2 0...... 56 J1(x) IIIA The Weighting Function W(T) with T Varying from T = 0.0001 to T = 0.0100 in Steps of 0.0001............ 41 IIIB The Weighting Function W(T) with T Varying from T = 0.0105 to T = 0.1000 in Steps of 0.0005.......... 42 IIIC The Weighting Function W(T) with T Varying from. T = 0.1010 to T = 0.2000 in Steps of 0.0010.......... 43 v

LIST OF FIGURES Figure Page 1. Head Fluctuation at Valve after Sudden Closure. Experimental Work of Holmboe. 11........................ 4 2. Ratio of Pressure Amplitudes Downstream/Upstream versus Frequency. Expeimental Work of A. F. D'Souza and R. Oldenburger..................................... 4 3. t and 5 as Functions of a o*oO........................ 31 4. Ratio of Unsteady and Steady Friction as Function of a.. 31 5. Weighting Function W as Function of T................. 40 6. Experimental System for Study of Pressure Amplification Due to Resonance....................5........o.......... 50 7. Ratio of the Amplitudes Between Downstream and Upstream Pressure versus Frequency................. o.............. 50 8. Hypothetical Pipe System for Numerical Study.............. 52 9. Assumed Head Variation at Point A of the System in Figure 8.............................. oo.............. 52 10. Head Variation at Point B Resulting from a Square Wave Input of Head at Point A..........O..................... 54 11. Discharge Variation at Point A Resulting from a Square Wave Input of Head at Point A........................ 55 12. Characteristic Lines on the x - t Plane.................. 58 13. Grid of Characteristics for Specified Time Intervals..... 58 14. Set-up for the Experimental Work of Holmboe(ll)....... 64 15. Head Fluctuation at Valve and Midstream after Sudden Valve Closure. Comparison with Calculated Results Using Steady State and Freluency Dependent Friction. Experimental Work of Holmb oe 11............................................ 66 16. Assumed Flow Variation in a Rigid Tube of 0.316 cm Diameter........................................... 72 vi

LIST OF FIGURES (CONT'D) Figure Page 17. Head Loss Due to Friction for the Flow Variation of Figure 16............................................... 72 18. Head Variation at Proximal and Distal End of a Tapered Artery....................................... 75 19. Flow Variation at Proximal End of Tapered Artery Resulting from the Head Variation of Figure 18..................... 75 vii

NOMENCLATURE A Pipe cross section area a Wave speed C Capacitance per unit length C1, C2 Constants D Pipe diameter e e = 2.71828 F = ^/p 3x f Darcy-Weisbach friction factor g Acceleration due to gravity H H = H(x,t), piezometric head hf Head loss due to fluid friction per unit length i i =~-1 i Integer JO Bessel function of first kind and zero order J, Bessel function of first kind and first order 71 pl (Z) = z Jo(z)/Jl(z), modified quotient of Bessel functions j Integer counter K Bulk modulus of fluid k Integer L Inertance per unit length 2 Length of pipe No Bessel Function of second kind and zero order P Amplitude of pressure fluctuation P Subscript denoting point P p p = p(x,t), pressure viii * * Vlll

NOMENCLATURE (CONT'D) Q Q = Q(x,t) discharge R Radius of pipe R Linearized resistance per unit length (only in Chapter V) R Subscript denoting point R Re Re = VD/v, Reynolds number r Distance from pipe axis S Subscript denoting point S s Complex variable s Roots of l (i R) - 2 = 0 sj = v R2 s Subscript denoting steady t Time u Time, used in convolution integral u Subscript denoting unsteady V V = V(x,t) = Q(x,t)/A, mean velocity v v = v(r,t), velocity W W = W(T), weights for past velocity changes x Distance along pipe axis x Real variable y Real variable z Complex variable Ca a=RTcn/v, dimensionless frequency 6D p ff= - aD-, Ytaper of pipe r Gamma function Y Specific weight of fluid ix

NOMENCLATURE (CONT'D) 5;See Equation (65) See Equation (22) Tj Roots of 1(Z) - 2 = 0 X See Equations (68) and (69) ]ii] = pv, dynamic viscosity v v = U/p, kinematic viscosity t See Equation (22) j Roots of JO(z) = 0 iT I = 3.14159 p Density of fluid V T T = - t, dimensionless time R To Shear stress at pipe wall See Equation (31) 3O Circular frequency x

I. INTRODUCTION Flow through pipes is always subject to energy loss due to fluid friction. While there is a well developed theory about the energy loss in steady pipe flow, little knowledge exists about the energy dissipation in unsteady flow systems. Although there is experimental and theoretical evidence that the fluid friction is highly dependent upon the rate of change of velocity, most investigators apply the steady state friction term to the analysis of unsteady flow and neglect the evidence that higher harmonics attenuate much faster than low frequency components of waves. The treatment of transient flow in pipes, commonly called waterhammer, begins historically with the well known works of Joukowsky(l5) and Allievi(l) who treated the propagation of pressure waves in frictionless systems, described by the wave equations 6H 1 V ax + - = O E gand aH a2 6V + g ~ - = 0 where H and V are dependent variables piezometric head and velocity, t and x are independent variables time and distance along the pipe axis, a is wave speed, and g is acceleration due to gravity. Later investigators included friction effects by applying the friction terms derived for incompressible steady laminar or turbulent flow, -1

-2 h = 32v gD2 and 12 hf = 1 V D 2g where hf is the friction loss per unit length of pipe, v is the kinematic viscosity, D is the pipe diameter, and f is the DarcyWeisbach friction factor. This approach assumes that the energy dissipation is only a function of the pipe constants and the sectional mean velocity. The wave equations including fluid friction are aH 1 av +- + g + hf = 0 and aH a2 av T + - = 0 Only when hf is a linear function of V, such as when flow is laminar or when the turbulent friction term is linearized, is an analytic solution of the equations possible. It has been studied in exhaustive detail in connection with the resistance-capacitance-inductance model of electric transmission lines and applied to waterhammer problems to calculate the transient and frequency response of fluid lines. A way to handle nonlinear friction in hydraulic lines is by lumping it into so called friction joints. However the graphical waterhammer theory permits the employment of only a very small number of joints, while the method of characteristics applied on a computer offers a procedure that can treat distributed nonlinear friction owing to the great number of sections which can be considered.

-5 The application of steady state friction to transients, both laminar and turbulent, implies that the damping effects are proportional only to the velocity. Since no dependency on the rate of change of velocity is taken into account, no change of shape of the waves can show up in the theoretical results. That means that discontinuities are propagated as discontinuities and higher frequency disturbances are subject to the same damping as low frequency disturbances. Two significant experimental results by Holmboe(11) and D'Souza and Oldenburger(8) are chosen to show the limitation of this assumption: 1. A pipe of one inch diameter and a length, 2, of 118.4 feet is connected upstream to a reservoir with constant head. The working fluid is an oil with a kinematic viscosity, v, of v = 0.000427 ft2/sec, almost 50 times that of water. The wave speed is 4345 f.p.s. The downstream valve is suddenly closed and the resulting head fluctuation is recorded in Figure 1o The step change in flow at the valve results in a sudden increase of the pressure at the valve and sends a steep pressure pulse through the system. After the relatively short wave travel time 22/a, when the reflected wave returns to the valve, a significant rxund-off effect of the leading edge of the pulse is evident and increases as time progresses. This effect cannot be predicted using a steady state friction term in the waterhammer analysis. Furthermore the attenuation of the fluctuation is greater then predicted in the conventional waterhammer theory.

-4 0 w I O H (\ KM I u <I <I TIME Figure 1. Head Fluctuation at Valve after Sudden Closure. Experimental Work of Holmboe. (11) 3.0 - 2.6 1 0 O O 0 2.2 I: 1.8: 1.4 C,) cn a-.6 0 0 0 0 0 0 0 0 0 0 0 00 0 0 II I 1 1 1 1 I I 0 10 20 30 40 50 60 70 80 90 100 FREQUENCY, CYCLES PER SECOND Figure 2. Ratio of Pressure Amplitudes Downstream/Upstream versus Frequency. Fperimental Work of A. F. D'Souza and R. Oldenburger.

-5 2. An experimental study of the frequency response of a pipe proves that the resonating pressure amplitudes of the higher harmonics are smaller than the lower harmonics. An experiment over a wide range of harmonics is not available, but Ao F. DISouza and R. Oldenburger (8) give the pressure ratio between downstream and upstream end of a 40o25 feet long steel pipe with an 0.495 in ID for a frequency range from 0 to 100 cycles/second. The speed of wave propagation is 4259 f.p.s., the density of the fluid 1.663 lb.sec2/ft4 and the kinematic viscosity, 0.000197 ft2/sec. Figure 2 demonstrates that the pressure amplification in the third harmonic is smaller than in the fundamental harmonic. Recently several investigators have attempted solutions for transient flow which include the effect of the varying velocity distribution over the cross-section. (354,811) Their work is limited to laminar flow and neglects all nonlinear effects. IThey succeeded in calculating the response of uniform fluid lines for sinusoidal, st ep and impulse inputs, building up an arbitrary input from these elements. The treatment of the nonuniform pipe is even more limited. Similar techniques as for the uniform pipe are used to derive transfer functions for pressure and flow of nonviscous and viscous fluids. ) They are valid only when the duration of the pulse is very small compared to the time required for the pulse to traverse the length of the pipe. After consideration of the limitations'of the present methods, this thesis was initiated with the objective of incorporating the influence of viscous dispersion effects into the one dimensional model of transient pipe flow.

-6 Following a review of noteworthy literature in Chapter II, the fundamental differential equations of waterhammer in uniform and nonuniform pipes and in tubes with very distensible walls are presented in Chapter IIII The equations for the friction head loss ir unsteady laminar flow are derived in the next chapter. If the flow varies sinusoidally the dissipation is a function of the instantaneous mean velocity and acceleration: IWEC 1 6+jV hf - 1 V + + aV f g T g d in which Tj and i are functions of the dimensionless frequency parameter a = R VI/v. If the flow varies arbitrarily with time the instantaneous value of the energy dissipation can be given as a function of the instantaneous velocity and the weighted past velocity changes: hf 8v V + 4v ft V,(u) W(t-u) du gR2 gR2 0 dt in which the weighting function W depends on the dimensionless time parameter T t. R2 Chapter V and VI present the application of the friction terms to waterhammer. For periodic flow cases use of the well established impedance method is proposed and modified by inclusion of unsteady friction which makes the inertance and resistance terms frequency dependent. Arbitrary flow cases are treated by use of the method of characteristics which again is altered in respect to the frictionr Both theories are verified experimentallyo A separate chapter describes the application of the method of characteristics including frequency dependent friction to the flow of blood in arteries.

It should be noted that frequency dependent friction is only one of several causes for dispersion effects in transient pipe flow. Inelastic behavior of the fluid and pipe walls may be of importance as well as a dependency of the wave speed on the frequency of the disturbance. However, for viscous fluids in small diameter pipes there is experimental evidence that these influences are not too importanto

II. LITERATURE The phenomenon of pressure waves in conduits has been studied very extensively during the last sixty years. For this literature review only some specific aspects are of particular interest. One aspect is not even directly concerned with waterhammer but deals with the friction loss due to the unsteady motion of an incompressible fluid, laminar and turbulent. Another one considers the response of uniform pipe to unsteady input taking into account the influence of the varying velocity profile. The existing literature about this topic is limited to laminar flow. Studies which try to extend the methods to nonuniform pipes are also noteworthy since the procedure offered in this thesis is universal in that respect. Finally some literature should be mentioned about the unsolved problem of the stability of unsteady laminar pipe flow. Wall Shear Stress in Unsteady Laminar Flow: E. G. Richardson(17) seems to have been the first who studied the velocity profile of oscillating pipe flow. In an experiment on sound waves in an Helmholtz oscillator he found that surprisingly the maximum velocity does not occur at the center of the pipe but close to the wall, and the higher the frequency of the periodic motion, the closer this maximum velocity is to the wall. T. Sexl(20) verified this phenomenon theoretically. Assuming that the pressure gradient along the pipe is varying as a sine wave, he found an exact solution of the Navier Stokes equations, which results in an expression for the velocity

-9 distribution throughout the cross section. This work was extended by S. Uchida(7) who superposed a periodic oscillation on an average flow and represented it by a sum of sine waves. He examined the velocity profile in detail and calculated the sectional mean velocity and the shear stress at the wall, pointing out the phase difference between them as well as the higher energy dissipation compared to steady flow. While for application of all these equations the pressure gradient had to be known Lambossy(14) eliminated it and derived the wall shear stress as a function of the sinusoidally varying mean velocity and acceleration. (25) For the nonperiodic case, F. Szymanski gave another exact solution of the Navier Stokes equations. With the pressure gradient varying as a Heaviside unit step function, he evaluated the velocity distribution. Due to the linearity of the equation any arbitrary motion can then be represented as a convolution integral, if the pressure gradient is known. This equation is not suitable however for inclusion in a one-dimensional model of transient flow and it is the main objective of Chapter IV of this thesis to derive an expression which allows the calculation of the wall shear stress if only the variation of the mean sectional velocity is known. Wall Shear Stress in Unsteady Turbulent Flow: This paper develops a method of considering the effects of unsteady friction on the propagation of waves through pipes with laminar flow. Since the author claims that an equivalent procedure is applicable to turbulent flow, a brief glance should be given to the little literature available on the subject of frequency dependent friction in unsteady turbulent flow.

-10 F. Schultz-Grunow(19) generated pulsations of flow in a smooth pipe of 11.5 feet length and 1.95 inch diameter. Dependent upon the frequency the instantaneous velocity varied between ^/3.8 f.p.s. and iv7.5 f.p.s. for a time period of the oscillation of T = 2.5 seconds;,v0.16 f.p.s. and ^v6.6 f.p.s. for T = 30 seconds; and /0.16 f.p.s. and v7.1 f.p.s. for T = 30 seconds. Schultz-Grunow noticed a strong similarity between velocity profiles of decelerated flow in a uniform pipe and steady flow in a divergent pipe: both flows are in the directtion of an adverse pressure gradient. A similar resemblance was evident between the velocity profiles of accelerated flow and flow in a convergent pipe. In experiments with a reversal of the flow direction, separation effects occurred and changed the profiles drastically. Although he concluded that the time average value of flow resistance is close to the steady state value, it should be noted that for transient flow the instantaneous value is of great importance since it causes the dispersion of waves. Moreover the frequencies in his experiments were relatively low. J. W. Daily et al. (6) used a one inch diameter pipe 8.25 feet in length. They established steady state flow and then subjected it to a slowly changing acceleration or deceleration; however, in the case of deceleration a rapid velocity change occurred in the start becoming more continuous afterwards. The range of velocities was between 9 and 18 f.p.s., and the range of accelerations or decelerations between 7 and 80 ft/s. When the authors determined the instantaneous DarcyWeisbach friction factor f = 2gDhf/V2 for both steady and unsteady flow they came to the following conclusions: if the fluid is accelerated

-11 the friction is slightly higher than under equivalent steady state conditions, but if the fluid is decelerated the friction is appreciably less than in steady state flow. Referring to the initial phase of -rapid velocity change, which seems to have influenced the friction value of the succeeding flow, the authors concluded that the particular state from which the flow is initiated affects the subsequent flow history. However, it seems reasonable to assume that for small accelerations it should be possible to express the wall shear stress by the instantaneous flow conditions. Following this idea, Mo R. Carstens and J. E. Roller(5) considered the velocity distribution to be the same in steady state flow and in unsteady flow. Then, using a power law of velocity distribution, they developed an expression for the ratio of the unsteady and steady friction factor which depends slightly on the acceleration of the fluid. The authors compared the theory with some of their own experimental results and with those of Daily but the agreement was not very convincing. As one of the possible reasons they too pointed out that unsteady friction may be a function of the history of the motion. If one summarizes these results, it becomes obvious that, up to the present time, no reliable prediction for the energy loss in unsteady turbulent flow is possible. Waterhammer in Uniform Pipes Including Effects of Frequency Dependent Friction: The application of unsteady friction effects to waterhammer problems is relatively new. The response of a pipe with laminar flow to an unsteady input of flow or pressure has been calculated for a

-12 sinusoidal, a step and an impulse input. The basic equations are in all cases the Navier Stokes equation av 1 lap 1 6av 5 + H- - v rA (r and the continuity equation 1 ap v+ a K;4 + -- in which all nonlinear and higher order viscous terms are neglected. In fact the equation of motion is identical with the basic equation for an incompressible fluid in a stiff cylindrical pipe, and therefore the velocity distribution and shear stress, which do not show up in the solution explicitly, are the ones of an incompressible fluid. In the continuity equation the bulk modulus K includes the compressibility of the water and the elasticity of the pipe walls. Based on these equations, A. F. D'Souza and R. Oldenburger(8) derived transfer functions for pressure or velocity between two points of a pipe for a sinusoidal input. Brown(3'4) calculated pressure and discharge at a pipe point due to a step change of pressure or discharge at another point. The analytical expressions of Brown are very complex. The response curve consists of three parts: the head of the curve must be calculated using a high frequency approximation for the propagation constant and characteristic impedance; the tail can be calculated using a low frequency approximation. In the remaining gap none of the approximations is accurate enough and no analytical expression available.

-13 A very complicated procedure is used to convert a more easily calcuable frequency response into a step response. In view of these mathematical difficulties, Brown's results are applicable in a practical sense only if one uses his carefully computed universal step response plots for semi-infinite lines. There remains the work of expressing transient responses by superposing step responses and of taking into account all reflections, an approach which can be very tedious in complicated piping systems. Waterhammer in Nonuniform Pipes Including Effects of Frequency Dependent Friction: Pressure pulses traveling through a tapered pipe are subject to a distortion which depends upon the geometry of the pipe and the viscous shear forces. The prediction of these effects is desirable for the design of pulse transformers in fluid control systems in which the taper is chosen in such a way as to change input amplitude to a desired output amplitude. The only theory of tapered fluid lines which takes into account viscous shear effects, is developed by W. T. Rouleau and F. J. Young. (18) The methods adopted are similar to these used for uniform fluid lines, discussed above, but are much more limited. The applied series converge only for high frequency components and therefore the procedure is only valid for short pulses, in fact only for pulses with a small duration compared to the wave travel time 2e/a. The underlying velocity distribution is assumed to be that of a cylindrical pipe. This assumption is justified only for very mild taper. Blasius(2) demonstrated that for steady laminar flow in a divergent tube a back flow at the wall starts when dR/dx > 12/Re

where Re denotes the Reynolds number. Back flow, however, is the starting phase of separation which causes a rapid increase of friction. The effect will be even higher if the fluid is decelerated since the adverse pressure gradient is steeper. F. J. Tarantine and W. T. Rouleau(26) derived a solution which is applicable to pulses of any duration but for a frictionless system only. Since a step-line-impedance technique is used the pipe can have any continuous taper function and arbitrary transient input must again be approximated by a sum of small step inputs. The method developed in this thesis will be applicable for any slightly tapered fluid line and the input can have any duration. Stability of Unsteady Laminar Flow in Pipes. The laminar flow theory for fluid transients in pipes is not valid if the flow is turbulent. The discrepancies are so large that it would be very important to have some criterion for the stability of laminar flow as a guideline for the upper limit of applicability of the theory. The research on this topic has been started but the reported results are not sufficient. D. A. Gilbrech and G. D. Combs (10) defined a Reynolds number for pulsating flow, relating it to the time average of the crosssectional mean velocity. The critical value of the Reynolds number was interpreted as the value at which the growth rate of turbulent plugs was zero, and then presented as a function of the ratio of the amplitude of the periodic components the time average of the velocity and as function of the dimensionless frequency R [O/v.

-15 The authors concluded that the critical value of the Reynolds number is greater than the one for steady flow and that it increases with increasing amplitudes of the pulsations, reaches a maximum and decreases again but is still greater than the steady state value. It seems that the velocity ratio at which the maximum Reynolds number occurred decreases with increasing values of dimensionless frequency. While these results are valuable for understanding certain aspects of the problem the range of the experiments is too limited to permit final conclusions. One serious problem in the application of the laminar flow theory is the separation which has been observedduring the reversal of the flow direction. Although the flow may still be laminar, the assumption of parallel streamlines is disturbed, which leads to higher energy loss than presented by the theory.

III. FUNDAMENTAL EQUATIONS OF WATERHAMMER The waterhammer phenomenon in pipes is commonly described by a one-dimensional model with pressure head and mean sectional velocity as dependent variables and with independent variables time and distance along the pipe axis. Applying momentum and continuity principles to a segment of dx of a pipe, several authors have derived the basic partial differential equations, which are the equation of motion +V + E + ghf =O gS Y ~+ a +g- T=~ (1) and the equation of continuity av +g vaH g aH 0 ax a2 Vx a2 at (2) in which a denotes the wave speed and hf friction which is related to the wall shear with 7 the specific weight of the fluid. taken to be steady state head loss the head stress by Head loss loss due to fluid hf = 4. T/7 D, is generally 1 V2 hf = f 1D 2g for turbulent flow and hf =32_ v 2 gDc for laminar flow.

-17 Concentrating now on laminar flow, it is of interest to see how the same: waterhammer equations starting with the Navier Stokes Equations can be derived since they will be used to develop a more general term for the fluid friction. The following derivation demonstrates that the problem of taking into account the varying velocity profile can be reduced to choosing an expression for the wall shear stress which is based on the true velocity distribution. If one uses cylindrical coordinates x, r, and assumes the flow to be parallel and axisymmetric in a horizontal circular pipe of constant cross section, the equation of motion is av 6v 1 [p 4 62v 1, _ +v - +- - V (r-) = (3) at ax p ax 3 x r r or and the continuity equation Op av 6p 5T+ p + v E = 0( and for the equation of state ap = ap (5) P K in which p = p(x,t) the pressure v = v(x,rt) the axial velocity v = the kinematic viscosity p =the fluid density K = the bulk modulus of the fluid.

-18 Neglecting the radial velocities implies that the pressure and therefore the density are constant across the cross section. The pipe walls are assumed to be rigid although the bulk modulus may be modified to take into account the elasticity of the wall. 4 62v The viscous term - 2 is small when compared with (r -) and may be neglected. r or or With these assumptions, the equations reduce to av av 1 lp 1 v a + v a + v 1 a (r ) = 0 (6) at ax p ax r or or K ax - (-~ + v P) + = 0 (7) K at ax ax If one integrates over the cross section, divides by the cross sectional area and defines the mean velocity R V -2 f 2T rv dr R R 0 one obtains aV 2 R r v 1 d2 v R a- + 22 f r v"I r + 1r+ P - 2 - (r bv) dr = 0 +R 0x P ax R2 0 (r s r K + V t x )+ K at ax ax Use of the Newton relation for the wall shear stress 6v TO = - pV - r(rR)

-19 simplifies further the integral R 6 -- R - v J a (r av)dr = r av R = R o0 r Tr r r (r=R) R to pv For the nonlinear term v -- an approximation is used ax 2 R 6v2 R | rv r dr - V f r dr Rf o o R2 0 ^ = vV ax then Equation (6) becomes aV av 1 ap 2To V + - + -- =0 at ox p ax Rp These equations are rewritten using the piezometric head H = p/y and the relation between wave speed and bulk modulus a2 K/p. 6v av 6H a-+ V a g a g hf = (1) H + V H a,2 _V With the exception of the approximation in the convective term V -, ax which is in itself small when compared to the remaining terms, the one-dimensional waterhammer Equations (1) and (2) represent every aspect of the Equations (6) and (7), if the friction is calculated from a velocity distribution which satisfies Equation (6).

-20 In the traditional transient flow treatment of laminar flow, To is calculated from the steady state flow velocity distribution which is a solution of av 1 ap 1 a av t + - v - (r -) = 0 at p ax r or or in which the pressure gradient ap = constant and therefore av = 0 ax at The velocity distribution is the well known Hagen-Poiseuille flow and the shear stress can be expressed as a function of the mean velocity To = 4pvV/R. To advance this approach in this thesis, the shear stress will be calculated from a velocity distribution which is a solution of Equation (8) for a time dependent pressure gradient. ax A special note shouldbe made concerning the basic equations for transient flow in nonuniform pipes. As derived in Reference 22 the differential equations are the momentum equation 6H 6V 6V gH V + V + + g hf = 0 g|H~v|Z+^ +gh, = 0 (1) ax ax at which is identical with the momentum equation for flow in uniform pipes and the continuity equation av+gva H+g-^. = o (9) ax +2 Vx * 2 a2t D = Compared with the continuity equation for a uniform pipe, Equation (2), the only difference is in the term, in which 2 designates the D

-21 aD taper of the pipe 3 = - x. The wave speed a is now a function of x due to the varying cross section of the pipe. As mentioned earlier, the underlying assumption of the onedimensional model of transient flow is that the pipe modulus of elasticity is small when compared to the bulk modulus of the fluid. A way to take into account the elasticity of the pipe walls is to modify the bulk modulus of the fluid, in other words its wave speed. The basic mathematical model is still that of a pipe with rigid walls. A further step forward is the application of the same concept to flow in tubes with highly flexible walls which was done in the (21,23) treatment of blood flow by Streeter et al. 1 In this case the compressibility of the fluid is almost negligible when compared with the elasticity of the tube walls. Nevertheless the one-dimensional model of the rigid tube is applied, which means Equations (1) and (9) are valid, and the tube is filled with a hypothetical fluid the bulk modulus or wave speed of which is a function of x and of the instantaneous pressure. This procedure has the advantage that a = a[x,H(t)] can be any nonlinear function, taking into account nonlinear wall properties for example. One specific relation of the (29) variation of the wave speed with the pressure was defined by Wylie. Because of the complexity of the equations, the wave speed is given as a function of the diameter, and, in a second equation, the diameter is related to the pressure head. If the modulus of elasticity of the pipe walls is constant, then a 1 - 2 n D (10) D p Po DDo 2o0 0 a

-22 H__ = 1 + H O ) (11) H ___i-.-2n + H( 2Eoeo D Do 2Eo0e Do, Eo, ao and eo denote diameter, modulus of elasticity, wave speed and wall thickness under an initial head of Ho. p and y are density and specific weight of the fluid. These modifications of the wave speed affect only the equation of continuity. This should be justification enough to neglect the time variation of the diameter in the calculation of the friction loss, which shows up only in the momentum equation, and is generally smaller than the pressure gradient and the inertia term. It is necessary for the derivation of an unsteady friction term, to hold the diameter at a cross section constant.

IV. WALL SHEAR STRESS IN UNSTEADY LAMINAR FLOW A. Fundamental Equations A constant pressure gradient along the axis of a circular pipe filled with an incompressible fluid causes the well known HagenPoiseuille flow which implies a parabolic velocity distribution. Under these conditions, the wall shear stress is customarily expressed as a function of the mean velocity. If the pressure gradient is time dependent, the velocity distribution is no longer parabolic but can be of very complex nature. Moreover the wall shear stress is not in phase with the mean velocity. This phenomenon can be understood by drawing a physical picture of the movement in the pipe. In laminar and in turbulent flow the viscous effect is concentrated in a thin layer close to the pipe walls, called the boundary layer. If the pressure gradient which is the driving force on a short section of the pipe is varying with the time, it affects the boundary layer and the fluid in the center of the pipe differently. Since, in the boundary layer, friction is prevailing and the inertial forces are small, the velocity near the pipe wall is in phase with the pressure gradient. However, in the inner region of the pipe cross section, the inertial forces are dominating, and therefore the pressure gradient is in phase with the acceleration of the fluid. In other words, a sudden change of the pressure gradient will first affect the boundary layer, and therefore change the shear stress at the wall, before the mean velocity is changed. The boundary layer may even change direction of flow before the mean velocity does, a possible cause of the separation observed during the reversal of -23

-24 flow direction. For several specific cases such as pressure varying sinusoidally, as a Heaviside unit step function or the Dirac delta function the velocity distribution has been calculated and discussed in detail. (927) But only for the sinusoidal variation has an equation been derived which relates the wall shear stress to the instantaneous values of mean velocity and acceleration. Some of these results are repeated in this chapter, and, in addition, a new equation is derived for the general case of the mean velocity varying arbitrarily. This equation relates the wall shear stress to the instantaneous mean velocity and past velocity changes. This approach is necessary in applying the friction term to a one-dimensional model of transient flow. The Navier Stokes equation for parallel axisymmetrical flow is restated 62v 1 6 6v 1 v 1 p/ax — ^ + _ -— >" - - > _(8) bra r Or v at v p and for simplification F(t) p Here the term v, which is px x small when compared with the remaining terms, is neglected. This means physically that the influence of the compressibility on the velocity distribution is ignored. Application of the Laplace transform(7) yields the subsidiary equation 2^ d2v 1 dv s l1. v = - F (12) dr r dr v v where v is the transform of v and F the transform of F, e.g.,

-25 C [v(r,t)] =v(r,s) = /v(r,t) e dt 0 and C [F(t)] F(s) - F(t) est dt 0 F(t) is assumed to be continuous or piecewise continuous and s is a complex number, the real part of which is positive and so large that the integral / F(t) e-st dt exists. 0 The partial differential Equation (8) is now reduced to the ordinary differential Equation (12), the solution of which is v(r,s) = C1 Jo(i i r) + C2 No(i r) F Jo and No denote the Bessel functions of first and second kind and of zero order respectively. C1 and C2 are constants which must be determined for the specific boundary conditions. These are v = 0 for r = R and v = finite for r = 0. Introducing them into the differential equation leads to C2 = 0 F C = _ Jo(i \/R) and therefore the transform of the velocity is v(r,s i r) i ) 1 (13) J(i' R)

-26 Furthermore, the transform of the mean velocity is ^. 2 R V(s) - f r dr R2 0 F 2 J (i _ R) -i 1R Jo(i R) F 2 = F f 2 l _ 1 (14) s Lwhc(i X R) J in which the function 7 1(Z) = z ~J() is defined as a modified quotient of Bessel functions of first order.l6) Then the transform of the shear stress at the wall is av To(S) = - [i r= 0 r(r=R) =- p RF Jl(i V R) i - R Jo(i - R) p R F -. pRF (15) 1(i X R) and the transform of the acceleration is _v (s) = s V(s) -= ir 2 - 1] (16) 9i(i R) J Some remarks are necessary about the function 1. An examination of. V5(z) in the complex plane leads to some significant results: 7 l(z) has no zeros and poles which are not real. On the real

axis there is an infinite number of zeros and poles which are interlaced with each other. On the imaginary axis, 71l(z) is real with lim 7 l(iy) = + o. 71(z) is an even function, e.g. l1(-z) = l(&z) y=+oo and the value at the origin is l,(z=O) = 2. The function 7l(aJ -i) with a real will be of particular interest and is discussed in the Appendix. There are some cases for which the inverse transforms of the Equations (13) to (16) can be found. This study is mainly concerned with the sinusoidal solution and a solution for the arbitrary transient case. B. Solution for Periodic Flow If for example the pressure gradient varies sinusoidally, e.g., 1 p iwt F(t) = -P e p ox then (,P - r P ic(t Jo -i) v(r,t) = -i- e 1 - (17) with a = R 4 which is a dimensionless number relating the three parameters radius, viscosity, and frequency. It is understood that in all complex equations only the real part gives the physical value. Furthermore, V(t) = -i P -it [1- 1 (18) To(t) = pRP eiWt(a -i) T (t) = pRP eicut 1 (19)

-28 Above equations are stated without derivation, which can be found in Reference 9. For details on the function 71 (aJ-i) see the Appendix and Chapter IV, Section A. To obtain a relation between the wall shear stress and the mean velocity the pressure gradient is eliminated by dividing Equation (19) by Equation (18). This leads to i pR. T (t) i P V(t) (20) and for the head loss due to friction h(t) = 2 = 2i V(t) (21) V(t) stands for the complex mean velocity V(t) = Re V(t) + i Im V(t) Since V(t) is varying sinusoidally, Im V(t) = - Re EV (t) 3 o-t Furthermore let T and t be two dimensionless numbers depending only on a 1 2i = Re I ( i(a d-i)-2 2i = Im 2 ~ (a $-i)-2 then hf = g l+ i (V- )t and the real part of hf is given by

-29 and hf = 1 v + t v h - lV+Q f g g t pRc 1 pR av T iV + 2 t o 2 a 2 t (22) (23) Numerical values of TM(a) and (a) are given in Table I, while Figure 3 and t (a) for a from 0.0 to 26.0. for a varying from 0.1 to 10.0 shows a graphical plot of ri(a) To examine the limiting value of hf for a -e 0 Equation (21) is used and. Sl(aJ-i) expressed by the first terms of a power series near the origin V1 (aOe -i) = 2 - i (c T-i)2 - (aJ -i)4 - ~4 96 h = C v — 2i (V - i 1 ) R2g 1 2 +1 4 at (^ 2) (V-i|R2 V) = V (8 + ) (V - i ) R2g 3 I2v at in which higher order terms of a are neglected. The real part of hf is h v + 1 aV (24) f Rg 5gt The first term of hf is identical with the expression for steady state shear stress, while the second represents a third of the inertance of the fluid in the pipe.

-30 TABLE I 1 AND E FOR a = 0.1, 0.2,... 10.0. o. *.30.40.50.60.70.80 9 C 1~00 1 * 10 1.20 1030 1.40 1.50 1.60 1.70 1.80 1.90; 2'.) 2. 1C?, 20 2.20 2.30 2.40 2.50 2.60 2.70 2.8C 2.90 3. OC 3.10 3.20 3.30 3.40 3.50 3.60 3.70 3.80 3.90 4.00 4.10 4.20 4.30 4.40 4.50 4.60 4.70 4*80 4.90 5.00 * 2 5,' 0 0.01125.0200C.03125,04499.06124.07997. 7'1 7.,15106 1750 6 21073. 2 4419.28003 31820, 3586('.40 13.44628.49323.54231 *59327.64607,70060 ~75674.81438.87338.93361,99494 1.05723 1.12035 1.18418 1.24859 1.31346 1.37870 1.44420 1.50989 1.57569 1.64154 1.70740 1.77324 1.83902 1.90474 1.97039 2.03597 2.10148 2.16695 2.23238 2.29779 2.36321 o 3 3 4.33333.33333.33333 33332.33331.33328. 33 24.3 31 8.33310.33300.332 35.33267.33245.33217.33183.33142 33094.33 037. 3971.32895.32808.32709.32599 ~32475.32337.32186.32020.31840.31644.31433.31208.30968.30714.30445.30164.29870.29565.29249.28923.28590.28248.27901.27549.27193.26835.26475.26114.25755.25396: * 1 0 5.70 f. o 0 6 * 0 5.20 6 O10 6740 6 70 7. 0 5.10 7. 20 7 C'0 7 50 7.60 7 70 7. 0 7.90 8.00 8.20 8.40 8.50 8.60 8.70 8.80 8 90 9.00 9.10 9.20 9.30 9.40 9.50 9.60 9.70 9.80 9.90 10.00 2.42 (i J4 2.49412 2,55965 2 62526 2.69096 2.75676 2 82267 2, 8 8 8 7 2. 54-F 3 0 21 1' 3 087 55 3.15410 3 * 2 2 0 7 8 3.22078 3 28758 3.35452 3.42157 3.488 75 3.5560' 3,6234. 3.69093 3.75852 3.82621 3 896 I1 3 0 9 71 4.02976 4.0977g 4.1 6:'6. 4 I 2 3 4, C +.3 02' 4.37047 4 50716 4.57559 4*64406 4*71258 4.78114 4.84975 4*91840 4.98709 5, 05582 5.12459 5.19340 5.26224 5.33112 5.40004 5 46898 5.53796 5.60698 5.67602 5.74509. 2 ) 4, 24687 239 2.23652 ~23317.229 6. *?O- 7 r-.. 7,'. 1 1'3.21733.21431.2113 9.20853.20573 * 20 307.20032 19771.10516.19267 19023 18785,18553.18325 *18103.17886 2 1 7 8576.17673.17465.17262.17063.16360.1657~.16491.16309.16130.15955.15784.15616.15451.15290.15131.14976.14824 e.14675.14529.14385.14245.14107.13971.13838

-31 77,__ X__ __________ <35( I0.30 )orA —— _- --.25 ~ ___ __^__ ___ __ ___,o _ _ -20 ~0 2 4 6 8 10 12 14 16 8 20 22 24 261 L10 0 2 4 6 8 10 12 14 16 18 20 22 24 26 a Figure 3. Tr and 5 as Functions of c. N Ici Ne4 9 8 7 6 5 4 3 2 1 __ _____ _____ _______ _ _ _ ______ _ _ _ _ 1 _ _ tL~i __ _ _ 450 40~ 35. 300 25~ 20~ 150 10 5~ 0~ 7 c 2 4 6 8 10 12 14 16 18 20 22 24 26 Figure 4. Ratio of Unsteady and Steady Friction as Function of ca s ( i,I + i t(a)). hf-S ~

To compare the frequency dependent friction with the steady state friction, one can calculate the ratio of the amplitudes and the phase shift between them as a function of the dimensionless parameter a = R No/v. For frequency dependent friction hf = + (22) u g rT g E t and for steady state friction hf - V 8 V s g gg The velocity is varying sinusoidally V = |IV ei1t hf = ( - 1 + - icn) IV eit u g~ g 8( |1 i t h = - 1VI e fs ga2 hfu C hf = 8 n s 2 = i tan-1 (.T) (25) 8 e+ The modulus -8 +; represents the ratio of the amplitudes and the phase angle tan-l(. rB) represents the phase shift between the head loss hf and hf Since hf is in phase with the mean velocity, the phase angle is identical with the phase shift between hfp and the mean velocity V. For numerical values see Figure 4. u

C. Solution for an Arbitrary Flow If F(t) = E is a nonperiodic time dependent function, the inverse transforms of Equations (13) to (15) can be given by the convolution integrals r 2R2 t 6F f ~ R 2 u v(r,t) = (t-u) 2 Z ~ -- e R v O at j=l 3J 1(j) (1 - 4)] du (26) - v J u V^t) = R J (t-U) [ j&-1 X4 e 8 1 - du (27) R2 2 R2 t F e 00 1 v (t) =p- 0 ~ (t- — j du (27) v o atT j=l 2 V kj R2 u 6F e 1 T0(t) = 2RP at (t-u) e du (28) The constants.j are the zeros of the Bessel function of first kind zero order, e.g. Jo(ij) = 0. They are tabulated in Reference 12. In the periodic case it was possible to eliminate the pressure gradient by combining the equations for the mean velocity and the shear stress. This procedure cannot be applied in the more general transient case. It is necessary to eliminate the pressure in the transformed equation and attempt an inverse transform of the resulting equation. Dividing Equation (15) by Equation (14) yields T (S) = — V(s) (29) "J(\1~

It is a necessary property of all Laplace transforms that they converge to zero when s tends through real values to + o. An examination of Equation (29) leads to the result that To as a function of v does not satisfy this necessity. Therefore To is expressed as a function A~ 0 of. Dividing Equation (15) by Equation (16) yields ot 0(S) 57= -pR V (s) (30) 71 (i Y R)-2 at Let p(s) = pR (31) J1 (i V R)-2 To find the inversion 0(t) of (i(s) the complex inversion formula can be used a+iD t 0(t) = lim ir 0(s) e ds (32) W-400 2i C -iw and evaluated by use of the residue theorem for analytic functions. 0(s) is analytic for all complex s except for singularities sj and has a zero at infinity. Thereforeo 1 st O(t) = _- ~(s) e ds 2~ri C = Res [0(s) e, s] (33) j=l J C is a path enclosing all the sj. The poles of 0(s) are the roots of 7l(i \ R) = 2.

-35A detailed consideration shows that the first root is s = 0 and an infinite number of roots are on the negative real axis at j2 s = - j v. The constants Tj had to be calculated, sj - R2 i Ill, 2, = 5.1556, 8.4172, 11.6198.... The interval between them is approaching t. The first 100 values are given in Table II. The function 0(t) is given by the sum of all residues of 0(t) e. They are s.t Res [,(s) et sj] e pR as (i Differentiation of 1 yields al g, R2 6as () 2v which is constant for all residues except the one at the origin. Therefore 2 and for the residue at the origin Res [0(s) est, ] -= R Finally, using Equation (33), 4vp 2vp 2 _j2 T (T) = +-.-.R Z. e (34a) (T) Rj=l

TABLE II I 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 tiE FIRST 100 ROOTS OF THE EUATION 5.13562 8.41724 11.61984 14.79595 17.95982 21.11700 24.27011 27.42057 30.56920 33.71652 36.86?R6 40C.00845 43.15345 46.29800 49*44216 52.58602 55.72963 58.87302 62.01622 65.15927 68.30219 71.44499 74.55769 77.73300 80.87283 84.01529 87.15768 90.3 003 93.44232 96.58456 99 *72677 102 86893 106 01107 109 15317 112.29524 115.43729 118.57931 121 72131 124X86329 128,00525 131 14720 134*28912 137*43104 140*57293 143.71482 146*85669 149.99855 153 14040 156.28224 159 42407 Jo(x) 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 162.56589 165 70770 168.84950 171 99130 175.13309 178.27487 181 41664 184 55841 187,70017 190.84193 193 98368 197 12543 200*26717 203,40891 206 55064 209.69237 212.83409 215 97581 219 11757 222*25924 225 40096 228.54266 231 68437 234 82607 237 96776 241 10946 244o25115 247*39284 250.53453 253 67622 256.81790 259.95958 263 10126 266 24294 269 38461 272.52628 275.66795 278 80962 281 95129 285.09296 288*23462 291 37629 294.51794 297.65961 300.80127 303.94292 307*08458 310 22623 313 36789 316.50954

-37 v in which T -= t is a dimensionless time equivalent to the dimenR2 sionless frequency a = R j/V in Chapter IV, Section B. The series in Equation (34a) is uniformly and absolutely convergent for T > 0. But since the convergence is relatively slow if T is smaller than approximately 0.02 an alternative solution for 0(T) will now be derived which is valid only for small values of T. The transform of 0(T) was given in Equation (31). Instead of applying the residue theorem, 0(s) is developed into a series for large values of s since in the theory of Laplace transforms the behavior of the image function at s -~oo is related to the behavior of the original function at t -0 0. When y is real and large then l,(iy) 1 51 1 63 1 27 1',J (iy) + I + 12 y 8y2 128 Y3 32 7 This series can be derived from the expansions for the Bessel functions JO and J1. Taking the first six terms will be sufficiently accurate. Division leads to 1 1 3 1 15 1 15 1 135 1 451 ^7(iy)-2 - + 2y+ 28 8 l2 1 Y 2 Y2Yy Now 0sA v1/2 1 3 v 1 15 V3/21 15 2 1 A(s) = pR ( - v ++.... R- 2 R2s 8 R3 s3/2 + s R ~ ^'i jr? ^ -^ 5/2 3 135 v5/2 1 45 v 1 128 R5 5/2 32 R6 s R~R s3

-38is tx-/r(x) for x real and x > 0 (6) The inverse transform of l/sX r(x) denotes the gamma function. Thus 0(t) = PR ( /2 t-1/2 C) + 3v 1 2 R2 r(l) 15 v3/2 8 R3 tl/2 3 r(2) 5 v2 t 8 r(2) 135 v5/2 128 R5 t3/2 r(2) 45 v3 352 R6 t2 r(3) ) v If the dimensionless time,T = T- t R is used and J r(1/2) = 1.772454 r () = 1.000000 r(3/2) = 0.886227 r(2) = 1.000000 r(5/2) = 1.329340 r(3) = 2.000000 then 1 = Vp. + (0.282095 T1 R R (T) - 1.250000 + 1.057855 T1/2 + 0.937500 T + 0.396696 T3/2 - 0.351563 T2) (34b) for T < 0.02. The terms have been arranged to show the similarity to Equation (34a).

-39 The convolution of 0(t) with the rate of change of velocity, Equation (30), leads to the desired equations for the wall shear stress and the head loss R ( 2pv t ~a r,0(t) = v V(t) + R (u) W(t-u) du 0R K0 5 a (35) 8v _47 hf(t) = v (t) g R g R T = R2 t av f (u) 0 t W(t-u) du (36) W(T) 00 2 - e- j T j=1 -5.1356 T = e -8. 1722 -11. 61982 + e -14.79602T -17.95982T +e 79oT+e (37a) T > 0.02 1 W(T) = 0.282095 T - 1.250000 + 1.057855 T 1/ + 0.937500 T + 0.396696 T/2 - 0.351563 T2 (37b) T < 0.02 Equations (35) and (36) mean that the instantaneous energy loss is the sum of the steady state value plus a term in which certain weights are given to the past velocity changes. The weights can be calculated from Equations (37a) and (37b) which deliver identical results in an overlapping region. Numerical values of W are given in Figure 5 and in Table III, A-C.

14 W(T) 12 C.04 T Figure 5. The Weighting Function W as Function of T.

TABLE IIIA E WEIGHTING FUNCTION W(T) WITH T VARYING FROM T - 0.0001 TO T = 0.0100 IN STEPS OF O.00OO T (WI Oi II *- )OO'~ 26.97014 10.29295 7.29161 5.84621 4.95630 4.33876 3.87844 3.51851 3 22726 2.98544 2.78058 2.60422 2.45037 2.31467 2.19386 2.08543 1 98745 1.89835 1 * 81691 1 *74210 9 * 4 4 3 6.9 3116 5.637:4 4.8 16 0 4 23 647 3.79969 3.45551 3. 1754 1 2.94182 2.74325 2 57180 2. 4219 2.28940 2 171 2 5 2 065 6 1 96896 1.88149 180144 W(T) 15.05535 8.75423 6.61328 5.44564 4.68501 4.13975 3.72457 3.39504 3.12541 2.89960 2.70699 2.54025 2.39411 2.26471 2.14912 2.04508 1 95082 1 86492 1 78624 12.87627 8.18575 6.33021 5.26962 4.5 36235 4.04811 3.65283 3.33694 3.07714 2.85870 2.67177 2 50951 2.36700 2.24057 2.12745 2.02550 1.93301 1 84864 1.77128 1.70004 11.380n7 7 70502 6.07606 5 10705 4.44718 3.6114 3.58421 3.28106 3.03051 2.81904 2.63752 2.47956 2.34053 2.21696 2.10623 2.00629 1.91552 1.83264 1.75657 1*68646. 0)6i.4 071,:' 91 1.72785 1.71384

-42 TABLE IIIB TrE WEIGHTING FUNCTION W(T) WITH T VARYING FROM T - 0.0105 TO T = 0.1000 IN STEPS OF 0.0005 TW(T).0105 1.62158 1 35744 1.16271 1.01220.0205.0305.89185.79316.71061.64049.58015.52768.48164.44094 1.56133 1 *31398 1.12963 *986C4.87059.77549.69569 *62771 q56908.518(C.47311.43337.39796 * 36625.33771.31191 1 50520 1,*2299 1.09816.96100.85013 *75844 *68124.61530,55831 *50856 *46478 *42596.34134.36030 *33234,30705 1.45274 1.23423 1.06816.93700.83044.74196.66724.60325.54783.49936.45664.41872.38487.35448.32708.30228 1.40360 1.19753 1.03954 *91397.81146.72603.65366.59154.53762.49039.44870.41165.37853.34877.32193.29760.0405.40473.37232.34318.31687.0505.29301.27132.25152.23341.0605.21681.20156 18752.17457.0705.16261.15155.14131.13182.0805.12301.11483.10722 *.10015.0905 *09356.08743.08171.07639;28851.26721.24777.22998 *21366.19866.18484.17210.16033.14944.13935 *13000.12132 11326 * 10576 e09879.09230 *08625 *08062.07536.28409.26318.24409.22660.21056 *19580 *18221.16967 15808.14736.13743.12821'11966.11172.10433.09746.09106,08509.07954.07436.27975.25922.24047.22328.20751.19300.17962.16728.15587.14531.13553.12645.11803.11020.10291.09614.08983.08395.07847.07336.27550.25534.23691.22002 2 0451.19024.17708.16493.15370.14330.13366.12472.11641.10870.10152 *09484.08862.08282.07742.07238

TABLE IIIC THE WEIGHTING FUNCTION W(T) WITH T VARYING FROM T = 0.1010 TO T - 0.2000 IN STEPS OF 0.0010 7TW(T).101.07C46 o.065.06678.06501,06329.C 6162.05929 05841.05687.05537.111. 5391.05249 * 0511 1.04977 *04846.44718.04594 *04474.04356.04242.121. 4131.04022 * 03917.03814.03714.03617 *03522 *03430.03340.03253.151.r3168 *03085 *03004.02926 *02849 o.2775 *027 v *02632.02563 *02496.141 *.2431.02367.02306.02245.02187 *02130.02074 *02020.01967.01916.151.01866.01817 *01770.01724.01679 *01635 *01592 *01551.01511 *01471.161.01433 *01395 301359.01324 0129O_ <.01256,01223.01191 *01160,01130.171.31100.01C72.01044.01017.00'99(.,0964.00939.00915.00891.00868.181.00845 *00823,00802.00781.00760.00741.00721 *00703.00684.00666.191.00649.00632.00616.00600.00584 -.0569.00554 *00540.00526.00512

V. APPLICATION OF FREQUENCY DEPENDENT FRICTION TO STEADY OSCILLATING FLOW A. Method of Solution Steady oscillations of pressure and discharge in fluid conduits are a specific problem of the waterhammer phenomenon. They occur if a periodically changing boundary condition exists long enough to establish flow conditions which are periodic at every point of the system, i.e., repeat themselves identically with every cycle. If the distribution of the fluid is taken into account, one is dealing with a steady state wave phenomenon the solution of which has been well known for a long time. A modern tool of representing the solutions of the wave equations is the method of impedances (2228) which is particularly convenient for handling complicated pipe systems. (8) A. F. D'Souza and R. Oldenburger have treated the dynamic response of fluid lines, taking into account the cross sectional velocity distribution. Their work results in a transfer function for pressure and mean velocity at two pipe locations. They do not derive an explicit expression for the friction term, but include the velocity distribution by dealing with the radial coordinate as a third independent variable, and later obtaining a cross sectional average. The present study is different insofar as it calculates the wall shear stress explicitly and includes it in the existing resistance-capacitanceinertance model which has been represented by the impedance model. The numerical results should be identical to those of D'Souza and Oldenburger. -44

-45 The equations of waterhammer are used in a simplified form since all nonlinear terms are neglected or linearized. Then the momentum equation and the continuity equation become jH + L a~ + R - 0 (38) C L +- 0 (39) at ax The notations R, C, L are equivalent to the ones used in the analysis of electric transmission lines. They describe the cross sectional parameters of a fluid line. C is the capacitance per unit length which represents the storage capability of the pipe L describes the inertia of the fluid R stands for the energy dissipation due to fluid friction. Equations (38) and (39) are modified 6H 1 6Q 4 -^ ^ — T+ o = o (40) ax gA at 7D 0 +H a2 (41) at gA ax and for the head loss due to friction Equation (22) is used

-46 4 hf = -D yD T = 0 Agir Q Ag at Ag at (22) v and. are dimensionless functions of the parameter see Chapter IV, Section B. Equation (40) becomes a R J)-/V,.. + ax - (1 + ) - + gA at Agn Ag r 1 and after comparison with Equation (38) one gets R = -cAgl (43) (44) 1 L = 1-g gA (1 + ) a2 C a gA (45) Use of the steady state friction term hf = 352 QR would lead to gD2A R = 32v gD2A 1 L = 1 gA 2 a C - gA (46) (47) (48) The limits for the inertance and resistance terms are

-47c 0 R 5 532v gD2A 4 1 L - 3 gA D- oo00 R oo L - 1 gA The steady state friction term h = 32 Q has been widely gD2A used to describe the energy dissipation in transient flow of a viscous fluid. However, it should be noted from Equation (24) that for -, 0 the head loss tends to hf = 1 + 2 Q, the second term of 3Ag at gD2A which is identical with the expression for steady state friction, while the first term amounts to one third of the inertance of the pipe. It is obvious that this relation is a more accurate low frequency approximation then the steady state friction term. While here the use of frequency dependent values for resistance R and inertance L is proposed, the capacitance has been held constant. It would be no major problem to include a frequency dependent wave speed and the resulting variable capacitance into the analysis. It seems, however, that at least for rigid pipes this influence is not of major importance. Equations (38) and (39) can be solved, assuming that the dependent variables H and V vary sinusoidally. Arbitrary periodic variations are handled using a Fourier analysis.

The derivations leading to the propagation characteristics of the pipe and transfer functions can be found in Reference 22. The following is a summary of the results. The complex propagation constant y = NCo (- L + iR) (49) and the characteristic impedance zc Y iC)C icuC (50) are functions of the physical properties of a cross section of the pipe and the frequency of the oscillation. The hydraulic impedance H Z = Q (51) includes the geometry of the piping system; it stands for the change of head necessary to cause a unit change of discharge. In the complex notation used here, the modulus of Z gives the ratio of the amplitude of head fluctuation to the amplitude of the discharge fluctuation, whereas the phase designates the phase difference between head and discharge. Given two of the values H, Q and x1, the values at another cross section x2 the following transfer functions: Z at one cross section can be calculated using

-49 H(x2) = H(xl) cosh 7 (x2-xl) - Q(xl)Z sinh 7(x2-x1) (52) H(xl) Q(x2) = - sinh 7 (x2-xl) + Q(xl) cosh y (x2-xl) (53) ZC Z(x2) = H(x2) Z(l) - ZC tanh y (2-xl) (4) Q(x2) Z(x2) 1 - Z tanh 7 (x2-xl) ZC Special boundary conditions and treatment of branches, series pipes etc. are discussed in Reference 22. B. Experimental Verification To demonstrate the significance of frequency dependent friction a simple experiment was carried out, which measured the pressure amplification due to resonance in a plastic tube. The system used is shown in Figure 6. The length of the tube was 72.3 ft., the diameter D = 0.248 in., and the wave speed was measured to be a = 1003 f.p.s. The fluid was water with a viscosity of 1.078 x 10-5 ft2/sec At one end of the tube discharge oscillations were forced by a piston pump, the frequency of which was varied between 0 and 50 radians per second, while the other end of the tube was closed. The pressure was measured at two points, using two pressure transducers. The ratio between downstream and upstream pressure was taken and compared with the results of the mathematical model. The transfer function for the head is H2 1 (55) Hi cosh 72 Figure 7 shows the numerical results obtained from this equation, using the frequency independent terms and the frequency dependent terms for L and R.

-50 DIAMETER OF PIPE.25 INCH WAVESPEED KINEMATIC 1003 fps. VISCOSITY 1.078 10-5 ft2/sec CLOSED END 72.3 FEET a. a. Figure 6. Experimental System for Study of Pressure Amplification Due to Resonance. ~12 ~ I I0 L to_ < I w 7^~~~~~~~~~~~~~ I I w a. w 1 o 0 < < u [RADIANS PER SECOND] Figure 7. Ratio of the Amplitude Between Downstream and Upstream Pressure versus Frequency. ------- Calculated Results Using Steady State Friction. Calculated Results Using Frequency Dependent Friction. 00 Experimental Results.

-51 Using steady state friction, the highest pressure amplification with resonanT- frequency is predicted to be 34.4, application of frequency dependent friction leads to 11.5 and the measured amplification was 7.5. The measured pressure build up, which is still smaller than predicted using frequency dependent friction, may be partially influenced by visco-elastic properties of the plastic tube. Of significance is the slight change of resonance frequency which can be seen in the experimental data and which is predicted using frequency dependent friction, but not using steady state friction. C. Numerical Example In the study of the pressure amplification versus frequency, Figure 7, one may observe that the prediction is fairly independent of the friction if the exciting frequency is not close to resonance conditions. The following numerical example will show that although the amplitudes may be only slightly different a significant change of the wave shape is nevertheless caused by frequency dependent effects. A hypothetical system is assumed with a one inch I.D. pipe, 1000 feet in length, Figure 8. The wave speed is 4000 f.p.s. and the kinematic viscosity of the fluid v = 0.000427 ft2/sec. The boundary conditions are a constant head of 50 feet downstream, and the head varying as a square wave upstream, Figure 9. This discontinuous variation causes waves with high frequency components which are subjected to higher damping than the low frequency components. The input frequency is C = 11.27 rad/sec. which is between the natural frequencies of the pipe for the first harmonic, C = 6.28 rad/sec., and the second harmonic, C = 12.56 rad/sec.

-52 DIAMETER WAVESPEED KINEMATIC HEAD IN OF PIPE I INCH 4000 fps VISCOSITY.000427 f?/sec RESERVOIR 50 FEET RESEVOIR RESERVOIR A L — 500 FEET Figure 8. B C 500 FEET -- Hypothetical Pipe System for Numerical Study. Uz 2 Z Iz a. 4 100L 50 m Q I T=.555 I I T=.555 Figure 9. Assumed Head Variation at Point A of the System in Figure 8.

-53 The square wave was approximated by 28 terms of a Fourier series, and the response of the system calculated using steady state friction and frequency dependent friction. Figures 10 and 11 show the discharge variation at the upstream point A, and pressure variation at the midstream point B. In both cases, the maximum fluctuations, using frequency dependent friction, are smaller than for steady state friction. The difference, however, would be much more significant if the input frequency were chosen closer to resonant conditions. The varied friction causes a phase shift between the results and a rounding effect of the waves. In the example, the frequencies of the 5th, 15th and 25th Fourier component of the square wave are very close to 9th, 27th and 45th harmonic of the pipe. This causes an amplification of the affected components in the calculation using steady state friction. The frequency dependent friction attributes more damping to the higher harmonics, and the scattering effect does not occur.

I A f-, I rLJ z LLJ m pz 0 iL LI 130 120 110 - 100 - 90 80 70 60 50 _ __'. 30 20 1 / - I - e___10 L i__I 0~~~ l^ >oo.1.2.3.4.5.6 TIME IN SECONDS Figure 10. Head Variation at Point B Resulting from a Square Wave Input of Head at Point A. - -.. — Using Steady State Fric-tion..-^ —- Using Frequency Dependent Friction.

-55-.006.oo005.004 --- j - t V-.o _... 1.,. CI I-.002... w O..00...001N SECON.0..2.3.4.5.6 TIME iN SECONDS Figure 11. Discharge Variationr at Point A Resulting from a Square Wave Input of Head at Point A. o —-o — Usirng -t;,eady St-ete Frti cition, -i ->0-o Using Frequentcy Dependent Friction.

VI. APPLICATION OF FREQUENCY DEPENDENT FRICTION TO TRANSIENT FLOW CASES A. Method of Solution The solution of the waterhammer equations given in Chapter V is limited to periodic oscillations. Although similar impedance techniques have been used to calculate the transient response of fluid systems, this thesis emphasizes the use of the method of characteristics(2224) because of the distinct advantages outlined in Chapter VIII. The governing equations for a tapered tube are two quasilinear hyperbolic partial differential equations: see Chapter III. 6H 6v 6v g + v + + g hf 0= (1) V 9+ vaH +g 6aH v0 O(9) 6x a2 V Tx a2 at D They can be transformed into a pair of total differential equations, the validity of which is restricted to certain lines in the x - t plane called the characteristics. This procedure leads to -++ghf-aV- a dt dt g g dH + dV h2 _ + +ghf+a (57) adt dt D -56

-57 and for the characteristic lines - = V +a (58) dt d = V - a (59) dt For most practical cases a first-order finite difference approximation has been found to be sufficiently accurate to approximate the four equations. Vp - VR + g- (Hp- ) + (tp-tR(ghfR- aR VR ) = (60) aR DR Vp - VS - - (Hp-HS) + (tp-tS)(ghfS+ aS VS D) = o (61) aS Xp - xR = (VR + aR)(tp -tR) (62) Xp - xS = (VS - aS)(tp - ts) (63) Having now four algebraic equations, the four unknowns Vp, Hp, xp, tp can be calculated from the known values at the points R and S, see Figure 12. In many cases, it is more convenient to use a grid with specified time intervals. This leads to an interpolation procedure described in Reference 22. The methods of treating different boundary conditions are also given there.

-5R t P S R Characterst Lines on the x -t Pane. Characteristic Lines on the x - t Plane. Figure 12. KI Pi, K I - j AtL 2 -7 -1 —--:~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ I i, i1i1 t i Ii 0 2 n I X IX rI Figure 13. Grid of Characteristics for Specified Time Intervals.

-59 If the velocity of the fluid is small when compared to the wave speed, the velocities VR and VS in the equations for the characteristic lines can be neglected. This approximation corresponds av 6H to a cancelling of the nonlinear convective terms V -- and V - ax ax in the differential Equations (1) and (9). No restriction have yet been made to the friction loss hf, and the wave speed a, which may be any linear or nonlinear functions of the instantaneous values V, H, t and x. But they certainly do affect the accuracy of the first order finite difference approximation, which is sufficient only if the friction term g hf, and the taper term a V 2, are small when compared to the other terms in the total D differential Equations (56) and (57). Henceforth the term"simplified equations" is used if in the basic Equations (1) and (9), the convective terms V v and V H ax ax are neglected, and the wave speed is assumed to be constant. Then specified time intervals can be used without interpolation, although for tapered tubes the Ax segments must be adjusted properly. It is customary to apply the steady state friction terms 1 V2 32v h = f - - for turbulent flow, and hf - V for laminar flow, D 2g D2g in the method of characteristics. As pointed out in the introduction, these terms cannot describe the energy dissipation with sufficient accuracy in fluid transients with higher frequency disturbances. In Chapter IV, Section C, a more accurate equation for the wall shear stress in unsteady laminar flow has been derived, which depends on the instantaneous velocity, and the past velocity changes. There is no way of expressing the unsteady fluid friction only by instantaneous

-60 values of the velocity or the acceleration except in sinusoidal flow. Since the method of characteristics is a time step method, in which the next step is calculated using the results of the previous time step the history of the velocity is known, and the equation for unsteady friction can be applied. The head loss at an instant t is given by Equation (36): hf(t) D2 V(t)+gD16- t hf(t) = 2v ) + 16v Wt V (u) W(t-u) du gD gD 0 at = h(t) + h'(t) (36) h is the steady state friction loss, h' the component due to the unsteady motion of the fluid. Using a grid with specified time intervals, as shown in Figure 13, the head loss at point Pik can i,k be calculated from a first-order approximation of Equation (36). h = h +h' i,ck i,k i,k - 532v hic k 2 Vi,k k-l hk 2 l iv j (Vi,j+l - Vi,j-l) W((k-j)At) i,k g j=l,5,... gD k-l h 16vw Z - (V j - V ) W(jAt) gD2 ik-j+l ik-j-l t) + (V6v [(Vi,k - Vi k-2) W(mt) + (Vi k-2 - Vi k-4)W(3At) gD2,,, +... + (Vi,2 - Vi,) w((k-l)At)] (64)

-61 Similar series can be set up for any grid in which the abscissas x of the grid points are constant. Also higher order approximation can be used to calculate the integral. The weights W are given in Table IIIa-c and Figure 5 as the function of the dimensionless parameter T = v- t. W(T) approaches R2 zero for T - o. If one choses a rmax so that all weights W(T) for T > Tma are negligible, then the number of terms in Equation (64) never surpasses Tmax R2/(2vAt). But since their number can still be high in many practical cases, a modified procedure to calculate hf will now be given. The idea is to employ the friction which has been calculated in the previous time step and calculate only its change within the time 2At. From Equation (36) one obtains hf(t-2At) = h(t-2At) + ht (t-2At) t-2At hi (t-2At): -2f (u) W(t-2At-u) du gD2 0 5t and therefore 16v t 6V h'(t) - h'(t-2At) g= D2 2t (u) W(t-u) du gD2 t-2At t 16v t-2At V + 72 f - (u) W(t-u) du gD O at t-2At 16v f V (u) W(t-2At-u) du gD2 o at

-62 and hf (t) = h(t) + h' (t) 32v V(t) + h' (t-2At) + i6 gD2 i6v ft-2At 6V + 16 /- (u) [W(t-u) gD2 0 t t f t-2At 6vU) W(t-u) du 6t - W(t-2At-u)] du (36a) and after using a first order approximation h. = h., + hi ik i,k ik - 32v i,k gD2 Vi,k hi,k, 16v i,7k-2 +gD2 k - Vik2) w(At) 16v gD2 k-3 j=,13,. = h.. 16v i,k-2 g+ 2 16v gD' z j =1,?3,. (Vi,j+l- Vi,j-_)[W((k-j)At)-W((k-j-2)At)] (Vik- Vi,k-2) w(At) (Vi,k-j-l- Vi,k-j-3)[W((j+2)At)-w(jAt)] (Vi,k - Vi,k-2) w(At) 2 - i,k-4)[W(3At) - w(At)], 16v i,k-2 gD2 + 6v {(v k gD2 ikgD + (vi,k-4 - Vi,k-6)[W(5At) - W(3At)] +.. + (Vi,2 - Vi,0) [W((k-1)At) - W((k-3)At)]} (64a)

-63The necessary number of terms in Equation (64a) is smaller than in Equation (64), because the velocity changes are no longer weighed with W(t) but with the differences W(t) - W(t-2At). Figure 5 demonstrates that they approach zero much more rapidly than W(t). One can choose a suitable Tmax and the maximum number of terms is again Tmax R2/(2vAt). B. Experimental Verification To verify the procedure of including frequency dependent friction in the method of characteristics, numerical results are compared with experimental results by E. L. Holmboe. He ran tests on a one inch I.D. copper tube connected upstream to a 60 gallon tank, which was maintained at constant pressure by compressed air. A quick closing valve was mounted to the downstream end. The tube had a length of 118.4 feet, was coiled in the form of a spiral about three feet in diameter, and embedded in concrete to reduce vibration; the wave speed was measured to be 4345 f.p.s. Pressure transducers were located at the valve and the midpoint of the tube. A sketch of the experimental set up is given in Figure 14. To produce a noticeable viscous shear effect, an oil was used as fluid having a viscosity of 0.427-10 3 ft /sec at 80~F which is almost 50 times that of water. The traces of the pressure fluctuation after a sudden closure of the valve are reproduced in Figure 15. Holmboe indicated a relatively high range of calibration error for the pressure of + 8%, and adjusted the pressure scale so that the sudden pressure step, immediately after the closure, corresponds to the well established value of a vo/g

DIAMETER OF PIPE I INCH WAVE SPEED 4345 fps 2 KINEMATIC VISCOSY.000427 ft/sec PRESSURE TRANSDUCERS QUICK CLOSING VALVE I Figure 14. Set-up for the Experimental Work of Holmboe.(11)

-65 Numerical results are obtained from the method of characteristics using laminar steady state friction and frequency dependent friction. For the mathematical model, the basic Equations (1) and (9) are simplified since V < a, a = constant, = 0. A grid with specified time intervals At = 0.00136 is used, the number of sections is 20 and the integral in Equation (36) is replaced by a first order approximation as done in Equation (64). Figure 15 shows the comparison between calculation and experiment. The application of steady state friction gives only a rough approximation of the measured pressure fluctuation. However, the inclusion of frequency dependent friction not only predicts the higher decay of the fluctuation, but ev very accurately details the disperson of the wave shape.

to I I 12 10 14 ir ^r ^ ~to — " EXPERIMENTAL ------- STEADY STATE FRICTION UNSTEADY FRICTION Figure 15. Head Fluctuation at Valve and Midstream after Sudden Valve Closure. Compnarison with Calculated Results Using Steady State and Frequency Dependent Friction. Ex-eri::: ental Work of Holmboe.(11)

VII. APPLICATION OF FREQUENCY DEPENDENT FRICTION TO THE FLOW OF BLOOD IN ARTERIES A. Basic Assumptions One important application of the method of characteristics including frequency dependent friction is in the calculation of blood flow in the arterial tree. As pointed out in the chapter on basic equations, the one-dimensional approach is approximate for a highly flexible tube like an artery. The flexibility of the walls is only taken care of by making the wave speed pressure dependent, which affects only the capacitance effect of the system. The balance of the forces, pressure gradient, inertia, friction is that of a rigid tube. An alternative to this approach is used by Wormersley,(30) who introduced the radial and axial velocities, and the pressure, as dependent variables of the cylindrical coordinates x and r and the time t. After some assumptions,the most important of which is the linearization of the equations, he formulates a relationship of the instantaneous mean flow to the pressure gradient. The one-dimensional approach is relatively crude compared to Wormerley's theory, but it has some distinct advantages in that it includes nonlinear effects, taper of the artery, branches, etc. The introduction of frequency dependent friction into the one-dimensional model moves it closer to Wormersley's model since the effects of the varying velocity profile on the equilibrium of forces are considered. The basic partial differential equations are Equations (1) and (9), which are transformed into total differential equations and approximated by first-order difference equations as shown in Chapter VI. -67

-68 The significant feature of the application of these equations to blood flow in arteries is that the wave speed depends upon the pressure. The relation used in the following example is the one of Equations (10) and (11). It implies that the modulus of elasticity of the wall is constant, and the arterial wall is incompressible, Poissonvs ratio = 0.5. The artery is tethered and the compressibility of blood is neglected. Two different equations for the friction loss have been derived in Chapter IV. They imply the assumption of treating blood as a Newtonian fluid, i.e. the viscosity is independent of the shear rate; 3 2v 16 t hf V + f2 f (u) W(t-u) du (36) gD2 1D2 0 and h = - V + _ (22) f g T g at Equation (36) is valid for any arbitrary flow case, while Equation (22) is true only for a sinusoidal motion, and,if using a Fourier analysis, for any periodic motion. Although the flow of blood is periodic, due to the rhythmic action of the heart, a direct application of Equation (22) on the method of characteristics is not possible because the Fourier components cannot be determined beforehand. In Chapter VI is a description given of how the friction can be calculated using Equation (36). For this purpose, the grid has to be chosen so that the x-coordinates of the grid points are fixed. This

-69 (21,23) demands an interpolation procedure which is quite time consuming. In the application to blood flow, a floating grid, as described by (29) Wylie, has been found to produce only minor shifts in the x-coordinates of the grid points. If a grid point P deviates a small distance 5x from the line x = constant, on which the integration is carried through, the velocity on the line may be calculated by V(x) = V(x+x) A(x+5x) and its value used for the integration. This v(x) = v(x+Gx) A(X) is true since Q is small compared to -Q A more approximate way of taking into account frequency dependent friction will now be given: for small values of a Equation (22) tends toward Equation (24) 3= 2v 1 +v h = 2v V + 1 a- (24) hf D2g 3g (t The proof is given in Chapter IV. Equation (24) may be written more generally: hf = hf + V (65) The unsteady friction term, hfu, is the sum of a steady friction term, hf, plus the acceleration multiplied with a constant,. The basic differential equations become: 6H 6v 6V g + V a+ (l+~ ) t + g hf = O (66) g (V H H ) +V 2- V - 0 (67) a -a x +t + x D

-70 For convenience let D = + ( VR 2 ( and XS = - Va| (1+) + (VS )2 +VS g Then the difference equations, equivalent to the Equations (60) to (63) in Chapter VI, are (68) (69) (l+~)(Vp-VR) + g (HP-R) aR2 + (g hf - XR D VR) (tP-t) R R D0 (70) (1+) (Vp-Vs) + gX3 __s (p-Hs) aS + (g hfS - XS D VS)(tp-t) 0 XR i= 1 (VR+X ) (tp tR) PX~ = D (vS(tt1 ) (71) (72) (73) These are four equations which can be solved to calculate the four unknowns Vp, H, Xp, tp. Instead of using Equation (24) one can in the same way apply Equation (22) with the constants q and t determined for the frequency of the first harmonic, if the time average of the flow is not too high.

-71 B. Computed Results A comparison of calculated results with in vivo measured blood flow data fails, due to the lack of reliable measurements. But there is no question that a mathematical model which includes the influence of the varying velocity profile, is more realistic than one which assumes the velocity distribution to be parabolic. However, more than in any technical application, approximations which simplify the resulting equations are justified in blood flow, since the assumption of a full developed velocity profile is very approximate in itself. Arteries are seldom long and straight enough, and the influence of entry disturbances and outflow at branches is too dominant to create truly axisymmetrical flow. But this should not effect the basic assumption that only when the history of the velocity at a cross section is known, an evaluation of the wall shear stress is possible. Two numerical studies have been made which are closely related to each other. In the first one, the head loss due to friction is calculated from a known flow curve. This task is somewhat theoretical but it demonstrates the different assumptions better than an actual model of an artery in which friction is only a secondary effect. Next, the flow in a tapered arteryis calculated, assuming the pressure variation at the proximal and the distal end to be known. For this the method of characteristics is applied as described. 1. Comparison of Computed Steady and Unsteady Friction The mean flow in a stiff, long tube with a 0.316 cm diameter is assumed to be periodic with a time period of 0.4 seconds, see Figure 16. The kinematic viscosity is v = 0.05 cm2/sec. Figure 17

3 2.8.7.6 Z z o w (C) w Q. 0 z LL z 4.4. n 0 ox w a.6 ".1 r) I 0.4 TIME IN SECONDS TIME IN SECONDS -I Figure 16. Assumed Flow Variation in a Rigid Tube of 0,316 cm Diameter. Figure 17. Head Loss Due to Friction for the Flow Variation of Figure 16. Curve a: Steady State Laminar Friction. Curve b: - -_ Steady State Turbulent Friction, f " 0.4. Curve c: ------ Unsteady Laminar Friction, Using Equations (36) or (22) Curve d:...... Approximation to Unsteady Laminar Friction Using Equation (24).

-73 shows the head loss due to friction calculated for one unit length of pipe under several different assumptions. Curve a: steady state laminar friction hf = 32 V gD 1 V2 Curve b: steady state turbulent friction h f = 0.5 D 2g Curve c: unsteady friction using Equation (36) and Equation (22) leading to almost identical results. This curve represents the exact solution of the Navier Stokes equation, Equation (8). For application of Equation (22) the flow curve is described by the first 18 terms of a Fourier series, hf calculated for every frequency and the results superposed. The integral in Equation (36) is evaluated, using a first order finite difference method, as described in Chapter VI. The calculation is carried through for two cycles to eliminate the transient solution, Tmax is 0.05. Curve d: approximation to unsteady friction based on Equation (65), 32v 1 +V hf _ 3v V + D2g 3g at Comparing the different curves, some important results can be noted: there is a significant phase. hift between unsteady and steady state friction, and therefore between unsteady friction and mean velocity. Furthermore, for accelerated flow, the unsteady friction, curve c, is greater, and for decelerated flow, smaller, than steady state friction, curve a. The ratio of the amplitudes is about 1.35 for the first positive peak but about 2.00 for the negative peak. Finally, the approximate expression of Equation (24) describes very closely the unsteady friction, at least for the physiological range of data of the example.

The Reynolds number never surpasses 230 during the cycle and therefore the probability is very small that turbulence develops. But even then it is interesting to add a friction curve which is calculated from the assumption that hf varies with the square of the velocity. A very high Darcy-Weisbach friction factor of f = 0. 5 had to be chosen to match the first positive peak, but the negative peak is still much smaller than calculated from a linear friction term. 2. Comparison of Computed Flow Curves The model consists of an elastic tapered tube 9.4 cm in length with a diameter of 0.316 cm and a wall thickness of 0.026 cm at the proximal end. The vessel tapers linearly to 0.228 cm at the distal end. The modulus of elasticity is 1.76-107 dynes/cm2 and the wave speed head relationship of Equations (10) and (11) is used for an initial pressure level of 11.5 cm Hg. The blood viscosity is chosen to be 0.05 cm2/sec. Proximal and distal pressure variations, Figure 18, are used as boundary conditions. Two different computer programs applying the method of characteristics were written. The first one uses Equations (60) to (65) and specified time intervals, and necessarily, an interpolation procedure as described in References 21 and 23. The second program uses Equations (70) through (73) and a free floating grid as described in Reference 29. In a first run, steady state laminar friction was employed in both programs, i.e. ~ = 0 in Program 2. This proved that the interpolation in Program 1 does not cause a significant change in the resulting flow curve, which is given as curve a in Figure 19. Then

-75 x 0 I z 0 w I Figure 18. 13 12.05.15.20 TIME IN SECONDS.25.30.35.40 Head Variation at Proximal and Distal End of a Tapered Artery. 9 3 a) 2 u - z IL.20 25 40 TIME IN SECONDS Figure 19. Flow Variation at Proximal End of Tapered Artery Resulting from the Head Variation of Figure 18. Curve a: Using Steady State Laminar Friction. Curve b: _.-_ Using Steady State Turbulent Friction, f=O. 5. Curve c: ----- Using Unsteady Laminar Friction, Equations (36) or (22). Curve d:..... Using Approximation to Unsteady Laminar Friction, Equation (24).

-76unsteady friction of Equation (56) was used in Program 1, curve c in Figure 19, and the approximate relation of Equation (22) in Program 2, curve d in Figure 19. The most important result of a comparison of the calculated flow curves is in the negative peak, which becomes significantly smaller if frequency dependent friction is employed in the analysis. Former comparisons of calculated flow curves with experimental blood flow data by Streeter et al. (21,23) demonstrated that the main disagreement between theory and experiment occurred in the prediction of the backflow portion of the flow which generally was much smaller in the in vivo data than in the calculation. Although the damping of the backflow still does not seem to be sufficient the trend of the results is significant. An effect which is neglected in the present study is the influence of the taper on the velocity profile. As pointed out earlier separation occurs even in steady flow in a divergent pipe. If in addition the fluid is decelerated, in other words the adverse pressure gradient increased, separation may occur even in very slightly tapered tubes and increase the energy dissipation during the backflow. Streeter(21) has taken into account such asymmetric friction by assigning different empirical friction factors for forward and backward flow, but did not consider the influence of the acceleration.

VIII. SUMMARY AND CONCLUSIONS The objective of this thesis was to include the effect of the varying velocity profile of unsteady laminar pipe flow into the one-dimensional model of transient flow. It has been shown that the problem can be reduced to choosing an expression for the wall shear stress which is based on the true velocity distribution. For sinusoidal flow, an equation which relates the wall shear stress to the instantaneous values of mean velocity and acceleration is applied to the RLC model of steady oscillations in pipes. The procedure results in frequency dependent terms for the resistance R and the inertance L. A solution of the equations with the established impedance method is possible without further modifications of the method. For a transient variation of the flow, an equation is derived which relates the wall shear stress to the instantaneous mean velocity and to the weighted past velocity changes. The term is applied to the method of characteristics in which the past velocities at every point are known. Experimental verification of the theory is given by comparing measured and calculated pressure amplification due to resonance in a tube, and by predicting the pressure fluctuation in a pipe due to a sudden valve closure. A numerical study demonstrates how frequency dependent friction alters the response of a pipe if the input contains high frequency components. -77

-78A separate chapter applies unsteady friction to the flow of blood in arteries and shows various results of flow curves using turbulent friction, steady and unsteady laminar friction, and gives a special approximate equation for unsteady friction. On the basis of experimental results and numerical calculations, one can conclude that the traditional application of steady state friction to unsteady flow problems is insufficient for viscous flow with high frequency components. Earlier methods of other authors, which consider frequency dependent effects, are all based on the classical linear theory of transfer operators for pressure and flow at two cross sections of the pipe. The method presented in this study uses the transfer technique only for periodic flow but applies the method of characteristics to the transient case. The employment of the developed theory is dependent on the use of a computer, however, the method has some distinct advantages: a. If transfer operators are used, all input functions must be expressed by a summation of step functions, and in addition all reflections have to be taken into account individually. This is somewhat tedious for complicated systems consisting of series, branches or parallel pipes. In the method of characteristics, however, all reflections are considered automatically, and complicated systems can be handled without significantly increasing the difficulty. b. Nonlinear boundary conditions which must be linearized in the transfer techniques can easily be handled.

-79 c. In the classical linear methods the equations had to be linearized and the convective terms, V - and H V, in ax ax particular were neglected. These terms are important if the fluid velocity and wave velocity are of the same order of magnitude, for example in pipes with very flexible walls. The method of characteristics includes these nonlinear terms. d. If an approximation for the wall shear stress of unsteady turbulent flow can be found, the method of characteristics offers the possibility of considering the nonlinear terms analytically. e. The treatment of an arbitrarily tapered tube is possible, if the taper is not so great that it affects the velocity distribution. This limitation was adopted by all earlier investigators. No restriction is made with respect to the duration of the pulse. The results of this thesis are valid for laminar flow and only a criterion of the laminar-turbulent transitions could describe the range of applicability of the theory. No reliable Reynolds number for unsteady flow has yet been developed, although some attempts are described in the literature review, Chapter II. Similarly, the problem of frequency dependent friction in fully developed turbulent flow is still unsolved but of highest importance. The literature available indicates that, even in turbulent flow, the friction is dependent upon the history of the motion. The way in which laminar friction is applied to the method of characteristics in

-80Chapter VI could be extended to turbulent friction if realistic weights for the velocity changes could be determined. Due to the strong lateral diffusive effect of the eddies, the time history of the velocity changes which actually has an influence on the instantaneous friction, seems to be much shorter for turbulent than for laminar flow. This makes the proposed method even more valuable, since it would reduce the necessary calculations. The topic should be subjected to an intensive experimental study.

APPENDIX DISCUSSION OF THE FUNCTION l(ac i-i), a REAL In Chapter IV l (a T-i) was defined as 1(aCT-i) = JO(aT-i) - a N-i. To calculate numerical values for G71 one can J (a 4-i) follow Kelvin and introduce the functions ber c, bei a, beri, beilea (Bessel-real and Bessel-imaginary) which stand for the real and imaginary parts of Jo and J1 and use T-i = -I ( - 1 + i) \T2 Then ( 1 (a -i a. lber a + i bei a 1= (-l+) berla + i beila ~71 can be separated into the real and imaginary parts, which, in most cases is not necessary, because the numerical calculation can be carried out in the complex plane. U( J-i-) = a ber a beila - ber a ber1a - bei a beila-bei a berla ber a + bei a 1 1 a ber a berla + bei abeila- bei a bera +ber a beila + i 2 2 NT2 berla + bei a The Kelvin series are absolutely and uniformly convergent series for any real a. They may be differentiated and integrated term by term. They are(15) -81

-82 () 4j ber a -= + 2T2 2 2 ( )!] 2 (-1! (_) J-i c bei a = (-1) Ji [(23-1)]2 2! 2 4 (22 ()4 1 22 1' 2! 2! 35 r f 2 (C4 beila = -- -! - + 1 1! 22 2'.1 3' 4 1-+ 2 2 4 22 a6 (a) + | 3-.,. 3! 4! (a)6 _!_ + 4. 3: 4! (8 4' 62 8' *0 If a is greater than 10 a progressive loss of accuracy occurs if the above series are used to calculate VL, since the first terms become very large. It is more precise to apply the asymptotic expansions for Kelvin functions with large arguments.(15) However, for most practical cases it is sufficient to calculate V 1 from the expansions Jo (a# -i) J1 (a4 -i) a 2 V a4W-i - 2 it a -i cos (a — i - ) cos (a T-i - 3 n) 4 which are the more accurate the more a increases. It follows that 7J l(a(-i) = - aT-i tan (a4-i - ) 1 4~~~~~~~~~~~ -= (1-i) s2 - sin (-'~T2 a + it) + i sinh (/T2 a) cos (,2 a) + cosh ( $2 a)

-83for a > 10 the fraction approaches i and therefore 7 (O-i) a, (l + i) 1 $2 If a is less than about 0.1, 71(C$ -i) can be calculated from the first terms of a power series near the origin 2 1(C i) =2 + i a4

REFERENCES 1. Allievi, L. "Teoria Generale Del Moto Perturbate Dell Acqua Nei Tubi in Pressione," Annali Della Societa Degli Ingegneri ed Architetti Italiani, Milan, 1903. English Translation: Theory of Waterhammer, Riccardo Garoni, Rome, Italy, 1925. 2. Blasius, H. "Laminare Stromung in Kanalen wechselnder Breite," Z. Math. u. Physik, 58, (1910), p. 225. 3. Brown, F. T. "The Transient Response of Fluid Lines," Journal of Basic Engineering, Trans. ASME, Series D, 84, (1962), p 547. 4. Brown, F. T., Nelson,S. E. "Step Responses of Liquid Lines with Frequency-Dependent Effects of Viscosity," Journal of Basic Engineering, Trans. ASME, Series D, 87, No. 2, (1965), pp5O04-510. 5. Carstens, M. R., Roller, J. E. "Boundary Shear Stress in Unsteady Turbulent Pipe Flow," Journal of the Hydraulic Division, ASCE, February 1959. 6. Daily, J. W., Hankey, W. L., Olive, R. W., Jordaun, J. M. "Resistance Coefficients for Accelerated and Decelerated Flows through Smooth Tubes and Orifices," ASME Paper, No. 55-SA78, 1955. 7. Doetsch, G. Guide to the Applications of Laplace Transforms, D. van Nostrand, London; Toronto, Canada, New York, N.Y., Princeton, N.J., 1961. 8. D'Souza, A. F., Oldenburger, R. "Dynamic Response of Fluid Lines," Journal of Basic Engineering, Trans. ASME, Series D, 86, No. 3, (1964), pp. 589-598. 9. Fan, C. Non-steady, Viscous, Incompressible Flow in Cylindrical and Rectangular Conduits (with Emphasis on Periodically Oscillating Flow), Doctoral Thesis, University of Illinois, 1958. 10. Gilbrech, D. A., Combs, G. D. Critical Reynolds Numbers for Incompressible Pulsating Flow in Tubes, Proceedings of the First Southeastern Conference on Theoretical and Applied Mechanics, (1962), p. 292. 11. Holmboe, E. L. Viscous Distortion in Wave Propagation as Applied to Waterhammer and Short Pulses, Doctoral Thesis, Carnegie Institute of Technology, 1964. 12. Jahnke, E., Emde, F. Tables of Functions, Dover Publications, New York, 1945. 13. Joukowsky, N. Uber den hydraulischen Stoss in Wasserleitungsrohren, Voss, Petersburg und Leipzig, 1900. English Translation: "Waterhammer", Proc. Amer. Water Works Assoc., 24, (1904), pp. 341-424. -84

-85 14. Lambossy, P. "OscillationsForcees d'un Liquide Incompressible et Visqueux dans un Tube Rigide et Horizontal. Calcul de la Force de Frottement," Helv. Phys. Acta, 25, 1952. 15. McLachlan, N. W. Bessel Functions for Engineers. Clarendon Press, Oxford, 1955. 16. Onoc, M. Tables of Modified Quotients of Bessel Functions of the First Kind for Real and Imaginary Arguments. Columbia University Press, New York, 1958. 17. Richardson, E. G., Tyler, E. "The Transverse Velocity Gradient Near the Mouth of Pipes in which an Alternating or Continuous Flow of Air is Established," Proc. Phys. Soc. London, 42, (1929), p. 1. 18. Rouleau, W. T., Young, F. J., "Distortion of Short Pulses in Tapered Tube Pulse Transformers. Part II - Viscid Liquid," Journal of Basic Engineering, Trans. ASME, Series D, 87, (1965), pp. 471-477. 19. Schultz-Grunow, Fo "Pulsierender Durchfluss durch Rohre," Forsch. Ing.-Wesen, 170, 1960. 20. Sexl, T.,'ber den von E. G. Richardson entdeckten Annulareffekt," Z. Physik, 61, (1930), p. 349. 21. Streeter, V. L., Keitzer, W. F., Bohr, D. F. "Energy Dissipation in Pulsatile Flow Through Distensible Tapered Vessels," in Attinger, F.P., ed., Pulsatile Blood Flow, Blakistan Div., McGraw-Hill, New York, 1964. 22. Streeter, V. L. Computer Solutions of Surge Problems, Symposium on Surges in Pipelines, Institution of Mechanical Engineers, Proc. 1965-66, Vol. 180, Part 3E. 23. Streeter, V. L., Keitzer, W. F., Bohr, D. F., "Pulsatile Pressure and Flow Through Distensible Vessels," Circulation Research, 13, (1963), pp. 3-20. 24. Streeter, V. L., Lai, Ch. "Waterhammer Analysis Including Fluid Friction," Journal of the Hydraulic Division, ASCE, May 1962. 25. Szymanski, F. "Quelques Solutions Exactes des Equations de 11Hydrodynamique de Fluide Visqueux dans le Cas d'un Tube Cylindrique," J. de math. pures et appliquees, series 9, 11, (1932), p. 67. 26. Tarantine, F. J., Rouleau, W. T. "Fluid Pressure Transients in a Tapered Transmission Line," Trans. ASME, Winter Annual Meeting, New York, Nov. 1966. 27. Uchida, S., "The Pulsating Viscous Flow Superimposed on the Steady Laminar Motion of Incompressible Fluid in a Circular Pipe," ZAMP VII, (1956), pp. 403-422.

-8628. Wylie, E. B. "Resonance in Pressurized Piping Systems," Journal of Basic Engineering, Trans. ASME, Series D, 87, No. 2, (1965) pp. 960-966. 29. Wylie, E. B. Flow Through Tapered Tubes with Nonlinear Wall Properties, Symposium on Biomechanics, ASME, New York, Nov. 1966. 30. Womersley, J. R. An Elastic Tube Theory of Pulse Transmission and Oscillatory Flow in Mammalian Arteries, Wright Air Development Center, WADC Report TR 56-614, 1957.