THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING PULSATILE PRESSURE AND FLOW THROUGH DISTENSIBLE VESSELS Victor L, Streeter, Sc.oD W, Ford Keitzer, MoDo David F, Bohr, MoDo December, 1962 IP-595

TABLE OF CONTENTS Page I, INTRODUCTION....1........................................ 1 II, THEORETICAL METHODS AND RESULTS........................... 6 A. Derivation of Basic Equations......................... 6 1. Basic Assumptions............................... 6 2. Elasticity Relationships......................... 7 3. Continuity Equation.............................. 10 4, Momentum Equation............................... 11 5. Development of Characteristic Equations........... 14 6. The Method of Specified Time Intervals............ 16 7. Boundary Conditions.............................. 19 8. Example.......................................... 22 B. Equations for Tapered Distensible Vessel With Distributed Outflow........................................ 27 III, COMPARISON OF MEASURED FLOW WITH FLOW COMPUTED FROM PRESSURE-TIME DATA IN A TAPERING VESSEL........................ 31 SUMMARY............................................... 34 ACKNOWLEDGEMENTSo.................... 35 NOTATION................................ 36 NOTATION FOR MAD PROGRAM.................................... 38 REFERENCES........................................ 40 ii

LIST OF FIGURES Figure Page 1 Relation Between Tensile Force Change in Wall dT and Pressure Change dP,.o.......................... 12 2 Material Balance Showing Mass Per Unit Time Entering and Leaving an Elemental Volume,.......... 12 3 Force Acting on Segment of Fluid (Solid Line), and Momentum Relationships (Dotted Line).............. 12 4 Characteristic Lines C+ and C_ Drawn Through Points R and S Respectively, Where V and H are Known, Their Intersection P is a Point Where x, t, V, and H May Be Determined From the Equations.... 17 5 By Specifying At and Ax, Which Locates Point P, Values of V and H at R and S May Be Found by Interpolation From Values at A, C, and B................o 17 b Proximal Boundary Relations. One Relation Between V and H at P is Obtained from Values of V and H at S................................................ 21 7 Distal Boundary Relations........................... 21 8 Computer Program and Sample of Calculations for Pulse Flow Into a Flexible Tube. Pulse Time 0,4 sec, Unstressed Tube i.do = 1.3 cm. t' = 0,l cm, Y = 5. x 106 dynes/cm2 Problem was started as a steady state flow with flow of 35.65 cc/sec. and head of 14, cm Hg at proximal end, Friction factor F = 0,3 and losses varying as square of the velocityo The calculations are printed out for the third pulse after steady flow, with values given for each 0o008 sec. time interval............... 23 9 Tapered Vessel with Distributed Outflow Along the Walls. Flow Per Unit Length Through i-th Reach Shown. HB is Terminal Bed Pressure................ 30 10 Pressure-Time Graphs for Two Cross Sections of the Femoral Artery of a Dog. Distance Between Sections 9.4 cm, Proximal Diameter 05316 cm. Distal Diameter 0,228 cm, Branch Arteries Between the Two Sections are Ligated....................................... 50 iii

LIST OF FIGURES CONT'D Figure Page 11 Electromagnetic Flow Meter Data For the Upstream Section of the Femoral Artery (Figure 10) and Computed Flow Using the Pressure-Time Data of Figure 10 for Boundary Conditionso In the computer solution the speed of pulse wave at H = 11,5 cm Hg was a = 1200 cm/sec. The friction factor was f = 0.4 and the energy dissipation was taken to vary as the square of the velocity.............. 32 iv

I INTRODUCTION The problems of characterizing pulsatile patterns of pressure and flow in the arterial system are intriguing and complex. Ejection of blood from the left ventricle initiates non-linear transients in pressure and flow at the root of the aorta. These transients initiate complex pulse patterns that are propagated throughout the arterial tree. Some of the factors that must be reckoned with in an analysis of these patterns are itemized in the list that follows: 1) The force initiating the transients is itself complex; the velocity of ventricular ejection increases rapidly with the opening of the aortic valve, then declines slowly to reach a negative nadir with the closure of the aortic valve. 2) The distensibility of the walls of the arteries receiving this positive increment of pressure and flow has an important influence on pulse patternso This physical property of the arterial wall is also responsible for changes in configuration and velocity of the transients as they pass over the arterial tree. The pulsatile patterns are further distorted by 3) fractional losses of both positive and negative flow of the blood, and by the 4) branching and tapering architecture of the arterial tree. -1

-25) The resistances to forward motion of blood through the distal arteriolar beds also have their effects on the contours of the arterial pulses observed upstream. Both experimental and theoretical methods have been used to analyze and to quantify the influence of these factors on the arterial pulse patterno The experimental approach has profited substantially from new instrumentation which permits a recording of the transients of both pressure and flow with a high degree of fidelity at various levels of the arterial tree(1l'2'34). The theoretical approach employs established mathematical relations between known physical parameters which permit the prediction of pressure and flow at specific points in the arterial tree and specific times in the cardiac cycleo If an adequate mathematical expression were available to describe these relations, the fit of these predictions to the actual measured values would then become a valuable tool both for assessing the accuracy of the selected values for the physical parameters and for checking the significance to the various factors that influence the transients in the arterial tree. The theoretical analysis of these transients then has two requirements: 1) realistic values for these factors that influence the transients, and 2) a mathematical statement of the inter-relationship of these factors which will define pressure and flow with respect to time and position in the arterial tree. Approximations of the required quantitative values for the physical factors involved can be found in the literature. However, one of the more unyielding problems has been the development of a mathematical expression for the inter-relation of these factors

