# Oscillatory Flow in a Rigid Circular Tube

1. Introduction

Pressure-driven flow in circular tubes is a phenomenon of enormous practical significance for scientists and engineers. In elementary treatments of this topic, the pressure gradient that maintains the fluid motion is assumed to be steady, with no appreciable time-based variation. This is one of the assumptions that lead to Hagen-Poiseuille flow, a key model introduced in elementary courses on fluid mechanics and transport phenomena. While convenient, pressure-gradient constancy cannot be assumed in some situations. The goal of this post is to do away with this assumption and develop a simple model of oscillatory (i.e., unsteady) flow in circular tubes with rigid walls, as presented in M. Zamir’s excellent volume The Physics of Pulsatile Flow (2000). Two additional noteworthy references are Truskey et al. (2009) and White & Majdalani (2022).

The assumptions under which the oscillatory flow at issue takes place are essentially the same as those that govern Hagen-Poiseuille flow. The assumptions are:

1. The cross-section of the tube is circular and uniform.
2. The tube wall is rigid, i.e., it does not deform in response to fluid motion.
3. The running fluid is Newtonian, with a dynamic viscosity that is constant and independent of shear rate, temperature, or other thermomechanical variables.
4. The flow is hydrodynamically fully developed.
5. Fluid motion takes place in the axial direction only; radial and tangential velocity components are both zero.

In view of these assumptions, only the x-component (x = axial coordinate) of the Navier-Stokes equations in cylindrical coordinates survives, and we express it as

$\displaystyle \rho \frac{{\partial u}}{{\partial t}}+\frac{1}{\rho }\frac{{\partial p}}{{\partial x}}=\mu \left( {\frac{{{{\partial }^{2}}u}}{{\partial {{r}^{2}}}}+\frac{1}{r}\frac{{\partial u}}{{\partial r}}} \right)\,\,\,\,\,\left( 1 \right)$

where $\displaystyle \rho$ and $\displaystyle \mu$ are the density and dynamic viscosity of the running fluid; p is pressure, which we assume to be a function of axial position x and time t; and u is the axial velocity, which we assume to be a function of radial position r and time t. The mere assumption that velocity and pressure vary with time already places us outside the realm of Poiseuille flow, otherwise we could cancel out the first term on the left-hand side of (1) and proceed with the solution that leads to the famous paraboloid-shaped velocity distribution. Pressure p is assumed to be oscillatory in time, varying as a trigonometric sine or cosine function. As the pressure rises to its peak, the flow increases gradually in response, and as the pressure falls, the flow follows again. There is a phase difference between the pressure wave and the flow that develops in response to it: as we’ll see shortly, the greater the frequency of the pressure wave, the larger this phase difference will be.

2. Setting up the equation

A crucial feature of partial differential equation (1) is that it is linear, so that we can separate it into a steady and an oscillatory part, solve them separately, and superpose them to obtain the final solution. In the present case, the steady part of PDE (1) turns out to be identical to the solution of Hagen-Poiseuille flow, which students reading this post have already mastered. Our focus, then, is on the oscillatory part.

If the steady and oscillatory parts of the pressure and velocity are identified by subscripts “s” and “$\displaystyle \phi$,” respectively, then to isolate the oscillatory flow problem we write

$\displaystyle p\left( {x,t} \right)={{p}_{s}}\left( x \right)+{{p}_{\phi }}\left( {x,t} \right)\,\,\,\,\,\left( {2.1} \right)$

$\displaystyle u\left( {r,t} \right)={{u}_{s}}\left( r \right)+{{u}_{\phi }}\left( {\phi ,t} \right)\,\,\,\,\,\left( {2.2} \right)$

Substituting these into (1), we obtain

$\displaystyle \frac{{d{{p}_{s}}}}{{dx}}-\mu \left( {\frac{{{{d}^{2}}{{u}_{s}}}}{{d{{r}^{2}}}}+\frac{1}{r}\frac{{d{{u}_{s}}}}{{dr}}} \right)=0\,\,\,\,\,\left( {3.1} \right)$

$\displaystyle \rho \frac{{\partial {{u}_{\phi }}}}{{\partial t}}+\frac{{\partial {{p}_{\phi }}}}{{\partial x}}-\mu \left( {\frac{{{{\partial }^{2}}{{u}_{\phi }}}}{{\partial {{r}^{2}}}}+\frac{1}{r}\frac{{\partial {{u}_{\phi }}}}{{\partial r}}} \right)=0\,\,\,\,\,\left( {3.2} \right)$

As stated, equation (3.1) leads to steady flow. Equation (3.2), in turn, is the governing equation for the oscillatory part of the flow.

Because of the independence of the steady- and oscillatory-flow equations, we can also break down the pressure gradients as

$\displaystyle k\left( t \right)={{k}_{s}}+{{k}_{\phi }}\left( t \right)\,\,\,\,\,\left( 4.1 \right)$

where

$\displaystyle k\left( t \right)=\frac{{\partial p}}{{\partial x}}\,\,\,\,\,;\,\,\,\,\,{{k}_{s}}=\frac{{d{{p}_{s}}}}{{dx}}\,\,\,\,\,;\,\,\,\,\,{{k}_{\phi }}=\frac{{\partial {{p}_{\phi }}}}{{\partial x}}\,\,\,\,\,\left( {4.2} \right)$

so that, substituting in (3.2), we arrive at

$\displaystyle \mu \left( {\frac{{{{\partial }^{2}}{{u}_{\phi }}}}{{\partial {{r}^{2}}}}+\frac{1}{r}\frac{{\partial {{u}_{\phi }}}}{{\partial r}}} \right)-\rho \frac{{\partial {{u}_{\phi }}}}{{\partial t}}={{k}_{\phi }}\left( t \right)\,\,\,\,\,\left( 5 \right)$

This is the desired form of the basic PDE for oscillatory flow, which we now attempt to solve for velocity component $\displaystyle {{u}_{\phi }}$.

3. The oscillatory velocity profile at moderate frequency

We assume that the flow is driven by an oscillatory pressure gradient described by

$\displaystyle {{k}_{\phi }}\left( t \right)={{k}_{s}}{{e}^{{i\omega t}}}\,\,\,\,\,\left( 6 \right)$

where $\displaystyle {{k}_{s}}$ is the steady-state pressure gradient, $\displaystyle i$ is the imaginary number, $\displaystyle \omega$ is the circular frequency of the oscillatory pressure gradient and t is time. As the reader may know, the rightmost complex exponential in (6) can be restated in the trigonometric form

$\displaystyle {{k}_{s}}{{e}^{{i\omega t}}}={{k}_{s}}\left( {\cos \omega t+i\sin \omega t} \right)\,\,\,\,\,\left( 7 \right)$

It follows that a solution to the partial differential equation (5) with the choice of oscillatory pressure gradient (6) will actually consist of two solutions, one for a flow driven by $\displaystyle {{k}_{\phi }}$(t) = $\displaystyle {{k}_{s}}$cos $\displaystyle \omega$t and the other for a flow driven by $\displaystyle {{k}_{\phi }}$(t) = $\displaystyle {{k}_{s}}$sin $\displaystyle \omega$t. The first is obtained by taking the real part of the solution and the second by taking the imaginary part. We’ve chosen to specify the oscillatory pressure gradient with an amplitude equal to its steady-state counterpart $\displaystyle {{k}_{s}}$ so that comparisons can be made between the oscillatory components of the solution and the corresponding results in Poiseuille flow.

The governing equation for oscillatory flow with the choice of pressure gradient (7) becomes

$\displaystyle \frac{{{{\partial }^{2}}{{u}_{\phi }}}}{{\partial {{r}^{2}}}}+\frac{1}{r}\frac{{\partial {{u}_{\phi }}}}{{\partial r}}-\frac{\rho }{\mu }\frac{{\partial {{u}_{\phi }}}}{{\partial t}}=\frac{{{{k}_{s}}}}{\mu }{{e}^{{i\omega t}}}\,\,\,\,\,\left( 8 \right)$

The main difficulty to carry out a solution for PDE (8) is that, as it stands, the equation is a function of two independent variables, namely the radial coordinate r and time t. (This sentence is a bit redundant because every PDE is a function of more than one independent variable.) This can be resolved if we propose a solution by separation of variables, such that $\displaystyle {{u}_{\phi }}$ is decomposed into a part that depends on r only and one that depends on t only. The proposed solution is then

$\displaystyle {{u}_{\phi }}\left( {r,t} \right)={{U}_{\phi }}\left( r \right){{e}^{{i\omega t}}}\,\,\,\,\,\left( 9 \right)$

Substituting this into (8), the factor $\displaystyle {{e}^{{i\omega t}}}$ cancels throughout and we’re left with an ordinary differential equation in $\displaystyle {{U}_{\phi }}$(r) only,

$\displaystyle \frac{{{{d}^{2}}{{U}_{\phi }}}}{{d{{r}^{2}}}}+\frac{1}{r}\frac{{d{{U}_{\phi }}}}{{dr}}-\frac{{i\alpha }}{R}{{U}_{\phi }}=\frac{{{{k}_{s}}}}{\mu }\,\,\,\,\,\left( {10} \right)$

where R is the tube radius and $\displaystyle \alpha$ is the Womersley number. The latter is one of the myriad dimensionless groups found in modelling of transport phenomena, and may be viewed as a dimensionless frequency (Zamir, 2000) or, equivalently, a characteristic time (Truskey et al., 2009) of the flow; it is expressed as

$\displaystyle \alpha =R\sqrt{{\frac{\omega }{\nu }}}\,\,\,\,\,\left( {11} \right)$

