Computational Modeling & Fluid Mechanics Universidad Carlos III de Madrid Bessel Navier-Stokes Solver

Pulsatile Hemodynamics & Womersley Flow Simulation

Numerical and analytical resolution of the unsteady Navier-Stokes equations for oscillatory blood flow in rigid cylindrical geometries. Evaluates how the Womersley dimensionless number (α = R√(ω/ν)) governs the transition from classical parabolic Poiseuille flow to blunted plug flow and high-frequency annular velocity overshoot (Richardson annular effect), validated against benchtop pulsatile flow experimental tracking.

Author & Simulation Lead: Álvaro García Casas Institution: Universidad Carlos III de Madrid

MATLAB R2024b Complex Bessel (J₀) Navier-Stokes Solvers Fourier Spectral Decomposition Particle Tracking Validation

Dynamic Velocity Profile Monitor • u(r, t)

Real-time solution of the unsteady momentum equation. Observe smooth profile morphing and wall shear deceleration across pulsatile cardiac phases.

Phase: 0°

α = 1

u(r, t) Radial Distribution • Boundary Layer δ

Peak Center Velocity: 1.00 u_max

[α = 1.0 — VISCOUS PARABOLIC]: Quasi-steady flow. Viscous diffusion acts rapidly across the lumen. Centerline velocity peaks at 1.00 with perfect convex parabolic curvature. Zero phase lag with driving pressure gradient. Characteristic of microvascular arterioles.

α ≤ 1.0

Parabolic (Poiseuille)

Viscous-dominated, in-phase flow

α ≈ 4–7

Plug-like Transition

Smooth shoulders, core blunting

α ≈ 14

Inertial Depression

Center dips to 0.65, dual peaks

α ≥ 20

Annular Overshoot

Deep bowl (0.52), wall peaks (0.89)

Theoretical Hydrodynamics

Fluid Mechanics: Transient Inertia vs. Viscous Shear

In cardiovascular hemodynamics, the heart drives flow via an oscillatory pressure gradient. The Womersley number (α) is the square root of the ratio between oscillatory inertial forces and viscous forces. It dictates how deeply wall friction penetrates into the vascular lumen over one cardiac cycle.

Womersley Parameter Formulation

α = R ⋅ √(ω / ν) = R ⋅ √(2πf ρ / μ)

R = vessel radius • ω = cardiac frequency • ν = kinematic blood viscosity

When α is small (< 2), momentum diffuses across the entire radius in a fraction of a cycle, preserving the classic parabolic Poiseuille shape. When α > 10, the core fluid has too much inertia to respond synchronously to the pressure gradient, compressing viscous dissipation into a thin Stokes boundary layer of thickness δ ∼ R/α.

Axisymmetric Navier-Stokes Equation

ρ ⋅ ∂u/∂t = −∂P/∂z + μ ⋅ (1/r) ⋅ ∂/∂r(r ⋅ ∂u/∂r)

Exact analytical solution via zeroth-order Bessel functions J₀(i^3/2α ⋅ r/R)

Because Bessel functions of complex arguments introduce real and imaginary components, the resulting velocity u(r, t) undergoes a continuous phase lag across the radius. Near the vessel wall, low inertia allows the fluid to reverse direction almost immediately when the pressure gradient inverts, while the core continues gliding forward.

Numerical Methodology & Benchtop Validation

Computational Pipeline & Experimental Validation

Resolution of the Womersley differential system in MATLAB using complex Bessel evaluation (besselj), compared directly against experimental particle tracking in physical pulsatile flow channels.

Empirical Validation

Physical Pulsatile Channel Tracking vs. Numerical Profiles

Laboratory Particle Tracking DAQ

Figure A • Empirical Profile at Wo = 7.09 (f = 0.5 Hz): Discrete particle image displacement mapped across a ±4 mm tube diameter. Confirms smooth plug-like core blunting with steep velocity gradient within 1 mm of the wall.

Figure B • Empirical Profile at Wo = 14.18 (f = 2.0 Hz): Under elevated pulsation frequency, the experimental tracking clearly displays the characteristic velocity overshoot lobes near ±2.4 mm, validating the M-profile prediction.

Wall Shear Stress τ_w(t) Over Cardiac Cycle

MATLAB Analytical

As α escalates, the wall shear stress exhibits a pronounced phase advance (up to 45°) relative to the driving pressure gradient, resulting in elevated oscillatory shear index (OSI).

Stokes Boundary Layer Thickness δ/R vs α

δ ∼ R / α

The viscous penetration depth decays hyperbolically. For α > 10, over 90% of the luminal cross-section behaves as an inviscid inertial plug.

Centerline Phase Lag φ_center vs α

Phase Retardation

Fluid in the center requires substantial momentum transfer to accelerate. As α increases, centerline velocity lags the driving pressure wave by asymptotically approaching 90°.

Oscillatory Flow Rate Attenuation |Q| / |Q_Poiseuille|

Impedance Damping

Rapid oscillation prevents full volumetric acceleration, causing vascular impedance to increase sharply with heart rate and vessel radius.

Regime Metrics

Morphological & Hemodynamic Characterization

Womersley Number Spectrum • Theoretical Values

Analytical Reference

Womersley α	Center Value u(0)	Peak Value u_max	Peak Location r/R	Profile Curvature & Morphology	Anatomical Site
1.0	1.00	1.00	0.00 (Center)	Convex parabolic (Poiseuille)	Terminal arterioles
4.0	0.88	0.88	0.00 (Center)	Convex arch with broader shoulders	Radial / Tibial arteries
7.0	0.755	0.755	0.00 (Center)	Plug-like: gentle plateau, rounded flanks	Femoral / Brachial arteries
14.0	0.645	0.760	±0.76	Concave center bowl, twin rounded lobes	Common Carotid artery
20.0	0.525	0.890	±0.84	Full M-profile (Richardson annular overshoot)	Ascending Aorta (resting/exercise)

Conclusions & Clinical Impact

Translational Significance in Vascular Biomechanics

Atherogenesis & Endothelial Shear Stress

The Richardson annular effect causes early flow deceleration and boundary layer reversal near the vessel wall while the core still moves forward. This creates low and oscillatory wall shear stress (OSI > 0.3), the proven mechanical trigger for focal atheromatous plaque formation in the carotid bifurcation.

Accuracy Limits of Poiseuille Assumptions

Clinical Doppler ultrasound algorithms that presume a parabolic Poiseuille velocity profile severely overestimate peak systolic volumetric flow in large conduits (α > 7) by up to 35%, proving the clinical necessity of Womersley correction filters in diagnostic hemodynamics.

Experimental Rig Validation

Physical benchtop tracking of tracer particles confirmed the exact numerical Bessel solution, demonstrating that non-dimensional Womersley scaling accurately predicts unsteady physiological velocity fields without requiring full 3D CFD discretization in straight conduits.

Author: Álvaro García Casas Institution: Universidad Carlos III de Madrid

← View All Projects