## Abstract

Many complications in physiology are associated with a deviation in flow in arteries due to a stenosis. The presence of stenosis may transition the flow to weak turbulence. The degree of stenosis as well as its configuration whether symmetric or non-symmetric to the parent artery influences whether the flow would stay laminar or transition to turbulence. Plenty of research efforts focus on investigating the role of varying degrees of stenosis in the onset of turbulence under steady and pulsatile flow conditions. None of the studies, however, have focused on investigating this under oscillatory flow conditions as flow reversal is a major occurrence in a number of physiologic flows, and is of particular relevance in cerebrospinal fluid flow research. Following up on the previous work in which a $75\%$ stenosis was studied, this contribution is a detailed investigation of the role of degrees of stenosis on transition in an oscillatory flow. A cylindrical pipe with $25\%$, $50\%$ and $60\%$ reductions in area in axisymmetric and eccentric configurations is studied for transition with 3 different pulsation frequencies of a purely oscillatory flow. Cycle averaged Reynolds numbers between 1800 and 2100 in steps of 100 are studied for each configuration resulting in 72 simulations each conducted on 76,800 CPU cores of a modern supercomputer. It is found that a higher degree of stenosis and eccentricity causes earlier transition to turbulence in oscillatory flow. The results further demonstrate that a higher frequency of oscillation results in larger hydrodynamic instability in the flow, which is more prominent in smaller degrees of stenosis.

### Supplementary Information

The online version contains supplementary material available at 10.1007/s10237-022-01579-0.

**Keywords:** Lattice Boltzmann method, Transitional flow, Turbulence, Stenosis, Hydrodynamic instability

## Introduction

A stenosis refers to obstruction in a blood vessel or an anatomical segment like the subarachnoid space, spinal canal or respiratory airway. A stenosis can induce turbulence and flow separation in the blood vessel or the anatomical segment as a result of adverse pressure gradients. The presence of turbulence can eventually lead to secondary complications to the physiology in question. The onset of turbulence itself is likely to be influenced by the degree of stenosis itself as well as the nature of the flow, whether steady, pulsatile or oscillatory. For example, stenoses in the carotid arteries are well known risk factor for ischemic stroke(Fairhead and Rothwell 2005). The degree of stenosis has been established as a factor for stroke risk and indication for intervention(Qureshi etal. 2007). A recent review byBrinjikji etal. (2016) delineates the current advances in plaque imaging and association between degree of stenosis and the progression of plaque, whileBerger and Jou (2000) reviewed studies of stenotic flows with a focus on atherosclerosis.

Previous studies have extensively investigated the complexity of fluid dynamics in stenosed vessels both numerically(Varghese etal. 2007a, b; Samuelsson etal. 2015) and experimentally(Ahmed and Giddens 1983, 1984). While studies of steady stenotic flow provide a general guidance about the nature of flow in such vessels, researchers try to assess the parts of cardiac cycles with most turbulent nature of the flow. In particular in the late systole when the flow decelerates, the flow has been shown to be highly disturbed, which stabilizes during the early diastole or acceleration phase of the cardiac cycle. The early studies ofYoganathan etal. (1986), Bluestein and Einav (1995) investigated in detail as to what degree of stenosis causes the flow to become turbulent.

Most of the studies of stenotic flows usually focus on uni-directional flows. This is natural because the areas where a stenosis are most relevant have a unidirectional steady or pulsatile flow. Earlier experiments ofAhmed and Giddens (1983) investigated post stenotic flow in a tube using LDA at mean Reynolds number of 600 and Womersley number ${\text{W}}_{\text{o}}=7.5$. They found strong turbulence and flow separation downstream of the $75\%$ stenosis. Ding etal. (2021) recently studied pulsatile flow in a 2D stenosed pipe using PIV experiments with a mean Reynolds number of 1750 and Womersley number of ${\text{W}}_{\text{o}}=6.15$. They employed $25\%$, $50\%$ and $75\%$ degree of stenosis and found increasing turbulent intensities with increasing constrictions. Xu etal. (2017) carried out experiments in a sinusoidal shaped pulsatile pipe flow and observed that the critical Reynolds number increased with an increase in the Womersley number until ${\text{W}}_{\text{o}}=2.5$. Then the critical Reynolds number reduced with a further increase in the Womersley number until it reached to ${\text{W}}_{\text{o}}=12$.

Some parts in the circulation however do experience a flow reversal. For example, the waveform of blood flow through arteries has been shown to change abruptly when passing from the thoracic aorta into the abdominal aorta in humans, exhibiting regions of flow reversal at the end of systole. Holenstein and Ku (1988) developed a lumped parameter model of such a flow and referred to it as *triphasic* flow with positive-negative-positive parts of the cardiac cycle. Another form of such a triphasic flow is a purely oscillatory flow, known to occur mainly in the central nervous system (CNS) namely the circulation of the cerebrospinal fluid (CSF). Another purely oscillatory flow is the airflow in the respiration system. In a purely oscillatory flow the underlying factors that trigger turbulence, the distribution of turbulence in various parts of the cycle, as well as the stabilization mechanisms are different than those from a steady or pulsatile flow with no reversal component. Earlier experiments byHino etal. (1976, 1983) have investigated detailed characteristics of turbulence in oscillatory pipe flow. Other similar studies(Sarpkaya 1966; Yellin 1966; Iguchi and Ohmi 1982) addressed such aspects but also in simple pipes without any distortion that would resemble an anatomical geometry.

