## Pulsatile Flow in Straight Cylindrical Tubes

The video demonstrates oscillatory axial flow velocities in a circular tube with a gradual transition from low \(\alpha\) to high. DOWNLOAD an avi of the video.

#### Where does it come from?

This article is largely a recapitulation of a former one on Poiseuille Flow. We are going to see exactly how and why fluid velocities do as they do in pulsatile flow. Unfortunately this is a purely mathematical/analytical endeavor; there is virtually no way anyone would figure this out experimentally. As in Poiseuille flow, the relationship comes not from experimental data but from a story problem or mathematical derivation. The solution is derived from 3 basic assumptions: 1) Newton's fundamental law of motion, i.e. \(\textbf{F} = m \textbf{a}\) ( force = mass times acceleration); 2) a geometric assumption about the nature of the flow, i.e. a circular cross-sectional geometry in which the fluid elements flow straight down the tube and that the flow is fully developed; and 3) an assumed relationship between relative fluid motion and shear stress (a Newtonian fluid is assumed where shear stress is proportional to *velocity gradient* through the Newtonian viscosity). While an understanding of these relationships is of great value, let us realize from the start that there is no location in the circulation where the last 2 assumptions are met. However, the first assumption ( Newton's laws of motion ) is more reliable in this vicinity of the universe than any medical or physiological "fact" (and that includes the one about death).

The derivation of the of pulsatile flow relationships is identical to the earlier Poiseuille solution except for the fact that we allow the fluid to be accelerated by a sinusoidally varying pressure gradient. Consequently the derivation of the solution is identical to the Poiseuille until about 3/4 of the way down the page. The derivation requires ruthless logic, scrupulous mathematical accounting of forces acting on fluid elements, and courage that we will be able to find our way to the result.

If you'd rather learn about this topic more without equations, you can skip to the next article which shows the results of our mathematical labors graphically.

#### Getting the facts on to "paper"

To this end, consider a straight circular tube of radius \( r_0\) and an annulus of fluid inside the tube. The annulus consists of a thin rim of fluid with mean radius \( r\) and thickness \( \Delta r\) in the radial direction and length \( \Delta z\) in the tube's axial direction, \( z\) .

We begin by actually assuming how the flow will occur, i.e. by moving in a straight line in the axial direction only. We'll also be assuming shortly that the forces on the fluid occur in an oscillatory manner, i.e. Sinusoidally. Newton's laws of motion tell us that force equals mass multiplied by acceleration, and that is the equation we are going to derive and solve.

Our next step is to enumerate and account for the forces acting on the fluid annulus (each and every one). Here we use a simple, but powerful conceptual tool called a *free body diagram* where the free body is the annulus of fluid. The free body diagram allows us to isolate an arbitrary part of the fluid and define rigorously the forces acting on the fluid. This process is equally valid for solid structures.

There are two basic kinds of forces when discussing mechanics of materials, fluid or solid: 1) Surface forces occur where the surface of the free body touches its surroundings, other fluid elements in this case. 2) Body forces are forces that seem to act from afar, to reach across space to affect an object. The most physiologically relevant of these in cardiovascular applications is gravity (until we start accelerating people in various kinds of transportation gadgetry). Even so, gravity is not very important to the discussion and cardiologists typically go out of their way to leave gravity out of the equation in a literal sense. Consequently, we only have surface forces in this problem.

\(\Large p A = p|_{\small{z-\Delta z/2}} 2 \pi r \Delta r\)

In the following figures, the surface under discussion is highlighted in color. The first surface force to consider is the pressure acting on the upstream surface of the annulus. Pressure is an example of a stress which has physical units of force/area. We calculate the value of this force by multiplying the pressure by the surface area it acts upon. Because of the way the problem is formulated, pressure has a different value at each axial location in the tube ( each value of \( z\) ). For the annulus in question, we are assuming its position at arbitrary axial location \( z\) with axial thickness \( \Delta z\) , so that the upstream surface is at axial coordinate \( z -\Delta z/2\) . Using common notation, we designate the pressure at this location as \( p|_{\small{z-\Delta z/2}}\) .

The surface area that this pressure acts upon is equal to a difference in areas between two circles; the outer one has radius \( r + \Delta r/2\) and the inner one has radius \( r - \Delta r/2\) . This comes out to \( 2 \pi r \Delta r\) . Without doing the math, you can think of this circular rim of surface as a rectangle that is been warped into a circle; the width of the rectangle is \( \Delta r\) and the length is \( 2 \pi r\) , the perimeter length of the circle.

The above pressure force is opposed by a similar one acting on the downstream rim of the annulus. The area of action is identical to the last calculation, but the pressure is different, occurring at a location \( \Delta z\) downstream of the previous pressure determination and represented by the notation \( p|_{\small{z+\Delta z/2}}\) .

