# Simulation of the FDA nozzle benchmark: A lattice Boltzmann study

Feng Huang, Romain Noël, Philipp Berg, Seyed Ali Hosseini

► To cite this version:

Feng Huang, Romain Noël, Philipp Berg, Seyed Ali Hosseini. Simulation of the FDA nozzle benchmark: A lattice Boltzmann study. Computer Methods and Programs in Biomedicine, 2022, 221, pp.106863. 10.1016/j.cmpb.2022.106863. hal-03881693

HAL Id: hal-03881693

https://inria.hal.science/hal-03881693

Submitted on 6 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Simulation of the FDA Nozzle Benchmark: A Lattice Boltzmann Study

Feng Huang<sup>a</sup>, Romain Noël<sup>c</sup>, Philipp Berg<sup>a,d</sup>, Seyed Ali Hosseini<sup>a,b,*</sup>

<sup>a</sup>Laboratory of Fluid Dynamics and Technical Flows, University of Magdeburg “Otto von Guericke”, D-39106 Magdeburg, Germany

<sup>b</sup>Department of Mechanical and Process Engineering, ETH Zürich, 8092 Zürich, Switzerland

<sup>c</sup>Univ. Gustave Eiffel, Inria, Cosys/SII, I4S, F-44344 Bouguenais, France

<sup>d</sup>Research Campus STIMULATE, University of Magdeburg “Otto von Guericke”, D-39106, Magdeburg, Germany

**Abstract**

*Background and objective:* Contrary to flows in small intracranial vessels, many blood flow configurations such as those found in aortic vessels and aneurysms involve larger Reynolds numbers and, therefore, transitional or turbulent conditions. Dealing with such systems require both robust and efficient numerical methods.

*Methods:* We assess here the performance of a lattice Boltzmann solver with full Hermite expansion of the equilibrium and central Hermite moments collision operator at higher Reynolds numbers, especially for under-resolved simulations. To that end the food and drug administration's benchmark nozzle is considered at three different Reynolds numbers covering all regimes: 1) laminar at a Reynolds number of 500, 2) transitional at a Reynolds number of 3500, and 3) low-level turbulence at a Reynolds number of 6500.

*Results:* The lattice Boltzmann results are compared with previously published inter-laboratory experimental data obtained by particle image velocimetry. Our results show good agreement with the experimental measurements throughout the nozzle, demonstrating the good performance of the solver even in under-resolved simulations.

*Conclusion:* In this manner, fast but sufficiently accurate numerical predictions can be achieved for flow configurations of practical interest regarding medical applications.

*Keywords:* Lattice Boltzmann method; FDA; Nozzle; Validation; Central Hermite multiple relaxation time.

*2010 MSC:* 00-01, 99-00

**1. Introduction**

Computational Fluid Dynamics (CFD) is known to be a demanding field in terms of computational resources and is sensitive to the numerical scheme employed. Traditional numerical methods used in CFD are computationally heavy and rely on relatively fastidious discrete operations, especially when the geometry is complex. The Lattice Boltzmann Method (LBM) is an alternative numerical technique that gained increasing attention, particularly so over the last two decades. Its popularity is due to the simplicity of the involved numerical operations, their locality – ensuring excellent portability to parallel computing architectures – and the relative ease at which complex geometries can be treated [1]. These strengths combined with its ability to model mechano-biological phenomena [2, 3], make the LBM a method of choice for the bio-medical community [4, 5, 6, 7, 8, 9, 10, 11, 12].

Nevertheless, computational costs are still often prohibitive in high Reynolds number flows for the classical

LBM, for instance in aortic aneurysms. That is mainly due to the fact that the standard LBM using a collision operator with a single parameter rate relaxation toward a second order expanded Equilibrium Distribution Function (EDF), suffers from shortcomings such as limited stability domain and lack of Galilean invariance [13, 14, 15].

Thus, many modifications have been proposed to extend the accuracy and stability of the LBM, most notably, through more advanced collision operators, e.g. the entropic formulation [16, 17], the Multiple Relaxation Times (MRT) [18], centered or cascaded [19], the Cumulant method [20], or the regularized LBM [21]. Furthermore, the inconsistency in the shear components can be readily overcome by using higher order expansions in the EDF [14]. The development of more consistent and robust schemes has also opened the door to the so-called under-resolved direct numerical simulations [15, 22]. A number of studies have shown that such simulations can match Large-Eddy Simulation (LES) with explicit sub-grid scale closures in terms of accuracy in many configurations, as demonstrated most recently in [23] for homogeneous isotropic turbulence. In the present work, we will further explore that avenue in the context of biological flows, through a well-established

<footnote>*Corresponding author
Email address: seyed.hosseini@ovgu.de (Seyed Ali Hosseini)</footnote>
Simulation of the FDA Nozzle Benchmark: A Lattice Boltzmann Study

Feng Huang<sup>a</sup>, Romain Noël<sup>c</sup>, Philipp Berg<sup>a,d</sup>, Seyed Ali Hosseini<sup>a,b,*</sup>

<sup>a</sup>Laboratory of Fluid Dynamics and Technical Flows, University of Magdeburg “Otto von Guericke”, D-39106 Magdeburg, Germany

<sup>b</sup>Department of Mechanical and Process Engineering, ETH Zürich, 8092 Zürich, Switzerland

<sup>c</sup>Univ. Gustave Eiffel, Inria, Cosys/SII, I4S, F-44344 Bouguenais, France

<sup>d</sup>Research Campus STIMULATE, University of Magdeburg “Otto von Guericke”, D-39106, Magdeburg, Germany

**Abstract**

*Background and objective:* Contrary to flows in small intracranial vessels, many blood flow configurations such as those found in aortic vessels and aneurysms involve larger Reynolds numbers and, therefore, transitional or turbulent conditions. Dealing with such systems require both robust and efficient numerical methods.

*Methods:* We assess here the performance of a lattice Boltzmann solver with full Hermite expansion of the equilibrium and central Hermite moments collision operator at higher Reynolds numbers, especially for under-resolved simulations. To that end the food and drug administration’s benchmark nozzle is considered at three different Reynolds numbers covering all regimes: 1) laminar at a Reynolds number of 500, 2) transitional at a Reynolds number of 3500, and 3) low-level turbulence at a Reynolds number of 6500.

*Results:* The lattice Boltzmann results are compared with previously published inter-laboratory experimental data obtained by particle image velocimetry. Our results show good agreement with the experimental measurements throughout the nozzle, demonstrating the good performance of the solver even in under-resolved simulations.

*Conclusion:* In this manner, fast but sufficiently accurate numerical predictions can be achieved for flow configurations of practical interest regarding medical applications.

*Keywords:* Lattice Boltzmann method; FDA; Nozzle; Validation; Central Hermite multiple relaxation time.

*2010 MSC:* 00-01, 99-00

**1. Introduction**

Computational Fluid Dynamics (CFD) is known to be a demanding field in terms of computational resources and is sensitive to the numerical scheme employed. Traditional numerical methods used in CFD are computationally heavy and rely on relatively fastidious discrete operations, especially when the geometry is complex. The Lattice Boltzmann Method (LBM) is an alternative numerical technique that gained increasing attention, particularly so over the last two decades. Its popularity is due to the simplicity of the involved numerical operations, their locality – ensuring excellent portability to parallel computing architectures – and the relative ease at which complex geometries can be treated [1]. These strengths combined with its ability to model mechano-biological phenomena [2, 3], make the LBM a method of choice for the bio-medical community [4, 5, 6, 7, 8, 9, 10, 11, 12].

Nevertheless, computational costs are still often prohibitive in high Reynolds number flows for the classical

LBM, for instance in aortic aneurysms. That is mainly due to the fact that the standard LBM using a collision operator with a single parameter rate relaxation toward a second order expanded Equilibrium Distribution Function (EDF), suffers from shortcomings such as limited stability domain and lack of Galilean invariance [13, 14, 15].

Thus, many modifications have been proposed to extend the accuracy and stability of the LBM, most notably, through more advanced collision operators, e.g. the entropic formulation [16, 17], the Multiple Relaxation Times (MRT) [18], centered or cascaded [19], the Cumulant method [20], or the regularized LBM [21]. Furthermore, the inconsistency in the shear components can be readily overcome by using higher order expansions in the EDF [14]. The development of more consistent and robust schemes has also opened the door to the so-called under-resolved direct numerical simulations [15, 22]. A number of studies have shown that such simulations can match Large-Eddy Simulation (LES) with explicit sub-grid scale closures in terms of accuracy in many configurations, as demonstrated most recently in [23] for homogeneous isotropic turbulence. In the present work, we will further explore that avenue in the context of biological flows, through a well-established

<footnote>*Corresponding author
Email address: seyed.hosseini@ovgu.de (Seyed Ali Hosseini)</footnote>
<table><caption>Table 1: Summary of numerical studies associated to the FDA benchmark nozzle.</caption><thead><tr><th>Study</th><th>Re= 500</th><th>Re= 2000</th><th>Re= 3500</th><th>Re= 5000</th><th>Re= 6500</th><th>Numerical method</th><th>Sub-grid model</th></tr></thead><tbody><tr><td>Sanchez Abad et al. [25]</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td></td><td>Spectral elements (Nek5000)</td><td>None</td></tr><tr><td>Fehn et al. [26]</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>High-order DG</td><td>None</td></tr><tr><td>Bergersen et al. [27]</td><td>✓</td><td></td><td>✓</td><td></td><td>✓</td><td>FEM (Oasis)</td><td>None</td></tr><tr><td>Pewowaruk et al. [28]</td><td>✓</td><td></td><td>✓</td><td></td><td>✓</td><td>FVM (Converge v2.4)</td><td>None/RANS</td></tr><tr><td>Stewart et al. [29]</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>various</td><td>various</td></tr><tr><td>Zmijanovic et al. [30]</td><td></td><td></td><td>✓</td><td></td><td></td><td>FVM (YALES2BIO)</td><td>σ-model</td></tr><tr><td>Manchester et al. [31]</td><td></td><td>✓</td><td></td><td></td><td></td><td>FVM (OpenFOAM)</td><td>WALE</td></tr><tr><td>Bhushan et al. [32]</td><td>✓</td><td></td><td>✓</td><td></td><td>✓</td><td>FVM (ANSYS Fluent)</td><td>URANS</td></tr><tr><td>Janiga [33]</td><td></td><td></td><td></td><td></td><td>✓</td><td>FVM (ANSYS Fluent)</td><td>Smagorinsky</td></tr><tr><td>Chabannes et al. [34]</td><td>✓</td><td>✓</td><td>✓</td><td></td><td></td><td>FEM (Feel++)</td><td>None</td></tr><tr><td>White et al. [35]</td><td>✓</td><td></td><td></td><td></td><td></td><td>regularized LBM (Palabos)</td><td>None</td></tr><tr><td>Jain [5]</td><td></td><td>✓</td><td>✓</td><td></td><td></td><td>raw moments MRT-LBM (Musubi)</td><td>None</td></tr><tr><td>Present work</td><td>✓</td><td></td><td>✓</td><td></td><td>✓</td><td>central Hermite MRT-LBM (ALBORZ)</td><td>None</td></tr></tbody></table>

benchmark.

The need for reliable, reproducible benchmarks confronted with experiments is urging in bio-medical applications where many variables may affect the results: temperature, in-vivo vs. ex-vivo vs. in-vitro, biological species, etc. The Food and Drug Administration (FDA), conscious of this reality, suggested a series of challenging benchmarks. The idealized medical nozzle device benchmark has been used by several laboratories to conduct Particle Image Velocimetry (PIV) measurements. Those measurements provided reference standards to test the accuracy, stability and efficacy for numerical solvers or alternative experimental measurement techniques [24]. Experimental data sets have been provided by independent laboratories at different Reynolds numbers Re= 500, 2000, 3500, 5000, and 6500, computed using the diameter of the throat (see later fig. 1). These values of Reynolds numbers are typically encountered in medical devices involving flow recirculations, contractions, or expansions. This range of Re covers very different flow regimes going from laminar (for Re = 500 and 2000) to transitional (for Re = 3500 and 5000) to low-level turbulence (for Re = 6500). The corresponding PIV experiments have been carried out completely independently by three different groups. It must be kept in mind that large inter-laboratory discrepancies have been observed, in particular regarding the jet breakdown location for transitional and turbulent flows. Conversely, a very good agreement was observed for fully laminar flows [29].

Many numerical studies of this benchmark nozzle have been performed over the years, employing various numerical schemes ranging from Finite Element Method (FEM) [27, 34] to Finite Volume Method (FVM) [30, 33], as well as a few using LBM [5, 35]. Table 1 summarizes the published numerical studies of the FDA benchmark nozzle with the associated Reynolds numbers, numerical methods and turbulence models employed. While many studies performed Direct Numerical Simulation (DNS) for the lower Reynolds numbers, most simulations in the turbulent regime were conducted using turbulence models like Unsteady Reynolds-Averaged Navier-Stokes (URANS) [28, 32] or LES [31, 33, 36].

Few studies simulated the FDA nozzle benchmark using LBM; none of those has covered all experimental conditions, more specifically the higher Reynolds numbers. White and Chong [35] performed LBM simulations at low Reynolds numbers (50 – 500) and studied the suitability of different lattice types. More recently, Jain [5] studied the FDA nozzle at Re=2000 and 3500 with a simple LBM scheme (raw moments MRT collision operator with second order EDF and D3Q19 lattice) focusing on fully-resolved simulations. Results were shown for the average axial velocity, pressure profiles, and shear stress in different planes. In the present study, we consider a wider range of Reynolds number, covering all flow regimes. The computations are done with an advanced LBM model relying on a central Hermite collision operator with a full expansion of the discrete equilibrium, implemented in the in-house LBM solver ALBORZ [37, 38, 39, 40]. The accuracy of this numerical model is assessed in under-resolved direct numerical simulations at different resolutions, akin to filtered simulations. A comprehensive comparison to experimental data covering all relevant fields is presented.

The article is organized as follows: Section 2 is dedicated to the numerical methods underlying our LBM model. Section 2.2 details geometries and conditions used. The obtained results are presented in Section 3 and discussed in Section 4. Section 5 brings concluding lights on this study.

## 2. Methods

### 2.1. Numerical model

