Numerical Simulations of Flows in a Cerebral Aneurysm Using the Lattice Boltzmann Method with the Half-Way and Interpolated Bounce-Back Schemes

: Lattice Boltzmann simulations and a velocity measurement of ﬂows in a cerebral aneurysm reconstructed from MRA (magnetic resonance angiography) images of an actual aneurysm were carried out and the numerical results obtained using the bounce-back schemes were compared with the experimental data to discuss the effects of the numerical treatment of the no-slip boundary condition of the complex boundary shape of the aneurysm on the predictions. The conclusions obtained are as follows: (1) measured data of the velocity in the aneurysm model useful for validation of numerical methods were obtained, (2) the numerical stability of the quadratic interpolated bounce-back scheme (QBB) in the ﬂow simulation of the cerebral aneurysm is lower than those of the half-way bounce-back (HBB) and the linearly interpolated bounce-back (LBB) schemes, (3) the ﬂow structures predicted using HBB and LBB are comparable and agree well with the experimental data, and (4) the ﬂuctuations of the wall shear stress (WSS), i.e., the oscillatory shear index (OSI), can be well predicted even with the jaggy wall representation of HBB, whereas the magnitude of WSS predicted with HBB tends to be smaller than that with LBB.


Introduction
A cerebral aneurysm is a disease in which the blood vessel of a cerebral artery is deformed into a saccular geometry. Rupture of cerebral aneurysms induces subarachnoid hemorrhage, which is an extremely serious disease with a high probability of death or poor prognosis, and therefore, it is important to understand the mechanisms of growth and rupture of the cerebral aneurysms. Although it has been pointed out that hemodynamic parameters, e.g., the wall-shear stress (WSS) acting on the artery wall, relate with the growth and rupture mechanisms [1,2], our knowledge on the mechanisms is still insufficient.
Computational fluid dynamics (CFD) is expected as a powerful tool to analyze fluid flows inside cerebral aneurysms. Foutrakis et al. [3] used artificial aneurysm models, which were two-dimensional channels with branches and convex geometries, to numerically investigate the characteristics of flows in the aneurysms. Valencia et al. [4] carried out three-dimensional simulations of flows in artificial aneurysm models and discussed the effects of the angle of connection between the aneurysm and the artery on the flow. These numerical studies were carried out using the artificial aneurysm models. On the other hand, numerical simulations with aneurysm models reconstructed from actual patient data have been carried out since the early 2000's [5][6][7]. Due to the complex shapes of actual cerebral aneurysms, it is desired to simply model the boundary condition at the artery wall. Since the lattice Boltzmann method (LBM) [8] can easily deal with such complex boundary shapes by using the bounce-back schemes, we developed a numerical method based on LBM for simulating blood flows in cerebral aneurysms in our previous studies [9][10][11].
The LBM predicts fluid flows based on the motion of imaginary fluid particles and the no-slip condition at the solid wall is simulated by the bounce-back motion of the particles. The half-way bounce-back scheme (HBB) [12] is the simplest, but the solid boundary is represented by a stepwise shape along the computational grid lines even for curved boundaries. On the other hand, the linearly interpolated bounce-back scheme (LBB) and the quadratic interpolated bounce-back scheme (QBB) deal with the boundary shape more accurately [13]. The HBB and LBB have been used in simulations of flows in cerebral aneurysms [14,15]. However the validation of HBB and LBB were usually carried out only for simple geometries [13] and the effects of the boundary model on the predicted flows in cerebral aneurysms have not been discussed in detail. The QBB uses the quadratic interpolation of the velocity distribution function. Therefore QBB is expected to be more accurate than LBB, which is based on the linear interpolation. Sanjeevi et al. [16] have however pointed out that QBB is less stable compared with LBB.
In our previous study [10], we investigated the effects of the collision models (single-, two-and multiple-relaxation time collision models) and the viscous stress model on numerical predictions of flows in cerebral aneurysms. However, the effects of the numerical treatment of the wall boundary condition on numerical predictions for the complex geometries were not discussed. For validation of numerical simulations with the complex shape of an actual cerebral aneurysm, experimental data must be obtained with the same channel geometry. In this study, we made a cerebral aneurysm model reconstructed from MRA (magnetic resonance angiography) images of an actual aneurysm for velocity measurement and lattice Boltzmann simulations and compared numerical results obtained using the bounce-back schemes with the experimental data to discuss the effects of the numerical treatment of the no-slip boundary condition on the predictions. The cerebral aneurysm model for the experiment was produced by making use of 3D printing from a STL (stereolithography) mesh of the cerebral aneurysm, which was also used to generate the scalar field (the signed-distance function) representing the boundary shape in the simulation to assure the same geometries in the experiment and the simulation.

