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. 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 The pulsatile flow profile is given in a straight pipe by
u ( r , t ) = R e { ∑ n = 0 N i P n ′ ρ n ω [ 1 − J 0 ( α n 1 / 2 i 3 / 2 r R ) J 0 ( α n 1 / 2 i 3 / 2 ) ] e i n ω t } , {\displaystyle u(r,t)=\mathrm {Re} \left\{\sum _{n=0}^{N}{\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^{in\omega t}\right\}\,,}
where:
Properties
Womersley number The pulsatile flow profile changes its shape depending on the Womersley number
α = R ( ω ρ μ ) 1 / 2 . {\displaystyle \alpha =R\left({\frac {\omega \rho }{\mu }}\right)^{1/2}.}
For α ≲ 2 {\displaystyle \alpha \lesssim 2} , viscous forces dominate the flow, and the pulse is considered quasi-static with a parabolic profile. For α ≳ 2 {\displaystyle \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.
Function limits
Lower limit The Bessel function at its lower limit becomes
lim z → ∞ J 0 ( z ) = 1 − z 2 4 , {\displaystyle \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 → 0 u ( r , t ) = − P 0 ′ 4 μ ( R 2 − r 2 ) , {\displaystyle \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 α → 0 u ( r , t ) = R e { − ∑ n = 0 N P n ′ 4 μ ( R 2 − r 2 ) e i n ω t } = − ∑ n = 0 N P n ′ 4 μ ( R 2 − r 2 ) cos ( n ω t ) . {\displaystyle \lim _{\alpha \to 0}u(r,t)=\mathrm {Re} \left\{-\sum _{n=0}^{N}{\frac {P'_{n}}{4\mu }}(R^{2}-r^{2})\,e^{in\omega t}\right\}=-\sum _{n=0}^{N}{\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 → ∞ J 0 ( z i ) = e z 2 π z , {\displaystyle \lim _{z\to \infty }J_{0}(z\,i)={\frac {e^{z}}{\sqrt {2\pi \,z}}}\,,}
which converges to
lim z → ∞ u ( r , t ) = R e { ∑ n = 0 N i P n ′ ρ n ω [ 1 − e α n 1 / 2 i 1 / 2 ( r R − 1 ) ] e i n ω t } = − ∑ n = 0 N P n ′ ρ n ω [ 1 − e α n 1 / 2 ( r R − 1 ) ] sin ( n ω t ) . {\displaystyle \lim _{z\to \infty }u(r,t)=\mathrm {Re} \left\{\sum _{n=0}^{N}{\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^{in\omega t}\right\}=-\sum _{n=0}^{N}{\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 {\displaystyle u(r=R,t)=0} , but the exponential term becomes negligible once α ( 1 − r / R ) {\displaystyle \alpha (1-r/R)} becomes large, and 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
ρ ∂ u ∂ t = − ∑ n = 0 N P n ′ . {\displaystyle \rho {\frac {\partial u}{\partial t}}=-\sum _{n=0}^{N}P'_{n}\,.}
However, close to the walls, in a layer of thickness O ( α − 1 ) {\displaystyle {\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:
The fluid is homogeneous, incompressible and Newtonian; The tube wall is rigid and circular; Motion is laminar, axisymmetric, and parallel to the tube's axis; The boundary conditions are axisymmetry at the centre, and no-slip on the wall; The pressure gradient is a periodic function that drives the fluid; and Gravitation has no effect on the fluid. Thus, the Navier-Stokes equation and the continuity equation are simplified as
ρ ∂ u ∂ t = − ∂ p ∂ x + μ ( ∂ 2 u ∂ r 2 + 1 r ∂ u ∂ r ) {\displaystyle \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
∂ u ∂ x = 0 , {\displaystyle {\frac {\partial u}{\partial x}}=0,}
respectively. The pressure gradient driving the pulsatile flow is decomposed in Fourier series,
∂ p ∂ x ( t ) = ∑ n = 0 N P n ′ e i n ω t , {\displaystyle {\frac {\partial p}{\partial x}}(t)=\sum _{n=0}^{N}P'_{n}e^{in\omega t},}
where i {\displaystyle i} is the imaginary number, ω {\displaystyle \omega } is the angular frequency of the first harmonic (i.e., n = 1 {\displaystyle n=1} ), and P n ′ {\displaystyle P'_{n}} are the amplitudes of each harmonic n {\displaystyle n} . P 0 ′ {\displaystyle P'_{0}} (standing for n = 0 {\displaystyle 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 ) = ∑ n = 0 N U n e i n ω t , {\displaystyle u(r,t)=\sum _{n=0}^{N}U_{n}e^{in\omega t},}
where U n {\displaystyle U_{n}} are the amplitudes of each harmonic of the periodic function, and the steady component ( n = 0 {\displaystyle n=0} ) is simply Poiseuille flow
U 0 = − P 0 ′ 4 μ ( R 2 − r 2 ) . {\displaystyle U_{0}=-{\frac {P'_{0}}{4\mu }}\left(R^{2}-r^{2}\right).}
Thus, the Navier-Stokes equation for each harmonic reads as
i ρ n ω U n = − P n ′ + μ ( ∂ 2 U n ∂ r 2 + 1 r ∂ U n ∂ r ) . {\displaystyle 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 ≥ 1 {\displaystyle n\geq 1} ) is
U n ( r ) = A n J 0 ( α r R n 1 / 2 i 3 / 2 ) + B n Y 0 ( α r R n 1 / 2 i 3 / 2 ) + i P n ′ ρ n ω , {\displaystyle 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 ( ⋅ ) {\displaystyle J_{0}(\cdot )} is the Bessel function of first kind and order zero, Y 0 ( ⋅ ) {\displaystyle Y_{0}(\cdot )} is the Bessel function of second kind and order zero, A n {\displaystyle A_{n}} and B n {\displaystyle B_{n}} are arbitrary constants, and α = R √ ( ω ρ / μ ) {\displaystyle \alpha =R\surd (\omega \rho /\mu )} is the dimensionless Womersley number. The axisymmetric boundary condition ( ∂ U n / ∂ r | r = 0 = 0 {\displaystyle \partial U_{n}/\partial r|_{r=0}=0} ) is applied to show that B n = 0 {\displaystyle B_{n}=0} for the derivative of above equation to be valid, as the derivatives J 0 ′ {\displaystyle J_{0}'} and Y 0 ′ {\displaystyle Y_{0}'} approach infinity. Next, the wall non-slip boundary condition ( U n ( R ) = 0 {\displaystyle U_{n}(R)=0} ) yields A n = − i P n ′ ρ n ω 1 J 0 ( α n 1 / 2 i 3 / 2 ) . {\displaystyle 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 {\displaystyle n} become
U n ( r ) = i P n ′ ρ n ω [ 1 − J 0 ( α n 1 / 2 i 3 / 2 r R ) J 0 ( α n 1 / 2 i 3 / 2 ) ] = i P n ′ ρ n ω [ 1 − J 0 ( Λ n r R ) J 0 ( Λ n ) ] , {\displaystyle 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 Λ n = α n 1 / 2 i 3 / 2 {\displaystyle \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 ) = P 0 ′ 4 μ ( R 2 − r 2 ) + R e { ∑ n = 1 N i P n ′ ρ n ω [ 1 − J 0 ( Λ n r R ) J 0 ( Λ n ) ] e i n ω t } . {\displaystyle u(r,t)={\frac {P'_{0}}{4\mu }}\left(R^{2}-r^{2}\right)+\mathrm {Re} \left\{\sum _{n=1}^{N}{\frac {i\,P'_{n}}{\rho \,n\,\omega }}\left[1-{\frac {J_{0}(\Lambda _{n}\,{\frac {r}{R}})}{J_{0}(\Lambda _{n})}}\right]e^{in\omega t}\right\}.}
Flow rate Flow rate is obtained by integrating the velocity field on the cross-section. Since,
d d x [ x p J p ( a x ) ] = a x p J p − 1 ( a x ) ⇒ d d x [ x J 1 ( a x ) ] = a x J 0 ( a x ) , {\displaystyle {\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\,xJ_{0}(a\,x)\,,}
then
Q ( t ) = ∬ u ( r , t ) d A = R e { π R 2 ∑ n = 1 N i P n ′ ρ n ω [ 1 − 2 Λ n J 1 ( Λ n ) J 0 ( Λ n ) ] e i n ω t } . {\displaystyle Q(t)=\iint u(r,t)\,dA=\mathrm {Re} \left\{\pi \,R^{2}\,\sum _{n=1}^{N}{\frac {i\,P'_{n}}{\rho \,n\,\omega }}\left[1-{\frac {2}{\Lambda _{n}}}{\frac {J_{1}(\Lambda _{n})}{J_{0}(\Lambda _{n})}}\right]e^{in\omega t}\right\}.}
Velocity profile
To compare the shape of the velocity profile, it can be assumed that
u ( r , t ) = f ( r ) Q ( t ) A , {\displaystyle u(r,t)=f(r)\,{\frac {Q(t)}{A}}\,,}
where
f ( r ) = u ( r , t ) Q ( t ) A = R e { ∑ n = 1 N [ Λ