pulsatile flow

{{Short description|A flow with periodic variations (fluid dyamics)}}

In fluid dynamics, a flow with periodic variations is known as pulsatile flow, or as Womersley flow. The flow profiles was first derived by John R. Womersley (1907–1958) in his work with blood flow in arteries.{{Cite journal

|author=Womersley, J.R.

|title=Method for the calculation of velocity, rate of flow and viscous drag in arteries when the pressure gradient is known

|journal=J. Physiol.

|volume=127

|issue=3

|pages=553–563

|date=March 1955

|pmid=14368548

|pmc=1365740

|doi=10.1113/jphysiol.1955.sp005276

}} The cardiovascular system of chordate animals is a very good example where pulsatile flow is found, but pulsatile flow is also observed in engines and hydraulic systems, as a result of rotating mechanisms pumping the fluid.

Equation

File:Womersley Flow Profile.svg

The pulsatile flow profile is given in a straight pipe by

: u(r, t) = \mathrm{Re}\left\{ \sum^N_{n=0} \frac{i \, P'_n}{\rho \, n \, \omega} \left[ 1 - \frac{J_0(\alpha \, n^{1/2} \, i^{3/2} \, \frac{r}{R})}{J_0(\alpha \, n^{1/2} \, i^{3/2})} \right] e^{i n \omega t} \right\} \, ,

where:

:

{{math|u}}is the longitudinal flow velocity,
{{math|r}}is the radial coordinate,
{{math|t}}is time,
{{math|α}}is the dimensionless Womersley number,
{{math|ω}}is the angular frequency of the first harmonic of a Fourier series of an oscillatory pressure gradient,
{{math|n}}are the natural numbers,
{{math|P'n}}is the pressure gradient magnitude for the frequency {{math|nω}},
{{math|ρ}}is the fluid density,
{{math|μ}}is the dynamic viscosity,
{{math|R}}is the pipe radius,
{{math|J0(·)}}is the Bessel function of first kind and order zero,
{{math|i}}is the imaginary number, and
{{math|Re{·}}}is the real part of a complex number.

Properties

= Womersley number =

The pulsatile flow profile changes its shape depending on the Womersley number

:\alpha = R \left( \frac{\omega \rho}{\mu} \right)^{1/2} \,.

For \alpha \lesssim 2, viscous forces dominate the flow, and the pulse is considered quasi-static with a parabolic profile.

For \alpha \gtrsim 2, the inertial forces are dominant in the central core, whereas viscous forces dominate near the boundary layer. Thus, the velocity profile gets flattened, and phase between the pressure and velocity waves gets shifted towards the core.{{cn|date=March 2021}}

= Function limits =

== Lower limit ==

The Bessel function at its lower limit becomes

{{cite web

| url =http://wwwf.imperial.ac.uk/~ajm8/BioFluids/lec1114.pdf

| title =Pulsatile flow in a long straight artery

| last =Mestel

| first =Jonathan

| date =March 2009

| website =

| publisher =Imperial College London

| access-date =6 January 2017

| quote = Bio Fluid Mechanics: Lecture 14

}}

: \lim_{z\to\infty} J_0(z) = 1 - \frac{z^2}{4} \,,

which converges to the Hagen-Poiseuille flow profile for steady flow for

: \lim_{n \to 0} u(r, t) = - \frac{P'_0}{4 \mu} \left(R^2 - r^2 \right) \, ,

or to a quasi-static pulse with parabolic profile when

: \lim_{\alpha \to 0} u(r, t) = \mathrm{Re} \left\{ - \sum^N_{n=0} \frac{P'_n}{4 \mu} (R^2 - r^2 ) \, e^{i n \omega t} \right\} = - \sum^N_{n=0} \frac{P'_n}{4 \mu} (R^2 - r^2 ) \, \cos(n \omega t) \, .

In this case, the function is real, because the pressure and velocity waves are in phase.

== Upper limit ==

The Bessel function at its upper limit becomes

: \lim_{z\to\infty} J_0(z \, i) = \frac{e^z}{\sqrt{2\pi\, z}} \,,

which converges to

: \lim_{z\to\infty} u(r, t) = \mathrm{Re} \left\{ \sum^N_{n=0} \frac{i \, P'_n}{\rho \, n \, \omega} \left[ 1 - e^{\alpha \, n^{1/2} \, i^{1/2} \left(\frac{r}{R} - 1 \right) } \right] e^{i n \omega t} \right\} = - \sum^N_{n=0} \frac{\, P'_n}{\rho \, n \, \omega} \left[ 1 - e^{\alpha \, n^{1/2} \left(\frac{r}{R} - 1 \right) } \right] \sin(n\,\omega\, t) \, .

This is highly reminiscent of the Stokes layer on an oscillating flat plate, or the skin-depth penetration of an alternating magnetic field into an electrical conductor.

On the surface u(r=R,t) = 0, but the exponential term becomes negligible once \alpha (1 - r/R) becomes large, the velocity profile becomes almost constant and independent of the viscosity. Thus, the flow simply oscillates as a plug profile in time according to the pressure gradient,

: \rho \frac{\partial u}{\partial t} = - \sum^N_{n=0} P'_n \,.

However, close to the walls, in a layer of thickness \mathcal{O}(\alpha^{-1}), the velocity adjusts rapidly to zero. Furthermore, the phase of the time oscillation varies quickly with position across the layer. The exponential decay of the higher frequencies is faster.

Derivation

For deriving the analytical solution of this non-stationary flow velocity profile, the following assumptions are taken:{{cite book

|author=Fung, Y. C.

|title=Biomechanics – Motion, flow, stress and growth

|year=1990

|page=569

|publisher=Springer-Verlag

|place=New York (USA)

|url=https://books.google.com/books?id=33qbOEKAWIwC&q=Biomechanics+-+Motion%2C+flow%2C+stress+and+growth|isbn=9780387971247

}}

{{cite journal

|last=Nield

|first=D.A.

|author2=Kuznetsov, A.V.

|title=Forced convection with laminar pulsating flow in a channel or tube

|journal=International Journal of Thermal Sciences

|year=2007

|volume=46

|issue=6

|pages=551–560

|doi=10.1016/j.ijthermalsci.2006.07.011|bibcode=2007IJTS...46..551N

}}

Thus, the Navier-Stokes equation and the continuity equation are simplified as

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

and

: \frac{\partial u}{\partial x} = 0 \, ,

respectively. The pressure gradient driving the pulsatile flow is decomposed in Fourier series,

: \frac{\partial p}{\partial x} (t) = \sum^N_{n=0}P'_n e^{in\omega t} \, ,

where i is the imaginary number, \omega is the angular frequency of the first harmonic (i.e., n = 1), and P'_n are the amplitudes of each harmonic n. Note that, P'_0 (standing for n = 0) is the steady-state pressure gradient, whose sign is opposed to the steady-state velocity (i.e., a negative pressure gradient yields positive flow). Similarly, the velocity profile is also decomposed in Fourier series in phase with the pressure gradient, because the fluid is incompressible,

: u(r,t) = \sum^N_{n=0}U_n e^{in\omega t} \, ,

where U_n are the amplitudes of each harmonic of the periodic function, and the steady component (n = 0) is simply Poiseuille flow

: U_0 = - \frac{P'_0}{4 \mu} \left(R^2 - r^2 \right) \, .

Thus, the Navier-Stokes equation for each harmonic reads as

: i\rho n\omega U_n = -P'_n +\mu \left(\frac{\partial^2 U_n}{\partial r^2} + \frac{1}{r} \frac{\partial U_n}{\partial r}\right) \, .

With the boundary conditions satisfied, the general solution of this ordinary differential equation for the oscillatory part (n \geq 1) is

: U_n(r) = A_n \, J_0 \left( \alpha \, \frac{r}{R} n^{1/2}\,i^{3/2} \right) + B_n \, Y_0 \left( \alpha \, \frac{r}{R} n^{1/2}\,i^{3/2} \right) + \frac{i\,P'_n}{\rho\, n\, \omega}\, ,

where J_0(\cdot) is the Bessel function of first kind and order zero, Y_0(\cdot) is the Bessel function of second kind and order zero, A_n and B_n are arbitrary constants, and \alpha = R \surd( \omega \rho / \mu ) is the dimensionless Womersley number. The axisymmetric boundary condition (\partial U_n/ \partial r|_{r=0} = 0) is applied to show that B_n = 0 for the derivative of above equation to be valid, as the derivatives J_0' and Y_0' approach infinity. Next, the wall non-slip boundary condition (U_n(R) = 0) yields A_n = - \frac{i \, P'_n}{\rho \, n \, \omega} \frac{1}{ J_0 \left( \alpha \, n^{1/2} \, i^{3/2} \right)}. Hence, the amplitudes of the velocity profile of the harmonic n becomes

: U_n(r) = \frac{i\,P'_n}{\rho \, n \, \omega} \left[ 1 - \frac{J_0(\alpha \, n^{1/2} \, i^{3/2} \, \frac{r}{R})}{J_0(\alpha \, n^{1/2} \, i^{3/2})} \right] = \frac{i\,P'_n}{\rho \, n \, \omega} \left[ 1 - \frac{J_0(\Lambda_n \, \frac{r}{R})}{J_0(\Lambda_n)} \right] \, ,

where \Lambda_n = \alpha \, n^{1/2} \, i^{3/2} is used for simplification.

The velocity profile itself is obtained by taking the real part of the complex function resulted from the summation of all harmonics of the pulse,

: u(r, t) = \frac{P'_0}{4 \mu} \left(R^2 - r^2 \right) + \mathrm{Re} \left\{ \sum^N_{n=1} \frac{i \, P'_n}{\rho \, n \, \omega} \left[ 1 - \frac{J_0(\Lambda_n \, \frac{r}{R})}{J_0(\Lambda_n)} \right] e^{i n \omega t} \right\} \, .

= Flow rate =

Flow rate is obtained by integrating the velocity field on the cross-section. Since,

: \frac{d}{dx} \left[ x^p J_p(a\,x) \right] = a\,x^p J_{p-1} (a\, x)

\quad \Rightarrow \quad

\frac{d}{dx} \left[ x\, J_1(a\,x) \right] = a\,x J_{0} (a\, x) \, ,

then

: Q(t) = \iint u(r, t) \, dA = \mathrm{Re} \left\{ \pi\,R^2 \, \sum^N_{n=1} \frac{i \, P'_n}{\rho \, n \, \omega} \left[ 1 - \frac{2}{\Lambda_n}\frac{J_1(\Lambda_n)}{J_0(\Lambda_n)} \right] e^{i n \omega t} \right\} \, .

=Velocity profile=

File:Comparison of velocity profile of Womersley flow.svg

To compare the shape of the velocity profile, it can be assumed that

:

u(r,t) = f(r)\,\frac{Q(t)}{A} \, ,

where

:

f(r) = \frac{u(r,t)}{\frac{Q(t)}{A}} = \mathrm{Re} \left\{ \sum^N_{n=1} \left[ \frac{\Lambda_n \, J_0(\Lambda_n) - \Lambda_n \, J_0(\Lambda_n \, \frac{r}{R})}{\Lambda_n \, J_0(\Lambda_n) - 2\,J_1(\Lambda_n)} \right] \right\}

is the shape function.

{{cite journal

|author=San, Omer

|author2=Staples, Anne E

|title=An improved model for reduced-order physiological fluid flows

|journal=Journal of Mechanics in Medicine and Biology

|year=2012

|volume=12

|issue=3

|pages=125–152

|doi=10.1142/S0219519411004666|arxiv=1212.0188

|s2cid=118525588

}}

It is important to notice that this formulation ignores the inertial effects. The velocity profile approximates a parabolic profile or a plug profile, for low or high Womersley numbers, respectively.

= Wall shear stress =

For straight pipes, wall shear stress is

:\tau_w = \mu \left. \frac{\partial u}{\partial r} \right|_{r=R} \, .

The derivative of a Bessel function is

:\frac{\partial}{\partial x}\left[ x^{-p} J_{-p} (a\,x) \right] = a\,x^{-p} J_{p+1}(a\,x)

\quad \Rightarrow \quad

\frac{\partial}{\partial x}\left[ J_0 (a\,x) \right] = -a\,J_1(a\,x) \, .

Hence,

:\tau_w = \mathrm{Re} \left\{ \sum^N_{n=1} P'_n \frac{R}{\Lambda_n} \frac{J_1(\Lambda_n)}{J_0(\Lambda_n)} e^{i n \omega t} \right\} \, .

= Centre line velocity =

If the pressure gradient P'_n is not measured, it can still be obtained by measuring the velocity at the centre line. The measured velocity has only the real part of the full expression in the form of

: \tilde{u}(t) = \mathrm{Re}(u(0, t)) \equiv \sum^N_{n=1} \tilde{U}_n \, \cos(n \, \omega \, t) \, .

Noting that J_0(0) = 1, the full physical expression becomes

: u(0, t) = \mathrm{Re} \left\{ \sum^N_{n=1} \frac{i \, P'_n}{\rho \, n \, \omega} \left[ \frac{J_0(\Lambda_n) - 1}{J_0(\Lambda_n)} \right] e^{i n \omega t} \right\}

at the centre line. The measured velocity is compared with the full expression by applying some properties of complex number. For any product of complex numbers (C = AB), the amplitude and phase have the relations |C| = |A||B| and \phi_C = \phi_A + \phi_B, respectively. Hence,

:

\tilde{U}_n = \left| \frac{i \, P'_n}{\rho \, n \, \omega} \left[ \frac{J_0(\Lambda_n) - 1}{J_0(\Lambda_n)} \right] \right|

\quad \Rightarrow \quad

P'_n = \tilde{U}_n \left| i \, \rho \, n \, \omega \left[ \frac{J_0(\Lambda_n)}{1 - J_0(\Lambda_n)} \right] \right|

and

:

\tilde{\phi} = 0 = \phi_{P'_n} + \phi_{U_n}

\quad \Rightarrow \quad

\phi_{P'_n} = \operatorname{phase} \left( \frac{i}{\rho \, n \, \omega} \left[ \frac{1 - J_0(\Lambda_n)}{J_0(\Lambda_n)} \right] \right) \, ,

which finally yield

:

\frac{1}{\rho} \frac{\partial p}{\partial x} = \sum^N_{n=1} \tilde{U}_n \left| i \, \rho \, n \, \omega \left[ \frac{J_0(\Lambda_n)}{1 - J_0(\Lambda_n)} \right] \right| \, \cos \left\{ n\,\omega\,t + \operatorname{phase} \left( \frac{i}{\rho \, n \, \omega} \left[ \frac{1 - J_0(\Lambda_n)}{J_0(\Lambda_n)} \right] \right) \right\} \, .

See also

References