\(\Large p A = p|_{\small{z+\Delta z/2}} 2 \pi r \Delta r\)

The first above pressure acts in the \( z\) direction whereas the second acts in the \( -z\) direction; these forces will have the opposite sign in the final accounting of forces.

Now we turn to the determination of shear stresses acting on the annulus. These are designated as \( \tau\) in general but we are specifically interested in shear stresses that translate to forces in the axial direction. Mechanical stress is a second order tensor. Vectors are first order tensors and have 3 components in a three-dimensional space, e.g. the \( x\) , \( y\) , and \( z\) directions. A second-order tensor has 9 components. The notation \( \tau_{rz}\) can be read as the "shear stress on the \( r\) surface in the \( z\) direction"; you can readily imagine how we would obtain 9 components in a three-dimensional space by having a stress on each of the 3 face orientations and in each of the 3 directions.

For the inner surface of the annulus, we must form an expression for the force \( \tau_{rz}A\) , i.e. the shear stress multiplied by the area it acts upon. We can designate the shear stress at the location \( (r-\Delta r/2)\) as \( \tau_{rz}|_{\small{r-\Delta r/2}}\) in keeping with previous notation. For a Newtonian fluid, shear stress is equal to the viscosity multiplied by the appropriate *velocity gradient* and this translates to \( \mu \Large\frac{\partial u}{\partial r}|_{\small{r-\Delta r/2}}\) for the problem at hand. We are using \( u\) to represent axial velocity and \( \Large\frac{\partial u}{\partial r}\) is the velocity gradient, mathematically the derivative of \( u\) with respect to the radial coordinate, \( r\) . \( \mu\) is the Newtonian viscosity.

\( \Large \tau_{rz}A = \tau_{rz}|_{\small{r-\Delta r/2}} 2 \pi r|_{\small{r-\Delta r/2}} \Delta z=\mu \Large\frac{\partial u}{\partial r}|_{\small{r-\Delta r/2}} 2 \pi r|_{\small{r-\Delta r/2}} \Delta z \)

The surface area that the shear stress acts upon at this location is slightly tricky to express because the area changes with radius, just as the stress itself does. Nevertheless, the area expression is simply that of a rectangle with thickness \( \Delta z\) and length \( 2 \pi (r - \Delta r/2)\) so that the area expression is \( A|_{r-\Delta r/2} = 2 \pi \Delta z (r - \Delta r/2)\) or \( 2 \pi \Delta z r|_{\small{r - \Delta r/2}}\) . The entire expression for the force is shown adjacent to the figure above. The shear stress acting on the outer rim of the annulus is derived similarly (below).

\(\Large \tau_{rz}A = \tau_{rz}|_{\small{r+\Delta r/2}} 2 \pi r|_{\small{r+\Delta r/2}} \Delta z =\mu \Large\frac{\partial u}{\partial r}|_{\small{r+\Delta r/2}} 2 \pi r|_{\small{r+\Delta r/2}} \Delta z \)

In case it has not yet occurred to you, it's (utterly) important to recognize that equations here are to express *physical* concepts. "Equations" used in cardiology are sometimes the result of statistical regressions, or they may express a concept or approximation in a general way. We can admit no conceptual or mathematical slack for the endevour at hand! We're deriving an equation whose physical units (thus far) are force. Pressure must be expressed in terms of force/area, e.g. dyne/cm^{2} (a dyne is a gram-cm/sec^{2} which is a mass multiplied by an acceleration). If you can only think of pressure in terms of mmHg, you've been spending too much time in the cath lab. When we multiply all of these things together, they have to have comparable physical units or we've got nothing but muck! The Newtonian viscosity, for example, is a physical quantity that multiplies the velocity gradient and yields a shear stress. If your shear stress has units of dyne/cm^{2} and your velocity gradient is in cm/sec/cm (i.e. sec^{-1}), then your viscosity had better have units of dyne-sec/cm^{2} (or gram-cm^{-1}-sec^{-1}) so that the physical units are correct! It just so happens that the physical unit of viscosity just named is called the Poise after the gentleman under discussion.

#### Putting it all together

We are now in a position to begin combining (mathematical ) terms into a coherent expression (have courage!). THIS IS WHERE THE SOLUTION DEPARTS from the Poiseuille situation. As noted previously, the upstream and downstream pressure forces oppose each other:

\( \Large 2\pi r\Delta r\left[p|_{\small{z-\Delta z/2}}-p|_{\small{z+\Delta z/2}}\right]\)

Similarly, shear stresses at the inner and outer surfaces of the annulus oppose each other:

\( \Large 2\pi \Delta z \mu\left[(r \Large\frac{\partial u}{\partial r})|_{\small{r+\Delta r/2}}-(r \Large\frac{\partial u}{\partial r})|_{\small{r-\Delta r/2}}\right]\)

