[gray]0.7 Accepted for Publication in Medical & Biological Engineering &
Computing 0.9

#  Efficacy of the FDA nozzle benchmark and the lattice Boltzmann method for
the analysis of biomedical flows in transitional regime

Kartik Jain Corresponding Author; E-mail<EMAIL_ADDRESS>Faculty of
Engineering Technology, University of Twente, P.O. Box 217, 7500AE Enschede,
The Netherlands

###### Abstract

Flows through medical devices as well as in anatomical vessels despite being
at moderate Reynolds number may exhibit transitional or even turbulent
character. In order to validate numerical methods and codes used for
biomedical flow computations, the U.S. food and drug administration (FDA)
established an experimental benchmark, which was a pipe with gradual
contraction and sudden expansion representing a nozzle. The experimental
results for various Reynolds numbers ranging from $500$ to $6500$ were
publicly released. Previous and recent computational investigations of flow in
the FDA nozzle found limitations in various CFD approaches and some even
questioned the adequacy of the benchmark itself. This communication reports
the results of a lattice Boltzmann method (LBM) based direct numerical
simulation (DNS) approach applied to the FDA nozzle benchmark for transitional
cases of Reynolds numbers $2000$ and $3500$. The goal is to evaluate if a
simple _off the shelf_ LBM would predict the experimental results without the
use of complex models or synthetic turbulence at the inflow. LBM computations
with various spatial and temporal resolutions are performed – in the
extremities of $44$ million to $2.8$ billion lattice cells – conducted
respectively on $32$ CPU cores of a desktop to more than $300^{\prime}000$
cores of a modern supercomputer to explore and characterize miniscule flow
details and quantify Kolmogorov scales. The LBM simulations transition to
turbulence at a Reynolds number $2000$ like the FDA’s experiments and
acceptable agreement in _jet breakdown locations, average velocity, shear
stress_ and _pressure_ is found for both the Reynolds numbers.

keywords— Lattice Boltzmann Method, transitional flow, turbulence, FDA,
nozzle, hydrodynamic instability

##  1 Introduction

Computational fluid dynamics (CFD) has seen a growing interest in the
biomedical community as flows within the cardiovascular system or in medical
devices can be computed with sufficient ease to analyze various clinically and
physiologically relevant details. However, the results from CFD can only be
relied if the methods and software tools have been comprehensively verified
and validated. The presence of turbulence in devices as well as in
physiological flows poses additional challenges for the validation of CFD
tools due to the complexity of the transitional flow physics itself.
Furthermore, an appropriate quantification of such a flow regime requires
fully resolved direct numerical simulation (DNS).

The U.S. food and drug administration (FDA) established a benchmark in $1999$
for this purpose 8. This benchmark was developed to contain features that
could closely resemble those in medical devices including regions of flow
contraction and expansion, flow recirculation, local high shear stresses as
well as flow regimes ranging from laminar to transitional and turbulent. In
addition, the benchmark was ensured to be simple enough such that CFD analyses
could be readily performed and compared with experiments. The benchmark was
thus a cylindrical pipe with a gradual contraction and sudden expansion.
Particle image velocimetry (PIV) experiments were conducted by several
laboratories on this device with Reynolds numbers $500$, $2000$, $3500$,
$5000$ and $6500$ to study all the flow regimes namely laminar, transitional
and fully developed turbulence. A large interlaboratory variability was found
in experiments especially in cases with transitional flow regimes
($Re_{th}=2000\,\&\,3500$) 26. The experimental datasets were publicly
released with the intention of serving as a benchmark for the validation of
CFD solvers thereof.

Author group |  $Re_{th}$ $500$ |  $Re_{th}$ $2000$ |  $Re_{th}$ $3500$ |  $Re_{th}$ $5000$ |  $Re_{th}$ $6500$ | Numerical method |  Turbulence model | $Re_{crit}$
---|---|---|---|---|---|---|---|---
Present work |  $\Theta$ |  ✓ |  ✓ |  $\Theta$ |  $\Theta$ | LBM DNS |  none | 2000
Fehn et al.,  7 |  ✓ |  ✓ |  ✓ |  ✓ |  ✓ | high order DG |  none |  $\sim$ 2400
Bergersen et al.,  1 |  $\Theta$ |  $\Theta$ |  ✓ |  $\Theta$ |  $\Theta$ | low order FEM |  none | 3500
Nicoud et al.,  19 |  ✓ |  ✓ |  ✓ |  ✓ |  $\Theta$ | fourth order FVM |  Sigma | 3500
Zmijanovic et al.,  34 |  ✓ |  $\Theta$ |  ✓ |  $\Theta$ |  $\Theta$ | fourth order FVM |  Sigma | 3500
Chabannes et al.,  4 |  ✓ |  ✓ |  ✓ |  $\Theta$ |  $\Theta$ | Taylor-Hood FEM |  none | 2000
Janiga,  15 |  $\Theta$ |  $\Theta$ |  $\Theta$ |  $\Theta$ |  ✓ | low order FVM |  Smagorinsky | $\Theta$
Passerini et al.,  20 |  ✓ |  ✓ |  ✓ |  $\Theta$ |  $\Theta$ | low order FEM |  none | 2000
Delorme et al.,  6 |  ✓ |  ✓ |  ✓ |  ✓ |  $\Theta$ | high order FD, IBM |  Vreman | 2000
White & Chong,  32 |  ✓ |  $\Theta$ |  $\Theta$ |  $\Theta$ |  $\Theta$ | LBM |  none | $\Theta$
Table 1: Computational studies of flow in the FDA nozzle geometry in reverse
chronological order. $Re_{th}$ refers to the Reynolds number at the throat of
the nozzle. The signs ✓and $\Theta$ respectively refer whether or not a
particular $Re_{th}$ was studied by the corresponding author group.
$Re_{crit}$ indicates the Reynolds number at which flow transitioned in the
corresponding study. Part of the information is extracted from Fehn et al.,  7

Various numerical schemes ranging from finite differences to finite element
and finite volume methods as well as high order discontinuous Galerkin (DG)
methods have been employed by researchers to evaluate if their numerical
method and the corresponding implementation can reproduce the results
published by the FDA. Table 1 lists major numerical studies conducted on the
FDA benchmark and the Reynolds number ($Re_{crit}$) at which the corresponding
study found flow transition. In particular researchers evaluate if their CFD
technique can predict the jet breakdown location and quantities like pressure
and velocity accurately. Many studies explore parameters and appropriate mesh
densities at which their results would match with the experimental data. Some
of these studies have questioned the suitability of the FDA benchmark itself,
inquiring if the benchmark should contain more information and be more robust
for comparison with CFD. More interestingly some groups have found flow to
transition to turbulence at $Re_{th}=2000$ 6, 20, 4 whereas others have found
flow transition only at $Re_{th}=3500$ 34, 1.

The study of Fehn et al.,  7 is the only exhaustive work to date in which
authors applied a high order DG method to study all the $Re_{th}$, and even
studied additional $Re_{th}$ to explore the $Re_{crit}$ i.e. the critical
Reynolds number at which the flow would transition in the nozzle. Several
others demonstrated different jet breakdown locations compared to experiments
20, 1, and some focused on outcomes like the _ever changing_ location of jet
breakdown with increasing spatial resolutions 1. Furthermore, it is
interesting to note that in some of the previous studies 34 adding synthetic
fluctuations to the inflow was necessary to make numerical simulations agree
with experiments for the $Re_{th}=3500$ case.

The lattice Boltzmann method (LBM), which is an alternative and relatively new
technique for the numerical solution of Navier-Stokes equations (NSE) has been
extensively used in the past decade by fluid dynamics researchers, and it has
found particular interest in the biomedical community 29, 33, 2, 13 as it can
represent complex anatomical geometries with ease and can enable simulations
on massively parallel computing architectures 13. Recent works 13, 12 applied
and validated LBM for complex transitional flows in anatomical geometries and
found the method efficient and valid for such flow regimes. Despite pressing
needs, however, no extensive effort to the author’s knowledge has been made to
evaluate its efficacy in the aforesaid FDA nozzle benchmark. On the other
hand, due to the diversity of results from different CFD studies and varied
opinions about the suitability of the FDA benchmark, it becomes imperative to
explore the suitability of the benchmark itself using a method that has not
been applied to it before. White & Chong,  32 did apply LBM to the FDA nozzle
benchmark but only for the laminar case and their study was more oriented
towards exploring the suitability of different lattice types. The previous
works of transitional flow computations using LBM 13, 12 were also for
moderate Reynolds numbers.

Motivated by the aforementioned needs, this work aims to evaluate if a simple
LB scheme, without the employment of complex collision models or synthetic
turbulence at the inflow, will accurately predict the results benchmarked by
the FDA. The focus is on transitional flow regime and thus Reynolds numbers
$2000$ and $3500$ are only analyzed. Physical quantities like _velocity, shear
stress_ and _pressure_ and observations like the _jet breakdown_ location are
compared from experiments and simulations. Furthermore, insight into questions
like _when (Re)_ , _where (locations)_ , _whether_ and _how_ of flow
transition is provided. These goals are achieved by conducting LBM based
direct numerical simulations (DNS) at a very high, and another extreme
resolution, which is found to be below the scales defined by Kolmogorov
theory. A comparison against Kolmogorov theory and assessment of Kolmogorov
scales is further shown to elucidate the role resolutions play in computation
of complex flow in such a medical device.

##  2 Methods

###  2.1 The FDA nozzle

Figure 1:  Cross sectional view of the FDA nozzle for the sudden expansion
case. Results were analyzed at various z-locations marked by dotted vertical
lines to compare against FDA’s experimental results:
$z_{1}=-0.088m,\,z_{2}=-0.064m,\,z_{3}=-0.048m,\,z_{4}=-0.02m,\,z_{5}=-0.008m,\,z_{6}=0.0m,\,z_{7}=0.008m,\,z_{8}=0.016m,\,z_{9}=0.024m,\,z_{10}=0.032m,\,z_{11}=0.06m,\,z_{12}=0.08m$.