In a recent workJain (2020b) we investigated transition to turbulence in a purely oscillatory flow with three pulsation frequencies in a pipe with $75\%$ area reduction in eccentric and axisymmetric configurations. While the role of symmetry of the geometry and the pulsation frequencies in the onset of turbulence was the focus of that study, the degree of stenosis itself was not investigated. The present contribution investigates the characteristics of turbulence in a stenosed pipe with area reductions of $60\%$, $50\%$ and $25\%$. All three degrees of stenosed pipes are studied in two configurations: one that is symmetric to the principal axis of the pipe while the other which is offset by 0.05 diameters of the pipe to introduce an eccentricity or a geometric aberration. Three different flow oscillation frequencies with Womersley numbers ranging between $\sim $ 5 and 8 are studied in all the degree and configurations of the stenosis to explore detailed characteristics of flow transition. Reynolds numbers between 1800 and 2100 are investigated in steps of 100 to compare the characteristics of flow transition with varying degrees of stenosis against the previous study(Jain 2020b).

The three main focal points of the study namely *degree of stenosis*, *pulsation frequencies* and *symmetry/asymmetry of the geometry* are representative of a broad range of CSF and respiratory flows. For example stenosis of the lumbar spine ranges from various degrees within the range studied here(Genevay and Atlas 2010) while the obstruction to CSF caused by Chiari I malformation also lies in this range(Linninger etal. 2016). Subglottic stenosis is a narrowing of a specific portion of the windpipe (trachea) known as the subglottis (just below the vocal cords). These are classified as grade 1 stenosis with a luminal obstruction $<50\%$, grade 2 and 3 stenosis with 51–70% and 71–99% obstruction respectively(Jefferson etal. 2016). The pulsation frequencies range from a Womersley number of approximately 5 (CSF) to 11 (respiratory airflow), which corresponds to the pulsation frequencies studied in this work. Occurrence of stenosis in a symmetric fashion is rare due to the complexity of the anatomy and most obstructions are found in a highly asymmetric fashion. Similarly, the Reynolds numbers, mainly due to the presence of stenosis can reach up to 2000 in CSF and 1500–8000 in respiratory airways. The present study is thus a generalization of stenotic flow in these two applications. The idea is to introduce representative anatomic and physiologic complexity while maintaining simplicity that could allow prospects for future studies in identifying the chances of turbulence. The goal for these canonical studies is to serve as benchmark for future comparisions and thus a patient specific anatomic case has not been studied.

For the absence of a clear definition for a flow regime that is neither laminar nor turbulent, we have addressed the flow in this study as transitional. Thus, if fluctuations occur in localized parts of the domain during parts of the cycle, we term the flow transitional. This is in accordance with several such studies in this direction(Varghese etal. 2007b; Samuelsson etal. 2015; Yellin 1966; Jain 2020a, b). The onset of a fully developed turbulence is not pursued and the focus is on exploring the critical Re at which the flow leaves a laminar regime for a particular degree of stenosis, and a specific pulsation frequency (Womersley number). The Reynolds number at which the flow transitions is judged on the basis of hydrodynamic instabilities in the velocity field and deflection of the jet downstream of the stenosis.

We have employed the lattice Boltzmann method (LBM) for the DNS of this study as the method scales well on parallel computers—allowing for accurate capture of transitional flow characteristics. The LBM when setup properly(Jain 2020a) is a very appropriate second order accurate method for the DNS of flow. While the order of accuracy is lower than competitive higher order methods, the strict control of numerical viscosity in LBM allows for a relatively easier capture of flow fluctuations. A comparison of LBM with other second-order accurate methods(Marié etal. 2009) has found that the numerical dissipation in LBM, even at the scales of grid spacing and the numerical dispersive effects are relatively smaller.

## Methods

The geometry of the stenosis was adopted from the experiments ofAhmed and Giddens (1984) (shown in Fig.1). Stenosis of varying degrees in axisymmetric and eccentric configurations were created using the Blender software previously byGericke (2016). The gray, green and red lines respectively demonstrate $60\%$, $50\%$ and $25\%$ reduction in area whereas there dotted counterparts indicate the eccentricity that was introduced.

Similar to previous studyJain (2020b) the pre and post-stenotic regions of the vessel were *equally* extended by 12 vessel diameters to account for the reversing nature of the flow, which resulted in a total pipe length of approximately 26 diameters. Due to flow reversal, the fluid travels equally on both sides of the stenosis throat, and the location of stenosis in the pipe is thus expected to influence the onset of transition as well as propagation of formed turbulent bursts.

### Oscillatory flow

The oscillatory nature of the flow can be described by the Womersley solution with a negative (flow reversal) component. The Womersley solution for laminar, pulsatile flow through rigid tubes is defined as:

$$\left(\begin{array}{c}\frac{{u}_{x}}{{u}_{\text{c}}}=A\left[1-\frac{{J}_{0}({i}^{3/2}\alpha 2r/D)}{{J}_{0}({i}^{3/2}\alpha )}\right]sin(\omega t),\hfill \\ \frac{{u}_{y}}{{u}_{\text{c}}}=0\hfill \\ \frac{{u}_{z}}{{u}_{\text{c}}}=0\hfill \end{array}\right\}$$ | 1 |

where ${u}_{\mathrm{c}}$ is the half-cycle-averaged inlet centerline velocity,* A* and $\omega $ respectively are the amplitude and angular frequency of pulsation, ${J}_{0}$ is the Bessel function of first kind and zeroth order, and $\alpha $ is the non-dimensional Womersley parameter ($=\frac{1}{2}D\sqrt{\omega /\nu}$, where $\nu $ is the kinematic viscosity). The Womersley parameter defines the extent to which the laminar profile departs from quasi-steadiness. This effect becomes significant when the Womersley parameter, $\alpha =3$. Three different frequencies of pulsation were studied where ${\omega}_{2}=2{\omega}_{1}$ and ${\omega}_{3}=4{\omega}_{1}$. The forward and backward flow are referred to as *blowing* and *suction* stages respectively in the following.

