Effect of Pulsatile Flow Waveform and Womersley Number on the Flow in Stenosed Arterial Geometry

Abstract

The salient hemodynamic flow features in a stenosed artery depend not only on the degree of stenosis, but also on its location in the circulatory system and the physiological condition of the body. The nature of pulsatile flow waveform and local Womersley number vary in different regions of the arterial system and at different physiological state, which affects the local hemodynamic wall parameters, for example, the wall shear stress (WSS) and oscillatory shear index (OSI). Herein, we have numerically investigated the effects of different waveforms and Womersley numbers on the flow pattern and hemodynamic parameters in an axisymmetric stenosed arterial geometry with 50% diametral occlusion. Temporal evolution of the streamlines and hemodynamic parameters are investigated, and the time-averaged hemodynamic wall parameters are compared. Presence of the stenosis is found to increase the OSI of the flow even at the far-downstream side of the artery. At larger Womersley numbers, the instantaneous flow field in the stenosed region is found to have a stronger influence on the flow profiles of the previous time levels. The study delineates how an approximation in the assumption of inlet pulsatility profile may lead to significantly different prediction of hemodynamic wall parameters.

1. Introduction

The investigation of the flow through a stenosed geometry is of interest because of its significance in connection to vascular diseases, like atherosclerosis which perpetrates to heart attack and stroke [1]. High level of cholesterol or triglycerides in blood results in continuous deposition of the fatty substances along the inner lining of the artery wall, constricting the flow passage and hardening the wall [2]. Numerous investigations have demonstrated that the altered shear stress activates mitogenic signaling pathways, promotes vascular cell proliferation and migration [3], enhances lipid deposition [4], and facilitates neointimal formation and atherogenesis [5]. Wootton and Ku [6] postulated a strong inverse correlation between wall shear stress and the atherosclerotic intimal thickening. Nerem [7] and Chiu et al. [8] concluded that the endothelial cells that line the arteries exhibit high cell division rates and low cell density near the reattachment point downstream to the stenosis, where the wall shear stress is low but the wall shear stress gradient (WSSG) is high. On the other hand, the increased pressure losses associated with flow separation at the stenosis reduce the flow rate, which can produce ischemia (localized absence of oxygen), needing surgical interventions such as stenting, angioplasty, and bypass operations [9]. Furthermore, in the presence of stenoses, hydrodynamic forces are involved in plaque fissure [10], poststenotic dilatation [11], thrombosis, and damage of the arterial wall in the poststenotic region [12, 13]. All these factors underscore the importance of investigating the hemodynamics in stenosed arteries to predict the severity of the atherosclerotic formation facilitating the pathological interpretation of the disease.

Blood flow and pressure in the vasculature are unsteady due to the pulsatile flow conditions in the arteries. Berger and Jou [15] presented a detailed review of the experimental and numerical studies in stenosed arterial geometries. Experimental studies (mostly in vitro) in stenosed arteries focused attention on the disturbances created downstream to the stenosis. Ahmed and Giddens [16] used ultrasound and laser Doppler techniques to analyze the flow disturbances and wall shear stress associated with the presence of stenoses. The photochromic tracer method used by Ojha et al. [14] enabled the description of coherent flow structures present in the pulsatile poststenotic flow field. Siouffi et al. [17] and Deplano and Siouffi [18] used the pulsed Doppler ultrasonic velocimetry for the measurement of velocity and determination of wall shear stress around the zone of stenosis, while Brunette et al. [19] used particle image velocimetry to find out the wall shear stress distribution and secondary flow associated with a mild stenosis in a phantom coronary artery.

Presently, joint clinical-computational studies have been extensively used in circulatory pathologies to plan surgical actions [20–22]. The rheology of the blood, expressed by the constitutive equation, is an important fluid property that affects the computational flow prediction. Various studies have been conducted [23–27] assuming different non-Newtonian constitutive equations (e.g., Herschel-Berkley, power law, Quemada, Casson models) for blood. On the other hand, blood behaves as a Newtonian fluid at shear rates higher than 100 s⁻¹ [15, 28]. Accordingly, many researchers have analyzed the blood flow (particularly at larger Reynolds numbers) considering the fluid to be Newtonian [9, 29]. Johnston et al. [26] analyzed the effect of blood rheology in four different human right coronary arteries through a transient study. They found that the Newtonian and non-Newtonian blood viscosity models predict only a little variation in the wall shear stress in the arteries though the particle paths are often predicted differently.