Lattice Boltzmann Method
The outline of the numerical method is described in the following. The details can be found in our previous paper [10]. The numerical method is based on LBM, in which fluid flows are simulated by the motion of imaginary fluid particles. The motion of imaginary fluid particles is governed by the following lattice Boltzmann equation: where f (x, t) is the velocity distribution function of particles with the velocity c at the position x and the time t, ∆t is the time step size, Ω is the collision term, and the subscript i represents the direction of the particle velocity. The D3Q19 discrete velocity model for isothermal incompressible flows is adopted, i.e., i = 0, 1, 2, . . . , 18. The multiple-relaxation time collision model (MRT) [17] is used for Ω: where M is the matrix transforming f into its moment vector m (= M f ), f is the vectorial representation of f i , and S is the relaxation parameter matrix. The m eq is the moment vector for the equilibrium distribution functions f eq , which are given by where w i is the weighting function, u is the fluid velocity, and ρ is the fluid density. The ρ and u are given by

Boundary Condition
The LBM often uses the bounce-back scheme, which assumes that particles collide with the solid wall and bounce back in the direction opposite to the original direction so as to impose the no-slip boundary condition. In HBB, particles travel from a cell center toward the wall during the time duration ∆t/2, and collide with the wall at t + ∆t/2. Then, the particles return to the cell center at t + ∆t. Therefore whereī is the opposite direction of the i direction. Figure 1 shows an example of the applications of HBB to a curved boundary. The curved wall is represented by the dashed line in HBB. The jaggy representation of the curved boundary possibly deteriorates the accuracy of predictions in the vicinity of the wall, especially when the spatial resolution is low. The interpolated bounce-back scheme (IBB), which interpolates f i by taking into account the ratio, q i , of the distance between the wall and the cell center to the grid spacing ∆x, was proposed to improve the accuracy of the boundary representation [13] as shown in Figure 2. The linearly interpolated bounce-back scheme (LBB) is IBB with the first-order interpolation, and f i directing towards the wall at x A and t ( Figure 2) is changed after the time duration ∆t as The IBB with the second-order interpolation is the following quadratic interpolated bounceback scheme (QBB): where the subscripts A, B and C denote the cells in Figure 2.

Cerebral Aneurysm Model
STL data of a cerebral aneurysm is generated from MRA images as shown in Figure 3, where the STL mesh is a collection of triangular elements. The STL data are imported into our in-house LBM code, AN2WER [9][10][11] , and are used to reconstruct the signed-distance function (level set function) [18], ψ, which is used to evaluate the distance fraction, q i , in IBB. The ψ is the negative distance from the wall in the solid region, whereas ψ is the positive distance in the fluid region, and therefore, the zero-level set, ψ = 0, represents the wall. The ψ is reconstructed as follows. For the Cth computational cell at x C , the distance from the wall given by the eth triangular element is given by where x e and n e are the position vector and the unit normal of the eth element, respectively. Contributions from all the elements near x C are collected to obtain ψ(x C ) using the following weighting average [19]: where the power, p, for the weight is set to seven. Figure 4 shows the level set function reconstructed from the STL data in Figure 3.
By making use of ψ, q i in Equations (7) and (8) are given by   Figure 5a shows the experimental setup, which consists of the test section, the reservoir tank (TRUSCO, F-3GR), the pump (IWAKI, MD-15R-N), the flow meter (NIPPON FROW CELL, FDT-H-5) and the laser light source (KENTEC, LDB2W-H (wavelength: 450 nm)). The liquid phase was supplied from the tank to the test section using the pump and returned to the tank after flowing the test section.

