Vortex dynamics under pulsatile flow in axisymmetric constricted tubes

Improved understanding of how vortices develop and propagate under pulsatile flow can shed important light on the mixing and transport processes occurring in such systems, including the transition to turbulent regime. For example, the characterization of pulsatile flows in obstructed artery models serves to encourage research into flow-induced phenomena associated with changes in morphology, blood viscosity, wall elasticity and flow rate. In this work, an axisymmetric rigid model was used to study the behaviour of the flow pattern with varying degrees of constriction (d0) and mean Reynolds (R̄e) and Womersley numbers (α). Velocity fields were obtained experimentally using Digital Particle Image Velocimetry, and generated numerically. For the acquisition of data, R̄e was varied from 385 to 2044, d0 was 1.0 cm and 1.6 cm, and α was varied from 17 to 33 in the experiments and from 24 to 50 in the numerical simulations. Results for the Reynolds numbers considered showed that the flow pattern consisted of two main structures: a central jet around the tube axis and a recirculation zone adjacent to the inner wall of the tube, where vortices shed. Using the vorticity fields, the trajectory of vortices was tracked and their displacement over their lifetime calculated. The analysis led to a scaling law equation for maximum vortex displacement as a function of a dimensionless variable dependent on the system parameters Re and α.


I. Introduction
Research on the dynamics of pulsatile flows through constricted regions has multiple applications in biomedical engineering and medicine. Cardiovascular diseases are the primary cause of death worldwide (accounting for 31% of total deaths in 2012) [1]. More than half of these deaths could have been avoided by prevention and early diagnosis. Atherosclerosis is characterized by the accumulation of fat, cholesterol and other substances in the intima layer, creating plaques which obstruct the arterial lumen. Atherosclerotic plaque fissuring and/or breaking are the major causes of cardiovas-cular stroke and myocardial infarct. Several factors leading to plaque complications have been reported: sudden increase in luminal pressure [2], turbulent fluctuations [3], hemodynamic shear stress [4], vasa vasorum rupture [5], material fatigue [6] and stress concentration within a plaque [7]. As most of these factors are flow dependent, understanding how a pulsatile flow behaves through a narrowed region can provide important insights for the development of reliable diagnostic tools.
Flow alterations due to arterial obstruction (i.e. stenosis) and aneurysms have been reported in the literature [8][9][10][11]. Simplified models of stenosed arteries have used constricted rigid tubes [12][13][14][15]. At moderate Reynolds numbers, the altered flow splits at the downstream edge of the constriction into a high-velocity jet along the centreline and a vortex shedding zone separating from the inner surface of the tube wall (recirculation flow). Increasing Reynolds numbers lead to a transition pattern and ultimately, turbulence.
Several experimental studies have been reported. Stettler and Hussain [16] studied the transition to turbulence of a pulsatile flow in a rigid tube using one-point anemometry. However, one-point anemometry fails to account for different flow structures which may play an important role at the onset of turbulence. Peacock [17] studied the transition to turbulence in a rigid tube using two-dimensional particle image velocimetry (PIV). Chua et al [18] implemented a three-dimensional PIV technique to obtain the volumetric velocity field of a steady flow. The work of Ahmed and Giddens [19,20] studied the transition to turbulence in a constricted rigid tube with varying degrees of constriction using twocomponent laser Doppler anemometry. For constriction degrees up to 50% of the lumen, the authors found no disturbances at Reynolds numbers below 1000. This study was not focused on describing vortex dynamics; however, the authors reported that turbulence, when observed, was preceded by vortex shedding.
Several authors carried out numerical studies. Long et al. [15] compared the flow patterns produced by an axisymmetric and a non-axisymmetric geometry, finding that flow instabilities through an axisymmetric geometry was more sensitive to changes in the degree of constriction. The work of Isler et al. [21] in a constricted channel found the instabilities that break the symmetry of the flow. Mittal et al. [22] and Sherwin and Blackburn [23] studied the transition to turbulence over a wide range of Reynolds numbers. Mittal et al. used a planar channel with a one-sided semicircular constriction. They found that downstream of the constriction the flow was composed of two shear layers, one originating at the downstream edge of the constriction and the other separating from the opposite wall. For Reynolds numbers above 1000, the authors reported transition to turbulence due to vortex shedding. Moreover, they found through spectral analysis that the characteristic shear layer frequency is associated with the frequency of vortex formation. Sherwin and Blackburn [23] studied the transition to turbulence using a three-dimensional axisymmetric geometry with a sinusoidal constriction. Based on the results of linear stability analysis, the authors reported the occurrence of Kelvin-Helmholtz instability. This reaffirms that instabil-ities involving vortex shedding take place in the transition to turbulence.
Finally, several studies compared experimental and numerical results. Ling et al. [24] compared numerical results with those obtained by hot-wire measurements. As mentioned above, onedimensional hot-wire measurements cannot be used to identify flow structures. Griffith et al. [25], using a rigid tube with a slightly narrowed section resembling stenosis, compared numerical results with experimental data. Using stability analysis by means of Floquet exponents, they demonstrated that the experimental flow was less stable than that of the simulated model. For low Reynolds numbers (50 − 700), the authors found that a ring of vortices formed immediately downstream of the stenosis and that its propagation velocity changed with the degree of stenosis. The work of Usmani and Muralidhar [26] compareed flow patterns in rigid and compliant asymmetric constricted tubes for a Reynolds range between 300 and 800 and Womersley between 6 and 8. The authors reported that the downstream flow was characterized by a high velocity jet and a vortex whose evolution was described qualitatively. The aforementioned studies highlight the significance of studying vortex dynamics, since vortex development precedes turbulence, and ultimately, contributes to the risk of a cardiovascular stroke.
The aim of this work is to characterize vortex dynamics under pulsatile flow in an axisymmetric constricted rigid tube. The study was carried out both experimentally and numerically for different constriction degrees and with mean Reynolds numbers varying from 385 to 2044 and Womersley numbers from 17 to 50. The results show that the flow pattern in these systems consists primarily of a central jet and a vortex shedding layer adjacent to the wall (recirculation flow), consistent with the literature. By tracking the vortex trajectory it was possible to determine the displacement of vortices over their lifetime. Vortex kinematics was described as a function of the system parameters in the form of a dimensionless scaling law for maximum vortex displacement.

