# Vortex dynamics under pulsatile flow in axisymmetric constricted tubes

Nicasio Barrere, Javier Brum, Alexandre L’Her, Gustavo L. Sarasúa, Cecilia
Cabeza  Instituto de Físca, Facultad de Ciencias, UdelaR, Uruguay

###### Abstract

An improved understanding of how vortices develop and propagate under
pulsatile flow can shed important light on the mixing and transport processes
including the transition to turbulent regime occurring in such systems. 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 constriction degree ($d_{0}$), mean Reynolds number
($\bar{Re}$) and Womersley number ($\alpha$). Velocity fields were acquired
experimentally using Digital Particle Image Velocimetry and generated
numerically. For the acquisition of data, $\bar{Re}$ was varied from 385 to
2044, $d_{0}$ was 1.0 cm and 1.6 cm, and $\alpha$ was varied from 17 to 33 in
the experiments and from 24 to 50 in the numerical simulations. Results for
the considered Reynolds number, 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 the maximum vortex displacement as a function of a dimensionless
variable dependent on the system parameters Re and $\alpha$.

PACS

47.20.Ib Instability of boundary layers; separation, 47.32.cb Vortex
interactions

###### pacs:

47.20.Ib ,47.32.cb

##  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) Naghavi  _et al._ (2003). 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
cardiovascular stroke and myocardial infarct. Several factors leading to
plaque complications have been reported: sudden increase in luminal pressure
Muller  _et al._ (1989), turbulent fluctuations Loree  _et al._ (1991),
hemodynamic shear stress Gertz and Roberts (1990), vasa vasorum rupture Barger
_et al._ (1984), material fatigue Falk  _et al._ (1995) and stress
concentration within a plaque Imoto  _et al._ (2005). As most of these factors
are flow-dependent, understanding how a pulsatile flow behaves through a
narrowed region can provide important insight toward the development of
reliable diagnostic tools.

Flow alterations due to arterial obstruction (i.e. stenosis) and aneurysms
have been reported in the literature Caro  _et al._ (1971); Ku  _et al._
(1985); Nerem and Levesque (1987); Gopalakrishnan  _et al._ (2014). Simplified
models of stenosed arteries have used constricted rigid tubes Arzani  _et al._
(2012); Ford  _et al._ (2005); Boussel  _et al._ (2009); Long  _et al._
(2001). 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 Stettler
and Hussain (1986) 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 Peacock  _et al._ (1998) studied the
transition to turbulence in a rigid tube using two-dimensional particle
imaging velocimetry (PIV). Chua et al Chua  _et al._ (2009) implement a three-
dimensional PIV technique to get a volumetrical velocity field of a steady
flow. The work of Ahmed and Giddens Ahmed and Giddens (1983, 1984) studied the
transition to turbulence in a constricted rigid tube with varying degree of
constriction using two-component 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. Long  _et al._
(2001) 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. Isler  _et al._ (2018) in an axisymmetric constricted channel
found the instabilities that breaks the symmetry of the flow. The works of
Mittal et al. Mittal  _et al._ (2003) and Sherwin and Blackburn Sherwin and
Blackburn (2005) 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 another 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 Sherwin and Blackburn (2005) 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 instabilities involving vortex shedding take place in the
transition to turbulence.

Finally, several studies compared experimental and numerical results. Ling et
al. Ling and Atabek (1972) compared 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.
Griffith  _et al._ (2009), 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 Usmani and Muralidhar (2016) compares flow patterns in rigid and
compliant asymmetric constricted tubes for Reynolds range between 300 and 800
and Womersley between 6 and 8. The authors report that the downstream flow is
characterized by a high velocity jet and a vortex whose evolution is
qualitatively described. The aforementioned studies highlight the significance
of study the 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 the 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 varying
mean Reynolds number from 385 to 2044 and Womersley number from 17 to 50.
These 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. The vortex kinematics was described as a function of the
system parameters in the form of a dimensionless scaling law for the maximum
vortex displacement.

##  II Methods and materials

###  II.1 Experimental setup