where $\displaystyle \nu$ is the kinematic viscosity of the running fluid. Differences in the magnitude of the Womersley number have a significant effect on the shape of the velocity distribution.

Equation (10) is a form of the Bessel equation, a type of second-order ODE that often arises in problems involving cylindrical geometries. Importantly, while the pressure gradient $\displaystyle {{k}_{\phi }}$(t) in equation (6) and the velocity $\displaystyle {{u}_{\phi }}$(r,t) in equation (9) appear to have the same oscillatory form in the time variable t because both are multiplied by $\displaystyle {{e}^{{i\omega t}}}$, this does not imply that the pressure gradient and velocity are in phase with each other. The reason for this is that the radial component $\displaystyle {{U}_{\phi }}$(r) in separated-variable solution (9) is itself a complex entity and, as such, modifies the phase angle of velocity $\displaystyle {{u}_{\phi }}$ when multiplied by the unsteady term $\displaystyle {{e}^{{i\omega t}}}$. As we’ve already noted, differences in phase between the pressure gradient wave $\displaystyle {{k}_{\phi }}$ and other mechanical variables are a distinctive feature of oscillatory flow.

Bessel equation (10) has a general solution of the form

$\displaystyle {{U}_{\phi }}\left( r \right)=\frac{{i{{k}_{s}}{{R}^{2}}}}{{\mu {{\alpha }^{2}}}}+A{{J}_{0}}\left( \zeta \right)+B{{Y}_{0}}\left( \zeta \right)\,\,\,\,\,\left( {12} \right)$

where A, B are constants and J0, Y0 are Bessel functions of order zero of the first and second kind, respectively. Zeta is a normalized radial distance given by

$\displaystyle \zeta \left( r \right)\equiv \Lambda \frac{r}{R}\,\,\,\,\,\left( {13} \right)$

where $\displaystyle \Lambda$ is a parameter that involves the Womersley number:

$\displaystyle \Lambda \equiv \left( {\frac{{i-1}}{{\sqrt{2}}}} \right)\alpha \,\,\,\,\,\left( {14} \right)$

The boundary conditions that the solution must satisfy are the same as those utilized in steady Poiseuille flow: first, we have the no-slip condition, which implies that velocity at the tube wall must equal zero,

$\displaystyle {{U}_{\phi }}\left( R \right)=0\,\,\,\,\,\left( {15.1} \right)$

Second, the velocity must be finite at the centerline of the tube,

$\displaystyle \left| {{{U}_{\phi }}\left( 0 \right)} \right|<\infty \,\,\,\,\,\left( {15.2} \right)$

It is known that Y0($\displaystyle \zeta$) becomes infinite when $\displaystyle \zeta$ $\displaystyle \to$ 0; therefore, boundary condition (15.2) implies that B = 0. Next, applying the first boundary condition (15.1) and noting that

$\displaystyle \zeta \left( {r=R} \right)=\Lambda \times \frac{R}{R}=\Lambda$

we have

$\displaystyle {{U}_{\phi }}\left( r \right)=\frac{{i{{k}_{s}}{{R}^{2}}}}{{\mu {{\alpha }^{2}}}}+A{{J}_{0}}\left( \zeta \right)$

$\displaystyle \therefore {{U}_{\phi }}\left( R \right)=0=\frac{{i{{k}_{s}}{{R}^{2}}}}{{\mu {{\alpha }^{2}}}}+A{{J}_{0}}\left( \Lambda \right)$

$\displaystyle \therefore A{{J}_{0}}\left( \Lambda \right)=-\frac{{i{{k}_{s}}{{R}^{2}}}}{{\mu {{\alpha }^{2}}}}$

$\displaystyle \therefore A=-\frac{{i{{k}_{s}}{{R}^{2}}}}{{\mu {{\alpha }^{2}}{{J}_{0}}\left( \Lambda \right)}}\,\,\,\,\,\left( {16} \right)$

Substituting in (12) ultimately yields

$\displaystyle {{U}_{\phi }}\left( r \right)=\frac{{i{{k}_{s}}{{R}^{2}}}}{{\mu {{\alpha }^{2}}}}+A{{J}_{0}}\left( \zeta \right)+B{{Y}_{0}}\left( \zeta \right)$

$\displaystyle \therefore {{U}_{\phi }}\left( r \right)=\frac{{i{{k}_{s}}{{R}^{2}}}}{{\mu {{\alpha }^{2}}}}-\frac{{i{{k}_{s}}{{R}^{2}}}}{{\mu {{\alpha }^{2}}}}\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}$

$\displaystyle \therefore {{U}_{\phi }}\left( r \right)=\frac{{i{{k}_{s}}{{R}^{2}}}}{{\mu {{\alpha }^{2}}}}\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right)\,\,\,\,\,\left( {17} \right)$

Substituting $\displaystyle {{U}_{\phi }}$ into (9), we find the oscillatory flow velocity

$\displaystyle {{u}_{\phi }}\left( {r,t} \right)={{U}_{\phi }}\left( r \right){{e}^{{i\omega t}}}$

$\displaystyle \therefore {{u}_{\phi }}\left( {r,t} \right)=\frac{{i{{k}_{s}}{{R}^{2}}}}{{\mu {{\alpha }^{2}}}}\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right){{e}^{{i\omega t}}}\,\,\,\,\,\left( {18} \right)$

This is the classical solution for oscillatory flow in a rigid circular tube, first published by Sexl in 1930 and further explored by Womersley (1955) and Uchida (1956).

The first element of the solution is a constant coefficient whose value depends on the amplitude of the pressure gradient $\displaystyle {{k}_{s}}$, the radius of the tube R, the viscosity $\displaystyle \mu$ of the fluid, and the Womersley number $\displaystyle \alpha$. The second element, inside the parentheses, is a function of normalized radial coordinate $\displaystyle \zeta$ and describes the velocity distribution in the circular cross-section. The third element is a function of time, which multiplies the other two terms and thereby describes the temporal evolution of velocity.

It is informative to normalize velocity profile (18) with respect to the centerline velocity in steady Poiseuille flow, which is given by the well-known result

$\displaystyle {{u}_{s}}=\frac{{-{{k}_{s}}{{R}^{2}}}}{{4\mu }}\,\,\,\,\,\left( {19} \right)$

so that

$\displaystyle \frac{{{{u}_{\phi }}\left( {r,t} \right)}}{{{{u}_{s}}}}=\frac{{{{i{{k}_{s}}{{R}^{2}}}}/{{\mu {{\alpha }^{2}}}}\;}}{{-{{{{k}_{s}}{{R}^{2}}}}/{{4\mu }}\;}}\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right){{e}^{{i\omega t}}}$

$\displaystyle \therefore \frac{{{{u}_{\phi }}\left( {r,t} \right)}}{{{{u}_{s}}}}=\frac{{-4i}}{{{{\alpha }^{2}}}}\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right){{e}^{{i\omega t}}}$

$\displaystyle \therefore \frac{{{{u}_{\phi }}\left( {r,t} \right)}}{{{{u}_{s}}}}=\frac{{-4i}}{{\left[ {{2}/{{{{{\left( {i-1} \right)}}^{2}}}}\;} \right]{{\Lambda }^{2}}}}\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right){{e}^{{i\omega t}}}$

$\displaystyle \therefore \frac{{{{u}_{\phi }}\left( {r,t} \right)}}{{{{u}_{s}}}}=\frac{{-4}}{{{{\Lambda }^{2}}}}\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right){{e}^{{i\omega t}}}\,\,\,\,\,\left( {20} \right)$

The nondimensional form of the oscillatory velocity profile is equal to unity when the oscillatory component of velocity, $\displaystyle {{u}_{\phi }}$, equals the maximum velocity in the corresponding Poiseuille flow case, $\displaystyle {{u}_{s}}$.

Since the pressure gradient wave (6) can be broken down into two components, it follows that oscillatory velocity profile (18) actually describes two solutions, one that describes flow in response to the real pressure gradient $\displaystyle {{k}_{{\phi R}}}$ and another that takes place in response to imaginary pressure gradient $\displaystyle {{k}_{{\phi I}}}$, where

$\displaystyle {{k}_{{\phi R}}}={{k}_{s}}\cos \omega t\,\,\,\,\,;\,\,\,\,\,{{k}_{{\phi I}}}={{k}_{s}}\sin \omega t\,\,\,\,\,\left( {21.1,2} \right)$

The corresponding velocities are the real and imaginary parts of $\displaystyle {{u}_{\phi }}$, that is,

$\displaystyle {{u}_{\phi }}={{u}_{{\phi R}}}+i{{u}_{{\phi I}}}={{U}_{\phi }}{{e}^{{i\omega t}}}$

$\displaystyle \therefore {{u}_{{\phi R}}}+i{{u}_{{\phi I}}}={{U}_{\phi }}\left( {\cos \omega t+i\sin \omega t} \right)\,\,\,\,\,\left( {22} \right)$

The real and imaginary parts of $\displaystyle {{U}_{\phi }}$ are denoted by $\displaystyle {{U}_{{\phi R}}}$ and $\displaystyle {{U}_{{\phi I}}}$, respectively, where

$\displaystyle \frac{{{{U}_{{\phi R}}}}}{{{{u}_{s}}}}=\Re \left[ {\frac{{-4}}{{{{\Lambda }^{2}}}}\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right)} \right]\,\,\,\,\,\left( {23.1} \right)$

$\displaystyle \frac{{{{U}_{{\phi I}}}}}{{{{u}_{s}}}}=\Im \left[ {\frac{{-4}}{{{{\Lambda }^{2}}}}\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right)} \right]\,\,\,\,\,\left( {23.2} \right)$

