A Simple Transient Poiseuille-Based Approach to Mimic the Womersley Function and to Model Pulsatile Blood Flow

Abstract

As it is known, the Womersley function models velocity as a function of radius and time. It has been widely used to simulate the pulsatile blood flow through circular ducts. In this context, the present study is focused on the introduction of a simple function as an approximation of the Womersley function in order to evaluate its accuracy. This approximation consists of a simple quadratic function, suitable to be implemented in most commercial and non-commercial computational fluid dynamics codes, without the aid of external mathematical libraries. The Womersley function and the new function have been implemented here as boundary conditions in OpenFOAM ESI software (v.1906). The discrepancy between the obtained results proved to be within 0.7%, which fully validates the calculation approach implemented here. This approach is valid when a simplified analysis of the system is pointed out, in which flow reversals are not contemplated.

1. Introduction

The Womersley function [1] models the transient/pulsatile velocity profile of blood through circular ducts. It was derived as an exact solution of viscous flow equations through a circular tube subjected to a periodic pressure gradient. Since then, the Womersley function has been widely used in hemodynamics. In the context of computational fluid dynamics modelization of arteries and blood vessels, several works employ the Womersley function in the inlet velocity boundary condition (see [2,3,4,5,6,7,8,9]).

Additional theoretical studies concerned with CFD have been conducted on the characterization of this function’s parameters. Shehada et al. [10], to improve the understanding of the pulsatile flow’s nature, pointed out the idealized three-dimensional velocity profiles for both the normal and femoral carotid arteries, considering the Fourier harmonics for each shape. The wave and the respective velocity profiles at each instant of the time were calculated with the Womersley equations.

Loudon et al. [11] have outlined the category of problems specific to internal flow by considering the variation of the dimensionless Womersley number, Wo (defined in Section 2), a very important parameter in the simulation of blood flow. The Womersley number is a dimensionless parameter that indicates the relationship between the pulse flow frequency and viscosity in biofluid mechanics. The exact analytical solution for transient flow between two parallel walls is based on the same fluid behavior pattern identified in the flow within cylinders. It has been demonstrated that when Wo < 1, the flow is expected to follow the oscillating pressure gradient faithfully, and the velocity profiles exhibit a parabolic shape such that the fluid oscillating with maximum amplitude is farther from the walls (“almost stable” behavior). When Wo > 1, the velocity profiles are no longer parabolic, and the flow is out of phase in time with respect to the oscillating pressure gradient. The amplitude of the oscillating fluid can increase or decrease as Wo > 1. Additional important studies regarding the Womersley number in pulsatile blood flow are recounted in [12,13,14,15].

It is worth mentioning that the Bessel functions of imaginary numbers represent the Womersley function. This work is aimed at the substitution of the above-mentioned Womersley function with the Poiseuille parabolic profile function [16], the latter being much simpler from the mathematical point of view, allowing its implementation without external mathematical libraries. The applicability and implementation difficulty of the Womersley function is widely known in the literature [17,18].

2. Womersley and Poiseuille Velocity Profile for a Pulsatile Blood Flow

Consider the fluid motion inside a tube of diameter D, subjected to a periodic pressure difference, P_in − P_out = Acos(ω t) (where A is the amplitude, ω is the angular frequency, and t is the time).

The solution to this problem is known as the Womersley velocity profile [1]. Several authors (i.e., [5]) adapt this function to the measured pulsed flow by writing the following function:

$$v_W\left(\frac{2r}{D},t\right)=\Re \left\{Wo\sqrt{i^3}\left[\frac{J_0\left(Wo\sqrt{i^3}\frac{2r}{D}\right)-J_0\left(Wo\sqrt{i^3}\right)}{2J_0\left(Wo\sqrt{i^3}\right)-Wo\sqrt{i^3}{J}_0\left(Wo\sqrt{i^3}\right)}\right]\frac{Q\left(t\right)}{Q_0}\right\},$$

(1)

where R{⋅} denotes the real part of the function defined in the complex plane, i = (−1)^0.5 is the imaginary unit, J₀(.) is the zero order modified Bessel function, (2r/D) is the dimensionless variable in which D is the tube diameter, r is the distance from the tube center line, Wo = 0.5D$\sqrt{\left(\omega \rho /\eta \right)}$ represents the Womersley number in which ω = 2π/T is the angular frequency determined from characteristic period T, ρ is the fluid density, η is the dynamic viscosity, Q(t) is the flow rate variable in time, and Q₀ is the average flow rate.

If in the same problem the inlet pressure difference is kept constant, the result is the Poiseuille velocity profile, illustrated in [11]:

$$v_P\left(\frac{2r}{D}\right)=\gamma \left[1-{\left(\frac{2r}{D}\right)}^2\right],$$

(2)

where γ denotes the non-dimensional amplitude.

Finally, it is worth mentioning that both velocity profiles refer to laminar and steady-state flow conditions.