Experimental Setup
The test section (Figure 5b) was the cerebral aneurysm model ( Figure 3) made by using 3D printing (Yasojima Proceed Co. Ltd., Kobe, Japan ). The aneurysm model had the single inlet and outlet, and its dimensions in the width, depth and height were approximately 50, 40 and 40 mm, respectively. The aneurysm models in the experiment and the simulation had the same geometry since the same STL data were used. The material was the transparent resin (stratasys, VeroClear). The resolution of the 3D printing was 30 µm per layer. The model was produced from two parts for smoothing the inner surface.
UV-cured coating was applied after smoothing, so that the model boundary was rigid. The test section was surrounded by the rectangular box of transparent acrylic resin and the box was filled with the liquid phase. For the liquid phase, silicone oil (Shin-Etsu Silicone, KF-56A (Newtonian fluid)) was used for index matching [20]. The liquid phase supplied by the pump was flowed into the aneurysm model through a horizontal circular pipe of 24 mm inner diameter, D, and 450 mm length (hydraulic entrance section). The temperature of the liquid phase was 24.8 ± 1 • C and the flow rate was 2.38 L/min. The Reynolds number defined by was 150, which is the typical value in cerebral arteries [21], where U is the mean velocity (0.088 m/s), and µ is the liquid viscosity. The velocity field was measured using the spatial filter velocimetry (SFV) [22] as described in the next section. The laser beam generated by the light source passed through the cylindrical lens to form a laser sheet for SFV. The position of the laser sheet was 12 mm above the bottom of the aneurysm. The thickness of the laser sheet was about 1.0 mm. Tracer particles, SiC (the specific gravity 1.69, the average particle diameter 3 µm), were added in the oil for the velocity measurement. The high-speed video camera (Redlake, Motion Pro X-3) and the camera lens (Nikon, 60 mm f/2.8D) were used to take particle images.

Spatial Filter Velocimetry
The SFV was used to measure the velocity in the aneurysm. The outline of SFV for the measurement of the velocity component in the x direction is as follows. Grayscale images of tracer particles are obtained by using a high-speed video camera. The images are divided into smaller pieces, which are called the interrogation areas (IA). The mean particle velocity in IA can be obtained by evaluating the change in the intensity of grayscale value due to the motion of particles. In order to evaluate the change in the intensity, the following spatial filter, F , is applied to the grayscale values: F (x, t) = cos(kx + ωt) (13) where k is the wave number of the spatial filter, and ω is the angular frequency. The integrated intensity, I(t n ), in IA of the nth frame image is given by where the subscripts l and m are indices of the pixels in the x and y directions, respectively, and g lm is the grayscale value of the pixel at (x l , y m ). The I is normalized as The m I and σ I are given by where N is the total number of images. The peak frequency, f p , of I is obtained by applying the wavelet transformation, and the mean velocity is then evaluated as where the shift frequency, f s , is given by (19) and the filter pitch, L, is given by The f s is introduced to detect particles moving in the −x direction. The velocity component in the y direction is obtained in the same manner. See Hosokawa et al. [22] for more details.