The flow of interest here, involving a fluid with properties similar to water at room temperature and therefore a sound speed of $c_s \approx 1450$ m/s, and a maximum characteristic speed of the order of 10 m/s, falls into the low-Mach regime, with Ma below $7 \times 10^{-3}$. As such it can readily be described with a low-Mach approximation to the Navier-Stokes (NS) equations:

$$ \partial_t(\rho u_\alpha) + \partial_\beta(\rho u_\alpha u_\beta) + \partial_\beta T_{\alpha\beta} = 0, \quad (1) $$
<table><caption>Table 1: Summary of numerical studies associated to the FDA benchmark nozzle.</caption><thead><tr><th>Study</th><th>Re= 500</th><th>Re= 2000</th><th>Re= 3500</th><th>Re= 5000</th><th>Re= 6500</th><th>Numerical method</th><th>Sub-grid model</th></tr></thead><tbody><tr><td>Sanchez Abad et al. [25]</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td></td><td>Spectral elements (Nek5000)</td><td>None</td></tr><tr><td>Fehn et al. [26]</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>High-order DG</td><td>None</td></tr><tr><td>Bergersen et al. [27]</td><td>✓</td><td></td><td>✓</td><td></td><td>✓</td><td>FEM (Oasis)</td><td>None</td></tr><tr><td>Pewowaruk et al. [28]</td><td>✓</td><td></td><td>✓</td><td></td><td>✓</td><td>FVM (Converge v2.4)</td><td>None/RANS</td></tr><tr><td>Stewart et al. [29]</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>various</td><td>various</td></tr><tr><td>Zmijanovic et al. [30]</td><td></td><td></td><td>✓</td><td></td><td></td><td>FVM (YALES2BIO)</td><td>σ-model</td></tr><tr><td>Manchester et al. [31]</td><td></td><td>✓</td><td></td><td></td><td></td><td>FVM (OpenFOAM)</td><td>WALE</td></tr><tr><td>Bhushan et al. [32]</td><td>✓</td><td></td><td>✓</td><td></td><td>✓</td><td>FVM (ANSYS Fluent)</td><td>URANS</td></tr><tr><td>Janiga [33]</td><td></td><td></td><td></td><td></td><td>✓</td><td>FVM (ANSYS Fluent)</td><td>Smagorinsky</td></tr><tr><td>Chabannes et al. [34]</td><td>✓</td><td>✓</td><td>✓</td><td></td><td></td><td>FEM (Feel++)</td><td>None</td></tr><tr><td>White et al. [35]</td><td>✓</td><td></td><td></td><td></td><td></td><td>regularized LBM (Palabos)</td><td>None</td></tr><tr><td>Jain [5]</td><td></td><td>✓</td><td>✓</td><td></td><td></td><td>raw moments MRT-LBM (Musubi)</td><td>None</td></tr><tr><td>Present work</td><td>✓</td><td></td><td>✓</td><td></td><td>✓</td><td>central Hermite MRT-LBM (ALBORZ)</td><td>None</td></tr></tbody></table>

benchmark.

The need for reliable, reproducible benchmarks confronted with experiments is urging in bio-medical applications where many variables may affect the results: temperature, in-vivo vs. ex-vivo vs. in-vitro, biological species, etc. The Food and Drug Administration (FDA), conscious of this reality, suggested a series of challenging benchmarks. The idealized medical nozzle device benchmark has been used by several laboratories to conduct Particle Image Velocimetry (PIV) measurements. Those measurements provided reference standards to test the accuracy, stability and efficacy for numerical solvers or alternative experimental measurement techniques [24]. Experimental data sets have been provided by independent laboratories at different Reynolds numbers Re= 500, 2000, 3500, 5000, and 6500, computed using the diameter of the throat (see later fig. 1). These values of Reynolds numbers are typically encountered in medical devices involving flow recirculations, contractions, or expansions. This range of Re covers very different flow regimes going from laminar (for Re = 500 and 2000) to transitional (for Re = 3500 and 5000) to low-level turbulence (for Re = 6500). The corresponding PIV experiments have been carried out completely independently by three different groups. It must be kept in mind that large inter-laboratory discrepancies have been observed, in particular regarding the jet breakdown location for transitional and turbulent flows. Conversely, a very good agreement was observed for fully laminar flows [29].

Many numerical studies of this benchmark nozzle have been performed over the years, employing various numerical schemes ranging from Finite Element Method (FEM) [27, 34] to Finite Volume Method (FVM) [30, 33], as well as a few using LBM [5, 35]. Table 1 summarizes the published numerical studies of the FDA benchmark nozzle with the associated Reynolds numbers, numerical methods and turbulence models employed. While many studies performed Direct Numerical Simulation (DNS) for the lower Reynolds numbers, most simulations in the turbulent regime were conducted using turbulence models like Unsteady Reynolds-Averaged Navier-Stokes (URANS) [28, 32] or LES [31, 33, 36].

Few studies simulated the FDA nozzle benchmark using LBM; none of those has covered all experimental conditions, more specifically the higher Reynolds numbers. White and Chong [35] performed LBM simulations at low Reynolds numbers (50 – 500) and studied the suitability of different lattice types. More recently, Jain [5] studied the FDA nozzle at Re=2000 and 3500 with a simple LBM scheme (raw moments MRT collision operator with second order EDF and D3Q19 lattice) focusing on fully-resolved simulations. Results were shown for the average axial velocity, pressure profiles, and shear stress in different planes. In the present study, we consider a wider range of Reynolds number, covering all flow regimes. The computations are done with an advanced LBM model relying on a central Hermite collision operator with a full expansion of the discrete equilibrium, implemented in the in-house LBM solver ALBORZ [37, 38, 39, 40]. The accuracy of this numerical model is assessed in under-resolved direct numerical simulations at different resolutions, akin to filtered simulations. A comprehensive comparison to experimental data covering all relevant fields is presented.

The article is organized as follows: Section 2 is dedicated to the numerical methods underlying our LBM model. Section 2.2 details geometries and conditions used. The obtained results are presented in Section 3 and discussed in Section 4. Section 5 brings concluding lights on this study.

## 2. Methods

### 2.1. Numerical model

The flow of interest here, involving a fluid with properties similar to water at room temperature and therefore a sound speed of $c_s \approx 1450$ m/s, and a maximum characteristic speed of the order of 10 m/s, falls into the low-Mach regime, with Ma below $7 \times 10^{-3}$. As such it can readily be described with a low-Mach approximation to the Navier-Stokes (NS) equations:

$$ \partial_t(\rho u_\alpha) + \partial_\beta(\rho u_\alpha u_\beta) + \partial_\beta T_{\alpha\beta} = 0, \quad (1) $$
where the stress tensor $\mathcal{T}$ is:

$$ \mathcal{T}_{\alpha\beta} = p\delta_{\alpha\beta} - \mu \left( \partial_{\beta}u_{\alpha} + \partial_{\alpha}u_{\beta} - \frac{2}{D}\partial_{\gamma}u_{\gamma}\delta_{\alpha\beta} \right) - \eta\partial_{\gamma}u_{\gamma}\delta_{\alpha\beta}, \quad (2) $$

where $p, \mu, \eta$ and $D$ are respectively the pressure, dynamic viscosity, bulk viscosity and number of dimensions in space; $\delta_{\alpha\beta}$ denotes the Kronecker delta function.

Low-Mach flows can be dealt with in a variety of ways: using compressible NS solvers would inevitably result in rather small time-steps and excessive computational costs due to the large discrepancy between the acoustic and convective modes in terms of propagation speeds; or incompressible NS solvers with an infinite sound speed (i.e. Ma = 0) which, while overcoming the stiffness stemming from acoustics, introduce a non-local elliptic equation for the pressure (Poisson equation). A third approach, inbetween these two, is a compressible solver with rescaled acoustics (still guaranteeing separation of scales between acoustics and hydrodynamics) with larger effective Mach number allowing for a fully hyperbolic/parabolic system of equations and larger time-steps compared to the fully compressible formulation. The now well-known and widely used LBM targets the latter approximation with a discrete set of evolution equations for discrete probability distribution functions, $f_i$ [41]:

$$ f_i(\mathbf{x} + \mathbf{c}_i \delta t, t + \delta t) - f_i(\mathbf{x}, t) = \frac{\delta t}{\bar{\tau}} (f_i^{\text{eq}}(\rho, \mathbf{u}) - f_i(\mathbf{x}, t)), \quad (3) $$

where $c_i$ are discrete particle velocities, $\bar{\tau}$ the relaxation time, $\delta t$ the time-step size and $f_i^{\text{eq}}$ the discrete equilibrium distribution function. This system of equations is readily obtained by integrating the discrete-in-phase-space Boltzmann equation with the Bhatnagar, Gross and Krook (BGK) collision operator along its characteristics [42, 43].

The discrete-in-phase-space Boltzmann equation is obtained by projecting the continuous space of particle velocities of dimension $D$ onto a set of ortho-normal base functions, e.g. Hermite polynomials, and keeping lower-order contributions needed to recover the correct hydrodynamic limit (i.e. NS level dynamics) [44]. As such the discrete equilibrium function is a truncated Hermite expansion of the Maxwell-Boltzmann equilibrium function [44]:

$$ f_i^{\text{eq}}(\rho, \mathbf{u}) = \sum_{n=0}^{N} \frac{w_i}{i! \theta_0^i} \mathbf{a}_n^{\text{eq}}(\rho, \mathbf{u}) : \mathcal{H}_n(\mathbf{c}_i), \quad (4) $$

where for a third-order quadrature the stencil reference temperature is $\theta_0 = \delta r^2/3\delta t^2$, $\mathcal{H}_n(\mathbf{c}_i)$ and $\mathbf{a}_n^{\text{eq}}(\rho, \mathbf{u})$ are the Hermite polynomials and corresponding equilibrium coefficients, while the weights $w_i$ and the set of discrete particle velocities $c_i$ are obtained from a Gauss-Hermite quadrature.

While most of the initial models and simulations relied on second order expansion of the equilibrium distribution, it has since been observed that this leads to Galilean invariance issues in the shear components of the stress ten-

sor. This problem can readily be removed by adding third-order components to the discrete equilibrium [45]. The same issue exists for normal components of the viscous tensor. However, this problem can only be removed via correction terms [13, 45, 46]. The present study targets the incompressible regime where bulk viscosity effects are negligible. Therefore, this correction is not necessary and will not be considered.

The relaxation time appearing in the evolution equation is tied to the local dynamic viscosity as:

$$ \bar{\tau} = \frac{\mu}{\rho \theta_0} + \frac{\delta t}{2}. \quad (5) $$

The original model, as described in eq. (3) is referred to as the single-relaxation time realization of the BGK collision operator, which has a limited range of stability. To enhance the domain of stability and have a more robust scheme allowing for under-resolved simulations we use a multiple relaxation time realization based on central Hermite polynomials. Then, eq. (3) changes into:

$$ f_i(\mathbf{x} + \mathbf{c}_i \delta t, t + \delta t) - f_i(\mathbf{x}, t) = \mathbf{T}^{-1} \mathbf{S} \mathbf{T} (f_i^{\text{eq}}(\rho, \mathbf{u}) - f_i(\mathbf{x}, t)) \quad (6) $$

where $\mathbf{T}$ and $\mathbf{T}^{-1}$ are the moments transform tensor and its inverse defined as:

$$ \mathbf{T} \cdot \mathbf{f} = \tilde{\mathbf{M}}, \quad (7) $$

and $\mathbf{S}$ is the diagonal tensor of relaxation coefficients. The central Hermite moments $\tilde{\mathbf{M}}$ are defined as [4, 39]:

$$ \tilde{\mathbf{M}}_n = \tilde{\mathbf{a}}_n = \sum_i \mathcal{H}_n(\mathbf{c}_i - \mathbf{u}) f_i. \quad (8) $$

While relaxation times appearing in $\mathbf{S}$ can be tuned individually for each moment, here all ghost relaxations are set to one, as this both minimizes computational costs and is very close to the optimal linear stability manifold in the space of ghost relaxation times [39]. For more details on the transforms and moments readers are referred to [4, 46].

## 2.2. Numerical simulations

**Cases description.** A cross-sectional view of the idealized medical device, subject of the FDA nozzle benchmark is shown in Figure 1. The geometry is composed of an axisymmetric nozzle with a convergent section on the flow inlet extremity, a constant throat section in the middle and a sudden expansion section toward the end. To eliminate a possible numerical effect of inlet and outlet boundary conditions, the nozzle has been extended by 0.052 m and 0.02 m at both ends, respectively. The final length of the nozzle is 0.24 m. To validate the numerical results, data is collected from both the horizontal center-line and eight distinct planes located at different coordinates along the horizontal x-axis. Figure 1 shows the locations of these planes. The simulations are run over a total period of
where the stress tensor $\mathcal{T}$ is:

$$ \mathcal{T}_{\alpha\beta} = p\delta_{\alpha\beta} - \mu \left( \partial_{\beta}u_{\alpha} + \partial_{\alpha}u_{\beta} - \frac{2}{D}\partial_{\gamma}u_{\gamma}\delta_{\alpha\beta} \right) - \eta\partial_{\gamma}u_{\gamma}\delta_{\alpha\beta}, \quad (2) $$

where $p, \mu, \eta$ and $D$ are respectively the pressure, dynamic viscosity, bulk viscosity and number of dimensions in space; $\delta_{\alpha\beta}$ denotes the Kronecker delta function.

Low-Mach flows can be dealt with in a variety of ways: using compressible NS solvers would inevitably result in rather small time-steps and excessive computational costs due to the large discrepancy between the acoustic and convective modes in terms of propagation speeds; or incompressible NS solvers with an infinite sound speed (i.e. Ma = 0) which, while overcoming the stiffness stemming from acoustics, introduce a non-local elliptic equation for the pressure (Poisson equation). A third approach, inbetween these two, is a compressible solver with rescaled acoustics (still guaranteeing separation of scales between acoustics and hydrodynamics) with larger effective Mach number allowing for a fully hyperbolic/parabolic system of equations and larger time-steps compared to the fully compressible formulation. The now well-known and widely used LBM targets the latter approximation with a discrete set of evolution equations for discrete probability distribution functions, $f_i$ [41]:

$$ f_i(\mathbf{x} + \mathbf{c}_i \delta t, t + \delta t) - f_i(\mathbf{x}, t) = \frac{\delta t}{\bar{\tau}} (f_i^{\text{eq}}(\rho, \mathbf{u}) - f_i(\mathbf{x}, t)), \quad (3) $$