A cross sectional view of the FDA nozzle is shown in figure 1. The flow
direction is from left to right implying a sudden expansion. If the flow is
applied from right to left the same geometry will act as a canonical diffuser,
which is not studied in this work. The 3D model of this nozzle was downloaded
from the FDA website111https://ncihub.org/wiki/FDA_CFD. The outlet of the
model was extended up to $z=0.253$ m and the inlet was extended up to
$z=-0.14$ m. The radial profiles were recorded at $12$ locations depicted in
figure 1 to enable a comparison against the particle image velocimetry (PIV)
experimental data released by the FDA.

Simulations were conducted with three different spatial and temporal
resolutions in order to explore the capabilities of an _off the shelf_ LBM on
a desktop machine to a federal supercomputer. Based on previous observations
on resolution requirements in LBM simulations of transitional flow 14 a
spatial resolution of about $80\mu m$ was predicted as a minimum requirement
for simulation of $Re_{th}=2000$ and accurate analysis of flow characteristics
thereof. This resolution resulted in $\sim 45$ million lattice sites
representing fluid inside the nozzle and $50$ lattice cells across the height
of the nozzle throat. This is referred to as normal resolution (NR) hereon.
For $Re_{th}=3500$ however, the flow was expected to transition and a high
resolution (HR) of $40\mu m$ was thus chosen, which resulted in about $357$
million lattice cells. As would be seen in section 3, HR resolution sufficed
for accurate flow simulations for both the Reynolds number. To further ensure
mesh independence, assess the flow in accordance with Kolmogorov theory, and
especially because it was not known whether LBM simulations will transition to
turbulence, an additional set of simulations were conducted with an _extreme_
resolution (XR) of $20\mu m$ resulting in about $2.88$ billion lattice cells
in the fluid domain.

| $\delta x(\times 10^{-6}m)$ | $\delta t(\times 10^{-6}s)$ |  nCells${}_{\text{th}}$ | nCells(M) | nCores | T(h)
---|---|---|---|---|---|---
NR | $80$ | $16$ | $50$ | $45$ | $32$ | $192$
HR | $40$ | $4$ | $100$ | $357$ | $38^{\prime}016$ | $3.5$
XR | $20$ | $1$ | $200$ | $2^{\prime}880$ | $304^{\prime}128$ | $12$
Table 2: Three different spatial and temporal resolutions and the
corresponding number of lattice cells inside the fluid domain. The NR
simulations were conducted on $32$ cores of a desktop while others on the
_SuperMUC-NG_. A total of $10$ physical seconds were simulated for all the
resolutions.

Table 2 lists the employed spatial and temporal resolutions, the resulting
number of lattice cells, the utilized CPUs and total execution time of the
simulation.

###  2.2 Setup of a simulation using the lattice Boltzmann method

To shed light on the aforementioned choice of resolutions and LB parameters a
brief explanation is provided here. The LBM is based on the mesoscopic
representation of movement of fictitious particles. The particles have
discrete velocities and they collide and stream to relax towards a
thermodynamic equilibrium. The LB equation recovers the NSE under the
continuum limits of low Mach and Knudsen numbers. Evolution of the particle
probability distribution functions over time is described by the Lattice
Boltzmann equation with the MRT collision matrix:

$\displaystyle f_{i}(\mathbf{r}+\mathbf{c}_{i}\delta t,t+\delta
t)=f_{i}(\mathbf{r},t)+\Omega_{ij}\left(f_{i}^{e}(\mathbf{r},t)-f_{i}(\mathbf{r},t)\right)$
(1)

where $f_{i}$ represents the density distributions of particles which are
moving with discrete velocity $\mathbf{c}_{i}$ at a position $\mathbf{r}$ at
time t. The indices which run from $i\\!=\\!1\\!\ldots\\!Q$ denote the links
per element i.e. the discrete directions, depending on the chosen stencil
(D3Q19 in our case). The collision matrix $\Omega_{ij}$ defines relaxation of
various modes of the distribution functions $f_{i}$ towards an equilibrium
$f_{i}^{e}$:

$\displaystyle
f_{i}^{e}=w_{i}\rho\Bigg{(}1+\frac{\mathbf{c}_{i}\cdot\mathbf{u}}{c_{s}^{2}}-\frac{\mathbf{u}^{2}}{2c_{s}^{2}}+\frac{1}{2}\frac{(\mathbf{c}_{i}\cdot\mathbf{u})^{2}}{c_{s}^{4}}\Bigg{)}$
(2)

where $w_{i}$ are the weights for each discrete link, $c_{s}$ is the reference
speed of sound in LBM obtained by integration of the discrete Boltzmann
equation along characteristics, and $\mathbf{u}$ is the fluid velocity. The
time step in LBM is coupled with the grid size by $\delta t\sim\delta x^{2}$
due to _diffusive_ scaling which is employed to recover the incompressible
NSE. Details on the computation of macroscopic quantities from LBM can be
found elsewhere 28; here we restrict ourselves on the process of setting up a
flow simulation.

####  2.2.1 Prescription of flow physics

Quantities like the mean velocity, $\bar{u}$, on which the Reynolds number is
based, fluid density $\rho$, and the kinematic viscosity $\nu$ are firstly
chosen for the fluid to be simulated. The Reynolds number is then calculated
as:

$\displaystyle\mathrm{Re}=\frac{\bar{u}\,d}{\nu}$ (3)

where d is the characteristic length equal to the throat diameter ($0.004$m)
in this case. In this study the Reynolds number is based on the nozzle throat.

The LB method requires the translation of these quantities into lattice units.
Under diffusive scaling, the principal relaxation parameter $\Omega$ is fixed
and can have a maximum value of $2.0$ for stability. The lattice viscosity,
$\nu_{lattice}$ is then calculated using the relation:

$\displaystyle\nu_{lattice}=c_{s}^{2}\bigg{(}\frac{1}{\Omega}-\frac{1}{2}\bigg{)}$
(4)

where $c_{s}^{2}$ is the speed of sound squared at reference state and is
equal to $\frac{1}{3}$ for the chosen stencil D3Q19. The time step is then
calculated as:

$\displaystyle\delta t=\frac{\nu_{lattice}\,\delta x^{2}}{\nu}$ (5)

Finally, the lattice velocity is obtained as:

$\displaystyle\bar{u}_{lattice}=\bar{u}\,\delta t/\delta x$ (6)

The MRT collision operator allows for different relaxation rates of various
moments of the distribution function, and it was employed throughout for the
DNS reported in this manuscript. The principal relaxation parameter that
defines the kinematic viscosity was set to $\Omega=1.90$. It is required that
the $\bar{u}_{lattice}$ in equation (6) is always less than 0.15 to enforce
the Mach number limit of the LBM. The $\Omega$ can be adjusted in order to
fine tune the lattice velocity and beyond its limit, the grid has to be
refined (reduce the $\delta x$), which is also the principle limitation of the
Lattice Boltzmann Method. Thus, the $\delta t$ is controlled by the $\delta x$
and the $\Omega$, and it is further constrained by the fluid velocity. To
assert correct prescription of parameters, if the characteristic length in
equation (3) is replaced by the number of fluid cells along that particular
characteristic length and $\bar{u}$ and $\nu$ are replaced by
$\bar{u}_{lattice}$ and $\nu_{lattice}$ respectively, the same Reynolds number
must be obtained as that obtained from equation (3), which uses physical
values.

####  2.2.2 The simulation framework

The employed simulation tool-chain is contained in the end-to-end parallel
framework APES (adaptable poly engineering simulator) 24,
17222https://apes.osdn.io. Meshes are created using the mesh generator
_Seeder_ 9 and computations are carried out using the LBM solver _Musubi_ 10.
_Musubi_ writes out binary files containing physics information to the disk.
These files are then converted to the visualization toolkit (VTK) format by
the post-processing tool _Harvester_ , which is contained within the APES
framework. The open source visualization tool
Paraview333https://www.paraview.org is then used to visualize the physics of
flow. The data for plots is written out by _Musubi_ as ASCII files that are
plotted using the Matplotlib plotting library within the Python programming
language.

The 3D model of the nozzle in STL format was read by the mesh generator
_Seeder_ and volume meshes for LBM computations were saved on the disk. A
higher order wall boundary condition described by Bouzidi et al.,  3 was
prescribed at the walls of the nozzle to reduce the influence of staircase
artifacts in LBM and ensure rotational symmetry of the setup. The _Musubi_ LBM
solver then computed these meshes (see table 2) on the _SuperMUC-NG_ petascale
system installed at the Leibniz Supercomputing Center in Munich, Germany. The
number of utilized cores on _SuperMUC-NG_ ranged from $38\,016$ to $304\,128$
(the whole system) based on our previous evaluations which have shown that
this choice of CPUs results in an optimal utilization of compute resources 17,
23.

####  2.2.3 Initial and boundary conditions

Initial conditions were set to zero pressure and velocity and a no-slip
boundary condition was maintained at the walls. The flow rates at inlet were
characterized based on the throat Reynolds number. The fluid itself was chosen
to be Newtonian, with a fluid density of $1056kg/m^{3}$ and dynamic viscosity
of $0.0035Ns/m^{2}$.

A parabolic velocity profile was prescribed at the inlet by quadratically
interpolating the velocity at lattice nodes, and a zero pressure was
maintained at the outlet. No synthetic turbulence was prescribed at the inlet
unlike some previous studies. The outflow in such flow regimes can be critical
resulting in back flow or instabilities in the LBM algorithm for which an
extrapolation boundary condition was employed 16.

Following Fehn et al.,  7 the simulated time interval was chosen on the basis
of a characteristic flow-through time. The length of the throat, $L_{th}$
section and the mean velocity, $\bar{u}$ were chosen as characteristic length
scale and the reference velocity respectively. This resulted in a flow through
time given by:

$\displaystyle T_{ft}=\frac{L_{th}}{\bar{u}}$ (7)

For the $Re_{th}=2000$ case this time was approximately $0.21$ s whereas for
the $Re_{th}=3500$ case it was $0.12$ seconds. Flow was allowed to develop for
initial $2$ seconds implying about $10$ and $20$ flow throughs for the
$Re_{th}=2000$ and $Re_{th}=3500$ case respectively. Within the computational
bounds a total of $10$ physical seconds were simulated for each case. A total
of $50000$ samples were gathered every second for the analysis of flow
characteristics.

###  2.3 Flow characterization

For the analysis of turbulent characteristics, the three dimensional velocity
field was decomposed into a mean and fluctuating component as:

$\displaystyle u_{i}(x,t)=\bar{u}_{i}(x)+u_{i}^{\prime}(x,t)$ (8)

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