Inflow Condition
A velocity measurement of the flow in the circular pipe (hydraulic entrance section) was carried out to obtain the velocity distribution at the inlet of the test section, which will also be used as the inlet condition in the numerical simulations of flows in the aneurysm. Table 1 shows the parameters for SFV. The measurement positions of each velocity data were corrected by taking into account the refraction at the pipe surface. Figure 6a shows a comparison between the measured data of the time-averaged streamwise velocity u and the Poiseuille distribution, i.e., u/u max = [1 − (r/R) 2 ], where u max is the maximum velocity at the pipe center, r is the radial coordinate, and R is the pipe radius. The data agree well with the laminar velocity profile, confirming that the flow was fully-developed in the entrance section and the inlet condition in the numerical simulations of flows in the aneurysm could be simply given by the Poiseuille distribution.  Laminar flows in the circular pipe were predicted using HBB, LBB and QBB. Table 2 shows the numerical conditions. The fluid was Newtonian and the Reynolds number was the same as in the experiment. The relaxation time, τ, was 0.517 and the liquid kinematic viscosity in the LBM simulation was given as ν = (τ − 1/2)/3 for ∆x = ∆t = 1. The number of cells in the streamwise direction was 400. The computational domain for the cross section was 48 × 48 cells, whereas D was 10% smaller than the domain size, i.e., D = 43.2 cells. Streamwise velocities predicted using HBB, LBB and QBB are shown in Figure 6b-d, respectively. In all the cases, the predictions agree well with the laminar velocity profile.
The spatial resolution of 43.2 for D is therefore sufficient to obtain a good result even with the jaggy boundary representation of HBB. However, in the cells closest to the pipe wall, the predicted velocity in HBB is slightly larger than the analytical value. The deviations of the predicted u in those cells from the analytical value (0.0087) were 0.027, 0.0052 and 0.00011 for HBB, LBB and QBB, respectively. The distance between the pipe axis and the wall on the line passing through the axis is 22.0 in HBB, which is slightly larger than the actual radius D/2 = 21.6 due to the stepwise representation of the wall. The fluid domain is therefore bit wider on this line, which resulted in the predicted velocity larger than the analytical value. On the other hand, the smaller errors in the velocities obtained with LBB and QBB clearly show that the interpolation calculation treats the position of the wall more accurately. Figure 7a shows the velocity field in the cerebral aneurysm model measured by using SFV, the parameters of which are shown in Table 3. The liquid flows from the left top region of the image into the aneurysm along the left-side wall and a vortical structure is formed inside it. The center of the vortical structure is close to that of the aneurysm. The magnitude of the velocity inside the aneurysm is much smaller than the mean velocity at the inlet of the test section, i.e., U = 0.088 m/s. Numerical simulations of flows in the cerebral aneurysm model were carried out. Table 4 shows the numerical conditions. Although the simulations were carried out using HBB, LBB and QBB, the simulation with QBB was unstable under the given conditions and a steady state solution could not be obtained. The numerical stability of QBB is therefore lower than the others in the present case. A similar trend was also confirmed by Sanjeevi et al. [16], i.e., for a high Re flow in a simpler channel geometry with a strong pressure gradient, QBB tends to be less stable than LBB.  Figure 7b,c show the velocity fields predicted using HBB and LBB, respectively. The number of the vectors was reduced in the figures for visualization purpose. Being similar to the data, the liquid enters into the aneurysm model along the left side wall and a vortical structure is formed inside the aneurysm in the predictions, which confirms that both HBB and LBB give qualitatively reasonable results. The differences between the results of HBB and LBB seem small, and therefore, from the qualitative point of view, the effects of the numerical treatment of the no-slip boundary condition on the predictions of the flow structure in the aneurysm are not significant. Figure 8 shows streamlines obtained using LBB to understand the flow structure in the whole aneurysm model. The liquid from the circular inlet flows almost rectilinearly in the main artery section, whereas the streamlines are largely bended due to the interaction with the aneurysm wall, which corresponds to the left side wall of the aneurysm shown in Figure 7. Swirling motion is then formed inside the aneurysm and the streamlines direct toward the artery outlet.