i. Experimental setup
The experimental setup ( Fig. 1 (a)) consisted of a circulating loop and a Digital Particle Image Velocimetry (DPIV) module capable of obtaining videos and processing velocity fields. The circulating loop consisted of a programmable pulsatile pump (PP), a reservoir (R), a flow development section (FDS) and a tube containing an annular constriction (CT).
The tube with the annular constriction (CT) consisted of a transparent acrylic rigid tube of length L = 51.0 ± 0.1 cm and inner diameter D = 2.6 ± 0.1 cm. The annular constriction consisted of a hollow cylinder 5 ± 1 mm in axial length, which was fitted inside the rigid tube. Keeping the outer diameter of the annular constriction fixed at 2.6 cm, i.e., equal to the inner diameter of the rigid tube, its inner diameter was changed to achieve different degrees of constriction. Tests were carried out using annular constrictions with inner diameters d 0 = 1.6 ± 0.1 cm and d 0 = 1.0 ± 0.1 cm, equivalent to degrees of constriction relative to D of 39% and 61%, respectively. A cross-sectional scheme of the constriction geometry is shown in Fig. 1 (b), along with the coordinate system used throughout the study. The radial direction is represented by r and the axial direction by z, with r = 0 coinciding with the tube axis and z=0 with the downstream edge of the constriction. Within this axis representation, the downstream region is defined by z >0 values and the upstream region is defined by z < 0 values. Finally, this constricted tube (CT) was placed inside a chamber filled with water, so that the refraction index of the fluid inside the tube matched that of the outside fluid.
The tube inlet was connected to the pulsatile pump via a flow development section (FDS), and its outlet was connected to the reservoir. The flow development section was designed to ensure a fully developed flow at the tube inlet and consisted of two sections: first, a conical tube of 35 cm in length whose inner diameter increased from 1 cm to 2.6 cm to provide a smooth transition from the outlet of the pump. Secondly, an acrylic rigid tube of inner diameter D and length of 48D connected to the constricted tube (CT). The reservoir (R) was used to set the minimum pressure inside the tube, which was set as equal to atmospheric pressure in the experiments.
The system was filled with degassed water and seeded with neutrally buoyant polyamide particles (0.13 g/l concentration, 50 µm diameter, DAN-TEC). The DPIV technique was used to obtain the velocity field [27,28]. A 1 W Nd:YAG laser was used to illuminate a 2 mm-thick section of the tube. Images were obtained at a frame rate of 180 Hz over a period equivalent to 16 cycles, using a CMOS camera (Pixelink, PL-B776F). The velocity field was finally computed using OpenPIV open-source software with 32 × 32 pixel 2 windows and an overlap of 8 pixels in both directions.
The region of interest was defined as -0.5< r/D <0.5, 0< z/D <1.5. Within this region no turbulence was observed, as confirmed by spectral analysis of the velocity fields. This observation is also consistent with previous observations [17,24,29]. Due to the pulsatility and the absence of turbulence, it was possible to consider each cycle as an independent experiment, using the ensemble average over all 16 cycles to obtain the final velocity fields.
For each degree of constriction different experiments were carried out, varying the flow velocity at the tube inlet. The pulsation period and the shape of the velocity as a function of t, in z = 0 and r = 0, shown in Fig The experiments were carried out for three different pulsatile periods, T : 0.96s, 2.39s and 3.58s, corresponding to α values of 33.26, 21.10 and 17.22 respectively. The peak Reynolds number upstream of the constriction is defined as Re u = Dv u /ν where v u is the peak velocity measured at the centreline. Four different values of Re u were used as inlet condition upstream of the constriction (see table 1). Similarly, the mean Reynolds at the downstream edge of the constriction is defined asRe = d 0v /ν, wherev is the mean velocity over a whole period, measured at (z, r) = (0, 0  summarized in table 1. The pulsatile pump (PP in Fig. 1 (a)) generated different flow conditions at the tube inlet, [30]. The pumping module consisted of a step-by-step motor driving a piston inside a cylindrical chamber. The control module contained the power supply for the entire system and an electronic board, programmable from a remote PC by means of custom software. The piston motion was programmed using the pump settings for pulsatile period, ejected volume, pressure and cycle shape. The piston motion of the programmable pump is shown in Fig. 2 in normalized units of displacement and time. Figure 2 shows that the piston motion, and hence the pulsatile waveform, are comparable among all Re u values, meaning that the inlet conditions are equivalent for all the experiments.
Finally, the condition of fully-developed flow at the upstream region was verified. The entrance length is usually estimated as L ≈ 0.05ReD. According to the work of Ku [31], the entrance length of a pulsatile flow with α >10 corresponds to the upstream mean Reynolds over a whole cycle,Re u , giving L ≈ 0.05Re u D. In our setup, the flow development section measures 48D,Re u <1000 and α >10 for all experiments, which makes it reasonable to assume a fully-developed flow. Nevertheless, velocity profiles were studied for Re u values of 820, 1187 and 1625 in order to ensure developed flow condition. Profiles are taken from a 1.5D length window upstream of the constriction, at fixed locations z 1 = −1.25D, z 2 = −0.75D and z 3 = −0.25D. Velocity values are normalized by the maximum velocity attained in a pulsatile cycle. Figure 3 shows such profiles at t = 0.25T , t = 0.5T , t = 0.75T and t = T . Profiles show that for each Re u , independently of the time, profiles are equal for z 1 =-1.25D, z 2 =-0.75D and z 3 =-0.25D, which demonstrates that the flow is developed.

ii. Numerical simulation
Three-dimensional numerical simulation using COMSOL software enabled comparison between the experimental data and the study of parameter configurations not achieved experimentally. This entailed solving time-dependent Navier-Stokes equations under the incompressibility condition.
The inlet condition was defined as normal inflow and taken from Re u =820 and Re u =1187 experimental data in the upstream region, where the flow is developed. In order to do this, a fourth-order Fourier decomposition of the experimental flow rate was carried out. The Fourier expansion yields v z (t) = a0 2 + 4 n=1 (a n cos (2πnf t) + b n sin (2πnf t)) with coefficients (×10 −5 , in cm/s) a 0 =196.32, a 1 =61.94, a 2 =8.22, a 3 =13.79, a 4 =3.87, b 1 =-27.66, b 2 =-3.81, b 3 =-5.13, b 4 =-2.84 for Re u =1187. Higher-order terms were not considered, as they do not contribute significantly to the shape of velocity profiles. In order to improve numerical stability, the velocity profile calculated by decomposition was multiplied by a ramp function of time increasing steadily from 0 to 1 in 0.6 seconds. Figure 4(a) shows the axial velocity at the centreline at the downstream edge of the constriction, and presents good agreement with its experimental analogue shown in Fig. 1 (c). Outlet conditions were set to p = 0, p being the difference between the outlet pressure and atmospheric pressure, disabling the normal flow and backflow suppression settings. No-slip conditions were imposed on the inner surface of the tube wall or the constriction. A scheme of the geometry of the constricted tube used in the numerical analysis is shown in Fig. 4(b). The simulations were initialized with the fluid at rest and were run for 20 periods. The first 4 periods were discarded in order to exclude transient effects. The step size for the data output was selected to be the same as for the experiments.
The mesh elements were tetrahedrons except those close to the no-slip boundaries, which were prismatic elements. This was done automatically by COMSOL, based on the boundary requirements. A COMSOL predefined physics-controlled mesh was built, enabling settings that allowed the element size to adapt to the physics at specific regions. A cross-sectional view of the tube and the mesh elements is shown in Fig. 4 (c). Figure 6 compares numerical and experimental velocity profiles in order to validate numerical results. Profiles are located upstream (z = −0.3D) and downstream (z = 0.75D) for the two inlet conditions Re u =820 and Re u =1187. The maximum difference between numerical and experimental profiles was less than 9%.
A mesh convergence analysis was carried out for three different meshes consisting of 180958, 351569 and 951899 elements at the highest Reynolds (Re u =1187) attained in simulation. In Fig.5 the axial velocity at r = 0, z = 0 vs. time (Fig.5(a)) and the velocity profiles at z = 0.5D (Fig.5(b)) are compared for the three different meshes, showing the results to be independent of the mesh chosen. The relative deviation at the centreline, r = 0, is less than 4.5%, although a slight deviation at t = T near the wall is observed. However, this does not affect vortex formation and propagation, which are equally predicted by the three meshes. Finally, a mesh of 351569 elements was chosen for this work in a trade-off between spatial resolution and computational time.

III. Results and discussion
i. Flow pattern The velocity at the centreline, r = 0, at the downstream edge of the constriction, z = 0, as a function of time (see Fig. 1(c)) was used to establish the time reference for all the experiments. The decelerating phase (diastolic phase) was defined as 0 < t < 0.5T and the accelerating phase (systolic phase) as 0.5T < t < T . Other features are dependent on the Reynolds number. For instance, forRe = 654 ( Fig. 7 (b)) the vortex had already developed at the first half of the decelerating phase and propagated until the beginning of the accelerating phase, as suggested by the dashed lines. At this stage, the thickness of the recirculation zone is approximately equal to h. During the accelerating phase, the thickness of the recirculation zone began to decrease as that of the central jet increased. Toward the end of this phase, the central jet was thicker and the previous vortex had almost entirely dissipated, while a new vortex had begun to develop in the vicinity of the constriction (z/D < 0.2). ForRe = 1106 (Fig. 7 (c)), the flow presented the same characteristics as forRe = 654. Here, the vortex travelled faster during the entire cycle and the recirculation zone was 27% thicker. This is ascribed to the increment in the peak velocity of the central jet and the consequent increase in shear stress, leading to enlargement of the recirculation zone and an increase in radial velocity. AtRe = 654 andRe = 1106 two vortices were observed in the region of interest over one pulsatile period (Figs. 7(b) and (c)). In the case ofRe = 1106 ( Fig. 8 (b)), only one vortex was observed over one period. This difference arose because the vortex propagates at lower velocities for lower Reynolds values, then it is expected that two vortices will be observed in the region of interest, which were shed in consecutive pulsatile periods. ForRe = 1106 at the beginning of the accelerating phase, the central jet decreased in thickness and the recirculation zone became enlarged to occupy almost the entire tube section at z/D ≥ 1. This can be attributed to a marked deceleration of the central jet at the start of the accelerating phase. Due to mass conservation, the radial velocity is expected to increase, leading to enlargement of the recirculation zone. This is also consistent with that reported previously by Sherwin [23]. AtRe = 1767 (Fig. 8 (c)), the vortex developed and propagated faster than in the previous cases, and a secondary vortex was observed. In the decelerating phase, the vortex reached its maximum size and almost escaped from the region of interest. When the velocity of the central jet at the con-striction was near its minimum, the vortex left the region of interest and the rest of the flow became disordered. In the accelerating phase, the vortex developed earlier than at lower Reynolds numbers.
Comparison of Figs. 7 and 8 illustrates how the flow pattern changed with the degree of constriction. While in Fig. 7 the vortex travelled without a substantial change in size, Fig. 8 (b) shows a sharp change in vortex size, measured radially, and Fig. 8 (c) shows that the vortex had almost entirely dissipated at t = 0.5T . Moreover, due to the reduction in d 0 (Fig. 8) the velocity of the central jet was higher, and the recirculation zone became enlarged with increasing z, which explains the increase in vortex size. This is consistent with that reported by Sherwin et al. [23] and Usmani [26].
Enlargement of the recirculation zone was present throughout the experiments. This could be explained in terms of circulation, as pointed out in the work of Gharib et al. [32] which studies the flow led by a moving piston into an unbounded domain. This work determined that a vortex forms up to a limited amount of circulation before shedding from the layer where it was created. The main difference from our work consists in the confinement that the walls impose on the vortex. The result is that a vortex which sheds with a certain amount of circulation enlarges in the axial direction. This mechanism can also explain the generation of a secondary vortex. Specifically, the vortex sheds after a precise amount of circulation is reached. Then any excess of circulation generated goes to a trail of vorticity behind the vortex, eventually identified as a secondary vortex.
Finally, experimental and numerical results were compared in order to validate the simulation. Figure 9 shows the evolution of the simulated flow at Re = 654 andRe = 1106. Comparison of Fig. 9 with Fig. 7 confirms similar behaviour for the experimental and simulated flows. The main structures, i.e. the central high-velocity jet and recirculation zone, were satisfactorily reproduced by numerical simulation.

ii. Vortex propagation
Vortex propagation was studied by measuring vortex displacement as a function of time along one pulsatile period. To this end, the study region was constrained to 0.3 < r/D < 0.5 in order to isolate the vortex fraction forming in the superior wall of the tube. In this region vorticity values of the vortex are positive. For each time frame, the vorticity field was used to extract the vortex position. Then, in each frame, all vorticity values below a threshold of 30% of the maximum vorticity value were disregarded in order to obtain the filtered vorticity field ζ(r, z). The vortex position was then calculated as (r,z) = (r, z)ζ(r, z). The vortex propagation was along z and is specified in figs. 7, 8 and 9 by a black dashed line. This aforementioned method enables us to track the vortex and finally to measure the vortex maximum displacement, which will be discussed in the next subsections.

iii. Numerical results for varying Womersley number
Maximum vortex displacement was defined as z * /D. Position z * is where vorticity becomes lower than 30% of the maximum vorticity and is measured from z 0 , the position where the vortex formed. Position z 0 was defined as the location of the vortex centre before it sheds. The dependence of the vortex displacement over its lifetime on α, or f , was studied numerically for fixed values ofRe = 654 andRe = 1106, and pulsatile frequency values of 0.5 f , 0.75 f , f , 1.5 f and 2 f , where f is the pulsatile frequency tested experimentally. Figure 10 clearly illustrates the dependence of the maximum vortex displacement on α. For a fixed value ofRe , as α increased, i.e., f increased, the flow behaviour tended to replicate within a smaller region of the tube over a shorter time span. With increasing α, vortices tended to have a shorter lifespan and their maximum displacement was therefore smaller. In other words, z * /D decreased with increasing pulsatile frequency.

iv. Scaling law
A full description of vortex displacement has been given in previous sections. Specifically, the maximum displacement of the vortex was studied; that is, the distance it travels before it vanishes. From our results we propose a scaling law which summarizes the behavior of z * /D as a function of the parameters involved,Re and α.
The physical dimensions for the relevant variables are ν, D, f andv. Based on the Vaschy-Buckingham theorem, it is possible to describe z * /D as a function of two independent dimensionless numbers,Re and α, which involve the relevant variables mentioned. Maximum displacement z * /D was found to be proportional to flow velocitȳ v (i.e.Re ), and inversely proportional to pulsation frequency f (i.e. α 2 ), suggesting that where K is the proportionality constant. Figure 11 shows z * /D as a function ofRe/α 2 , showing that both the experimental data and the numerical results are a good fit to Eq. 1. The fitted slope was K = 1.15 ± 0.19 at a confidence level of 95%. This result highlights the dependence of z * /D onRe and 1/α 2 , as could be inferred from the data of Fig. 10.
The work of Gharib et al. [32] shows that the developing vortex reaches a threshold of circulation before it sheds. If an excess of circulation is generated, this could lead to the formation of other structures such as secondary vortices. This means that pulsatile frequency f coincides with the shedding frequency of the primary vortex, the one that was tracked. Hence, the parameterRe/α 2 can be related to the Strouhal number through Re/α 2 = 2 π d0 D 1 Sr . A discussion of the results in terms of the Strouhal number clarifies the relation between the oscillatory component of the flow and the stationary component of the flow. For lower values ofRe/α 2 , that is, Sr close to 1, stationary and oscillatory components are comparable, which explains a vortex with low to null displacement. The extreme case occurs when the flow oscillates but no vortices are shed. For instance, this could be attained for Reynolds numbers below those used in this work. AsRe/α 2 increases, that is for Sr 10 −1 , vortices are shed and travel further before vanishing as a consequence of the stationary component of the flow being greater than the oscillatory component. ParameterRe/α 2 also states the kinematic nature of vortex displacement, and by rewritingR e α 2 = 2d0v πD 2 f we conclude that z * /D depends on kinematic variables and not on viscosity.
Finally, this shows that Eq. 1 can be considered as a scaling law that describes the vortex displacement over its lifetime for any combination of the relevant parameters Re and α within the range studied and the constriction shape used. Specifically for the dependence on constriction shape, preliminary results with a Guassian-shaped constriction were obtained numerically, giving a difference below 20% for the maximun vortex displacement for Re u =1187. Vortex dynamics depend on inlet condition shape [23,33]. However, studies are not conclusive regarding the dependence of the vortex maximum displacement on the shape of the inlet condition. Further research is being carried out on this point.

IV. Conclusions
In this work, a pulsatile flow in an axisymmetric constricted geometry was studied experimentally and results were compared with those obtained numerically for the same values ofRe and α. After validation of the numerical results, simulations were run over a range of α values that could not be tested experimentally, in order to identify trends in flow behaviour with varying α.
The flow structure was found to consist of a cen-tral jet around the centreline and a recirculation zone adjacent to the wall, with vortices shedding in the latter. The analysis addressed how the vortex trajectory and size changed with the system parameters. Specifically, for a fixed value of α, vortex size grew with decreasing d 0 , and vortex displacement was larger for increasingRe values. The dependence on α was such that as α increased, the behavior of the flow was reproduced in a shorter extension of the tube. The vortex trajectory was tracked and its displacement over its lifetime deter- mined. Results showed that vortex displacement over its lifetime decreased with increasing α. The analysis led to a scaling law establishing linear dependence of the vortex displacement on a dimensionless parameter combiningRe and α, namelȳ Re/α 2 . Moreover, this parameter was found to be proportional to the inverse of the Strouhal number. This directly relates the vortex behavior to the ratio between the pulsatile component and the stationary component of the flow.
As seen from the medical perspective on the issue of stenosed arteries, these results provide insight into vortex shedding and displacement in a simplified model of stenosed arteries. This becomes crucial since vortex shedding precedes turbulence, and turbulence is related to plaque complications which may finally lead to cardiovascular stroke. The authors encourage further research into the behaviour of vortices over a wider range of parameters, including different constriction shapes and sizes.