Figure2 depicts the flow waveform and the 6 time intervals of* T*/6 where flow quantities were analyzed. Like the studies ofJain (2020b), the flow profile here consists of a negative component which accounts for the flow reversal, and the mean flow in this case is thus zero describing a purely oscillatory flow. The Reynolds number was based on the main vessel diameter, *D*, and the mean inlet centerline velocity of *half the cycle*, ${u}_{\mathrm{c}}$. The Re was varied from 600 to 2100 in steps of $\delta \text{Re}=100$ for both the configurations of the stenosis and three pulsation frequencies. This resulted in a total 90 sets of simulations.

The value of A and $\alpha $ in Eq.1 were 0.667 and 5.4 respectively. For the lowest pulsation frequency ${\omega}_{1}$ the Womersley parameter ${\alpha}_{1}$ was chosen as 5.4 keeping the frequency of oscillations of the cerebrospinal fluid (CSF) in view. To attain the pulsation frequencies of ${\omega}_{2}=2{\omega}_{1}$ and ${\omega}_{3}=4{\omega}_{1}$, the corresponding Womersley parameters became ${\alpha}_{2}=7.636$ and ${\alpha}_{3}=10.8$ respectively. Zero pressure and velocity were set at every lattice cell as initial conditions. A zero pressure gradient was maintained at the outlet of the pipe which translates to an extrapolation of incoming populations at the outlet(Junk and Yang 2011). The flow waveform shown in Fig.2 was prescribed at the inlet. Flow was allowed to develop for two initial cycles before analysis was performed.

### Direct numerical simulation

The adaptable poly engineering simulator (APES) framework(Roller etal. 2012; Klimach etal. 2014; Qi etal. 2016) was chosen for the simulations which contains a full tool-chain of simulation software ranging from the mesh generator *Seeder*(Harlacher etal. 2012) and the LBM solver *Musubi*(Hasert etal. 2014). Modules within APES are available as an open source software tool for academic research.1 The *Musubi* LBM solver is mainly managed at the Institute of Software Methods for Product Virtualization within the German Aerospace Center (DLR) in Dresden, Germany. Extensive development and use of the solver for applications in physiology is done at the University of Twente, The Netherlands.

Extensive verification of the LBM solver *Musubi* has been performed in prevous works(Hasert etal. 2014). About 30 verification test cases are executed automatically every night that ensure that the physics computation and parallel as well as serial performance of the solver compares well with the benchmarks. Several efforts in the direction of validation have been placed, see for example a comparison against experiments reported inJohannink etal. (2015). A comprehensive validation for the simulation of physiologic flows in transitional and turbulent regime was conducted using the benchmark set out by the US Food and Drug Administration (FDA)(Jain 2020a). For each new biomedical application, comparisons against in vitro and in vivo measurements are performed during the continuous use and development of *Musubi*. See for example, a recent work conducted in parallel in which simulation results were compared against *in vitro* experiments of nasal airflows(Hebbinketal. 2022).

Due to the simplicity of the geometry studied in this work, the employment of the Bhatnagar–Gross–Krook (BGK) relaxation scheme was sufficient. The BGK relaxation parameter was set to $\mathrm{\Omega}=1.84$ in the present study that keeps the lattice Mach number within the stability limits of the LBM(Junk etal. 2005). Further details on these aspects can be referred elsewhere(Jain 2020a).

The spatial and temporal resolutions were chosen as $\delta x=32\times {10}^{-4}$ and $\delta t=7.5\times {10}^{-6}$ units respectively. Based in the degree of stenosis, this resulted in $\sim $ 700 million cells. When the pipe had the minimum degree of stenosis ($25\%$), there was more volume of the fluid resulting in approximately 1 million more lattice cells than the most severe case. Presuming the main vessel diameter, *D*, of 1 unit, the spatial resolution $\delta x$ resulted in 124, 156 and 232 cells along the diameters of the throat for $60\%$, $50\%$ and $25\%$ degree stenoses respectively. The number of cells along the main diameter of the parent pipe was 312. This ensured the symmetry of the mesh. It may be remarked that the number of lattice cells here does not translate directly to the *order* of polynomial or *resolution in space* for other higher order methods like the spectral methods. We employed the D3Q19 stencil of the LBM in this study, which means that each lattice cell had 19 degrees of freedom. In higher order methods, the order of the polynomial allows for a higher accuracy within a cell, which is not the case with LBM. Thus these resolutions cannot be compared with those employed in similar studies(Varghese etal. 2007b). A detailed comparison of methods is left for future efforts.

Our previous worksKlimach etal. (2014), Hasert etal. (2014), Qi etal. (2016) have shown that the *Musubi* LBM solver exhibits excellent weak and strong scaling as well as parallel efficiency of more than $90\%$ on supercomputers when every core has more than about 4000 lattice nodes. Simulations were thus executed on 76,800 cores of the *SuperMUC-NG* supercomputer installed at the Leibniz Supercomputing Center in Munich, Germany. With this setup each core (or MPI rank) had about 9000 lattice nodes for computations thereby reducing the communication to computation ratio to an optimal minimum, and ensuring fully efficient utilization of compute resources.

Simulations at Reynolds of 2000 and 2100 are only reported in this work due to the laminar nature of the flow at lower Reynolds numbers. The simulations for all the cases were first conducted for only $n=4$ cycles as in case of laminar flow, computation of more cycles was redundant. Some of these simulations were conducted using 7200 CPUs of the Snellius system, which is the national supercomputer in The Netherlands. In case of a transitional or a turbulent flow, it is always challenging to decide how many total cycles should be used for ensemble averaging until the quantities are statistically converged as the compute costs and physical details need to be balanced. Ideally, the ensemble average should attain the shape of a laminar flow, similar to inflow or the turbulent fluctuations should wash out from the averaged plot.