The net sum of all the forces acting on the annulus is NOT equal to zero, but equal to the mass of the annulus multiplied by its rate of * acceleration*. The mass of the annulus is its volume, \(2\pi r\Delta r \Delta z\), multiplied by the fluid density \(\rho\). The acceleration is \(\partial u/\partial t\).

\( \Large 2\pi r\Delta r\left[p|_{\small{z-\Delta z/2}}-p|_{\small{z+\Delta z/2}}\right]+2\pi \Delta z \mu\left[(r \Large\frac{\partial u}{\partial r})|_{\small{r+\Delta r/2}}-(r\Large\frac{\partial u}{\partial r})|_{\small{r-\Delta r/2}}\right]=\rho 2\pi r\Delta r \Delta z \frac{\partial u}{\partial t}\)

We divide the entire equation by \( 2\pi \Delta r \Delta z\) and simplify:

\( \Large r \frac{p|_{\small{z-\Delta z/2}}-p|_{\small{z+\Delta z/2}}}{\Delta z}+\mu \Large\frac{(r \Large\frac{\partial u}{\partial r})|_{\small{r+\Delta r/2}}-(r \Large\frac{\partial u}{\partial r})|_{\small{r-\Delta r/2}}}{\Delta r}=\rho \frac{\partial u}{\partial t}\)

(The equation just acquired new physical units : force/cm2) The next step is allow \( \Delta r\) and \( \Delta z\) approach zero and to recognize the definition of the derivative as defined in calculus:

\( \Large\frac{df}{dx} \equiv \Large\frac{f|_{x+\Delta x} - f|_{x}}{\Delta x}_{lim\Delta x \rightarrow 0}\)

This allows us to express the work thus far as a differential equation:

\( -r\Large\frac{\partial p}{\partial z}+\mu \Large\frac{\partial }{\partial r} (r \Large\frac{\partial u}{\partial r})=\rho \frac{\partial u}{\partial t}\)

or

\( -\Large\frac{\partial p}{\partial z}+\Large\frac{\mu}{r}\Large\frac{\partial}{\partial r}(r \Large\frac{\partial u}{\partial r})=\rho \frac{\partial u}{\partial t}\)

This is a linear, ordinary differential equation whose solution will tell us how the velocity changes as a function of radius. It may look like there are 2 unknown functions in the equation, \( p(z)\) and \( u(r)\) , but there's really only one. \( -\partial p/\partial z\) is a given in this problem.

THE NEXT STEP IN THE SOLUTION IS TO REPLACE ALL OF THE TIME VARYING QUANTITIES IN THE EQUATION WITH SINUSOIDS. RECOGNIZING THAT THE EQUATION IS LINEAR ALLOWS US TO SEEK A PARTICULAR SOLUTION TO THE EQUATION FOR AN OSCILLATORY PRESSURE GRADIENT OF ANY FREQUENCY YOU'D LIKE TO CHOOSE. IF WE CAN DO THAT, WE CAN ADD ANY COMBINATION OF OSCILLATIONS TO OBTAIN A SOLUTION FOR ANY INPUT.

Subsequently, \(u(r)\) and \(p\) are no longer time-domain signals but repressent sinusoids of angular frequency \(\omega\). Symbols are subsequently represented by capitals, i.e. \(U(r)\) and \(P\). To be perfectly explicit, we should write \(U(r,j\omega)\) and \(P(j\omega)\) where \(\omega\) is the * angular frequency* of the oscillation; \(\omega = 2 \pi f\) where \(f\) is the frequency (Hertz) of the oscillation.

\(\Large \rho j\omega U-\Large\frac{\mu}{r}\Large\frac{\partial}{\partial r}(r \Large\frac{\partial U}{\partial r})= -\frac{\partial P}{\partial z}\)

Because we've just switched to sinusoids, the \(\partial u/\partial t\) (acceleration) term has been replaced by \(j \omega U\). In mathematical terms, we have switched from the Time Domain to the Frequency or ** Fourier Domain**. \(j\) is the imaginary constant \(\sqrt{-1}\); we found earlier that this seemingly silly number is a very useful thing to have around. It allows us to REDUCE THE PROBLEM FROM A PARTIAL DIFFERENTIAL EQUATION TO AN ORDINARY DIFFERENTIAL EQUATION WITH THE DEPENDENT VARIABLE \(U(r)\). \(-\partial P/\partial z\) is NOT a variable in this equation; it's the pressure gradient supplied for the problem of interest (known) and it's just a sinusoid. It's an oscillating pressure force that drives the motion of each fluid annulus. The equation is written equivalently as:

\(\Large \rho j\omega U-\mu(\frac{\partial^2 U}{\partial r^2}+\frac{1}{r}\frac{\partial U}{\partial r})= -\frac{\partial P}{\partial z}\)