The experimental setup (Fig. 1 (a)) consisted of a circulating loop and a
Digital Particle Imaging Velocimetry (DPIV) module capable of acquiring 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 with $L=51.0\pm 0.1$ cm in length and inner diameter $D=2.6\pm 0.1$
cm. The annular constriction consisted of a hollow cylinder $5\pm 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\pm 0.1$ cm and $d_{0}=1.0\pm
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.

Figure 1: : (a) Schematic drawing of the experimental setup; (b) lengthwise
view of a tube of diameter $D=2.6\pm 0.1$ cm and an annular constriction of
diameter $d_{0}$; (c) experimental axial velocity at $z=0$ and $r=0$ for
$\bar{Re}$ = 1002, with $d_{0}$ = 1.6 cm .

The tube inlet was connected to the pulsatile pump via a flow development
section (FDS), whereas 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
with inner diameter increasing 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 48$D$ connected to the constricted tube (CT).
The reservoir (R) was used to set the minimum pressure inside the tube, which
was set 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 $\mu$m diameter, DANTEC). DPIV
technique was used to acquire the velocity field Westerweel (1997); Raffel
_et al._ (2007). A 1 W Nd:YAG laser was used to illuminate a 2 mm-thick
section of the tube. Images were acquired 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\times 32$ pixel2 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 Peacock  _et al._ (1998); Ling and Atabek (1972); Cassanova and
Giddens (1977). 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 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, $\alpha$, is defined as
$\alpha=\frac{D}{2}\sqrt{\frac{2\pi f}{\nu}}$, where $f$ is the pulsatile
frequency and $\nu$ the kinematic viscosity of water ($\nu=1.0\times 10^{-6}$
m2/s).

