Papers in Physics, vol. 12, art. 120002 (2020) Received: 20 February 2020, Accepted: 25 May 2020 Edited by: S. Rafai Reviewed by: J. A. Isler Licence: Creative Commons Attribution 4.0 DOI: https://doi.org/10.4279/PIP.120002 www.papersinphysics.org ISSN 1852-4249 Vortex dynamics under pulsatile flow in axisymmetric constricted tubes Nicasio Barrere1*, Javier Brum1, Alexandre L’her1, Gustavo L. Sarasúa1, Cecilia Cabeza1 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, in- cluding the transition to turbulent regime. For example, the characterization of pulsatile flows in obstructed artery models serves to encourage research into flow-induced phenom- ena 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 experi- ments 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. Cardiovascu- lar diseases are the primary cause of death world- wide (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 accumula- tion 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- * nbarrere@fisica.edu.uy [1] Instituto de F́ısica, Facultad de Ciencias, UdelaR, Mon- tevideo 11400, Uruguay cular stroke and myocardial infarct. Several fac- tors leading to plaque complications have been re- ported: sudden increase in luminal pressure [2], tur- bulent 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, under- standing 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–11]. Simplified models of stenosed ar- teries have used constricted rigid tubes [12–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 120002-1 Papers in Physics, vol. 12, art. 120002 (2020) / Nicasio Barrere et al. 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 us- ing one-point anemometry. However, one-point anemometry fails to account for different flow struc- tures 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 two- component laser Doppler anemometry. For con- striction degrees up to 50% of the lumen, the au- thors found no disturbances at Reynolds numbers below 1000. This study was not focused on describ- ing 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 pro- duced 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 Black- burn [23] studied the transition to turbulence over a wide range of Reynolds numbers. Mittal et al. used a planar channel with a one-sided semicircu- lar 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 constric- tion. Based on the results of linear stability analy- sis, 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] com- pared numerical results with those obtained by hot-wire measurements. As mentioned above, one- dimensional hot-wire measurements cannot be used to identify flow structures. Griffith et al. [25], using a rigid tube with a slightly narrowed section resem- bling stenosis, compared numerical results with ex- perimental 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 vor- tices formed immediately downstream of the steno- sis 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 Womers- ley between 6 and 8. The authors reported that the downstream flow was characterized by a high veloc- ity 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 ulti- mately, 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 num- bers varying from 385 to 2044 and Womersley num- bers from 17 to 50. The results show that the flow pattern in these systems consists primarily of a cen- tral jet and a vortex shedding layer adjacent to the wall (recirculation flow), consistent with the litera- ture. By tracking the vortex trajectory it was pos- sible 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. 120002-2 Papers in Physics, vol. 12, art. 120002 (2020) / Nicasio Barrere et al. Figure 1: : (a) Schematic drawing of the experimental setup; (b) lengthwise view of a tube of diameter D = 2.6 ± 0.1 cm and an annular constriction of diameter d0; (c) experimental axial velocity at z = 0 and r = 0 for R̄e = 1106, with d0 = 1.6 cm . II. Materials and methods 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 circu- lating 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) con- sisted 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 di- ameter was changed to achieve different degrees of constriction. Tests were carried out using annular constrictions with inner diameters d0 = 1.6 ± 0.1 cm and d0 = 1.0 ± 0.1 cm, equivalent to degrees of constriction relative to D of 39% and 61%, re- spectively. A cross-sectional scheme of the constric- tion 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 ax- ial 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 val- ues. Finally, this constricted tube (CT) was placed inside a chamber filled with water, so that the re- fraction 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 in- ner 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. Im- ages were obtained at a frame rate of 180 Hz over a period equivalent to 16 cycles, using a CMOS cam- era (Pixelink, PL-B776F). The velocity field was finally computed using OpenPIV open-source soft- ware with 32 × 32 pixel2 windows and an overlap of 8 pixels in both directions. 120002-3 Papers in Physics, vol. 12, art. 120002 (2020) / Nicasio Barrere et al. 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 spec- tral analysis of the velocity fields. This observa- tion 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 experi- ments 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.1 (c), remained unchanged throughout the experiments. Pulsation period and peak velocity can be related to the Womersley and Reynolds numbers, respectively. The Womersley number, α, is defined as α = D 2 √ 2πf ν , where f is the pulsatile frequency and ν the kinematic viscos- ity of water (ν = 1.0 × 10−6 m2/s). The experiments were carried out for three dif- ferent 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 Reu = Dvu/ν where vu is the peak velocity measured at the centreline. Four different values of Reu were used as inlet con- dition upstream of the constriction (see table 1). Similarly, the mean Reynolds at the downstream edge of the constriction is defined as R̄e = d0v̄/ν, where v̄ is the mean velocity over a whole period, measured at (z,r) = (0, 0). Finally, experiments are labeled through the R̄e values, since this rep- resents a unique combination of Reu, d0 and α, as α R̄e Reu Constriction d0 (cm) 33.26 654 820 39% 1.6 1002 820 61% 1.0 1106 1187 39% 1.6 1767 1187 61% 1.0 2044 1625 61% 1.0 21.10 385 505 39% 1.6 963 820 39% 1.6 17.22 585 820 39% 1.6 Table 1: Experimental degrees of constriction and R̄e numbers. Figure 2: Normalized displacement of piston yp for Reu=505 (diamond green), Reu=820 (blue circles), Reu=1187 (red squares), Reu=1625 (black plus sign). 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 mo- tor driving a piston inside a cylindrical chamber. The control module contained the power supply for the entire system and an electronic board, pro- grammable from a remote PC by means of custom software. The piston motion was programmed us- ing the pump settings for pulsatile period, ejected volume, pressure and cycle shape. The piston mo- tion of the programmable pump is shown in Fig. 2 in normalized units of displacement and time. Fig- ure 2 shows that the piston motion, and hence the pulsatile waveform, are comparable among all Reu values, meaning that the inlet conditions are equiv- alent 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. Ac- cording 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, R̄eu, giving L ≈ 0.05R̄euD. In our setup, the flow de- velopment section measures 48D, R̄eu <1000 and α >10 for all experiments, which makes it rea- sonable to assume a fully-developed flow. Never- theless, velocity profiles were studied for Reu val- ues of 820, 1187 and 1625 in order to ensure de- veloped flow condition. Profiles are taken from a 120002-4 Papers in Physics, vol. 12, art. 120002 (2020) / Nicasio Barrere et al. Figure 3: Upstream velocity profiles for (a) Reu=820, (b) Reu=1187 and (c) Reu=1625. Profiles are located at z1=-1.25D (blue circles), z2=-0.75D (red triangles) and z3=-0.25D (black plus sign). 1.5D length window upstream of the constriction, at fixed locations z1 = −1.25D, z2 = −0.75D and z3 = −0.25D. Velocity values are normalized by the maximum velocity attained in a pulsatile cy- cle. Figure 3 shows such profiles at t = 0.25T , t = 0.5T , t = 0.75T and t = T . Profiles show that for each Reu, independently of the time, profiles are equal for z1=-1.25D, z2=-0.75D and z3=-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 param- eter 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 Reu=820 and Reu=1187 experi- mental 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 vz(t) = a0 2 + ∑4 n=1 (ancos (2πnft) + bnsin (2πnft)) with coefficients (×10−5, in cm/s) a0 =196.32, a1 =61.94, a2 =8.22, a3 =13.79, a4 =3.87, b1 =-27.66, b2 =-3.81, b3 =-5.13, b4 =-2.84 for Reu=1187. Higher-order terms were not consid- ered, as they do not contribute significantly to the shape of velocity profiles. In order to improve nu- merical stability, the velocity profile calculated by decomposition was multiplied by a ramp function of 120002-5 Papers in Physics, vol. 12, art. 120002 (2020) / Nicasio Barrere et al. time increasing steadily from 0 to 1 in 0.6 seconds. Figure 4(a) shows the axial velocity at the centre- line at the downstream edge of the constriction, and presents good agreement with its experimen- tal 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 sur- face 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 sim- ulations 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 ele- ment size to adapt to the physics at specific re- gions. 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 re- sults. Profiles are located upstream (z = −0.3D) and downstream (z = 0.75D) for the two inlet con- ditions Reu=820 and Reu=1187. The maximum difference between numerical and experimental pro- files 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 (Reu=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 com- putational time. III. Results and discussion i. Flow pattern The velocity at the centreline, r = 0, at the down- stream edge of the constriction, z = 0, as a func- tion of time (see Fig. 1(c)) was used to establish the time reference for all the experiments. The de- celerating phase (diastolic phase) was defined as 0 < t < 0.5T and the accelerating phase (systolic phase) as 0.5T < t < T. Figures 7 and 8 show the temporal flow patterns obtained from the experimental data correspond- ing to R̄e=654, 1002, 1106, 1767 and 2044, while Fig. 9 shows the patterns obtained via simulation. Figures 7(a), 8(a) and 9(a) show the flow evolution at times 0.25 T , 0.5 T , 0.6 T and 0.9 T. Some features were common to all the experiments and numerical results. During the accelerating phase two main flow structures can be distinguished: a central, high-velocity central jet and a recirculation zone between the former and the inner surface of the tube wall, with vortices developing in the latter. In Figs. 7, 8 and 9, the central jet is shown in blue and the recirculation zone in red. The central jet is distinguished from the recirculation zone by having vorticity values below the threshold of 30% of the maximum of the absolute value of vorticity in each timestep. The aim of this criterion is to separate regions approximately rather than give their pre- cise location. Within the recirculation zone, vor- tices separated from the wall and travelled along the tube. Other features are dependent on the Reynolds number. For instance, for R̄e = 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). For R̄e = 1106(Fig. 7 (c)), the flow presented the same characteristics as for R̄e = 654. Here, the vortex travelled faster during the entire cycle 120002-6 Papers in Physics, vol. 12, art. 120002 (2020) / Nicasio Barrere et al. Figure 4: (a) Axial inlet velocity at r = 0; corresponds to Reu = 1187; (b) axial velocity at r = 0 and z = 0; corresponds to R̄e = 1106; (c) schematic drawing of the flow development section and the axisymmetric constricted tube; (d) cross sectional view of the mesh. 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. At R̄e = 654 and R̄e = 1106 two vortices were ob- served in the region of interest over one pulsatile pe- riod (Figs. 7(b) and (c)). In the case of R̄e = 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. For R̄e = 1106 at the beginning of the accelerating phase, the central jet decreased in thickness and the re- circulation 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 recircula- tion zone. This is also consistent with that reported previously by Sherwin [23]. At R̄e = 1767 (Fig. 8 (c)), the vortex developed and propagated faster than in the previous cases, and a secondary vortex was observed. In the de- celerating 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 constric- tion. 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 en- tirely dissipated at t = 0.5T. Moreover, due to the reduction in d0 (Fig. 8) the velocity of the central jet was higher, and the recirculation zone became enlarged with increasing z, which explains the in- crease 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 do- main. 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 dif- ference 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 sec- ondary vortex. Specifically, the vortex sheds after 120002-7 Papers in Physics, vol. 12, art. 120002 (2020) / Nicasio Barrere et al. Figure 5: Behavior comparison for three meshes: 180958 elements (Mesh 1), 351569 elements (Mesh 2) and 951899 elements (Mesh 3). (a) Axial velocity at r = 0, z = 0 vs time and (b) Velocity profiles in down- stream location z = 0.5D at 0.25T, 0.5T, 0.75T and T. 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. Fig- ure 9 shows the evolution of the simulated flow at R̄e = 654 and R̄e = 1106. Comparison of Fig. 9 with Fig. 7 confirms similar behaviour for the ex- perimental and simulated flows. The main struc- tures, i.e. the central high-velocity jet and recircu- lation zone, were satisfactorily reproduced by nu- merical simulation. ii. Vortex propagation Vortex propagation was studied by measuring vor- tex 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 vor- tex 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 disre- garded 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 mea- sure the vortex maximum displacement, which will be discussed in the next subsections. iii. Numerical results for varying Womers- ley 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 z0, the position where the vortex formed. Position z0 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 of R̄e = 654 and R̄e = 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 of R̄e , 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 de- creased with increasing pulsatile frequency. iv. Scaling law A full description of vortex displacement has been given in previous sections. Specifically, the maxi- mum displacement of the vortex was studied; that is, the distance it travels before it vanishes. From 120002-8 Papers in Physics, vol. 12, art. 120002 (2020) / Nicasio Barrere et al. Figure 6: Experimental (blue circles) and numerical (red solid line) velocity profiles at upstream location z = −0.3D for (a) Reu=820, (b) Reu=1187 and at downstream location z = 0.75D for (c) Reu=820, (d) Reu=1187 our results we propose a scaling law which sum- marizes the behavior of z∗/D as a function of the parameters involved, R̄e and α. The physical dimensions for the relevant vari- ables are ν, D, f and v̄. Based on the Vaschy- Buckingham theorem, it is possible to describe z∗/D as a function of two independent dimension- less numbers, R̄e and α, which involve the rele- vant variables mentioned. Maximum displacement z∗/D was found to be proportional to flow velocity v̄ (i.e. R̄e ), and inversely proportional to pulsation frequency f (i.e. α2), suggesting that z∗/D = K R̄e α2 (1) where K is the proportionality constant. Figure 11 shows z∗/D as a function of R̄e/α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 on R̄e 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 circula- tion 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 parameter R̄e/α2 can be related to the Strouhal number through R̄e/α2 = 2 π d0 D 1 Sr . A discussion of the results in terms of the Strouhal number clarifies the rela- tion between the oscillatory component of the flow and the stationary component of the flow. For lower values of R̄e/α2, that is, Sr close to 1, sta- tionary and oscillatory components are compara- ble, which explains a vortex with low to null dis- placement. The extreme case occurs when the flow oscillates but no vortices are shed. For instance, this could be attained for Reynolds numbers below 120002-9 Papers in Physics, vol. 12, art. 120002 (2020) / Nicasio Barrere et al. Figure 7: (a) Axial velocity profile at z = 0 and r = 0 and time at which the velocity field was obtained (red circles); (b) streamlines derived from experimental data for R̄e = 654, with d0 = 1.6 cm; (c) streamlines derived from experimental data for R̄e = 1106, with d0 = 1.6 cm. In both cases, red streamlines represent the recirculation zone and blue streamlines the central jet (colours available online only). those used in this work. As R̄e/α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 os- cillatory component. Parameter R̄e/α2 also states the kinematic nature of vortex displacement, and by rewriting R̄e α2 = 2d0v̄ πD2f we conclude that z∗/D depends on kinematic variables and not on viscos- ity. Finally, this shows that Eq. 1 can be consid- ered as a scaling law that describes the vortex dis- placement over its lifetime for any combination of the relevant parameters Re and α within the range studied and the constriction shape used. Specifi- cally for the dependence on constriction shape, pre- liminary results with a Guassian-shaped constric- tion were obtained numerically, giving a difference below 20% for the maximun vortex displacement for Reu=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. 120002-10 Papers in Physics, vol. 12, art. 120002 (2020) / Nicasio Barrere et al. Figure 8: (a) Axial velocity profile at z = 0 and r = 0 and time at which the velocity field was obtained (red circles); (b) streamlines derived from experimental data for R̄e = 1002, with d0 = 1.0 cm; (c) streamlines derived from experimental data for R̄e = 1767, with d0 = 1.0 cm. In both cases, red streamlines represent the recirculation zone and blue streamlines the central jet (colours available online only). IV. Conclusions In this work, a pulsatile flow in an axisymmetric constricted geometry was studied experimentally and results were compared with those obtained nu- merically for the same values of R̄e and α. Af- ter 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 vor- tex trajectory and size changed with the system parameters. Specifically, for a fixed value of α, vor- tex size grew with decreasing d0, and vortex dis- placement was larger for increasing R̄e 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- 120002-11 Papers in Physics, vol. 12, art. 120002 (2020) / Nicasio Barrere et al. Figure 9: (a) Axial velocity profile at z = 0 and r = 0 and time at which the velocity field was obtained (red circles); (b) streamlines derived from simulation results for R̄e = 654; (c) streamlines derived from simulation results for R̄e = 1106. In both cases, red streamlines represent the central jet layer and blue streamlines represent the central jet (colours available online only). mined. Results showed that vortex displacement over its lifetime decreased with increasing α. The analysis led to a scaling law establishing linear de- pendence of the vortex displacement on a dimen- sionless parameter combining R̄e and α, namely R̄e/α2 . Moreover, this parameter was found to be proportional to the inverse of the Strouhal num- ber. 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 simpli- fied model of stenosed arteries. This becomes cru- cial since vortex shedding precedes turbulence, and turbulence is related to plaque complications which may finally lead to cardiovascular stroke. The au- thors encourage further research into the behaviour of vortices over a wider range of parameters, includ- ing different constriction shapes and sizes. 120002-12 Papers in Physics, vol. 12, art. 120002 (2020) / Nicasio Barrere et al. Figure 10: Dependence of the dimensionless maximum vortex displacement on α. Experimental results in blue squares and numerical results in red circles (R̄e = 654) and red triangles (R̄e = 1106). Figure 11: Dimensionless maximum vortex displace- ment, z∗/D, as a function of the dimensionless param- eter R̄e/α2 for experimental data (blue squares) and numerical data (red circles). Acknowledgements This research was supported by CSIC I+D Uruguay, through the I+D project 2016 “Estu- dio dinámico de un flujo pulsátil y sus impli- caciones hemodinámicas vasculares”, ANII (Doc- toral scholarship POSNAC-2015-1-109843), ECOS- SUD (project reference number U14S04) and PEDECIBA, Uruguay. The manuscript was edited by Eduardo Speranza. References [1] M Naghavi et al., From Vulnerable Plaque to Vulnerable Patient: a call for new definitions and risk assessment strategies: Part I, Circu- lation 108, 1664 (2003). [2] J E Muller, G H Tofler, P H Stone, Circadian variation and triggers of onset of acute cardio- vascular disease, Circulation 79, 733 (1989). [3] H M Loree, R D Kamm, C M Atkinson, R T Lee, Turbulent pressure fluctuations on surface of model vascular stenosis, Am. J. of Physiol. 261, 644 (1991). [4] S D Gertz, W C Roberts, Hemodynamic shear force in rupture of coronary arterial atherso- clerosis plaques, Am. J. of Cardiol. 66, 1368 (1990). [5] A C Barger, R Beeuwkes 3rd, L L Lainey, K J Silverman, Hypothesis: vasa vasorum and neo- vascularization of human coronary arteries. A possible role in the pathophysiology of athe- riosclerosis, N. Engl. J. Med 19. 175 (1984). [6] E Falk, P K Shah, V Fuster, Coronary plaque disruption, Circulation 92, 657 (1995). [7] K Imoto textitet al., Longitudinal Structural Determinants of Atherosclerotic Plaque Vul- nerability, J. Am. Coll. Cardiol. 46, 1507 (Oct. 2005). [8] C G Caro, Atheroma and arterial wall shear observations, correlation and proposal of a shear dependent mass transfer mechanism for atherogenesis, Proc. R. Soc. Lond. B 17, 109 (1971). [9] D N Ku, D P Giddens, C K Zarins, S Glagov, Pulsatile flow and atherosclerosis in the hu- man carotid bifurcation, Arteriosclerosis 5, 293 (1985). [10] R M Nerem, M J Levesque, Fluid mechanics in atherosclerosis, In: Handbook of bioengineer- ing, Eds. R Skalak, S Chien, Chap. 21, Pag. 21.1, McGraw-Hill, New York (1987). [11] S S Gopalakrishnan, B Pier, A Biesheuvel, Dy- namics of pulsatile flow through model abdom- inal aortic aneurysms, J. Fluid Mech. 758, 150 (2014). 120002-13 https://doi.org/10.1161/01.CIR.0000087480.94275.97 https://doi.org/10.1161/01.CIR.0000087480.94275.97 https://doi.org/10.1161/01.CIR.79.4.733 https://doi.org/10.1152/ajpheart.1991.261.3.H644 https://doi.org/10.1152/ajpheart.1991.261.3.H644 https://doi.org/10.1016/0002-9149(90)91170-B https://doi.org/10.1016/0002-9149(90)91170-B https://doi.org/10.1056/NEJM198401193100307 https://doi.org/10.1161/01.CIR.92.3.657 https://doi.org/10.1016/j.jacc.2005.06.069 https://doi.org/10.1016/j.jacc.2005.06.069 https://doi.org/10.1098/rspb.1971.0019 https://doi.org/10.1098/rspb.1971.0019 https://doi.org/10.1161/01.ATV.5.3.293 https://doi.org/10.1161/01.ATV.5.3.293 https://doi.org/10.1017/jfm.2014.535 https://doi.org/10.1017/jfm.2014.535 Papers in Physics, vol. 12, art. 120002 (2020) / Nicasio Barrere et al. [12] A Arzani, P Dyverfeldt, T Ebbers, S C Shad- den, In vivo validation of numerical prediction for turbulence intensity in aortic coarctation, Ann. Biomed. Eng. 40, 860 (2012). [13] M D Ford et al., Virtual angiography for visu- alization and validation of computational mod- els of aneurysm hemodynamics, IEEE T. Med. Imaging 24, 1586 (2005). [14] L Boussel et al., Phase contrast magnetic res- onance imaging measurements in intracranial aneurysms in vivo of flow patterns, velocity fields, and wall shear stress: Comparison with computational fluid dynamics, Magn. Reson. Med. 61, 409 (2009). [15] Q Long, X Y Xu, K V Rammarine, P Hoskins, Numerical investigation of physiologically real- istic pulsatile flow through arterial stenosis, J. Biomech. 34, 1229 (2001). [16] J C Stettler, A K M Fazle Hussain, On transi- tion of pulsatile pipe flow, J. Fluid Mech. 170, 169 (1986). [17] J Peacock, T Jones, R Lutz, The onset of turbulence in physiological pulsatile flow in a straight tube, Exp. Fluids 24, 1 (1998). [18] C S Chua et al., Particle image of non axisym- metic stenosis model, 8th International Sym- posium on Particle Image Velocimetry (2009). [19] S A Ahmed, D P Giddens, Flow disturbance measurements through a constricted tube at moderate Reynolds numbers, J. Biomech. 16, 955 (1983). [20] S A Ahmed, D P Giddens, Pulsatile post- stenotic flow studies with laser Doppler- anemometry, J. Biomech. 17, 695 (1984). [21] J A Isler, R S Gioria, B S Carmo, Bifurcations and convective instabilities ofsteady flows in a constricted channel, J. Fluid Mech. 849, 777 (2018). [22] R Mittal, S P Simmons, F Najjar, Numerical study of pulsatile flow in a constricted channel, J. Fluid Mech. 485, 337 (2003). [23] S J Sherwin, H M Blackburn, Three dimen- sional instabilities and transition of steady and pulsatile axisymmetric stenotic flows, J. Fluid Mech. 533, 297 (2005). [24] S C Ling, H B Atabek, A nonlinear analysis of pulsatile flow in arteries, J. Fluid Mech. 55, 493 (1972). [25] M D Griffith, T Leweke, M C Thompson, K Hourigan, Pulsatile flow in stenotic geome- tries: flow behaviour and stability, J. Fluid Mech. 622, 291 (2009). [26] A Y Usmani, K Muralidhar, Pulsatile flow in a compliant stenosed asymmetric model, Exp. Fluids 57, 186 (2016). [27] J Westerweel, Fundamentals of digital particle velocimetry, Meas. Sci.Technol. 8, 1379 (1997). [28] M Raffel, C Willert, S Wereley, J Kompen- hans, Particle image velocimetry. A practical guide, Second edition, Springer (2007). [29] R A Cassanova, D P Giddens, Disorder dis- tal to modeled stenoses in steady and pulsatile flow, J. Biomech. 11, 441 (1977). [30] G Balay, Elasticidad en tejidos arteriales, diseno de un corazon artficial in vitro y nuevo metodo ultrasonico de determinacion de elas- ticidad arterial, MA thesis, Universidad de la Republica, Uruguay, 2012. [31] D N Ku, Blood flow in arteries, Ann. Rev. Fluid Mech. 29, 399 (1997). [32] M Gharib, E Rambod, K Shariff, A universal time scale for vortex ring formation, J. Fluid Mech. 360, 121 (1998). [33] H Liu, T Yamaguchi, Waveform Dependence of Pulsatile Flow in a Stenosed Channel, J. Biomech. Eng. 123, 88 (2000). 120002-14 https://doi.org/10.1007/s10439-011-0447-6 https://doi.org/10.1109/TMI.2005.859204 https://doi.org/10.1109/TMI.2005.859204 https://doi.org/10.1002/mrm.21861 https://doi.org/10.1002/mrm.21861 https://doi.org/10.1016/S0021-9290(01)00100-2 https://doi.org/10.1016/S0021-9290(01)00100-2 https://doi.org/10.1017/S0022112086000848 https://doi.org/10.1017/S0022112086000848 https://doi.org/10.1007/s003480050144 https://doi.org/10.1016/0021-9290(83)90096-9 https://doi.org/10.1016/0021-9290(83)90096-9 https://doi.org/10.1016/0021-9290(84)90123-4 https://doi.org/10.1017/jfm.2018.410 https://doi.org/10.1017/jfm.2018.410 https://doi.org/10.1017/S002211200300449X https://doi.org/10.1017/S0022112005004271 https://doi.org/10.1017/S0022112005004271 https://doi.org/10.1017/S0022112072001971 https://doi.org/10.1017/S0022112072001971 https://doi.org/10.1017/S0022112008005338 https://doi.org/10.1017/S0022112008005338 https://doi.org/10.1007/s00348-016-2274-x https://doi.org/10.1007/s00348-016-2274-x https://doi.org/10.1088/0957-0233/8/12/002 https://doi.org/10.1007/978-3-540-72308-0 https://doi.org/10.1016/0021-9290(78)90056-8 https://doi.org/10.1146/annurev.fluid.29.1.399 https://doi.org/10.1146/annurev.fluid.29.1.399 https://doi.org/10.1017/S0022112097008410 https://doi.org/10.1017/S0022112097008410 https://doi.org/10.1115/1.1339818 https://doi.org/10.1115/1.1339818 Introduction Materials and methods Experimental setup Numerical simulation Results and discussion Flow pattern Vortex propagation Numerical results for varying Womersley number Scaling law Conclusions