Misra and coworkers [30, 31] presented approximate analytical solutions to pulsatile flow in mild stenoses, taking deformability of the wall into consideration following the thin-wall theory. Stergiopulos et al. [32] and Moayeri and Zendehbudi [29] considered the elastic properties of the undisturbed artery, but the stenosis zone is presumed to be rigid. The simulation with deformable wall predicted an increased time-averaged pressure drop, while the maximum shear stress on the artery wall was found to decrease [29]. Though the effect of artery wall compliance is potentially important, it is of secondary significance in the hemodynamic study through stenosed arteries, as atherosclerosis considerably reduces the elastic property of the wall [2]. On the other hand, Nerem [7] propounded that the primary effect of wall deformability is on the shape of the pulsatile waveform in flow and pressure, so that when the actual physiological flow waveform is considered, the wall can be considered as rigid. Therefore, several investigations ignored the effect of wall motion and considered the tube to be rigid [9, 33, 34] for the flow study through stenosed artery.

A variety of pulsatile flow waveforms are observed in the arteries (see e.g., Zendehbudi and Moayeri [35]; McDonald [36]), differing from individual to individual, with position in the arterial system and in response to physical exertion. However, most experimental and numerical studies assumed a simple harmonic pulse upstream to the stenosis [14, 16, 25, 37]. Johnston et al. [26] and Wiwatanapataphee et al. [38] performed their studies with physiological waveforms observed in the human right coronary artery, while Long et al. [39] considered the volumetric flow waveform in the carotid artery, measured by the Doppler ultrasound. Mills et al. [34] compared the flows with simple pulsatile waveform and physiological waveform in the canine femoral artery. Sherwin and Blackburn [9] studied the instability and transition to turbulence with different flow waveforms. The flow with two-harmonic waveforms having a higher peak to mean velocity ratio was found to be more stable than the flow with single harmonic waveform. Koo et al. [40] uses the CFD technique to measure FFRCT non-invasively from the reconstructed CCTA images with invasive FFR as a standard index. Plaque accumulation is predominant at bends and bifurcations where blood flow is disturbed, where a secondary flow from the main vessel similar to collaterals is created, where flow separation regions are found [41]. Govindaraju et al. [42] explains the role of fluid dynamics in measuring the severity of CAD.

The time period of the flow oscillation, expressed in terms of Womersley number (𝛼=𝑅₀𝜔/𝜈), is influenced by the physiological condition of the body and affects the hemodynamic parameters in the blocked artery. Following the recommendations of Caro et al. [43] and Tutty [44], 𝛼 can vary from 1 to 12.5. Buchanan et al. [25] observed a measurable effect of the Womersley number on the wall hemodynamic parameters where an increase in Womersley number resulted in an increasing disturbance in the fluid element pathline. Sherwin and Blackburn [9] reported that the blood flow remains stable and laminar up to a higher Reynolds number at lower Womersley number.

The decision of surgical intervention on a patient having atherosclerosis is taken based on the degree of stenosis [24]. However, irrespective of the degree, the location of the stenosis in the arterial system and the physiological condition of the body may have different implications on the hemodynamic parameters since the waveform and the local Womersley number of the pulsatile flow differ in different sections of the circulatory system. Wall shear stress and pressure drops are two important hemodynamic parameters based on which the severity of the disease on the patient’s health can be interpreted. In addition, the altered flow features, like jet expansion and contraction due to the presence of the stenosis, have adverse pathological implications. The local pulsatile waveform and time period of oscillation greatly affect the flow pattern and hemodynamic parameters. The main objective of the present work is to numerically study the effects of different waveforms and time periods of oscillation on the flow pattern and hemodynamic parameters. Presented here are the obtained pulsatile flow patterns and wall shear stress distributions in stenosed tubes with three typical flow waveforms found in human common carotid arteries. The laminar flow of blood has been modeled in an axisymmetric stenosed artery with 50% diametral occlusion. The artery wall is considered to be noncompliant, and the blood is assumed to be Newtonian.

2. Physical Model

Figure 1 shows the half-section of the axisymmetric stenosed artery in which the flow has been analysed. The artery is treated as a rigid tube of diameter 𝑑 (radius 𝑅₀), which is divided into three segments, namely, the inlet segment, upstream to the stenosis, having length 𝑙_𝑢=4𝑑, the stenosed segment of length 𝑙_𝑠=𝑑, and the downstream segment of length 𝑙_𝑑=45𝑑. The degree of occlusion (𝑆) is expressed as 𝑆=((𝑑−𝑑_𝑜)/𝑑)×100%, where 𝑑₀ denotes the diameter at the throat of the stenosis, and the study has been performed with 50% occlusion. The shape of the stenosis is expressed by the following equation [45]: 𝑅=𝑅₀11−2𝑆1−cos𝜋𝑧𝑑²,for𝑙_𝑢≤𝑧≤𝑙_𝑢+𝑙_𝑠. (1)

Figure 1 

Geometrical configurations of the computational domain with cosine-shaped stenosis.

