Central moments multiple relaxation time LBM for hemodynamic simulations in intracranial aneurysms: An in-vitro validation study using PIV and PC-MRI
S.A. Hosseini, P. Berg, F. Huang, C. Roloff, G. Janiga, D. Thévenin
CCentral moments multiple relaxation time LBM for hemodynamic simulations inintracranial aneurysms: An in-vitro validation study using PIV and PC-MRI
Seyed Ali Hosseini a,b, ∗ , Philipp Berg a,c , Feng Huang a , Christoph Roloff a , G´abor Janiga a , Dominique Th´evenin a a Laboratory of Fluid Dynamics and Technical Flows, University of Magdeburg “Otto von Guericke”, D-39106 Magdeburg, Germany b Department of Mechanical and Process Engineering, ETH Z¨urich, 8092 Z¨urich, Switzerland c Research Campus STIMULATE, University of Magdeburg “Otto von Guericke”, D-39106 Magdeburg, Germany
Abstract
The lattice Boltzmann method (LBM) has recently emerged as an efficient alternative to classical Navier-Stokes solvers.This is particularly true for hemodynamics in complex geometries. However, in its most basic formulation, i.e. withthe so-called single relaxation time (SRT) collision operator, it has been observed to have a limited stability domain inthe Courant/Fourier space, strongly constraining the minimum time-step and grid size. The development of improvedcollision models such as the multiple relaxation time (MRT) operator in central moments space has tremendously widenedthe stability domain, while allowing to overcome a number of other well-documented artifacts, therefore opening the doorfor simulations over a wider range of grid and time-step sizes. The present work focuses on implementing and validatinga specific collision operator, the central Hermite moments multiple relaxation time model with the full expansion of theequilibrium distribution function, to simulate blood flows in intracranial aneurysms. The study further proceeds witha validation of the numerical model through different test-cases and against experimental measurements obtained viastereoscopic particle image velocimetry (PIV) and phase-contrast magnetic resonance imaging (PC-MRI). For a patient-specific aneurysm both PIV and PC-MRI agree fairly well with the simulation. Finally, low-resolution simulations wereshown to be able to capture blood flow information with sufficient accuracy, as demonstrated through both qualitativeand quantitative analysis of the flow field while leading to strongly reduced computation times. For instance in the caseof the patient-specific configuration, increasing the grid-size by a factor of two led to a reduction of computation timeby a factor of 14 with very good similarity indices still ranging from 0.83 to 0.88.
Keywords:
Lattice Boltzmann method; Particle Image Velocimetry; Intracranial aneurysm; Magnetic ResonanceImaging; Validation; Single relaxation time; Central Hermite multiple relaxation time.
1. Introduction
Intracranial aneurysms (IA) are local malformations ofthe cerebral vasculature, which occur with an estimatedprevalence of approximately 3% in the western popula-tion [1]. It is known that aneurysms can grow and de-velop into a stable state or rupture, if the hemodynamicforces exceed the vascular resistance. Since a reliable andclinically applicable rupture risk assessment procedure re-mains challenging until now, increasing research effort isbeing put into developing one. Owing to their non-invasivenature resulting in risk-free (for the patient) assessmentof rupture probability and the possibility to access highlyresolved velocity fields (in both space and time) computa-tional fluid dynamics (CFD) studies are becoming increas-ingly popular [2, 3].While valuable insights have been gained through such nu-merical studies [4, 5], a number of challenges and gaps ∗ Corresponding author
Email address: [email protected] (Seyed Ali Hosseini) remain to be bridged, before predictive and reliable mod-els based on CFD for aneurysm dynamics can be devel-oped [6, 7]. To put these challenges into perspective, Berget al. [8] summarized the individual working steps involvedin image-based blood flow simulations and provided cor-responding recommendations to avoid simulation inaccu-racies. As for any numerical simulation, the choice of thenumerical solver, resolution and grid configuration are ofthe utmost importance to ensure convergence of the ob-tained solutions. The latter two also being consequencesof the choice of the numerical method, the former is adetermining factor concerning both convergence and effi-ciency of the solver.The performances of the employed solvers can only be as-sessed through a systematic benchmarking/validation pro-cedure against reliable experimental data. One of the mostfrequently used sources of experimental reference data isthe particle image velocimetry (PIV) relying on laser-basedhigh-speed camera flow measurements [9, 10, 11, 12, 13].PIV measurements allow for more accurate measurementsand higher resolutions as compared to other data acqui-
Preprint submitted to ... January 28, 2021 a r X i v : . [ phy s i c s . c o m p - ph ] J a n omenclature H α,n Hermite polynomial of order n T Central Hermite moments transform matrix a ( eq ) n Hermite equilibrium coefficient of order n c α Discrete particle velocity in [m/s] u Velocity in [m/s] δ t Time-step size in [s] δ x Grid size in [m]Ω α Collision operator ρ Density in [kg/m ] τ Relaxation time in [s] (cid:101) a ( eq ) n Central Hermite equilibrium coefficient of order nc s lattice sound speed in [m/s] f α Discrete probability distribution function f ( eq ) α Discrete equilibrium distribution function w α Weight associated to discrete velocity c α sition tools such as magnetic resonance imaging (MRI)[14] or cerebral angiography [15]. A number of studieshave presented qualitative comparisons of velocity fieldsfor patient-specific aneurysms as obtained from CFD andPIV measurements, and reported good qualitative agree-ment between the two [16, 17].While most of the early publications relied on discretesolvers for the Navier-Stokes-Poisson (NSP) equations, e.g.using finite volume (FV) or finite element (FE) methods,emergent numerical methods such as smoothed particlehydrodynamics (SPH) [18, 19] and the lattice Boltzmannmethod (LBM) [20, 21, 22, 23, 24, 25, 26, 27] are increas-ingly applied to such flows. Contrary to classical NSPsolvers they are not strictly incompressible, as they ap-proximate the low Mach flow regime through appropri-ate isothermal flow manifolds, and as such rely on a sys-tem of purely hyperbolic equations, making the algorithminherently local. The absence of the elliptic componentof the NSP equations allows for dramatic computationcost reduction while the parabolic nature of the partialdifferential equation (PDE) describing pressure evolutionmakes the formulation suitable for unsteady simulations –contrary to other weakly compressible formulations suchas the artificial compressibility method (ACM). Further-more, the compressible nature of the formulation (involv-ing a thermodynamic pressure) along with kinetic heuris-tic boundary closures (such as the bounce-back rule) pro-vide for an efficient and consistent implementation of wallboundaries [28, 29, 30]. The drastic decrease in compu-tation time makes the method more viable compared toapproaches based on the NSP equations. Based on theseobservations, a number of dedicated LBM-based solvers forblood flow simulations have been developed over the pastdecade, e.g. [21, 31, 32]. For example a research prototypewas developed by the Siemens Healthineers AG (Erlangen,Germany), which is increasingly used to address clinicalquestions [33, 34, 35]. However, it must be noted that theLBM in its simplest form, i.e. with a single relaxationtime (SRT) collision operator and a second-order discreteequilibrium distribution function (EDF) has a rather lim- ited stability domain. The SRT formulation also leads toa number of other well-documented numerical artifacts in-cluding, but not limited to, the (non-dimensional) viscosity-dependence of the solid wall position when used with bounce-back-type boundary treatments [36].The aim of the present work is to assess the performancesof a specific class of multiple relaxation time (MRT) op-erator based on Hermite central moments and a fully ex-panded EDF. This collision model is shown to alleviatesome of the traditional shortcomings of the classical SRT-based solvers, while – through its much wider stabilitydomain – allowing for stable low-resolution simulations.The performances of the central Hermite multiple relax-ation time (CHMRT) model are assessed, via our in-housesolver ALBORZ [37, 38], through a variety of test-casesspanning ideal and patient-specific geometries consider-ing steady and pulsatile flows. The numerical results arecompared with high-resolution in-vitro measurements forvalidation. They are also compared to results from aSRT solver with second-order EDF to further showcasethe added value of the CHMRT collision operator. Theissue of under-resolved simulations is also considered bysystematically conducting simulations at different (lower)resolutions. It is shown that at low resolution (inaccessibleto the SRT collision operator), the proposed scheme is stillable to capture properly the flow dynamics.
2. Materials and methods
The LBM is a solver for the Boltzmann equation in thelimit of the hydrodynamic regime [39]. The time-evolutionof the discrete probability distribution functions is writtenas [30]: f α ( x + c α δ t , t + δ t ) = f α ( x , t ) + δ t Ω α , (1)where f α are the discrete distribution functions, c α thecorresponding particle velocities, δ t the time-step size and2 α the discretized particle collision operator. The colli-sion operator Ω α is usually approximated via a Batnagar-Gross-Krook (BGK) linear relaxation model defined as [40]:Ω α = 1 τ (cid:16) f ( eq ) α − f α (cid:17) , (2)where f ( eq ) α is the EDF. The EDF used in LBM solvers isa truncated approximation (using a Taylor-McLaurin ex-pansion in the limit of vanishing Mach numbers [41] or aHermite expansion [42]) to the Maxwell-Boltzmann distri-bution. It can be written as [43]: f ( eq ) α = w α N (cid:88) n =0 n ! c ns a ( eq ) n : H α,n , (3)where w α are weights associated to each discrete popula-tion, c s is the non-dimensional sound speed at the referencetemperature tied to the time-step and grid sizes, H α,n theHermite polynomial of order n and a ( eq ) n the Hermite co-efficient of the corresponding order. Conserved field vari-ables appearing in the EDF, i.e. density and momentum,are computed from distribution functions as: ρ = (cid:88) α f α , (4) ρ u = (cid:88) α c α f α , (5)where ρ and u are the fluid density and velocity. As briefly mentioned in the introduction, the form ofthe LBM most commonly used in medical applications,i.e. with the SRT collision operator and a second-orderEDF, is subject to a number of well-documented short-comings. It can be readily shown through a multi-scaleperturbation analysis, that at the macroscopic scale, theviscous stress tensor recovered by the second-order EDFadmits deviations from the target Newtonian stress scal-ing as O (Ma ) [30]. In the LBM literature this deviationis usually referred to as the Galilean invariance problem(rightly so as it ties the effective kinematic viscosity to lo-cal velocity) of the stress tensor. It can readily be shownthat the Galilean invariance of the viscous stress tensorcan be restored by extending the Hermite expansion ofthe EDF to higher orders (at least order three). For thebulk viscosity however, given the bias between the firstand third-order moments introduced by the limited num-ber of discrete velocities, one must introduce an additionalcorrection term. The bulk viscosity not being of much in-terest nor relevance in the near-incompressible regime, itwill not be discussed further. Detailed numerical and the-oretical proofs can be found in [44, 45, 46, 47, 48].One of the major advantages of the LBM is the way bound-ary conditions can be implemented. The bounce-back ruleis one of the most popular approaches to enforce wall, ve-locity and pressure boundary conditions. However, it has been shown that when used with the SRT operator it issubject to numerical artifacts such as the (non-dimensional)viscosity-dependence of the solid wall position [36]. Thisshortcoming can be readily proven through asymptoticanalysis of the corresponding system of discrete equationsor simple permeability studies of flow in porous media [38,30]. The seminal work of I. Ginzburg and subsequent de-velopment of a collision operator with a wall position in-dependent from the relaxation coefficient led to the tworelaxation time (TRT) operator [49]. The so-called gen-eralized collision operator first proposed in [50, 51] waslater parametrized to fix its shortcomings [36]. It hasbeen shown that in order to control the wall position in-dependently from the viscosity, at least two different re-laxation coefficients are needed (one for odd, one for evenmoments).Last but not least, it is common knowledge that the classi-cal SRT model is very sensitive to the Fourier number de-fined as Fo = νδ t δ x . It is practically unusable for Fo < . The classical BGK collision operator with one relax-ation coefficient having stability issues for vanishing non-dimensional viscosities, in the context of the present studywe use the CHMRT model defined as:Ω α = T − S T (cid:16) f ( eq ) − f (cid:17) , (6)where T is the moments transform matrix, and S themoments relaxation rate diagonal matrix. In the chosenmoments space, assuming a D3Q27 stencil, the followingequilibrium moments are recovered: T f ( eq ) = (cid:101) a ( eq ) , (7)3here (cid:101) a ( eq ) are the equilibrium Hermite coefficients in cen-tral moments space: (cid:101) a ( eq )0 = ρ, (8a) (cid:101) a ( eq ) n = 0 , ∀ n (cid:54) = 0 . (8b)It must be noted that these equilibrium moments are onlyrecovered using the full Hermite expansion supported bythe stencil, and not the classical second-order EDF usuallyemployed in LBM solvers. Given the added benefits of afully expanded EDF [43, 47], the present work will use asixth-order Hermite-expansion for all simulations.Apart from the relaxation rates of second-order moments,i.e. a xx , a yy , a zz , a xy , a xz and a yz , tied to the fluid viscos-ity, the remainder of the relaxation rates are freely-tunableparameters. Appropriate choices of these free parametershave been shown to widen the linear stability domain ofthe solver and reduce dispersion errors [46, 47]. Given the complex nature of geometries considered inaneuryms, appropriate treatment of wall boundaries andtheir curvature is an important point. In the context ofthe present study the curved-boundary bounce-back for-mulation proposed in [59] is used. At a given boundarynode x f , the missing incoming populations are computedas: f α ( x f , t + δ t ) = 2 qf ¯ α ( x f + c ¯ α , t + δ t )+ (1 − q ) f ¯ α ( x f , t + δ t ) , ∀ q < , (9a) f α ( x f , t + δ t ) = 12 q f ¯ α ( x f + c ¯ α , t + δ t )+ 2 q − q f ¯ α ( x f , t + δ t ) , ∀ q ≥ , (9b)where ¯ α designates the direction opposite α and q : q = || x f − x s |||| c α || , (10)with x s denoting the wall position in direction α . To illustrate the points brought forward in the previ-ous section and showcase the performances of the CHMRTcollision operator for the applications of interest, three dif-ferent test-cases are considered: (a) ideal aneurysm understeady flow, (b) ideal aneurysm subject to pulsatile flow,and finally (c) patient-specific geometry. While the lat-ter two are intended as detailed validations of the solveragainst experimental data, the first test-case is presentedto briefly illustrate the difference between the classicalSRT operator and the proposed model.
The first configuration is a rather simple one consistingof an ideal spherical aneurysm of radius 20 mm and vesselsof diameter 6 mm positioned with an angle of 90 ◦ . Thedistance between the inlet and the opposite point in theaneurysm sac is 60 mm while for the outlet it is 50 mm.The domain is subject to a constant inflow at the inlet andan open boundary at the outflow. The case is modeled us-ing both the CHMRT and second-order SRT collision oper-ators at different resolutions. The obtained results are fur-ther validated against the commercial finite-volume codeSTAR-CCM+ 14.04 (Siemens Product Lifecycle Manage-ment Software Inc., Plato, TX, USA). For the purpose of validation, an idealized sphericalsidewall aneurysm with a diameter of 20 mm was virtu-ally created. The parent vessel featured a diameter of4 mm and was bended under an angle of 120 ◦ (see Fig. 1).The 3-D vessel geometry served as basis for the numer-ical simulation as well as for the phantom manufactur-ing prior to the in-vitro PIV measurements. Basedon the 3-D CAD model, a transparent phantom modelwas manufactured using a lost-core technique that resultedin a silicone block incorporating the hollow vessel struc-ture. The refractive index of the silicone was measuredto be n silicone = 1 . ◦ C (Abbemat 200, AntonPaar, Ostfildern, Germany). As blood analogue liquid(BAL), a mixture of distilled water, glycerin, sodium io-dide, and sodium thiosulphate was used to match the re-fractive index of the silicone block ( n BAL = 1 . ρ BAL = 1221 kg/m . The kinematic viscos-ity of the fluid used for this configuration was ν BAL =3 . × − m / s. This corresponds to a dynamic viscos-ity of µ blood = 4 . × − Pa.s (assuming a density of ρ blood = 1222 kg/m ).During the PIV measurements, the phantom was placedinside a transparent acrylic box with two inclined walls andfilled with index-matching fluid. A laser light sheet wasdirected to illuminate the sagittal plane of the aneurysm.As seeding for the PIV measurements, small resin micro-spheres doped with Rhodamine B (diameter d=10 . ± . µ m, density ρ = 1510 kg/m ) were used. The twostereoscopic PIV high-speed cameras (sCMOS, 2560 × max = 1025 and Womersley number Wo = 2 . y z xyz Figure 1: (left) Idealized aneurysm geometry as used both for the phantom manufacturing and subsequent PIV measurement and in LBsimulation [60], and (right) Segmentation of the patient-specific intracranial aneurysm model used for the LBM simulations and the in-vitrovalidation [61]. (PC-MRI) – Siemens Magnetom. Further details regard-ing the flow acquisition can be found in [62]. To ensurea fully developed laminar flow profile the inlet connectorto the phantom consisted of a 500 mm straight tubing.The PIV double frame recordings were conducted at a fre-quency of 500 Hz (triggered by the pump control), wherethe interframe time was set to 200 µ s. In total, 36 peri-odic cycles were recorded, resulting in 36,000 double framepairs. Velocity processing was conducted via DaVis 8.4.0(LaVision, G¨ottingen, Germany) using a multi-pass stereocross-correlation with a final interrogation window size of32 ×
32 px and 50 % overlap, resulting in one velocityvector every 141 µ m. The velocity fields were then phase-averaged. Additional details on the PIV data acquisitionand processing procedure can be found in [60]. Figure 2: Average velocity profile at the inlet over one cardiac cyclefor the idealized pulsatile aneurysm.
A pulsatile hemodynamic simulation was carried outusing identical conditions compared to the experiment.Specifically, the same geometry and boundary conditionswere applied, with a higher temporal and spatial resolu-tion. The LBM simulation was performed using time-stepand grid sizes of δ t = 1 . × − s and δ x = 1 . × − m,respectively. Two additional lower resolution simulationswere also conducted to provide qualitative/quantitativecharacterizations of convergence and under-resolution ef-fects. Furthermore, the simulations were ran for multiplecycles and data sampling was carried out at the fourth cy-cle. The state of the flow, and the absence of any artifactstied to initial conditions were assessed by monitoring thevelocity changes at three points around the aneurysm sac.The obtained results (from the highest resolution simula-tion) are shown in Fig. 3. In order to increase the anatomical complexity and en-able a qualitative and quantitative comparison in a morerealistic scenario, a patient-specific aneurysm was furtherconsidered (see Fig. 1). Here, an internal carotid arteryaneurysm with a spherical shape was chosen, since thelocation represents one of the most frequent sites of oc-currence. Further, a small aneurysm is located oppositethe larger one leading to the formation of interesting flowstructures.To account for the most comprehensive situation from ahemodynamic perspective, peak systolic flow conditionswere chosen for the validation. Specifically, a flow rate of430 mL/min was defined resulting in a Reynolds numberof Re peak = 563 and a corresponding mean velocity in theproximal vessel of 0 .
512 m/s, respectively.Equivalent to the experiments described in the previous5
Figure 3: Velocity magnitude changes over four cycles at three dif-ferent monitoring points around the aneurysm sac for the idealizedpulsatile aneurysm. section, stereoscopic PIV measurements were carried outto obtain in-vitro flow data. Specifically, the micro-gearpump (HNP Mikrosysteme, Schwerin, Germany) providedthe desired flow rate of Q = 430 . ± . ×
32 px (correspond-ing to 208 × µ m) and 50 % overlap. This resulted ina velocity vector every 104 µ m. The final results wereobtained by averaging all processed recordings.In addition to the PIV acquisitions, phase-contrast PC-MRI measurements were carried out in a 7 Tesla whole-body system (Siemens Healthineers, Forchheim, Germany)with a voxel size of 570 × × µ m . The measurementswere repeated six times and afterwards the data were av-eraged. Before evaluation, eddy currents and other back-ground phase effects were corrected using a second datasetwith identical scan parameters. For further details regard-ing the scan parameters and the data processing the inter-ested reader is referred to [17, 63].For the corresponding LBM simulation, similar to theidealized geometry, blood is assumed to be a Newtonianfluid with density ρ = 1221 kg/m and kinematic viscosity ν = 3 . × − m / s, a valid choice since shear-thinningeffect in the considered neurovasculature are negligible.The time-step and grid sizes were set to δ t = 5 × − sand δ x = 1 × − m, respectively. This choice of pa-rameters led to a non-dimensional relaxation coefficient of τδ t = 0 .
506 and a maximum non-dimensional velocity of u max δ t δ x ≈ . δ t = 1 × − s and δ x = 2 × − m, leading to τδ t = 0 .
3. Results
In the following, qualitative and quantitative compar-isons between the hemodynamic simulations based on theLB approach and the in-vitro phantom measurement arepresented.
As a first step, before going into detailed validation ofthe solver against experimental data, we present a briefcomparative study of the solver (against a classical SRTsolver) through a simple test-case. The steady-state ve-locity profiles along the incoming and outgoing vessels,as obtained from the different simulations are shown inFig. 4. The simulations using the classical SRT modelwere conducted using three different resolutions, i.e. δ x =3 . × − , 2 . × − and 1 . × − m (R1, R2, andR3 in Fig. 4). For the CHMRT only one simulation at δ x = 3 . × − was conducted. It can clearly be ob-served that at the lowest resolution, the CHMRT solvermatches results from the SRT simulation with the highestresolution, and even surpasses it in some regions. The low-resolution SRT simulation exhibits clear underestimationof the maximum velocity at the center of the inflow vessel.The CHMRT simulation on the other hand, even throughrelying on a grid-size three times that of the highest res-olution SRT, is the closest to the reference STAR-CCM+simulation, showing the superior performances of the for-mer. It is worth noting that the differences observed nearthe inlet are to be expected, due the different ways bound-ary conditions are applied. The length of the inlet pipe, asobserved in the velocity distribution plots, guarantees thatat the inlet of the sac the velocity profiles are establishedfor all simulations. The angle of the model for the second case was chosenin a way that the flow can enter into the sac through thedistal part of the aneurysm ostium. After impinging onthe aneurysm wall, the flow aligns along the perfect spher-ical shape. This forms a large vortex within the aneurysmwith a stagnation zone approximately at the center of thesac.The comparison between the experimental (PIV) and thenumerical (LBM) flow acquisition demonstrates an excel-lent agreement within the considered plane (perpendicularto the z -axis passing through the center of the aneurysmsac). Results for 11 different times spanning a cycle areshown in Fig. 5. Almost identical flow structures are cap-tured and slight deviations are only visible in the out-flow region of the aneurysm close to the peak-systolic time6 [m/s]0.190 xz y x[m] u [ m / s ] Star-CCMSRT-R1SRT-R2SRT-R3CH-MRT y[m] u [ m / s ] Star-CCMSRT-R1SRT-R2SRT-R3CH-MRT
Figure 4: (from left to right) Geometrical configuration of the steady flow in the idealized aneurysm, velocity distribution in the x − direction,and velocity distribution in the y − direction (along the dashed lines shown in the left subfigure). point (0.4 s). Additionally, minor differences can be ob-served most notably in the form of, less pronounced, un-steady structures not observed in the experimental data.The absence of these structure can partially be explainedby the data acquisition modes in experiments and simula-tions. While the time-dependent flow field from the PIVare phase-averaged over a number of cycles, those from thesimulations represent instantaneous shots from the fourthcycle after initialization. In addition to this qualitative ob-servation, a quantitative comparison is presented in Fig. 6.Here, the velocity components ( u x and u y ) are comparedalong two orthogonal lines on the plane perpendicular tothe z -axis for three different resolutions, confirming theinitial findings: Both experiment and simulation lead toalmost identical profiles throughout the cardiac cycle andthe numerical approach nearly always occurs within thepresented error bars of the measurement. It is worth not-ing that two additional simulations were performed by set-ting δ x to 2 × − m and 3 × − m; the correspond-ing time-steps were set to respectively 1 . × − s and2 . × − s. At the lowest resolution, as compared tothe highest resolution simulation, the computational costis reduced by a factor of 25, while velocity profiles arestill in very good agreement (and mostly within the errorbars) with the experimental data. Overall, based on bothqualitative and quantitative comparisons with PIV data,it is shown that the solver is able to correctly capture theflow structure and match experimental observations. Thegood agreement is still well maintained even for relativelyunder-resolved simulations. The last part of the validation was carried out in apatient-specific aneurysm model enabling a clearly morerealistic flow scenario. The obtained velocity fields (bothnumerical and experimental) on three perpendicular planescutting through the main aneurysm sac are shown in Fig. 7.It must be noted that both PIV and PC-MRI measure-ments were conducted here (as opposed to the previous configuration where only PIV measurements were made).Similar to the idealized test-case, a vortex forms withinthe sac aligning along the luminal surface and forming astagnation zone at the center. Due to the helical vesselanatomy proximal to the aneurysm as well as the differentostium size, this vortex appears less stable.Regarding the qualitative comparison of the three indepen-dent methods of flow acquisition one can notice that theoverall flow structure is captured by all approaches. Theentering flow jet of the simulation appears to be slightlynarrower and contains accordingly higher velocity valuescompared to the experimental techniques. A look at theconvergence indicator in the simulations and the flow fieldevolution over time showed that this configuration resultsin an unsteady flow although the inlet velocity is not pul-satile. This numerical observation was also confirmed bystudying PIV instantaneous velocity fields. As such thedata averaging procedure in time, i.e. sampling frequencyand time interval, can have non-negligible effect on out-coming velocity fields from both experiments and numer-ical simulations. The third column in Fig. 7, representingthe instantaneous fields from the LBM simulation are in-cluded to point out the unsteady nature of the flow. Itmust be noted that due to the computational costs andlarge sizes of associated files, an interval of 0.25 s andsampling frequency of 400 Hz were used to get the averagevelocity fields for the simulation while as noted earlier, forthe PIV the data acquisition was operated at a frequencyof 5 Hz and spanned 100 s. For the PC-MRI results on theother hand, the final velocity fields are the average of sixsnapshots taken at a repetition/echo time of 6.9/3.586 ms.These differences in the averaging process stemming fromlimitations of each one of these approaches can explain,to a great extent, the discrepancies observed in Fig. 7. Inpractice, the low data sampling frequency in PIV and PC-MRI measurements operates as a low-pass filter while thelimited sampling interval of the simulation is equivalent toa high-pass filter. This in turn can explain the thickerand smoother shear layers observed in the former two.7 y x 𝑢 [m/s] P I V L B M A B C D E F G H I J K t[s] A B C D E
F G
H I J K 𝑢 [ m / s ] Figure 5: Qualitative comparison of the transient flow field between PIV experiments (top row) and LB simulation (bottom row) at elevendifferent time points covering an entire cardiac cycle (from 0 s to 1 s with ∆ t =0.1 s). The considered plane is shown on the upper left figure.The positions of the snapshots in time (relative to the cardiac cycle) are shown with red dashed lines on the upper right plot. The subsequent quantification presented in Fig. 8 providesa better insight in occurring differences and further con-firms the previously discussed effect regarding the discrep-ancies coming from data acquisition. Here, velocity mag-nitudes are shown along three orthogonal lines through theaneurysm. Overall, the courses of all approaches are in agood agreement and the values show no strong deviationsfrom each other. It must be noted that error bars in Fig. 8only represent the 99.7% confidence interval. The simula-tion results fall within the 99.7% confidence interval at allpoints.Matching coordinates and positions of the correspondingvelocity fields, together with possible distortions in spaceresulting from the different experimental procedures, con-stitute also a major challenge for this comparison; it cer-tainly explains to some extent the observed discrepancy.This effect is particularly visible on the line along the y -axis in Fig. 8. While the velocity profiles from the sim-ulation tend towards zero near the aneurysm sac walls,PIV results do not show such a clear behavior. Interest-ingly, the agreement between all methods is excellent forplane III (see Fig. 8, right), while larger discrepancies arenoticeable in the other two (e.g., smoother courses of themeasurement compared to the simulation). Overall, whilea good qualitative agreement between the three differentapproaches can be observed, non-negligible discrepanciespersist. In-depth analysis of the results showed that theunsteady nature of the configuration, along with the dif-ferent data acquisition modes (i.e. frequency and time-span of velocity field averaging process) affect the finalvelocity fields. This effect was illustrated by comparinginstantaneous and averaged velocity fields in the simula-tion. The distortions and manufacturing artifacts in thefinal aneurysm geometry can also contribute to these dis-crepancies. Further reasons are discussed in Section 4. Finally, to showcase the ability of the collision operatorto deal with under-resolved cases, and provide a more in-depth analysis, the patient-specific configuration was alsomodeled using a coarse grid. The time-averaged velocityfields obtained using the high- and the low-resolution LBMsimulation are shown in Fig. 9. The orthogonal planes re-veal that the predicted flow structures, while admittingdiscrepancies, agree very well. The quantitative compar-ison of both simulations leads to a domain-averaged sim-ilarity of 0.8294, 0.8765 and 0.8359 for planes I, II andIII, respectively. This similarity index on each plane iscomputed as the average similarity index over all discretepoints on the plane [60]:SI = (cid:18) u HR · u LR || u HR |||| u LR || (cid:19)(cid:18) − | || u HR |||| u HR || max − || u LR |||| u LR || max | (cid:19) , (11)where u HR and u LR are the velocity vectors from thehigh- and low-resolution simulations, respectively. Thedeviations are further illustrated in Fig. 10, where thelow-resolution velocity magnitudes are plotted against thehigh-resolution ones on the three considered planes. Thenormalized frequency of the similarity index on all threeplanes is shown in Fig. 11. Despite these minor differenceswith respect to the velocity distribution, it is importantto mention that the low-resolution simulation reduced thenumber of grid-points by a factor of eight, and – giventhe acoustic scaling – the number of time-steps needed forconvergence by a factor of two. This in turn reduced thecomputation time by a factor of 14 on the same machineand using the same number of processing units.8
10 20-0.8-0.20 -0.2 0 01020
Figure 6: Velocity profiles along the (in red) x - and (in blue) y -directions at the center of the aneurysm sac on the central z -plane at fourdifferent times: (from left to right) 0.16, 0.27, 0.51 and 0.61 s from the start of the fourth cycle, and (top to bottom) decreasing resolutions inthe lattice Boltzmann simulations. Velocity profiles obtained from the LBM simulation are shown with red and blue plain lines while blackplain lines with error bars (corresponding to the 99.7% confidence interval) represent data obtained from PIV measurements. z y x 𝑢 [m/s] IIIIII PC-MRI PIV LBM Averaged LBM
Figure 7: From left to right: Qualitative comparison of the velocity fields acquired using PC-MRI (1st column), stereoscopic PIV (2ndcolumn), instantaneous (3rd column) and time-averaged LBM (last column) simulation in the aneurysm sac on three different planes, i.e.from top to bottom: planes I, II and III. -3 -3 Figure 8: Comparison of velocity magnitude profiles along three lines in the aneurysm sac as obtained from PIV (red plain lines with errorbars corresponding to one standard deviation), PC-MRI (blue dashed lines with square symbols), instantaneous LB velocity field (black plainlines) and time-averaged (over 0.25 s) LB velocity field (grey plain lines). z y x 𝑢 [m/s] IIIII I Figure 9: Qualitative comparison of the velocity fields in theaneurysm sac on three different planes, i.e. from top to bottom:planes I, II and III, as obtained from the high-resolution (left) andlow-resolution (right) simulation.
4. Discussion
With an improved linkage of medical and engineeringdisciplines as well as increasing computational resourcesnew potentials in interdisciplinary research are revealed.Specifically, the challenging question of rupture risk assess-ment for intracranial aneurysm remains unanswered untilnow; however, more and more insights into patient-specifichemodynamics become accessible [64]. On the other hand,insufficient data processing, the absence of reliable bound-ary conditions and mostly long simulation times lead to areduced translation into clinical practice [65].To highlight the potential of image-based blood flow simu-lations and demonstrate its reliability with respect to flowprediction, a multimodal validation study was carried out.This included different configuration, flow scenarios, mea-surement methodologies and simulation settings, respec-
Figure 10: Quantitative comparison of the velocity magnitude dis-tribution in the high- and low-resolution simulations on planes I (topleft), II (top right) and III (bottom). The color bar represents thenormalized (via the total number of sampled points) number of grid-points with the corresponding velocities in the low-resolution ( x -axis)and high-resolution ( y -axis) simulation. tively.The first part of the validation relied on flow assessmentsunder well-controlled conditions. A brief comparative studyshowed that the CHMRT collision operator with a full ex-tension of the EDF performed much better than the classi-cal SRT model. To further validate the solver an idealizedmodel of an intracranial aneurysm was created leading tothe development of an organized flow field. The quali-tative and quantitative comparison of both independenttechniques (LBM versus PIV) showed an excellent agree-ment. Due to the unsteady inflow, minor fluctuations oc-curred throughout the cardiac cycle, especially near peak10 igure 11: Histograms displaying the normalized frequency of theSI on planes I (blue), II (red) and III (black). The correspondingoverall SIs are shown using dashed lines of the same color. systole. Still, both simulation and measurement were invery good agreement during the whole period.In a second part of this study, complexity was increasedby considering a patient-specific anatomy. Although un-steady phenomena occurred at peak-systolic inflow con-ditions, flow predictions by simulation and measurementwere comparable. Beside the validation using PIV, PC-MRI acquisitions were carried out. Here, minor differencesespecially in high-velocity regions close to the aneurysmallumen were present, which could be a result of vessel wallcompliance [66]. However, taking into account all under-lying differences with respect to temporal and spatial res-olution as well as the uncertainties associated to the mea-surement procedures, the validation of the LBM simula-tion was successfully conducted.Apart from the validation purpose, additional simulationswith low spatial resolutions were performed. The compar-ison of low- and high-resolution LBM computation showedthat similar velocity fields were calculated, although only1 /
14 of the simulation time was required. Hence, depend-ing on the clinical research question and the desired ac-curacy of the hemodynamic description the low-resolutionapproach, made possible by the CHMRT collision opera-tor, can be a promising development direction, especiallygiven that the computational overhead associated to thechange of collision operator is very minor [67].Beside the presented findings, it is important to mentionthat this study has several limitations: First, the numberof cases considered is relatively small for an extensive val-idation study. Only two idealized and one patient-specificgeometries were selected. However, it is important to pointout that intracranial aneurysms show a high variability re-garding size (small versus giant), location (lateral versusterminal), morphology (spherical versus complex). Conse-quently, this can lead to increased challenges for the under-lying numerical schemes. Nevertheless, multiple numericaland experimental processing steps are involved.Second, precise qualitative and quantitative comparisonsare only feasible, when identical conditions are warranted. However, due to the manufacturing of the phantom mod-els for the in-vitro investigations or the subsequent regis-trations processes, minor misalignment between the idealCAD geometries and the silicone models can occur.Third, blood was considered as an incompressible New-tonian fluid in the simulations. There is an ongoing de-bate about the necessity of a more detailed considerationof this suspension [68, 69, 70]. Nevertheless, blood rheol-ogy still appears to have a secondary impact compared toprimary factors such as segmentation or boundary condi-tions [71, 8]. Finally, although PC-MRI is a technique toobtain flow data in-vivo, within this study in-vitro mea-surements in a silicone phantom were carried out. There-fore, the MR sequence might lead to slightly different re-sults when applied to a real aneurysm.Since this study demonstrates the feasibility of applyingthe introduced LBM-based solver to clinically relevant re-search problems, future work includes the considerationof an increased number of aneurysm patients. Further-more, specific questions with respect to rupture risk as-sessment, thrombus formation or treatment support willbe addressed, while accompanying experimental validationmeasurements are carried out to ensure the reliability ofour numerical models.
5. Conclusions
This validation study comprising multimodal measure-ment techniques and various experimental scenarios demon-strates that the presented LBM solver relying on a CHMRTcollision operator is able to generate reliable and valid sim-ulation results. These results are in line with observationsmade in the LBM literature on the effect of such modelsin eliminating a number of numerical artefacts associatedto the classical SRT collision operator with second-orderEDF while maintaining a low computational cost low –around 0.4 Mega lattice updates per second (MLUPS) perprocessor. As shown in [47] the added cost of the modifiedcollision operator is negligible. Furthermore, the widerstability domain and enhanced spectral properties of thecollision model (as compared to the SRT) allowed for sim-ulations relying on lower resolutions (not stable with theSRT) yielding comparable numerical predictions with anextensive reduction of computational load (grid-points andtime-steps). Given the difficulties of applying explicit sub-grid model-based large eddy simulations (LES) to pulsatileflows in complex geometries, tools and numerical formula-tions allowing for model-free under-resolved simulations –tantamount to implicit LES – can be an interesting pathtowards a rapid evaluation of flow parameters in complexbiomedical settings.
Acknowledgements
S.A.H. would like to acknowledge the financial supportof the Deutsche Forschungsgemeinschaft (DFG, German11esearch Foundation) in TRR 287 (Project-ID 422037413).F.H. was supported by the State Scholarship Fund of theChina Scholarship Council (grant number 201908080236).Furthermore, P.B. acknowledges the funding by the Fed-eral Ministry of Education and Research within the Re-search Campus
STIMULATE (13GW0473A) and the Ger-man Research Foundation (BE 6230/2-1).The authors further thank Franziska Gaidzik and DanielStucht (both University of Magdeburg, Germany) for theirassistance during the PC-MRI flow measurements.
References [1] M. H. Vlak, A. Algra, R. Brandenburg, G. J. Rinkel, Prevalenceof unruptured intracranial aneurysms, with emphasis on sex,age, comorbidity, country, and time period: a systematic reviewand meta-analysis, The Lancet. Neurology 10 (7) (2011) 626–636.[2] L. Liang, D. A. Steinman, O. Brina, C. Chnafa, N. M. Can-celliere, V. M. Pereira, Towards the clinical utility of CFD forassessment of intracranial aneurysm rupture – a systematic re-view and novel parameter-ranking tool, Journal of Neurointer-ventional Surgery 11 (2) (2019) 153–158.[3] K. M. Saqr, S. Rashad, S. Tupin, K. Niizuma, T. Hassan,T. Tominaga, M. Ohta, What does computational fluid dynam-ics tell us about intracranial aneurysms? a meta-analysis andcritical review, Journal of Cerebral Blood Flow & Metabolism40 (5) (2020) 1021–1039.[4] H. Meng, V. M. Tutino, J. Xiang, A. Siddiqui, High WSS or lowWSS? complex interactions of hemodynamics with intracranialaneurysm initiation, growth, and rupture: toward a unifyinghypothesis, American Journal of Neuroradiology 35 (7) (2014)1254–1262.[5] J. Cebral, E. Ollikainen, B. J. Chung, F. Mut, V. Sippola, B. R.Jahromi, R. Tulamo, J. Hernesniemi, M. Niemel¨a, A. Robert-son, J. Fr¨osen, Flow conditions in the intracranial aneurysm lu-men are associated with inflammation and degenerative changesof the aneurysm wall, American Journal of Neuroradiology38 (1) (2017) 119–126.[6] D. Fiorella, C. Sadasivan, H. H. Woo, B. Lieber, Regard-ing “Aneurysm rupture following treatment with flow-divertingstents: computational hemodynamics analysis of treatment”,American Journal of Neuroradiology 32 (5) (2011) E95–7; au-thor reply E98–100.[7] D. F. Kallmes, Point: CFD–computational fluid dynamics orconfounding factor dissemination, American Journal of Neuro-radiology 33 (3) (2012) 395–396.[8] P. Berg, S. Saalfeld, S. Voß, O. Beuing, G. Janiga, A re-view on the reliability of hemodynamic modeling in intracra-nial aneurysms: why computational fluid dynamics alone can-not solve the equation, Neurosurgical Focus 47 (1) (2019) E15.[9] M. D. Ford, H. N. Nikolov, J. S. Milner, S. P. Lownie, E. M. De-mont, W. Kalata, F. Loth, D. W. Holdsworth, D. A. Steinman,PIV-measured versus CFD-predicted flow dynamics in anatomi-cally realistic cerebral aneurysm models, Journal of Biomechan-ical Engineering 130 (2) (2008) 021015.[10] M. Raschi, F. Mut, G. Byrne, C. M. Putman, S. Tateshima,F. Vi˜nuela, T. Tanoue, K. Tanishita, J. R. Cebral, CFDand PIV analysis of hemodynamics in a growing intracra-nial aneurysm, International Journal for Numerical Methodsin Biomedical Engineering 28 (2) (2012) 214–228.[11] T. Yagi, A. Sato, M. Shinke, S. Takahashi, Y. Tobe, H. Takao,Y. Murayama, M. Umezu, Experimental insights into flow im-pingement in cerebral aneurysm by stereoscopic particle imagevelocimetry: transition from a laminar regime, Journal of theRoyal Society, Interface 10 (82) (2013) 20121031.[12] P. Bouillot, O. Brina, R. Ouared, K.-O. Lovblad, M. Farhat,V. M. Pereira, Particle imaging velocimetry evaluation of in- tracranial stents in sidewall aneurysm: hemodynamic transitionrelated to the stent design, PloS one 9 (12) (2014) e113762.[13] P. Berg, C. Roloff, O. Beuing, S. Voss, S.-I. Sugiyama, N. Aris-tokleous, A. S. Anayiotos, N. Ashton, A. Revell, N. W. Bressloff,A. G. Brown, B. J. Chung, J. R. Cebral, G. Copelli, W. Fu,A. Qiao, A. J. Geers, S. Hodis, D. Dragomir-Daescu, E. Nor-dahl, Y. Bora Suzen, M. Owais Khan, K. Valen-Sendstad,K. Kono, P. G. Menon, P. G. Albal, O. Mierka, R. M¨unster,H. G. Morales, O. Bonnefous, J. Osman, L. Goubergrits, J. Pal-lares, S. Cito, A. Passalacqua, S. Piskin, K. Pekkan, S. Ra-malho, N. Marques, S. Sanchi, K. R. Schumacher, J. Sturgeon,H. ˇSvihlov´a, J. Hron, G. Usera, M. Mendina, J. Xiang, H. Meng,D. A. Steinman, G. Janiga, The computational fluid dynamicsrupture challenge 2013–phase II: variability of hemodynamicsimulations in two intracranial aneurysms, Journal of Biome-chanical Engineering 137 (12) (2015) 121008.[14] P. Van Ooij, J. Schneiders, H. Marquering, C. Majoie,E. Van Bavel, A. Nederveen, 3-D cine phase-contrast MRI at 3Tin intracranial aneurysms compared with patient-specific com-putational fluid dynamics, American Journal of Neuroradiology34 (9) (2013) 1785–1791.[15] J. R. Cebral, R. S. Pergolizzi Jr, C. M. Putman, Computationalfluid dynamics modeling of intracranial aneurysms:: Qualita-tive comparison with cerebral angiography, Academic Radiology14 (7) (2007) 804–813.[16] N. Paliwal, R. J. Damiano, N. A. Varble, V. M. Tutino, Z. Dou,A. H. Siddiqui, H. Meng, Methodology for computational fluiddynamic validation for medical use: application to intracra-nial aneurysm, Journal of Biomechanical Engineering 139 (12)(2017).[17] C. Roloff, D. Stucht, O. Beuing, P. Berg, Comparison of in-tracranial aneurysm flow quantification techniques: standardPIV vs stereoscopic PIV vs tomographic PIV vs phase-contrastMRI vs CFD, Journal of Neurointerventional Surgery 11 (3)(2019) 275–282.[18] M. M¨uller, S. Schirm, M. Teschner, Interactive blood simulationfor virtual surgery based on smoothed particle hydrodynamics,Technology and Health Care 12 (1) (2004) 25–31.[19] J. Qin, W.-M. Pang, B. P. Nguyen, D. Ni, C.-K. Chui, Particle-based simulation of blood flow and vessel wall interactions invirtual surgery, in: Proceedings of the 2010 Symposium on In-formation and Communication Technology, 2010, pp. 128–133.[20] G. Z´avodszky, G. Pa´al, Validation of a lattice Boltzmannmethod implementation for a 3D transient fluid flow in an in-tracranial aneurysm geometry, International Journal of Heatand Fluid Flow 44 (2013) 276–283.[21] M. D. Mazzeo, P. V. Coveney, HemeLB: A high performanceparallel lattice-Boltzmann code for large scale fluid flow in com-plex geometries, Computer Physics Communications 178 (12)(2008) 894–914.[22] D. Groen, J. Hetherington, H. B. Carver, R. W. Nash, M. O.Bernabeu, P. V. Coveney, Analysing and modelling the per-formance of the hemeLB lattice-Boltzmann simulation environ-ment, Journal of Computational Science 4 (5) (2013) 412–422.[23] H. Anzai, M. Ohta, J.-L. Falcone, B. Chopard, Optimization offlow diverters for cerebral aneurysms, Journal of ComputationalScience 3 (1-2) (2012) 1–7.[24] C. Huang, B. Shi, Z. Guo, Z. Chai, Multi-GPU based lat-tice Boltzmann method for hemodynamic simulation in patient-specific cerebral aneurysm, Communications in ComputationalPhysics 17 (4) (2015) 960–974.[25] R. Ouared, B. Chopard, Lattice Boltzmann simulations of bloodflow: non-Newtonian rheology and clotting processes, Journalof Statistical Physics 121 (1-2) (2005) 209–221.[26] M. Finck, D. H¨anel, I. Wlokas, Simulation of nasal flow by lat-tice Boltzmann methods, Computers in Biology and Medicine37 (6) (2007) 739–749.[27] Y. Chen, L. Navarro, Y. Wang, G. Courbebaisse, Segmentationof the thrombus of giant intracranial aneurysms from CT an-giography scans with lattice Boltzmann method, Medical ImageAnalysis 18 (1) (2014) 1–8.
28] Y. Ma, R. Mohebbi, M. Rashidi, Z. Yang, Numerical simula-tion of flow over a square cylinder with upstream and down-stream circular bar using lattice Boltzmann method, Interna-tional Journal of Modern Physics C 29 (04) (2018) 1850030.[29] Y. Ma, R. Mohebbi, M. Rashidi, Z. Yang, M. A. Sheremet,Numerical study of mhd nanofluid natural convection in a baf-fled u-shaped enclosure, International Journal of Heat and MassTransfer 130 (2019) 123–134.[30] T. Kr¨uger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva,E. M. Viggen, The lattice Boltzmann method, Springer Inter-national Publishing 10 (2017) 978–3.[31] J. Latt, O. Malaspinas, D. Kontaxakis, A. Parmigiani, D. La-grava, F. Brogi, M. B. Belgacem, Y. Thorimbert, S. Leclaire,S. Li, et al., Palabos: Parallel lattice Boltzmann solver, Com-puters & Mathematics with Applications (2020).[32] M. Hasert, K. Masilamani, S. Zimny, H. Klimach, J. Qi,J. Bernsdorf, S. Roller, Complex fluid simulations with theparallel tree-based lattice Boltzmann solver musubi, Journal ofComputational Science 5 (5) (2014) 784–794.[33] T. Suzuki, C. Ioan Nita, S. Rapaka, H. Takao, V. Mi-halef, S. Fujimura, C. Dahmani, P. Sharma, H. Mamori,T. Ishibashi, T. Redel, M. Yamamoto, Y. Murayama, Verifica-tion of a research prototype for hemodynamic analysis of cere-bral aneurysms, Conference proceedings : Annual InternationalConference of the IEEE Engineering in Medicine and BiologySociety. IEEE Engineering in Medicine and Biology Society. An-nual Conference 2016 (2016) 2921–2924.[34] P. Berg, S. Saalfeld, S. Voß, T. Redel, B. Preim, G. Janiga,O. Beuing, Does the dsa reconstruction kernel affect hemody-namic predictions in intracranial aneurysms? An analysis ofgeometry and blood flow variations, Journal of Neurointerven-tional Surgery 10 (3) (2018) 290–296.[35] T. Suzuki, H. Takao, S. Rapaka, S. Fujimura, C. Ioan Nita,Y. Uchiyama, H. Ohno, K. Otani, C. Dahmani, V. Mihalef,P. Sharma, A. Mohamed, T. Redel, T. Ishibashi, M. Yamamoto,Y. Murayama, Rupture risk of small unruptured intracranialaneurysms in japanese adults, Stroke 51 (2) (2020) 641–643.[36] C. Pan, L. S. Luo, C. T. Miller, An evaluation of lattice Boltz-mann schemes for porous medium flow simulation, Computers& Fluids 35 (8-9) (2006) 898–909.[37] S. A. Hosseini, H. Safari, N. Darabiha, D. Th´evenin,M. Krafczyk, Hybrid lattice Boltzmann-finite difference modelfor low Mach number combustion simulation, Combustion andFlame 209 (2019) 394–404.[38] S. A. Hosseini, N. Darabiha, D. Th´evenin, Theoretical and nu-merical analysis of the lattice kinetic scheme for complex flowsimulations, Physical Review E 99 (2) (2019) 023305.[39] S. Succi, The lattice Boltzmann equation: for fluid dynamicsand beyond, Oxford University Press, 2001.[40] P. L. Bhatnagar, E. P. Gross, M. Krook, A model for collisionprocesses in gases. I. Small amplitude processes in charged andneutral one-component systems, Physical Review 94 (3) (1954)511.[41] X. He, L.-S. Luo, Theory of the lattice Boltzmann method:From the Boltzmann equation to the lattice Boltzmann equa-tion, Physical Review E 56 (6) (1997) 6811.[42] X. Shan, X.-F. Yuan, H. Chen, Kinetic theory representationof hydrodynamics: a way beyond the Navier–Stokes equation,Journal of Fluid Mechanics 550 (2006) 413–441.[43] S. A. Hosseini, C. Coreixas, N. Darabiha, D. Th´evenin, Stabilityof the lattice kinetic scheme and choice of the free relaxationparameter, Physical Review E 99 (6) (2019) 063305.[44] N. I. Prasianakis, I. V. Karlin, Lattice Boltzmann method forsimulation of compressible flows on standard lattices, PhysicalReview E 78 (1) (2008) 016704.[45] I. Karlin, P. Asinari, Factorization symmetry in the latticeBoltzmann method, Physica A: Statistical Mechanics and itsApplications 389 (8) (2010) 1530–1548.[46] S. A. Hosseini, N. Darabiha, D. Th´evenin, Compressibility inlattice Boltzmann on standard stencils: effects of deviation fromreference temperature, Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engi-neering Sciences 378 (2175) (2020) 20190399.[47] S. A. Hosseini, Development of a lattice Boltzmann-based nu-merical method for the simulation of reacting flows, Ph.D.thesis, Universit´e Paris-Saclay/Otto von Guericke University(2020).[48] M. H. Saadat, F. B¨osch, I. V. Karlin, Lattice Boltzmann modelfor compressible flows on standard lattices: Variable prandtlnumber and adiabatic exponent, Physical Review E 99 (1)(2019) 013306.[49] I. Ginzburg, F. Verhaeghe, D. d’Humieres, Two-relaxation-timelattice Boltzmann scheme: About parametrization, velocity,pressure and mixed boundary conditions, Communications inComputational Physics 3 (2) (2008) 427–478.[50] F. J. Higuera, S. Succi, R. Benzi, Lattice gas dynamics withenhanced collisions, EPL (Europhysics Letters) 9 (4) (1989) 345.[51] D. d’Humieres, Generalized lattice Boltzmann equations, Rar-efied Gas Dynamics (1992) 450–458.[52] P. Lallemand, L. S. Luo, Theory of the lattice Boltzmannmethod: Dispersion, dissipation, isotropy, Galilean invariance,and stability, Physical Review E 61 (6) (2000) 6546.[53] D. d’Humieres, Multiple–relaxation–time lattice Boltzmannmodels in three dimensions, Philosophical Transactions of theRoyal Society of London. Series A: Mathematical, Physical andEngineering Sciences 360 (1792) (2002) 437–451.[54] M. Geier, A. Greiner, J. G. Korvink, Cascaded digital latticeBoltzmann automata for high Reynolds number flow, PhysicalReview E 73 (6) (2006) 066705.[55] J. Latt, B. Chopard, Lattice Boltzmann method with regular-ized pre-collision distribution functions, Mathematics and Com-puters in Simulation 72 (2-6) (2006) 165–168.[56] S. Ansumali, I. V. Karlin, Single relaxation time model for en-tropic lattice Boltzmann methods, Physical Review E 65 (5)(2002) 056312.[57] I. V. Karlin, F. B¨osch, S. S. Chikatamarla, Gibbs’ principle forthe lattice-kinetic theory of fluid dynamics, Physical Review E90 (3) (2014) 031302.[58] F. B¨osch, S. S. Chikatamarla, I. V. Karlin, Entropic multire-laxation lattice Boltzmann models for turbulent flows, PhysicalReview E 92 (4) (2015) 043309.[59] M. Bouzidi, M. Firdaouss, P. Lallemand, Momentum transferof a Boltzmann-lattice fluid with boundaries, Physics of Fluids13 (11) (2001) 3452–3459.[60] F. Gaidzik, D. Stucht, C. Roloff, O. Speck, D. Th´evenin,G. Janiga, Transient flow prediction in an idealized aneurysmgeometry using data assimilation, Computers in Biology andMedicine 115 (2019) 103507.[61] C. Roloff, D. Stucht, O. Beuing, P. Berg, Comparison of in-tracranial aneurysm flow quantification techniques: standardPIV vs stereoscopic PIV vs tomographic PIV vs phase-contrastMRI vs CFD, Journal of Neurointerventional Surgery 11 (3)(2019) 275–282.[62] P. Berg, D. Stucht, G. Janiga, O. Beuing, O. Speck,D. Th´evenin, Cerebral blood flow in a healthy Circle of Willisand two intracranial aneurysms: computational fluid dynam-ics versus four-dimensional phase-contrast magnetic resonanceimaging, Journal of Biomechanical Engineering 136 (4) (2014).[63] C. Roloff, P. Berg, T. Redel, G. Janiga, D. Thevenin, Tomo-graphic particle image velocimetry for the validation of hemo-dynamic simulations in an intracranial aneurysm, Conferenceproceedings Annual International Conference of the IEEE En-gineering in Medicine and Biology Society. IEEE Engineering inMedicine and Biology Society. 2017 (2017) 1340–1343.[64] Y. Murayama, S. Fujimura, T. Suzuki, H. Takao, Computa-tional fluid dynamics as a risk assessment tool for aneurysmrupture, Neurosurgical Focus 47 (1) (2019) E12.[65] P. Berg, S. Vos, M. Becker, S. Serowy, T. Redel, G. Janiga,M. Skalej, O. Beuing, Bringing hemodynamic simulationscloser to the clinics: a CFD prototype study for intracranialaneurysms, Conference proceedings: Annual International Con-ference of the IEEE Engineering in Medicine and Biology Soci- ty. 2016 (2016) 3302–3305.[66] S. Tupin, K. M. Saqr, M. Ohta, Effects of wall complianceon multiharmonic pulsatile flow in idealized cerebral aneurysmmodels: comparative PIV experiments, Experiments in Fluids61 (7) (2020) 164.[67] C. Coreixas, B. Chopard, J. Latt, Comprehensive comparisonof collision models in the lattice Boltzmann framework: Theo-retical investigations, Physical Review E 100 (3) (2019) 033305.[68] K. M. Saqr, O. Mansour, S. Tupin, T. Hassan, M. Ohta, Ev-idence for non-Newtonian behavior of intracranial blood flowfrom doppler ultrasonography measurements, Medical & Bio-logical Engineering & Computing 57 (5) (2019) 1029–1036.[69] U. Y. Lee, J. Jung, H. S. Kwak, D. H. Lee, G. H. Chung, J. S.Park, E. J. Koh, Patient-specific computational fluid dynam-ics in ruptured posterior communicating aneurysms using mea-sured non-Newtonian viscosity : A preliminary study, Journalof Korean Neurosurgical Society 62 (2) (2019) 183–192.[70] S. Mahrous, N. A. C. Sidik, K. M. Saqr, Newtonian and non-Newtonian CFD models of intracranial aneurysm: A review,CFD Letters 12 (1) (2020) 62–86.[71] M. O. Khan, D. A. Steinman, K. Valen-Sendstad, Non-Newtonian versus numerical rheology: Practical impact ofshear-thinning on the prediction of stable and unstable flowsin intracranial aneurysms, International Journal for NumericalMethods in Biomedical Engineering 33 (7) (2017).ty. 2016 (2016) 3302–3305.[66] S. Tupin, K. M. Saqr, M. Ohta, Effects of wall complianceon multiharmonic pulsatile flow in idealized cerebral aneurysmmodels: comparative PIV experiments, Experiments in Fluids61 (7) (2020) 164.[67] C. Coreixas, B. Chopard, J. Latt, Comprehensive comparisonof collision models in the lattice Boltzmann framework: Theo-retical investigations, Physical Review E 100 (3) (2019) 033305.[68] K. M. Saqr, O. Mansour, S. Tupin, T. Hassan, M. Ohta, Ev-idence for non-Newtonian behavior of intracranial blood flowfrom doppler ultrasonography measurements, Medical & Bio-logical Engineering & Computing 57 (5) (2019) 1029–1036.[69] U. Y. Lee, J. Jung, H. S. Kwak, D. H. Lee, G. H. Chung, J. S.Park, E. J. Koh, Patient-specific computational fluid dynam-ics in ruptured posterior communicating aneurysms using mea-sured non-Newtonian viscosity : A preliminary study, Journalof Korean Neurosurgical Society 62 (2) (2019) 183–192.[70] S. Mahrous, N. A. C. Sidik, K. M. Saqr, Newtonian and non-Newtonian CFD models of intracranial aneurysm: A review,CFD Letters 12 (1) (2020) 62–86.[71] M. O. Khan, D. A. Steinman, K. Valen-Sendstad, Non-Newtonian versus numerical rheology: Practical impact ofshear-thinning on the prediction of stable and unstable flowsin intracranial aneurysms, International Journal for NumericalMethods in Biomedical Engineering 33 (7) (2017).