$\displaystyle k=\frac{1}{2}\Big{(}{u_{x}^{\prime 2}+u_{y}^{\prime
2}+u_{z}^{\prime 2}}\Big{)}$ (9)

The dimensionless Strouhal number, St is defined as:

$\displaystyle St=\frac{fd}{\bar{u}}$ (10)

where f is the frequency of flow fluctuations, d and $\bar{u}$ are the
characteristic length and mean velocities. A power spectral density of the TKE
was computed as a function of the Strouhal number at various centerline
locations using the Welch’s periodogram method to analyze the inertial and
viscous subranges.

###  2.4 Kolmogorov microscales

The smallest structures that can exist in a turbulent flow can be estimated on
the basis of Kolmogorov’s theory 22 – used in this work to assess the quality
of the resolutions in DNS. Viscosity dominates and the TKE is dissipated into
heat at the Kolmogorov scale 22. The Kolmogorov theory, in general, applies to
fully developed turbulent flows with Reynolds numbers much higher than those
studied in this work. The reference to Kolmogorov scales and theory here is
thus indicative for the assessment of mesh independence as has been done in
various such studies at low Reynolds numbers 12, 11.

The Kolmogorov scales, non-dimensionalized with respect to the velocity scale
$\bar{u}$ and the length scale d (throat diameter) are computed from the
_fluctuating_ component of the non-dimensional strain rate defined as:

$\displaystyle s^{\prime}_{ij}=\frac{1}{2}\bigg{(}\frac{\partial
u_{i}^{\prime}}{\partial x_{j}}+\frac{\partial u_{j}^{\prime}}{\partial
x_{i}}\bigg{)}$ (11)

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

$\displaystyle\eta=\bigg{(}\frac{1}{Re^{2}}\frac{1}{2s^{\prime}_{ij}s^{\prime}_{ij}}\bigg{)}^{1/4}$
(12)
$\displaystyle\tau_{\eta}=\bigg{(}\frac{1}{2s^{\prime}_{ij}s^{\prime}_{ij}}\bigg{)}^{1/2}$
(13) $\displaystyle
u_{\eta}=\bigg{(}\frac{2s^{\prime}_{ij}s^{\prime}_{ij}}{Re^{2}}\bigg{)}^{1/4}$
(14)

Based on these scales, the quality of the spatial and temporal resolutions
employed in a simulation can be estimated as the ratio of $\delta x$ and
$\delta t$ against the corresponding Kolmogorov scales i.e.

$\displaystyle l^{+}=\frac{\delta x}{\eta}\;\hskip 10.0ptt^{+}=\frac{\delta
t}{\tau_{\eta}}$ (15)

The Kolmogorov scales in this study were computed along various locations at
the centerline. The fluctuating strain rate was averaged between locations
$z_{8}$ and $z_{12}$ (figure 1) due to variations in TKE caused by jet
breakdown in these locations, and the Kolmogorov scales were computed from
this averaged strain rate.

###  2.5 Comparison to experiments

The radial velocity profiles in the streamwise direction and the mean
centerline velocity at $12$ locations shown in figure 1 were plotted against
corresponding data from $5$ PIV experiments. Instantaneous centerline
velocities and pressure were analyzed at $120$ locations each $0.002$ m apart
along the length of the nozzle. The $5$ FDA experiments were conducted by $3$
independent laboratories where one laboratory ran three trials resulting in
$5$ data sets 8. The axial component of the velocity $u_{z}$ was plotted
directly against the data from experiments444Note that other works 7, 20 have
normalized the $u_{z}$ with respect to $\bar{u}_{in}$.

Pressure was probed and averaged at $12$ circular cross sections across the
length of the nozzle ($z_{1}\,-\,z_{12}$). To enable comparison against
experiments, the mean pressure difference was normalized with respect to the
average velocity at the throat $\bar{u}$:

$\displaystyle\Delta p^{norm}=\frac{p_{z}-p_{z_{0}}}{0.5\rho\bar{u}^{2}}$ (16)

where $\rho$ is the fluid density and $\bar{u}$ is the mean velocity at the
nozzle throat, on which the $Re_{th}$ was based (see equation (3)).

The shear stress was computed from the LBM simulations as it is one of the
quantities of interest to predict hemolysis, and was directly compared to the
experimental data.

The differences between a particular LBM simulation case and an experimental
dataset were quantified on the basis of simple relative percentage errors,
computed as:

$\displaystyle\delta=\left|~{}\frac{U_{ref}-U_{h}}{U_{ref}}~{}\right|\times
100$ (17)

where $U_{ref}$ denotes the reference solution (experiments in our case) and
$U_{h}$ denotes the LBM computed solution; the corresponding solutions can be
for pressure, velocity or shear stress.

##  3 Results

The flow transitioned to turbulence at $Re_{th}=2000$ and it became fully
turbulent at $Re_{th}=3500$. The NR simulation was stable for $Re_{th}=2000$
during the whole $10$ seconds that were simulated. For $Re_{th}=3500$,
however, the NR simulation crashed as soon as the flow reached the outflow due
to the instabilities in that region, and local fluctuations in the Mach number
limit of the LB algorithm (equation 6).

#### Kolmogorov microscales

To assess the mesh independence of the simulations we start with an analysis
of the Kolmogorov microscales, shown for different resolutions and $Re_{th}$
in table 3.

| $Re_{th}=2000$ | $Re_{th}=3500$
---|---|---
|  $\eta(\mu m)$ |  $\tau_{\eta}(\mu s)$ |  $u_{\eta}(m/s)$ |  $l^{+}$ |  $t^{+}$ |  $\eta(\mu m)$ |  $\tau_{\eta}(\mu s)$ |  $u_{\eta}(m/s)$ |  $l^{+}$ |  $t^{+}$
NR |  $4.94$ |  $8.37$ |  $1.19$ |  $16.19$ |  $1.91$ |  $\Theta$ |  $\Theta$ |  $\Theta$ |  $\Theta$ |  $\Theta$
HR |  $16.32$ |  $8.69$ |  $1.97$ |  $2.45$ |  $0.46$ |  $12.57$ |  $4.08$ |  $3.80$ |  $3.18$ |  $0.98$
XR |  $21.97$ |  $7.14$ |  $2.10$ |  $0.91$ |  $0.14$ |  $18.86$ |  $4.34$ |  $4.10$ |  $1.06$ |  $0.23$
Table 3: Kolmogorov scales computed from NR, HR and XR resolutions for
$Re_{th}=2000$ and $Re_{th}=3500$.

The $l^{+}$ from NR resolution for $Re_{th}=2000$ is $16.19$ thus clearly
indicating this resolution as insufficient for accurate assessment of
turbulent characteristics. These scales from the HR simulations for both
$Re_{th}$ are sufficiently close to the corresponding Kolmogorov scales while
from XR simulations they attain a value equal to the Kolmogorov scales. As
observed later, the XR resolutions, being resolved exactly at the Kolmogorov
length scales result in minor deviances in flow characteristics if compared to
HR resolutions. The $\tau_{\eta}$ are $<1$ in most of the cases due to the
small time step of LB simulations.

#### Jet breakdown location

(a) $Re_{th}=2000$ NR

(b) $Re_{th}=2000$ HR

(c) $Re_{th}=2000$ XR

(d) $Re_{th}=3500$ HR

(e) $Re_{th}=3500$ XR

Figure 2: Snapshots of the instantaneous vorticity during t=10 second of the
simulation across a bisecting plane in the FDA nozzle from HR and XR LBM
simulations for $Re_{th}=2000$ and $Re_{th}=3500$. The vorticity magnitude is
scaled according to $Re_{th}$ and ranges from $0-2Re_{th}$.

(a) $Re_{th}=2000$

(b) $Re_{th}=3500$

Figure 3: Spectral density of the turbulent kinetic energy at $6$ locations
along the centerline after the sudden expansion ($z_{7}-z_{12}$ in figure 1)
along the centerline is plotted against the Strouhal number (St). The PSD is
computed during the last $8$ seconds of the simulations using $4\times 10^{5}$
time steps in total. Dark and dotted gray lines show the PSD respectively from
HR and XR LBM simulations. For $Re_{th}=2000$ PSD from NR simulations is shown
in green line. The solid gray line shows the Kolmogorov’s $\frac{-5}{3}$
decay.

The regions of flow breakdown and maximum turbulent activity are highlighted
in figure 2, which shows the instantaneous vorticity magnitude at the $t=10$
second of simulations. Jet breakdown location identified from NR simulations
of $Re_{th}=2000$ lies between $z_{9}$ and $z_{10}$. A secondary flow jet
however starts emanating right after the expansion. It may be concluded that
the NR resolution captures the onset of turbulence but with a largely
compromised accuracy as the location of jet breakdown is different from all
the previously reported experiments 26. On the other hand the vorticity plots
for $Re_{th}=2000$ from both HR and XR LBM simulations demonstrate that the
location of jet breakdown is more or less identical (between $z_{11}$ and
$z_{12}$) and resolutions do not seem to influence the location of jet
breakdown. The jet breakdown location is between $z_{9}$ and $z_{10}$ at
$Re_{th}=3500$ from HR simulations while it is shifted further downstream by
approximately half of the throat diameter from XR simulations. It is
emphasized that these are _instantaneous_ vorticity plots during the t=10
second of the simulation as due to immense memory requirements an ensemble
averaging was not possible. Furthermore, it is expected that an ensemble
average of at least $10\,000$ snapshots would be needed to have an accurate
assessment of jet breakdown location.

In the supplementary material, animations of vorticity and velocity field for
both $Re_{th}$ from HR simulations are provided over the last one second
($9-10$) of the simulations. It may be observed from the animations that at
$Re_{th}=2000$ the jet itself experiences distortion over time course of the
simulation. If we take a closer look at the figure 2b we can see development
of _discontinuities_ over the jet after $3$ nozzle diameters downstream of
expansion. When the flow breaks down, vortices merge, annihilate and interact
with the jet to _distort_ the jet itself over the course of time and
resolution does not seem to play a role here. This circumstance is in fact
natural for a transitional flow as it is not fully developed turbulence but
the jet is between laminar and transitional regime. For $Re_{th}=3500$, on the
other hand, the jet does not get distorted over the course of time but it
shifts further downstream at XR resolution (figure 2e). In this case, a shear
layer develops around the jet itself as can be clearly seen from the
animations as well as the instantaneous snapshots of vorticity. A further grid
refinement for the case of $Re_{th}=3500$ (figure 2e) resolves this shear
layer better to shift the jet breakdown location further downstream by about
half a throat diameter ($d/2=0.002m$) as this is a fully developed turbulent
flow field.