Due to the lack of this information a priori, we first simulated the flow with pulsation frequency ${\omega}_{1}$. In the cases where flow transition was identified, we restarted those simulations from the 4th cycle onwards and allowed the computation of further $n=20$ cycles resulting in a total of $n=24$ cycles. For higher pulsation frequencies, due to the same time step, the number of cycles automatically became $n=48$ for ${\omega}_{2}$ and $n=96$ for ${\omega}_{3}$. Each cycle required $\sim $ 4.3 min of execution time. An *abort criteria* was set in *Musubi* to stop the simulation upon achievement of a steady ensemble average. Thus, each subsequent cycle was ensemble averaged with the previously computed cycles and when the velocity fluctuations in ensemble averages vanished, the simulations were stopped automatically. This analysis revealed that a total of about $n=20$ cycles were enough for the most turbulent simulation, which further implicated the choice of $n=20$ cycles for all the cases. Based on this a total of $n=22$ cycles were computed for each simulation and the last 20 were used for the analysis of turbulent flow characteristics. It may be noted that 20 cycles are relatively less than what would be needed for statistical convergence. A recent study by Andersson and Karlsson (2021) reports that employment of much larger number of cycles is required for statistical convergence. The 20 cycles in this case are sufficient because of the zero mean nature of the flow, which is absent in a pulsatile blood flow reported in the study ofAndersson and Karlsson (2021).

### Flow characterization

A total of $n=22$ (where initial 2 cycles are discarded from analysis) cycles were computed for each Re simulation. Velocity over these cycles was ensemble averaged for the analysis of turbulent characteristics as:

$$\begin{array}{c}\hfill \overline{u}(x,t)=\frac{1}{n}\sum _{k=0}^{n-1}u(x,t+kT)\end{array}$$ | 2 |

where *u*(*x*,*t*) is the instantaneous velocity at a particular lattice site, **x** denotes the spatial coordinates, **t** is the time and **T** is the period of cycle.

We assumed that turbulence was the only source of fluctuations in the flow and composed the velocity field into a mean and a fluctuating component. This is known as Reynolds’ decomposition and is commonly used for steady flows. The oscillatory flows have a reversal component in addition to the pulsating nature of the flow. Triple decomposition(Hussain and Reynolds 1970) is thus considered a more suitable method for the analysis of such flows:

$$\begin{array}{c}\hfill {u}_{i}(x,t)=\overline{{u}_{i}}(x)+{u}_{i}^{\prime}(x,t)+\stackrel{~}{{u}_{i}}(x,\varphi )\end{array}$$ | 3 |

Thus, a periodic component denoted by $\stackrel{~}{{u}_{i}}(x,\varphi )$ gets added to Reynolds’ decomposition, which is a function of the time within the cycle. The phase is defined as $\varphi =mod(t/T,1)$ where* T* is the total length of the cycle. Cycle-to-cycle variations are also intrinsically taken into account by the method of triple decomposition.

The Turbulent Kinetic Energy (TKE) is derived from the fluctuating components of velocity in 3 directions as:

$$\begin{array}{c}\hfill k=\frac{1}{2}({u}_{x}^{\prime 2}+{u}_{y}^{\prime 2}+{u}_{z}^{\prime 2})\end{array}$$ | 4 |

The frequency components present in a transitional or turbulent flow can be analyzed from the power spectral density (PSD) of the TKE. The PSD was computed in this work using the Welch’s periodogram method, which was related with Kolmogorov’s energy decay to observe the inertial and viscous ranges.

### Kolmogorov microscales

In the previousJain (2020b) as well as in related workHelgeland etal. (2014) the Kolmogorov theory has been taken as a reference to estimate the quality of the mesh employed for DNS. It may be noted that the degrees of stenosis considered in this article are relatively lower than those in the previous work where stenosis with a $75\%$ area reduction was chosen. This implies that the amount of turbulence (if any) should be lower in the present work, which in turn suggests that the resolution employed in simulations should be enough for the accurate capture of transitional characteristics. It was however not a priori known whether turbulence will occur in these cases at Reynolds of up to 2100 or not. This means that simulations at even higher Reynolds number could have been required and thus Kolmogorov microscales were still calculated,

The Kolmogorov scales are defined as spatial and temporal scales when viscosity dominates and the TKE is dissipated into heat(Pope 2000).

The Kolmogorov scales, non-dimensionalized with respect to the velocity scale ${u}_{\mathrm{c}}$ and the length scale *D* are computed from the *fluctuating* component of the non-dimensional strain rate defined as:

$$\begin{array}{c}\hfill {s}_{\mathit{ij}}^{\prime}=\frac{1}{2}(\frac{\partial {u}_{i}^{\prime}}{\partial {x}_{j}}+\frac{\partial {u}_{j}^{\prime}}{\partial {x}_{i}})\end{array}$$ | 5 |

The Kolmogorov length, time and velocity scales are then respectively computed as:

$$\begin{array}{cc}\hfill \eta & =(\frac{1}{{\text{Re}}^{2}}\frac{1}{2{s}_{\mathit{ij}}^{\prime}{s}_{\mathit{ij}}^{\prime}}{)}^{1/4}\hfill \end{array}$$ | 6 |

$$\begin{array}{cc}\hfill {\tau}_{\eta}& =(\frac{1}{2{s}_{\mathit{ij}}^{\prime}{s}_{\mathit{ij}}^{\prime}}{)}^{1/2}\hfill \end{array}$$ | 7 |

$$\begin{array}{cc}\hfill {u}_{\eta}& =(\frac{2{s}_{\mathit{ij}}^{\prime}{s}_{\mathit{ij}}^{\prime}}{{\text{Re}}^{2}}{)}^{1/4}\hfill \end{array}$$ | 8 |

Based on these scales, the quality of the spatial and temporal resolution of a simulation is estimated by computing the ratio of $\delta x$ and $\delta t$ against the corresponding Kolmogorov scales i.e.

$$\begin{array}{c}\hfill {l}^{+}=\frac{\delta x}{\eta}\phantom{\rule{0ex}{0ex}}{t}^{+}=\frac{\delta t}{{\tau}_{\eta}}\end{array}$$ | 9 |