OR

\(\Large \frac{\partial^2 U}{\partial r^2}+\frac{1}{r}\frac{\partial U}{\partial r}-\frac{\rho j\omega}{\mu} U-=\frac{1}{\mu}\frac{\partial P}{\partial z}\)

It turns out to be a good idea at this point in the proceedings to substitute the radius ** fraction** \(y=r/r_0\) in place of the radial coordinate (\(r_0\) is the tube radius). The equation could be solved in that form but it's also a good place to introduce the oscillatory flow parameter, \(\alpha = r_0 \sqrt{\rho \omega/\mu}\). The pressure gradient ( \(\partial P/\partial z\) ) is just a constant in the equation, albeit a complex one. Then the equation is:

\(\Large \frac{\partial^2 U}{\partial y^2}+\frac{1}{y}\frac{\partial U}{\partial y} -j \alpha^2 U=\frac{\partial P}{\partial z}\frac{r_0^2}{\mu}\)

Solving the Equation of Motion

The \(\overline{F}=m\overline{a}\) equation in a physics problem is often denoted the *equation of motion* and is typically a differential equation like the one above. This is an example of Bessel's equation, a differential equation that drops out of the sky frequently when working with circular geometries.

\(\Large x^2 \frac{d^2 y}{dx^2} + x \frac{dy}{dx} + (x^2 - a^2)y = 0\)

In our particular example the value of \(a\) is zero and the solution involves a Bessel Function of order 0. No doubt solving this equation could be pretty scary to some (me too), but the solution is as follows (Witzig 1914, revived and popularized by Womersley in 1955):

\( \Large U(y)= -\frac{\partial P}{\partial z} \frac{r_0^2}{j \mu \alpha^2} [1-\frac{J_0(j^{3/2}\alpha y)}{J_0(j^{3/2}\alpha)}] \)

Equivalently:

\( \Large U(y)= -\frac{\partial P}{\partial z} \frac{1}{j \rho \omega} [1-\frac{J_0(j^{3/2}\alpha y)}{J_0(j^{3/2}\alpha)}] \)

I prefer to keep \(-\partial P/\partial z\) as an intact entity, i.e. with the minus sign firmly attached. This is to continuously remind me that it is the * negative* pressure gradient that drives fluid flow. Reminding you how this comes about, an upstream pressure that is greater than the downstream ( positive force driving flow downstream) constitutes a negative pressure gradient. A shorthand I'll also introduce is that \(P_z \equiv \partial P/\partial z\) so that:

\( \Large U(y)= -P_z\frac{r_0^2}{j \mu \alpha^2} [1-\frac{J_0(j^{3/2}\alpha y)}{J_0(j^{3/2}\alpha)}] \)

and

\( \Large U(y)= -P_z \frac{1}{j \rho \omega} [1-\frac{J_0(j^{3/2}\alpha y)}{J_0(j^{3/2}\alpha)}] \)

The \(J_0\) part means the "Bessel function of order 0". While Bessel functions may not be as familiar as some others, just imagine for the moment that there's a button on your calculator next to the "Sin" or "Cos" button that will spit out "\(J_0\)" of whatever number you decide to put in; we'll see what the numbers look like on the next page. Remember that \(y\) is the radius fraction so it's going to equal 1.0 at the wall of the tube. So at the wall the whole thing in the square brackets [] is going to be 0; the velocity is 0.0 at the wall of the tube. The part in front of the brackets is just a thing that multiplies the solution (e.g. bigger negative pressure gradient, greater velocities). \(\alpha\) is the thing inside the Bessel function that affects the SHAPE of the velocity profile. IT'S THE ONLY THING that affects the shape of the velocity profile and that's why it's been singled out as a separate, specific entity.

The solution for the axial velocity \(U(y)\) above is a function, a number that varies with the radius fraction, \(y\) ( \(y\) is 0 at the tube centerline, 1 at the wall). But \(U(y)\) is a ** complex** number. It has a magnitude or

*at each value of \(y\) but also a*

**modulus****- an**

*phase*

**angle.****\(-P_z\) is a**

**complex number (magnitude and phase) that you supply to this function when you want to find the velocity. Velocity is not (in general) in phase with the pressure force (negative pressure gradient). -Pressure gradient constitutes a force that acts to accelerate the fluid, to change it's velocity. We'll see that the force due to pressure gradient increases before the velocity does.**

*single*Having these complex numbers lying around can be a little daunting, but let's see how we use these numbers to actually see what the velocities are as a function of time. You've supplied \(-P_z\) as a modulus and phase; let's call this \(|-P_z|\), the absolute value or magnitude (modulus) of \(-P_z\), and \(\angle(-P_z)\) the angle (phase) of \(-P_z\). Then the negative pressure gradient as a function of time is :