where $c_i$ are discrete particle velocities, $\bar{\tau}$ the relaxation time, $\delta t$ the time-step size and $f_i^{\text{eq}}$ the discrete equilibrium distribution function. This system of equations is readily obtained by integrating the discrete-in-phase-space Boltzmann equation with the Bhatnagar, Gross and Krook (BGK) collision operator along its characteristics [42, 43].

The discrete-in-phase-space Boltzmann equation is obtained by projecting the continuous space of particle velocities of dimension $D$ onto a set of ortho-normal base functions, e.g. Hermite polynomials, and keeping lower-order contributions needed to recover the correct hydrodynamic limit (i.e. NS level dynamics) [44]. As such the discrete equilibrium function is a truncated Hermite expansion of the Maxwell-Boltzmann equilibrium function [44]:

$$ f_i^{\text{eq}}(\rho, \mathbf{u}) = \sum_{n=0}^{N} \frac{w_i}{i! \theta_0^i} \mathbf{a}_n^{\text{eq}}(\rho, \mathbf{u}) : \mathcal{H}_n(\mathbf{c}_i), \quad (4) $$

where for a third-order quadrature the stencil reference temperature is $\theta_0 = \delta r^2/3\delta t^2$, $\mathcal{H}_n(\mathbf{c}_i)$ and $\mathbf{a}_n^{\text{eq}}(\rho, \mathbf{u})$ are the Hermite polynomials and corresponding equilibrium coefficients, while the weights $w_i$ and the set of discrete particle velocities $c_i$ are obtained from a Gauss-Hermite quadrature.

While most of the initial models and simulations relied on second order expansion of the equilibrium distribution, it has since been observed that this leads to Galilean invariance issues in the shear components of the stress tensor.

This problem can readily be removed by adding third-order components to the discrete equilibrium [45]. The same issue exists for normal components of the viscous tensor. However, this problem can only be removed via correction terms [13, 45, 46]. The present study targets the incompressible regime where bulk viscosity effects are negligible. Therefore, this correction is not necessary and will not be considered.

The relaxation time appearing in the evolution equation is tied to the local dynamic viscosity as:

$$ \bar{\tau} = \frac{\mu}{\rho \theta_0} + \frac{\delta t}{2}. \quad (5) $$

The original model, as described in eq. (3) is referred to as the single-relaxation time realization of the BGK collision operator, which has a limited range of stability. To enhance the domain of stability and have a more robust scheme allowing for under-resolved simulations we use a multiple relaxation time realization based on central Hermite polynomials. Then, eq. (3) changes into:

$$ f_i(\mathbf{x} + \mathbf{c}_i \delta t, t + \delta t) - f_i(\mathbf{x}, t) = \mathbf{T}^{-1} \mathbf{S} \mathbf{T} (f_i^{\text{eq}}(\rho, \mathbf{u}) - f_i(\mathbf{x}, t)) \quad (6) $$

where $\mathbf{T}$ and $\mathbf{T}^{-1}$ are the moments transform tensor and its inverse defined as:

$$ \mathbf{T} \cdot \mathbf{f} = \tilde{\mathbf{M}}, \quad (7) $$

and $\mathbf{S}$ is the diagonal tensor of relaxation coefficients. The central Hermite moments $\tilde{\mathbf{M}}$ are defined as [4, 39]:

$$ \tilde{\mathbf{M}}_n = \tilde{\mathbf{a}}_n = \sum_i \mathcal{H}_n(\mathbf{c}_i - \mathbf{u}) f_i. \quad (8) $$

While relaxation times appearing in $\mathbf{S}$ can be tuned individually for each moment, here all ghost relaxations are set to one, as this both minimizes computational costs and is very close to the optimal linear stability manifold in the space of ghost relaxation times [39]. For more details on the transforms and moments readers are referred to [4, 46].

## 2.2. Numerical simulations

**Cases description.** A cross-sectional view of the idealized medical device, subject of the FDA nozzle benchmark is shown in Figure 1. The geometry is composed of an axisymmetric nozzle with a convergent section on the flow inlet extremity, a constant throat section in the middle and a sudden expansion section toward the end. To eliminate a possible numerical effect of inlet and outlet boundary conditions, the nozzle has been extended by 0.052 m and 0.02 m at both ends, respectively. The final length of the nozzle is 0.24 m. To validate the numerical results, data is collected from both the horizontal center-line and eight distinct planes located at different coordinates along the horizontal x-axis. Figure 1 shows the locations of these planes. The simulations are run over a total period of
<figure><img src="image_1.png" /><figcaption>Figure 1: Schematic representation of the FDA nozzle with all dimensions.</figcaption></figure>

time covering 300 $t_c$, where $t_c$ is the flow-through period defined as,

$$t_c = \frac{L_c}{u_c}, \qquad (9)$$

where the characteristic speed $u_c$ and length $L_c$ are set to the velocity in the throat $u_{thr}$ and to the length of the throat section, respectively. Only data from the second half, i.e. 150–300 $t_c$ is used for the time-averaging process to get average flow fields. The PIV data used for validation and the 3-D model geometry are available from the FDA official website<fnref n="1" />.

*Flow conditions.* Following experimental settings and benchmark requirements, the fluid was considered to be Newtonian, with a density of 1056 kg/m³ and a dynamic viscosity of 3.5 mPa s. The Reynolds number is computed as:

$$Re = Re_{thr} = \frac{\rho \bar{u}_{thr} d_{thr}}{\mu}, \qquad (10)$$

where $d_{thr}$ is the throat diameter and $\bar{u}_{thr}$ is defined as:

$$\bar{u}_{thr} = \left( \frac{d_0}{d_{thr}} \right)^2 \bar{u}_0. \qquad (11)$$

The average inlet velocity $\bar{u}_0$ is obtained dividing the volume flow rate by the inlet surface area $S_0 = \pi d_0^2/4$. To cover all flow regimes, the present study considers Re=500, 3500, and 6500. Based on critical pipe Reynolds numbers, the first would fall in the laminar regime, the second in the transitional, and the third in the low-level turbulent regime, respectively.

*Simulation details.* The inlet boundary condition for the velocity is prescribed in the form of a laminar parabolic profile:

$$u(r) = 2\bar{u}_0 \left( 1 - 4 \frac{r^2}{d_0^2} \right). \qquad (12)$$

Constant pressure boundary conditions are implemented at the outlet while no-slip boundary conditions are employed along the walls. The no-slip conditions are enforced with the curved treatment of the bounce-back method [47]

while the outflow constant pressure is applied using the non-equilibrium extrapolation approach [48].

The computational domain is tridimensional ($D = 3$) and has a volume of 0.012 × 0.012 × 0.24 m³. Simulations were conducted with different grid resolutions, $\delta r$, listed in table 2. The respective time-step sizes, $\delta t$, are defined by fixing the maximum predicted convective Courant–Friedrichs–Lewy (CFL) number in the domain:

$$CFL = \frac{2\bar{u}_{thr}}{\delta r / \delta t}, \qquad (13)$$

in turn leading to:

$$\delta t = \frac{\rho d_{thr} \delta r \text{CFL}}{2 \mu \text{Re}_{\text{thr}}}. \qquad (14)$$

For the analysis of the turbulent statistics, the velocity is decomposed into a mean and a fluctuating part as:

$$\mathbf{u}(\mathbf{x}, t) = \bar{\mathbf{u}}(\mathbf{x}) + \mathbf{u}'(\mathbf{x}, t), \qquad (15)$$

where the mean part $\bar{\mathbf{u}}$ is defined by:

$$\bar{\mathbf{u}}(\mathbf{x}) = \frac{1}{\Delta t} \int_{t_0}^{t_0+\Delta t} \mathbf{u}(\mathbf{x}, t) dt, \qquad (16)$$

$t_0$ indicates the start of averaging and $\Delta t$ the period over which it is done. Averaging the convective momentum flux term in the NS equations using this decomposition, one gets:

$$\overline{(u_{\alpha} + u'_{\alpha})(u_{\beta} + u'_{\beta})} = \overline{u_{\alpha}u_{\beta}} + \overline{u'_{\alpha}u'_{\beta}}, \qquad (17)$$

The last term of this equation appears in the Reynolds stress tensor, defined as:

$$\Pi'_{\alpha\beta} = \overline{\rho u'_{\alpha} u'_{\beta}}, \qquad (18)$$

where $u'(\mathbf{x}, t)$ and $\rho$ indicate velocity fluctuations and fluid density, respectively. The Reynolds stress tensor is symmetrical and consists of three normal stress and six shear stress components. The normal stress can be calculated as $\rho u'_{\alpha} u'_{\alpha}$ where the term $u'_{\alpha} u'_{\alpha}$ indicates the time-averaged product of normal velocity fluctuations. Finally, the viscous shear stress $\sigma$ is calculated, in this study, as below:

$$\sigma_{\alpha\beta} = \rho\nu(\partial_{\beta}u_{\alpha} + \partial_{\alpha}u_{\beta}). \qquad (19)$$

Four different resolutions will be considered in the study, from R1 (coarsest) to R4 (finest), listed in table 2.

<footnote n="1">¹https://ncihub.org/wiki/FDA_CFD</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Schematic representation of the FDA nozzle with all dimensions.</figcaption></figure>

time covering 300 $t_c$, where $t_c$ is the flow-through period defined as,

$$t_c = \frac{L_c}{u_c}, \qquad (9)$$

where the characteristic speed $u_c$ and length $L_c$ are set to the velocity in the throat $u_{thr}$ and to the length of the throat section, respectively. Only data from the second half, i.e. 150–300 $t_c$ is used for the time-averaging process to get average flow fields. The PIV data used for validation and the 3-D model geometry are available from the FDA official website<fnref n="1" />.

*Flow conditions.* Following experimental settings and benchmark requirements, the fluid was considered to be Newtonian, with a density of 1056 kg/m³ and a dynamic viscosity of 3.5 mPa s. The Reynolds number is computed as:

$$Re = Re_{thr} = \frac{\rho \bar{u}_{thr} d_{thr}}{\mu}, \qquad (10)$$

where $d_{thr}$ is the throat diameter and $\bar{u}_{thr}$ is defined as:

$$\bar{u}_{thr} = \left( \frac{d_0}{d_{thr}} \right)^2 \bar{u}_0. \qquad (11)$$

The average inlet velocity $\bar{u}_0$ is obtained dividing the volume flow rate by the inlet surface area $S_0 = \pi d_0^2/4$. To cover all flow regimes, the present study considers Re=500, 3500, and 6500. Based on critical pipe Reynolds numbers, the first would fall in the laminar regime, the second in the transitional, and the third in the low-level turbulent regime, respectively.

*Simulation details.* The inlet boundary condition for the velocity is prescribed in the form of a laminar parabolic profile:

$$u(r) = 2\bar{u}_0 \left(1 - 4\frac{r^2}{d_0^2}\right). \qquad (12)$$

Constant pressure boundary conditions are implemented at the outlet while no-slip boundary conditions are employed along the walls. The no-slip conditions are enforced with the curved treatment of the bounce-back method [47]

while the outflow constant pressure is applied using the non-equilibrium extrapolation approach [48].

The computational domain is tridimensional ($D = 3$) and has a volume of 0.012 × 0.012 × 0.24 m³. Simulations were conducted with different grid resolutions, $\delta r$, listed in table 2. The respective time-step sizes, $\delta t$, are defined by fixing the maximum predicted convective Courant–Friedrichs–Lewy (CFL) number in the domain:

$$CFL = \frac{2\bar{u}_{thr}}{\delta r / \delta t}, \qquad (13)$$

in turn leading to:

$$\delta t = \frac{\rho d_{thr} \delta r \text{CFL}}{2 \mu \text{Re}_{\text{thr}}}. \qquad (14)$$

For the analysis of the turbulent statistics, the velocity is decomposed into a mean and a fluctuating part as:

$$\mathbf{u}(\mathbf{x}, t) = \bar{\mathbf{u}}(\mathbf{x}) + \mathbf{u}'(\mathbf{x}, t), \qquad (15)$$

where the mean part $\bar{\mathbf{u}}$ is defined by:

$$\bar{\mathbf{u}}(\mathbf{x}) = \frac{1}{\Delta t} \int_{t_0}^{t_0+\Delta t} \mathbf{u}(\mathbf{x}, t) dt, \qquad (16)$$

$t_0$ indicates the start of averaging and $\Delta t$ the period over which it is done. Averaging the convective momentum flux term in the NS equations using this decomposition, one gets:

$$\overline{(u_{\alpha} + u'_{\alpha})(u_{\beta} + u'_{\beta})} = \overline{u_{\alpha}u_{\beta}} + \overline{u'_{\alpha}u'_{\beta}}, \qquad (17)$$

The last term of this equation appears in the Reynolds stress tensor, defined as:

$$\Pi'_{\alpha\beta} = \overline{\rho u'_{\alpha} u'_{\beta}}, \qquad (18)$$

where $\mathbf{u}'(\mathbf{x}, t)$ and $\rho$ indicate velocity fluctuations and fluid density, respectively. The Reynolds stress tensor is symmetrical and consists of three normal stress and six shear stress components. The normal stress can be calculated as $\rho u'_{\alpha} u'_{\alpha}$ where the term $\overline{u'_{\alpha}u'_{\alpha}}$ indicates the time-averaged product of normal velocity fluctuations. Finally, the viscous shear stress $\sigma$ is calculated, in this study, as below:

$$\sigma_{\alpha\beta} = \rho\nu(\partial_{\beta}u_{\alpha} + \partial_{\alpha}u_{\beta}). \qquad (19)$$

Four different resolutions will be considered in the study, from R1 (coarsest) to R4 (finest), listed in table 2.