so that

$\displaystyle {{U}_{\phi }}={{U}_{{\phi R}}}+i{{U}_{{\phi I}}}\,\,\,\,\,\left( {23.3} \right)$

which, when substituted into (22), yields

$\displaystyle {{u}_{{\phi R}}}+i{{u}_{{\phi I}}}=\left( {{{U}_{{\phi R}}}+i{{U}_{{\phi I}}}} \right)\left( {\cos \omega t+i\sin \omega t} \right)$

$\displaystyle \therefore {{u}_{{\phi R}}}={{U}_{{\phi R}}}\cos \left( {\omega t} \right)-{{U}_{{\phi I}}}\sin \left( {\omega t} \right)\,\,\,\,\,\left( {24.1} \right)$

$\displaystyle \therefore {{u}_{{\phi I}}}={{U}_{{\phi I}}}\cos \left( {\omega t} \right)+{{U}_{{\phi R}}}\sin \left( {\omega t} \right)\,\,\,\,\,\left( {24.2} \right)$

It is clear from solution (18) that the shape of the oscillatory velocity profile depends critically on the frequency of oscillation, which appears explicitly in $\displaystyle {{e}^{{i\omega t}}}$ but also implicitly in $\displaystyle \Lambda$, $\displaystyle \alpha$, and $\displaystyle \zeta$. The best way to appreciate this effect is to rely on the Womersley number $\displaystyle \alpha$, which, after all, may be regarded as a dimensionless frequency. Figure 1 shows a set of oscillatory velocity profiles corresponding to $\displaystyle \alpha$ = 3 and to the real part of the pressure gradient, namely $\displaystyle {{k}_{s}}$cos $\displaystyle \omega$t. It is observed that velocity profiles oscillate between the forward direction and the backward direction, but neither the phase nor amplitude of these peak profiles correspond with the peaks of the oscillatory pressure. In particular, the maximum velocity in the peak velocity profile is less than 1.0, which means that it is less than the maximum velocity in the corresponding Poiseuille profile.

Figure 1. Oscillatory velocity profiles $\displaystyle {{u}_{\phi }}$ (normalized by $\displaystyle {{u}_{s}}$) for Womersley number $\displaystyle \alpha$ = 3.0.

3.1. Oscillatory flow rate and shear stress at moderate frequency

Oscillatory flow rate. The volumetric flow rate $\displaystyle {{q}_{\phi }}$ in oscillatory flow through a tube is obtained by integrating the oscillatory velocity profile over a cross-section of the tube. Since the oscillatory velocity $\displaystyle {{u}_{\phi }}$(r,t) is a function of r and t, the result is a function of time given by

$\displaystyle {{q}_{\phi }}\left( t \right)=\int_{0}^{R}{{{{u}_{\phi }}\left( {r,t} \right)\times 2\pi rdr}}$

$\displaystyle \therefore {{q}_{\phi }}\left( t \right)=\frac{{2\pi i{{k}_{s}}{{R}^{2}}}}{{\mu {{\alpha }^{2}}}}{{e}^{{i\omega t}}}\int_{0}^{R}{{\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right)\times rdr}}\,\,\,\,\,\left( {25} \right)$

To evaluate the integral on the right-hand side, we first use (13) to replace radial coordinate r with its normalized form $\displaystyle \zeta$,

$\displaystyle \int_{0}^{R}{{\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right)\times rdr}}=\int_{0}^{\Lambda }{{\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right)\times \frac{{{{R}^{2}}}}{{{{\Lambda }^{2}}}}\zeta d\zeta }}$

$\displaystyle \therefore \int_{0}^{R}{{\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right)\times rdr}}=\frac{{{{R}^{2}}}}{{{{\Lambda }^{2}}}}\int_{0}^{\Lambda }{{\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right)\zeta d\zeta }}$

$\displaystyle \therefore \int_{0}^{R}{{\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right)\times rdr}}=\frac{{{{R}^{2}}}}{{{{\Lambda }^{2}}{{J}_{0}}\left( \Lambda \right)}}\int_{0}^{\Lambda }{{\left( {{{J}_{0}}\left( \Lambda \right)-{{J}_{0}}\left( \zeta \right)} \right)\zeta d\zeta }}\,\,\,\,\,\left( {26} \right)$

Then, we carry out the integration as usual,

$\displaystyle \int_{0}^{R}{{\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right)\times rdr}}=\frac{{{{R}^{2}}}}{{{{\Lambda }^{2}}{{J}_{0}}\left( \Lambda \right)}}\int_{0}^{\Lambda }{{\left( {{{J}_{0}}\left( \Lambda \right)\zeta -{{J}_{0}}\left( \zeta \right)\zeta } \right)d\zeta }}$

$\displaystyle \therefore \int_{0}^{R}{{\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right)\times rdr}}=\frac{{{{R}^{2}}}}{{{{\Lambda }^{2}}{{J}_{0}}\left( \Lambda \right)}}\left. {\left[ {\frac{{{{J}_{0}}\left( \Lambda \right){{\zeta }^{2}}}}{2}-{{J}_{1}}\left( \zeta \right)\zeta } \right]} \right|_{{\zeta =0}}^{{\zeta =\Lambda }}$

$\displaystyle \therefore \int_{0}^{R}{{\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right)\times rdr}}=\frac{{{{R}^{2}}}}{{{{\Lambda }^{2}}{{J}_{0}}\left( \Lambda \right)}}\left( {\frac{{{{J}_{0}}\left( \Lambda \right){{\Lambda }^{2}}}}{2}-{{J}_{1}}\left( \Lambda \right)\Lambda } \right)$

$\displaystyle \therefore \int_{0}^{R}{{\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right)\times rdr}}=\frac{{{{R}^{2}}}}{{{{\Lambda }^{2}}{{J}_{0}}\left( \Lambda \right)}}\left( {\frac{{{{J}_{0}}\left( \Lambda \right){{\Lambda }^{2}}}}{2}-\frac{{{{J}_{1}}\left( \Lambda \right){{\Lambda }^{2}}}}{\Lambda }} \right)$

$\displaystyle \therefore \int_{0}^{R}{{\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right)\times rdr}}=\frac{{{{R}^{2}}}}{{{{J}_{0}}\left( \Lambda \right)}}\left( {\frac{{{{J}_{0}}\left( \Lambda \right)}}{2}-\frac{{2{{J}_{1}}\left( \Lambda \right)}}{{2\Lambda }}} \right)$

$\displaystyle \therefore \int_{0}^{R}{{\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right)\times rdr}}=\frac{{{{R}^{2}}}}{2}\left( {1-\frac{{2{{J}_{1}}\left( \Lambda \right)}}{{\Lambda {{J}_{0}}\left( \Lambda \right)}}} \right)\,\,\,\,\,\left( {27} \right)$

where J1 is a Bessel function of first order and first kind, related to J0 by

$\displaystyle \int{{\zeta {{J}_{0}}\left( \zeta \right)d\zeta }}=\zeta {{J}_{1}}\left( \zeta \right)\,\,\,\,\,\left( {28} \right)$

Substituting integration result (27) into (26) gives the oscillatory flow rate

$\displaystyle {{q}_{\phi }}\left( t \right)=\frac{{2\pi i{{k}_{s}}{{R}^{2}}}}{{\mu {{\alpha }^{2}}}}\int_{0}^{R}{{\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right)\times rdr}}$

$\displaystyle \therefore {{q}_{\phi }}\left( t \right)=\frac{{2\pi i{{k}_{s}}{{R}^{2}}}}{{\mu {{\alpha }^{2}}}}\times \frac{{{{R}^{2}}}}{2}\left( {1-\frac{{2{{J}_{1}}\left( \Lambda \right)}}{{\Lambda {{J}_{0}}\left( \Lambda \right)}}} \right)$

$\displaystyle \therefore {{q}_{\phi }}\left( t \right)=\frac{{i\pi {{k}_{s}}{{R}^{4}}}}{{\mu {{\alpha }^{2}}}}\left( {1-\frac{{2{{J}_{1}}\left( \Lambda \right)}}{{\Lambda {{J}_{0}}\left( \Lambda \right)}}} \right){{e}^{{i\omega t}}}\,\,\,\,\,\left( {29} \right)$

The net flow rate $\displaystyle {{Q}_{\phi }}$ over one oscillatory cycle is given by

$\displaystyle {{Q}_{\phi }}=\int_{0}^{T}{{{{q}_{\phi }}\left( t \right)dt}}$

$\displaystyle \therefore {{Q}_{\phi }}=\frac{{i\pi {{k}_{s}}{{R}^{4}}}}{{\mu {{\alpha }^{2}}}}\left( {1-\frac{{2{{J}_{1}}\left( \Lambda \right)}}{{\Lambda {{J}_{0}}\left( \Lambda \right)}}} \right)\int_{0}^{T}{{{{e}^{{i\omega t}}}dt}}$

$\displaystyle \therefore {{Q}_{\phi }}=\frac{{i\pi {{k}_{s}}{{R}^{4}}}}{{\mu {{\alpha }^{2}}}}\left( {1-\frac{{2{{J}_{1}}\left( \Lambda \right)}}{{\Lambda {{J}_{0}}\left( \Lambda \right)}}} \right)\int_{0}^{T}{{\left( {\cos \omega t+i\sin \omega t} \right)dt}}$

$\displaystyle \therefore {{Q}_{\phi }}=\frac{{i\pi {{k}_{s}}{{R}^{4}}}}{{\mu {{\alpha }^{2}}}}\left( {1-\frac{{2{{J}_{1}}\left( \Lambda \right)}}{{\Lambda {{J}_{0}}\left( \Lambda \right)}}} \right)\int_{0}^{{2\pi }}{{\left( {\cos \theta +i\sin \theta } \right)d\theta }}$