3. Approximation Method and Validation Model

The Womersley velocity profile (Equation (1)) is a function of r and t, while the Poiseuille velocity profile (Equation (2)) is a function of r. In order to compare the mentioned two curves, the Womersley profile was considered in the following version:

$$v_W\left(\frac{2r}{D},t\right)=f\left(\frac{2r}{D}\right)g\left(t\right),$$

(3)

where the function f(2r/D) is composed from the real part of Equation (1) and g(t) is equal to R(Q(t)/Q₀) (the remaining term of Equation (1)). In this modality, it is now possible to compare Equation (2) with the function f(2r/D) defined in Equation (3).

Therefore, the problem is now represented by the equivalence:

$$\Re \left\{Wo\sqrt{i^3}\left[\frac{J_0\left(Wo\sqrt{i^3}\frac{2r}{D}\right)-J_0\left(Wo\sqrt{i^3}\right)}{2J_0\left(Wo\sqrt{i^3}\right)-Wo\sqrt{i^3}{J}_0\left(Wo\sqrt{i^3}\right)}\right]\right\}=\gamma \left[1-{\left(\frac{2r}{D}\right)}^2\right],$$

(4)

γ denotes the unknown variable, which will be a function of (2R/D). From Equation (4), the following calculation is obtained:

$$\gamma \left(\frac{2r}{D}\right)=\frac{\Re \left\{Wo\sqrt{i^3}\left[\frac{J_0\left(Wo\sqrt{i^3}\frac{2r}{D}\right)-J_0\left(Wo\sqrt{i^3}\right)}{2J_0\left(Wo\sqrt{i^3}\right)-Wo\sqrt{i^3}{J}_0\left(Wo\sqrt{i^3}\right)}\right]\right\}}{\left[1-{\left(\frac{2r}{D}\right)}^2\right]}.$$

(5)

In order to validate the approach pointed out here, typical values of the blood flow inside an artery were considered (data are taken from Vimmr et al. [5]): D = 0.003 m, blood density ρ = 1060 kg/m³, blood dynamic viscosity η = 3.45 × 10⁻³ Pa s and cardiac cycle period T = 1.68 s.

The plot of Equation (5) is represented in Figure 1 together with the average value of this calculated according to the formula:

$$M\left(\gamma \left(\frac{2r}{D}\right),\left[a,b\right]\right)=\frac{1}{b-a}{{\displaystyle \int }}_a^b\gamma \left(\frac{2r}{D}\right)d\left(\frac{2r}{D}\right),$$

(6)

where [a, b] = [0, 1] is the integral range. The average value M of Equation (6) is equal to 1.99923.

The calculated γ average value can be used as an approximate amplitude value in the Poiseuille Equation (2). Therefore, the function R{⋅} has been approximated in a transient Poiseuille velocity profile. For the sake of clarity, the “transient” Poiseuille velocity profile has been referred to as the Poiseuille velocity profile multiplied by the function that adjusts the amplitude over time; that is, the following equation is intended:

$$v_P\left(\frac{2r}{D},t\right)=\gamma \left[1-{\left(\frac{2r}{D}\right)}^2\right]g\left(t\right),$$

(7)

where g(t) is defined in Equation (3).

What varies during the transient flow is, therefore, the function g(t), which has an influence on the width of the profile.

The relative percentage error is evaluated according to the relation:

$$\epsilon \%=\frac{\left|f\left(\frac{2r}{D}\right)-v_P\left(\frac{2r}{D}\right)\right|}{f\left(\frac{2r}{D}\right)}100.$$

(8)

Figure 2 shows the comparison between function f(2r/D) and v_P(2r/D) and the corresponding percentage error along the dimensionless axis. In addition, all the corresponding numerical values have been reported in

Table 1. Numerical value of the Womersley f(2r/D) profile and Poiseuille profile considering the example values: ρ = 1060 kg/m³, η = 3.45 × 10⁻³ Pa s, D = 0.003 m, T = 1.68 s, and Wo = 1.61.

2r/D	f(2r/D)	vP(2r/D)	ε%	2r/D	f(2r/D)	vP(2r/D)	ε%
0.00	1.98845	1.99923	0.542061	0.55	1.39521	1.39446	0.053681
0.05	1.98360	1.99423	0.536241	0.60	1.28133	1.27951	0.142441
0.10	1.96902	1.97924	0.518870	0.65	1.15722	1.15456	0.230149
0.15	1.94472	1.95425	0.490221	0.70	1.02283	1.01961	0.314656
0.20	1.91065	1.91926	0.450745	0.75	0.87812	0.87466	0.393649
0.25	1.86679	1.87428	0.401073	0.80	0.72308	0.71972	0.464652
0.30	1.81310	1.81930	0.342012	0.85	0.55772	0.55479	0.525025
0.35	1.74952	1.75433	0.274547	0.90	0.38204	0.37985	0.571956
0.40	1.67601	1.67936	0.199834	0.95	0.19611	0.19493	0.602458
0.45	1.59249	1.59439	0.119200	1.00	0.00000	0.00000	-
0.50	1.49891	1.49942	0.034141