<footnote n="1">¹https://ncihub.org/wiki/FDA_CFD</footnote>
<table><caption>Table 2: Spatial and temporal resolutions of simulations.</caption><thead><tr><th>Re</th><th>δr</th><th>×10<sup>-4</sup>m</th><th>δt</th><th>×10<sup>-6</sup>s</th></tr></thead><tbody><tr><td>R1</td><td>500</td><td>2.3</td><td></td><td>49.03</td></tr><tr><td>R2</td><td>500</td><td>1.725</td><td></td><td>36.75</td></tr><tr><td>R3</td><td>500</td><td>1.15</td><td></td><td>24.5</td></tr><tr><td>R1</td><td>3500</td><td>2.3</td><td></td><td>9.86</td></tr><tr><td>R2</td><td>3500</td><td>1.725</td><td></td><td>7.38</td></tr><tr><td>R3</td><td>3500</td><td>1.15</td><td></td><td>4.93</td></tr><tr><td>R1</td><td>6500</td><td>2.3</td><td></td><td>5.31</td></tr><tr><td>R2</td><td>6500</td><td>1.725</td><td></td><td>3.98</td></tr><tr><td>R3</td><td>6500</td><td>1.15</td><td></td><td>2.65</td></tr><tr><td>R4</td><td>6500</td><td>0.776</td><td></td><td>1.7</td></tr></tbody></table>

## 3. Results

Numerical results are compared to experimental data from three independent laboratories, presented in the form of five different data-sets [24], i.e. three different measurements from the first laboratory, and one each for the second and third laboratories. Throughout this section the three independent sets of measurements from the first laboratory are shown with diamonds, circular, and square markers. Data from the second and third laboratories are represented by asterisks and upward-pointing triangular markers, respectively.

### 3.1. Laminar flow, Re = 500

At Re = 500 (inlet Reynolds number of 167), the flow regime is expected and has been observed to be fully laminar from both PIV and numerical results, meaning that no jet break-down is observed after the sudden expansion (see Figure 2).

<figure><img src="image_2.png" /><figcaption>Figure 2: Results for Re=500. Snapshot of the steady velocity field as obtained from R2 simulations at time t = 2.4 s.</figcaption></figure>

**Axial velocity.** The results as obtained from R1, R2 and R3 resolutions are compared to experimental data in Figures 3 and 4. First, we look into the time-averaged axial velocity (note that, in this laminar configuration, time-averaged fields are equivalent to the steady fields) along the center-line. All three numerical profiles closely match

each other – lying often on top of each other –, indicating that the chosen grid-sizes allow to correctly resolve the flow features.

<figure><img src="image_3.png" /><figcaption>Figure 3: Results for Re = 500. Comparison of mean axial velocity between LBM results and PIV experiments along the center-line of the nozzle. The results from R1, R2, and R3 are indicated by solid black, dashed black and solid red lines respectively.</figcaption></figure>

It can be noted from the axial velocity plot that the flow is constant upstream of the contraction section ($x < -0.063$ m) and accelerates when flowing into the contraction region, reaching a peak value of 0.76 m/s at the entrance of the expansion section, plane $X_2$ (corresponding to $x = 0$). The flow starts decelerating with an approximately linear rate in the sudden expansion section ($x > 0$) due to the increase in cross-section area. In both Figures 3 and 4, large discrepancies are observed between the different measurements. Stewart et al. [29] reported that the experimental data-sets are strongly affected by normalization errors, which could partly explain the observed variations.

**Radial profiles of axial velocity on sampling planes.** The profiles of the predicted axial velocity are illustrated at eight cross-sections in Figure 4. A good quantitative agreement is observed between PIV data and numerical results throughout the nozzle.

However, the peak values of the axial velocity appear to be slightly overestimated after the sudden expansion section. Hariharan et al. [24] reported that there is a less than 1% fluctuation in flow rate during the experiments, which can lead to this difference. Given that discrepancies between experimental and numerical results are at the same level as measurement uncertainties, it can still be concluded that numerical simulations agree well with the experiments for this case.

**Viscous shear stress.** The magnitude of the mean viscous shear stress component $\sigma_{xy}$ along the same cross-sections is illustrated in Figure 5. Overall, numerical viscous shear stresses are in good agreement with the values reported from the PIV measurements. Peak stresses between the central jet and the recirculation zones show a decreasing trend in downstream direction after entering the sudden
<table><caption>Table 2: Spatial and temporal resolutions of simulations.</caption><thead><tr><th>Re</th><th>δr</th><th>×10<sup>-4</sup>m</th><th>δt</th><th>×10<sup>-6</sup>s</th></tr></thead><tbody><tr><td>R1</td><td>500</td><td>2.3</td><td></td><td>49.03</td></tr><tr><td>R2</td><td>500</td><td>1.725</td><td></td><td>36.75</td></tr><tr><td>R3</td><td>500</td><td>1.15</td><td></td><td>24.5</td></tr><tr><td>R1</td><td>3500</td><td>2.3</td><td></td><td>9.86</td></tr><tr><td>R2</td><td>3500</td><td>1.725</td><td></td><td>7.38</td></tr><tr><td>R3</td><td>3500</td><td>1.15</td><td></td><td>4.93</td></tr><tr><td>R1</td><td>6500</td><td>2.3</td><td></td><td>5.31</td></tr><tr><td>R2</td><td>6500</td><td>1.725</td><td></td><td>3.98</td></tr><tr><td>R3</td><td>6500</td><td>1.15</td><td></td><td>2.65</td></tr><tr><td>R4</td><td>6500</td><td>0.776</td><td></td><td>1.7</td></tr></tbody></table>

## 3. Results

Numerical results are compared to experimental data from three independent laboratories, presented in the form of five different data-sets [24], i.e. three different measurements from the first laboratory, and one each for the second and third laboratories. Throughout this section the three independent sets of measurements from the first laboratory are shown with diamonds, circular, and square markers. Data from the second and third laboratories are represented by asterisks and upward-pointing triangular markers, respectively.

### 3.1. Laminar flow, Re = 500

At Re = 500 (inlet Reynolds number of 167), the flow regime is expected and has been observed to be fully laminar from both PIV and numerical results, meaning that no jet break-down is observed after the sudden expansion (see Figure 2).

<figure><img src="image_3.png" /><figcaption>Figure 2: Results for Re=500. Snapshot of the steady velocity field as obtained from R2 simulations at time t = 2.4 s.</figcaption></figure>

**Axial velocity.** The results as obtained from R1, R2 and R3 resolutions are compared to experimental data in Figures 3 and 4. First, we look into the time-averaged axial velocity (note that, in this laminar configuration, time-averaged fields are equivalent to the steady fields) along the center-line. All three numerical profiles closely match

each other – lying often on top of each other –, indicating that the chosen grid-sizes allow to correctly resolve the flow features.

<figure><img src="image_4.png" /><figcaption>Figure 3: Results for Re = 500. Comparison of mean axial velocity between LBM results and PIV experiments along the center-line of the nozzle. The results from R1, R2, and R3 are indicated by solid black, dashed black and solid red lines respectively.</figcaption></figure>

It can be noted from the axial velocity plot that the flow is constant upstream of the contraction section ($x < -0.063$ m) and accelerates when flowing into the contraction region, reaching a peak value of 0.76 m/s at the entrance of the expansion section, plane $X_2$ (corresponding to $x = 0$). The flow starts decelerating with an approximately linear rate in the sudden expansion section ($x > 0$) due to the increase in cross-section area. In both Figures 3 and 4, large discrepancies are observed between the different measurements. Stewart et al. [29] reported that the experimental data-sets are strongly affected by normalization errors, which could partly explain the observed variations.

**Radial profiles of axial velocity on sampling planes.** The profiles of the predicted axial velocity are illustrated at eight cross-sections in Figure 4. A good quantitative agreement is observed between PIV data and numerical results throughout the nozzle.

However, the peak values of the axial velocity appear to be slightly overestimated after the sudden expansion section. Hariharan et al. [24] reported that there is a less than 1% fluctuation in flow rate during the experiments, which can lead to this difference. Given that discrepancies between experimental and numerical results are at the same level as measurement uncertainties, it can still be concluded that numerical simulations agree well with the experiments for this case.

**Viscous shear stress.** The magnitude of the mean viscous shear stress component $\sigma_{xy}$ along the same cross-sections is illustrated in Figure 5. Overall, numerical viscous shear stresses are in good agreement with the values reported from the PIV measurements. Peak stresses between the central jet and the recirculation zones show a decreasing trend in downstream direction after entering the sudden
expansion. The viscous shear stresses are consistently zero along the center-line.

### 3.2. Transitional flow, Re = 3500

At $Re = 3500$ (inlet Reynolds number of 1169), the flow regime has been observed to be transitional with a jet breakdown downstream of the sudden expansion region observed both in the PIV experiments and numerical simulations (see Figure 6).

**Axial velocity.** Looking at Figure 7, the axial velocity reaches a peak of 4 m/s on the finest grid (R3) at $x = 0.02$ m, but then decreases afterwards rapidly to 0.5 m/s at $x = 0.06$ m, which indicates jet break-down. The velocity distributions obtained from all three simulations match perfectly up to the smooth contraction zone. In the throat, small discrepancies (of the order of 3% relative differences) are observed in the predicted maximum velocities. Overall, the time-averaged axial velocities along the center-line predicted by LBM simulations match the PIV data very well, particularly so for resolution R3. The differences observed between R1, R2 and R3 are in the same range as the variations found in the measurements.

**Radial profiles of axial velocity on sampling planes.** Figure 8 shows axial velocity profiles along the chosen cross-section planes. In contrast to the parabolic profiles found at $Re = 500$, the axial velocity profiles are now plug-like in the throat section, as expected. Starting at $X_4$, negative axial velocities appear close to the walls indicating a large, recirculating toroidal zone. Further downstream, the profiles slowly grow into a parabolic distribution with a decreasing peak value, down to 0.5 m/s at the late stages, i.e. the last plane $X_8$. Visible differences between different

grid resolutions are observed most notably at $X_6$, where the highest levels of turbulence (and hence smallest structures) are present, as will be discussed later. The results with the finest mesh R3 (red lines in Figure 8) lead everywhere to an excellent agreement with measurement data.

Even strongly under-resolved simulations, like R1, lead to fully correct and stable qualitative predictions, and stay even in fair quantitative agreement with measurement data outside of the transitional region, demonstrating the robustness of the employed numerical solver.

**Viscous shear stress.** The magnitude of the mean viscous shear stress component $\sigma_{xy}$ is shown in Figure 9. The peak of the viscous shear stress is found at around 25 Pa in the throat section. The stress decreases gradually from plane $X_2$ to the end of the nozzle. The peak is shifted toward the shear layer between the jet and the recirculation zone during jet breakdown. The shear stress remains zero in the boundary layers of the recirculation region at planes $X_3$ and $X_4$, while it starts increasing there at plane $X_5$, before progressively going back to zero further downstream (relaminarization).

**Reynolds stress.** Figures 10 and 11 show the normal Reynolds normal in axial direction at selected cross-sections and along the center-line, respectively. Ahead of the sudden expansion section, the Reynolds stress is nearly zero, indicating laminar conditions.

It experiences a sharp increase around plane $X_4$, reaching maximum values in planes $X_5$ and $X_6$. The peak of the normal Reynolds stress is found to be around 800 Pa for the finest resolution at $X_6$, where strong flow instabilities can be observed from Figure 6. Further downstream

<figure><img src="image_4.png" /><figcaption>Figure 4: Results for $Re = 500$. Axial velocity profiles compared with PIV experiments at selected cross-sections. The radial distance is normalized by the radius of the nozzle, $R_0$. The results from R1, R2, and R3 are indicated by solid black, dashed black and solid red lines respectively. Cross-section planes $X_1$ to $X_8$ are represented from top left to bottom right.</figcaption></figure>
expansion. The viscous shear stresses are consistently zero along the center-line.

### 3.2. Transitional flow, $Re = 3500$

At $Re = 3500$ (inlet Reynolds number of 1169), the flow regime has been observed to be transitional with a jet breakdown downstream of the sudden expansion region observed both in the PIV experiments and numerical simulations (see Figure 6).

**Axial velocity.** Looking at Figure 7, the axial velocity reaches a peak of 4 m/s on the finest grid (R3) at $x = 0.02$ m, but then decreases afterwards rapidly to 0.5 m/s at $x = 0.06$ m, which indicates jet break-down. The velocity distributions obtained from all three simulations match perfectly up to the smooth contraction zone. In the throat, small discrepancies (of the order of 3% relative differences) are observed in the predicted maximum velocities. Overall, the time-averaged axial velocities along the center-line predicted by LBM simulations match the PIV data very well, particularly so for resolution R3. The differences observed between R1, R2 and R3 are in the same range as the variations found in the measurements.

**Radial profiles of axial velocity on sampling planes.** Figure 8 shows axial velocity profiles along the chosen cross-section planes. In contrast to the parabolic profiles found at $Re = 500$, the axial velocity profiles are now plug-like in the throat section, as expected. Starting at $X_4$, negative axial velocities appear close to the walls indicating a large, recirculating toroidal zone. Further downstream, the profiles slowly grow into a parabolic distribution with a decreasing peak value, down to 0.5 m/s at the late stages, i.e. the last plane $X_8$. Visible differences between different

grid resolutions are observed most notably at $X_6$, where the highest levels of turbulence (and hence smallest structures) are present, as will be discussed later. The results with the finest mesh R3 (red lines in Figure 8) lead everywhere to an excellent agreement with measurement data.

Even strongly under-resolved simulations, like R1, lead to fully correct and stable qualitative predictions, and stay even in fair quantitative agreement with measurement data outside of the transitional region, demonstrating the robustness of the employed numerical solver.

**Viscous shear stress.** The magnitude of the mean viscous shear stress component $\sigma_{xy}$ is shown in Figure 9. The peak of the viscous shear stress is found at around 25 Pa in the throat section. The stress decreases gradually from plane $X_2$ to the end of the nozzle. The peak is shifted toward the shear layer between the jet and the recirculation zone during jet breakdown. The shear stress remains zero in the boundary layers of the recirculation region at planes $X_3$ and $X_4$, while it starts increasing there at plane $X_5$, before progressively going back to zero further downstream (relaminarization).

**Reynolds stress.** Figures 10 and 11 show the normal Reynolds normal in axial direction at selected cross-sections and along the center-line, respectively. Ahead of the sudden expansion section, the Reynolds stress is nearly zero, indicating laminar conditions.

It experiences a sharp increase around plane $X_4$, reaching maximum values in planes $X_5$ and $X_6$. The peak of the normal Reynolds stress is found to be around 800 Pa for the finest resolution at $X_6$, where strong flow instabilities can be observed from Figure 6. Further downstream