5-3 in pulsatile flow in a distensible vessel. Owing to the difficulties encountered, greatly simplifying assumptions have been made, such as lumped parameters, laminar flow, linear frictional resistance, and steady oscillatory flow(l12) The resulting equations, despite the restrictive assumptions, are very complicated and do not lend themselves to ease of computation when practical boundary conditions are introduced. The current study presents a new mathematical approach which permits a more realistic solution of specific equations describing these inter-relations. Specific values for pressure and flow with respect to time and position in the arterial tree can be computed. This approach has not previously been applied to a study of the transients in blood vesselso It starts with two established equations dealing with these transients and the physical factors that influence them. These simultaneous equations are: 1) the continuity equation which equates the net influx of blood entering a small segment of the arterial tree with the increase of volume of that segment, and 2) the momentum equation (Newton's second law) which equates the forward force acting on this segment of blood with the backward force plus the force exerted by the arterial wall and the force required to overcome friction in the artery (Figure 3)5 The inclusion in this equation of the statement for friction makes it a non-linear, partial differential equation which could not be solved by previous methodso The new aspect of the current approach is the application of the method of characteristics which permits the solution of these two simultaneous equations. Specific values can be computed for the unknown functions of pressure and velocity along

-4characteristic lines relating the independent variables of time and distance (Figures 4 and 5)~ With the aid of a high speed computer these unknowns have been determined for frequent periods of time (1/400'th of a physiological pulse cycle) and for short segments along the arterial tree (1 cm)o The problem of reflections and their interpretation in an analysis of unsteady flow become less significant with this approach. A pressure pulse is transmitted through the vessel at a speed that depends upon the tube properties and the pressure within the tubeo As the pressure varies both with distance along the tube and with time, the speed of the pulse wave changes continuously with the independent variables, distance along the tube and time. For any change in speed of the pulse wave, reflections are set up which move in the reverse direction, and the transmitted wave is affected was well. As reflections also occur at boundaries, ioe., entrances, exits, branches and obstructions, previous methods of keeping track of all these reflections became hopelessly complicatedo The new method outlined here automatically takes all these reflections into account, not by keeping track of them separately, but by satisfying all of the basic equations at closely-spaced sections along the vessel at frequent time intervals, and by satisfying the boundary equations. Owing to the relative ease of handling the computations for interior sections and for boundaries, solutions may be programmed simulating branching arteries with satisfactory accuracyo

-5The theory of characteristics, which applies to the solution of hyperbolic partial differential equations, first gained prominence in solution of supersonic flow problems by Courant and Friedricks.(5) These methods were extended to applications of free-surface flow cases later, by Stoker.(6) More recently they have been applied to water hammer situations by Lai(7), Streeter and Lai(8), and Streeter(9), in which non-linear terms for wall expansion and for wall friction have been retained in the equations. This paper presents an extension of the theory to flexible vessels, tapering vessels, and vessels with distributed outflow along their length, together with in vivo experimental confirmation of the method.

II. THEORETICAL METHODS AND RESULTS In this section the mathematical and physical relationships for flow through flexible tubes are derived, including equations for tapering vessels with distributed outflow. A. Derivation of Basic Equations The assumptions required for analytical handling of pulsatile flow are first discussed, followed by development of elasticity relations and the continuity and momentum equations for vessels of constant initial diameter. From these relationships the characteristic equations are obtained and finite difference methods applied to develop the equations for the method of specified time intervalso After discussion of boundary conditions, an example is presented. 1. Basic Assumptions The basic assumptions required in developing the working equations are: ao One-dimensional flow; the velocity at a cross section is given by the average velocity at the section at a given instant. b. The vessel walls are elastic, with a Poisson ratio of 0,5, and they are tethered; hence the volume of elastic vessel wall per unit length remains constant. c. The fluid density is constant. Compressibility of blood is small compared with expansion of the vessel walls under increased pressure, and -6

-7d. Pressure losses due to wall friction may be expressed as proportional to some power of the velocity at a cross section. 20 Elasticity Relationships Although arteries have viscoelastic properties, their effect seems to be minor and previous investigators(l0) have concluded that the stress-strain curve may be linearized without introducing appreciable error. If D is the inside diameter of the vessel, and t' the effective wall thickness, then as a consequence of the constant volume of elastic wall material resulting from the assumption of Poisson's ratio of 0.50, Dt' = Doto (1) Do is the unstressed diameter and to is the corresponding wall thickness. If the tube is subjected to an increment of pressure dP internally, the tensile force dT resisting this pressure increment per unit length of vessel, Figure 1, is dT - d(PD) 2 Since the diameter change, on a percentage basis, is much less than the pressure change, the term dD/2 is neglected in expanding the right hand side, leaving dT = dP (2) 2 By dividing through by the wall thickness the change in tensile stress dS (force per unit area) is obtained dT DdP dS - = - t' 2t'

-8Now, by dividing through by the elastic modulus of the vessel wall, Y, the unit strain is obtained, which is the change in length per unit length caused by dPo Since circumference changes are proportional to diameter changes, the unit strain is dD/D, dD D dP (3) D 2t' Y After use of Equation (1) to eliminate t, and after separating variables dP = 2Yt D (4) 0 D3 Integrating, 1 1 P= YtD - (5) 0 o D0 D o in which the condition has been used that D = Do when P = 0. By multiplying and dividing the right-hand side by nt/4 to introduce the vessel cross sectional area A, and correspondingly Ao, A 1 1 (6) 0 1 -DO t' Y 0 The pressure P has been expressed as force per unit area (dynes/cm2). It is customary to express it in terms of the height of a liquid column (the fluid flowing). These are related by P = pgH in which p is the mass density (gm/cc), g is gravity (980 cm/sec2), H is the pressure in height of fluid flowing (cm)o Two additional substitutions are made, let a2 = to Y a2 t 0 P D P D7

-9Then A D2 1 A D2 21I! (8) A 0 o 1 - gH/a2 By use of Equations (1) and (7) a = to D D2 A a2 t' Do D2 Ao After equating expressions (8) and (9) 2 2 a = a gH (10) 0 and from Equation (9) Do ao D = ~a (11) a a is the speed of the pressure pulse wave through the vesselo It changes both with time t and with distance x along the vessel, and whenever a changes reflections are produced. The partial derivatives of A with respect to the two independent variables, time t and distance x are needed later and result from Equations (8) to (11) (a variable subscript x or t represents partial differentiation with respect to that variable, i.e., Ax = A/6xo A a2 02 2 A = 2 a At = a2 0 a5 o a t Ao a3 Ao a3 In these equations A is considered to be constant. Then, from Equation (10) 2a a =gHX 2a at = gH x -gHx = -grt

-10Now, by eliminating ax, at, and Ao from the last three sets of equations Ax - gHx A gHt ( (12) A a2 A a2 Also, from the relation P = pgH P = PgHx (1 These elastic relationships are used in the continuity and momentum equations that are now developed. 3. Continuity Equation The continuity eq-ation is a material balance for a small segment of vessel, which states that the net mass inflow per unit time is just equal to the time rate of increase of mass within the segment, Figure 2. If V is the average velocity at the entrance to the element, the rate of mass inflow is pAV, and the rate of mass outflow is pAV + (pAV)xdx. The net mass inflow per unit time is then - (pAV)xdx and must just equal the time rate of increase of mass within the segment (pAdx)t. Equating these expressions (pAV)dx = (pAdx)t (14) p is constant for all practical purposes when considering flow in a distensible vessel. A and V are dependent variables. After expanding Equation (14), remembering that x is independent of t, and then dividing through by the mass of the fluid segment, p Adx, VA + At 15 - + - (15) A x A which is the continuity equation and must hold throughout the vessel.

-114. Momentum Equation The momentum equation when applied in the x-direction to the fluid in a segmental volume, Figure 3, is a statement that the resultant x-component of force on the segment of fluid is just equal to the net efflux of x-momentum plus the time rate of increase of x-momentum within the segment. The resultant force component, Figure 35 is dx F = PA + (P + Px 2 ) Axdx - [PA + (PA)xdx] - ToDDdx The first term is due to pressure within the fluid acting over the cross section at x, the second term is the force component in the x-direction due to the tube wall pushing against the fluid (zero for constant A, as A dx is the increase in cross section in the length dx)o The term in x brackets is the force pushing against the element on the distal side. The action of fluid friction at the tube wall is given by the product of shear stress 0 at the wall and area of wall surface tDdxo This force is assumed to act wholly in the x-directionn By expanding the expression for F it becomes F =-PxAdx - ToDdx + P A (2d Since the term with square of dx becomes of a higher order of smallness as dx approaches zero, it may be dropped from the equation. The momentum influx at x is oAV2 and the momentum efflux at x + dx is pAV2 + (pAV2)xdx, with a net efflux of (pAV2)xdxo The time rate of increase of momentum within the segment is given by (pAdxV)to After combining the force and momentum terms9

-12T +dT D ) PD+d(PD) d ^ i T+dT Figure 1. Relation Between Tensile Force Change in Wall dT and Pressure Change dP. PAV --' - PAV+ (pAV)xdX x ldx - figure 2. Material Balance Showing Mass Per Unit Time Entering and Leaving an Elemental Volume. /-TIME RATE OF INCREASE (P+PxdX/2)AXdX OF MOMENTUM I (PAVdX) PA (P PA +(PA)dX - - MOMENTUM IN I MOMENTUM OUT PAV2 PAV2+ (pAV2)xdX ToTDdX X dX Figure 3. Force Acting on Segment of Fluid (Solid Lines), and Momentum Relationships (Dotted Lines).

-135- PxAdx - To0TDdx = (pAV2)xdx + (pAdxV)t By expanding the partial derivatives, then dividing through by the mass of the segment pAdx, and after replacing Px by pgHx, TotD A 2 At gHx + t + + 2VV + V +V+ V 0 (16) pA A t A The wall shear stress To may be written in the form T = k2 (17) o 2 For established, steady laminar flow k = 16/R, with R the Reynolds number VDp/p, in which p is the fluid viscosity (Poise). For turbulent flow k = f/4, with f the commonly used Darcy-Weisbach friction factor. (11) By inserting Equation (17) into Equation (16), using f, gH + + x V2 + 2V Vx + Vt + V At = 0 (18) xy 2D A t A This is the momentum equation for flow through a distensible vessel. By multiplying Equation (15) by V and subtracting it from Equation (18), substantial simplification results* fV2 gHx + V Vx + Vt + - = 0 (19) 2D Equations (15) and (19) contain the continuity and momentum principles. After substituting the elastic relationships given by Equations (12) into Equation (15) it becomes, upon simplification, a2 VHx+ Ht +-V = (20) g X *By writing the V2 term in the friction expression as V V, it will reverse sign if the flow reverses, and hence act in the proper direction at all times o

-14Equations (19) and (20) are used to develope the final equations. 5. Development of Characteristic Equations By calling Equation (19) L1 and Equation (20) L2, they may be combined linearly using an unknown multiplier X, as follows: g1 k XL2 = X[H,( ha2 fV2 L L1 + L = [Hx( + V) + Ht] + [Vx(V + — ) + Vt] + = 0 (21) X -) V] 2D If two distinct values of X are taken, two equations result which contain the momentum and continuity principleso The theory of characteristics determines two special values of X which result in great simplification of the equations. To review some fundamental relations in calculus, if H = H(x,t) and V = V(x,t), then the total derivatives of H and V with respect to t are dH dH x t V V dx + V dt x dt dt x dt where in these relationships H and V are pressure and velocity of a particle as it moves (x becomes a function of t). By examining Equation (21) it is seen that the first bracket contains dH/dt if g +v = d (22) X dt and the second bracket contains dV/dt if Xa2 dx V + Xa2 = (23) g dt These expressions must be the same. By equating them and solving for X g v a2 + V = V+ -- x g

-15and = a (25) -a Now, by restricting the applicability of Equation (21) to those characteristic lines for which Equations (22) and (23) are satisfied, it may be written in the simple form dH dV fl2 L = dH + + -— = O (25) dt dt 2D By applying X = +g/a to Equations (22) and (25) g dH + V + fV(26) a dt dt 2D ])C d x = V + a (27) dt j and by applying X = -g/a to the same equations -g dH dV + T28) -~ ^ ^+l ~v+ ~_i= 0 (28) a dt dt 2D =! C dx = V- a I (29) dt i Equation (26), a total differential equation, is valid only along a line defined by Equation (27), if a plot is made on an x - t plane, as shown in Figure 4. If C in the figure represents the characteristic line defined by dx/dt = V + a and passing through a point R where V, H, x and t are known then Equation (26) may be used as one relation between pressure H and velocity V along this line. Similarly Equation (28) is valid only along a line defined by Equation (29), Figure 4o Now with V and H known at the known points R and S, four

-16equations are available at the intersection P of the two characteristic curves, for computing V, H, x and to 6. The Method of Specified Time Intervals For computation purposes, Equations (26) to (29) are written as finite difference equations. For use with a high-speed digital computer, the theory of characteristics method may be extended to a system in which the time and distance intervals are preselectedo(12) This is called the method of specified time intervals, and entails interpolation of values of V and H at unknown points R and S, Figure 5, such that P occurs at equally spaced distances Ax along the vessel for equal time increments. In Figure 5 consider that V and H have been computed for the first two rows at the equally spaced sectionso Values of V and H are to be computed next for point P at time tp, which is tC + At, and values of V and H are known at A, B, and Co Equations (26) to (29) are written as finite difference equations, for points R and P on C+ and for points S and P on C_, with At = tp - tR tp - tSo The mesh (Ax, At) is assumed to be so fine that the velocity in the friction term may be evaluated at known conditions at Co Also aR and aS are replaced by ac, as well as DR and DS by DCo g (Hp HR) + Vp VR + fcVC VC At 0 (3) a 2DC xp - xR = (V + a)C At (31) -& (Hp - HS) + Vp - VS + fcC Ivcl At = (32) aC 2DC

-17t P CC+ ~/ ~ S R l _ X Figure 4. Characteristic Lines C+ and C Drawn Through Points R and S Respectively, Where V and H are Known. Their Intersection P is a Point Where x, t, V, and H May Be Determined From the Equations. P At B -- -- - ---- --- At A R C S B ~~~~~~~~~~0t' 0 Figure 5. By Specifying At and Ax, Which Locates Point P, Values of V and H at R and S May Be Found by Interpolation From Values at A, C, and B.

-18xp - XS = (V - a)C At (33) The V2 friction term has been replaced by VI V| VR, HR, VS, and HS are now to be computed from the equations, with reference to Figure 5, by linear interpolation. The mesh (values of Ax, At) is assumed to be fine enough so that the slopes of the characteristics lines are given adequately by evaluating V + a and V - a at C. From Figure 5 VR- VA XR - XR A (x4) Vc - VA xc - XA Ax By remembering that xp = Xc, from Equation (31) XC - xR = (V + a)C At Then xR - xA = Ax - (xC - XR) = x - (V + a)C At (35) For convenience the grid mesh ratio At/Ax is called 9o Substituting into Equation (35), the first of the following four equations is obtained, VR = VA + (VC - VA)(1 - @(V + a)C) (36) HR = HA + (HC - HA)(1 - G(V + a)) (37) VS = VB + (VC - Vg)(l + G(V - a)C) (38) HS = HB + (HC - Hg)(l + G(V - a)C) (39) The last three equations are found in a manner similar to that used in obtaining Equation (36). In

-19= At (40) Ax At must be selected so that R and S lie within the reach defined by points A and Bo With the interpolated values of VR, HR, VS, HS known, Equations (30) and (32) may now be solved for Hp and Vp, Hp= R S+a (VR -Vs) (41) 2 2aC V + V. o + g fCVC lVC At PVPF - -a + - (HR - Hs) - (42) 2 2aC 2DC Equations (36)through (41) permit all interior points of the grid to be computed, i.eo, all points where corresponding points, A, B, and C are known. At the end points an additional condition is needed to solve for Vpand Hp, since only one of the Equations (30) and (32) is available at a given boundary. This additional condition is called the boundary condition. 7. Boundary Conditions Boundary conditions may take on many forms. Some examples are given here to illustrate procedures in developing them~ At the proximal end of the vessel Equation (32) applies, Figure 6. fcVc VVc At Vp VS+ (HP- HS) - (43) aC 2DC in which Vp and Hp are the unknowns, as VS and HS are given by the interpolation Equations (38) and (39). If a known pulse inflow Qp into the vessel is known at this instant, then a second equation becomes available, as follows:

-20 - Qp Vp Dp 2 (44) But from Equations (10) and (11) 2 D2a2 Do2a2 Dp 2- - = -p —-- (45) a a07 - gHp By substituting Equation (45) into Equation (44) to eliminate Dp, and after solving for Vp, 4Qp 2 Vp — 2 (ao - gHp) (46) rtDo ao which is the second relationship required, With the pulse inflow given for each time increment, Vp and Hp at the proximal end may be computed progressively as the rest of the solution proceeds. Another example is to have the pressure pulse specified as a function of time at x = 0. In this case the known Hp is inserted into Equation (43) and Vp may be computed. At the distal end, Figure 7, the outflow may be expressed in terms of the pressure difference between the vessel Hp and the terminal bed HB as Qp = k1 (Hp - HB)m Vp Dp2 (47) Equation (30), solved for Vp is VP - VR - (Hp - HR) - fCvc IC (48) V VR aC 2DC

-21P CC 0 0 S B Figure 6. Proximal Boundary Relations. One Relation Between V and H at P is Obtained from Values of V and H at S. P / C A R Figure 7. Distal Boundary Relations.

-22By use of Equations (10) and (11), Equation (47) may be solved for Vp in terms of Hp 4k 2 V = 2 2 (Hp - HB) (a - gHp) (49) ~Do ao A trial solution may be required, depending upon the value of exponent m in the terminal bed resistance relationship. Care must be taken in applying a terminal bed Q - H relation. Within the vessel proximal to the bed, there is no simple Q - H relation, primarily due to the complex reflections. If a segment of an arterial system is to be analyzed, the boundary conditions should generally be expressed as H versus t, Q versus t, or Vp versus t. For known pressure as a function of time, Vp is found directly from Equation (48). 8. Example To illustrate the application of the characteristics theory to a flow situation, a known pulse flow is injected into a distensible vessel, with a linear terminal bed at the distal end. The computer program, Figure 8, is in the MAD (Michigan Algorithmic Decoder) language and an IBM 709 is used in performing the calculations. The pulse flow from the heart, measured at the ascending aorta of a dog using the Square Wave electrio magnetic flowmeter,* has been expressed by a series of empirical formulas. The average flow is 35565 cc/sec over the 0o4 sec. pulse cycle, with a maximum inflow of 160 cc/seco A friction factor **Carolina Medical Electronics, Inc., Winston-Salem, North Carolina.

-23PULSATILE FLOW THROUGH A DISTENSIBLE TUBE DIMENSION V(20),VP(20),H(20),HP(20),D(20),Q(20),AA(20) *001 INTEGER IU,NPKK *OG2 PRiNT COMMEN$iS.. GIVEN DATA FOR PROBLEMS... OU3 Al READ DATA _004 PRINT RESULTS Y,THICKDO,L;HBHOQAVERHOAV' R O GSG,FNDELTPTIME 005 2,PKK *005 PRELIMINARY CONSIDERATIONS TM=.12~PTIME ___006 PMT~.25*PTIME "007 PPT~.32*PTIME.008 PTT.34MPTIME *-009 OM=160..010 AO.SORT.(Y*THICK/(7RHODO)).011 VO=.7854*DO*UO*AO.AO >012 TH=DELT*N/L *013 A=SQRT.(AO~AO-G*SG*HO) 0l14 DUDe0AO/A t015 VP=QAVE/(.7854D0~D ) -016 DHF=F*L*VP.VP/ (N**2.*G*SG) *017 THROUGH A2,FOR I=0,1,I.G.N *018 V(I)=VP *019 H(I)=HO-UHF.I *020 0(11=0.021 0(I)=OAVE *022 A-2.....A IGS..A( -I)=-SQRT.IA(- -AO-G-S'G'H(I )... —---—. —— 3 — ---------------------- --- 23 —K=QAVE/(HO-HB) -024 T=-5.aDELT,025 U=-5 ~026 C1=K/(.7854.DODDO*AO*AO) -027 C2=C.l(AO*AO+G*SG*HB) *028 """"" C- -- ---------------------------------------- ------ - --- -------- 02...... C3=C1.G.SG.029 C4=C1*HB*AO*AO _030 ACN=SQRT.(AO*AO-G*SG~H(N)) 31 PRINT COMMENT.O HEADS(CM HG)~ VELOCIIIES(CM/.032 2SEC), DIAMETERS(CMI, AND FLOW(CC/SEC)$ -032 PRINT COMMENTSO TIME FLOW X/L=.0.1.2 *033 " -. —- -. ^ ^. -6 -- -- --------—. —-- ^^ —.2.3.4.5.6,7.8.9...033 3$ *033 A3 PRINT FORMAT BlTIOH(O),Hi2);H(4),H( 6)4)H(6H(8)t,H( 0),H(12)H(14 *034 2),H(16),H(18),H(20).034 -PRINI FORMAT B2,V(O) V(,V2),V (4),V(6), V(8 l),V10), V(2),VV035 216),V(18),V(20) *035. —--- """ FfiNTT" g k Vi —-----— ff( VD( 2)''D. 6) D( 8 );D( 1' 0'2);.- -----------------— 0 —.1.-3-6 216)vD(18),U(20) *036 PRINT FORMAT 84,Q(0)(Q(2)(2,Q(4)Q(6)Q(8)QQ(10),Q(12),Q(14),Q( ~037 216),U(18),Q(20) *037 VECTOR VALUES Bl=$IHO,F8.3,F9.2',S2',3H H=,llF8.3*.-038 VECTOR VALUES 82=$1H,S20,2HV=,11F8.3$ -*039 VECTOR VALUES 83=SIH,S20,2HD=,11F8.3t$ *040 VECTOR VALUES B4=$1H ~S20,2HQ=,llF8.3*S *041 A4 T=T+DELT -042 U=U+l -043 WHENEVER U.G.PTRANSFER TO Al *044 CALCULATION OF INTERIOR POINTS THROUGH A5,FOR I=t,1,I.E.N.045 HR=H(I-1)+(H(I)-H(I-1))*(1.-TH*(V(I)+AA(I))) -046 VR=V(I-'l )+(V(i)-V(I-1))*(l.-TH(V(II)+AA(I))). 047 VS=V(I+1)+(V(I)-V(I+1))~(1.+TH*(V(I)-AA(I))) -048 HS=H(i+1)+(H(I)-H(IY+1))(I.+THe(V(I)-AA(I))).049 HP(I)=(HR+HS)/2.+AA(I)~(VR-VS)/(2.*G*SG) 00__ VP(I)=(VR+VS)/2.+G*SG*(HR-HS)/(2.*AA(I))-F'V(I)*.ABS.V(I)DEL *051 2T/(2.0D(l)) -Ol5 D(I)=DO.AO/AA(I) ~053 A5 0(I)=.7854*Db(I).(I)*VP(I) *054 CALCULATION OF PROXIMAL BOUNDARY CONDIlION WHENEVER T.L.O.Q=OAVE *055 WHENEVER TGE.O..AND.T.LE.PMT,Q=QM*T/(TM*EXP.I(T/TM-i.)) ---.056 WHENEVER T.G.PMT.AND.T.LE.PPT,Q=QM*T/(TM.EXP.(T/TM-1.))-.56.0 *057 2M(iT-PMT)/(PPT-PMr)'...... —.........057 WHENEVER T.G.PPT.AND.T.LE.PTTtQ=QM.(T-PTT)/PPT -.058 WHENEVER T.G.PTT.AND.T.LE.PTIMEQ=O..059 WHENEVER.ABS.(T-PTIME).L..001~T=0. *060 AC=SQRT. (AO-AO-G-SG-H)- 061 VS=V(1)+(V-V(1))*(1.+TH*(V-AC)) -062 "-""""""HS=H(i-)+H-HI (1)).1. +TH.IV-AC ) 1- -- -063 HP=(Q/(.7854*DO.DO)-VS+G.SG*HS/AC+F*.V.ABS.V.DELT/(2.*D))/(G. *064 2SG-(1./AC+Q/VO)) *064 A=AOAO-GSGHP -- --- 65 VP=Q.A/VO -066 D=00*AO/SQRT.(A) _067 DISTAL BOUNDARY CONDITION VR= VIN-L)+ V(N)-V(N-1))I1. —TH_(VN)+AN ------------- - *068 HR=H(N-1)+(H(N)-H{N-1 )I 1.-THVINiN) +A CN)) - 069 VEE=VR+G.SG*HR/ACN-FV(N)*.ABS.VN)DELT/2..(N)) 0 —--------- 10 C5=G*SG/ACN *071 C6=(C5+C2)/C3 *072 C7=(VEE+C4)/C3 *073 HPIN)=C6/2.-SQRT.(C6BC6/4.-C7) *074 VP(N)=VEE-CS-HP(N) -075 ACN=SQRT.(AO.AO-GUSG*HP(N)) 076 D(N)=DO*AO/ACN - -77 Q(N)=VP(N)*.7854.D(N).D(N) ______ *U78 ___ THROUGH A6,FOR I=0,1,I.G.N -079 VII )=VPI 1) *080 A6 H(I)=HP(I) -081 WHENEVER U/KK.KK.E.UTRANSFER TO A3.082 TRANSFER TO A4 -083 END OF PROGRAM.__ "84 Figure 8 (Cont'd)

-24GIVEN DAIA FOR PROBLEM _............... Y_= 5.OOOOOOE 06,_._ THICK =.100000t, D = 1.300000,, = 50:.0000C HB. =...0... ].._O.O000000, _ HO = 14.000000, QAVE = 35.650000, RHO = 1.000001,............... G = 980.000000, SG = 13.600000, F =.300000, N = _____DELT = 4.00000E-03# PTIME =.400000, P = 300, KK = 2 HEADS(CM HG)L, VELOCITIES(CM/SEC), DIAMETERS(CM) AND FLOW(CC/SEC) TIME FLOW X/L=.0.1.2.3.4.5.6.7.8.9 1..000.00 H= 12.945 12.959 13.009 13.11Z 13.289 13.543 13.846 14.138 14.364 14.496 — 14.536 ___ __ _ ____V=.000 2.891 5.573 7.863 9.600 10.716 11.390 12.038 12.946 14.033 15.116 0= 1.751 1.751 1.754 1.760 1.770 1.785 1.802 1.820 1.834 1.843 1.845.........._ —=._._00_ 0 6.964 13.469 19.128 23.619 26.803 29.062 31.325 34.214 371428 40.421.008 61.36 H= 13.713 12.899 12.962 13.083 13.275 13.533 13.823 14.096 14.309 14.442 14.485.-.... -.................. -' ——.-;-_..... —-,672 4.748 4.542 4.910 6.45'3.9.029 12.022 15..000 _..... ____ __ 0= 1.794 1.748 1.752 1.758 1.769 1.784 1.801 1.818 1.831_ 1.839_J.82.._ O= 61.359 5.125 9.164 11.344 11.670 11.352 12.510 16.743 23.772 31.942 39.975.016 1603.88" H 14. —f4-.- -203- I-3.458 12.948 13.093 13.297 13.5407 3.793 02i 14.21i4 —-.34- 5 14.390 V= 39.743 18.362 1.591 1.083 -.154 -1.161 -.783 1.510 5.350 9.965 14.719.D 182.... 1.780 1.751 1.759 1.770 1.784 1.799 1.813 1.825 1.833 1.836 Q= 103.879 45.673 3.831 2.630 -.380 -2.903 -1.991 3.900 13.999 26.300 39.130.024 131.90 H= 14.531 14.047 13.366 13.134 13.339 13.553 13.756 13.940 14.098 14.210 14.253 V...... -33'-43-4381 10.433 -2.881 -4.875 -6.044 -5.389 -2.563 2.121 8.015 14.453 D= 1.845 1.815 1.774 1.761 1.773 1.785 1.797 1.808 1.818 1.825 1.827. — _1'318-^8} -8 25.196 -7.017 -12.032 -15.125 -13.669 -6.580 5.505 20.958 37.906.032 148.87 H= 14.760 14.468 14.004 13.446 13.384 13.558 13.710 13.843 13.958 14.044 -14.080 V= 54.790 44.065 24.657.401 -9.090 -9.899 -8.809 -5.621 -.420 6.379 14.028 D= 1.8.60 1.841 1.812 1.779 1.775 1.785 1.794 1.802 1.809 1.814 1.817 Q= 148.865 117.297 63.581.998 -22.500 -24.782 -22.274 -14.340 -1.080 16.493 36.361.040 157.51 H= 14.931 14.772 14.517 14.034 13.568 13.550 13.654 13.738 13.806 13.851 13.882 V= 57.272 48.927 34.603 12.226 -8.132 -12.697 -11.097 -1.590 -2.077 5.225 13.52 7 D= 1.871 1.861 1.844 1.814 1.786 1.785 1.791 1.796 1.800 1.803 1.805.............................. -5-5t-. 8-._-5 —=3;- _ _- 45- 92.424 31.591 -20.373 -31.771 -27.957 -19.228 -5.285 13.341 34.598.048 160.00 H= 15.661 -' 15.001 14.877 14.570 14.009 13.610 13.587 13.624 13.646 13.661 13.674 V_ 57.605 50.660 39.903 23.704 1.939 -12.017 -12.362 -8.488 -2.770 4.625 12.980 )U- 1.881 1.876 1.868 1.848 1.812 1.788 1.787 1.789 1.791 1.791 1.792 Q= 160.000 140.041 109.327 63.548 5.002 -30.188 -31.006 -21.343 -6.975 11.651 32.746.056 158.01 H= 15.182 15.179 15.127 14.941 14.502 13.888 13.554 13.503 13.484 13.468 13.471 V. 56.44- 50.487 42.366 31.446 15.210 -3.545 -11.313 -8.421 -2.595 4.545 12.426 0= 1.888 1.888 1.885 1.872 1.843 1.805 1.785 1.782 1.781 1.180 1.1780 Q= 158.010 141.377 118.190 86.553 40.583 -9.070 -28.315 -21.007 -6.465 11.313 30.932.-_-064 1 —52.-86~ H= 15.'28~0 15.:316- 1-5=,299 i5-176 4.-871- 14..3. 02 13.'688 — 13'00-. 135 —-i32.2-9-.2-. V= 54.184 49.269 43.332 36.230 26.047 10.333 -4.154 -6.794 -1.773 4.855 11.898 0= 1.895 1.898 1-897-1.'888 1.86"7 1..830. 1].793 1.176.1.9772. 1.170 1.70 o= 152.860 139.357 122.410 101.431 71.330 27.191 -10.489 -16.835 -4.371 11.945 29.262.072 145.57 H= 15.358 15.415 15.411 15.318 15.095 14.662 14.000 13.427 13.190 13.129 13.119 V= 51.303 47.573 43.567 39.271 33.548 23.801 9.307 -.746 -.129 5.3721 11.422 U= 1.901 1.905 1.904 1.898 1.882 1.854 1.812 1.778 1.'64 1. 7b1 1.760............= 145.567 1 —-.-35.552 124, 100 111.094 93.368 64.222 23.994 -1.852 -.316 13.061 21,.799.080 136.91 H= 15.412 15.478. 5476 15.395 - 15.215 14.889 14.343 1 3.650 13.162 12.999 12.960 V= 48.061 45.716 43.475 41.353 38.732 33.984 24.347 11.306 4.416 6.193 11.009 0= 1.904 1.909 1.909 1.903 1.891 1.868 1.833 1.791 1.763 1.754 1.753 9= 136.911 130.871 124.443 117.655 108.741 93.185 64.250 28.476 10.776 14.957 26.558.088 127.48 H= 15.436 15.505 15.502 15.425 15.264 15.001 14.588 13.960 13.332 12.958 12.813 V= 44.671 43.833 43.242 42.915 42.536 41.148 36.540 26.253 14.109 8.925 10.666 0= 1.906 1.911 1. 9,11 1.905 1.894 1.876 1.849 1.811 1.772 1.751 1_.747 Q= 127.482 125.737 124.021 122.370 119.858 113.744 98.084 67.589 34.807 21.501 25.607.096 117.72 -H ~i5 428 1i5-. 498 15.494 15.416 15.266 15.040 14.722 14.266 13.656 13.096~-1-2.837 V= 41.273 41.964 42.935 44.175 45.489 46.240 44.981 39.174 27.549 15.984 10.686 D= 1.906 1.911 1.910 1.905 1.894 1.879 1.857 1.828 1.791 1.759 1.747 Q= 117.721 120.312 123.059 125.883 128.195 128.181 121.886 102.836 69.414 38.844 25.607 Figure 8 (Cont'd)

-25-.104 95.15 H= 15.257 15.458 15.453 15.376 15.235 15.035 14.781 14.449 13.986 13.431 13.123 V= 33..'86.40.119 42.558 45.195 47.817 49.912 50.529 48.052 40.219 26.690 11.434 D= 1.894 1.908 1.907 1.902 1.892 1.878 1.861 1.840.1.8.1.L.L78_..... t= 95.153 114.678 121.604 128.410 134.445 138.312 137.489 127.740 103.589 66.269 27.834.112 60.01 H= 14.900 15.314 15.382 15.309 15.179 15.005 14.797 14.549 14.236 13.881 13.667 V= 21.867 36.025 42.089 45.947 49.574 52.501 54.100 53.496 48.837 35.952 12.961 o= 1.869 1.898 1.902 1.897 1.888 1.876 1.862 1.846 1.826 1.805 1.792 0= 60.010 101.884 119.630 129.900 138.821 145.168 147.373 143.201 127.934 91.951 32.681.120 25.25 H= 14.485 14.971 15.243 15.220 15.107 14.959 14.790 14.604 14.434 14.356 14.335 V= 9.476 26.831 40.230 46.356 50.728 54.168 56.279 56.477 52.54U 39.295 14.648 D= 1.842 1.874 1.893 1.891 1.883 1.873 1.862 1.850 1.839 1.834 1.833 Q= 25.252 74.006 113.185 130.193 141.301 149.286 153.232 151.768 139.530 103.791 38.632.128 -9.01 H= 14.008 14.549 14.930 15.089 15.021 14.903 14.771 14.656 14.648 14.801 14.906 V= -3.494 16.400 33.612 45.628 51.233 54.999 57.331 5/.222 51.515 36.964 15.927 0= 1.812 1.846 1.871 1.882 1.877 1.869 1.861 1.853 1.853 1.863 1.8/0 Q= -9.013 43.901 92.440 126.930 141.825 150.970 155.894 154.335 138.860 100.722 43.728.136 -.00 H= 14.045 14.067 14.521 14.821 14.912 14.841 14.757 14.745 14.903 15.173 15.301 V= -.000 5.165 24.715 40.711 50.651 55.024 57.123 55.345 46.689 32.113 16.7Z1 D= 1.814 1.816 1.844 1.864 1.870 1.865 1.860 1.859 1.869 1.888 1.897 Q= -.000 13.375 66.032 111.092 139.118 150.361 155.169 150.211 128.151 89.887 47.244.144.00 H= 13.940 13.900 14.053 14.450 14.700 14.769 14.769 14.891 15.1/3 15.450 15.548 V=.000 4.109 14.555 32.732 46.535 53.872 55.139 50.627 39.861 27.417 11.182 0= 1.808 1.806 1.815 1.840 1.856 1.861 1.861 1.869 1.888 1.907 1.914 9=.000 10.522 37.658 87.019 125.900 146.466 149.910 138.848 111.567 78.325 49.446.152.00 H= 13.815 13.799 13.774 14.020 14.388 14.641 14.825 15.083 15.409 15.634 15.698 V=.00o 5.199 10.237 23.014 38.955 49.252 50.696 43.687 32.988 23.8/5 17.448 0= 1.801 1.800 1.798 1.813 1.836 1.852 1.864 1.882 1.904 1.920 1.925'=.000 13.223 25.995 59.410 103.121 132.698 138.379 121.484 93.953 69.158 50.787.160.0-.= 1 - 3. 674 13.654 13.698 14.026 14.456 14.881 15.277 15.582 15.747 15.7u8 V=.000 5.567 10.400 16.326 28.922 40.102 42.365 35.793 27.232 21.403 11.602 D= 1.792 1.792 1.791 1.794 1.813 1.840 1.868 1.895 1.917 1.929 1.932 Q=.000 14.044 26.201 41.250 74.691 106.661 116.103 100.950'78.5..6.530 6.53 51.587.360.00 H= 13.470 13.507 13.601 13.709 13.7'94 13.840 13.857 13.8/4 13.913 13.960 13.9(1 V=.000 7.545 15.155 22.433 28.587 32.857 34.743 33.861 29.862 22.823 13.754 D= 1.780 1.782 1.788 1.794 1.799 1.802 1.803 1.804 1.806 1.8u9 1.810 3=.000 18.826 38.047 56.723 72.690 83.803 88.716 86.562 76.536 58.680 35.389.368.00 H=-13.308 13.338 13.421 13.535 13.652 13.753 13.841 13.934 14.042 14.138 14.169 V=.000 6.033 12.558 19.440 25.834 30.596 32.738 31.724 27.6-8 21.464 14.248 D= 1.771 1.773 1.777 1.784 1.791 1.797 1.802 1.808 1.814 1.820 1.822 Q=.000 14.890 31.161 48.596 65.077 77.587 83.503 81.422 71.557 55.855 37.15 3.376.00 H= 13.177 13.199 13.264 13.370 13.506 13.659 13.821 13.991 14.154 14.273 14.310 V=.000 4.891 10.333 16.378 22.351 27.027 29.159 28.147 24.521 19.658 14.589 0= 1.764 1.765 1.768 1.775 1.782 1.791 1.801 1.811 1.821 1.829 1.831 Q=.o000 11.965 25.382 40.506 55.766 68.111 74.276 72.518 63.891 51.628 38.414.384.00 H= 13.069 13.084 13.133 13.227 13.374 13.570 13.799 14.035 14.235 14.362 4.398 V=.000 4.094 8.552 13.465 18.385 22.302 24.108 23.387 20.858 17.7/0 14.797.= 1.758 1.758 1.761 1.766 1.775 1.786 1.800 1.814 1.826 1.834 1.837 O=.000 9.942 20.830 32.997 45.481 55.878 61.323 60.435 54.637 46.953 39.197.392.00 H= 12.975 12.987 13.029 13.118 13.275 13.502 13.778 14.054 14.273 14.402 14.439 V=.000 3.503 7.075 10.702 14.114 16.729 17.988 17.872 16.968 15.901.14.892 0= 1.752 1.753 1.755 1.760 1.769 1.782 1.798 1.815 1.829 1.837 1. 639.=.-........000 8.453 17.120 26.046 34.692 41.729 45.691 46.244 44.564 42.133 39.561.000.....00.H=.12.'894 12.907 12.952 13.049 13.217 13.463 13.757 14.043 14.265 14.396 14.436 V=.000 2.939 5.644 7.919 9.622 10.726 11.432 12.121 13.012 13.967 14.8o6.D.=. 1,-748 - 1..749 1.751 1.756 1.766 1.780 1.797 1.814 1.828 1.83t 1.839 9=.000 7.058 13.592 19.187 23.565 26.686 28.997 31.338 34.155 31.046 39.536 Figure 8. Computer Program and Sample of Calculations for Pulse Flow Into a Flexible Tube. Pulse Time 0.4 sec. Unstressed Tube i.d. =1.3 cm. t' - 0.1 cm, Y - 5. x 10~ dynes/cm2. Problem was started as a steady state flow with flow of 35.65 cc/sec. and head of 14. cm Hg at proximal end. Friction factor F - 0.3 and losses varying as square of the velocity. The calculations are printed out for the third pulse after steady flow, with values given for each 0.008 sec. time interval.

-26f = 0.3 was used and a time increment of At = 0O004 sec. taken with twenty equal reaches of Ax 2.5 cm. The pressure H, velocity V, diameter D and flow rate Q were calculated, and values printed out for every.008 sec. time interval and at 5 cm distances along the tube. The terminal bed pressure was taken as 100 mm Hg, and the problem was iniated by first setting up a steady flow equal to the average flow (QAVE). The sequence of main calculations in the program are as follows (see Figure 8): 1. Calculation of steady state problem, storing of H, V, D, and Q for the 21 sections 2.5 cm apart for time t = 0. 2o Calculation of interior points, statements 45 through 54, for the next time increment. 3. Calculation of the proximal boundary condition, statements 55 through 67. 4o Calculation of the distal boundary condition equations, statements 68 through 780 5o Print out of every second set of results, incrementation of T and U, and check on determination of end of solution. The program was run through 3 pulse cycles. The end of the second and third pulses showed rather close agreement, indicating that steady oscillatory flow has almost been established. A gross check on continuity was made as a measure of the accuracy of the finite difference

-27method. The total inflow in the 3rd pulse is 35.65 x 0.4 = 14.26 cc. The total outflow is 14.34 cc, and the volume of fluid within the tube has decreased by 0.29 cc. This indicates an error in continuity of about 1.5%, which could be further reduced by taking shorter reaches and smaller time increments, or by going to a method of calculation(8) having second order accuracy. B. Equations for Tapered Distensible Vessel with Distributed Outflow Since the vascular system is so complex, and several computer statements are generally required for any boundary condition, a special program has been developed for flow through a vessel having its unstressed diameter varying with length along the tube, Figure 9, and with flow leaving the vessel in a distributed manner. Such a vessel could represent the aorta, with branches of various sizes along its length. The outflow along the vessel is set up as a flow per unit length of vessel, with rate proportional to the head difference inside and outside the tube. If q represents the outflow per unit length, then q = ki (H - HB) (50) in which ki has an assigned value for each of the N sections of the vessel. The continuity equation, Equation (14), now has an extra term - (pAV)xdx - pqdx = (pAdx)t (51) which simplifies to Ax At q 5 V + V + A + 0 A A A

-28By assuming that the fluid leaving through the walls has its axial momentum reduced by contact with the branch, the momentum equation becomes - PxAdx - To-rDdx = (pAV2)x dx + (pAVdx)t + pqVdx After reducing this equation in a manner similar to Equations (18) and (19), one obtains Equation (19) as before. In order to take the tapering effect into account Ax is evaluated as follows: A = a0 with A a function of x. Then a 2 A = A o - 2A ~ a Ox a2 Oa3 x But a2 = a02 -gH and 2a ax = - gHx Hence A_ A + gH = + gH (53) A Ao a2 Ao a2 in which a is the rate of change of unstressed vessel area per unit length. At/A is obtained as beforeo Following the previous procedures in developing the finite difference equations, the two controlling equations become (the interpolation equations VR, HR, VS, HS are unchanged):

-29Vp = V - g (H - HR)l - V t - (54) VP R a )A 2D t' Aoao2 Vp VS C + g (Hp HS) FVc IVcAt fVC / VC I At a ap3g3t P = Vs + -- (Hp HS ) -VC'At +?_ (55) aC 2D Ao A0 By addition, and then by subtraction, the two equations yield values of Vp and Hp at interior points: V V fV IVs At Vp = + VS + (H - HS) - C A (56) 2 2aC R 2D Ha + HS ac + 2agcVAt 2ac qgt (57) p + 2g - Aa 2 0 The boundary conditions are handled exactly as before, except that either Equation (59) or Equation (55) is used, depending upon whether it is a right-end or a left-end boundary, respectively. For actual computation Equation (50) is used to eliminate q from the working equations. As very small time increments are generally used (.001 sec.), it is a satisfactory approximation to allow H to equal HC in Equation (50), although this is not absolutely necessary. With this type of distributed outflow (Equation 50), each branch is treated as if it were a terminal bed, i.e., a definite relation between pressure H and flow q is established for each branch.

-30Ki ( H - HB) I I I V A I IH, Figure 9. Tapered Vessel with Distributed Outflow Along the Walls. Flow Per Unit Length Through i-th Reach Shown. EB is Terminal Bed Pressure. 17. 17.'\ _____._i PROXIMAL / \ __ _ - DISTAL 16.\ i{~ \ \ ~ Figure 10. Pressure-Time Graphs for Two Cross 1l\(P~~~~~ \ ~Sections of the Femoral Artery of a Dog. /. \ xDistance Between Sections 9.4 cm, Proximal'9 \ \ Diameter 0.316 cm, Distal Diameter 0.228 cm. 0 _1____U5. _ \ \ Branch Arteries Between the Two Sections are Ligated. I \ H I/ \ Cm 14. Hg I 12. II / ___________________ _________ _________ 0.1.2.3.4 TIME, SEC.

III. COMPARISON OF MEASURED FLOW WITH, FLOW COMPUTED FROM PRESSURE-TIME DATA IN A TAPERING VESSEL To compare a calculated flow with experimentally measured flow an in vivo experiment was performed. Pressure and flow data were obtained from the femoral artery of a dog. Pressure was determined by inserting short needles at two points in the artery 9.4 cm apart. Small branches between were ligated at the vessel wall. Suitable catheters and fittings were connected between the needles and two Sanborn differential transducers. Flow was measured by means of a non-cannulating electromagnetic flow probe placed immediately upstream from the proximal pressure needle. Data was recorded on a four-channel Sanborn (350) system. Figure 10 shows the proximal and distal pressure obtained for one cycle. The solid line in Figure 11 is the flow obtained for the same cycle. The dashed line in Figure 11 is the flow calculated by using the experimental pressure data as boundary conditions. In computing the flow from the pressure-time data values of pressure were read off the strip chart for each OoG1 sec. By parabolic interpolation these values were replaced by 0.001 seco interval values over the pulse cycle. These computed values then provided the boundary conditions needed to compute flow during the cycle. The diameter of artery was assumed to decrease linearly along the length under consideration, and the length was split into seven equal reaches. Initially the flow was assumed to be steady at about the average flow, with the head (pressure) constant along the artery. After two cycles have been computed the flow has almost achieved its steady oscillatory charactero -531-.

3. i / ---- \ ELECTROMAGNETIC FLOWMETER 1/ \\ 2 L \\ --— CALCULATED: 2. \\ /I \ I LLI 0 0.1 Q2 0.3 0.4 Q5 TIME, SEC Figure 11. Electromagnetic Flow Meter Data For the Upstream Section of the Femoral Artery (Figure 10) and computed Flow Using the Pressure-Time Data of Figure 10 for Boundary Conditions. In the Computer Solution the Speed of Pulse Wave at H = 11.5 cm Hg was a = 1200 cm/sec. The Friction Factor was f = 0.4 and the Energy Dissipation was Taken to Vary as the Square of the Velocity.

-33The calculated flow does not coincide exactly with the experimental: several factors contribute to this discrepancy. Pressure measurements are extremely critical; errors of 1 mm in locating the zero pressure datum line can result in a computation that reverses the direction of average flow. Physical dimensions of the transducer fittings and lines may well introduce their own values and transients that alter the recorded pressure pulse. Dampening of the recorded flow from amplifier filtering factors in the flowmeter must also be considered. Several of these factors have been made an object of special study using programming and the theoretical model. Results of these studies, though not a subject of this paper, have shown that more refined raw data will produce better correspondence between calculated and experimental results. In making the calculations the frictional effects are quite unknown, but very important. By a gradient method, values of f and n in the friction term f CIVC n-l 2Dc were obtained that gave the best fit (least square method) with the flowmeter data. These results yielded values of f of about 0.4 and n about 2.0.

SUMMARY The basic differential equations for elastic wall material and for continuity and momentum are derived, including fluid frictional resistance of the wall of the tubes, based on one-dimensional flow. These partial differential equations are transformed into four ordinary differential equations using the theory of characteristicso Then difference equations are developed and by an interpolation method (method of specified time intervals) equations are obtained for computation of velocity and pressure at equally-spaced sections along the vessel at specified equal time intervals. The equations are first applied to a flexible tube of initial constant diameter, with a pulse flow taken from in vivo experimentso Equations are then developed for tapering tubes with distributed outflow along their lengths (to simulate branches)o Pressure-time data from femoral artery measurements are then used to compute flow through. the artery, and the results are compared with electromagnetic flowmeter datae 534

ACKNOWLEDGEMENTS The authors are indebted to the University of Michigan Computing Center for use of its digital computer, and to the NIH for its support through project No. HEO7443-01 on "Pulsalite Flow Through Distensible Tubes". The experimental work was completed while Dr. Keitzer was supported in part by projects U.S.P.H.S.B-2290C1 and U.S.P.H.S.H-4260C3 -35

NOTATION a = speed of pressure pulse through vessel a = speed of pressure pulse through vessel at zero pressure A = area of vessel cross section, point on x-t plot A = area of vessel cross section at zero pressure o B = point on x-t plot C = designation of characteristic curve, point on x-t plot D = diameter of vessel D = diameter of vessel at pressure zero f = Darcy-Weisbach friction factor F = force on fluid element g = acceleration due to gravity H = pressure expressed in length of column of fluid flowing K = constant L = label for a differential equation m = exponent P = pressure, or point to be computed q = outflow per unit length Q = flow through vessel R known point on characteristic curve R Reynolds number S = tensile stress in vessel wall; known point on characteristic curve t = -time -36

-57t' = thickness of vessel wall t' = thickness of vessel wall at pressure zero 0o T = tensile force per unit length in tube wall V - velocity in vessel X = distance along vessel Y = elastic modulus of vessel wall = rate of change of unstressed vessel area per unit length @ = ratio At/Ax X = undetermined multiplier = dynamic viscosity p = density of fluid T = frictional wall shear stress o

NOTATION FOR MAD PROGRAM A speed of pressure pulse at pressure HO AC = speed of pressure pulse at pressure H AO = speed of pressure pulse at zero pressure ACM = speed of pressure pulse at pressure H(N) C1, C2, C3, C4, C5, C6, C7 = constants D = diameter DELT = time increment DHF = steady state pressure drop in reach Ax F = friction factor G = gravity H = pressure, from previous calculation HO = bed pressure HB initial pressure HP = pressure to be calculated HR - interpolated head HS = interpolated head I = integer K = constant in terminal bed relation KK = constant N = number of reaches P constant PMT, PPT, PTT = constants in computing pulse PTIME = pulse period Q = flow -58

-59QAVE = steady-state flow QM = maximum pulse flow IHO = density SG = specific gravity of mercury in terms of fluid flowing T = time TM = constant in computing pulse U = integer V = velocity, from previous calculation VP = velocity to be calculated VR interpolated velocity VS = interpolated velocity

REFERENCES 1. Womersley, J. R. "An Elastic Tube Theory of Pulse Transmission and Oscillatory Flow in Mammalian Arteries." WADC Technical Report Tr 56-614, January 1957. 2. McDonald, D. A. "Blood Flow Through Arteries," The Williams and Wilkins Co., Baltimore, 1960. 3. Warner, Homer R. "Use of Analogue Computers in the Study of Control Mechanisms in the Circulation" Federation Proceedings Vol. 21, No. 1, January - February, 1962 Symposium on Computer Techniques in Physiology. 4. Remington, J. W. "Contour Changes of the Aortic Pulse During Propagation" Am. J. Phys. Vol. 199, No. 2, August 1960, PP. 331-334. 5. Courant, R., and Friedrichs, K. 0. "Supersonic Flow and Shock Waves," Interscience Publishers, New York, 1948. 6. Stoker, J. J. "The Formation of Breakers and Bores," Communications on Applied Mathematics, Vol. 1. No. 1, 1948. 7. Lai, Chintu, "A Study of Waterhammer Including Effect of Hydraulic Losses." Doctoral Dissertation, The University of Michigan, 1961. 8. Streeter, V. L., and Lai, Chintu. "Water-Hammer Analysis Including Fluid Friction." Journal of the Hydraulics Division, Proc. ASCE, Paper 3135, HY3, pp. 79-112, May 1962. 9. Streeter, V. L. "Waterhammer Analysis with Nonlinear Frictional Resistance." Conference on Hydraulics and Fluid Mechanics, University of Western Australia (In press, presented December 1962). 10. Peterson, L. H., Jensen, R. E., and Parnell, John. "Mechanical Properties of Arteries in Vivo" Circulation Research, Vol. VIII, No. 3, pp. 622-639, 1960. 11. Streeter, V. L. "Fluid Mechanics," 3rd ed. McGraw-Hill Book Co. Inc., New York, 1962, pp. 210-222. 12. Lister, M. "The Numerical Solution of Hyperbolic Partial Differential Equations by the Methods of Characteristics." in Mathematical Methods for Digital Computers, John Wiley and Sons, New York, 1960. -40