An accurate assessment of turbulent activity at various locations, and from
different resolutions can be obtained from the PSD plots of figure 3, computed
at locations downstream of the expansion ($z_{7}-z_{12}$ from figure 1). The
PSD is computed over $8$ seconds of the simulation to obtain abundant
statistics and overcome previous issues of insufficient averaging. The dark
and dotted gray lines respectively show PSD from HR and XR LBM simulations
whereas the solid gray line indicates Kolmogorov’s $\frac{-5}{3}$ decay. For
$Re_{th}=2000$ the green line shows the PSD from NR simulations, which
corroborates previous observations from vorticity plots and Kolmogorov scales
and confirms the inadequacy of this resolution. For the $Re_{th}=2000$ case it
is seen that the flow transitions at locations $z_{11}$ and $z_{12}$ and the
turbulent activity captured by HR and XR simulations is the same with a few
frequencies in the inertial subrange. The third plot of figure 3b is the most
enlightening about the jet breakdown location for $Re_{th}=3500$ from HR and
XR simulations and corroborates the observations from plot 4b. Clearly the jet
does not break down at location $z_{9}$ in XR simulations as the spectrum
falls off rapidly whereas in HR simulations the spectrum tends to attain a
Kolmogorov like decay 12 at this location. The following figures then
substantiate this observation whereby at $z_{10}$ the jet is already turbulent
from both HR and XR simulations and turbulent activity becomes similar beyond
this (last two plots of figure 3b).

#### Comparison of velocity, pressure and shear stress against experiments

(a) $Re_{th}=2000$

(b) $Re_{th}=3500$

Figure 4: Time averaged centerline velocity at $86$ locations along the length
of the nozzle for two different Reynolds numbers and all the employed LBM
resolutions. The centerline velocity is compared with $12$ ($z_{1}$ – $z_{12}$
from figure 1) experimental data points along the centerline available from
$5$ distinct PIV experiments. The number in the legend corresponds to the
experiment ID from the FDA website. The gray bars highlight the regions with
largest discrepancy between experiments and simulations. At $Re_{th}=3500$ the
NR simulation was unsuccessful.

(a) $Re_{th}=2000$

(b) $Re_{th}=3500$

Figure 5: Time averaged centerline pressure normalized by the mean throat
velocity $u_{th}$ versus the axial distance computed from two sets of LBM
resolutions is compared against pressure data from $5$ PIV experiments of the
FDA for two different Reynolds numbers. Experimental data is plotted at $12$
locations along the centerline ($z_{1}$ – $z_{12}$ in figure 1) for
experiments whereas $86$ points between these locations are plotted from
simulations. Note that for $Re_{th}=2000$ case the pressure data from
experiment $\\#999$ is not available. At $Re_{th}=3500$ the NR simulation was
unsuccessful.

Figure 4 shows the averaged centerline velocity computed at $86$ axial
locations between $z_{1}$ and $z_{12}$ for both the Reynolds numbers from HR
and XR LBM simulations and additionally from NR simulations for the
$Re_{th}=2000$ case. The corresponding data from $5$ experiments at $12$
locations is plotted for comparison. The shaded bars in the plots highlight
the regions of largest discrepancy between experiments and simulations. As can
be seen, for $Re_{th}=2000$ there is a significant difference between
experiments and simulations between locations $z=-0.04$ and $z=-0.02$, which
is the region just after the end of gradual contraction, or the beginning of
the nozzle throat. The velocity matches well at the latter part of the throat
and this difference is seen again in the jet breakdown location, which matches
reasonably with one experiment only ($e243$) for both the HR and XR
resolutions.  The difference in velocity in regions between $z=0.06$ and
$z=0.07$ is about $10\%$ between few experiments and simulations, computed
using equation (17).  From NR simulations the jet breakdown location is
completely different from experiments and HR and XR simulations as was seen in
previous plots.

For the $Re_{th}=3500$ case similar trends in the throat area are seen while
the jet breakdown location from HR simulations matches more closely to the
experiments. This location from the XR simulations is comparatively further
away by half the throat diameter.  The velocities estimated by the LBM at the
throat are also higher by about $6-9\%$ than those from the experiments. For
$Re_{th}=3500$, the velocities computed by LBM right after the expansion at
locations around $z=0.025$ are overestimated by about $10\%$ compared to the
experiments.

The cause of discrepancies in the starting of the throat area (shaded regions)
cannot be unequivocally stated as the data points available from the
experiments are very few. In the plot from the simulations, there are $14$
samples between locations $z_{3}$ and $z_{4}$ (figure 1). If the data points
between these two locations were omitted, the curves would look like a
straight line as in the experiments and previously reported studies 7. The
pressure at the throat right after the contraction reduces considerably and
the minor discontinuities at the corners might explain these differences.

This is further elaborated in the centerline average pressures plotted in
figure 5. Qualitatively, similar trends from experiments and simulations are
observed, and the difference between HR and XR simulations for both the
Reynolds number is much smaller.  Interestingly, the pressure drop for
$Re_{th}=3500$ at locations before the expansion ($z<0.0$) is overestimated by
about $5-6\%$ and matches only one experiment.  The data points for pressure
from PIV experiments are relatively few disabling a better and insightful
comparison. Furthermore because the outflow length in LBM simulations is
larger than that in experiments, it may account for minor differences in
pressure changes.

(a) $Re_{th}=2000$

(b) $Re_{th}=3500$

Figure 6: Radial velocity profiles at $10$ locations ($z_{3}-z_{12}$ shown in
figure 1) along the length of the nozzle for two different Reynolds numbers
averaged over $8$ seconds or $400000$ time steps. The velocities are
normalized to enable an intuitive comparison. The black lines show the
averaged velocities computed from XR LBM simulations whereas red, blue, green
and olive colored lines are the data from experiment ids $\\#243$, $\\#297$,
$\\#763$ and $\\#999$ respectively. The gray dotted lines at locations
$z_{11}$ and $z_{12}$ for $Re_{th}=2000$ and at $z_{10}$, $z_{11}$ and
$z_{12}$ for $Re_{th}=3500$ show instantaneous radial velocities from LBM
simulations to depict the location of jet breakdown.

(a) $Re_{th}=2000$

(b) $Re_{th}=3500$

Figure 7: Radial shear stress profiles at $10$ locations ($z_{3}-z_{12}$ shown
in figure 1) along the length of the nozzle for two different Reynolds numbers
averaged over $8$ seconds or $400000$ time steps. The shear stresses are
normalized to enable an intuitive comparison. The black lines show the
averaged velocities computed from XR LBM simulations whereas red, blue, green
and olive colored lines are the data from experiment ids $\\#243$, $\\#297$,
$\\#763$ and $\\#999$ respectively.

Figure 6 shows the radial profiles of mean streamwise velocity from XR LBM
simulations at $10$ locations after the contraction ($z_{3}-z_{12}$ from
figure 1). The ensemble averages are gathered over the last $8$ seconds of the
simulation or for $4\times 10^{5}$ time steps.  In this case the velocities
have been normalized by $\bar{u}$ for a direct intuition about the differences
in experiments and simulations. The red, blue, green and olive colored lines
are the data from experiment ids $\\#243$, $\\#297$, $\\#763$ and $\\#999$
respectively555Plots for experiment $\\#468$ have been omitted as several
radial profiles are not available from the data of this experiment..  At
locations $z_{11}$ and $z_{12}$ for $Re_{th}=2000$ and at $z_{10}-z_{12}$ for
$Re_{th}=3500$, the instantaneous radial profiles at t=10 second of the
simulation are plotted in gray dotted lines to depict regions where the flow
jet broke down. For the $Re_{th}=3500$ case the velocity from LBM in the post
expansion region is higher than the experiments mainly in the jet core region
as was also observed in figure 4b. Except for the jet breakdown locations, a
reasonable agreement in radial profiles can be seen against all the
experiments. As has been observed throughout this study the interlaboratory
variations in the experiments are quite large 26 making a quantitative
comparison more difficult.

Shear stress from experiments and computations is compared in figure 7. Here
the shear stress has been normalized by mean shear stress so that a direct
comparison can be enabled. Excellent agreement for both the Reynolds number
can be observed in mean radial velocities as well as the shear stress666The
shear stress at $z_{3}$ seems to _go out_ of the nozzle. This is a visual
artifact as the shear rises dramatically at the walls.. The shape of the shear
stress profile at locations before the expansion ($z<0.0$) is flatter in the
simulations whereas in experiments it exhibits a skewed shape. Other minor
differences lie in the regions of jet breakdown and are a consequence of the
fact that $4\times 10^{5}$ sample points are used for gathering statistics,
while good, are not perfect for the reproduction of a converged state. Despite
the computational resources deployed for this study, sampling more statistics
was considered prohibitive. For $Re_{th}=3500$, LB overestimates the velocity
and the shear stress at the location of jet breakdown, which was also observed
in figure 4.  These differences are about $3-4\%$ are mostly confined around
the jet breakdown locations.

##  4 Discussion

It is the first time that LBM has been applied to study transitional flows in
the two decade old FDA benchmark. The study has found that LBM can predict the
jet breakdown location and compute macroscopic quantities to a reasonable
accuracy for a transitional flow in a biomedically relevant device. In
addition to a comparison the main physical insight is the distortion in jet
for transitional flow case of $Re_{th}=2000$ and assessment of the Kolmogorov
microscales. Here we discuss the agreement and discrepancy in results as well
as implications of this study.

#### Analysis of the flow characteristics and comparison to previous works and
experiments