The experiments, were carried out for three different pulsatile period $T$:
0.96s, 2.39s and 3.58s corresponding to $\alpha$ values of 33.26, 21.10 and
17.22 respectively. The peak Reynolds number upstream of the constriction is
defined as $Re_{\mathrm{u}}=Dv_{\mathrm{u}}/\nu$ where $v_{\mathrm{u}}$ is the
peak velocity measured at the centreline. Four different values of
$Re_{\mathrm{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 as $\bar{Re}=d_{0}\bar{v}/\nu$, where $\bar{v}$ is the
mean velocity over a whole period, measured at $(z,r)=(0,0)$. Finally,
experiments are labeled through its $\bar{Re}$ values since it represents a
unique combination of $Re_{\mathrm{u}}$, $d_{0}$ and $\alpha$ as summarized in
table 1.

$\bm{\alpha}$ | $\bm{\bar{Re}}$ | $Re_{\mathrm{u}}$ | Constriction |  $d_{0}$ (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 $\bar{Re}$ numbers.

The pulsatile pump (PP in Fig. 1 (a)) generated different flow conditions at
the tube inlet, Balay (2012). 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 all the system and an electronic board,
programmable from a remote PC by means of a 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 showed
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_{\mathrm{u}}$ values, meaning that the inlet conditions are equivalent for
all the experiments.

Figure 2: Normalized displacement of piston $y_{p}$ for $Re_{\mathrm{u}}$=505
(diamond green), $Re_{\mathrm{u}}$=820 (blue circles), $Re_{\mathrm{u}}$=1187
(red squares), $Re_{\mathrm{u}}$=1625 (black plus sign).

Finally, the condition of fully developed flow at the upstream region was
verified. The entrance length is usually estimated as $L\approx 0.05ReD$.
According to the work of Ku Ku (1997), the entrance length of a pulsatile flow
with $\alpha>$10, corresponds to the upstream mean Reynolds over a whole
cycle, $\bar{Re}_{\mathrm{u}}$, giving $L\approx 0.05\bar{Re}_{\mathrm{u}}D$.
In our setup, the flow development section measures 48D,
$\bar{Re}_{\mathrm{u}}<$1000 and $\alpha>$10 for all experiments, which makes
reasonable assuming developed flow. Nevertheless, velocity profiles were
studied for $Re_{\mathrm{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 shows that for each
$Re_{\mathrm{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.

Figure 3: Upstream velocity profiles for (a) $Re_{\mathrm{u}}$=820, (b)
$Re_{\mathrm{u}}$=1187 and (c) $Re_{\mathrm{u}}$=1625. Profiles are located at
$z_{1}$=-1.25D (blue circles), $z_{2}$=-0.75D (red triangles) and
$z_{3}$=-0.25D (black plus sign).

###  II.2 Numerical simulation

Three-dimensional numerical simulation using COMSOL™ software enabled a
comparison with 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 experimental
data in the upstream region, where the flow is developed. In order to do so, a
fourth-order Fourier decomposition of the experimental flow rate was carried
out. Higher-order terms were neglected, 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 a good agreement with its experimental analogue
shown in Fig. 1 (c). Outlet conditions were set to $p=0$, being $p$ 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 neglect
transient effects. The step size for the data output was selected to be the
same of the experiments.

Figure 4: For $\bar{Re}$ = 1002: (a) axial velocity at $r$ = 0 and $z$ = 0;
(b) schematic drawing of the flow development section and the axisymmetric
constricted tube; (c) cross sectional view of the mesh.

The mesh used was a free tetrahedral mesh composed of 351569 elements in
total, with an average element quality of 0.6244 and an element volume ratio
of $6.358\times 10^{-4}$. It was built using a COMSOL predefined physics-
controlled mesh, enabling the normal element size settings so as to allow the
mesh to adapt to the physics at specific regions. The mesh was found to be
adequate by comparing results to those obtained in a mesh of 951899 elements
for the inlet conditions given by $Re_{\mathrm{u}}$=820 and
$Re_{\mathrm{u}}$=1187. A cross-sectional view of the tube and the mesh
elements is shown in Fig. 4 (c).

##  III Results and Discussion

###  III.1 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$.

Figure 5: (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 $\bar{Re}$ = 654, with $d_{0}$ = 1.6 cm; (c) streamlines
derived from experimental data for $\bar{Re}$ = 1002, with $d_{0}$ = 1.6 cm.
In both cases, red streamlines represent the recirculation zone and blue
streamlines the central jet (colours available online only).

Figure 6: (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 $\bar{Re}$ = 1106, with $d_{0}$ = 1.0 cm; (c)
streamlines derived from experimental data for $\bar{Re}$ = 1767, with $d_{0}$
= 1.0 cm. In both cases, red streamlines represent the recirculation zone and
blue streamlines the central jet (colours available online only).

Figures 5 and 6 show the temporal flow patterns obtained from the experimental
data corresponding to $\bar{Re}$=654, 1002, 1106, 1767 and 2044, while Fig. 7
shows the patterns obtained via simulation. Figures 5(a), 6(a) and 7(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. 5, 6 and 7, the central jet is shown in blue and the recirculation zone
in red. The central jet is distinguished from recirculation zone for having
vorticity values below the threshold of 30% of the maximum of the absolute
value of vorticity in each timestep. The aim of this criteria, is to
approximately separate regions rather than giving a precise location of such
regions. Within the recirculation zone, vortices separated from the wall and
travelled along the tube.

Other features are dependent on the Reynolds number. For instance, for
$\bar{Re}$ = 654 (Fig. 5 (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 nearly entirely
dissipated, while a new vortex had begun to develop in the vicinity of the
constriction ($z/D<0.2$).

For $\bar{Re}$ = 1002(Fig. 5 (c)), the flow presented the same characteristics
as for $\bar{Re}$ = 654. Here, the vortex travelled faster during the whole
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 an enlargement of the recirculation zone and an
increase in radial velocity.

At $\bar{Re}$ = 654 and $\bar{Re}$ = 1002 two vortices were observed in the
region of interest over one pulsatile period, Figs. 5(b) and (c). In the case
of $\bar{Re}$ = 1106 (Fig. 6 (b)), only one vortex was observed over one
period. This difference arise because vortex propagates at lower velocities
for lower Reynolds values, then is expected to observe two vortex in the
region of interest which were shed in consecutive pulsatile periods. For
$\bar{Re}$ = 1106 at the beginning of the accelerating phase, the central jet
decreased in thickness and the recirculation zone has become enlarged to
occupy nearly the entire tube section at $z/D\geq 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 the enlargement of the recirculation zone. This is also consistent
with that reported previously by Sherwin Shadden  _et al._ (2005).

At $\bar{Re}$ = 1767 (Fig. 6 (c)), the vortex developed and propagated faster
than in the previous cases, and a secondary vortex is 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
constriction 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.

A comparison of Figs. 5 and 6 illustrates how the flow pattern changed with
the degree of constriction. While in Fig. 5 the vortex travelled without a
substantial change in size, Fig. 6 (b) shows a sharp change in vortex size,
measured radially, and Fig. 6 (c) shows that the vortex had nearly entirely
dissipated at $t=0.5T$. Moreover, due to the reduction in $d_{0}$ (Fig. 6) 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. Sherwin and Blackburn
(2005) and Usmani Usmani and Muralidhar (2016).

The enlargement of the recirculation zone is present through the experiments.
This could be explained in terms of circulation as pointed out in the work of
Gharib et al. Gharib  _et al._ (1998) which studies the flow led by a moving
piston into an unbounded domain. Such work states that vortex formed up to a
limited amount of circulation before it sheds from the layer where is created.
The main difference with our work consist on the confinement that the walls
impose to the vortex. Giving as result that vortex which sheds with a certain
amount of circulation enlarges in the axial direction. This mechanism can also
explain the generation of secondary vortex. Specifically, 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 secondary vortex.

Finally, experimental and numerical results were compared in order to validate
the simulation. Figure 7 shows the evolution of the simulated flow at
$\bar{Re}$ = 654 and $\bar{Re}$ = 1002. Comparison of Fig. 7 with Fig. 5
confirms a similar behaviour for the experimental and simulated flows. The
main structures, i.e. central high-velocity jet and recirculation zone, were
satisfactorily reproduced by numerical simulation. As expected, the flow
pattern had a sharper definition in the simulation. Nevertheless, the relative
error between the numerical and experimental velocity fields was less than 9%,
showing a satisfactory degree of agreement.

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
simulation results for $\bar{Re}$ = 654; (c) streamlines derived from
simulation results for $\bar{Re}$ = 1002. In both cases, red streamlines
represent the central jet layer and blue streamlines represent the central jet
(colours available online only).

###  III.2 Vortex Propagation

Vortex propagation was studied by measuring the 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 to obtain the
filtered vorticity field $\bar{\zeta}(r,z)$. Then, the vortex position was
calculated as $(\bar{r},\bar{z})=\sum(r,z)\bar{\zeta}(r,z)$. The vortex
propagation was along $z$ and is specified in figs. 5, 6 and 7 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 on
the next subsections.

###  III.3 Numerical results for varying Womersley number

The maximum vortex displacement was defined as $z^{*}/D$. The 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. The 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 $\alpha$, or $f$,
was studied numerically for fixed values of $\bar{Re}$ = 654 and $\bar{Re}$ =
1002 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 8 clearly illustrates the dependence of the maximum vortex displacement
on $\alpha$. For a fixed value of $\bar{Re}$ , as $\alpha$ 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 $\alpha$, vortices
tended to have a shorter lifespan and their maximum displacement was therefore
smaller. In other words, $z^{*}/D$ decreased with increasing pulsatile
frequency.

Figure 8: Dependence of the dimensionless maximum vortex displacement on
$\alpha$ (numerical results). Data for $\bar{Re}$ = 654 (red circles) and
$\bar{Re}$ = 1002(blue squares)

###  III.4 Scaling law

A full description of vortex displacement was made on previous sections.
Specifically, it was studied the maximum displacement of the vortex, that is
the distance it travels before it vanishes. Then, we propose a scaling law
which summarizes the behavior of $z^{*}/D$ as function of the involved
parameters, $\bar{Re}$ and $\alpha$.

The physical dimensions for the relevant variables are $\nu$, $D$, $f$ and
$v_{0}$. Based on the Vaschy-Buckingham theorem, it is possible to describe
$z^{*}/D$ as a function of two independent dimensionless numbers, $\bar{Re}$
and $\alpha$ which involve the mentioned relevant variables. The maximum
displacement $z^{*}/D$ was found to be proportional to the flow velocity
$v_{0}$ (i.e. $\bar{Re}$ ) and inversely proportional to the pulsation
frequency $f$ (i.e. $\alpha^{2}$), suggesting that

$z^{*}/D=K\frac{\bar{Re}}{\alpha^{2}}$ (1)

where $K$ is the proportionality constant.

Figure 9: Dimensionless maximum vortex displacement, $z^{*}/D$, as a function
of the dimensionless parameter $\bar{Re}/\alpha^{2}$ for experimental data
(blue squares) and numerical data (red circles).

Figure 9 shows $z^{*}/D$ as a function of $\bar{Re}/\alpha^{2}$ , showing that
both the experimental data and the numerical results provide a good fit to Eq.
1. The fitted slope was $K=1.15\pm 0.19$ at a confidence level of 95%. This
result, remarks the dependence of $z^{*}/D$ with $\bar{Re}$ and $1/\alpha^{2}$
as could be inferred from data of figure 8.

The work of Gharib et al. Gharib  _et al._ (1998) shows that in the developing
process of the vortex, it reaches a threshold of circulation before it sheds.
In case an excess of circulation is generated, this could derive in the
formation of other structures such as secondary vortex. This means that
pulsatile frequency $f$ coincides with the shedding frequency of the primary
vortex, the one that was tracked. Hence, the parameter $\bar{Re}/\alpha^{2}$
can be related to the Strouhal number through
$\bar{Re}/\alpha^{2}=\frac{2}{\pi}\frac{d_{0}}{D}\frac{1}{\mathrm{Sr}}$. A
discussion of the results in terms of Strouhal number clarifies the relation
between the oscillatory component of the flow and the stationary component of
the flow. For lower values of $\bar{Re}/\alpha^{2}$ , that is $\mathrm{Sr}$
close to 1, stationary and oscillatory components are comparable, which
explain vortex with low to null displacement. The extreme case occurs when the
flow oscillates but no vortex are shed. For instance, this could be attained
for Reynolds below the used. As $\bar{Re}/\alpha^{2}$ increases, that is for
$\mathrm{Sr}\lesssim 10^{-1}$, vortex are shed and travel further before it
vanishes as consequence of a stationary component of the flow greater than the
oscillatory component. Parameter $\bar{Re}/\alpha^{2}$ also states the
kinematic nature of vortex displacement, by rewriting
$\frac{\bar{Re}}{\alpha^{2}}=\frac{2d_{0}\bar{v}}{\pi D^{2}f}$ we conclude
that $z^{*}/D$ depends on kinematics 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 $\bar{Re}$ and $\alpha$ on the range studied for this
particular shape of constriction. Specifically for the dependence on the
constriction shape, preliminary results with a Guassian-shaped constriction
were obtained numerically, giving a difference below 20% on the maximun vortex
displacement for $Re_{\mathrm{u}}=$1187\. We encourage further research
regarding the dependence on constriction shape.

##  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 of $\bar{Re}$ and $\alpha$. After a validation
of numerical results, simulations were run over a range of $\alpha$ values
that could not be tested experimentally in order to identify trends in the
flow behaviour with varying $\alpha$.

The flow structure was found to consist of a central 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 $\alpha$, vortex
size grows while decreasing $d_{0}$, and vortex displacement was larger for
increasing $\bar{Re}$ values. The dependence with $\alpha$, stated that as
$\alpha$ increased, the behavior of the flow was reproduced in a shorter
extension of the tube, that means for lower $z/D$ values. The vortex
trajectory was tracked and the vortex displacement over its lifetime
determined. Results showed that the vortex displacement over its lifetime
decreased with increasing $\alpha$. The analysis led to a scaling law
establishing a linear dependence of the vortex displacement on a dimensionless
parameter combining $\bar{Re}$ and $\alpha$, namely $\bar{Re}/\alpha^{2}$ .
Moreover, this parameter was found to be proportional to the inverse of
Strouhal number. This relates directly the vortex behavior with 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, this
results provide insight of the vortex shedding and displacement on a
simplified model of stenosed arteries. This becomes crucial since vortex
shedding precedes turbulence and turbulence is related to plaque complications
finally leading 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.

##  V Acknowledgements

This research was supported by CSIC I+D Uruguay, through the I+D project 2016
“Estudio dinámico de un flujo pulsátil y sus implicaciones hemodinámicas
vasculares”, ANII (Doctoral scholarship POSNAC-2015-1-109843), ECOS-SUD
(project reference number U14S04) and PEDECIBA, Uruguay. The manuscript was
edited by Eduardo Speranza.

## References

  * Naghavi  _et al._ (2003) M. Naghavi _et al._ , “From vulnerable plaque to vulnerable patients. a call for new definitions and risk assesment strategies part 1,” Circulation 108, 1664–1772 (2003).
  * Muller  _et al._ (1989) J.E. Muller, G.H Tofler, and P.H. Stone, “Circadian variation and triggers of onset of acute cardiovascular disease,” Circulation 79, 733–743 (1989).
  * Loree  _et al._ (1991) H.M. Loree, R.D. Kamm, C.M. Atkinson, and R.T. Lee, “Turbulent pressure fluctuations on surface of model vascular stenosis,” American Journal of Physiology 261, 644–650 (1991).
  * Gertz and Roberts (1990) S. D. Gertz and W.C. Roberts, “Hemodynamic shear force in rupture of coronary arterial athersoclerosis plaques,” American Journal of Cardiology 66, 1368–1372 (1990).
  * Barger  _et al._ (1984) A.C. Barger, R. Beeuwkes 3rd, L.L. Lainey, and K. J. Silverman, “Hypothesis: vasa vasorum and neovascularization of human coronary arteries. a possible role in the pathophysiology of atheriosclerosis,” N. Engl. J. Med 19, 175–177 (1984).
  * Falk  _et al._ (1995) E. Falk, P.K. Shah, and V. Fuster, “Coronary plaque disruption,” Circulation 92, 657–371 (1995).
  * Imoto  _et al._ (2005) K. Imoto _et al._ , “Longitudinal structural determinants of atherosclerotic plaque vulnerability: a computational analysis of stress distribution using vessel models and three-dimensional intravascular ulrasound imaging,” Journal of the American College of Cardiology 18, 1507–1515 (2005).
  * Caro  _et al._ (1971) C.G. Caro _et al._ , “Atheroma and arterial wall shear observations, correlation and proposal of a shear dependent mass transfer mechanism for atherogenesis,” Proc. R. Soc. London Ser B 17, 109–159 (1971).
  * Ku  _et al._ (1985) D. N. Ku _et al._ , “Pulsatile flow and atherosclerosis in the human carotid bifurcation,” Arteriosclerosis 5, 293–302 (1985).
  * Nerem and Levesque (1987) R. M. Nerem and M.J. Levesque, “Fluid mechanics in atherosclerosis,” (McGraw-Hill, 1987) Chap. 21, pp. 21.1–21.22.
  * Gopalakrishnan  _et al._ (2014) Shyam Sunder Gopalakrishnan, Benoît Pier, and Arie Biesheuvel, “Dynamics of pulsatile flow through model abdominal aortic aneurysms,” Journal of Fluid Mechanics 758, 150–179 (2014).
  * Arzani  _et al._ (2012) A. Arzani _et al._ , “In vivo validation of numerical prediction for turbulence intensity in aortic coarctation,” Ann Biomed Eng 40(4), 860–870 (2012).
  * Ford  _et al._ (2005) M.D. Ford _et al._ , “Virtual angiography for visualization and validation of computational models of aneurysm hemodynamics,” IEEE Transactions on Medical Imaging 24, 1586–1592 (2005).
  * Boussel  _et al._ (2009) L. Boussel _et al._ , “Phase-contrast magnetic resonance imaging measurements in intracranial aneurysms in vivo of flow patterns, velocity ffields and wall shear stress: comparisson with computational fluid dynamics,” Magn. Reson. Med 61, 409–417 (2009).
  * Long  _et al._ (2001) Q. Long _et al._ , “Numerical investigation of physiologically realistic pulsatile flow through arterial stenosis,” J.Biomechanics 34, 1229–1242 (2001).
  * Stettler and Hussain (1986) J.C. Stettler and A.K.M. Fazle Hussain, “On transition of pulsatile pipe flow,” J. Fluid Mech. 170, 169–197 (1986).
  * Peacock  _et al._ (1998) J. Peacock, T. Jones, and R. Lutz, “The onset of turbulence in physiological pulsatile flow in a straight tube,” Experiments in Fluids 24, 1–9 (1998).
  * Chua  _et al._ (2009) C.S. Chua _et al._ , “Particle image of non axisymmetic stenosis model,” 8th international symposium on particle image velocimetry (2009).
  * Ahmed and Giddens (1983) S.A. Ahmed and D.P. Giddens, “Flow disturbance measurements through a constricted tube at moderate reynolds numbers,” J.Biomechanics 16, 955–963 (1983).
  * Ahmed and Giddens (1984) S.A. Ahmed and D.P. Giddens, “Pulsatile poststenotic flow studies with laser doppler anemometry,” J.Biomechanics 17, 695–705 (1984).
  * Isler  _et al._ (2018) João A. Isler, Rafael S. Gioria, and Bruno S. Carmo, “Bifurcations and convective instabilities of steady flows in a constricted channel,” Journal of Fluid Mechanics 849, 777–804 (2018).
  * Mittal  _et al._ (2003) R. Mittal _et al._ , “Numerical study of pulsatile flow in a constricted channel,” J.Fluid Mech. 485, 337–378 (2003).
  * Sherwin and Blackburn (2005) S.J. Sherwin and H.M. Blackburn, “Three dimensional instabilities and transition of steady and pulsatile axisymmetric stenotic flows,” J. Fluid Mech. 533, 297–327 (2005).
  * Ling and Atabek (1972) S.C. Ling and H.B. Atabek, “A nonlinear analysis of pulsatile flow in arteries,” J. Fluid Mech. 55, 493–511 (1972).
  * Griffith  _et al._ (2009) M.D. Griffith _et al._ , “Pulsatile flow in stenotic geometries: flow behaviour and stability,” J.Fluid Mech. 622, 291–320 (2009).
  * Usmani and Muralidhar (2016) A.Y. Usmani and K. Muralidhar, “Pulsatile flow in a compliant stenosed asymmetric model,” Exp. Fluids 57(12), 186 (2016).
  * Westerweel (1997) J. Westerweel, “Fundamentals of digital particle velocimetry,” Measurement science and technology 8, 1379–1392 (1997).
  * Raffel  _et al._ (2007) M. Raffel, C. Willert, S. Wereley, and J. Kompenhans,  _Particle image velocimetry, a practical guide._ , 2nd ed. (Springer, 2007).
  * Cassanova and Giddens (1977) R. A. Cassanova and D. P. Giddens, “Disorder distal to modeled stenoses in steady and pulsatile flow,” J. Biomechanics 11, 441–453 (1977).
  * Balay (2012) G. Balay,  _Elasticidad en tejidos arteriales, diseño de un corazón artficial in vitro y nuevo método ultrasónico de determinación de elasticidad arterial_ , Master’s thesis, Universidad de la República, Uruguay (2012).
  * Ku (1997) D.N. Ku, “Blood flow in arteries,” Ann. Rev. Fluid Mech 29, 399–434 (1997).
  * Shadden  _et al._ (2005) S.C. Shadden, F. Leiken, and J.E. Marsden, “Definition and properties of lagrangian coherent structures from finite-time lyapunov exponents in two-dimendional aperiodic flows,” Physica D 212, 271–304 (2005).
  * Gharib  _et al._ (1998) M. Gharib, E. Rambod, and K. Shariff, “A universal time scale for vortex ring formation,” Journal of Fluid Mechanics 360, 121–140 (1998).