<figure><img src="image_5.png" /><figcaption>Figure 4: Results for $Re = 500$. Axial velocity profiles compared with PIV experiments at selected cross-sections. The radial distance is normalized by the radius of the nozzle, $R_0$. The results from R1, R2, and R3 are indicated by solid black, dashed black and solid red lines respectively. Cross-section planes $X_1$ to $X_8$ are represented from top left to bottom right.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 5: Results for Re = 500. Viscous shear stress compared with PIV experiments at selected cross-sections along the center-line. The results from R1, R2, and R3 are indicated by solid black, dashed black and solid red lines respectively. Cross-section planes X₁ to X₈ are represented from top left to bottom right. Note the different vertical scales between the two rows.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 6: Results for Re=3500. Snapshots of instantaneous velocity field as obtained from R3 simulation at time $t = 4.2$ s. The top-left inlay is a zoom over the sudden expansion.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 7: Results for Re = 3500. Comparison of mean axial velocity between LBM results and PIV experiments along the center-line of the nozzle. The results from R1, R2 and R3 are indicated by solid black, dashed black and solid red lines respectively.</figcaption></figure>

(cross-section planes X₇ & X₈) the turbulent stress goes back to zero, indicating again relaminarization of the flow. While noticeable differences are observed among the three different resolutions (especially around the jet breakdown region), they are observed to converge toward experimental data, as seen in Figure 11. The numerical results from the finest mesh R3 fall within the PIV data and follow closely the corresponding trends, leading to an accurate description of jet breakdown. It must also be kept in mind that an accurate measurement of the Reynolds stress tensor based on PIV remains extremely challenging due to the needed resolution in space and time. As such, the experimental data used for the present comparison is certainly associated with a non-negligible level of experimental uncertainty, already visible in the widely different results obtained.
<figure><img src="image_6.png" /><figcaption>Figure 5: Results for Re = 500. Viscous shear stress compared with PIV experiments at selected cross-sections along the center-line. The results from R1, R2, and R3 are indicated by solid black, dashed black and solid red lines respectively. Cross-section planes X₁ to X₈ are represented from top left to bottom right. Note the different vertical scales between the two rows.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Results for Re=3500. Snapshots of instantaneous velocity field as obtained from R3 simulation at time $t = 4.2$ s. The top-left inlay is a zoom over the sudden expansion.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Results for Re = 3500. Comparison of mean axial velocity between LBM results and PIV experiments along the center-line of the nozzle. The results from R1, R2 and R3 are indicated by solid black, dashed black and solid red lines respectively.</figcaption></figure>

(cross-section planes X₇ & X₈) the turbulent stress goes back to zero, indicating again relaminarization of the flow. While noticeable differences are observed among the three different resolutions (especially around the jet breakdown region), they are observed to converge toward experimental data, as seen in Figure 11. The numerical results from the finest mesh R3 fall within the PIV data and follow closely the corresponding trends, leading to an accurate description of jet breakdown. It must also be kept in mind that an accurate measurement of the Reynolds stress tensor based on PIV remains extremely challenging due to the needed resolution in space and time. As such, the experimental data used for the present comparison is certainly associated with a non-negligible level of experimental uncertainty, already visible in the widely different results obtained.
<figure><img src="image_8.png" /><figcaption>Figure 8: Results for Re = 3500. Time-averaged axial velocity profiles compared with PIV experiments at selected cross-sections. The radial distance is normalized by the radius of the nozzle, $R_0$. The results from R1, R2 and R3 are indicated by solid black, dashed black and solid red lines respectively. Cross-section planes X₁ to X₈ are represented from top left to bottom right.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 9: Results for Re = 3500. Time-averaged viscous shear stress compared with PIV experiments at selected cross-sections. The results from R1, R2 and R3 are indicated by solid black, dashed black and solid red lines respectively. Cross-section planes X₁ to X₈ are represented from top left to bottom right. Note the different vertical scales between the two rows.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Results for Re = 3500. Time-averaged axial velocity profiles compared with PIV experiments at selected cross-sections. The radial distance is normalized by the radius of the nozzle, $R_0$. The results from R1, R2 and R3 are indicated by solid black, dashed black and solid red lines respectively. Cross-section planes X₁ to X₈ are represented from top left to bottom right.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Results for Re = 3500. Time-averaged viscous shear stress compared with PIV experiments at selected cross-sections. The results from R1, R2 and R3 are indicated by solid black, dashed black and solid red lines respectively. Cross-section planes X₁ to X₈ are represented from top left to bottom right. Note the different vertical scales between the two rows.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 10: Results for Re = 3500. Time-averaged normal Reynolds stress in axial direction compared with PIV experiments at selected cross-sections. The results from R1, R2 and R3 are indicated by solid black, dashed black and solid red lines respectively. Cross-section planes X₁ to X₈ are represented from top left to bottom right. Note the different vertical scales between the two rows.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 11: Results for Re= 3500. Comparison between LBM results and experimental data for normal Reynolds stress in axial direction along the center-line. The results from R1, R2 and R3 are indicated by solid black, dashed black and solid red lines respectively.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 12: Results for Re = 6500. Instantaneous axial velocity field for mesh R1 at time t = 11.925 s.</figcaption></figure>

tained by the three experimental groups involved in the
study.

3.3. Mildly turbulent flow, Re = 6500

At Re = 6500 (inlet Reynolds number of 2171), the flow regime is observed to be locally turbulent downstream of the sudden expansion region in the PIV and in the numerical simulations (see Figure 12).

**Axial velocity.** The profiles of time-averaged axial velocity for the PIV data-sets and for the four LBM meshes (R1, R2, R3 and R4) are shown along the center-line in Figure 13. As in Figure 7 the LBM axial velocity profiles show very good agreement with PIV experimental data-sets. The results with mesh resolution R3 and R4 fall within the experimental symbols. It is worth noting

that the jet from the coarsest mesh, R1, starts breakdown
slightly later than at higher resolution. The axial velocity
increases to a peak value around 7.5 m/s at x = 0.01 m,
before decreasing very sharply to less than 0.5 m/s. This
abrupt decrease marks the onset of jet breakdown. All
four resolutions are able to accurately predict the axial
velocity distribution along the center-line throughout the
FDA nozzle.

Radial profiles of axial velocity along sampling planes. The time-averaged axial velocity profiles are illustrated at the selected cross-sections along the flow direction in Figure 14. It can be seen from the first two plots that the flow remains at the same peak velocity in the throat section; the numerical results from the four resolutions match well with the PIV experiments. Although the flow enters the
<figure><img src="image_11.png" /><figcaption>Figure 10: Results for Re = 3500. Time-averaged normal Reynolds stress in axial direction compared with PIV experiments at selected cross-sections. The results from R1, R2 and R3 are indicated by solid black, dashed black and solid red lines respectively. Cross-section planes X₁ to X₈ are represented from top left to bottom right. Note the different vertical scales between the two rows.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11: Results for Re = 3500. Comparison between LBM results and experimental data for normal Reynolds stress in axial direction along the center-line. The results from R1, R2 and R3 are indicated by solid black, dashed black and solid red lines respectively.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: Results for Re = 6500. Instantaneous axial velocity field for mesh R1 at time $t = 11.925$ s.</figcaption></figure>

tained by the three experimental groups involved in the study.

### 3.3. Mildly turbulent flow, Re = 6500

At Re = 6500 (inlet Reynolds number of 2171), the flow regime is observed to be locally turbulent downstream of the sudden expansion region in the PIV and in the numerical simulations (see Figure 12).

**Axial velocity.** The profiles of time-averaged axial velocity for the PIV data-sets and for the four LBM meshes (R1, R2, R3 and R4) are shown along the center-line in Figure 13. As in Figure 7 the LBM axial velocity profiles show very good agreement with PIV experimental data-sets. The results with mesh resolution R3 and R4 fall within the experimental symbols. It is worth noting

that the jet from the coarsest mesh, R1, starts breakdown slightly later than at higher resolution. The axial velocity increases to a peak value around 7.5 m/s at $x = 0.01$ m, before decreasing very sharply to less than 0.5 m/s. This abrupt decrease marks the onset of jet breakdown. All four resolutions are able to accurately predict the axial velocity distribution along the center-line throughout the FDA nozzle.

**Radial profiles of axial velocity along sampling planes.** The time-averaged axial velocity profiles are illustrated at the selected cross-sections along the flow direction in Figure 14. It can be seen from the first two plots that the flow remains at the same peak velocity in the throat section; the numerical results from the four resolutions match well with the PIV experiments. Although the flow enters the
<figure><img src="image_13.png" /><figcaption>Figure 13: Results for Re = 6500. Comparison between LBM results and PIV experiments for mean axial velocity along the center-line of the nozzle. The results from solid black, dashed black, solid red and dashed red lines respectively.</figcaption></figure>

expansion section at location X₃, the peak axial velocity is still the same as in the previous locations, in good agreement with PIV data. It is worth noting that a small negative velocity appears at cross-section X₄ near the boundary layers, while the peak axial velocity decreases only slightly. On the other hand, the peak axial velocity decreases drastically from section X₅ to X₆ in association with jet breakdown. At later stages, in the cross-sections X₇ and X₈, it can be observed that the peak values are relatively small (around 0.75 m/s). The numerical predictions match very well with the PIV data in this region where the jet collapses almost completely and the flow reattaches. Overall, the LBM results are in very good agreement with the experiments concerning jet breakdown. The mesh R3 captures correctly the axial velocity profile and the jet breakdown process throughout the nozzle.

**Viscous shear stress.** The magnitude of the time-averaged viscous shear stress component $\sigma_{xy}$ is shown at selected cross-sections in Figure 15. It can be noted that the viscous shear stress reaches a peak around 50 Pa in the throat section, the largest velocity gradients being found as expected in the boundary layer. In the sudden expansion section, the peak of the viscous shear stress is predicted to be less than 8 Pa and is found in the shear layers between the central, high-velocity jet and the recirculating fluid. The peak of the viscous shear stress decreases further downstream. Later on, the effect of viscous shear stress can be fully neglected due to the uniform and low velocity found in the nozzle. The largest (but low) shear stress is found close to the boundary layer in cross-sections X₇ and X₈, which can be attributed to the jet reattachment to the wall. As expected, the shear stress is zero along the center-line. The height of the boundary layer is estimated to be around half of the throat diameter in both experimental and numerical results.

Overall, viscous shear can be neglected as the peak viscous stress is very small throughout the nozzle. Such a value should not be able to damage red blood cells. De-

pending on the studies, a shear stress exceeding about 100 Pa is necessary to observe any hemolysis [49].

**Reynolds stress.** Figures 16 and 17 illustrate the normal Reynolds stress for selected cross-planes and along the center-line, respectively. It can be observed that the results are in good agreement with the experimental data throughout the nozzle, keeping in mind the large variations also found in the experiments. In the throat section (planes X₁& X₂), the numerical results are close to zero, as in most measurements (while few experimental sets show local peaks up to 1500 Pa; the experimental results vary very widely there). After flowing through the sudden expansion region, the Reynolds stress increases largely until the late stage of jet breakdown. The peak of Reynolds stress appears in this region in the shear layers between the central high-speed jet and the surrounding recirculating fluid, reaching around 3000 Pa between cross-sections X₅ and X₆. Further downstream, the Reynolds stress starts decreasing sharply towards zero. The Reynolds stress is negligible in the throat as well as at the end of the nozzle.

## 4. Discussion

In this section, further information is provided concerning agreement and/or discrepancies between numerical results and PIV measurements in this study.

**Laminar flow.** The in-house numerical tool ALBORZ successfully predicts the nature of the laminar flow. The jet induced after leaving the throat does not break down in the sudden expansion section. The numerical simulations accurately capture the fields of velocity and viscous shear stress.

**Transitional flow.** Both numerical simulations and experiments show a transitional regime with a jet breakdown in the sudden expansion section at this Reynolds number of 3500. The location of the jet breakdown is accurately predicted by the simulations, around $x = 0.02$ m, which is later than for Re = 6500. The effect of the viscous shear stress would be negligible regarding hemolysis, due to its low peak value. Concerning the Reynolds stress, used to quantify the strength of turbulent velocity fluctuations and to assess indirectly the potential damage to blood cells. The obtained values are much higher and would be sufficient to trigger hemolysis. However, as explained in Yu et al. [49], there is still no clear consensus regarding critical values and unsteady effects must be taken into account to clearly assess the impact of stress on hemolysis.

**Mildly turbulent flow.** The turbulent flow exhibits an earlier onset of jet break-down, already around $x = 0.01$ m. Ahead of the sudden expansion section, the value of the Reynolds stresses are small, showing that there is no notable flow disturbance produced before and within the throat, as expected. Some discrepancies are found between numerical results and PIV measurements regarding
<figure><img src="image_14.png" /><figcaption>Figure 13: Results for Re = 6500. Comparison between LBM results and PIV experiments for mean axial velocity along the center-line of the nozzle. The results from solid black, dashed black, solid red and dashed red lines respectively.</figcaption></figure>

expansion section at location X₃, the peak axial velocity is still the same as in the previous locations, in good agreement with PIV data. It is worth noting that a small negative velocity appears at cross-section X₄ near the boundary layers, while the peak axial velocity decreases only slightly. On the other hand, the peak axial velocity decreases drastically from section X₅ to X₆ in association with jet breakdown. At later stages, in the cross-sections X₇ and X₈, it can be observed that the peak values are relatively small (around 0.75 m/s). The numerical predictions match very well with the PIV data in this region where the jet collapses almost completely and the flow reattaches. Overall, the LBM results are in very good agreement with the experiments concerning jet breakdown. The mesh R3 captures correctly the axial velocity profile and the jet breakdown process throughout the nozzle.

**Viscous shear stress.** The magnitude of the time-averaged viscous shear stress component $\sigma_{xy}$ is shown at selected cross-sections in Figure 15. It can be noted that the viscous shear stress reaches a peak around 50 Pa in the throat section, the largest velocity gradients being found as expected in the boundary layer. In the sudden expansion section, the peak of the viscous shear stress is predicted to be less than 8 Pa and is found in the shear layers between the central, high-velocity jet and the recirculating fluid. The peak of the viscous shear stress decreases further downstream. Later on, the effect of viscous shear stress can be fully neglected due to the uniform and low velocity found in the nozzle. The largest (but low) shear stress is found close to the boundary layer in cross-sections X₇ and X₈, which can be attributed to the jet reattachment to the wall. As expected, the shear stress is zero along the center-line. The height of the boundary layer is estimated to be around half of the throat diameter in both experimental and numerical results.