$\displaystyle \therefore {{Q}_{\phi }}=0\,\,\,\,\,\left( {30} \right)$

where T = 2$\displaystyle \pi$/$\displaystyle \omega$ is the period of integration. This result indicates that in oscillatory flow the fluid moves back and forth, with no net flow in either direction.

Analysis of oscillatory flow rate (29) is made more instructive if we normalize it with respect to the flow rate for steady Poiseuille flow, which is given by the well-known expression

$\displaystyle {{q}_{s}}=\frac{{-{{k}_{s}}\pi {{R}^{4}}}}{{8\mu }}\,\,\,\,\,\left( {31} \right)$

so that

$\displaystyle \frac{{{{q}_{\phi }}\left( t \right)}}{{{{q}_{s}}}}=\frac{{{{i\pi {{k}_{s}}{{R}^{4}}}}/{{\mu {{\alpha }^{2}}}}\;}}{{-{{{{k}_{s}}\pi {{R}^{4}}}}/{{8\mu }}\;}}\left( {1-\frac{{2{{J}_{1}}\left( \Lambda \right)}}{{\Lambda {{J}_{0}}\left( \Lambda \right)}}} \right){{e}^{{i\omega t}}}$

$\displaystyle \therefore \frac{{{{q}_{\phi }}\left( t \right)}}{{{{q}_{s}}}}=\frac{{{i}/{{{{\alpha }^{2}}}}\;}}{{-{1}/{8}\;}}\left( {1-\frac{{2{{J}_{1}}\left( \Lambda \right)}}{{\Lambda {{J}_{0}}\left( \Lambda \right)}}} \right){{e}^{{i\omega t}}}$

$\displaystyle \therefore \frac{{{{q}_{\phi }}\left( t \right)}}{{{{q}_{s}}}}=\frac{{-8i}}{{{{\alpha }^{2}}}}\left( {1-\frac{{2{{J}_{1}}\left( \Lambda \right)}}{{\Lambda {{J}_{0}}\left( \Lambda \right)}}} \right){{e}^{{i\omega t}}}\,\,\,\,\,\left( {32} \right)$

But, squaring (14),

$\displaystyle {{\Lambda }^{2}}={{\left( {\frac{{i-1}}{{\sqrt{2}}}} \right)}^{2}}{{\alpha }^{2}}$

$\displaystyle \therefore {{\Lambda }^{2}}=-i{{\alpha }^{2}}$

$\displaystyle \therefore \frac{{{{\Lambda }^{2}}}}{{-i}}={{\alpha }^{2}}\,\,\,\,\,\left( {33} \right)$

Substituting in (32),

$\displaystyle \frac{{{{q}_{\phi }}\left( t \right)}}{{{{q}_{s}}}}=\frac{{-8}}{{{{\Lambda }^{2}}}}\left( {1-\frac{{2{{J}_{1}}\left( \Lambda \right)}}{{\Lambda {{J}_{0}}\left( \Lambda \right)}}} \right){{e}^{{i\omega t}}}\,\,\,\,\,\left( {34} \right)$

This ratio represents the oscillatory flow rate scaled in terms of the corresponding flow rate  in Poiseuille flow, so that a value of 1.0 represents a flow rate equal to that of a Poiseuille flow in the same tube and under a constant pressure gradient equal to $\displaystyle {{k}_{s}}$. It is apparent that the flow rate depends heavily on the frequency of oscillation. Flow rate (34) is complex, with a real component corresponding to a cosenoidal oscillatory pressure gradient and an imaginary component corresponding to a sinusoidal pressure gradient.

Figure 2 shows the imaginary component of normalized oscillatory flow rate (34) for a Womersley number $\displaystyle \alpha$ = 3.0. It is seen that the flow rate oscillates between a peak in the forward direction and a peak in the backward direction, but neither peak is particularly close to positive or negative unity; that is, the oscillatory flow falls short of reaching the rate associated with steady Poiseuille flow. Also shown is the imaginary pressure wave $\displaystyle {{k}_{{\phi I}}}$ = $\displaystyle {{k}_{s}}$ sin $\displaystyle \omega$t; it is seen that peak flow occurs later than peak pressure, that is, the flow wave lags the pressure wave.

Figure 2. Oscillatory flow rate $\displaystyle {{q}_{\phi }}$ (normalized by steady-state value $\displaystyle {{q}_{s}}$) for Womersley number $\displaystyle \alpha$ = 3.0.

Oscillatory shear stress. In oscillatory flow, as fluid moves back and forth in response to the oscillatory pressure gradient, the shear stress exerted by the fluid on the tube wall varies correspondingly as a function of time given by

$\displaystyle {{\tau }_{\phi }}\left( t \right)=-\mu {{\left. {\left( {\frac{{\partial {{u}_{\phi }}\left( {r,t} \right)}}{{\partial r}}} \right)} \right|}_{{r=R}}}\,\,\,\,\,\left( {35} \right)$

It is important to note that, since this shear is produced by only the oscillatory velocity $\displaystyle {{u}_{\phi }}$, it occurs in addition to that produced by the steady velocity $\displaystyle {{u}_{s}}$, when present. Substituting (18) into (35) and carrying out the differentiation, we obtain

$\displaystyle {{\tau }_{\phi }}\left( t \right)=-\mu \frac{\partial }{{\partial r}}{{\left. {\left[ {\frac{{i{{k}_{s}}{{R}^{2}}}}{{\mu {{\alpha }^{2}}}}\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right){{e}^{{i\omega t}}}} \right]} \right|}_{{r=R}}}$