\(\Large -p_z(t)=|-P_z| \cos[\omega t + \angle(-P_z)]\)

Similarly:

\(\Large u(y,t)=|U(y)| \cos[\omega t + \angle(U(y))]\)

The small letters for \(-p_z\) and \(u(y,t)\) indicate we're back to the ** time domain** – each of these functions oscillates sinusoidally at angular frequency \(\omega\). Be aware that \(U(y)\) has a

*modulus and phase at each radial coordinate. We already saw that the modulus of the velocity is going to be zero at the wall, but will also find that the velocity is ( relatively )*

**different***with the pressure gradient force near the wall and ( relatively )*

**in phase***near the tube centerline.*

**out of****phase**The solution shown is for a particular sinusoidal angular frequency ( \(\omega\) ) or, more specifically, for a particular value of \(\alpha\). However we can obtain a solution for *any* value of \(\omega\) you like (** every** value that applies to the situation at hand). Because the equations are linear, we can add them together to produce ** any** pressure gradient ( or flow ) desired. Each frequency is independent of every other because of the linearity of the equations. Whether this linearity applies to the system of interest is a separate issue.

The quantity \(\alpha\) is worth taking a closer look at.

\( \Large \alpha = r_0 \sqrt{\frac{\rho \omega}{\mu}} \)

Since this is physical stuff, we'd better figure out the physical units of this quantity to get a handle on its meaning.

\(\mu\) is viscosity and has physical units of stress / strain rate. It's the thing that multiplies velocity gradient to get a force per unit area. Using metric system units, stress units of force/area are * gm-cm/s^{2}* (a unit of force called a dyne) divided by

*or*

**cm**^{2}**. Strain rate is**

*gm/cm-s*^{2}**(velocity) per cm; this is a velocity gradient with units**

*cm/s**. Then dividing stress by strain rate gives*

**1/s***divided by*

**gm/cm-s**^{2}*which yields*

**1/s***; that's a unit of viscosity called a*

**gm/cm-s***.*

**Poise**Angular frequency has units of radians per second and radians are dimensionless; so \(\omega\) has units of * 1/s*. Density is in

*so the bit under the square root sign resulting from \(\rho \omega /\mu\) is*

**gm/cm**^{3}*which boils down to*

**(gm/cm**^{3})(1/s)/(gm/cm-s)*. Take the square root and it's just*

**1/cm**^{2}*. Multiply that by the radius \(r_0\) (*

**1/cm***) and you get ... \(\alpha\)*

**cm**

**is a dimensionless number -- no physical units.**

Without going into the gore of the Buckingham \(\Pi\) theory I'll just state that this is what happens in physics if you do things correctly. We find that physical behavior depends on **ratios,** non-dimensional parameters like \(\alpha\), and not the absolute size, speed, or mass of things per se. \(\alpha\) simply shows up of its own accord as a result of solving the equation of motion (that Bessel equation thing). Let's look at it another way:

\(\Large \alpha^2 = r_0^2 \frac{\rho \omega}{\mu}= r_0^2 \frac{\rho 2\pi f}{\mu}\)

** \(\alpha^2\) is proportional to the oscillation frequency.** In essence it's a nondimensional

**, a frequency that's been rescaled to suit the physics of the problem. Similarly, \(1/\alpha^2\) is a nondimensional time scale. Everything in this universe has an appropriate time, length, and mass scale. (You wouldn't use a micrometer to measure the distance from here to the sun.)**

*frequency*Now we'll take another look at the solution for velocity, particularly the big scary fraction within the square brackets. We're going to supply a number, \(j^{3/2}\alpha y\) (called the argument of the function), to the Bessel function gadget, \(J_0\), and it spits back another number. Yes they're complex numbers, but don't let it get to you. The thing we are giving to the Bessel function consists of 3 * pure* ( nondimensional ) numbers multiplied together; \(j\) is the square root of -1, \(y\) is the fractional radial coordinate, and we just figured out that \(\alpha\) is nondimensional also. That turns out to be a good thing; a Bessel function actually ends up raising the \(j^{3/2}\alpha y\) argument to fractional powers. Have you ever tried to take the square root of a centimeter? While you may find something like that in a clinical journal, it's nonsensical when we're talking physics. So the whole thing in the square brackets is a pure number and is of Order (1); it varies across the radius but doesn't shoot off to infinity. So the answer has a nondimensional part ( inside the square brackets ) that describes the shape of the velocity profile, and a dimensional part (with tube radius, pressure gradient, fluid density and viscosity) that scales (sizes) the part between the [] for our particular problem (tube).

\(\alpha\) also has the interpretation of a ratio of * forces* affecting the fluid. Notice the presence of things having to mass (\(\rho\)) and oscillation (acceleration, \(\omega\) ) in the numerator, friction (viscosity, \(\mu\)) in the denominator. We don't see anything about \(\alpha\) having to do with velocity per se because the analysis ends up with velocity in both numerator and denominator; they cancel out in the end.