Overall, viscous shear can be neglected as the peak viscous stress is very small throughout the nozzle. Such a value should not be able to damage red blood cells. De-

pending on the studies, a shear stress exceeding about 100 Pa is necessary to observe any hemolysis [49].

**Reynolds stress.** Figures 16 and 17 illustrate the normal Reynolds stress for selected cross-planes and along the center-line, respectively. It can be observed that the results are in good agreement with the experimental data throughout the nozzle, keeping in mind the large variations also found in the experiments. In the throat section (planes X₁& X₂), the numerical results are close to zero, as in most measurements (while few experimental sets show local peaks up to 1500 Pa; the experimental results vary very widely there). After flowing through the sudden expansion region, the Reynolds stress increases largely until the late stage of jet breakdown. The peak of Reynolds stress appears in this region in the shear layers between the central high-speed jet and the surrounding recirculating fluid, reaching around 3000 Pa between cross-sections X₅ and X₆. Further downstream, the Reynolds stress starts decreasing sharply towards zero. The Reynolds stress is negligible in the throat as well as at the end of the nozzle.

## 4. Discussion

In this section, further information is provided concerning agreement and/or discrepancies between numerical results and PIV measurements in this study.

**Laminar flow.** The in-house numerical tool ALBORZ successfully predicts the nature of the laminar flow. The jet induced after leaving the throat does not break down in the sudden expansion section. The numerical simulations accurately capture the fields of velocity and viscous shear stress.

**Transitional flow.** Both numerical simulations and experiments show a transitional regime with a jet breakdown in the sudden expansion section at this Reynolds number of 3500. The location of the jet breakdown is accurately predicted by the simulations, around $x = 0.02$ m, which is later than for Re = 6500. The effect of the viscous shear stress would be negligible regarding hemolysis, due to its low peak value. Concerning the Reynolds stress, used to quantify the strength of turbulent velocity fluctuations and to assess indirectly the potential damage to blood cells. The obtained values are much higher and would be sufficient to trigger hemolysis. However, as explained in Yu et al. [49], there is still no clear consensus regarding critical values and unsteady effects must be taken into account to clearly assess the impact of stress on hemolysis.

**Mildly turbulent flow.** The turbulent flow exhibits an earlier onset of jet break-down, already around $x = 0.01$ m. Ahead of the sudden expansion section, the value of the Reynolds stresses are small, showing that there is no notable flow disturbance produced before and within the throat, as expected. Some discrepancies are found between numerical results and PIV measurements regarding
<figure><img src="image_14.png" /><figcaption>Figure 14: Results for Re = 6500. Time-averaged axial velocity compared with PIV experiments at the chosen cross-sections along the center-line. The radial distance is normalized by the radius of the nozzle, R₀. The results from R1, R2, R3 and R4 are indicated by solid black, dashed black, solid red and dashed red lines respectively. Cross-section planes X₁ to X₈ are represented from top left to bottom right.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 15: Results for Re = 6500. Time-averaged viscous shear stress compared with PIV experiments at selected cross-sections along the center-line. The results from R1, R2, R3 and R4 are indicated by solid black, dashed black, solid red and dashed red lines respectively. Cross-section planes X₁ to X₈ are represented from top left to bottom right. Note the different vertical scales between the two rows.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 14: Results for Re = 6500. Time-averaged axial velocity compared with PIV experiments at the chosen cross-sections along the center-line. The radial distance is normalized by the radius of the nozzle, R₀. The results from R1, R2, R3 and R4 are indicated by solid black, dashed black, solid red and dashed red lines respectively. Cross-section planes X₁ to X₈ are represented from top left to bottom right.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 15: Results for Re = 6500. Time-averaged viscous shear stress compared with PIV experiments at selected cross-sections along the center-line. The results from R1, R2, R3 and R4 are indicated by solid black, dashed black, solid red and dashed red lines respectively. Cross-section planes X₁ to X₈ are represented from top left to bottom right. Note the different vertical scales between the two rows.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 16: Results for Re=6500. Comparison between LBM results and experimental data for axial Reynolds normal stress at different cross-sections. The results from R1, R2, R3 and R4 are indicated by solid black, dashed black, solid red and dashed red lines respectively. Cross-section planes X₁ to X₈ are represented from top left to bottom right.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 17: Results for Re=6500. Comparison between LBM results and experimental data for axial Reynolds normal stress along the center-line. The results from R1, R2, R3 and R4 are indicated by solid black, dashed black, solid red and dashed red lines respectively.</figcaption></figure>

Reynolds normal stress for cross-section planes X₅ and X₆. However, the differences are below the variations observed among the different experimental data-sets, the orders of magnitude are well predicted, and the shapes of the profiles coincide. It is therefore impossible to decide if the disagreements come from the numerical model or from experimental limitations.

## 5. Conclusions

Lattice Boltzmann simulations were performed on the FDA's benchmark nozzle, mimicking a simplified blood-contacting medical device, to assess the capability of the in-house solver ALBORZ in modelling such flows under laminar, transitional, and mildly turbulent conditions. For this purpose, Reynolds numbers Re = 500, 3500, and 6500

have been considered, for which experimental data are available. Velocity, viscous shear stress, and turbulent Reynolds stress were compared to the inter-laboratory PIV experiments documented in the literature [24]. Overall, numerical results are in very good agreement with the experimental data, demonstrating that the solver is capable of predicting accurately the underlying flow physics. This applies even for under-resolved simulations (grid R1), allowing therefore very fast estimates, since results on grid R1 are typically obtained almost 20 times faster than on grid R3. When a very high accuracy is required, by increasing the resolution, results converge toward experimental profiles within the corresponding uncertainty. This shows that the solver can be used in a stable manner either 1) for a fast estimate of the flow conditions with a fair accuracy, or 2) to get very accurate flow predictions at the price of a much longer computational time. To our knowledge, this work is the first LBM study of the FDA's benchmark nozzle considering simultaneously laminar, transitional, and mildly turbulent flow conditions, providing comprehensive quantitative comparisons for the flow fields.

## Acknowledgements

F.H. would like to acknowledge the financial support of China Scholarship Council (grant number 201908080236). S.A.H. acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) in TRR 287 (Project-ID 422037413). P.B. acknowledges the funding by the Federal Ministry of Education and Research in Germany within the Forschungscampus STIMULATE (grant number 13GW0473A) and the German Research Foundation (grant number BE 6230/6 -1). The authors
<figure><img src="image_17.png" /><figcaption>Figure 16: Results for Re=6500. Comparison between LBM results and experimental data for axial Reynolds normal stress at different cross-sections. The results from R1, R2, R3 and R4 are indicated by solid black, dashed black, solid red and dashed red lines respectively. Cross-section planes X₁ to X₈ are represented from top left to bottom right.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 17: Results for Re=6500. Comparison between LBM results and experimental data for axial Reynolds normal stress along the center-line. The results from R1, R2, R3 and R4 are indicated by solid black, dashed black, solid red and dashed red lines respectively.</figcaption></figure>

Reynolds normal stress for cross-section planes X₅ and X₆. However, the differences are below the variations observed among the different experimental data-sets, the orders of magnitude are well predicted, and the shapes of the profiles coincide. It is therefore impossible to decide if the disagreements come from the numerical model or from experimental limitations.

## 5. Conclusions

Lattice Boltzmann simulations were performed on the FDA's benchmark nozzle, mimicking a simplified blood-contacting medical device, to assess the capability of the in-house solver ALBORZ in modelling such flows under laminar, transitional, and mildly turbulent conditions. For this purpose, Reynolds numbers Re = 500, 3500, and 6500

have been considered, for which experimental data are available. Velocity, viscous shear stress, and turbulent Reynolds stress were compared to the inter-laboratory PIV experiments documented in the literature [24]. Overall, numerical results are in very good agreement with the experimental data, demonstrating that the solver is capable of predicting accurately the underlying flow physics. This applies even for under-resolved simulations (grid R1), allowing therefore very fast estimates, since results on grid R1 are typically obtained almost 20 times faster than on grid R3. When a very high accuracy is required, by increasing the resolution, results converge toward experimental profiles within the corresponding uncertainty. This shows that the solver can be used in a stable manner either 1) for a fast estimate of the flow conditions with a fair accuracy, or 2) to get very accurate flow predictions at the price of a much longer computational time. To our knowledge, this work is the first LBM study of the FDA's benchmark nozzle considering simultaneously laminar, transitional, and mildly turbulent flow conditions, providing comprehensive quantitative comparisons for the flow fields.

## Acknowledgements

F.H. would like to acknowledge the financial support of China Scholarship Council (grant number 201908080236). S.A.H. acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) in TRR 287 (Project-ID 422037413). P.B. acknowledges the funding by the Federal Ministry of Education and Research in Germany within the Forschungscampus STIMULATE (grant number 13GW0473A) and the German Research Foundation (grant number BE 6230/6 -1). The authors
also gratefully acknowledge the Gauss centre for providing computation time under grants "pn29du" and "pn73ta" on the GCS supercomputer superMUC-NG at Leibniz supercomputing centre, Munich, Germany.

## Declaration of competing interest

None.

## Appendix A. Central Hermite moments collision operator

In the present work the continuous particle-speed space is discretized via a third-order Gauss-Hermite quadrature. The third-order quadrature leads to stencils of the form $DdQ^3d$, where $d$ is the number of physical dimensions, i.e. $d \in \{1, 2, 3\}$. Furthermore the third-order quadrature yields the following weights in 1-D:

$$w_{i,\alpha} = \begin{cases} 2/3 & \text{for } c_{i,\alpha} = 0, \\ 1/6 & \text{for } c_{i,\alpha} = \delta r/\delta t, \end{cases} \quad (A.1)$$

which in turn for the 3-D stencil is merely the tensorial product of the 1-D weights:

$$w_i = \prod_{\alpha=x,y,z} w_{i,\alpha}, \quad (A.2)$$

and $\theta_0 = \delta r^2/3\delta t^2$. The Hermite polynomials appearing in Eq. 4 are defined as:

$$\mathcal{H}_0 = 1, \quad (A.3a)$$

$$\mathcal{H}_{\alpha} = c_{i,\alpha}, \quad (A.3b)$$

$$\mathcal{H}_{\alpha\beta} = c_{i,\alpha}c_{i,\beta} - \delta_{\alpha\beta}, \quad (A.3c)$$

$$\mathcal{H}_{\alpha\beta\gamma} = c_{i,\alpha}c_{i,\beta}c_{i,\gamma} - \text{cyc}(\delta_{\alpha\beta}c_{i,\gamma}), \quad (A.3d)$$

$$\mathcal{H}_{\alpha\beta\gamma\zeta} = c_{i,\alpha}c_{i,\beta}c_{i,\gamma}c_{i,\zeta} - \text{cyc}(\delta_{\alpha\beta}c_{i,\gamma}c_{i,\zeta}) \\ + \text{cyc}(\delta_{\alpha\beta}\delta_{\gamma\zeta}), \quad (A.3e)$$

$$\mathcal{H}_{\alpha\beta\gamma\zeta\iota} = c_{i,\alpha}c_{i,\beta}c_{i,\gamma}c_{i,\zeta}c_{i,\iota} - \text{cyc}(\delta_{\zeta\iota}c_{i,\alpha}c_{i,\beta}c_{i,\gamma}) \\ + \text{cyc}(c_{i,\alpha}\delta_{\beta\gamma}\delta_{\zeta\iota}), \quad (A.3f)$$

$$\mathcal{H}_{\alpha\beta\gamma\zeta\iota\kappa} = c_{i,\alpha}c_{i,\beta}c_{i,\gamma}c_{i,\zeta}c_{i,\iota}c_{i,\kappa} - \text{cyc}(c_{i,\alpha}c_{i,\beta}c_{i,\gamma}c_{i,\zeta}\delta_{\iota\kappa}) \\ + \text{cyc}(c_{i,\alpha}c_{i,\beta}\delta_{\gamma\zeta}\delta_{\iota\kappa}) - \text{cyc}(\delta_{\alpha\beta}\delta_{\gamma\zeta}\delta_{\iota\kappa}), \quad (A.3g)$$

with cyc a cyclic permutation and resulting equilibrium coefficients:

$$\mathbf{a}_0^{\text{eq}} = \rho, \quad (A.4a)$$

$$\mathbf{a}_{\alpha}^{\text{eq}} = \rho u_{\alpha}, \quad (A.4b)$$

$$\mathbf{a}_{\alpha\beta}^{\text{eq}} = \rho u_{\alpha} u_{\beta}, \quad (A.4c)$$

$$\mathbf{a}_{\alpha\beta\gamma}^{\text{eq}} = \rho u_{\alpha} u_{\beta} u_{\gamma}, \quad (A.4d)$$

$$\mathbf{a}_{\alpha\beta\gamma\zeta}^{\text{eq}} = \rho u_{\alpha} u_{\beta} u_{\gamma} u_{\zeta}, \quad (A.4e)$$

$$\mathbf{a}_{\alpha\beta\gamma\zeta\iota}^{\text{eq}} = \rho u_{\alpha} u_{\beta} u_{\gamma} u_{\zeta} u_{\iota}, \quad (A.4f)$$

$$\mathbf{a}_{\alpha\beta\gamma\zeta\iota\kappa}^{\text{eq}} = \rho u_{\alpha} u_{\beta} u_{\gamma} u_{\zeta} u_{\iota} u_{\kappa}. \quad (A.4g)$$

In the central moments space, assuming a D3Q27 stencil, the following equilibrium moments are recovered:

$$T f^{\text{eq}} = \tilde{a}^{\text{eq}}, \quad (A.5)$$

where $\tilde{a}^{\text{eq}}$ are the equilibrium Hermite coefficients in central moments space:

$$\tilde{a}_{0}^{\text{eq}} = \rho, \quad (A.6a)$$

$$\tilde{a}_{n}^{\text{eq}} = 0, \forall n \neq 0. \quad (A.6b)$$

It must be noted that all equilibrium moments of order $n > 0$ become zero in central space only using the full Hermite expansion supported by the stencil. In the present work we use a sixth-order Hermite-expansion of the discrete equilibrium for all simulations.