. The maximum approximation error is slightly greater than 0.60% in the boundary points. This discrepancy validates the calculation approach developed here.

The Womersley number Wo, defined in Equation (1), is a very important parameter that affects the validity of the approximation explained above. As shown in Figure 3, as the Womersley number increases, the behavior of the Womersley function moves further and further away from the parabolic profile.

Due to this, the approximation of the Womersley profile with the Poiseuille profile is valid for little values of the Womersley number Wo (Wo < 2).

It is worth mentioning that the limit of this approximation mainly relays on the neglection of possible flow inversions. More specifically, the above-mentioned approximation is valid for code implementation, when an input flow value is imposed or a fully developed parabolic velocity profile is considered, provided that Wo < 2.

4. Case Study

Vimmr et al. [5] have applied the Womersley function as a velocity input condition for the study of a bypass in different geometric conditions. This study’s peculiarity was to adapt the Womersley function to the specific flow conditions of the heartbeat.

A time-dependent inlet flow rate Q(t) was considered that corresponds to flow rate values measured in the right coronary artery during rest. The flow rate is prescribed in the form of the following Fourier series [5]:

$$Q\left(t\right)=Q_0+{\displaystyle \sum }_{k=1}^5{Q}_k\cos \left(k\omega t-{\phi }_k\right),$$

(9)

where the cardiac cycle period is T = 1.68 s, Q₀ = 65.07 mL/min represents the average inlet flow rate, and Q_k, and φ_k, k = 1, …, 5 are the amplitude and phase angle, respectively. The values used for Equation (9) are coherent with the ones found by Vimmr et al. [5]: Q₁ = 18.149 mL/min, Q₂ = 34.828 mL/min, Q₃ = 12.329 mL/min, Q₄ = 9.107 mL/min, Q₅ = 2.944 mL/min, φ₁ = 1.944 rad, φ₂ = 2.836 rad, φ₃ = −2.124 rad, φ₄ = −1.875 rad and φ₅ = −0.447 rad.

Thus, the Womersley function used by Vimmr et al. [5] to define the velocity shape inlet was the same as Equation (1). The same properties defined in Section 3 are reported: Newtonian blood with density ρ = 1060 kg/m³ and dynamic viscosity η = 3.45 × 10⁻³ Pa s, artery diameter D = 0.003 m and Wo = 1.61.

The approximation of Equation (1), according to the Poiseuille profile (Equation (2)) is characterized by a γ value equal to 1.99923, as shown in Section 3.

To compare the difference between the profiles, two simulations using OpenFOAM ESI (v.1906) have been performed: the first considering the Womersley profile according to Equation (1) as the velocity input, the second considering the transient Poiseuille equation with a γ value equal to 1.99923.

A rigid-walled cylinder with length L = 25D and pressure outlet P_out equal to 0 Pa was considered for both simulations.

To reduce the computational effort and the simulation time without losing numerical accuracy, an axially anisotropic with constant spacing mesh was adopted. A prism layer along the wall was employed to capture the boundary layer’s steepest effects. The resulting mesh consists of 81,600 hexahedral cells, as depicted in Figure 4.

Second-order numerical schemes for spatial and temporal discretization have been selected. For the velocity field, a preconditioned bi-conjugate gradient (PBiCG) coupled solver algorithm has been preferred, due to its numerical stability and convergence ratio. The pressure field has been solved using a generalized algebraic multi-grid (GAMG) algorithm with a Gauss–Seidel smoother.

Figure 5 shows the velocity comparisons between the Womersley velocity profile (shown in solid line) and the Poiseuille function profile (indicated with dots) in different sections (denoted with A, B, C, D, and E) for different times relating to the cardiac cycle defined with Equation (8). Based on pure observation, the results clearly show that there is no evident difference between the two studied velocity profiles.

5. Conclusions

In fluid flows inside blood vessels, or in general, in laminar conditions inside channels of relatively small sections, the motion generated by pressure variation over time of the wave type induces a velocity profile, which has been modelled in the literature as the Womersley function. This function has proved to be difficult to be implemented numerically. The numerical investigation conducted here demonstrates that a simplification by means of the Poiseuille function does not generate detected errors on the speed profile since the results agree approximately within 0.6%. Furthermore, the choice to substitute the Womersley equation with the Poiseuille function is advantageous even if only one unknown variable must be determined. Finally, since the non-invasive estimation of arterial blood volume flow (BVF) has become a central issue in the assessment of cardiovascular risk, it is possible to affirm that, based on the results obtained here, the Poiseuille and Womersley approaches are practically equivalent and can be commonly used to assess the BVF from the centerline velocity. Finally, it is important to specify that the range of applicability of the model developed here is for Womersley number values not exceeding 2.