### Flow Rate and Longitudinal Impedance from the Velocity Profile

Our next problem is to find the flow rate stemming from the pressure gradient. To do that we need to integrate the velocity function over the cross-section of the tube:

\( \Large Q(j\omega) = 2\pi r_0^2 \int_0^1 U(y) y dy \)

\( \Large Q(j\omega) = 2\pi r_0^2 \int_0^1 -P_z\frac{r_0^2}{j \mu \alpha^2} \left[1-\frac{J_0(j^{3/2}\alpha y)}{J_0(j^{3/2}\alpha)}\right] y dy\)

\( \Large Q(j\omega) =-P_z \frac{j \pi r_0^2 }{\rho \omega}\left[1-\frac{2 J_1(j^{3/2}\alpha)}{j^{3/2}\alpha J_0(j^{3/2}\alpha)}\right]\)

The above is the oscillatory flow rate due to oscillatory pressure gradient \(P_z\). Once again it has a dimensional part that multiplies a non-dimensional shape factor within the square brackets []. Oscillatory flow increases with pressure gradient and the square of the radius; it is inversely proportional to fluid density and frequency. The part in the square brackets contains Bessel functions of the first kind (\(J\)), order 0 (\(J_0\)) and order 1 (\(J_1\)). The \(J_1\) happens in the process of integrating the function obtained previously for the velocity. Again think of this as a button on your calculator that spits out a number when you supply an input.

The stuff in the brackets is just a number, albeit a complex one. To simplify the math, Womersley called the fraction inside the square brackets \(F_{10}\), the 10 having to do with Bessel functions of order 1 and 0:

\(\Large \left[1-\frac{2 J_1(j^{3/2}\alpha)}{j^{3/2}\alpha J_0(j^{3/2}\alpha)}\right] = 1-F_{10}\)

Now the equation appears as :

\( \Large Q(j\omega) =-P_z \frac{j \pi r_0^2 }{\rho \omega}\left[1-F_{10} \right]\)

We're now in a position to write down an expression for the longitudinal impedance:

\(\Large \frac{-P_z}{Q} = \frac{\rho \omega}{j \pi r_0^2} \frac{1}{\left[1-F_{10} \right]}\)

\( \frac{-P_z}{Q}\) is the **definition** of the longitudinal impedance; the previous equation shows how we might compute it given physical characteristics of the tube and fluid. Womersley also decided to express \(1-F_{10}\) in terms of a modulus and phase:

\(\Large M'_{10} e^{j\varepsilon_{10}} =1- F_{10} = \left[1-\frac{2 J_1(j^{3/2}\alpha)}{j^{3/2}\alpha J_0(j^{3/2}\alpha)}\right]\)

\(M'_{10}\) is the modulus and \(\varepsilon\) the phas of \(1-F_{10}\). It's just a (complex) number! Here's what it looks like plotted:

And HERE is a table of values. We'll look at this again later.

### Effect of \(\alpha\) on Velocity Profile Shape

Now we'll look at a graphical representation of all that math. The next figure is a plot of the * Modulus* of the velocity (y-axis) across half the tube (x-axis) for several values of \(\alpha\). The x-coordinate is the radius fraction; it goes from 0 at the centerline of the tube (left) to 1.0 at the wall (right), regardless of the physical radius of the tube. The average velocity over the cross section (spatial average) has also been normalized to an value of 1.0 for

**frequency (value of \(\alpha\)); they all have the same peak value of flow. For \(\alpha = 0\), there is no oscillation and we have the Poiseuille solution - a parabola with a velocity that peaks out at**

*each***the average flow velocity. As has been stated previously, the velocity has to go to 0.0 at the wall.**

*twice*

As you can see, the velocity profile changes quite a bit with increasing oscillation frequency or * specifically* as \(\alpha\) changes. Changing tube radius, oscillation frequency, fluid density, or fluid viscosity, will all affect the shape of the velocity profile

*they affect \(\alpha\). The velocity profile flattens out considerably with increasing \(\alpha\); consequently fluid velocities do not vary as much over cross section. However the velocity must always go to 0 at the wall and that means there has to be a region near the wall where the velocity transitions to 0 from whatever it's doing in the middle of the tube. This means that there will be a steeper velocity gradient with increasing value of \(\alpha\). The next figure is just more of the same but for a different distribution of \(\alpha\) ( higher values ) just to show off how the velocity gradient becomes steeper and steeper with higher values of \(\alpha\).*

**because**

DOWNLOAD the Linear Hemodynamics software which is an executable file. (Software and Scenerio Files subject to change at any time and without provocation. Download the files afresh if it's not working for you.) Your computer will warn you that I am NOT a trusted source which is certainly true. The software comes up with a default model and settings that start when you **double-click the executable file (****LinearHemodynamics.exe)*** ,* but will also load

*– model and display presets that set the system up quickly to demonstrate particular aspects of hemodynamics. Individual scenario files can be downloaded as you follow through the webpages to set the software up for the topic at hand. None of the software on this website is designed to be malicious - I apologize in advance if the software causes EITHER of you any problems.*

**"Scenario" files**DOWNLOAD an executable zip file (HemodynamicScenerios.exe) containing ALL currently available Scenerio Files for this software. **Extract the files to a location of your choice by double-clicking on ****HemodynamicScenerios.exe****Double-click the ****LinearHemodynamics.exe file **** **to start the software. Click on the "SCENERIOS" tab, click the "LOAD" button, and **select the *** VelocityDemo1.xml* to set up the software to show velocities like the images displayed here.

This demonstration is to show the effect of alpha on the velocity profile, particularly in the frequency or Fourier domain. The scenario loads with several frequency components for the input, all with 1.0 average velocity over the cross section of the tube.

Look at the "Fourier Flow Velocity" Display Form which is set to show only the modulus (magnitude) of the velocity from the tube centerline (left-hand side of the plot at x = 0) to the wall (right-hand side of the plot where x = 1). Velocity always goes to 0 at the wall.

The plots show the velocity distribution (velocity "profile") across the tube for steady flow (alpha = 0) along with 4 harmonics, i.e. 4 oscillation frequencies. The zeroeth harmonic (alpha = 0) corresponds to Poiseuille flow and a parabolic velocity profile. The velocity at the center line is twice the average velocity across the tube. At higher frequencies, more specifically at higher values of alpha, the velocity profile flattens out – becomes more uniform across the tube. Nevertheless there must always be a region close to the wall where the velocity transitions to 0 at the wall. Higher values of alpha coincides with a steeper velocity gradient near the wall.

Alpha is proportional to tube radius (diameter), SQUARE ROOT of fluid density and oscillation frequency, and the RECIPROCAL SQUARE ROOT of fluid viscosity. Click on the "MODEL" tab of the Application Control Form to adjust the model parameters and observe how the Fourier Domain velocity profile changes.

Hint: "INC" and "DEC" buttons increase or decrease the parameter selected from the list by a fraction that you determine; 10% in this case. Nothing happens when you change the fraction, only when you push one of the buttons. You're not allowed to change some of the parameters, ones that are determined from other information. Diameter can be changed for example, not Radius.

Click on the "INPUTS" tab and change the base "FREQUENCY" of the input to see the effect of higher/lower oscillation frequency on the velocity profile. Change the value of \(\alpha\) through these parameters and see how the velocity profiles approach the parabolic Poiseuille profile at low \(\alpha\), or tend towards a completely flat profile at very high \(\alpha\). You'll see also on this tab that the inputs have been set for the same average velocity ("measured velocity " for every frequency so that we can compare shape of the profiles at the same flow rate. Of course you can change the magnitude of the flow velocities here, but that only affects the size of the velocity profile, not the shape.

Hint: The "RESCALE" button multiplies ALL the inputs by the factor that you select. So the button won't do anything with the number 1.0 in there. You can keep the shape of the input the same ( size relationships of all the harmonics ) but scale the magnitude up or down. There is also a "SET ALL" that sets all the harmonics to the SAME value ( both modulus and phase), whatever is shown in the numeric selection window.

Hint: It's very important that you understand what aspect of the input you are specifying. "Velocity Measured" is shown in the drop-down window here, but you are allowed to specify Flow Rate, Pressure, Pressure Gradient, OR Velocity ( the specified velocity is the * average* velocity over the cross section i.e. flow rate divided by cross-sectional area ). We haven't talked about wave reflections yet, but you're allowed to indicate whether you are specifying the antegrade, retrograde, or measured component of these quantities.

*: This can get you into trouble. If you try to specify a nonzero value for a reflected variable any system that doesn't have any reflections, you're going to crash the program ( until I get some safeguards installed).*

**WARNING**

### The Effect of \(\alpha\) on Velocity

In the last section, we kept the average velocity constant for all the different values of \(\alpha\), specifically to evaluate the effect of \(\alpha\) on velocity profile shape. However the equation shows us that velocity is also affected by pressure gradient, oscillation frequency, density, and viscosity. The last three of these also affect the velocity profile shape (through \(\alpha\)).

The next plot shows both velocity and * phase angle* for a range of frequencies, i.e. \(\alpha\)). The thing that's being held constant between the frequencies in this plot is the magnitude ( and phase ) of the

**; each flow oscillation is driven by a sinusoidal pressure gradient and all the sinusoids have the same magnitude.**

*pressure gradient*The next plot is simply for a greater range of \(\alpha\) (higher base or fundamental frequency selected for the input).

DOWNLOAD an executable zip file (HemodynamicScenerios.exe) containing ALL currently available Scenerio Files for this software. **Extract the files to a location of your choice by double-clicking on ****HemodynamicScenerios.exe****Double-click the ****LinearHemodynamics.exe file **** **to start the software. **Click on the "SCENERIOS" tab, click the "LOAD" button, and select the ****VelocityDemo2****.xml** to set up the software to show velocities like the images displayed here.

This demonstration is to show the effect of oscillation frequency of the pressure gradient (input) on the velocity profile, particularly in the frequency or Fourier domain. The scenario loads with several frequency components for the input, all with the same magnitude for the oscillatory pressure GRADIENT. Look at the "Fourier Flow Velocity" Display Form which is set to show both the modulus (magnitude) and phase of the velocity from the tube centerline (left-hand side of the plot at x = 0) to the wall (right-hand side of the plot where x = 1). Velocity always goes to 0 at the wall.

The plots show the velocity distribution (velocity "profile") across the tube for steady flow (alpha = 0) along with several harmonics, i.e. oscillation frequencies. The zeroeth harmonic (alpha = 0) corresponds to Poiseuille flow and a parabolic velocity profile. The velocity at the center line is twice the average velocity across the tube. * With higher frequency of the pressure gradient, the average velocity of flow diminishes*. This of course will coincide with a lower value for the oscillatory flow rate. In steady flow (Poiseuille), the pressure gradient is exactly balanced by friction and all of the pressure change is due to a

**loss of energy.****As soon as the flow starts to oscillate, a portion of the pressure gradient goes into**

*the fluid elements. They are being*

**accelerating**

**accelerated throughout****the oscillation cycle. Sometimes velocity is increasing and sometimes it is decreasing, but it is**

*in a sinusoidal oscillation. For a given*

**always changing***of oscillatory force (modulus of the pressure gradient), fluid velocity*

**magnitude***with increasing oscillation frequency ( greater accelerations ) and/or with greater fluid density ( greater mass ).*

**decreases**A good way to explore the physics in this situation is to change the fundamental frequency on the "INPUTS" tab. That will change the value of \(alpha\) without doing anything to the Poiseuille flow solution, i.e. the blue parabolic velocity profile will stay put and only the non-zero frequency components will change, getting smaller as you keep increasing the fundamental frequency.

Alternatively you can change the other parameters that affect \(\alpha\):

1) Rescale the pressure gradient for all harmonics by entering a number other than 1.0 in the numeric window and pushed the "RESCALE" button.

2) Change the value of \(\mu\) ( viscosity ), \(\rho\) ( density ), or diameter for parameters on the "MODEL" tab.