We have expressed the pulsatility in the flow through a variation in Reynolds number [14, 45] around a mean Reynolds number (Re_mean) of 200. The Reynolds number (Re=𝜌𝑢𝑑/𝜇) is defined on the basis of the diameter of the unstenosed artery, time-averaged mean velocity of flow, and the density and viscosity of blood, which are taken as 𝜌=1056Kg/m³ and 𝜇=.00345 PaS, respectively. Three different types of pulsatile waveforms (shown in Figure 2(a)) are considered for the present analysis, representing the occurrence of the stenosis at different locations of the arterial system:  Profile I:Re(𝑡)=Re_mean1+⁵_𝑛=1𝐴_𝑛𝑡Cos2𝑛𝜋𝑇−𝐵_𝑛. (2a) (see Table 1).  Profile II: Re(𝑡)=Re_mean𝑡1+𝐴Sin2𝜋𝑇𝑡+𝐵Cos4𝜋𝑇, (2b) where  𝐴=0.75,𝐵=−0.75.   Profile III: Re(𝑡)=Re_mean𝑡1+𝐴Sin2𝜋𝑇, (2c)where  𝐴=1.

Table 1 

(a)

(b)

Figure 2 

Details of pulsatile flow waveform considered for the analysis: (a) three profiles together, and (b) important time planes considered for analysis of Profile I.

The first profile represents the physiological waveform in the human right coronary artery [38], the second one is the physiological waveform adopted by Sherwin and Blackburn, [9], which is probable at a downstream location in the arterial system, and the third velocity profile is of a simple sinusoidal waveform adopted by many researchers for the blood flow studies.

The particular input pulse evaluates the Womersley number (𝛼=𝑅₀𝜔/𝜈), where 𝑅₀ is the radius of the artery,  𝜈 is the kinematic viscosity of the blood, and 𝜔 is the fundamental frequency, based on the time period of oscillation. Following the recommendations of Caro et al. [43] and Tutty [44], we adopted three different values of 𝛼 as 4, 7.5, and 12.5 following different physiological conditions of the body. In a flow distribution network that progresses from a large tube to many small tubes (e.g., a blood vessel network), the Womersley number is also changing progressively from a higher value to the lower value. The Womersley numbers tend to 1 at the level of the terminal arteries. In the arterioles, capillaries, and venules, the Womersley numbers are less than one. In these regions the inertia force becomes less important, and the flow is determined by the balance of viscous stresses and the pressure gradient. Some typical values for the Womersley number in the cardiovascular system under resting condition is given by Fung [46], which shows that in the ascending aorta the value is around 13.2, and it drops to 4.4 in the carotid artery, which can further reduced to 0.005 inside the capillaries.

3. Numerical Model

3.1. Theory

The Navier-Stokes equations along with the continuity equation for incompressible fluid flow have been considered for the solution of axial and radial velocities (𝑢 and 𝑣 resp.) and pressure (𝑝). The governing equations are as follows: continuity:1𝑟𝜕(𝑟𝑣)+𝜕𝑟𝜕𝑢𝜕𝑧=0, (3) axial momentum:𝜕𝑢+1𝜕𝑡𝑟𝜕()+𝜕𝑟𝑟𝑢𝑣𝜕𝑢²1𝜕𝑧=−𝜌𝜕𝑝+1𝜕𝑧𝑟𝜕𝜇𝜕𝑟𝜌𝑟𝜕𝑢+𝜕𝑟𝜕𝑣+𝜕𝜕𝑧𝜕𝑧2𝜇𝜌𝜕𝑢,𝜕𝑧 (4a) Radial Momentum:𝜕𝑣+1𝜕𝑡𝑟𝜕𝜕𝑟𝑟𝑣²+𝜕()1𝜕𝑧𝑢𝑣=−𝜌𝜕𝑝+1𝜕𝑟𝑟𝜕𝜕𝑟2𝜇𝜌𝑟𝜕𝑣−𝜕𝑟2𝜇𝜌𝑣𝑟²+𝜕𝜇𝜕𝑧𝜌𝜕𝑢+𝜕𝑟𝜕𝑣.𝜕𝑧 (4b)