The ${l}^{+}$ and ${t}^{+}$ were computed using Eq.9 for only the $60\%$ stenosis case at the highest frequency as the fluctuations were highest in this case. The values of these scales were 1.18 and 0.62 respectively.2

## Results

The flow at $\text{Re}=1800$ was laminar and it started to show valleys and peaks in the centerline velocities at $\text{Re}=1900$ while at $\text{Re}=2000$ these valleys and peaks became more pronounced and minor fluctuations started to appear mostly in the $60\%$ case. The TKE components up to this Reynolds number were depictive of a laminar flow with minor vortex shedding in the $60\%$ eccentric stenosis case. The flow only transitioned to weak turbulence at $\text{Re}=2100$ commensurate with the previous study(Jain 2020b) of $75\%$ stenosis. Figures3 and4 shows traces of instantaneous centerline velocity along the streamwise direction for axisymmetric and eccentric cases at $\text{Re}=2000$ and $\text{Re}=2100$ respectively.

Clearly the fluctuations at $\text{Re}=2100$ are markedly higher than those at $\text{Re}=2000$ for the $60\%$ and $50\%$ degrees of stenosis. The flow in both stenosis configurations is essentially laminar for the $25\%$ stenosis. At a higher oscillation frequency of ${\omega}_{3}$, minor distortions in the sinusoidal profile however start to appear even in the $25\%$ case. For the $50\%$ stenosis case the flow shows minor fluctuations during the peak flow and the deceleration phase. In the eccentric case these fluctuations are relatively higher at this degree of stenosis. The flow attains a transitional character only at a stenosis degree of $60\%$ for both the stenosis configurations. The fluctuations are relatively higher in the eccentric case. The regions downstream of the stenosis where the flow transitions and thereafter stabilizes appears to shift closer to the stenosis throat with increasing oscillation frequency. For the eccentric stenosis, the regions of maximum chaos are mostly confined in the locations $2D<x<4D$, $1D<x<3D$, $0D<x<2D$ respectively for oscillation frequencies ${\omega}_{1}$, ${\omega}_{2}$ and ${\omega}_{3}$. In the axisymmetric case, however, the regions of jet breakdown are shifted further downstream by one diameter whereas the regions of restabilization are the same as in eccentric stenosis for all the studied frequencies. At the highest studied frequency, ${\omega}_{3}$, there are signs of fluctuations even at the stenosis throat ($x=0D$), which are a result of fast oscillation of the fluid that transfers the vortices on both sides of the stenosis rapidly.

These findings are further elaborated in the power spectral density (PSD) plots of Fig.5. Each column from top to bottom shows the PSD of the TKE at the stenosis throat and 4 locations downstream of it during the blowing stage whereas the black, green and red lines respectively demonstrate the PSD for $60\%$, $50\%$ and $25\%$ degrees of stenoses. In the axisymmetric case for oscillation frequency ${\omega}_{1}$, the spectrum at $x=2D$ and $x=3D$ begins to take a shape with few frequencies in the inertial range depicted by the $-5/3$ decay. At frequencies higher than ${10}^{3}$ the spectrum rolls off to $-10/3$ and subsequently to $-7$. The spectrum at $x=4D$ is similar for $60\%$ and $50\%$ degree of stenosis while the flow in $25\%$ case is laminar. Large valleys and peaks (red line) are seen in $50\%$ and $25\%$ stenoses at higher oscillation frequencies of ${\omega}_{1}$ and ${\omega}_{2}$ depictive of the minor disturbances that were seen in instantaneous traces as well. Similar observations can be drawn from the PSD plots of the eccentric case with a clear indication of larger turbulent activity. There are larger number of frequencies in the $-5/3$ range and the roll of to $-10/3$ and $-7$ happens later in case of eccentric stenosis.

### Vortex structures

Figures6 and7 show the profiles of vorticity magnitude across a bisecting plane in the axisymmetric and the *x*–*z* bisecting plane in the eccentric stenoses cases of all variants and oscillation frequencies. The observation point P3 (see Fig.2) has been chosen as the flow exhibits maximum chaos at this observation point in the sinusoidally oscillating flow. The suction stages (points P4–P6) have been omitted because more or less similar profiles are seen on the left side of the stenosis throat.3 The previous observation that the regions of jet breakdown and restabilization shift closer to the stenosis throat with increasing oscillation frequencies are clearly seen in the vorticity plots.

In the axisymmetric case the jet that emanates from the throat has enough time to travel downstream at oscillation frequency ${\omega}_{1}$ before its breakdown and transition of flow to weak turbulence. At higher frequencies, the flow breaks down closer to the throat as the deceleration phase commences sooner in these cases. It may be specifically seen that at $60\%$ stenosis and oscillation frequency ${\omega}_{3}$ (Fig.6c) the jet is deflected towards the lower part downstream of the stenosis in this particular cycle. This is a consequence of the fact that when the flow reverses its direction, the vortices merge just downstream of the throat. The frequency (${\omega}_{3}$) is relatively high that breaks the symmetry that is otherwise present at lower frequencies. The momentum of the flow is suddenly lost when it reverses its direction and this happens at a rapid speed at this frequency resulting in this loss of symmetry even in the axisymmetric case. This is also seen to a certain extent in $50\%$ stenosis with ${\omega}_{3}$ (Fig.6f) albeit at a lower extent due to the lower intensity of turbulence in the first place due to less constriction in the pipe. For the $25\%$ stenosis case this phenomenon is never observed as the flow is clearly laminar for all the oscillation frequencies.

The eccentricity of the stenosis causes the jet to deflect earlier than the axisymmetric case resulting in a jet breakdown location closer to the stenosis throat. Here a higher oscillation frequency amplifies the chaos and factors that influence the onset of transition, in addition to eccentricity, are faster blowing and suction as was stated above. For the $50\%$ stenosis case at oscillation frequency ${\omega}_{1}$ (Fig.7f), the jet does not appear to be deflected as a result of eccentricity but an inflexion point is created about 2.5 diameters downstream of stenosis where the turbulent activity occurs and extends beyond to up to about 4 diameters. At higher frequencies in this degree of stenosis, interesting patterns are observed with a similar fluid dynamical behavior as in axisymmetric stenosis.