Qualitatively the LB computed flow transitioned to turbulence at
$Re_{th}=2000$ as it did in most of the experiments 26, and flow field at
$Re_{th}=3500$ was reminiscent of a partially developed turbulence indicated
by the immense chaotic activity. Quantitatively, a satisfactory match between
experiments and computations was observed in averaged velocities, shear
stress, and pressure as well as the jet breakdown locations. This agreement
has been obtained without parameter tuning in the method, and in particular
without adding synthetic fluctuations as has been necessary for other studies
34. The centerline velocity and pressures computed from LBM had marked
differences, in particular in the throat and the expansion area. In particular
for the $Re_{th}=3500$ case the LB computed velocities were higher by about
$10\%$ in comparison to the experiments. While an exact reason for this is not
known it may be observed that the LB computed radial velocity profiles (figure
6b) had a more pronounced recirculation region and a higher velocity only in
the jet core – not seen in the experiments.

A mesh sensitivity analysis revealed that resolutions within an order of
magnitude of the Kolmogorov microscales are necessary to accurately capture
the flow field and a further refinement down to the Kolmogorov scales results
in slightly distal jet breakdown of flow downstream of expansion for
$Re_{th}=3500$ while the averaged macroscopic quantities are not much
affected. The NR resolution failed for the higher $Re_{th}=3500$ and the
results were erroneous for $Re_{th}=2000$. Previous studies have found a much
larger dependence of results on the mesh density. Delorme et al.,  6 used LES
model and found good match except that their simulations found relatively less
turbulent nature of the flow attributed to grid resolutions. Passerini et al.,
20 performed a number of simulations for each $Re$ they studied (table 1) to
arrive at an optimal grid resolution. Zmijanovic et al.,  34 used three mesh
resolutions in their LES simulations and found excellent match with
experiments from the highest resolution of $50M$ cells. In the present LBM
simulations, going from HR to XR increased the computational effort remarkably
and brought the resolutions right down to the Kolmogorov microscales. The
benefit of this was pronounced at $Re_{th}=3500$ while at $Re_{th}=2000$
almost no advantage in the results was seen from HR to XR. The HR setup is
thus adequate in this configuration as suggested by previous studies 14 as
well. The scalability of LBM to XR scale may be leveraged in future to
incorporate physics of red blood cells or other particles in the blood 29 to
answer relevant questions in physiology.

In the present work the role of boundary conditions has not been explored in
detail and the prescription of inflow and outflow boundary conditions is
likely to influence the results. Passerini et al.,  20 extended the in- and
outflow according to the $Re_{th}$ and the recent precursor approach adopted
by Fehn et al.,  7 seems to overcome most of the effects of inflow pipe
length. When a simulation is refined, the accuracy of the boundary conditions
also increases, which explains the _eternal_ change in solution upon
increasing resolutions. In particular the disturbances that emanate from the
outflow are reduced upon refining the mesh and time step, which would explain
the relatively further breakdown of jet downstream of the expansion.

#### Jet breakdown location

A number of previous studies have found the jet breakdown location sensitive
to parameters used in CFD 20, 34. Here, as noted above there were no
parameters that were varied in the LBM simulations except for the grid
refinement and the corresponding time step size. For $Re_{th}=2000$, we noted
that the jet breakdown location was largely the same from HR and XR
resolutions. A new observation, seemingly overlooked in previous studies was
the propagating distortion in the jet at $Re_{th}=2000$ (animations in
supplementary information). As the flow was just transitioning, the jet lost
its momentum at times, tended to restabilize the flow, and then gained
momentum again resulting in shifts in the breakdown location. For
$Re_{th}=3500$ the jet broke down half a nozzle diameter further downstream of
the expansion in XR simulations compared to HR. The reason for that is the
better accuracy at the sudden expansion, which stabilized the perturbances in
that area, and their propagation thereof. Nonetheless, the breakdown location
from XR simulations was just half the throat diameter, $d/2$, downstream, and
its relevance can be queried. More interestingly the flow restabilization
regions from both the resolutions were the same for $Re_{th}=3500$, and the XR
resolution thus narrowed down the length of chaotic activity (figure 3). It
would be biased and impractical to extrapolate these observations of jet
breakdown location to other numerical studies 20, 34, 7 as the LBM essentially
solves the Boltzmann equation to recover NSE. Also, the parameters that have
been of discussion in related studies 7, 6 like Courant number cannot be
directly related with LBM. It is however important to remark that numerical
dissipation in LBM, even at the scales of grid spacing, and the numerical
dispersive effects are smaller compared to other second-order accurate methods
18, which to a certain extent explains the consistent jet breakdown locations
with increasing resolutions.

In this study, we upsurged from NR to HR and XR directly without exploring a
mid resolution range. It is noted that these are not a minimum resolution
requirement for LBM simulations, and it is thus re-emphasized that this study
was designed to explore LBM’s suitability in performing fully resolved DNS,
and thus resolutions, _as high as possible_ , within the confines of present
computational paradigms were employed. It may be noted that the LB literature
has grown considerably in the past and several ways to incorporate turbulence
models within LB have been devised 27. Evaluations of these models using the
FDA nozzle is left for future efforts.

#### Onset of flow transition

If we focus on the FDA nozzle only, previous studies have found conflicting
$Re_{crit}$, and why the numerically computed flow field did, or did not
transition at $Re_{th}=2000$ in one particular study or the other is still
unknown. It is fair to state that a flow in a perfectly symmetric setup as
that of the FDA benchmark should not transition to turbulence in a simulation,
and such an event must be a consequence of the numerics. Different numerical
methods, parameters, stabilization techniques and resolutions may thus result
in suppression or amplification of a turbulence triggering mechanism. Intense
discussions have curtailed in the past about the role of resolutions to
capture transitional phenomena 31, 6, 34, 14 while the physics, non-linear
dynamics of a transitional flow, and its mechanobiological significance, if
any, must be viewed with equal attention. A closer look at figure 2b and the
animations shows disturbance in the flow jet that emanates from the nozzle
throat discussed above. These disturbances in the jet were also seen in the
work of Fehn et al.,  7 despite the flow remained laminar in their simulations
at this $Re_{th}$. Upon slight increase in $Re_{th}$ to $\sim 2400$ the jet
did breakdown in their simulations to transition the flow. It may be inferred
that at $Re_{th}=2000$, the flow is on the verge of breakdown, which,
depending on the numerics and inflow conditions used may or may not
quantitatively breakdown. A complementary question is the circumstances that
perturb the flow in the first place to trigger turbulence. It may be
hypothesized that the perfect symmetry of the mesh in higher order methods
suppresses the onset of transition but a conclusive statement on that cannot
be made. Despite the fact that symmetry was ensured in the LB setup, there
could have been artifacts from the boundaries that manifested as instabilities
in the flow thereby triggering turbulence. White & Chong,  32 demonstrated the
inferior rotational invariance of the $D3Q19$ lattice and found $D3Q27$
superior whereas Dellar,  5 advocated that multi time relaxation (MRT) model
of the LBM recovers the rotational invariance. Peng et al.,  21 also recently
found both these types of lattice to yield accurate turbulent flow statistics.
None of the study to the author’s knowledge has investigated the role of
lattice types in combination with the higher order wall boundary condition
that was employed in this study 3 and also none has used these high
resolutions. It may be inferred that a combination of MRT, higher order wall
representation, and high resolutions must have overcome the minor deficiencies
of the $D3Q19$ lattice but a detailed investigation of that is left for future
efforts. It is also emphasized that the LBM is inherently a transient scheme,
which might explain a closer match to the experiments. In a setup like the FDA
nozzle it is clear that the sudden expansion resulted in adverse gradients of
pressure, which, at a sufficiently high $Re_{th}$ departed the flow from its
laminar regime.

It is finally remarked that the FDA nozzle is essentially a device to study
blood flow and the non-Newtonian affects due to blood cells should be
accounted for in future 25. In LBM simulations non-Newtonian models that
account for shear thinning behavior of fluid can be easily incorporated, and
such a model is expected to _delay_ the onset of flow transition. On the other
hand, however, Tupin et al.,  30 have recently demonstrated a unique inverse
energy cascade in blood flow and have found the turbulence of non-Kolmogorov
type. The studies of transitional bioflows in future may study explicit
transport of red blood cells, which requires incorporation of, and coupling
with, other methods with LBM 29. The present work has prospects to serve as
reference for future extensions and comparisons.

##  5 Conclusions and implications

  1. 1.

The LBM is an adequate numerical technique for the DNS of biomedical flows in
transitional regime, and can reproduce flow characteristics without much
parameter tuning even at higher Reynolds numbers. The FDA benchmark, for
transitional and turbulent flow regimes is suitable but not fully vigorous for
a detailed quantitative comparison of CFD codes. The definition of the
benchmark should have a reliable and quantifiable turbulence triggering
mechanism, which may be incorporated into CFD models.

  2. 2.

The practice of _quantitatively_ and _qualitatively_ comparing experiments and
simulations for a transitional flow is not entirely appropriate. Previously
conducted experiments can only provide guidance on the setup of simulation and
help in the analysis of results thereof. A comparison can only be performed by
extensive joint efforts of experimentalists and computational researchers to
_tweak_ experimental aspects like noise at inflow according to the
simulations, or adjusting simulation aspects like initial conditions or inflow
distortions. Previously reported discrepancies between experimental and
computational results are non quantifiable and their source, while can be
conjectured, it cannot be ascertained.

  3. 3.

A transitional flow is characterized by chaotic eddies and vortices with rapid
annihilation and merger of vortices, and their interaction with the flow and
each other, which results in distortions within the jet that emanates from the
throat, continous restabilization and re-disruption of the jet at regular
intervals. This results in shifts in the jet breakdown location, which are
more pronounced at lower Reynolds numbers close to the $Re_{crit}$.

  4. 4.

Questions like _jet breakdown location_ have received major attention in the
literature. The $Re_{crit}$, however, should be given equal attention and
causes of discrepancies in its identification should be scrutinized to have
better understanding of mechanobiological aspects like, for example, hemolysis
or endothelial dysfunction.

  5. 5.

The LBM being a second order method requires relatively higher mesh density
compared to other numerical methods but it can compute flows accurately, and
in particular it can predict the onset of transition accurately in a
relatively lesser time. On the other hand the method scales impeccably on
massively parallel computer architectures thereby allowing simulations on
complex geometries at any scale. Whether these facets of the LBM are assets or
liabilities would depend on the perspective, the problem in hand and the
research question itself.

  6. 6.

LBM based DNS approach would likely become _impractical_ for turbulent flows
in complex geometries at Reynolds numbers higher than $6500$, and complex
collision models or turbulence models should be explored in future.

### acknowledgements