No slip boundary condition has been imposed on the wall, while a linearised convective out flow boundary condition [47] is considered at the outlet. The prescribed inlet velocity profile is generated performing simulations on sufficiently long tubes [48] with oscillating plug flow, given by (2a)–(2c), at the inlet. The length of the tube is chosen so that the flow attains its spatially fully developed profile (similar to the Womersley velocity profile for Profile III, given by [49], corresponding to inlet flow as given in (2c) at the exit.

The conservation equations of momentum and mass are solved simultaneously with their appropriate boundary conditions by an explicit finite difference-computing scheme as adopted by Banerjee et al. [50] and Nag and Datta [51]. The variables are defined following a staggered grid arrangement. The axial and radial velocities are specified at the cell faces while the pressure and the fluid properties are calculated at the cell centers. The advection and diffusion terms of the momentum equations are discretized following the power law scheme [52]. At each time step, the velocity variables are advanced explicitly using the axial and radial momentum equations. Thereafter, the pressure and velocity corrections are enforced through the continuity equation to obtain the corresponding temporal values. All the hemodynamic parameters, like wall shear stress, and pressure drop, are calculated based on the temporal values before advancing the variables in time. The time increment for the explicit advancement is calculated satisfying both the stability conditions stipulated by the grid Fourier number and the Courant-Freidrichs-Lewy (CFL) criteria. The solution is advanced in time till one complete time cycle is completed.

3.2. Hemodynamic Wall Parameters

One objective in characterizing flows in the stenosed artery is to quantify the hemodynamic forces that the artery wall experiences due to the system’s input pulse under the changed geometrical scenario. In order to clearly elucidate the effects on the flow field, certain hemodynamic wall parameters have been calculated, which are indicator functions for the onset and progression of arterial diseases, for example, atherosclerosis and thrombosis [25].

We have used the time-averaged wall shear stress, to describe an important hemodynamic wall parameter as follows: 1WSS=𝑇1𝜌𝑢^2𝑇₀𝜏_wall𝑑𝑡. (5)

The oscillatory shear index (OSI) is used to describe the unsteady nature of blood flow through the artery and is defined as 1OSI=2∫1−^𝑇0𝜏_wall𝑑𝑡∫^𝑇₀𝜏_wall=1𝑑𝑡21−WSS_meanWSS_mag. (6)

This index quantifies the pulsatility and the main direction of the flow. It ranges from 0 (steady, forward flow throughout the cardiac cycle) to 1 (fully reversed steady flow). An OSI of 0.5 corresponds to a pure oscillating flow with a WSS_mean of 0. It is also effective at locating points of time-averaged separation and reattachment points [25].

4. Result and Discussions

4.1. Grid Sensitivity

A grid sensitivity analysis has been performed to optimize the number of grids. Five different mesh sizes (Table 2), represented by 𝑀₁, 𝑀₂, 𝑀₃, 𝑀₄, and 𝑀₅, arranged according to their increasing degree of fineness, have been considered. The pulsatile blood flow is simulated considering simple harmonic waveform (Profile III) and with Re_mean=200 and 𝛼=12.5. Four different degrees of stenosis, varying from 𝑆 = 30% to 65%, are considered for the simulation, and the hemodynamic parameters at each time step have been evaluated. Table 2 compares the peak wall shear stress at 𝑡/𝑇=0.5 for each of the cases. It has been found that refining the grid from 𝑀₄ to 𝑀₅ does not change the results (up to second digit after decimal) in any of the cases. Therefore, we adopted grid 𝑀₄ to simulate the results in our computation.

Table 2 

Variation of peak wall shear stress at different degree of stenosis for different mesh size: test for grid-independent result.

4.2. Model Validation

The present model has been validated by comparing the numerical predictions against the experimental results of Ojha et al. [14]. Flow geometry similar to the experimental one has been considered for the validation with a constriction of trapezoidal profile and a 45% occlusion based on area. The fluid properties and flow conditions are also considered to be the same as in the experiments (𝜌=755 kg/m³, 𝜇=0.00143 N/(m² s) and 𝛼=7.75). In absence of a definite expression of the input waveform considered in the experiments, we have considered an inlet pulse for the analysis as [53] (𝑡)𝑡Re=575+360Sin2𝜋𝑇+2.23. (7)

The above waveform results in a similar velocity variation with time as considered in the experiment. This has been represented in Figure 3(a), which shows that the assumed inlet centerline velocities are similar to those obtained experimentally throughout one cardiac cycle. The temporal variations of the centerline velocity at three axial locations following the stenosis (at distances of 𝑑, 2.5𝑑, and 4.3𝑑 from the stenosis throat) are used to validate the numerical predictions against the experimental data (Figure 3(b)). A generally good agreement (maximum relative error = 2.6%) has been achieved with small differences that can be attributed to the difference in the input wave pulse from the experimental values.

(a)

(b)

Figure 3 

Comparison of the present numerical computation with that carried out by Ojha et al. [14]: (a) inlet centerline velocity profile and (b) centerline axial velocity at different poststenotic region. Here, 𝑧 denotes the distance from the stenosis throat.

Finally a comparison between the measured [14] and computed velocity profiles at different time stamps (𝑇₁=0.175,𝑇₂=0.525,𝑇₃=0.875) and at axial distances of 𝑑, 2.5𝑑, and 4.3𝑑 from the throat of the stenosis is shown in Figure 4. The figure shows an excellent agreement between the predicted and measured values at all times and positions. The above validations ensure the predictive capability of the code.

(a)

(b)

(c)

Figure 4 

Model validation: axial velocity profile at three different locations measured from the center of stenosis for three different time planes: (a) at 𝑇1, (b) at 𝑇2, and (c) at 𝑇3. Here 𝑧 represents distance from the stenosis throat.

4.3. Study of Predicted Flow Pattern and Hemodynamic Parameters

The validated numerical code is used to simulate the flow with Re_mean=200 through a stenosed artery with 50% diametral occlusion. Three different input pulsatile waveforms, as given by (2a)–(2c), and three Womersley numbers (𝛼=4.0, 7.5, and 12.5) have been considered in the study. We have first presented the flow pattern and hemodynamic parameters for one case (considered as the base case) with input waveform given by Profile I and 𝛼=4.0. This waveform profile is compatible to that in the human coronary artery (about 4 mm in diameter), and the Womersley number pertains to the resting heartbeat of 74 minute⁻¹.

The temporal flow pattern in the stenosed artery, particularly in the poststenotic zone, is of interest due to its pathological implications. Considerable variations in the instantaneous streamline pattern are observed throughout the cycle (Figure 5).

Figure 5 

Instantaneous stream lines for the base case.

Figures 6(a) and 6(b) show the axial variation of the instantaneous wall shear stress (WSS) at different time of the cardiac cycle, around the stenosis and in the post stenotic zone immediately following the stenosis, respectively. The peak WSS at 𝑡/𝑇=0.20 is almost 10.8 times higher than that at the start of the cardiac cycle and this value reduced to 6.06 for 𝑡/𝑇=0.67. From Figure 6(b), it is quite clear that the recirculation length increases as the flow accelerates during the systole, and the variation in the recirculation length experiences a time lag compared to the variation in the flow rate.

(a)

(b)

Figure 6 

Time-dependent behavior of WSS for base case: (a) at the stenosis zone and (b) at the recirculation zone.

Figure 7 presents the axial variations of the time-averaged wall shear stress (WSS) and the OSI for the base case. The time-averaged separation (𝑧^∗=4.83) and reattachment (𝑧^∗=12.81) points are clearly delineated in the WSS profile (i.e., where WSS=0) which also correspond to local maxima values of the OSI. The OSI value drops down to nearly zero at 𝑧^∗=4.47, that is, just before the throat, indicating that the effect of flow pulsatility is minimum there. The poststenotic region is characterized by a high OSI value, indicating that the flow consists of a high-inertia core jet that does not exhibit flow reversal, surrounded by a recirculation zone that exhibits greater pulsatility with lower net flow.

Figure 7 

Time-averaged WSS and OSI distribution for the base case. (Inset shows the drop and rise of OSI near the throat.)

The time-average nondimensional pressure drop (Δ𝑝^∗=WP−WP_inlet) at various locations is plotted for the base case (Figure 8). The results show that the unsteady flow produces a 1.48 times higher pressure drop than the steady flow. Beyond the location of minimum pressure, a recovery of pressure takes place due to the Bernoullian conversion of the kinetic energy following the expansion.

Figure 8 

Time average dimensionless pressure drop for the base case.

4.4. The Effects of Different Pulsatile Flow Waveforms

The waveform of the pulsatile flow in the coronary artery has five harmonics in it. However, the waveform changes as the blood flows away from the heart in the arterial system. In fact, in the small arterioles far away from the heart the flow becomes almost steady [13]. It is therefore of relevance to study how much the variation in waveform affects the hemodynamic parameters.

We have examined the details of the flow fields for three different pulsatile waveforms (as represented by (2a), (2b), and (2c)) for 𝛼=4, without varying the mean Reynolds number. The comparisons have been presented at four different times of the cardiac cycle, namely, at 𝑡/𝑇=0.25, 0.50, 0.75 and 1.0, through the instantaneous streamlines at these times for the three different waveforms (refer to Figures 9(a)–9(c)).

(a)

(b)

(c)

Figure 9 

Instantaneous streamline plots for different pulsatile flow waveforms at four different time planes: (a) Profile I, (b) Profile II, and (c) Profile III.

The flow field with Profile I is already described Section 4.3 (albeit the four time instants in Figure 9(a) differ from those in Figure 5). The flow in the proximal side of the stenosis remains parallel till it comes close to the stenosis at all the four instants. The length of the recirculating vortex distal to the stenosis, changes depending upon the jet inertia following the constriction as discussed earlier.

Figure 9(b) represents the flow fields at the four instants with Profile II. For this case, the peak flow is achieved at 𝑡/𝑇=0.25. The flow approaches the stenosis as parallel flow and takes a turn close to the proximal wall of the stenosis to confront the constriction. Distal to the stenosis a long recirculation vortex is established. Two counter-rotating vortices are seen at the core of the vortex. A major change in the flow field is observed at 𝑡/𝑇=0.5. A recirculation adjacent to the wall is now seen proximal to the occlusion, while the recirculation after the occlusion grows in size. Figure 2 clearly shows that at 𝑡/𝑇=0.5 the flow rate is close to the minimum flow rate of the cycle in case of Profile II. The inertia of the flow is so low at the instant that the flow cannot confront the occlusion which close to the wall. The recirculation length following the stenosis grows as the jet created by the constriction contracts during the deceleration period beyond 𝑡/𝑇=0.25. At 𝑡/𝑇=0.5 the jetting action has almost subsided and the distal recirculation reaches close to its maximum. As the jet expands again following the increase in flow rate in the diastole, a new recirculating vortex is created distal to the stenosis. However, the recirculation length at 𝑡/𝑇=0.75 (when the flow rate reaches its peak diastolic value) is shorter than that at 𝑡/𝑇=0.25. This is because the peak flow rate during diastole is much less than the peak flow rate in systole. The recirculation grows again during the deceleration of the flow after 𝑡/𝑇=0.75. The distal recirculating vortex at 𝑡/𝑇=1.0 is reduced further than that at 𝑡/𝑇=0.75. However, though the flow rate at 𝑡/𝑇=1.0 is close to the minimum again, no recirculation is observed proximal to the stenosis at this time.

Finally, Figure 9(c) presents the flow field during the four time instants with the sinusoidal profile representing simple harmonic flow. The flow fields with this profile are clearly different from the flow fields of the corresponding time instants with physiological profiles. The differences are in the lengths of the distal recirculation zones and formation of the proximal recirculation vortex. The growth in the distal recirculation during the decelerating phase is observed. The proximal recirculation is observed at 𝑡/𝑇=0.75, when the flow rate is zero. Due to that there is no net flow, the recirculating vortices on either side of the stenosis extend almost to the centerline.

The comparison of instantaneous WSS distributions for the three waveforms of flow at time 𝑡/𝑇=0.50 is presented in Figure 10(a). The maximum WSS is observed at the throat of the stenosis for each of the profiles. The peak WSS reached with the simple harmonic profile (Profile III) is much higher than the peak WSS with either of the physiological profiles at this time instant. The peak WSS of Profile I is 1.75 times less than the peak value with Profile III, while the peak WSS of Profile II is even smaller. The wall shear stress depends on the axial velocity gradient on the wall. It is evident from Figure 2(a) that at time 𝑡/𝑇=0.5 the instantaneous flow rate is the maximum with Profile III, while the second highest flow rate is achieved with Profile I and the flow rate of Profile II is the minimum of the three profiles. The higher flow rate increases the maximum velocity and the wall velocity gradient in the artery. Accordingly, the shear stress at the throat of the artery reaches a higher value with the higher instantaneous flow rate. Figure 10(a) further reveals that in each of the three cases the wall shear stress reaches from a low value to the maximum at the stenosis throat. After reaching the maximum value at the throat the wall shear stress drops even faster in the distal side of the stenosis. As a result a high wall shear stress gradient is established on the proximal and distal sides of the stenosis.

(a)

(b)

Figure 10 

Axial variation of different hemodynamic parameters for different pulsatile flow waveforms: (a) WSS at 𝑡/𝑇=0.50 and (b) OSI at 𝛼=4.

However, Figure 10(a) compares the wall shear stress distribution with the three flow profiles at a particular time instant of the complete cardiac cycle. The relative magnitudes of the wall shear stresses at the other time steps will be different. Table 3 lists the maximum and the minimum wall shear stress at the four different time instants of 𝑡/𝑇=0.25, 0.5, 0.75 and 1.0 for the three different wave profiles considered. The maximum value of the shear stress at any point for all the three profiles is achieved with Profile II at 𝑡/𝑇=0.25. It is because the maximum instantaneous Reynolds number, indicating the maximum flow rate, is achieved in that case. Therefore, the wall shear stress distribution and the distribution of wall shear stress gradient are highly dependent on the actual flow waveform through the artery. It is therefore not prudent to study the flow with an arbitrarily assumed profile as that may lead to wrong results.

Table 3 

Maximum and minimum wall shear stress at different time instants for different pulsatile flow waveform when 𝛼=4.

In order to get the time-averaged results, we study the effect of OSI with the different pulsatile flow waveform considered for the analysis (see Figure 10(b)). The figure clearly indicates that the location of time-averaged reattachment point is greatly influenced by the nature of the input flow waveform. For the Profile I the time-averaged point of reattachment is located at 𝑧^∗=12.89𝑑 measured from the inlet; on the other hand, it is almost 16.06% and 10.4% more for Profile II and III,respectively. This clearly suggests the significance of choosing the most physiological waveform, (i.e., Profile I) over the other two waveforms. Also, the far-downstream OSI values of the Profile III show considerable deviation from that due to Profile I. Therefore, choice of the correct waveform is important not only for predicting the correct hemodynamic parameters for the stenosed regions, but also the far-downstream regions. From this perspective, Profile II yields a better far-downstream OSI values (i.e., closer to that due to Profile I).

Figure 11 shows the time-average axial distribution of Δ𝑝^∗ for the three different profiles. It is observed that the time-averaged irreversible pressure drop (Δ𝑝_𝐼), as a result of the occlusion in the artery, is the minimum when the physiological waveform of the coronary artery (Profile I) is considered for the simulation. On the other hand, the maximum time average irreversible pressure drop occurs with the sinusoidal profile (Profile III).

Figure 11 

Time average dimensionless pressure drop for different pulsatile flow waveforms at = 4.

4.5. The Effects of the Womersley Number

We have already shown that results of the physiological waveform differ considerably from the sinusoidal waveform both in terms of characteristics and values. It is therefore worthy to study the effect of Womersley number with the physiological wave profile (Profile I) for a better pathological interpretation. The Womersley numbers of 4.0 (whose results are already presented as the base case in Figures 5–8), 7.5, and 12.5 are considered for the study.

The instantaneous streamlines for 𝛼=7.5 and 12.5 at the selected time levels are shown in Figures 12 and 13 respectively. Referring to the streamlines in Figures 5, 12, and 13, the difference in flow patterns can be noticed for different Womersley numbers. For a comparatively low Womersley number (𝛼=4.0), the flow pattern looks similar to those of the quasiunsteady flow, and the flow field is largely influenced by the instantaneous flow rate. There is always a recirculating vortex formed downstream to the constriction for nonzero instantaneous flow rate at low Womersley number.

Figure 12 

Instantaneous streamline plots for Profile I at 𝛼=7.5.

Figure 13 

Instantaneous streamline plots for Profile I at 𝛼=12.5.

In contrast, for higher Womersley numbers, the flow field shows a stronger dependence on the flow rate at the previous time levels. The shape of the recirculating vortex shows larger deviations from their quasi-unsteady counterparts (see Figure 12). For example, with 𝛼=7.5 and at 𝑡/𝑇=0.13, the flow remains attached to the distal wall of the stenosis, and no recirculation zone is predicted (which is not the case for 𝛼=4.0). At 𝑡/𝑇=0.20, the flow separates from the stenosis wall, and the vortex continues to grow in size. Finally, at 𝑡/𝑇=0.53, the large vortex breaks into two pieces with a secondary vortex following the primary one. As the time increases, the secondary vortex continues to grow in size and the primary vortex gets completely washed away (at 𝑡/𝑇=1.0). Clearly, the flow field appears to evolve at a more sluggish pace as compared to that in Figure 5. For example, the vortex breakdown for 𝛼=7.5 appears to occur at 𝑡/𝑇=0.50, while the same for 𝛼=4.0 occurs at an earlier phase (𝑡/𝑇=0.20).

Consistent to the above trend, the flow field for 𝛼=12.5 shows an even larger transient inertial effects relative to the viscous effects. In Figure 13, the vortex detachment is seen to occur at an even later phase (𝑡/𝑇=0.53) than that observed in Figure 12. A remarkable difference in flow pattern is also observed at 𝑡/𝑇=0.13 compared to the case when 𝛼=7.5. A strong recirculation zone, which is detached from the stenosis wall, is formed. There occurs a strong inward velocity of the fluid at the occlusion, which essentially tries to attach the vortex to the distal wall of the stenosis. At the peak of the systolic pulse (time level 𝑡/𝑇=0.20), the recirculation zone attaches the distal wall of the stenosis but remains short in length. With the increase in time (𝑡/𝑇=0.27), the core of the vortex is shifted towards the downstream direction. During the diastolic phase, at 𝑡/𝑇=0.53, the vortex grows to a very large size, and a fissure appears in it. One of the vortices (primary) travels downstream with the flow, while the other (secondary) remains attached to the stenosis distal wall. When the flow reaches the peak diastolic time, that is, at 𝑡/𝑇=0.67, the secondary vortex grows in size, while the strength of the primary vortex decreases continuously. As the time further increases to 𝑡/𝑇=1.0, the primary vortex is no longer visible, and thereafter the flow repeats itself.

The axial WSS distribution at a specified time level (𝑡/𝑇=0.20) for 𝛼=4.0,7.5 and 12.5 follows almost the same distribution up to 𝑧=4.65𝑑 as shown in Figure 14. Beyond that, a variation in the distribution of WSS is observed. This is due to the varying flow field in the poststenotic region under different Womersley number. At 𝛼=4.0, the wall shear stress sharply decreases from the peak value at the throat. The wall shear stress then remains negative in the wall region shown in the figure. However, for the other two cases (𝛼=7.5 and 12.5), an oscillation in the shear stress is observed in the poststenotic zone. The oscillation even results in the multiple changes in sign of the wall shear stress. Table 4 presents the numerical values of maximum and minimum WSS at two specified time levels (peak systole and peak diastole) for different values of 𝛼. It is clear from the table that there is only a very little change in the maximum value of WSS with 𝛼. However, there is a sharp change in the values of the minimum WSS as reported in the table.

Table 4 

Maximum and minimum wall shear stress at peak systolic and peak diastolic instants for different values of 𝛼 and for Profile I.

Figure 14 

Axial variation of WSS for the base case for different Womersley numbers at 𝑡/𝑇=0.20.

Figures 15 and 16 present the hemodynamic wall parameters for each Womersley number. The main features of the flow field, that is, the separation and reattachment of the flow distal to the stenosis, are distinct despite the time averaging. Figure 15 compares the axial distribution of the OSI values for 𝛼=4.0, 7.5 and 12.5. The OSI values are larger (in most of the sections of the artery studied here) for the higher values of Womersley number. At larger 𝛼, the high-inertia core of the flow gets little time for a complete flow reversal, which in turn renders the flow in the low-inertia peripheral regions more fluctuating in nature. All the profiles show a near-zero value right at the upstream of the throat, where there is hardly any flow reversal at any radial location. The time-averaged separation and reattachment points are clearly defined in the distribution of OSI (corresponding to the maximum values of OSI, see Figure 15), which is again justified by the time-averaged WSS (see Figure 16(b)). From Figure 16(b) it is quite clear that the reattachment point locations from the stenosis are not monotonic, which is similar to the observations made by Buchanan et al. [25]. The time-averaged point of reattachment is observed to be maximum (𝑧^∗=13.19) for the intermediate Womersley number of 7.5, while it is least (𝑧^∗=12.12) for the maximum value of 12.5.

Figure 15 

Spatial distribution of OSI for three different Womersley numbers for Profile I.

(a)

(b)

Figure 16 

Spatial distribution of time-averaged WSS for three different Womersley numbers for Profile I: (a) full length view and (b) enlarged view, showing only the recirculation zone.

From the Figure 16(a), the nature of the time-averaged WSS is almost similar for all the Womersley number considered for the analysis. But a decrease in peak time-averaged WSS is observed with the increase in Womersley number. For 𝛼=4.0 the peak value is 1.36, while it is almost 2% less for 𝛼=12.5.

5. Biomedical Implication

Magnetic resonance (MR)- [54] and intravascular ultrasound (IVUS)- [55] based techniques have been described to provide valuable diagnostic and anatomical information; however, their temporal and spatial resolution is usually not sufficient for deriving information based on the differentiated velocity vector field such as the wall shear stress distribution [56]. CFD techniques, on the other hand, can provide a significantly larger amount of information on coronary flow dynamics, with greater spatial and temporal resolution, noninvasively, on a patient-specific basis and within a completely virtual environment.

The ability to perform such flow simulations on a patient-specific basis in a totally virtual environment, noninvasively and thus repetitively, could lend these techniques an important role as a predictive tool, which means that such simulations could predict potential flare-up scenarios where preventive measures could be instituted. Several investigations that are underway, reveal the important contributions that CFD can offer when combined with presently available diagnostic techniques [57–59].

6. Conclusion

Pulsatile blood flow in stenosed artery has been studied to analyze the influence of different flow waveform and Womersley number on the key hemodynamic parameters. The following inferences can be drawn from the investigation.(1)The instantaneous streamlines of the physiological flow show a considerable variation throughout the cycle, characterized by the formation of a vortex just downstream of the stenosis zone during the accelerating phase and the vortex elongation and formation of a second corotating vortex during the decelerating phase. (2)The magnitude of WSS is the largest at the peak systolic time, that is, when the flow is the maximum in the cycle.(3)The recirculation length increases as the flow accelerates during the systole. The variation in the recirculation length experiences a time lag compared to the variation in the flow rate. The lag depends on the strength of the transitional inertia effect given by the Womersley number. (4)The unsteady flow produces a higher-pressure drop than the steady flow.(5)Presence of the stenosis leads to a net increase in the OSI at the far-downstream side of an artery.(6)The instantaneous flow field and the hemodynamic wall parameters are highly dependent on the actual flow waveform through the artery. The study shows how an arbitrarily assumed profile may lead to different connotation.(7)At higher Womersley number, the instantaneous flow field shows a stronger influence of the flow profiles of previous time levels. With an increase in Womersley number the peak value of time-averaged WSS decreases, while the OSI values exhibit an increasing trend.