## Discussion

The direct numerical simulations of oscillatory flow in varying degrees of stenosis in axisymmetric and eccentric configurations have found a correlation between severity of the stenosis and the intensity of turbulence in an oscillatory flow. This influence is a well established fact but none of the previous studies have characterized this for oscillatory flow.

In comparison to the previous study with $75\%$ stenosis(Jain 2020b), this study has contrastingly found that the increasing frequencies of oscillation result in earlier breakdown of the flow, a phenomenon most pronounced in lower degrees of stenosis. In a constriction as large as $75\%$ the onset of turbulence was predominated by the geometry rather than the complexity of the flow field. Upon reducing the constriction in the present study the frequencies of the oscillatory flow field resulted in the onset and sustainment of turbulent flow. At high frequencies of oscillation, the instant at which the downstream flow reversed direction had a lag in different radial directions resulting in merger of vortices and creation of newer ones.

The axisymmetric or the eccentric nature of the geometry here also had less pronounced effect in comparison to the previous study, mainly because of the aforementioned factors. The most interesting distinction between eccentric and axisymmetric cases is seen in $25\%$ degree of stenosis and pulsation frequency ${\omega}_{3}$. As can be seen from Fig.7h, i miniature patches of vortices are created downstream of the stenosis. These vortices that are shed are a consequence of jet deflection, which sows the seeds of flow transition already at this small constriction in the pipe. These seeds are not intense enough to cause a breakdown of the flow at this Reynolds number, and the flow reversal further adds to their incapability in bringing about flow transition.

The intermittent nature of turbulence along the flow cycle demonstrates how flow reversal stabilizes the flow field—commensurate with the previous study (Jain 2020b). The centerline velocities however demonstrate different characteristics with different oscillation frequencies. At larger pulsation frequencies ${\omega}_{2}$ and ${\omega}_{3}$ the turbulence has more a decaying nature because the deceleration phase has a reduced time with increasing frequencies.

One of the most surprising observations stems from a comparison of Figs. 6c and 7c. While the former figure is for axisymmetric case and the latter for eccentric case, the shape of the flow jet seems to suggest as if these figures have been interchanged. The time lapse between events like flow jet breakdown, reversal and stabilization is extremely small in these cases due to the high frequency of the sinusoidal waveform. Due to that the vortices that were created during the deceleration phase are immediately disturbed by the flow reversal that tries to stabilize the flow field and push the flow in the other direction. A fast merger and annihilation of vortices breaks the shape of the post stenotic jet giving rise to this behavior. If we take a closer look at Fig.7c we can observe that the jet is actually deflected upwards due to the eccentricity.4

From the PSD plots (Fig.5) it may be seen that whatever turbulent nature of the flow there is, it is less intense than that for $75\%$ stenosis. At higher degrees of stenoses the PSD attains a $\frac{-5}{3}$ decay albeit for a very small range of frequencies. A roll of to $\frac{-10}{3}$ and $-7$ happens much sooner demonstrating the decaying nature of turbulence.

The degree of constriction, in almost all pathologies indicates the severity of the situation. Bluestein and Einav (1995) focused their study with background of heart valves and for pulsating flow their general observation, that a higher constriction results in larger turbulence, was consistent with the results of the present study. It may be queried how relevant a zero mean oscillatory flow is in physiology. Holenstein and Ku (1988) has demonstrated the relevance of reverse flow in infrarenal vessels as well albeit the flow there only has a reversed component and the mean of the flow is not zero. This study will thus have a main relevance in CSF and respiratory flow research that are purely oscillatory.

This study has several limitations, most of which are similar to the previous study(Jain 2020b). In the context of degrees of stenosis it may be mentioned that despite the study of three different pulsation frequencies, the waveform has always been sinusoidal. In both CSF and respiratory flows, the waveforms are very complex due to a number of physiologic factors. Such waveforms might result in alternate turbulence mechanisms. Also, despite eccentricity the stenosis shape is still symmetric, which is unlikely to happen in an anatomic case. This study may thus be inferred as one where complexity of the anatomy has been incorporated in the form of major factors that are present, while maintaining the constitutive simplicity so that these results can act as a benchmark.

## Supplementary Information

Below is the link to the electronic supplementary material.

Download video file^{ (14MB, mp4) }

Supplementary material 1: Vorticity in the axisymmetric case with 60% stenosis and pulsation frequency w3 (MP4 14381 kb)

Download video file^{ (12.8MB, mp4) }

Supplementary material 2: Vorticity in the axisymmetric case with 60% stenosis and pulsation frequency w2 (MP4 13101 kb)

Download video file^{ (11.7MB, mp4) }

Supplementary material 3: Vorticity in the axisymmetric case with 60% stenosis and pulsation frequency w1 (MP4 11951 kb)

Download video file^{ (9.5MB, mp4) }

Supplementary material 4: Vorticity across the x-z plane in the eccentric case with 60% stenosis and pulsation frequency w1 (MP4 9778 kb)

Download video file^{ (14.2MB, mp4) }

Supplementary material 5: Vorticity across the x-z plane in the eccentric case with 60% stenosis and pulsation frequency w3 (MP4 14491 kb)

Download video file^{ (13MB, mp4) }

Supplementary material 6: Vorticity across the x-z plane in the eccentric case with 60% stenosis and pulsation frequency w2 (MP4 13308 kb)

## Acknowledgements

Compute resources on *SuperMUC* and *SuperMUC-NG* were provided by the Leibniz Supercomputing Center (LRZ), Munich, Germany. The author is grateful for the able support of colleagues at LRZ. Some simulations were performed on the Dutch national supercomputer Snellius. Compute resources on Snellius were provided by the SURFSara through NWO grant 2019/NWO/00768083.