Compute resources on the _SuperMUC-NG_ were provided by the Leibniz
Supercomputing Center (LRZ), Munich, Germany. The author is grateful for the
able support of colleagues at LRZ.

### Conflict of Interest Statement

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

## References

  * Bergersen et al., 2019 Bergersen, A. W., Mortensen, M., & Valen-Sendstad, K. (2019).  The FDA nozzle benchmark:”in theory there is no difference between theory and practice, but in practice there is”.  International journal for numerical methods in biomedical engineering, 35(1), e3150.
  * Bernsdorf et al., 2008 Bernsdorf, J., Harrison, S. E., Smith, S. M., Lawford, P. V., & Hose, D. R. (2008).  Applying the lattice Boltzmann technique to biofluids: A novel approach to simulate blood coagulation.  Computers & Mathematics with Applications, 55(7), 1408–1414.
  * Bouzidi et al., 2001 Bouzidi, M., Firdaouss, M., & Lallemand, P. (2001).  Momentum transfer of a Boltzmann-lattice fluid with boundaries.  Physics of Fluids, 13, 3452.
  * Chabannes et al., 2017 Chabannes, V., Prud’Homme, C., Szopos, M., & Tarabay, R. (2017).  High order finite element simulations for fluid dynamics validated by experimental data from the FDAbenchmark nozzle model.  arXiv preprint arXiv:1701.02179.
  * Dellar, 2014 Dellar, P. J. (2014).  Lattice Boltzmann algorithms without cubic defects in Galilean invariance on standard lattices.  Journal of Computational Physics, 259, 270–283.
  * Delorme et al., 2013 Delorme, Y. T., Anupindi, K., & Frankel, S. H. (2013).  Large eddy simulation of FDA’s idealized medical device.  Cardiovascular engineering and technology, 4(4), 392–407.
  * Fehn et al., 2019 Fehn, N., Wall, W. A., & Kronbichler, M. (2019).  Modern discontinuous Galerkin methods for the simulation of transitional and turbulent flows in biomedical engineering: A comprehensive les study of the FDA benchmark nozzle model.  International Journal for Numerical Methods in Biomedical Engineering, 01(1), e3228.
  * Hariharan et al., 2011 Hariharan, P., Giarra, M., Reddy, V., Day, S. W., Manning, K. B., Deutsch, S., Stewart, S. F., Myers, M. R., Berman, M. R., Burgreen, G. W., et al. (2011).  Multilaboratory particle image velocimetry analysis of the FDA benchmark nozzle model to support validation of computational fluid dynamics simulations.  Journal of biomechanical engineering, 133(4), 041002.
  * Harlacher et al., 2012 Harlacher, D. F., 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 et al., 2014 Hasert, M., Masilamani, K., Zimny, S., Klimach, H., Qi, J., Bernsdorf, J., & Roller, S. (2014).  Complex fluid simulations with the parallel tree-based lattice Boltzmann solver Musubi.  Journal of Computational Science, 5(5), 784–794.
  * Helgeland et al., 2014 Helgeland, A., Mardal, K.-A., Haughton, V., & Anders Pettersson Reif, B. (2014).  Numerical simulations of the pulsating flow of cerebrospinal fluid flow in the cervical spinal canal of a chiari patient.  Journal of Biomechanics.
  * Jain, 2019 Jain, K. (2019).  Transition to turbulence in an oscillatory flow through stenosis.  Biomechanics and Modeling in Mechanobiology, 1, 1–19.  PMID:31359287.
  * Jain et al., 2017 Jain, K., Ringstad, G., Eide, P.-K., & Mardal, K.-A. (2017).  Direct numerical simulation of transitional hydrodynamics of the cerebrospinal fluid in chiari I malformation: The role of cranio-vertebral junction.  International journal for numerical methods in biomedical engineering, 33(9), e02853.
  * Jain et al., 2016 Jain, K., Roller, S., & Mardal, K.-A. (2016).  Transitional flow in intracranial aneurysms–a space and time refinement study below the Kolmogorov scales using lattice Boltzmann method.  Computers & Fluids, 127, 36–46.
  * Janiga, 2014 Janiga, G. (2014).  Large eddy simulation of the FDA benchmark nozzle for a Reynolds number of 6500.  Computers in biology and medicine, 47, 113–119.
  * Junk & Yang, 2011 Junk, M. & Yang, Z. (2011).  Asymptotic analysis of lattice Boltzmann outflow treatments.  Communications in Computational Physics, 9(5), 1117–1127.
  * Klimach et al., 2014 Klimach, H., Jain, K., & Roller, S. (2014).  End-to-end parallel simulations with apes.  In Parallel Computing: Accelerating Computational Science and Engineering (CSE), volume 25 (pp. 703–711).
  * Marié et al., 2009 Marié, S., Ricot, D., & Sagaut, P. (2009).  Comparison between lattice Boltzmann method and Navier–Stokes high order schemes for computational aeroacoustics.  Journal of Computational Physics, 228(4), 1056–1070.
  * Nicoud et al., 2018 Nicoud, F., Chnafa, C., Siguenza, J., Zmijanovic, V., & Mendez, S. (2018).  Large-eddy simulation of turbulence in cardiovascular flows.  (pp. 147–167).
  * Passerini et al., 2013 Passerini, T., Quaini, A., Villa, U., Veneziani, A., & Canic, S. (2013).  Validation of an open source framework for the simulation of blood flow in rigid and deformable vessels.  International journal for numerical methods in biomedical engineering, 29(11), 1192–1213.
  * Peng et al., 2018 Peng, C., Geneva, N., Guo, Z., & Wang, L.-P. (2018).  Direct numerical simulation of turbulent pipe flow using the lattice Boltzmann method.  Journal of Computational Physics, 357, 16–42.
  * Pope, 2000 Pope, S. B. (2000).  Turbulent flows.  Cambridge university press.
  * Qi et al., 2016 Qi, J., Jain, K., Klimach, H., & Roller, S. (2016).  Performance evaluation of the LBM solver Musubi on various HPC architectures.  In Advances in Parallel Computing: On the Road to Exascale, volume 27 of Advances in Parallel Computing (pp. 807–816).: IOS Press.
  * Roller et al., 2012 Roller, S., Bernsdorf, J., Klimach, H., Hasert, M., Harlacher, D., Cakircali, M., Zimny, S., Masilamani, K., Didinger, L., & Zudrop, J. (2012).  An adaptable simulation framework based on a linearized octree.  In High Performance Computing on Vector Systems 2011 (pp. 93–105).
  * Saqr et al., 2019 Saqr, K. M., Mansour, O., Tupin, S., Hassan, T., & Ohta, M. (2019).  Evidence for non-newtonian behavior of intracranial blood flow from doppler ultrasonography measurements.  Medical & biological engineering & computing, 57(5), 1029–1036.
  * Stewart et al., 2012 Stewart, S. F., Paterson, E. G., Burgreen, G. W., Hariharan, P., Giarra, M., Reddy, V., Day, S. W., Manning, K. B., Deutsch, S., Berman, M. R., et al. (2012).  Assessment of cfd performance in simulations of an idealized medical device: results of fda’s first computational interlaboratory study.  Cardiovascular Engineering and Technology, 3(2), 139–160.
  * Succi, 2008 Succi, S. (2008).  Lattice Boltzmann across scales: from turbulence to DNA translocation.  The European Physical Journal B, 64(3-4), 471–479.
  * Succi et al., 1991 Succi, S., Benzi, R., & Higuera, F. (1991).  The lattice Boltzmann equation: a new tool for computational fluid-dynamics.  Physica D: Nonlinear Phenomena, 47(1), 219–230.
  * Sun & Munn, 2008 Sun, C. & Munn, L. L. (2008).  Lattice-Boltzmann simulation of blood flow in digitized vessel networks.  Computers & Mathematics with Applications, 55(7), 1594–1600.
  * Tupin et al., 2020 Tupin, S., Saqr, K. M., Rashad, S., Niizuma, K., Ohta, M., & Tominaga, T. (2020).  Non-kolmogorov turbulence and inverse energy cascade in intracranial aneurysm: Near-wall scales suggest mechanobiological relevance.  arXiv preprint arXiv:2001.08234.
  * Ventikos, 2014 Ventikos, Y. (2014).  Resolving the issue of resolution.  American Journal of Neuroradiology, 35(3), 544–545.
  * White & Chong, 2011 White, A. T. & Chong, C. K. (2011).  Rotational invariance in the three-dimensional lattice Boltzmann method is dependent on the choice of lattice.  Journal of Computational Physics, 230(16), 6367–6378.
  * Zhang et al., 2008 Zhang, J., Johnson, P. C., & Popel, A. S. (2008).  Red blood cell aggregation and dissociation in shear flows simulated by lattice Boltzmann method.  Journal of biomechanics, 41(1), 47–55.
  * Zmijanovic et al., 2017 Zmijanovic, V., Mendez, S., Moureau, V., & Nicoud, F. (2017).  About the numerical robustness of biomedical benchmark cases: interlaboratory FDA’s idealized medical device.  International journal for numerical methods in biomedical engineering, 33(1), e02789.


Case c): $p\in(\varrho,1)$. (This corresponds to part (c)). We can set
$\mathcal{u}=\mathcal{u}_{3}:=\frac{1-p}{2\theta-1}$ to attain the optimal
rate. Taking $\nu=1$, (57) can be written as

$y_{k+1}\leqslant\left[1+\frac{\ell_{1}\alpha^{\tau}}{(k+\gamma)^{\tau
p}}-\frac{2\theta\ell_{2}\alpha}{k+\gamma}\right]y_{k}+\left[\frac{(2\theta-1)\ell_{2}\alpha}{(k+\gamma)^{1+\mathcal{u}_{3}}}+\frac{\ell_{3}\alpha^{\tau}}{(k+\gamma)^{\tau
p}}\right].$

If $\alpha\geqslant\frac{2\mathcal{u}_{3}}{\theta\ell_{2}}$ and
$\gamma\geqslant\max\\{(\frac{\alpha^{\tau-1}\ell_{1}}{\theta\ell_{2}})^{\frac{1}{\tau
p-1}},(\frac{\alpha^{\tau-1}\ell_{3}}{\ell_{2}})^{\frac{1}{\tau
p-\mathcal{u}_{3}-1}},\alpha\theta\ell_{2}\\}$, it follows