Steady Flow in Cerebral Aneurysm
The measured and predicted velocities at the positions indicated by the blue squares on the lines, H1 (horizontal), H2, V1 (vertical) and V2, shown in Figure 7a are compared in the following to discuss the validity of the predictions in more detail. Figure 9 shows the measured velocities (the open circles) on H1, where v x and v y are the x and y components of the velocities normalized by the mean velocity at the inlet. On H1 in Figure 7a, the flow directs downward (+y direction) for x < 0.02 m due to the liquid inflow from the main artery, whereas for larger x the liquid in the vortical structure flows from the right to the left (−x direction). Therefore the values of v y for x < 0.02 m are positive and large and, for larger x, v x has negative values and its magnitude is much larger than that of v y . The v x and v y of HBB and LBB agree well with those of SFV though some deviations appear, especially for x < 0.025 m. No significant differences are found in the predictions of HBB and LBB, that is, the sum of the absolute differences between the velocity components of HBB and LBB on each line are less than several percent of the sum of the absolute value of the velocity components.  Figure 10 shows the velocities on H2, which passes the center of the vortical structure (x ∼ 0.03 m) and therefore the velocity components change their signs near that position. The liquid flow along the left side wall exhibits large velocities and passes H2 almost diagonally. Therefore the magnitudes of v x and v y are large in the small x region. The predictions of HBB and LBB reproduce well the velocity profile though they are somewhat smaller than those of SFV.    Figures 11 and 12 show v x and v y on V1 and V2, respectively. Being similar to the velocity components on H1 and H2, the velocity profiles, in other words, the vortical structure in the cerebral aneurysm model, are well predicted using HBB and LBB, and the predictions of HBB and LBB are not so different.
As discussed above, the velocities in the aneurysm model of HBB and LBB agree with the data and the differences between the predictions of HBB and LBB are small, which confirms that both numerical treatments of the wall can give reasonable predictions for the present complex geometry and the effects of the wall boundary treatment on the predicted flow structure in the aneurysm are small. We should however examine also the effects of the wall boundary treatment at lower spatial resolutions since the predictions of HBB would be worse than those of LBB as the spatial resolution decreases due to the stepwise representation of the wall. We therefore carried out the simulation with lower spatial resolutions. The spatial resolutions additionally tested were N x × N y × N z = 128 × 172 × 178 and 81 × 109 × 114. The other conditions were the same as in the previous simulation. The numbers of cells for D were 31 and 20 for the former and the latter, respectively, and τ were 0.5125 and 0.5079. It should be noted that a spatial resolution lower than N x = 81 was also tested, but a numerical instability took place and a steady state solution could not be obtained. Figure 13 shows the comparison of the velocity components, v x and v y , between N x = 162, 128 and 81. The flow characteristics in the experiment are reproduced well with HBB and LBB even at the lowest spatial resolution although the decrease in the spatial resolution slightly deteriorates the prediction. Interestingly the prediction of HBB is not so different from that of LBB even at the lowest spatial resolution, and therefore, the jaggy representation of the aneurysm wall is still acceptable to reproduce the flow structure.