## Declarations

### Conflict of interest

The authors declare that they have no conflict of interest. No funding was received for this specific work.

## Footnotes

^{1}

^{2}

A detailed analysis of Kolmogorov microscales has been omitted because the turbulence intensity in these cases is lower than the previous study(Jain 2020b). Please see the results section for more details.

^{3}

Please refer to the animations in supplementary material for more details.

^{4}

Such observations are more clear in the animations provided as supplementary information.

**Publisher's note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Change history**

4/23/2022

The missing Electronic Supplementary Material is included

## References

- Ahmed SA, Giddens DP. Velocity measurements in steady flow through axisymmetric Stenoses at moderate Reynolds numbers. J Biomech. 1983;16(7):505–516. doi: 10.1016/0021-9290(83)90065-9. [DOI] [PubMed] [Google Scholar]
- Ahmed SA, Giddens DP. Pulsatile poststenotic flow studies with laser doppler anemometry. J Biomech. 1984;17(9):695–705. doi: 10.1016/0021-9290(84)90123-4. [DOI] [PubMed] [Google Scholar]
- Andersson M, Karlsson M. Model verification and error sensitivity of turbulence-related tensor characteristics in pulsatile blood flow simulations. Fluids. 2021;6(1):11. doi: 10.3390/fluids6010011. [DOI] [Google Scholar]
- Berger S, Jou L-D. Flows in stenotic vessels. Ann Rev Fluid Mech. 2000;32(1):347–382. doi: 10.1146/annurev.fluid.32.1.347. [DOI] [Google Scholar]
- Bluestein D, Einav S. The effect of varying degrees of stenosis on the characteristics of turbulent pulsatile flow through heart valves. J Biomech. 1995;28(8):915–924. doi: 10.1016/0021-9290(94)00154-V. [DOI] [PubMed] [Google Scholar]
- Brinjikji W, Huston J, Rabinstein AA, Kim G-M, Lerman A, Lanzino G. Contemporary carotid imaging: from degree of stenosis to plaque vulnerability. J Neurosurg. 2016;12(4):127–142. doi: 10.3171/2015.1.JNS142452. [DOI] [PubMed] [Google Scholar]
- Ding G, Choi K-S, Ma B, Kato T, Yuan W. Transitional pulsatile flows with stenosis in a two-dimensional channel. Phys Fluids. 2021;33(3):034115. doi: 10.1063/5.0042753. [DOI] [Google Scholar]
- Fairhead J, Rothwell PM. The need for urgency in identification and treatment of symptomatic carotid stenosis is already established. Cerebrovasc Dis. 2005;19(6):355–358. doi: 10.1159/000085201. [DOI] [PubMed] [Google Scholar]
- Genevay S, Atlas SJ. Lumbar spinal stenosis. Best Pract Res Clin Rheumatol. 2010;24(2):253–265. doi: 10.1016/j.berh.2009.11.001. [DOI] [PMC free article] [PubMed] [Google Scholar]
- Gericke J (2016) Direct numerical simulation of transition to turbulence in pipe flows with varying degree of stenosis using the lattice Boltzmann method (Bachelor’s Thesis). Universität Siegen, Germany
- Harlacher DF, Hasert M, Klimach H, Zimny S, Roller S (2012) Tree based voxelization of STL Data. In: High performance computing on vector systems, 2011, pp 81–92
- Hasert M, Masilamani K, Zimny S, Klimach H, Qi J, Bernsdorf J, Roller S. Complex fluid simulations with the parallel tree-based lattice Boltzmann solver Musubi. J Comput Sci. 2014;55:784–794. doi: 10.1016/j.jocs.2013.11.001. [DOI] [Google Scholar]
- Hebbink R, Wessels B, Hagmeijer R, Jain K (2022) Simulations of unassisted breathing in a nasal cavity. Med Biol Eng Comput (submitted)
- Helgeland A, Mardal K-A, Haughton V, Anders PetterssonReif B (2014) Numerical simulations of the pulsating flow of cerebrospinal fluid flow in the cervical spinal canal of a Chiari patient. J Biomech [DOI] [PubMed]
- Hino M, Kashiwayanagi M, Nakayama A, Hara T. Experiments on the turbulence statistics and the structure of a reciprocating oscillatory flow. J Fluid Mech. 1983;13(1):363–400. doi: 10.1017/S0022112083001378. [DOI] [Google Scholar]
- Hino M, Sawamoto M, Takasu S. Experiments on transition to turbulence in an oscillatory pipe flow. J Fluid Mech. 1976;75(02):193–207. doi: 10.1017/S0022112076000177. [DOI] [Google Scholar]
- Holenstein R, Ku DN. Reverse flow in the major infrarenal vessels—a capacitive phenomenon. Biorheology. 1988;25(6):835–842. doi: 10.3233/BIR-1988-25604. [DOI] [PubMed] [Google Scholar]
- Hussain AKMF, Reynolds WC. The mechanics of an organized wave in turbulent shear flow. J Fluid Mech. 1970;41(02):241–258. doi: 10.1017/S0022112070000605. [DOI] [Google Scholar]
- Iguchi M, Ohmi M. Transition to turbulence in a pulsatile pipe flow: part 2, characteristics of reversing flow accompanied by relaminarization. Bull JSME. 1982;252(08):1529–1536. doi: 10.1299/jsme1958.25.1529. [DOI] [Google Scholar]
- Jain K. Efficacy of the FDA nozzle benchmark and the lattice Boltzmann method for the analysis of biomedical flows in transitional regime. Med Biol Eng Comput. 2020;58:1817–1830. doi: 10.1007/s11517-020-02188-8. [DOI] [PMC free article] [PubMed] [Google Scholar]
- Jain K. Transition to turbulence in an oscillatory flow through stenosis. Biomech Model Mechanobiol. 2020;19:113–131. doi: 10.1007/s10237-019-01199-1. [DOI] [PubMed] [Google Scholar]
- Jefferson ND, Cohen AP, Rutter MJ (2016) Subglottic stenosis. 25(3):138–143 [DOI] [PubMed]
- Johannink M, Masilamani K, Mhamdi A, Roller S, Marquardt W. Predictive pressure drop models for membrane channels with non-woven and woven spacers. Desalination. 2015;37(6):41–54. doi: 10.1016/j.desal.2015.07.024. [DOI] [Google Scholar]
- Junk M, Klar A, Luo L-S. Asymptotic analysis of the lattice Boltzmann equation. J Comput Phys. 2005;21(02):676–704. doi: 10.1016/j.jcp.2005.05.003. [DOI] [Google Scholar]
- Junk M, Yang Z (2011) Asymptotic analysis of lattice Boltzmann outflow treatments. Commun Comput Phys 9(5):1117–1127
- Klimach H, Jain K, Roller S. End-to-end parallel simulations with apes. Parallel Comput Accel Comput Sci Eng (CSE) 2014;25:703–711. [Google Scholar]
- Linninger A, Tangen K, Hsu C-Y, Frim D. Cerebrospinal fluid mechanics and its coupling to cerebrovascular dynamics. Ann Rev Fluid Mech. 2016;48(1):219–257. doi: 10.1146/annurev-fluid-122414-034321. [DOI] [Google Scholar]
- Marié S, Ricot D, Sagaut P. Comparison between lattice Boltzmann method and Navier–Stokes high order schemes for computational aeroacoustics. J Comput Phys. 2009;22(84):1056–1070. doi: 10.1016/j.jcp.2008.10.021. [DOI] [Google Scholar]
- Pope SB. Turbulent flows. Cambridge: Cambridge University Press; 2000. [Google Scholar]
- Qi J, Jain K, Klimach H, Roller S (2016) Performance evaluation of the LBM solver Musubi on various HPC architectures. In Parallel Computing: On the Road to Exascale: 807–816. IOS Press
- Qureshi AI, Alexandrov AV, Tegeler CH, Hobson RW, Dennis Baker J, Hopkins LN. Guidelines for screening of extracranial carotid artery disease: a statement for healthcare professionals from the multidisciplinary practice guidelines committee of the american society of neuroimaging; cosponsored by the society of vascular and interventional neurology. J Neuroimaging. 2007;17(1):19–47. doi: 10.1111/j.1552-6569.2006.00085.x. [DOI] [PubMed] [Google Scholar]
- Roller S, Bernsdorf J, Klimach H, Hasert M, Harlacher D, Cakircali M, $\dots $ Zudrop J (2012) An adaptable simulation framework based on a linearized octree. In: High performance computing on vector Systems 2011 High performance computing on vector systems 2011, pp93–105
- Samuelsson J, Tammisola O, Juniper M. Breaking axi-symmetry in stenotic flow lowers the critical transition Reynolds number. Phys Fluids (1994-present) 2015;27(10):104103. doi: 10.1063/1.4934530. [DOI] [Google Scholar]
- Sarpkaya T. Experimental determination of the critical Reynolds number for pulsating Poiseuille flow. J Fluids Eng. 1966;88(3):589–598. [Google Scholar]
- Varghese SS, Frankel SH, Fischer PF. Direct numerical simulation of stenotic flows. Part 1. Steady flow. J Fluid Mech. 2007;58(2):253–280. doi: 10.1017/S0022112007005848. [DOI] [Google Scholar]
- Varghese SS, Frankel SH, Fischer PF. Direct numerical simulation of stenotic flows. Part 2. Pulsatile flow. J Fluid Mech. 2007;58(2):281–318. doi: 10.1017/S0022112007005836. [DOI] [Google Scholar]
- Xu D, Warnecke S, Song B, Ma X, Hof B. Transition to turbulence in pulsating pipe flow. J Fluid Mech. 2017;83(1):418–432. doi: 10.1017/jfm.2017.620. [DOI] [Google Scholar]
- Yellin EL. Laminar–turbulent transition process in pulsatile flow. Circ Res. 1966;19(4):791–804. doi: 10.1161/01.RES.19.4.791. [DOI] [PubMed] [Google Scholar]
- Yoganathan A, Ball J, Woo Y-R, Philpot E, Sung H-W, Franch R, Sahn D. Steady flow velocity measurements in a pulmonary artery model with varying degrees of pulmonic stenosis. J Biomech. 1986;19(2):129–146. doi: 10.1016/0021-9290(86)90143-0. [DOI] [PubMed] [Google Scholar]

## Associated Data

*This section collects any data citations, data availability statements, or supplementary materials included in this article.*

### Supplementary Materials

Download video file^{ (14MB, mp4) }

Supplementary material 1: Vorticity in the axisymmetric case with 60% stenosis and pulsation frequency w3 (MP4 14381 kb)

Download video file^{ (12.8MB, mp4) }

Supplementary material 2: Vorticity in the axisymmetric case with 60% stenosis and pulsation frequency w2 (MP4 13101 kb)

Download video file^{ (11.7MB, mp4) }

Supplementary material 3: Vorticity in the axisymmetric case with 60% stenosis and pulsation frequency w1 (MP4 11951 kb)

Download video file^{ (9.5MB, mp4) }

Supplementary material 4: Vorticity across the x-z plane in the eccentric case with 60% stenosis and pulsation frequency w1 (MP4 9778 kb)

Download video file^{ (14.2MB, mp4) }

Supplementary material 5: Vorticity across the x-z plane in the eccentric case with 60% stenosis and pulsation frequency w3 (MP4 14491 kb)

Download video file^{ (13MB, mp4) }

Supplementary material 6: Vorticity across the x-z plane in the eccentric case with 60% stenosis and pulsation frequency w2 (MP4 13308 kb)