$\frac{\ell_{1}\alpha^{\tau}}{(k+\gamma)^{\tau
p}}-\frac{2\theta\ell_{2}\alpha}{k+\gamma}\leqslant-\frac{\theta\ell_{2}\alpha}{k+\gamma},\quad\frac{\ell_{3}\alpha^{\tau}}{(k+\gamma)^{\tau
p}}\leqslant\frac{\ell_{2}\alpha}{(k+\gamma)^{1+\mathcal{u}_{3}}},\quad\frac{\theta\ell_{2}\alpha}{k+\gamma}\leqslant
1,$

$\alpha\theta\ell_{2}\geqslant 2\mathcal{u}_{3}$ for all $k\geqslant 0$.
Consequently, we can infer

$y_{k+1}\leqslant\left[1-\frac{\theta\ell_{2}\alpha}{k+\gamma}\right]y_{k}+\frac{2\theta\ell_{2}\alpha}{(k+\gamma)^{1+\mathcal{u}_{3}}}.$

Thus, applying Theorem 24 (b) with $c=\theta\ell_{2}\alpha$,
$d=2\theta\ell_{2}\alpha$, $q=\mathcal{u}_{3}$ and utilizing
$c-q\geqslant\frac{c}{2}$, we obtain

$y_{k+1}\leqslant
4(k+1+\gamma)^{-\mathcal{u}_{3}}+y_{0}[\gamma^{-1}(k+1+\gamma)]^{-\theta\ell_{2}\alpha}.$

Sub-Case 2-2: $p=1$. Here, we introduce the auxiliary sequence
$u_{k}:=\log(k+\gamma)^{-\frac{1}{2\theta-1}}$. Mimicking the earlier cases,
there exists $\gamma_{0}>0$, which solely depends on
$\alpha,\ell_{1},\ell_{2},\ell_{3},\tau$, such that for any
$\gamma\geqslant\gamma_{0}$, we have
$\gamma\log(\gamma)\geqslant\alpha\theta\ell_{2}$ and (57) can be simplified
to

$y_{k+1}\leqslant\left[1-\frac{\alpha\theta\ell_{2}}{(k+\gamma)\log(k+\gamma)}\right]y_{k}+\frac{2\alpha\theta\ell_{2}}{(k+\gamma)\log(k+\gamma)^{\frac{2\theta}{2\theta-1}}}.$

(We will omit an explicit construction of such $\gamma_{0}$ here). Let us now
define $s(x):=\frac{1}{\alpha\theta\ell_{2}}x\log(x)$,
$t(x):=\frac{1}{2\alpha\theta\ell_{2}}x\log(x)^{\frac{2\theta}{2\theta-1}}$,
$r(x):=2\log(x)^{-\frac{1}{2\theta-1}}$, and $b_{k}=k+\gamma$. Then, it holds
that

$\displaystyle u(x)$ $\displaystyle:=s^{\prime}(x)-r(x)t^{\prime}(x)$
$\displaystyle=\frac{1}{\alpha\theta\ell_{2}}(1+\log(x))-\frac{1}{\alpha\theta\ell_{2}}\cdot\frac{1}{\log(x)^{\frac{1}{2\theta-1}}}\left(\frac{2\theta}{2\theta-1}\log(x)^{\frac{1}{2\theta-1}}+\log(x)^{\frac{2\theta}{2\theta-1}}\right)$
$\displaystyle=-\frac{1}{\alpha\theta(2\theta-1)\ell_{2}}.$

Thus, if $\alpha\geqslant\frac{2}{\theta(2\theta-1)\ell_{2}}$, it follows
$(b_{k+1}-b_{k})u(b_{k})\geqslant-\frac{1}{2}$. This the rate function $r$ is
clearly convex on $I=[\gamma,\infty)$, all assumptions in Theorem 1 (b) are
satisfied with $\lambda_{k}\equiv 2$. Hence, applying Theorem 1 (b), we can
infer

$\displaystyle y_{k+1}\leqslant\leavevmode\nobreak$ $\displaystyle
4\log(k+1+\gamma)^{-\frac{1}{2\theta-1}}+y_{0}\exp\Big{(}-{\sum}_{i=0}^{k}\frac{\alpha\theta\ell_{2}}{(i+\gamma)\log(i+\gamma)}\Big{)}$
$\displaystyle\leqslant\leavevmode\nobreak$ $\displaystyle
4\log(k+1+\gamma)^{-\frac{1}{2\theta-1}}+y_{0}\Big{(}\frac{\log(k+1+\gamma)}{\log(\gamma)}\Big{)}^{-\alpha\theta\ell_{2}},$

where we utilize the integral test
$\sum_{i=0}^{k}\frac{1}{(k+\gamma)\log(k+\gamma)}\geqslant\int_{0}^{k+1}\frac{1}{(x+\gamma)\log(x+\gamma)}\,\mathrm{d}x=\log(\log(k+1+\gamma))-\log(\log(\gamma))$
in the last line.

###  E.2 Results for Polynomial Step Sizes

Proof of Theorem 20 The condition $\gamma\geqslant({\sf
L}\alpha)^{{1}/{\varrho_{1}}}$ ensures $\alpha_{k}\leqslant\frac{1}{{\sf L}}$
for all $k$ and hence, the descent property (15) is applicable. The first
convergence rate (44) in Theorem 20 can be obtained by applying Lemma 19 (b)
to (15) with $\ell_{1}=\frac{{\sf A}{\sf L}}{2}$, $\ell_{2}=\mu$,
$\ell_{3}=\frac{{\sf L}\sigma^{2}}{2}$, and $\tau=2$.

To establish (45), we substitute
$\alpha=(\frac{2\omega_{1}}{\zeta_{1}})^{\varrho_{1}}$ in
$\gamma=\max\\{(\frac{{\sf A}{\sf
L}}{2\zeta_{1}})^{\frac{1}{\omega_{1}}}\alpha^{\frac{1}{\varrho_{1}}},\zeta_{1}\alpha^{\frac{1}{\varrho_{1}}},({\sf
L}\alpha)^{\frac{1}{\varrho_{1}}}\\}$ and (44). Then, using
$\frac{1}{\varrho_{1}}=2-\frac{1}{2\theta}$, $\frac{1}{\omega_{1}}=4\theta-1$,
and $\zeta_{1}=\theta\mu^{\frac{1}{2\theta}}(\frac{{\sf
L}\sigma^{2}}{2})^{1-\frac{1}{2\theta}}$, we have

$\displaystyle\gamma=2\omega_{1}\max\Big{\\{}\Big{(}\frac{{\sf A}{\sf
L}}{2\zeta_{1}}\Big{)}^{\frac{1}{\omega_{1}}}\frac{1}{\zeta_{1}},1,\frac{{\sf
L}^{{1}/{\varrho_{1}}}}{\zeta_{1}}\Big{\\}}$
$\displaystyle=2\omega_{1}\max\Big{\\{}\frac{({\sf A}{\sf
L})^{4\theta-1}}{2^{4\theta-1}\zeta_{1}^{4\theta}},1,\frac{2^{1-\frac{1}{2\theta}}{\sf
L}}{\theta\mu^{\frac{1}{2\theta}}\sigma^{2-\frac{1}{\theta}}}\Big{\\}}$
$\displaystyle=2\omega_{1}\max\\{1,C_{1},C_{2}\\}$

and $(\frac{{\sf
L}\sigma^{2}}{2\mu})^{\frac{1}{2\theta}}\alpha^{\frac{1}{2\theta}}=(\frac{{\sf
L}\sigma^{2}}{2\mu})^{\frac{1}{2\theta}}(\frac{2\omega_{1}}{\zeta_{1}})^{\omega_{1}}=(\frac{\omega_{1}}{\theta})^{\omega_{1}}(\frac{{\sf
L}\sigma^{2}}{\mu^{2}})^{\omega_{1}}$.

Proof of Theorem 21 The inequality $\gamma\geqslant(2{\sf
L}\alpha)^{\frac{1}{\varrho_{2}}}$ ensures that
$\alpha_{k}\leqslant\frac{1}{2{\sf L}}$, which means that (16) is applicable.
The first inequality (46) can be obtained by applying Lemma 19 to (16) with
$\ell_{1}=\frac{{\sf A}{\sf L}^{2}}{2N}$, $\ell_{2}=\frac{\mu}{2}$, and
$\ell_{3}=\frac{{\sf L}^{2}\sigma^{2}}{2N}$. To establish (47), we substitute

$\displaystyle\gamma$ $\displaystyle=\max\Big{\\{}(\tfrac{{\sf A}{\sf
L}^{2}}{2N\zeta_{2}})^{\frac{1}{\omega_{2}}}\alpha^{\frac{1}{\varrho_{2}}},\zeta_{2}\alpha^{\frac{1}{\varrho_{2}}},(2{\sf
L}\alpha)^{\frac{1}{\varrho_{2}}}\Big{\\}}=2\omega_{2}(1+\log(N))\max\cdot\Big{\\{}\tfrac{C_{3}}{\sqrt{N}},C_{4}N^{1-\frac{1}{2\theta}},1\Big{\\}}$
$\displaystyle\alpha$
$\displaystyle=\Big{[}\tfrac{2\omega_{2}(1+\log(N))}{\zeta_{2}}\Big{]}^{\varrho_{2}}\quad\text{and}\quad
K\geqslant 2\gamma$

into (46), then it holds that