Each of these changes however, will cause ALL of the velocity profiles to change including the Poiseuille profile. You can have the software automatically adjust the plot range by clicking on "AutoRange" for each of the modulus and phase displays on the "DISPLAY" tab, You'll see the SCALE of the plots change as you adjust model parameters or inputs. You can remove frequencies or an entire plot using the checkboxes in the Display Tree on that tab.

If you DECREASE the fluid viscosity on the "MODEL" tab, you'll see the velocities become more disparate; the higher frequencies tend towards very low velocity. This is apparent from one of the equations for the velocity above showing that velocity is inversely proportional to viscosity and \(\alpha^2\) which turns out to be equivalent to saying it's inversely proportional to oscillation frequency.

Focus now on the * phase* of the velocity. Because the input in all cases is a sinusoidal pressure gradient ( with phase set to 0.0 ), this plot shows the phase relationship of the velocity to the force ( negative pressure gradient ) applied. For the Poiseuille flow, we would say that the negative pressure gradient is all

*with the velocity although you could argue that this doesn't make that much sense since there is no oscillation at all. It is clear that \(\alpha\) doesn't have to increase very much before we see that the phase becomes a negative. This means that the velocity starts to*

**in phase***the oscillatory pressure force; the peak of the ( negative ) pressure gradient occurs*

**lag***the velocity peak. This is of course what we expect. Whenever a net force is applied it causes an acceleration - a*

**before***in velocity. It takes a while after the pressure force is applied before the velocity actually responds so we have the peak value of the force occurring before the peak value of the velocity. Higher accelerations ( specifically greater value of \(\alpha\)) coincides with a greater lag in the velocity. Furthermore we can see that the lag is more pronounced towards the center of the tube. Near the wall of the tube, frictional forces are much more significant due to the high velocity gradient there. So a greater fraction of the pressure force there only acts against friction, NOT to accelerate the fluid. Again the Poiseuille solution is the situation where ALL of the pressure force is applied against friction (no acceleration) which is why the velocity is entirely*

**change***with the pressure force, all the way across the tube.*

**in phase**