Apart from the relaxation rates of second-order moments tied to the fluid viscosity, i.e. $a_{xx}$, $a_{yy}$, $a_{zz}$, $a_{xy}$, $a_{xz}$ and $a_{yz}$, all other relaxation rates are set to 1. A set of Matlab code are included as supplementary material to compute the transform and inverse tensors along with corresponding moments for the collision model used in the present work.

It is also interesting to note that the viscous stress tensor obtained as the second-order moment of the non-equilibrium part of the distribution function [50] can be readily obtained from Hermite coefficients:

$$\frac{1}{2} (\partial_{\alpha} u_{\beta} + \partial_{\beta} u_{\alpha}) = \frac{1}{2\rho\theta_{0}\bar{\tau}} (a_{\alpha\beta} - a_{\alpha\beta}^{\text{eq}}), \quad (A.7)$$

or central Hermite coefficients:

$$\frac{1}{2} (\partial_{\alpha} u_{\beta} + \partial_{\beta} u_{\alpha}) = \frac{1}{2\rho\theta_{0}\bar{\tau}} (\tilde{a}_{\alpha\beta}). \quad (A.8)$$

## References

[1] S. Succi, The lattice Boltzmann equation: for fluid dynamics and beyond, Oxford university press, 2001.

[2] R. Noël, F. Ge, Y. Zhang, L. Navarro, G. Courbebaisse, Lattice Boltzmann method for modelling of biological phenomena, in: 25th European Signal Processing Conference (EUSIPCO), IEEE, (2017), pp. 2654–2658.

[3] R. Ouared, B. Chopard, B. Stahl, D. Rüfenacht, H. Yilmaz, G. Courbebaisse, Thrombosis modeling in intracranial aneurysms: a lattice Boltzmann numerical algorithm, Computer Physics Communications 179 (1-3) (2008) 128–131.

[4] S. A. Hosseini, P. Berg, F. Huang, C. Roloff, G. Janiga, D. Thévenin, Central moments multiple relaxation time LBM for hemodynamic simulations in intracranial aneurysms: An in vitro validation study using PIV and PC-MRI, Computers in Biology and Medicine 131 (2021) 104251.

[5] K. Jain, Efficacy of the FDA nozzle benchmark and the lattice Boltzmann method for the analysis of biomedical flows in transitional regime, Medical & Biological Engineering & Computing 58 (2020) 1817–1830.

[6] M. Hasert, K. Masilamani, S. Zimny, H. Klimach, J. Qi, J. Bernsdorf, S. Roller, Complex fluid simulations with the parallel tree-based Lattice Boltzmann solver Musubi, Journal of Computational Science 5 (5) (2014) 784–794.

[7] O. Malaspinas, A. Turjman, D. Ribeiro de Sousa, G. Garcia-Cardena, M. Raes, P.-T. T. Nguyen, Y. Zhang, G. Courbebaisse, C. Lelubre, K. Zouaoui Boudjeltia, B. Chopard, A spatio-temporal model for spontaneous thrombus formation in cerebral aneurysms, Journal of Theoretical Biology 394 (2016) 68–76.
also gratefully acknowledge the Gauss centre for providing computation time under grants "pn29du" and "pn73ta" on the GCS supercomputer superMUC-NG at Leibniz supercomputing centre, Munich, Germany.

## Declaration of competing interest

None.

## Appendix A. Central Hermite moments collision operator

In the present work the continuous particle-speed space is discretized via a third-order Gauss-Hermite quadrature. The third-order quadrature leads to stencils of the form $DdQ^3d$, where $d$ is the number of physical dimensions, i.e. $d \in \{1, 2, 3\}$. Furthermore the third-order quadrature yields the following weights in 1-D:

$$w_{i,\alpha} = \begin{cases} 2/3 & \text{for } c_{i,\alpha} = 0, \\ 1/6 & \text{for } c_{i,\alpha} = \delta r/\delta t, \end{cases} \quad (A.1)$$

which in turn for the 3-D stencil is merely the tensorial product of the 1-D weights:

$$w_i = \prod_{\alpha=x,y,z} w_{i,\alpha}, \quad (A.2)$$

and $\theta_0 = \delta r^2/3\delta t^2$. The Hermite polynomials appearing in Eq. 4 are defined as:

$$\mathcal{H}_0 = 1, \quad (A.3a)$$

$$\mathcal{H}_{\alpha} = c_{i,\alpha}, \quad (A.3b)$$

$$\mathcal{H}_{\alpha\beta} = c_{i,\alpha}c_{i,\beta} - \delta_{\alpha\beta}, \quad (A.3c)$$

$$\mathcal{H}_{\alpha\beta\gamma} = c_{i,\alpha}c_{i,\beta}c_{i,\gamma} - \text{cyc}(\delta_{\alpha\beta}c_{i,\gamma}), \quad (A.3d)$$

$$\mathcal{H}_{\alpha\beta\gamma\zeta} = c_{i,\alpha}c_{i,\beta}c_{i,\gamma}c_{i,\zeta} - \text{cyc}(\delta_{\alpha\beta}c_{i,\gamma}c_{i,\zeta}) \\ + \text{cyc}(\delta_{\alpha\beta}\delta_{\gamma\zeta}), \quad (A.3e)$$

$$\mathcal{H}_{\alpha\beta\gamma\zeta\iota} = c_{i,\alpha}c_{i,\beta}c_{i,\gamma}c_{i,\zeta}c_{i,\iota} - \text{cyc}(\delta_{\zeta\iota}c_{i,\alpha}c_{i,\beta}c_{i,\gamma}) \\ + \text{cyc}(c_{i,\alpha}\delta_{\beta\gamma}\delta_{\zeta\iota}), \quad (A.3f)$$

$$\mathcal{H}_{\alpha\beta\gamma\zeta\iota\kappa} = c_{i,\alpha}c_{i,\beta}c_{i,\gamma}c_{i,\zeta}c_{i,\iota}c_{i,\kappa} - \text{cyc}(c_{i,\alpha}c_{i,\beta}c_{i,\gamma}c_{i,\zeta}\delta_{\iota\kappa}) \\ + \text{cyc}(c_{i,\alpha}c_{i,\beta}\delta_{\gamma\zeta}\delta_{\iota\kappa}) - \text{cyc}(\delta_{\alpha\beta}\delta_{\gamma\zeta}\delta_{\iota\kappa}), \quad (A.3g)$$

with cyc a cyclic permutation and resulting equilibrium coefficients:

$$\mathbf{a}_0^{\text{eq}} = \rho, \quad (A.4a)$$

$$\mathbf{a}_{\alpha}^{\text{eq}} = \rho u_{\alpha}, \quad (A.4b)$$

$$\mathbf{a}_{\alpha\beta}^{\text{eq}} = \rho u_{\alpha} u_{\beta}, \quad (A.4c)$$

$$\mathbf{a}_{\alpha\beta\gamma}^{\text{eq}} = \rho u_{\alpha} u_{\beta} u_{\gamma}, \quad (A.4d)$$

$$\mathbf{a}_{\alpha\beta\gamma\zeta}^{\text{eq}} = \rho u_{\alpha} u_{\beta} u_{\gamma} u_{\zeta}, \quad (A.4e)$$

$$\mathbf{a}_{\alpha\beta\gamma\zeta\iota}^{\text{eq}} = \rho u_{\alpha} u_{\beta} u_{\gamma} u_{\zeta} u_{\iota}, \quad (A.4f)$$

$$\mathbf{a}_{\alpha\beta\gamma\zeta\iota\kappa}^{\text{eq}} = \rho u_{\alpha} u_{\beta} u_{\gamma} u_{\zeta} u_{\iota} u_{\kappa}. \quad (A.4g)$$

In the central moments space, assuming a D3Q27 stencil, the following equilibrium moments are recovered:

$$Tf^{\text{eq}} = \tilde{a}^{\text{eq}}, \quad (A.5)$$

where $\tilde{a}^{\text{eq}}$ are the equilibrium Hermite coefficients in central moments space:

$$\tilde{a}_{0}^{\text{eq}} = \rho, \quad (A.6a)$$

$$\tilde{a}_{n}^{\text{eq}} = 0, \forall n \neq 0. \quad (A.6b)$$

It must be noted that all equilibrium moments of order $n > 0$ become zero in central space only using the full Hermite expansion supported by the stencil. In the present work we use a sixth-order Hermite-expansion of the discrete equilibrium for all simulations.

Apart from the relaxation rates of second-order moments tied to the fluid viscosity, i.e. $a_{xx}$, $a_{yy}$, $a_{zz}$, $a_{xy}$, $a_{xz}$ and $a_{yz}$, all other relaxation rates are set to 1. A set of Matlab code are included as supplementary material to compute the transform and inverse tensors along with corresponding moments for the collision model used in the present work.

It is also interesting to note that the viscous stress tensor obtained as the second-order moment of the non-equilibrium part of the distribution function [50] can be readily obtained from Hermite coefficients:

$$\frac{1}{2} (\partial_{\alpha} u_{\beta} + \partial_{\beta} u_{\alpha}) = \frac{1}{2\rho\theta_{0}\bar{\tau}} (a_{\alpha\beta} - a_{\alpha\beta}^{\text{eq}}), \quad (A.7)$$

or central Hermite coefficients:

$$\frac{1}{2} (\partial_{\alpha} u_{\beta} + \partial_{\beta} u_{\alpha}) = \frac{1}{2\rho\theta_{0}\bar{\tau}} (\tilde{a}_{\alpha\beta}). \quad (A.8)$$

## References

[1] S. Succi, The lattice Boltzmann equation: for fluid dynamics and beyond, Oxford university press, 2001.

[2] R. Noël, F. Ge, Y. Zhang, L. Navarro, G. Courbebaisse, Lattice Boltzmann method for modelling of biological phenomena, in: 25th European Signal Processing Conference (EUSIPCO), IEEE, (2017), pp. 2654–2658.

[3] R. Ouared, B. Chopard, B. Stahl, D. Rüfenacht, H. Yilmaz, G. Courbebaisse, Thrombosis modeling in intracranial aneurysms: a lattice Boltzmann numerical algorithm, Computer Physics Communications 179 (1-3) (2008) 128–131.

[4] S. A. Hosseini, P. Berg, F. Huang, C. Roloff, G. Janiga, D. Thévenin, Central moments multiple relaxation time LBM for hemodynamic simulations in intracranial aneurysms: An in vitro validation study using PIV and PC-MRI, Computers in Biology and Medicine 131 (2021) 104251.

[5] K. Jain, Efficacy of the FDA nozzle benchmark and the lattice Boltzmann method for the analysis of biomedical flows in transitional regime, Medical & Biological Engineering & Computing 58 (2020) 1817–1830.

[6] M. Hasert, K. Masilamani, S. Zimny, H. Klimach, J. Qi, J. Bernsdorf, S. Roller, Complex fluid simulations with the parallel tree-based Lattice Boltzmann solver Musubi, Journal of Computational Science 5 (5) (2014) 784–794.

[7] O. Malaspinas, A. Turjman, D. Ribeiro de Sousa, G. Garcia-Cardena, M. Raes, P.-T. T. Nguyen, Y. Zhang, G. Courbebaisse, C. Lelubre, K. Zouaoui Boudjeltia, B. Chopard, A spatio-temporal model for spontaneous thrombus formation in cerebral aneurysms, Journal of Theoretical Biology 394 (2016) 68–76.
[8] J. McCullough, P. Coveney, An efficient, localised approach for the simulation of elastic blood vessels using the lattice Boltzmann method, arXiv preprint arXiv:2108.08783 (2021).

[9] J. McCullough, P. Coveney, High fidelity physiological blood flow in patient-specific arteriovenous fistula for clinical applications, arXiv preprint arXiv:2012.04639 (2020).

[10] O. E. Kadri, V. D. Chandran, M. Surblyte, R. S. Voronov, In vivo measurement of blood clot mechanics from computational fluid dynamics based on intravital microscopy images, Computers in biology and medicine 106 (2019) 1–11.

[11] A. Khosravanian, M. Rahmanimanesh, P. Keshavarzi, S. Mozafarari, Fast level set method for glioma brain tumor segmentation based on superpixel fuzzy clustering and lattice boltzmann method, Computer Methods and Programs in Biomedicine 198 (2021) 105809.

[12] H. H. Afrouzi, M. Ahmadian, M. Hosseini, H. Arasteh, D. Toghraie, S. Rostami, Simulation of blood flow in arteries with aneurysm: Lattice boltzmann approach (lbm), Computer methods and programs in biomedicine 187 (2020) 105312.

[13] N. I. Prasianakis, I. V. Karlin, Lattice Boltzmann method for simulation of compressible flows on standard lattices, Physical Review E 78 (1) (2008) 016704.

[14] P. J. Dellar, Incompressible limits of lattice Boltzmann equations using multiple relaxation times, Journal of Computational Physics 190 (2) (2003) 351–370.

[15] S. A. Hosseini, C. Coreixas, N. Darabiha, D. Thévenin, Stability of the lattice kinetic scheme and choice of the free relaxation parameter, Physical Review E 99 (6) (2019) 063305.

[16] S. Ansumali, I. V. Karlin, Single relaxation time model for entropic lattice Boltzmann methods, Physical Review E 65 (5) (May 2002).

[17] I. V. Karlin, F. Bösch, S. S. Chikatamarla, Gibbs’ principle for the lattice-kinetic theory of fluid dynamics, Physical Review E 90 (3) (2014) 031302.

[18] D. d’Humieres, I. Ginzburg, M. Krafczyk, P. Lallemand, L.-S. Luo, Multiple-relaxation-time lattice Boltzmann models in three dimensions, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 360 (1792) (2002) 437–451.

[19] M. Geier, A. Greiner, J. G. Korvink, Cascaded digital lattice boltzmann automata for high reynolds number flow, Phys. Rev. E 73 (2006) 066705.

[20] M. Geier, M. Schönherr, A. Pasquali, M. Krafczyk, The cumulant lattice Boltzmann equation in three dimensions: Theory and validation, Computers & Mathematics with Applications 70 (4) (2015) 507–547.

[21] O. Malaspinas, Increasing stability and accuracy of the lattice Boltzmann scheme: Recursivity and regularization, arXiv preprint arXiv:1505.06900 (2015).

[22] M. Geier, A. Pasquali, M. Schönherr, Parametrization of the cumulant lattice Boltzmann method for fourth order accurate diffusion part II: Application to flow around a sphere at drag crisis, Journal of Computational Physics 348 (2017) 889–898.