$\displaystyle y_{K}$ $\displaystyle\leqslant
4\Big{(}\frac{\ell_{3}}{\ell_{2}}\Big{)}^{\frac{1}{2\theta}}\Big{(}\frac{2{\omega_{2}}}{\zeta}\Big{)}^{\omega_{2}}\Big{[}\frac{1+\log(N)}{K+\gamma}\Big{]}^{{\omega_{2}}}+y_{0}[\gamma^{-1}(K+\gamma)]^{-2{\omega_{2}}(1+\log(N))}$
$\displaystyle\leqslant
4\Big{(}\frac{4{\omega_{2}}}{\theta}\Big{)}^{\omega_{2}}\Big{(}\frac{{\sf
L}\sigma}{\mu^{{3}/{2}}}\Big{)}^{\omega_{2}}\Big{[}\frac{1+\log(N)}{\sqrt{N}(K+\gamma)}\Big{]}^{{\omega_{2}}}+y_{0}\Big{[}{\frac{\gamma}{K+\gamma}}\Big{]}^{2{\omega_{2}}}[\gamma^{-1}(K+\gamma)]^{\frac{-2{\omega_{2}}\log(N)}{\log[(K+\gamma)/\gamma]}}$
$\displaystyle=4\Big{(}\frac{4\omega_{2}}{\theta}\Big{)}^{\omega_{2}}\Big{(}\frac{{\sf
L}\sigma}{\mu^{{3}/{2}}}\Big{)}^{\omega_{2}}\Big{[}\frac{1+\log(N)}{\sqrt{N}(K+\gamma)}\Big{]}^{\omega_{2}}+y_{0}\Big{[}{\frac{2\omega_{2}(1+\log(N))\max\\{\frac{C_{3}}{\sqrt{N}},C_{4}N^{1-\frac{1}{2\theta}},1\\}}{N(K+\gamma)}}\Big{]}^{2\omega_{2}}.$

This completes the proof.

## References

  * Abadi et al. (2015) M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, et al.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015.
  * Agarwal et al. (2021) A. Agarwal, S. M. Kakade, J. D. Lee, and G. Mahajan.  On the theory of policy gradient methods: Optimality, approximation, and distribution shift.  _The Journal of Machine Learning Research_ , 22(1):4431–4506, 2021.
  * Ahn et al. (2020) K. Ahn, C. Yun, and S. Sra.  SGD with shuffling: optimal rates without component convexity and large epoch requirements.  _Advances in Neural Information Processing Systems_ , 33:17526–17535, 2020.
  * Bolte et al. (2014) J. Bolte, S. Sabach, and M. Teboulle.  Proximal alternating linearized minimization for nonconvex and nonsmooth problems.  _Mathematical Programming_ , 146(1-2, Ser. A):459–494, 2014.
  * Bottou et al. (2018) L. Bottou, F. E. Curtis, and J. Nocedal.  Optimization methods for large-scale machine learning.  _SIAM Review_ , 60(2):223–311, 2018.
  * Chung (1954) K. L. Chung.  On a stochastic approximation method.  _The Annals of Mathematical Statistics_ , pages 463–483, 1954.
  * Derman (1956) C. Derman.  An application of chung’s lemma to the Kiefer-Wolfowitz stochastic approximation procedure.  _The Annals of Mathematical Statistics_ , pages 532–536, 1956.
  * Fabian (1968) V. Fabian.  On asymptotic normality in stochastic approximation.  _The Annals of Mathematical Statistics_ , pages 1327–1332, 1968.
  * Fatkhullin et al. (2022) I. Fatkhullin, J. Etesami, N. He, and N. Kiyavash.  Sharp analysis of stochastic optimization under global Kurdyka-łojasiewicz inequality.  In _Advances in Neural Information Processing Systems_ , volume 35, pages 15836–15848, 2022.
  * Fontaine et al. (2021) X. Fontaine, V. De Bortoli, and A. Durmus.  Convergence rates and approximation results for SGD and its continuous-time counterpart.  In _Conference on Learning Theory_ , pages 1965–2058. PMLR, 2021.
  * Ge et al. (2019) R. Ge, S. M. Kakade, R. Kidambi, and P. Netrapalli.  The step decay schedule: A near optimal, geometrically decaying learning rate procedure for least squares.  _Advances in neural information processing systems_ , 32, 2019.
  * Goffin (1977) J.-L. Goffin.  On convergence rates of subgradient optimization methods.  _Mathematical programming_ , 13:329–347, 1977.
  * Gower et al. (2021) R. Gower, O. Sebbouh, and N. Loizou.  SGD for structured nonconvex functions: Learning rates, minibatching and interpolation.  In _International Conference on Artificial Intelligence and Statistics_ , pages 1315–1323. PMLR, 2021.
  * Gower et al. (2019) R. M. Gower, N. Loizou, X. Qian, A. Sailanbayev, E. Shulgin, and P. Richtárik.  SGD: General analysis and improved rates.  In _International conference on machine learning_ , pages 5200–5209. PMLR, 2019.
  * Gürbüzbalaban et al. (2019) M. Gürbüzbalaban, A. Ozdaglar, and P. A. Parrilo.  Convergence rate of incremental gradient and incremental Newton methods.  _SIAM Journal on Optimization_ , 29(4):2542–2565, 2019.
  * Gürbüzbalaban et al. (2021) M. Gürbüzbalaban, A. Ozdaglar, and P. Parrilo.  Why random reshuffling beats stochastic gradient descent.  _Mathematical Programming_ , 186(1-2):49–84, 2021.
  * HaoChen and Sra (2019) J. Z. HaoChen and S. Sra.  Random shuffling beats SGD after finite epochs.  In _International Conference on Machine Learning_ , pages 2624–2633, 2019.
  * Huang et al. (2021) X. Huang, K. Yuan, X. Mao, and W. Yin.  An improved analysis and rates for variance reduction under without-replacement sampling orders.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Karimi et al. (2016) H. Karimi, J. Nutini, and M. Schmidt.  Linear convergence of gradient and proximal-gradient methods under the Polyak-Łojasiewicz condition.  In _Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2016, Riva del Garda, Italy, September 19-23, 2016, Proceedings, Part I 16_ , pages 795–811. Springer, 2016.
  * Li et al. (2021) X. Li, Z. Zhuang, and F. Orabona.  A second look at exponential and cosine step sizes: Simplicity, adaptivity, and performance.  In _International Conference on Machine Learning_ , pages 6553–6564. PMLR, 2021.
  * Li et al. (2023) X. Li, A. Milzarek, and J. Qiu.  Convergence of random reshuffling under the Kurdyka–Łojasiewicz inequality.  _SIAM Journal on Optimization_ , 33(2):1092–1120, 2023.
  * Łojasiewicz (1959) S. Łojasiewicz.  Sur le problème de la division.  _Studia Mathematica_ , 18:87–136, 1959.
  * Łojasiewicz (1963) S. Łojasiewicz.  Une propriété topologique des sous-ensembles analytiques réels.  _Les équations aux dérivées partielles_ , 117:87–89, 1963.
  * Loshchilov and Hutter (2016) I. Loshchilov and F. Hutter.  SGDR: Stochastic gradient descent with warm restarts.  In _International Conference on Learning Representations_ , 2016.
  * Mei et al. (2020) J. Mei, C. Xiao, C. Szepesvari, and D. Schuurmans.  On the global convergence rates of softmax policy gradient methods.  In _International Conference on Machine Learning_ , pages 6820–6829. PMLR, 2020.
  * (26) K. Mishchenko, A. Khaled, and P. Richtarik.  Random reshuffling: Simple analysis with vast improvements.  In _Advances in Neural Information Processing Systems_ , pages 17309–17320.
  * Moulines and Bach (2011) E. Moulines and F. Bach.  Non-asymptotic analysis of stochastic approximation algorithms for machine learning.  In _Advances in Neural Information Processing Systems_ , volume 24. Curran Associates, Inc., 2011.
  * Needell et al. (2014) D. Needell, R. Ward, and N. Srebro.  Stochastic gradient descent, weighted sampling, and the randomized Kaczmarz algorithm.  _Advances in neural information processing systems_ , 27, 2014.
  * Nemirovski et al. (2009) A. Nemirovski, A. Juditsky, G. Lan, and A. Shapiro.  Robust stochastic approximation approach to stochastic programming.  _SIAM Journal on optimization_ , 19(4):1574–1609, 2009.
  * Nguyen et al. (2021) L. M. Nguyen, Q. Tran-Dinh, D. T. Phan, P. H. Nguyen, and M. van Dijk.  A unified convergence analysis for shuffling-type gradient methods.  _Journal of Machine Learning Research_ , 2021.
  * Paszke et al. (2019) A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  _Advances in neural information processing systems_ , 32, 2019.
  * Polyak (1963) B. T. Polyak.  Gradient methods for the minimisation of functionals.  _USSR Computational Mathematics and Mathematical Physics_ , 3(4):864–878, 1963.
  * Polyak (1987) B. T. Polyak.  _Introduction to optimization_.  Translations Series in Mathematics and Engineering. Optimization Software, Inc., Publications Division, New York, 1987.  Translated from the Russian, With a foreword by Dimitri P. Bertsekas.
  * Robbins and Monro (1951) H. Robbins and S. Monro.  A stochastic approximation method.  _The annals of mathematical statistics_ , pages 400–407, 1951.
  * Sacks (1958) J. Sacks.  Asymptotic distribution of stochastic approximation procedures.  _The Annals of Mathematical Statistics_ , 29(2):373–405, 1958.
  * Shamir (2016) O. Shamir.  Without-replacement sampling for stochastic gradient methods.  In _Advances in Neural Information Processing Systems_ , pages 46–54, 2016.
  * Shamir and Zhang (2013) O. Shamir and T. Zhang.  Stochastic gradient descent for non-smooth optimization: Convergence results and optimal averaging schemes.  In _International conference on machine learning_ , pages 71–79. PMLR, 2013.
  * Sohl-Dickstein et al. (2014) J. Sohl-Dickstein, B. Poole, and S. Ganguli.  Fast large-scale optimization by unifying stochastic gradient and quasi-newton methods.  In _International Conference on Machine Learning_ , pages 604–612. PMLR, 2014.
  * Touvron et al. (2023) H. Touvron, L. Martin, K. Stone, P. Albert, A. Almahairi, Y. Babaei, N. Bashlykov, S. Batra, P. Bhargava, S. Bhosale, et al.  Llama 2: Open foundation and fine-tuned chat models.  _arXiv preprint arXiv:2307.09288_ , 2023.
  * Tran et al. (2021) T. H. Tran, L. M. Nguyen, and Q. Tran-Dinh.  Smg: A shuffling gradient-based method with momentum.  In _International Conference on Machine Learning_ , pages 10379–10389. PMLR, 2021.
  * Wang et al. (2021) X. Wang, S. Magnússon, and M. Johansson.  On the convergence of step decay step-size for stochastic optimization.  _Advances in Neural Information Processing Systems_ , 34:14226–14238, 2021.
  * Wang and Zou (2022) Y. Wang and S. Zou.  Policy gradient method for robust reinforcement learning.  In _International Conference on Machine Learning_ , pages 23484–23526. PMLR, 2022.