$\displaystyle \therefore {{\tau }_{\phi }}\left( t \right)=-\frac{{i{{k}_{s}}{{R}^{2}}}}{{{{\alpha }^{2}}}}\frac{\partial }{{\partial r}}{{\left. {\left[ {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right]} \right|}_{{r=R}}}{{e}^{{i\omega t}}}$

$\displaystyle \therefore {{\tau }_{\phi }}\left( t \right)=-\frac{{i{{k}_{s}}{{R}^{2}}}}{{{{\alpha }^{2}}}}\frac{\partial }{{\partial \zeta }}{{\left. {\left[ {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right]} \right|}_{{\zeta =\Lambda }}}\frac{\Lambda }{R}{{e}^{{i\omega t}}}$

$\displaystyle \therefore {{\tau }_{\phi }}\left( t \right)=-\frac{{i{{k}_{s}}{{R}^{2}}}}{{{{\alpha }^{2}}}}\frac{\Lambda }{R}\left( {\frac{{{{J}_{1}}\left( \Lambda \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right){{e}^{{i\omega t}}}$

$\displaystyle \therefore {{\tau }_{\phi }}\left( t \right)=-\frac{{i{{k}_{s}}{{R}^{2}}}}{{{{{{\Lambda }^{2}}}}/{{-i}}\;}}\frac{\Lambda }{R}\left( {\frac{{{{J}_{1}}\left( \Lambda \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right){{e}^{{i\omega t}}}$

$\displaystyle \therefore {{\tau }_{\phi }}\left( t \right)=-\frac{{{{k}_{s}}R}}{\Lambda }\left( {\frac{{{{J}_{1}}\left( \Lambda \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right){{e}^{{i\omega t}}}\,\,\,\,\,\left( {36} \right)$

Note that we’ve used the Bessel function property

$\displaystyle \frac{{d{{J}_{0}}\left( \zeta \right)}}{{d\zeta }}=-{{J}_{1}}\left( \zeta \right)\,\,\,\,\,\left( {37} \right)$

As before, it is convenient to put (36) into nondimensional form by scaling the oscillatory shear stress by the corresponding shear stress in Poiseuille flow, which is given by

$\displaystyle {{\tau }_{s}}=\frac{{-{{k}_{s}}R}}{2}\,\,\,\,\,\left( {38} \right)$

so that

$\displaystyle {{\tau }_{\phi }}\left( t \right)=\frac{{{{-{{k}_{s}}R}}/{\Lambda }\;}}{{{{-{{k}_{s}}R}}/{2}\;}}\left( {\frac{{{{J}_{1}}\left( \Lambda \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right){{e}^{{i\omega t}}}$

$\displaystyle \therefore {{\tau }_{\phi }}\left( t \right)=\frac{2}{\Lambda }\left( {\frac{{{{J}_{1}}\left( \Lambda \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right){{e}^{{i\omega t}}}\,\,\,\,\,\left( {39} \right)$

Oscillatory shear (39) is complex, with a real component corresponding to a real oscillatory pressure gradient and an imaginary component corresponding to an imaginary pressure gradient. Figure 3 shows the imaginary component of normalized wall shear $\displaystyle {{\tau }_{\phi }}$ for a Womersley number $\displaystyle \alpha$ = 3.0, along with the sinusoidal pressure wave. It is seen that the shear lags the pressure, somewhat like the imaginary flow rate in Figure 2. The shear oscillates back and forth in each direction, peaking at a maximum value close to one half the steady Poiseuille flow value.

In pulsatile flow consisting of steady- and oscillatory-flow components, the oscillatory shear stress adds to and subtracts from the steady shear stress. The results in Figure 3 show that in pulsatile flow at $\displaystyle \alpha$ = 3.0 the shear stress would oscillate between a high of approximately 1.5 to a low of approximately 0.5 times its constant value in steady Poiseuille flow.

Figure 3. Oscillatory wall shear stress $\displaystyle {{\tau }_{\phi }}$ (normalized by steady-state value $\displaystyle {{\tau }_{s}}$) for Womersley number $\displaystyle \alpha$ = 3.0.

4. Oscillatory velocity profile at low frequency

When subjected to an unsteady pressure gradient of low frequency, flow is better able to keep pace with the changing pressure. Indeed, in the limit of near-zero frequency, the relation between flow and pressure becomes instantaneously the same as in steady Poiseuille flow. That is, at each point in time within the oscillatory cycle, the velocity profile is what it would be in steady Poiseuille flow under a pressure gradient equal to the value of the oscillatory pressure gradient at that instant. Zamir (2000) suggests the term oscillatory Poiseuille flow for this type of behavior.

A series expansion of the zero-order Bessel function of first kind for small z reads as

$\displaystyle {{J}_{0}}\left( z \right)=1-\frac{{{{z}^{2}}}}{{{{2}^{2}}}}+\frac{{{{z}^{4}}}}{{{{2}^{2}}\times {{4}^{2}}}}-\frac{{{{z}^{6}}}}{{{{2}^{2}}\times {{4}^{2}}\times {{6}^{2}}}}+\ldots \,\,\,\,\,\left( {40} \right)$

This expansion is valid for complex values of the independent variable z as required for application to the complex solution of the pulsatile flow equation. Thus, an approximation for the quotient term in equation (18) for the velocity profile, using only the first three terms of the series, and recalling that $\displaystyle \zeta$ = $\displaystyle \Lambda$r/R, is given by

$\displaystyle \frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}=\frac{{{{J}_{0}}\left( {{{\Lambda r}}/{R}\;} \right)}}{{{{J}_{0}}\left( \Lambda \right)}}$

$\displaystyle \therefore \frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}\approx \left( {1-\frac{{{{\Lambda }^{2}}{{r}^{2}}}}{{4{{R}^{2}}}}+\frac{{{{\Lambda }^{4}}{{r}^{4}}}}{{64{{R}^{4}}}}} \right)\times {{\left( {1-\frac{{{{\Lambda }^{2}}}}{4}+\frac{{{{\Lambda }^{4}}}}{{64}}} \right)}^{{-1}}}$

$\displaystyle \therefore \frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}\approx \left( {1-\frac{{{{\Lambda }^{2}}{{r}^{2}}}}{{4{{R}^{2}}}}+\frac{{{{\Lambda }^{4}}{{r}^{4}}}}{{64{{R}^{4}}}}} \right)\times \left( {1+\frac{{{{\Lambda }^{2}}}}{4}+\frac{{3{{\Lambda }^{4}}}}{{64}}} \right)$

$\displaystyle \therefore \frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}\approx 1+\frac{{{{\Lambda }^{2}}}}{4}\left[ {\left( {1-\frac{{{{r}^{2}}}}{{{{R}^{2}}}}} \right)+\frac{{{{\Lambda }^{2}}}}{{16}}\left( {3-\frac{{4{{r}^{2}}}}{{{{R}^{2}}}}+\frac{{{{r}^{4}}}}{{{{R}^{4}}}}} \right)} \right]\,\,\,\,\,\left( {41} \right)$

where in each step we’ve retained only terms of order $\displaystyle {{\zeta }^{4}}$ or lower. Substituting the expression above into equation (20) for the dimensionless velocity profile, we obtain the approximate expression

$\displaystyle \frac{{{{u}_{\phi }}\left( {r,t} \right)}}{{{{u}_{s}}}}\approx \left[ {\left( {1-\frac{{{{r}^{2}}}}{{{{R}^{2}}}}} \right)-\frac{{i{{\alpha }^{2}}}}{{16}}\left( {3-\frac{{4{{r}^{2}}}}{{{{R}^{2}}}}+\frac{{{{r}^{4}}}}{{{{R}^{4}}}}} \right)} \right]{{e}^{{i\omega t}}}\,\,\,\,\,\left( {42} \right)$

The real and imaginary parts of the velocity are then given by

$\displaystyle \frac{{{{u}_{{\phi R}}}\left( {r,t} \right)}}{{{{u}_{s}}}}\approx \left( {1-\frac{{{{r}^{2}}}}{{{{R}^{2}}}}} \right)\cos \omega t+\frac{{{{\alpha }^{2}}}}{{16}}\left( {3-\frac{{4{{r}^{2}}}}{{{{R}^{2}}}}+\frac{{{{r}^{4}}}}{{{{R}^{4}}}}} \right)\sin \omega t\,\,\,\,\,\left( {43.1} \right)$

$\displaystyle \frac{{{{u}_{{\phi I}}}\left( {r,t} \right)}}{{{{u}_{s}}}}\approx \left( {1-\frac{{{{r}^{2}}}}{{{{R}^{2}}}}} \right)\sin \omega t-\frac{{{{\alpha }^{2}}}}{{16}}\left( {3-\frac{{4{{r}^{2}}}}{{{{R}^{2}}}}+\frac{{{{r}^{4}}}}{{{{R}^{4}}}}} \right)\cos \omega t\,\,\,\,\,\left( {43.2} \right)$

These expressions are easier to use than equation (20) because they do not involve Bessel functions and can be used in place of that equation when the frequency is small (i.e, the Womersley number is small). Furthermore, substituting $\displaystyle {{u}_{s}}$ from equation (19) and using equation (7) for the real and imaginary parts of the pressure gradient, we have

$\displaystyle \frac{{{{u}_{{\phi R}}}\left( {r,t} \right)}}{{{{-{{k}_{s}}{{R}^{2}}}}/{{4\mu }}\;}}\approx \left( {1-\frac{{{{r}^{2}}}}{{{{R}^{2}}}}} \right)\cos \omega t+\frac{{{{\alpha }^{2}}}}{{16}}\left( {3-\frac{{4{{r}^{2}}}}{{{{R}^{2}}}}+\frac{{{{r}^{4}}}}{{{{R}^{4}}}}} \right)\sin \omega t$

$\displaystyle \therefore {{u}_{{\phi R}}}\left( {r,t} \right)\approx -\frac{{{{k}_{s}}{{R}^{2}}}}{{4\mu }}\left[ {\left( {1-\frac{{{{r}^{2}}}}{{{{R}^{2}}}}} \right)\cos \omega t+\frac{{{{\alpha }^{2}}}}{{16}}\left( {3-\frac{{4{{r}^{2}}}}{{{{R}^{2}}}}+\frac{{{{r}^{4}}}}{{{{R}^{4}}}}} \right)\sin \omega t} \right]$

$\displaystyle \therefore {{u}_{{\phi R}}}\left( {r,t} \right)\approx -\frac{{{{k}_{s}}\cos \left( {\omega t} \right){{R}^{2}}}}{{4\mu }}\left[ {\left( {1-\frac{{{{r}^{2}}}}{{{{R}^{2}}}}} \right)+\frac{{{{\alpha }^{2}}}}{{16}}\left( {3-\frac{{4{{r}^{2}}}}{{{{R}^{2}}}}+\frac{{{{r}^{4}}}}{{{{R}^{4}}}}} \right)\tan \omega t} \right]$

$\displaystyle \therefore {{u}_{{\phi R}}}\left( {r,t} \right)\approx -\frac{{{{k}_{{\phi R}}}{{R}^{2}}}}{{4\mu }}\left[ {\left( {1-\frac{{{{r}^{2}}}}{{{{R}^{2}}}}} \right)+\frac{{{{\alpha }^{2}}}}{{16}}\left( {3-\frac{{4{{r}^{2}}}}{{{{R}^{2}}}}+\frac{{{{r}^{4}}}}{{{{R}^{4}}}}} \right)\tan \omega t} \right]\,\,\,\,\,\left( {44.1} \right)$

$\displaystyle \frac{{{{u}_{{\phi I}}}\left( {r,t} \right)}}{{{{-{{k}_{s}}{{R}^{2}}}}/{{4\mu }}\;}}\approx \left( {1-\frac{{{{r}^{2}}}}{{{{R}^{2}}}}} \right)\sin \omega t-\frac{{{{\alpha }^{2}}}}{{16}}\left( {3-\frac{{4{{r}^{2}}}}{{{{R}^{2}}}}+\frac{{{{r}^{4}}}}{{{{R}^{4}}}}} \right)\cos \omega t$

$\displaystyle \therefore {{u}_{{\phi I}}}\left( {r,t} \right)\approx -\frac{{{{k}_{s}}{{R}^{2}}}}{{4\mu }}\left[ {\left( {1-\frac{{{{r}^{2}}}}{{{{R}^{2}}}}} \right)\sin \omega t-\frac{{{{\alpha }^{2}}}}{{16}}\left( {3-\frac{{4{{r}^{2}}}}{{{{R}^{2}}}}+\frac{{{{r}^{4}}}}{{{{R}^{4}}}}} \right)\cos \omega t} \right]$

$\displaystyle \therefore {{u}_{{\phi I}}}\left( {r,t} \right)\approx -\frac{{{{k}_{s}}\sin \left( {\omega t} \right){{R}^{2}}}}{{4\mu }}\left[ {\left( {1-\frac{{{{r}^{2}}}}{{{{R}^{2}}}}} \right)-\frac{{{{\alpha }^{2}}}}{{16}}\left( {3-\frac{{4{{r}^{2}}}}{{{{R}^{2}}}}+\frac{{{{r}^{4}}}}{{{{R}^{4}}}}} \right)\cot \omega t} \right]$

$\displaystyle \therefore {{u}_{{\phi I}}}\left( {r,t} \right)\approx -\frac{{{{k}_{{\phi I}}}{{R}^{2}}}}{{4\mu }}\left[ {\left( {1-\frac{{{{r}^{2}}}}{{{{R}^{2}}}}} \right)-\frac{{{{\alpha }^{2}}}}{{16}}\left( {3-\frac{{4{{r}^{2}}}}{{{{R}^{2}}}}+\frac{{{{r}^{4}}}}{{{{R}^{4}}}}} \right)\cot \omega t} \right]\,\,\,\,\,\left( {44.2} \right)$

At low frequency, the term that multiplies $\displaystyle \alpha$ in (44.1) and (44.2) can be neglected, and the relation between velocity and pressure becomes

$\displaystyle {{u}_{{\phi R}}}\left( {r,t} \right)\approx \frac{{{{k}_{{\phi R}}}}}{{4\mu }}\left( {{{R}^{2}}-{{r}^{2}}} \right)\,\,\,\,\,\left( {45.1} \right)$

$\displaystyle {{u}_{{\phi I}}}\left( {r,t} \right)\approx \frac{{{{k}_{{\phi I}}}}}{{4\mu }}\left( {{{R}^{2}}-{{r}^{2}}} \right)\,\,\,\,\,\left( {45.2} \right)$

These expressions are identical to the parabolic velocity profile for steady Poiseuille flow, but with instantaneous values of velocity and pressure, thereby justifying the term oscillatory Poiseuille flow. Some representative oscillatory velocity profiles corresponding to $\displaystyle \alpha$ = 1.0 and to the real part of pressure gradient (7) are shown in Figure 4. Here we have two main differences relatively to the profiles illustrated earlier in Fig. 1, which referred to a Womersley number of 3.0. First, the velocity profiles in Fig. 4 reach their peak form at approximately 0o and 180o, that is, at the same phase angles as the pressure gradient cos $\displaystyle \omega$t; there is essentially no phase lag. Second, the normalized velocity in Fig. 4 peaks very close to unity, indicating that the maximum velocity in this low-frequency limit is almost the same as the peak velocity observed in steady Poiseuille flow.

Figure 4. Oscillatory velocity profiles $\displaystyle {{u}_{\phi }}$ (normalized by $\displaystyle {{u}_{s}}$) for Womersley number $\displaystyle \alpha$ = 1.0.

4.1. Oscillatory flow rate and shear stress at low frequency

Oscillatory flow rate. For the flow rate in ‘oscillatory Poiseuille flow,’ we need the series expansion for J1($\displaystyle \zeta$) for small values of $\displaystyle \zeta$, which reads as

$\displaystyle {{J}_{1}}\left( \zeta \right)\approx\frac{\zeta }{2}-\frac{{{{\zeta }^{3}}}}{{{{2}^{2}}\times 4}}+\frac{{{{\zeta }^{5}}}}{{{{2}^{2}}\times {{4}^{2}}\times 6}}+\frac{{{{\zeta }^{7}}}}{{{{2}^{2}}\times {{4}^{2}}\times {{6}^{2}}\times 8}}+\ldots \,\,\,\,\,\left( {46} \right)$

Using only the first three terms of the series for the quotient term in equation (32), we find

$\displaystyle \frac{{{{J}_{1}}\left( \Lambda \right)}}{{{{J}_{0}}\left( \Lambda \right)}}\approx \frac{\Lambda }{2}\left( {1-\frac{{{{\Lambda }^{2}}}}{8}+\frac{{{{\Lambda }^{4}}}}{{192}}} \right)\times {{\left( {1-\frac{{{{\Lambda }^{2}}}}{4}+\frac{{{{\Lambda }^{4}}}}{{64}}} \right)}^{{-1}}}$

$\displaystyle \therefore \frac{{{{J}_{1}}\left( \Lambda \right)}}{{{{J}_{0}}\left( \Lambda \right)}}\approx \frac{\Lambda }{2}\left( {1-\frac{{{{\Lambda }^{2}}}}{8}+\frac{{{{\Lambda }^{4}}}}{{192}}} \right)\times \left( {1+\frac{{{{\Lambda }^{2}}}}{4}+\frac{{3{{\Lambda }^{4}}}}{{64}}} \right)$

$\displaystyle \therefore \frac{{{{J}_{1}}\left( \Lambda \right)}}{{{{J}_{0}}\left( \Lambda \right)}}\approx \frac{\Lambda }{2}\left( {1+\frac{{{{\Lambda }^{2}}}}{8}+\frac{{{{\Lambda }^{4}}}}{{48}}} \right)\,\,\,\,\,\left( {47} \right)$

where, again, only terms of order $\displaystyle {{\Lambda }^{4}}$ or lower were retained in each step. Substituting this result into (32),

$\displaystyle \frac{{{{q}_{\phi }}\left( t \right)}}{{{{q}_{s}}}}=\frac{{-8i}}{{{{\alpha }^{2}}}}\left( {1-\frac{{2{{J}_{1}}\left( \Lambda \right)}}{{\Lambda {{J}_{0}}\left( \Lambda \right)}}} \right){{e}^{{i\omega t}}}$

$\displaystyle \therefore \frac{{{{q}_{\phi }}\left( t \right)}}{{{{q}_{s}}}}=\frac{{-8i}}{{{{\alpha }^{2}}}}\left[ {1-\frac{2}{\Lambda }\times \frac{\Lambda }{2}\left( {1+\frac{{{{\Lambda }^{2}}}}{8}+\frac{{{{\Lambda }^{4}}}}{{48}}} \right)} \right]{{e}^{{i\omega t}}}$

$\displaystyle \therefore \frac{{{{q}_{\phi }}\left( t \right)}}{{{{q}_{s}}}}=\frac{{-8i}}{{{{\alpha }^{2}}}}\left[ {1-1-\frac{{{{\Lambda }^{2}}}}{8}-\frac{{{{\Lambda }^{4}}}}{{48}}} \right]{{e}^{{i\omega t}}}$

$\displaystyle \therefore \frac{{{{q}_{\phi }}\left( t \right)}}{{{{q}_{s}}}}\approx \frac{{8i}}{{{{\alpha }^{2}}}}\left( {\frac{{{{\Lambda }^{2}}}}{8}+\frac{{{{\Lambda }^{4}}}}{{48}}} \right){{e}^{{i\omega t}}}\,\,\,\,\,\left( {48} \right)$

Noting that $\displaystyle \Lambda$2 = –i$\displaystyle \alpha$2 and $\displaystyle \Lambda$4 = –$\displaystyle \alpha$4, we can restate the above relation as

$\displaystyle \frac{{{{q}_{\phi }}\left( t \right)}}{{{{q}_{s}}}}=\frac{{8i}}{{{{\alpha }^{2}}}}\left( {\frac{{-i{{\alpha }^{2}}}}{8}+\frac{{-{{\alpha }^{4}}}}{{48}}} \right){{e}^{{i\omega t}}}$

$\displaystyle \therefore \frac{{{{q}_{\phi }}\left( t \right)}}{{{{q}_{s}}}}=8i\left( {\frac{{-i}}{8}+\frac{{-{{\alpha }^{2}}}}{{48}}} \right){{e}^{{i\omega t}}}$

$\displaystyle \therefore \frac{{{{q}_{\phi }}\left( t \right)}}{{{{q}_{s}}}}=\left( {-\frac{{8{{i}^{2}}}}{8}-\frac{{{{\alpha }^{2}}i}}{6}} \right){{e}^{{i\omega t}}}$

$\displaystyle \therefore \frac{{{{q}_{\phi }}\left( t \right)}}{{{{q}_{s}}}}=\left( {1-\frac{{{{\alpha }^{2}}i}}{6}} \right){{e}^{{i\omega t}}}$

$\displaystyle \therefore \frac{{{{q}_{\phi }}\left( t \right)}}{{{{q}_{s}}}}\approx \left( {1-\frac{{{{\alpha }^{2}}i}}{6}} \right)\left( {\cos \omega t+i\sin \omega t} \right)\,\,\,\,\,\left( {49} \right)$

Thus, the real and imaginary parts of the flow rate are given by

$\displaystyle \frac{{{{q}_{{\phi R}}}\left( t \right)}}{{{{q}_{s}}}}\approx \cos \omega t+\frac{{{{\alpha }^{2}}}}{6}\sin \omega t\,\,\,\,\,\left( {50.1} \right)$

$\displaystyle \frac{{{{q}_{{\phi I}}}\left( t \right)}}{{{{q}_{s}}}}\approx \sin \omega t-\frac{{{{\alpha }^{2}}}}{6}\cos \omega t\,\,\,\,\,\left( {50.2} \right)$

Substituting $\displaystyle {{q}_{s}}$ from equation (31),

$\displaystyle \frac{{{{q}_{{\phi R}}}\left( t \right)}}{{{{q}_{s}}}}=\frac{{{{q}_{{\phi R}}}\left( t \right)}}{{{{-{{k}_{s}}\pi {{R}^{4}}}}/{{8\mu }}\;}}\approx \cos \omega t+\frac{{{{\alpha }^{2}}}}{6}\sin \omega t$

$\displaystyle \therefore {{q}_{{\phi R}}}\left( t \right)\approx \frac{{-{{k}_{s}}\pi {{R}^{4}}}}{{8\mu }}\left( {\cos \omega t+\frac{{{{\alpha }^{2}}}}{6}\sin \omega t} \right)$

$\displaystyle \therefore {{q}_{{\phi R}}}\left( t \right)\approx \frac{{-{{k}_{s}}\cos \omega t\pi {{R}^{4}}}}{{8\mu }}\left( {1+\frac{{{{\alpha }^{2}}}}{6}\tan \omega t} \right)$

$\displaystyle \therefore {{q}_{{\phi R}}}\left( t \right)\approx \frac{{-{{k}_{{\phi R}}}\pi {{R}^{4}}}}{{8\mu }}\left( {1+\frac{{{{\alpha }^{2}}}}{6}\tan \omega t} \right)\,\,\,\,\,\left( {51.1} \right)$

A similar expression can be obtained for the imaginary component of flow,

$\displaystyle {{q}_{{\phi I}}}\left( t \right)\approx \frac{{-{{k}_{{\phi I}}}\pi {{R}^{4}}}}{{8\mu }}\left( {1+\frac{{{{\alpha }^{2}}}}{6}\cot \omega t} \right)\,\,\,\,\,\left( {51.2} \right)$

At low frequency, such that the second terms in parentheses in (51.1) and (51.2) can be neglected, the relation between flow rate and pressure gradient becomes the same as that in steady flow at each instant, that is,

$\displaystyle {{q}_{{\phi R}}}\left( t \right)\approx \frac{{-{{k}_{{\phi R}}}\pi {{R}^{4}}}}{{8\mu }}\,\,\,\,\,\left( {52.1} \right)$

$\displaystyle {{q}_{{\phi I}}}\left( t \right)\approx \frac{{-{{k}_{{\phi I}}}\pi {{R}^{4}}}}{{8\mu }}\,\,\,\,\,\left( {52.2} \right)$

Figure 5 shows the imaginary component of oscillatory flow for a Womersley number $\displaystyle \alpha$ = 1.0. In contrast to the moderate Womersley number scenario shown in Figure 2, here the flow is nearly in phase with the pressure wave. Moreover, the amplitude is close to one, indicating that the flow rate reaches an intensity nearly equal to that of steady Poiseuille flow.

Figure 5. Oscillatory flow rate $\displaystyle {{q}_{\phi }}$ (normalized by steady-state value $\displaystyle {{q}_{s}}$) for Womersley number $\displaystyle \alpha$ = 1.0.

Oscillatory shear stress. At low frequency, we can introduce series expansion (47) into oscillatory shear stress (36) and manipulate to obtain

$\displaystyle {{\tau }_{\phi }}\left( t \right)\approx -\frac{{{{k}_{s}}R}}{\Lambda }\left( {\frac{{{{J}_{1}}\left( \Lambda \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right){{e}^{{i\omega t}}}$

$\displaystyle \therefore {{\tau }_{\phi }}\left( t \right)\approx -\frac{{{{k}_{s}}R}}{\Lambda }\times \left[ {\frac{\Lambda }{2}\times \left( {1+\frac{{{{\Lambda }^{2}}}}{8}+\frac{{{{\Lambda }^{4}}}}{{48}}} \right)} \right]\times {{e}^{{i\omega t}}}$

$\displaystyle \therefore {{\tau }_{\phi }}\left( t \right)\approx -\frac{{{{k}_{s}}R}}{2}\times \left( {1+\frac{{{{\Lambda }^{2}}}}{8}+\frac{{{{\Lambda }^{4}}}}{{48}}} \right)\times {{e}^{{i\omega t}}}$

$\displaystyle \therefore {{\tau }_{\phi }}\left( t \right)\approx -\frac{{{{k}_{s}}R}}{2}\times \left( {1+\frac{{-i{{\alpha }^{2}}}}{8}+\frac{{-{{\alpha }^{4}}}}{{48}}} \right)\times {{e}^{{i\omega t}}}$

$\displaystyle \therefore {{\tau }_{\phi }}\left( t \right)\approx -\frac{{{{k}_{s}}R}}{2}\times \left( {1-\frac{{{{\alpha }^{2}}i}}{8}-\frac{{{{\alpha }^{4}}}}{{48}}} \right)\times {{e}^{{i\omega t}}}$

$\displaystyle \therefore {{\tau }_{\phi }}\left( t \right)\approx -\frac{{{{k}_{s}}R}}{2}\times \left( {1-\frac{{{{\alpha }^{2}}i}}{8}-\frac{{{{\alpha }^{4}}}}{{48}}} \right)\times \left( {\cos \omega t+i\sin \omega t} \right)\,\,\,\,\,\left( {53} \right)$

Expanding the product of the two terms in parentheses and retaining only terms of order $\displaystyle \alpha$2 or less, we get, for the real and imaginary components of shear,

$\displaystyle {{\tau }_{{\phi R}}}\left( t \right)\approx \frac{{-{{k}_{{\phi R}}}R}}{2}\left( {1+\frac{{{{\alpha }^{2}}}}{8}\tan \omega t} \right)\,\,\,\,\,\left( {54.1} \right)$

$\displaystyle {{\tau }_{{\phi I}}}\left( t \right)\approx \frac{{-{{k}_{{\phi I}}}R}}{2}\left( {1-\frac{{{{\alpha }^{2}}}}{8}\cot \omega t} \right)\,\,\,\,\,\left( {54.2} \right)$

Again, if the frequency is low enough for the trigonometric term inside the parentheses to be neglected, the wall shear becomes instantaneously equal to that of steady Poiseuille flow. Figure 6 shows the imaginary component of normalized wall shear for a Womersley number $\displaystyle \alpha$ = 1.0.

Figure 6. Oscillatory wall shear stress $\displaystyle {{\tau }_{\phi }}$ (normalized by steady-state value $\displaystyle {{\tau }_{s}}$) for Womersley number $\displaystyle \alpha$ = 1.0.

5. Oscillatory velocity profile at high frequency

At high frequency, oscillatory flow in a tube is less able to keep pace with the changing pressure, thus reaching less than the fully developed Poiseuille flow profile at the peak of each cycle. The higher the frequency, the lower the peak velocity the flow is able to reach. In the limit of infinite frequency, the velocity reached at the peak of each cycle is zero, that is, the fluid does not move at all.

An approximate expression for J0($\displaystyle \zeta$) when $\displaystyle \zeta$ is large is given by

$\displaystyle {{J}_{0}}\left( \zeta \right)\approx \frac{{\sin \zeta +\cos \zeta }}{{\sqrt{{\pi \zeta }}}}\,\,\,\,\,\left( {55} \right)$

For the purpose of algebraic manipulation, we write $\displaystyle \zeta$ = i$\displaystyle {{\zeta }_{1}}$, so that

$\displaystyle {{J}_{0}}\left( \zeta \right)\approx \frac{{\sin \left( {i{{\zeta }_{1}}} \right)+\cos \left( {i{{\zeta }_{1}}} \right)}}{{\sqrt{{\pi i{{\zeta }_{1}}}}}}$

$\displaystyle \therefore {{J}_{0}}\left( \zeta \right)\approx \frac{{i\sinh \left( {{{\zeta }_{1}}} \right)+\cosh \left( {{{\zeta }_{1}}} \right)}}{{\sqrt{{\pi i{{\zeta }_{1}}}}}}$

$\displaystyle \therefore {{J}_{0}}\left( \zeta \right)\approx \frac{{\left( {1+i} \right)}}{2}\frac{{{{e}^{{{{\zeta }_{1}}}}}}}{{\sqrt{{\pi i{{\zeta }_{1}}}}}}\,\,\,\,\,\left( {56} \right)$

and similarly, writing $\displaystyle \Lambda$ = i$\displaystyle {{\Lambda }_{1}}$,

$\displaystyle {{J}_{0}}\left( \Lambda \right)\approx \frac{{\left( {1+i} \right)}}{2}\frac{{{{e}^{{{{\Lambda }_{1}}}}}}}{{\sqrt{{\pi i{{\Lambda }_{1}}}}}}\,\,\,\,\,\left( {57} \right)$

Near the tube wall where r/R $\displaystyle \approx$ 1: Inserting the above approximations into equation (20) for the velocity profile, and recalling from (13) that $\displaystyle \zeta$ = $\displaystyle \Lambda$r/R, hence $\displaystyle {{\zeta }_{1}}$ = $\displaystyle {{\Lambda }_{1}}$r/R, we find

$\displaystyle \frac{{{{u}_{\phi }}\left( {r,t} \right)}}{{{{u}_{s}}}}=\frac{{-4}}{{{{\Lambda }^{2}}}}\left( {1-\frac{{{{J}_{0}}\left( \zeta \right)}}{{{{J}_{0}}\left( \Lambda \right)}}} \right){{e}^{{i\omega t}}}$

$\displaystyle \therefore \frac{{{{u}_{\phi }}\left( {r,t} \right)}}{{{{u}_{s}}}}\approx \frac{{-4}}{{{{\Lambda }^{2}}}}\left( {1-\sqrt{{\frac{{{{\Lambda }_{1}}}}{{{{\zeta }_{1}}}}}}{{e}^{{\left( {{{\zeta }_{1}}-{{\Lambda }_{1}}} \right)}}}} \right){{e}^{{i\omega t}}}$

$\displaystyle \therefore \frac{{{{u}_{\phi }}\left( {r,t} \right)}}{{{{u}_{s}}}}\approx \frac{{-4}}{{{{\Lambda }^{2}}}}\left( {1-\sqrt{{\frac{R}{r}}}{{e}^{{{{\Lambda }_{1}}\left( {\frac{r}{R}-1} \right)}}}} \right){{e}^{{i\omega t}}}$

$\displaystyle \therefore \frac{{{{u}_{\phi }}\left( {r,t} \right)}}{{{{u}_{s}}}}\approx \frac{{4i}}{\Lambda }\left( {1-\frac{r}{R}} \right){{e}^{{i\omega t}}}\,\,\,\,\,\left( {58} \right)$

the real and imaginary parts of which are given by

$\displaystyle \frac{{{{u}_{{\phi R}}}\left( {r,t} \right)}}{{{{u}_{s}}}}\approx \frac{{2\sqrt{2}}}{\alpha }\left( {1-\frac{r}{R}} \right)\left( {\cos \omega t+\sin \omega t} \right)\,\,\,\,\,\left( {59.1} \right)$

$\displaystyle \frac{{{{u}_{{\phi I}}}\left( {r,t} \right)}}{{{{u}_{s}}}}\approx \frac{{2\sqrt{2}}}{\alpha }\left( {1-\frac{r}{R}} \right)\left( {\sin \omega t-\cos \omega t} \right)\,\,\,\,\,\left( {59.2} \right)$

Near the center of the tube, where r/R $\displaystyle \approx$ 0, but $\displaystyle \Lambda$ is large because of high frequency: Here the ratio J0($\displaystyle \zeta$)/J0($\displaystyle \Lambda$) in the expression for velocity (eq. (20)) requires an approximation for J0($\displaystyle \zeta$) for small $\displaystyle \zeta$ and an expansion of J0($\displaystyle \Lambda$) for large $\displaystyle \Lambda$. The ratio ultimately vanishes, with the result

$\displaystyle \frac{{{{u}_{\phi }}\left( {r,t} \right)}}{{{{u}_{s}}}}\approx \frac{{-4}}{{{{\Lambda }^{2}}}}{{e}^{{i\omega t}}}\,\,\,\,\,\left( {60} \right)$

the real and imaginary parts of which are given by

$\displaystyle \frac{{{{u}_{{\phi R}}}\left( {r,t} \right)}}{{{{u}_{s}}}}\approx \frac{4}{{{{\alpha }^{2}}}}\sin \omega t\,\,\,\,\,\left( {61.1} \right)$

$\displaystyle \frac{{{{u}_{{\phi I}}}\left( {r,t} \right)}}{{{{u}_{s}}}}\approx \frac{{-4}}{{{{\Omega }^{2}}}}\cos \omega t\,\,\,\,\,\left( {61.2} \right)$

Some representative oscillatory velocity profiles corresponding to $\displaystyle \alpha$ = 10 and to the real part of pressure gradient (7) are shown in Figure 7. This graph illustrates the striking behavior of high-frequency oscillatory flow: because the pressure gradient oscillates so fast, the running fluid has little time to adjust itself and to develop a meaningful flow velocity. As a result, the velocity profiles are nearly flat. (The profiles are so close together that I didn’t even bother to include a legend for the phase angles.)

Figure 7.Oscillatory velocity profiles $\displaystyle {{u}_{\phi }}$ (normalized by $\displaystyle {{u}_{s}}$) for Womersley number $\displaystyle \alpha$ = 10.0.

5.1. Oscillatory flow rate and shear stress at high frequency

For the flow rate, similarly, we use the approximation for large $\displaystyle \zeta$ of the Bessel function of first order, namely

$\displaystyle {{J}_{1}}\left( \zeta \right)\approx \frac{{\sin \zeta -\cos \zeta }}{{\sqrt{{\pi \zeta }}}}\,\,\,\,\,\left( {62} \right)$

As before, writing $\displaystyle \zeta$ = i$\displaystyle {{\zeta }_{1}}$, $\displaystyle \Lambda$ = i$\displaystyle {{\Lambda }_{1}}$ gives

$\displaystyle {{J}_{1}}\left( \zeta \right)\approx \frac{{\left( {i-1} \right){{e}^{{{{\zeta }_{1}}}}}}}{{2\sqrt{{\pi i{{\zeta }_{1}}}}}}$

$\displaystyle {{J}_{1}}\left( \Lambda \right)\approx \frac{{\left( {i-1} \right){{e}^{{{{\Lambda }_{1}}}}}}}{{2\sqrt{{\pi i{{\Lambda }_{1}}}}}}\,\,\,\,\,\left( {63} \right)$

so that, substituting into equation (34), we get

$\displaystyle \frac{{{{q}_{\phi }}\left( t \right)}}{{{{q}_{s}}}}=\frac{{-8i}}{{{{\alpha }^{2}}}}\left( {1-\frac{{2{{J}_{1}}\left( \Lambda \right)}}{{\Lambda {{J}_{1}}\left( \Lambda \right)}}} \right){{e}^{{i\omega t}}}$

$\displaystyle \therefore \frac{{{{q}_{\phi }}\left( t \right)}}{{{{q}_{s}}}}\approx \frac{{-8i}}{{{{\alpha }^{2}}}}\left[ {1-\frac{2}{{i{{\alpha }^{2}}}}\frac{{\left( {i-1} \right)}}{{\left( {i+1} \right)}}} \right]{{e}^{{i\omega t}}}$

$\displaystyle \therefore \frac{{{{q}_{\phi }}\left( t \right)}}{{{{q}_{s}}}}\approx \frac{8}{{{{\alpha }^{2}}}}\left( {\sin \omega t-i\cos \omega t} \right)\,\,\,\,\,\left( {64} \right)$

where the term containing 1/$\displaystyle \alpha$2 is neglected in the last step, and recalling that $\displaystyle \Lambda$2 = –i$\displaystyle \alpha$2. The real and imaginary parts are given by

$\displaystyle \frac{{{{q}_{{\phi R}}}\left( t \right)}}{{{{q}_{s}}}}\approx \left( {\frac{8}{{{{\alpha }^{2}}}}} \right)\sin \left( {\omega t} \right)\,\,\,\,\,\left( {65.1} \right)$

$\displaystyle \frac{{{{q}_{{\phi I}}}\left( t \right)}}{{{{q}_{s}}}}\approx \left( {\frac{{-8}}{{{{\alpha }^{2}}}}} \right)\cos \left( {\omega t} \right)\,\,\,\,\,\left( {65.2} \right)$

These results are similar to those for velocity near the center of the tube, therefore flow rate behaves like velocity near the center of the tube. The equations also show how flow rate achieves ever decreasing values for increasing Womersley $\displaystyle \alpha$ (that is, for increasing frequency). Figure 8 shows the imaginary component of normalized oscillatory flow for a Womersley number $\displaystyle \alpha$ = 10.0. As can be seen, the flow rate barely deviates from zero. Although this is not apparent in the figure, there is an approximately 90o phase difference between flow rate and pressure gradient.

Figure 8. Oscillatory flow rate $\displaystyle {{q}_{\phi }}$ (normalized by steady-state value $\displaystyle {{q}_{s}}$) for Womersley number $\displaystyle \alpha$ = 10.0.

The corresponding results for shear stress, omitting the details, are given by

$\displaystyle \frac{{{{\tau }_{{\phi R}}}\left( {r,t} \right)}}{{{{\tau }_{s}}}}\approx \left( {\frac{{\sqrt{2}}}{\Omega }} \right)\left( {\sin \omega t+\cos \omega t} \right)\,\,\,\,\,\left( {66.1} \right)$

$\displaystyle \frac{{{{\tau }_{{\phi I}}}\left( {r,t} \right)}}{{{{\tau }_{s}}}}\approx \left( {\frac{{\sqrt{2}}}{\Omega }} \right)\left( {\sin \omega t-\cos \omega t} \right)\,\,\,\,\,\left( {66.2} \right)$

Figure 9 shows the imaginary component of normalized wall shear for a Womersley number $\displaystyle \alpha$ = 10.0. Like flow rate, shear stress is very low at high frequency. Again, there is an approximately 90o phase difference between shear stress and pressure gradient.

Figure 9.Oscillatory wall shear stress $\displaystyle {{\tau }_{\phi }}$ (normalized by steady-state value $\displaystyle {{\tau }_{s}}$) for Womersley number $\displaystyle \alpha$ = 10.0.

6. Summary

Our quick investigation of oscillatory flow in a rigid tube showed the following:

1. The behavior of the flow is mainly governed by the frequency of oscillation, as represented by the Womersley number $\displaystyle \alpha$.
2. The velocity distribution moves back and forth in time, yielding profiles with shape governed by the value of $\displaystyle \alpha$. For moderate Womersley numbers, the flow rate and wall shear stress both lag the pressure wave somewhat. This is the case when the oscillatory pressure gradient is a well-defined sine or cosine wave, but other types of response may follow from other kinds of unsteady pressure; see Truskey et al. (2009).
3. For small $\displaystyle \alpha$, the frequency of oscillation is small enough that the flow may be modeled as though it were instantaneously steady. This enables us to pursue an analogy with Hagen-Poiseuille flow and momentarily do away with the Bessel functions that appear in Sexl’s (1930) general solution.
4. For large $\displaystyle \alpha$, the unsteady pressure gradient oscillates so fast that fluid motion can barely be established, that is, velocity profiles are nearly flat. As the Womersley number increases, the phase lag between flow and pressure approaches 90 degrees.

References

Books

• Truskey, G.A., Yuan, F. and Katz, D.F. (2009). Transport Phenomena in Biological Systems. 2nd edition. Pearson.
• White, F.M. and Majdalani, J. (2022). Viscous Fluid Flow. 4th edition. McGraw-Hill.
• Zamir, M. (2000). The Physics of Pulsatile Flow. Springer.

Papers

• Sexl, T. (1930). Über den von EG Richardson entdeckten “Annulareffekt“. Zeitschrift für Physik, 61, 349–362.
• Uchida, S. (1956). The pulsating viscous flow superimposed on the steady laminar motion of incompressible fluid in a circular pipe. Zeitschrift für angewandte Mathematik und Physik, 7, 377.
• Womersley, J.R. (1955). XXIV. Oscillatory motion of a viscous liquid in a thin-walled elastic tube – I: The linear approximation for long waves. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 46(373), 199–221.