Wall Shear Stress and Oscillatory Shear Index
Some hemodynamic parameters, e.g., the time-averaged wall shear stress (TAWSS), the oscillatory shear index (OSI) and so on, have been proposed to describe the growth and rupture mechanisms of cerebral aneurysms so far, where TAWSS is the magnitude of WSS acting on an aneurysm wall averaged for the single cardiac cycle, T, and OSI is an indicator representing the degree of directional change of WSS during T. In our previous study [11], we reported that thin aneurysm wall regions of un-ruptured cerebral aneurysms are correlated with OSI at a high probability, i.e., a low OSI region in numerical predictions corresponds to a thin wall region observed in an endovascular surgery. The TAWSS also showed correlations with the thickness of the aneurysm wall, i.e., a high TAWSS region corresponds to a thin wall region, whereas OSI would be more effective than TAWSS.
We therefore carried out numerical simulations of flows in the cerebral aneurysm model with an unsteady inlet condition and evaluated TAWSS and OSI to discuss the effects of the wall boundary treatment on these hemodynamic parameters. The Casson model [23] was used to account for the shear-thinning characteristics [10,24]. The inlet condition was time-dependent to roughly mimic the change in the flow rate during T, i.e., the flow rate increased at the early stage of T until it reached the peak value (the contraction period), and then, it decreased and returned to the initial value (the dilation period) (see Kimura et al. [9] for the detail).
The instantaneous wall shear stress vector, τ W , is given by where x is the position of the center of a cell closest to the wall, I is the unit tensor and n is the unit normal to the wall given by n = ∇ψ |∇ψ| (22) The TAWSS is given by The OSI is a hemodynamic parameter that describes the fluctuation of the WSS vector during T and is given by [25] The spatial resolutions were N x = 162, N x = 128 and N x = 81, Re was 150, and the Bingham number, Bi, was 0.1, where Here τ 0 is the yield stress and µ P is the reference viscosity. The time average of U was set to 0.02. Figure 14a,b show TAWSS predicted using HBB and LBB, respectively, with N x = 162. The distributions of TAWSS of HBB and LBB are similar although TAWSS of HBB is somewhat smaller than that of LBB in some areas, e.g., the region indicated by the dashed oval in the figure. Table 5 shows TAWSS averaged for the aneurysm. The TAWSSs of HBB are about 11% smaller than those of LBB at all the spatial resolutions. The TAWSS depending on the wall boundary model implies that the velocity vectors near the aneurysm wall are different in HBB and LBB though the velocity profiles in the entire aneurysm are almost the same as shown in the previous section. Table 6 shows the mean values of the magnitude of the velocity component tangent to the wall, v t = |(I − nn) · u(x)|, in the computational cells closest to the aneurysm wall. Being similar to TAWSS, the mean tangential velocity of HBB is smaller than that of LBB regardless of the spatial resolution. Therefore the smaller TAWSS in HBB is due to the smaller velocity magnitude in the near wall region.  Figure 15 shows the predicted OSI of N x = 162, where the maximum value of OSI for visualization is set to 0.01 since regions of OSI < 0.01 were confirmed to have correlations with thin wall regions [11]. The distributions of OSI predicted with HBB and LBB show only a little difference. The mean OSIs for the aneurysm are shown in Table 5; the values of HBB and LBB are very similar. The same trend was also seen for N x = 128 and N x = 81. Hence HBB can predict fluctuations of WSS as well as LBB in spite of the jaggy representation of the curved boundary.
Thus the effects of the wall boundary treatment on predicted flow structure are expected to be small and fluctuations of WSS can be well predicted even with the jaggy wall representation of HBB, whereas the magnitude of WSS predicted HBB can be smaller than that of LBB.

Conclusions
Numerical simulations of a flow in a cerebral aneurysm were carried out using the lattice Boltzmann method with bounce-back schemes, i.e., the half-way bounce-back (HBB), the linearly interpolated bounce-back (LBB) and the quadratic interpolated bounce-back (QBB) schemes, to investigate the effects of the numerical treatment of the complex boundary shape of the aneurysm on the predicted velocity field and hemodynamic parameters. The cerebral aneurysm for the simulation was reconstructed from STL (stereolithography) mesh data of an actual aneurysm. For validation of the numerical method, experimental data of the velocity field in a cerebral aneurysm model, which was made by 3D printing with the STL mesh, were obtained. The conclusions obtained for the present aneurysm case are as follows:

1.
Measured data of the velocity in the cerebral aneurysm model useful for validation of numerical methods were obtained. The present experimental method is to be of use for further validation with different types of cerebral aneurysm.

2.
The numerical stability of QBB in the flow simulation of the cerebral aneurysm is lower than HBB and LBB. 3.
The flow structures predicted using HBB and LBB are comparable and agree well with the experimental data. 4.
The fluctuations of the wall shear stress (WSS), i.e., the oscillatory shear index (OSI), can be well predicted even with the jaggy wall representation of HBB, whereas the magnitude of WSS predicted with HBB tends to be smaller than that with LBB.