[23] M. Geier, S. Lenz, M. Schönherr, M. Krafczyk, Under-resolved and large eddy simulations of a decaying Taylor–Green vortex with the cumulant lattice Boltzmann method, Theoretical and Computational Fluid Dynamics 35 (2) (2021) 169–208.

[24] P. Hariharan, M. Giarra, V. Reddy, S. W. Day, K. B. Manning, S. Deutsch, S. F. Stewart, M. R. Myers, M. R. Berman, G. W. Burgreen, et al., 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) (2011).

[25] N. Sánchez Abad, R. Vinuesa, P. Schlatter, M. Andersson, M. Karlsson, Simulation strategies for the Food and Drug Administration nozzle using Nek5000, AIP Advances 10 (2) (2020) 025033.

[26] N. Fehn, W. A. Wall, M. Kronbichler, 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 35 (12) (2019) e3228.

[27] A. W. Bergersen, M. Mortensen, K. Valen-Sendstad, 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) (2019) e3150.

[28] R. Pewowaruk, Y. Li, D. Rowinski, A. Roldán-Alzate, Solution adaptive refinement of cut-cell cartesian meshes can improve FDA nozzle computational fluid dynamics efficiency, International Journal for Numerical Methods in Biomedical Engineering 37 (4) (2021) e3432.

[29] S. F. Stewart, E. G. Paterson, G. W. Burgreen, P. Hariharan, M. Giarra, V. Reddy, S. W. Day, K. B. Manning, S. Deutsch, M. R. Berman, et al., 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) (2012) 139–160.

[30] V. Zmijanovic, S. Mendez, V. Moureau, F. Nicoud, About the numerical robustness of biomedical benchmark cases: Interlaboratory FDA’s idealized medical device, International Journal for numerical methods in biomedical engineering 33 (1) (2017) e02789.

[31] E. L. Manchester, X. Y. Xu, The effect of turbulence on transitional flow in the FDA’s benchmark nozzle model using large-eddy simulation, International Journal for Numerical Methods in Biomedical Engineering 36 (10) (2020) e3389.

[32] S. Bhushan, D. K. Walters, G. W. Burgreen, Laminar, turbulent, and transitional simulations in benchmark cases with cardiovascular device features, Cardiovascular Engineering and Technology 4 (4) (2013) 408–426.

[33] G. Janiga, Large eddy simulation of the FDA benchmark nozzle for a Reynolds number of 6500, Computers in Biology and Medicine 47 (2014) 113–119.

[34] V. Chabannes, C. Prud’Homme, M. Szopos, R. Tarabay, High order finite element simulations for fluid dynamics validated by experimental data from the FDA benchmark nozzle model, arXiv preprint arXiv:1701.02179 (2017).

[35] A. T. White, C. K. Chong, Rotational invariance in the three-dimensional lattice Boltzmann method is dependent on the choice of lattice, Journal of Computational Physics 230 (16) (2011) 6367–6378.

[36] Y. T. Delorme, K. Anupindi, S. H. Frankel, Large eddy simulation of FDA’s idealized medical device, Cardiovascular Engineering and Technology 4 (4) (2013) 392–407.

[37] S. A. Hosseini, N. Darabiha, D. Thévenin, Theoretical and numerical analysis of the lattice kinetic scheme for complex-flow simulations, Physical Review E 99 (2) (2019) 023305.

[38] S. A. Hosseini, H. Safari, N. Darabiha, D. Thévenin, M. Krafczyk, Hybrid Lattice Boltzmann-finite difference model for low mach number combustion simulation, Combustion and Flame 209 (2019) 394–404.

[39] S. A. Hosseini, Development of a lattice Boltzmann-based numerical method for the simulation of reacting flows, Ph.D. thesis, Université Paris-Saclay; Otto-von-Guericke-Universität Magdeburg (2020).

[40] S. Hosseini, A. Abdelsamie, N. Darabiha, D. Thévenin, Low-mach hybrid lattice Boltzmann-finite difference solver for combustion in complex flows, Physics of Fluids 32 (7) (2020) 077105.

[41] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E. M. Viggen, The lattice Boltzmann method, Springer International Publishing 10 (978-3) (2017) 4–15.

[42] X. He, L.-S. Luo, Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation, Physical Review E 56 (6) (1997) 6811.

[43] X. He, L.-S. Luo, A priori derivation of the lattice boltzmann equation, Physical Review E 55 (6) (1997) R6333.

[44] X. Shan, X.-F. Yuan, H. Chen, Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation, Journal of Fluid Mechanics 550 (2006) 413–441.

[45] P. J. Dellar, Lattice Boltzmann algorithms without cubic defects
[8] J. McCullough, P. Coveney, An efficient, localised approach for the simulation of elastic blood vessels using the lattice Boltzmann method, arXiv preprint arXiv:2108.08783 (2021).

[9] J. McCullough, P. Coveney, High fidelity physiological blood flow in patient-specific arteriovenous fistula for clinical applications, arXiv preprint arXiv:2012.04639 (2020).

[10] O. E. Kadri, V. D. Chandran, M. Surblyte, R. S. Voronov, In vivo measurement of blood clot mechanics from computational fluid dynamics based on intravital microscopy images, Computers in biology and medicine 106 (2019) 1–11.

[11] A. Khosravanian, M. Rahmanimanesh, P. Keshavarzi, S. Mozafarari, Fast level set method for glioma brain tumor segmentation based on superpixel fuzzy clustering and lattice boltzmann method, Computer Methods and Programs in Biomedicine 198 (2021) 105809.

[12] H. H. Afrouzi, M. Ahmadian, M. Hosseini, H. Arasteh, D. Toghraie, S. Rostami, Simulation of blood flow in arteries with aneurysm: Lattice boltzmann approach (lbm), Computer methods and programs in biomedicine 187 (2020) 105312.

[13] N. I. Prasianakis, I. V. Karlin, Lattice Boltzmann method for simulation of compressible flows on standard lattices, Physical Review E 78 (1) (2008) 016704.

[14] P. J. Dellar, Incompressible limits of lattice Boltzmann equations using multiple relaxation times, Journal of Computational Physics 190 (2) (2003) 351–370.

[15] S. A. Hosseini, C. Coreixas, N. Darabiha, D. Thévenin, Stability of the lattice kinetic scheme and choice of the free relaxation parameter, Physical Review E 99 (6) (2019) 063305.

[16] S. Ansumali, I. V. Karlin, Single relaxation time model for entropic lattice Boltzmann methods, Physical Review E 65 (5) (May 2002).

[17] I. V. Karlin, F. Bösch, S. S. Chikatamarla, Gibbs’ principle for the lattice-kinetic theory of fluid dynamics, Physical Review E 90 (3) (2014) 031302.

[18] D. d’Humieres, I. Ginzburg, M. Krafczyk, P. Lallemand, L.-S. Luo, Multiple-relaxation-time lattice Boltzmann models in three dimensions, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 360 (1792) (2002) 437–451.

[19] M. Geier, A. Greiner, J. G. Korvink, Cascaded digital lattice boltzmann automata for high reynolds number flow, Phys. Rev. E 73 (2006) 066705.

[20] M. Geier, M. Schönherr, A. Pasquali, M. Krafczyk, The cumulant lattice Boltzmann equation in three dimensions: Theory and validation, Computers & Mathematics with Applications 70 (4) (2015) 507–547.

[21] O. Malaspinas, Increasing stability and accuracy of the lattice Boltzmann scheme: Recursivity and regularization, arXiv preprint arXiv:1505.06900 (2015).

[22] M. Geier, A. Pasquali, M. Schönherr, Parametrization of the cumulant lattice Boltzmann method for fourth order accurate diffusion part II: Application to flow around a sphere at drag crisis, Journal of Computational Physics 348 (2017) 889–898.

[23] M. Geier, S. Lenz, M. Schönherr, M. Krafczyk, Under-resolved and large eddy simulations of a decaying Taylor–Green vortex with the cumulant lattice Boltzmann method, Theoretical and Computational Fluid Dynamics 35 (2) (2021) 169–208.

[24] P. Hariharan, M. Giarra, V. Reddy, S. W. Day, K. B. Manning, S. Deutsch, S. F. Stewart, M. R. Myers, M. R. Berman, G. W. Burgreen, et al., 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) (2011).

[25] N. Sánchez Abad, R. Vinuesa, P. Schlatter, M. Andersson, M. Karlsson, Simulation strategies for the Food and Drug Administration nozzle using Nek5000, AIP Advances 10 (2) (2020) 025033.

[26] N. Fehn, W. A. Wall, M. Kronbichler, 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 35 (12) (2019) e3228.

[27] A. W. Bergersen, M. Mortensen, K. Valen-Sendstad, 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) (2019) e3150.

[28] R. Pewowaruk, Y. Li, D. Rowinski, A. Roldán-Alzate, Solution adaptive refinement of cut-cell cartesian meshes can improve FDA nozzle computational fluid dynamics efficiency, International Journal for Numerical Methods in Biomedical Engineering 37 (4) (2021) e3432.

[29] S. F. Stewart, E. G. Paterson, G. W. Burgreen, P. Hariharan, M. Giarra, V. Reddy, S. W. Day, K. B. Manning, S. Deutsch, M. R. Berman, et al., 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) (2012) 139–160.

[30] V. Zmijanovic, S. Mendez, V. Moureau, F. Nicoud, About the numerical robustness of biomedical benchmark cases: Interlaboratory FDA’s idealized medical device, International Journal for numerical methods in biomedical engineering 33 (1) (2017) e02789.

[31] E. L. Manchester, X. Y. Xu, The effect of turbulence on transitional flow in the FDA’s benchmark nozzle model using large-eddy simulation, International Journal for Numerical Methods in Biomedical Engineering 36 (10) (2020) e3389.

[32] S. Bhushan, D. K. Walters, G. W. Burgreen, Laminar, turbulent, and transitional simulations in benchmark cases with cardiovascular device features, Cardiovascular Engineering and Technology 4 (4) (2013) 408–426.

[33] G. Janiga, Large eddy simulation of the FDA benchmark nozzle for a Reynolds number of 6500, Computers in Biology and Medicine 47 (2014) 113–119.

[34] V. Chabannes, C. Prud’Homme, M. Szopos, R. Tarabay, High order finite element simulations for fluid dynamics validated by experimental data from the FDA benchmark nozzle model, arXiv preprint arXiv:1701.02179 (2017).

[35] A. T. White, C. K. Chong, Rotational invariance in the three-dimensional lattice Boltzmann method is dependent on the choice of lattice, Journal of Computational Physics 230 (16) (2011) 6367–6378.

[36] Y. T. Delorme, K. Anupindi, S. H. Frankel, Large eddy simulation of FDA’s idealized medical device, Cardiovascular Engineering and Technology 4 (4) (2013) 392–407.

[37] S. A. Hosseini, N. Darabiha, D. Thévenin, Theoretical and numerical analysis of the lattice kinetic scheme for complex-flow simulations, Physical Review E 99 (2) (2019) 023305.

[38] S. A. Hosseini, H. Safari, N. Darabiha, D. Thévenin, M. Krafczyk, Hybrid Lattice Boltzmann-finite difference model for low mach number combustion simulation, Combustion and Flame 209 (2019) 394–404.

[39] S. A. Hosseini, Development of a lattice Boltzmann-based numerical method for the simulation of reacting flows, Ph.D. thesis, Université Paris-Saclay; Otto-von-Guericke-Universität Magdeburg (2020).

[40] S. Hosseini, A. Abdelsamie, N. Darabiha, D. Thévenin, Low-mach hybrid lattice Boltzmann-finite difference solver for combustion in complex flows, Physics of Fluids 32 (7) (2020) 077105.

[41] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E. M. Viggen, The lattice Boltzmann method, Springer International Publishing 10 (978-3) (2017) 4–15.

[42] X. He, L.-S. Luo, Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation, Physical Review E 56 (6) (1997) 6811.

[43] X. He, L.-S. Luo, A priori derivation of the lattice boltzmann equation, Physical Review E 55 (6) (1997) R6333.

[44] X. Shan, X.-F. Yuan, H. Chen, Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation, Journal of Fluid Mechanics 550 (2006) 413–441.

[45] P. J. Dellar, Lattice Boltzmann algorithms without cubic defects
in Galilean invariance on standard lattices, Journal of Computational Physics 259 (2014) 270–283.

[46] S. A. Hosseini, N. Darabiha, D. Thévenin, Compressibility in lattice Boltzmann on standard stencils: effects of deviation from reference temperature, Philosophical Transactions of the Royal Society A 378 (2175) (2020) 20190399.

[47] M. Bouzidi, M. Firdaouss, P. Lallemand, Momentum transfer of a Boltzmann-lattice fluid with boundaries, Physics of fluids 13 (11) (2001) 3452–3459.

[48] G. Zhao-Li, Z. Chu-Guang, S. Bao-Chang, Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method, Chinese Physics 11 (4) (2002) 366.

[49] H. Yu, S. Engel, G. Janiga, D. Thévenin, A review of hemolysis prediction models for computational fluid dynamics, Artificial Organs 41 (2017) 603–621.

[50] T. Krüger, F. Varnik, D. Raabe, Shear stress in lattice boltzmann simulations, Physical Review E 79 (4) (2009) 046704.
in Galilean invariance on standard lattices, Journal of Compu-
tational Physics 259 (2014) 270–283.

[46] S. A. Hosseini, N. Darabiha, D. Thévenin, Compressibility in
lattice Boltzmann on standard stencils: effects of deviation from
reference temperature, Philosophical Transactions of the Royal
Society A 378 (2175) (2020) 20190399.

[47] M. Bouzidi, M. Firdaouss, P. Lallemand, Momentum transfer of a Boltzmann-lattice fluid with boundaries, Physics of fluids 13 (11) (2001) 3452–3459.

[48] G. Zhao-Li, Z. Chu-Guang, S. Bao-Chang, Non-equilibrium ex-
trapolation method for velocity and pressure boundary condi-
tions in the lattice Boltzmann method, Chinese Physics 11 (4)
(2002) 366.

[49] H Yu, S. Engel, G. Janiga, D. Thévenin, A review of hemolysis prediction models for computational fluid dynamics, Artificial Organs 41 (2017) 603–621.

[50] T. Krüger, F. Varnik, D. Raabe, Shear stress in lattice boltzmann simulations, Physical Review E 79 (4) (2009) 046704.