Central moments multiple relaxation time LBM for hemodynamic simulations in intracranial aneurysms: An in-vitro validation study using PIV and PC-MRI

Highlights

• Lattice Boltzmann simulation of intracranial aneurysm flows.

• Use of Galilean invariant collision operator based on central Hermite moments.

• Cross-validation with stereoscopic particle image velocimetry and phase-contrast magnetic resonance imaging.

• Under-resolved simulations via extended numerical properties of collision operator.

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. with the so-called single relaxation time (SRT) collision operator, it has been observed to have a limited stability domain in the Courant/Fourier space, strongly constraining the minimum time-step and grid size. The development of improved collision models such as the multiple relaxation time (MRT) operator in central moments space has tremendously widened the stability domain, while allowing to overcome a number of other well-documented artifacts, therefore opening the door for simulations over a wider range of grid and time-step sizes. The present work focuses on implementing and validating a specific collision operator, the central Hermite moments multiple relaxation time model with the full expansion of the equilibrium distribution function, to simulate blood flows in intracranial aneurysms. The study further proceeds with a validation of the numerical model through different test-cases and against experimental measurements obtained via stereoscopic 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 were shown to be able to capture blood flow information with sufficient accuracy, as demonstrated through both qualitative and quantitative analysis of the flow field while leading to strongly reduced computation times. For instance in the case of the patient-specific configuration, increasing the grid-size by a factor of two led to a reduction of computation time by a factor of 14 with very good similarity indices still ranging from 0.83 to 0.88.

Introduction

Intracranial aneurysms (IA) are local malformations of the cerebral vasculature, which occur with an estimated prevalence of approximately 3% in the western population [1]. It is known that aneurysms can grow and develop into a stable state or rupture, if the hemodynamic forces exceed the vascular resistance. Since a reliable and clinically applicable rupture risk assessment procedure remains challenging until now, increasing research effort is being put into developing one. Owing to their non-invasive nature resulting in risk-free (for the patient) assessment of rupture probability and the possibility to access highly resolved velocity fields (in both space and time) computational fluid dynamics (CFD) studies are becoming increasingly popular [2,].

While valuable insights have been gained through such numerical studies [3,4], a number of challenges and gaps remain to be bridged, before predictive and reliable models based on CFD for aneurysm dynamics can be developed [5,6]. To put these challenges into perspective, Berg et al. [7] summarized the individual working steps involved in image-based blood flow simulations and provided corresponding recommendations to avoid simulation inaccuracies. As for any numerical simulation, the choice of the numerical solver, resolution and grid configuration are of the utmost importance to ensure convergence of the obtained solutions. The latter two also being consequences of the choice of the numerical method, the former is a determining factor concerning both convergence and efficiency of the solver.

The performances of the employed solvers can only be assessed through a systematic benchmarking/validation procedure against reliable experimental data. One of the most frequently used sources of experimental reference data is the particle image velocimetry (PIV) relying on laser-based high-speed camera flow measurements [[8], [9], [10], [11], [12]]. PIV measurements allow for more accurate measurements and higher resolutions as compared to other data acquisition tools such as magnetic resonance imaging (MRI) [13] or cerebral angiography [14]. A number of studies have presented qualitative comparisons of velocity fields for patient-specific aneurysms as obtained from CFD and PIV measurements, and reported good qualitative agreement between the two [15,16].

While most of the early publications relied on discrete solvers for the Navier-Stokes-Poisson (NSP) equations, e.g. using finite volume (FV) or finite element (FE) methods, emergent numerical methods such as smoothed particle hydrodynamics (SPH) [17,18] and the lattice Boltzmann method (LBM) [[19], [20], [21], [22], [23], [24], [25], [26]] are increasingly applied to such flows. Contrary to classical NSP solvers they are not strictly incompressible, as they approximate the low Mach flow regime through appropriate isothermal flow manifolds, and as such rely on a system of purely hyperbolic equations, making the algorithm inherently local. The absence of the elliptic component of the NSP equations allows for dramatic computation cost reduction while the parabolic nature of the partial differential equation (PDE) describing pressure evolution makes the formulation suitable for unsteady simulations – contrary to other weakly compressible formulations such as the artificial compressibility method (ACM). Furthermore, the compressible nature of the formulation (involving a thermodynamic pressure) along with kinetic heuristic boundary closures (such as the bounce-back rule) provide for an efficient and consistent implementation of wall boundaries [[27], [28], [29]]. The drastic decrease in computation time makes the method more viable compared to approaches based on the NSP equations. Based on these observations, a number of dedicated LBM-based solvers for blood flow simulations have been developed over the past decade, e.g. Refs. [20,30,31]. For example a research prototype was developed by the Siemens Healthineers AG (Erlangen, Germany), which is increasingly used to address clinical questions [[32], [33], [34]]. However, it must be noted that the LBM in its simplest form, i.e. with a single relaxation time (SRT) collision operator and a second-order discrete equilibrium distribution function (EDF) has a rather limited stability domain. The SRT formulation also leads to a number of other well-documented numerical artifacts including, but not limited to, the (non-dimensional) viscosity-dependence of the solid wall position when used with bounce-back-type boundary treatments [35].

The aim of the present work is to assess the performances of a specific class of multiple relaxation time (MRT) operator based on Hermite central moments and a fully expanded EDF. This collision model is shown to alleviate some of the traditional shortcomings of the classical SRT-based solvers, while – through its much wider stability domain – allowing for stable low-resolution simulations. The performances of the central Hermite multiple relaxation time (CHMRT) model are assessed, via our in-house solver ALBORZ [36,37], through a variety of test-cases spanning ideal and patient-specific geometries considering steady and pulsatile flows. The numerical results are compared with high-resolution in-vitro measurements for validation. They are also compared to results from a SRT solver with second-order EDF to further showcase the added value of the CHMRT collision operator. The issue of under-resolved simulations is also considered by systematically conducting simulations at different (lower) resolutions. It is shown that at low resolution (inaccessible to the SRT collision operator), the proposed scheme is still able to capture properly the flow dynamics.