Volumetric lattice Boltzmann method for wall stresses of image-based pulsatile flows

Image-based computational fluid dynamics (CFD) has become a new capability for determining wall stresses of pulsatile flows. However, a computational platform that directly connects image information to pulsatile wall stresses is lacking. Prevailing methods rely on manual crafting of a hodgepodge of multidisciplinary software packages, which is usually laborious and error-prone. We present a new computational platform, to compute wall stresses in image-based pulsatile flows using the volumetric lattice Boltzmann method (VLBM). The novelty includes: (1) a unique image processing to extract flow domain and local wall normality, (2) a seamless connection between image extraction and VLBM, (3) an en-route calculation of strain-rate tensor, and (4) GPU acceleration (not included here). We first generalize the streaming operation in the VLBM and then conduct application studies to demonstrate its reliability and applicability. A benchmark study is for laminar and turbulent pulsatile flows in an image-based pipe (Reynolds number: 10 to 5000). The computed pulsatile velocity and shear stress are in good agreements with Womersley's analytical solutions for laminar pulsatile flows and concurrent laboratory measurements for turbulent pulsatile flows. An application study is to quantify the pulsatile hemodynamics in image-based human vertebral and carotid arteries including velocity vector, pressure, and wall-shear stress. The computed velocity vector fields are in reasonably well agreement with MRA (magnetic resonance angiography) measured ones. This computational platform is good for image-based CFD with medical applications and pore-scale porous media flows in various natural and engineering systems.

Formulation of volumetric lattice Boltzmann method. The VLBM 28 was developed specifically for treating arbitrarily oriented boundaries with or without boundary movement. In the VLBM, fluid particles are uniformly distributed in lattice cells. By introducing a volumetric parameter P( x, t) , defined as the occupation of solid volume in the cell, i.e.,P(� x, t) ≡ �V s (� x, t)/�V (� x, t) , we distinguish three types of lattice cells in the simulation domain: solid cell (pure solid occupation, P = 1), fluid cell (pure fluid occupation, P = 0), and boundary cell (partial solid and partial fluid, 0 < P < 1), as illustrated in Fig. 2. Formulation of the VLBM is self-regularized to deal with the complex flow geometries through P . It consists of three operations: (1) collision considering the momentum exchange between the moving boundary and the flow; (2) streaming accompanying a volumetric bounce-back procedure in boundary cells; and (3) migrating volumetrically moving the residual fluid particles into the flow domain when the boundary swipes over a boundary cell toward a solid cell. The VLBM strictly satisfies mass conservation and can handle any irregular boundary orientation and motion with respect to the mesh. In this work, we consider fixed flow geometries only, thus P(� x, t) = P(� x), which is a one-time calculation from STL image data. Thus, the collision reverts to the conventional operation in LBM and the migrating operation is omitted.
We use the D3Q19 lattice model to formulate the volumetric LBEs. The 19 discrete molecular velocities are where c(= δx δt with δx and δt the lattice width and time interval respectively) represents the lattice speed in lattice units. We select c = 1 , i.e., δx = δt = 1 , meaning that the particles stream one lattice unit per time step. The sound speed of the D3Q19 model in the lattice units is c s = c √ 3 ; and the weighting factors ω i are ω 0 = 1 3 , ω 1−6 = 1 18 , and ω 7−18 = 1 36 .   www.nature.com/scientificreports/ The VLBM deals with the time evolution of the particle population in each lattice cell, and the volumetric LBE with BGK (Bhatnagar, Gross, and Krook) approximation 52 for the collision operation reads where n i ( x, t) , n eq i ( x, t) , τ , and F i ( x, t) are the particle population, the equilibrium particle population, the relaxation time due to fluid-particle collisions, and the forcing term, respectively, in cell x at time t. The forcing term is formulated 44,53 as with a the acceleration due to an external force. The equilibria for incompressible flows are where N(� x, t) ≡ 18 i=0 n i (� x, t) ≡ 18 i=0 n eq i (� x, t) and N(� x, t)� u(� x, t) ≡ 18 i=0 � e i n i (� x, t) ≡ 18 i=0 � e i n eq i (� x, t). The relation between the particle density distribution function f i ( x, t) in a node-based LBM and the particle population The implementation of Eq. (1) includes two operations. The first involves the collision operator, � i ( x, t) , resulting in n ′ i (� x, t) that is referred to as 'post-collision' particle population.
The second operation is for streaming, which reflects the particles moving from the current cell to the neighboring cells. Since a boundary cell is partially occupied by fluid, only a partial volume fraction of fluid particles can stream to a neighboring cell in the direction of the molecular velocity. In the original paper 28 , the streaming equation is formulated as follows.
where i * is the particle velocity direction opposite to the i-th direction, i.e., � e i * = −� e i . In Eq. (5), the first term is the streaming term. All the particles in the upwind neighboring cell, n ′ i (� x − � e i δt, t), stream to the current cell but the current cell can only accept . This part assumes that the upwind cell is a fluid cell. The second term is the bounce-back term. [1 − P(� x + � e i * δt)] · n ′ i * (� x, t) particles stream to the upwind neighboring cell from the current cell and P(� x + � e i * δt) · n ′ i * (� x, t) particles bounce back to the current cell. This part assumes that the current cell is a fluid cell. However, either the current or the neighboring cell can be a boundary cell, with which Eq. (5) becomes problematic because the capability to stream and/or receive particles is determined by the relative fluid volume between the two cells. Also, Eq. (5) does not apply for the scenario that both current and upwind cells are boundary cells with an equal fluid fraction. To address these deficiencies, we consider two scenarios in the streaming operation as illustrated in Fig. 3. Assume cell B, located at x , is receiving particles during the streaming, and cell A is cell B's upwind neighboring cell located at � x + � e i * δt . x + � e i * δt are both boundary cells. A parameter G is introduced to distinguish these two cases. (a) P B > P A , G = 1. The receiving particles are streamed from the upwind cell. Because the fluid fraction of B is smaller than that of cell A, only 1−P B 1−P A n ′ Ai particles of cell A can stream into cell B and no bounced-back particles. (b) P B < P A , G=0. The receiving particles consists of two parts: streamed from the upwind cell A (n′ Ai ) and the bounced-back particles from the current cell P A −P B 1−P B n ′ Bi * in the opposite direction.  Fig. 3a, the receiving particles are purely streamed from cell A. Because the fluid fraction of cell B is smaller than that of cell A, only 1−P B 1−P A n ′ Ai particles among the total particles in cell A in the i-th direction can stream into cell B and the remaining particles P B −P A 1−P A n ′ Ai will bounce back to cell A in the direction of � e i * . Figure 3b illustrates another scenario, P B < P A , in which the receiving particles are from two sources: n ′ Ai (streamed from cell A) and P A −P B 1−P A n ′ Bi * (bounced back from n ′ Bi * ). We now introduce parameter G as and generalize Eq. (5) as In the case of P B = P A , Eq. (7) becomes n i (� x, t + δt) = n i (� x + � e i * δt, t) , which Eq. (5) cannot achieve. The resulting density, velocity, and pressure of the flow system are obtained as follows: and where p 0 and ρ 0 (= 1) are reference pressure and density in lattice units, respectively. The LBM can handle arbitrary oriented boundaries with or without motion, and it satisfies mass conservation strictly. The formulation taking into consideration of bounce-back boundary conditions is self-regulated by the volumetric parameter P( x) that distinguishes fluid, boundary, and solid cells. The implementation of the volumetric LBEs only occurs in fluid and boundary cells. Thus, the key to the VLBM is the determination of P( x) for a flow domain.
Geometrical processing for volumetric parameter P x and local wall normality N w ( x). Geometrical processing is an important step in ICFD needed to connect image segmentation to CFD. Usually, it requires extra effort for geometry rescaling and mesh generation. When using the VLBM as the CFD solver, a seamless connection between the image information and CFD can be achieved by quantifying the volumetric parameter P( x) for the entire simulation domain. The normal direction of each boundary cell, N w ( x) , which is required for calculating WSS and WNS, is calculated concurrently. The geometrical process from STL data to the volumetric parameter P( x) and wall normality N w ( x) is seen in Fig. 1. The green part is the output of the STL image processing (blue part) and the input of VLBM (yellow part), enabling a seamless connection between image information and CFD. The geometrical process consists of three main components described in the following subsections.
Level set method for signed-distance field f ( x) from STL data. The level set method is a conceptual framework for using level sets as a tool for numerical analysis of surfaces and shapes. One can perform numerical computations involving curves and surfaces on a fixed Cartesian grid without having to parameterize these objects 54 . A signed-distance function φ( x ) determines the distance of a given point x from an interface, with the sign determined by whether x is inside or outside the interface. When a flow domain is closed and bounded by an interface Γ with inward normal N w , as shown in Fig. 4, the signed distance field φ( x) represents the geometric space using a level set of with The STL data describes a raw, unstructured triangulated surface by the unit normal and vertices (ordered by the right-hand rule) of the triangles using a 3-D Cartesian coordinate system. This file format is supported by many other software packages, e.g., MATLAB. We use existing open-source packages of READ_stl.m, VOX-ELISE.m 55 , and ac-reinit.m 56 to read the STL file, generate a 3-D uniform mesh with the specified resolution, and calculate the signed distance field φ( x) in MATLAB.

Volumetric parameter P x and wall-normal direction
For an illustration, we use a 2-D lattice, in Fig. 5, to describe how to determine the P value for a cell from the For a boundary cell, we use the following steps to calculate the P value.
• Uniformly divide the boundary cell into q 2 sub-cells. Each cell is represented by the square center x j s with a solid volume V j s , j = 1, . . . , q 2 .
• Calculate the φ value at its x j s using the known φ values at 4 vertices of the cell through linear interpolation for each sub-cell.
• Obtain the total solid volume through a summation of the solid volume of all the q 2 sub-cells: • Obtain the solid volume fraction, i.e., the P value, as P = V s /q 2 .
Note that more sub-cells, equivalent to larger q, result in a more accurate P value, but more computation will be needed to calculate the volume fraction. Thus, q should be appropriately chosen to balance between accuracy and computational cost.
The inward normal direction at the volumetric center of a boundary cell can be easily determined. The normal direction to the boundary wall is needed to calculate the WSS and WNS, see Eqs. (17) and (18) . Then based on the volumetric center    ∂x α , where α(= 1, 2, 3) and β(= 1, 2, 3) represent the indices of the axial directions in a 3-D Cartesian coordinate system. We use these notations in the LBM, it can be calculated from the non-equilibrium distribution function as 57 It has been demonstrated that Eq. (13) is more robust than the conventional finite difference method (FDM) to calculate the strain rate from the velocity field 44 .
In the VLBM, Eq. (13) becomes Observing that shear stress tensor σ αβ = 2µS αβ , pressure p = c 2 s N 1−P , and the viscosity µ = c 2 s τ − 1 2 , we can formulate the total stress tensor as where δ αβ is the Kronecker unit tensor.
The Cauchy formula gives the overall stress on the wall with a normal vector N w : in which Einstein summation convention with index notation has been used. Then its projection onto the normal direction yields the WNS Here γ (= 1, 2, 3) is also an index for the axial direction in the 3-D Cartesian coordinate system. The WSS is then computed as the difference between the overall stress and its projection onto the normal: With the determination of the local wall normality N w in Eq. (12), both WSS and WNS (when needed) can be obtained en-route during the VLBM implementation.

Benchmark study of oscillating flows in a 3-D pipe
We now apply the VLBM to simulate pulsatile laminar and turbulent flows in an image-based 3-D pipe, focusing on demonstrating the reliability and applicability of the computational platform by comparing the computational results with analytical and experimental results.
Computation set-up. A 3-D rigid, right-circular straight pipe is generated in Solidworks in STL format with a length L = 19.05 mm and a radius R = 9.525 mm. The pulsatile flow along the pipe in the z-direction is driven by either a pressure gradient, P(t) or a velocity profile, u in (x, y, t) , which will be specified below. A no-slip boundary condition is applied on the pipe wall, realized by the bounce-back boundary condition. It is noted that in VLBM, the bounce-back boundary condition is involved in the streaming operation, see Eqs. (6) and (7), regulated by the volumetric parameter P( x) . The initial conditions are constant pressure and zero velocity in the entire flow domain. A uniform mesh in the VLBM is used with the cell number, N D , across the pipe diameter to represent the spatial resolution.
The volumetric parameter P( x) of each lattice cell and the normal direction of each boundary cell N w ( x) are calculated based on the STL file by our in-house MATLAB code. As mentioned in "Volumetric parameter P( x) and wall-normal direction N w ( x) from signed distance field φ( x)", when the volumetric parameter P( x) of a boundary cell at location x is calculated, the boundary cell is divided into q 3 sub-cells and q should be appropriately chosen to balance between accuracy and computation cost. We perform an analysis on this theme. First, we evaluate the influence of the spatial resolution represented by N D through the relative error, δ V (%), of the P( x)-based fluid volume (summation of 1 − P(� x) over the entire computation domain) from the specified pipe volume ( 2πR 2 L ). Table 1 shows, when q = 8 is selected, the error reduces to below 1% when N D =127. We then stick to N D = 127 to study the effects of q on the computation time and accuracy of P( x) in Table 2. It is seen that when q increases, the computation time increases significantly but the accuracy remains almost the same. Since our general goal for accuracy is less than 1%, we use q = 8 to compute the P( x) field with N D > 127 in what follows. The wall-clock time was recorded on a computer with Inter®Xeon® CPU E5-1650 v3 @ 3.50 GHz, RAM 32 GB, 64-bit Operating System, ×64-based processor, and Matlab R2020a.  (2) where � a(t) = P(t) ρ � k along the z-direction. With the body force to drive the flow, we impose a periodic boundary condition along z-direction. The density and kinematic viscosity of the fluid is 1.0e3 kg/m 3 and 1.0e-6 m/s 2 , respectively.
We first perform a grid-function convergence check to determine the appropriate spatial resolution. Seven N D s of 31, 59, 87, 117, 153, 181, 209, and 241 are used. Figure 6a shows the effects of N D on the spatial mean downstream velocity u and mean shear stress σ on a cross-section. It is seen that each of the curves becomes flat as N D increases. The quantitative measurement of the convergence is shown in Fig. 6b, in which the convergence rate δ is calculated as the relative difference of the quantity of interest-velocity or shear stress-between two successive spatial resolutions. It indicates that when N D ≥150, δ < 0.5% for both velocity and shear stress. For the simulations of laminar and turbulent flows below, N D is in the range of 150-240. Correspondingly, the relaxation time τ is in the range of 0.55-0.95.
Order of accuracy. The order of accuracy of a finite-difference-related numerical scheme is usually evaluated by the discretization error. As pointed out in Ferziger's textbook 58 , the exact solution of the discretized equations on the grid h(= 2R N D ) , ϕ h , differs from the exact solution of the partial differential equation,Φ , by the discretization error δ h , i.e., δ h = Φ − ϕ h . For sufficiently fine grids, this discretization error is proportional to the leading term in the Taylor series: δ h ≈ αh q + H where H stands for higher-order terms and α depends on the derivatives at the given point but is independent of h . The discretization error can be estimated if solutions on systematically refined grids are compared through the expression Φ = ϕ h + αh b + H = ϕ 3h + α(3h) b + H with 3 times (for the cell-based cubic mesh) for the increase of reference spacing from h . The exponent b represents the order of the numerical scheme and can be estimated as follows We simulated the flow with three successive resolutions with N D = 31, 93, 279. Using the velocity magnitude at three locations r/R = 0, 0.5, and 0.86, we calculate the orders of accuracy for the three locations using Eq. (19) and obtained b = 1.2, 1.2, and 1.1, respectively, which indicates that the spatial accuracy of VLBM is 1.2 in interior and 1.1 at the boundary. It is well-known that the spatial accuracy of node-based LBM is theoretically no more than the second order in the interior. The accuracy reduces to the first order near the boundary with a bounce-back boundary condition. When the wall lies exactly at the nodes or middle of two neighboring nodes (half-way bounce-back). where u is the downstream velocity, J 0 (·) is the Bessel function of the first kind and order zero, Re(·) is the real part of a complex function, α(= R ω n ) is the Womersley number, and r(= x 2 + y 2 ) is the distance to the center of a cross-section. The strain-rate, S = ∂u ∂r , shear stress σ = 2µS , and WSS σ w = σ | wall can be calculated from Eq. (20). Keeping the same pipe configuration, we use P s = 0.3Pa , P o = 1.8Pa , and ω = 1.57s −1 to simulate the flow in this part. The Womersley number α and the base flow Reynolds number Re 0 = 2Ru 0 n , with u 0 the spatial mean velocity and ν(= µ ρ ) the kinematic viscosity, of the flow, are 6.89 and 9.4, respectively. Figure 7 compares the simulation results (symbols) with the analytical solutions (lines) at 8 representative time instants in a pulsation for the downstream velocity u (a) and shear stress σ (b) and good agreements between them are shown. More quantitative comparisons are presented in Table 3 with the relative difference between the numerical

Computational results vs. analytical solutions of Womersley flow.
where f is for u and σ w (WSS) on a cross-section of the pipe at the same representative time instants as embedded in Fig. 7b. Besides, we compared the strain rate (not shown) calculated from en-route VLBM implementation and post-processing of velocity field using FDM. The VLBM is more accurate than the 3-point FDM but nearly equal accurate between VLBM and 5-point FDM.
Last but not the least, we demonstrate how Eq. (7) improves the modeling of the generalized streaming operation in the VLBM. Keeping everything else the same, we use Eqs. (5) and (7) to solve the velocity field. The accuracy of the spatial mean velocity comparing with the analytical solutions from Eq. (20) at the 8 time-points in one cycle is presented in Table 4. The average relative error is reduced to 0.653% (Eq. (7)) from 0.721% (Eq. (5)).
where f is for u and σ w on a cross-section of the pipe at the same representative time instants as embedded in Fig. 7b. The average relative errors at the 8 time points are 1.53% for velocity and 3.98% for WSS. The analytical solution of σ w at t 8 tends to zero thus the relative error is not available. www.nature.com/scientificreports/ Turbulent pulsatile flows. We now validate the reliability and accuracy of the numerical simulation for more realistic pulsatile flows, relatively high base flow Reynolds number Re 0 , by comparing velocity and WSS waveforms with laboratory measurements 60 . Time-resolved particle image velocimetry (PIV) techniques were used to acquire two components of velocity vector data on the plane of laser illumination that included the longitudinal axis of an acrylic pipe. The optical system consisted of a single frequency continuous-wave AR + laser, model Spectra-Physics Millennia Vs of 5.5 W with continuous output power at 532 nm. A fast frame-rate CMOS camera, Vision Research v710 Phantom, was used to acquire flow images with a 1200 × 800 pixel resolution. In each experiment, 20,000 images were acquired with a rate of 1 kHz. Hollow glass particles of 9-13 μm diameter (type borosilicate glass spheres LaVision 110P8) were used as tracers to visualize the flow. The PIV images were processed using an in-house modified code based on the open-source PIVLab software for MATLAB platforms. This code is a multi-pass PIV FFT (fast Fourier transform)-correlation based algorithm. The initial interrogation size was 128 × 128 pixels. It is reduced to a final window size of 16 × 16 pixels during three iterations in order to improve the signal-to-noise ratio. Five different experiments were carried out with Re 0 in the range 535 to 4825 and α from 11.91 to 23.82. We simulated all the five flows along the downstream region from location 183 to 184.55 in the pipe via direct numerical simulation. The diameter of the pipe (19.05 mm) and the fluid properties, i.e., density ( 998.008kg/m 3 ) and kinematic viscosity ( 1.004 × 10 −6 m/s 2 ), are identical to the experiment. In the simulating region, velocity has been measured at 61 downstream locations. A velocity profile at each location contains velocity information at 39 spatial locations along the pipe diameter and 400 temporal samples during a pulsatile period, i.e., 0.25 Hz. We use the velocity profiles at downstream locations at 183 and 184.55 as the inflow and outflow boundary conditions, respectively. They are introduced in the VLBM via a non-equilibrium extrapolation scheme 61 . We use N d = 243 and τ = 0.506 for direct numerical simulations of turbulence. The spatial-temporal velocity information at each location is introduced in the D3Q19 lattice cell model on the corresponding cross-section via piecewise cubic Hermite interpolation. The comparisons of (a) velocity waves at the center (r/R = 0), middle location (r/R = 0.5), and near the boundary (r/R = 0.98) and (b) WSS wave between numerical simulation (lines) and experimental measurement (symbols) for the case of Re 0 = 4735 and α = 11.91 are shown in Fig. 8. The WSS representing the experimental results in Fig. 8b was calculated from the experimentally measured velocity fields via a 5-point forward FDM. Figure 8a shows that the simulated velocity waves in the inner pipe at r/R = 0 and 0.5 , agree with the experimental measurements very well. Whereas at the inner wall of the pipe ( r/R = 0.98 ), both velocity (Fig. 8a) and WSS (Fig. 8b) waves from simulation still agree with the experimental measurements, but the simulation results are consistently less noisy than those from experiments.

Application study
Beyond the benchmark study, we now apply the presented computational method to a medical research project to explore the capability of non-invasive quantification of 4-D hemodynamics in image-based human vessels. Our focus is the pulsatile velocity, pressure, and WSS in human vertebral and carotid arteries based on MRA. Cardiovascular disease is known to be strongly influenced by hemodynamic WSS. Although the advanced diagnostic imaging modalities such as digital subtraction angiography, CTA, ultrasound, and MRA can provide very good information of the vessel morphology and the total blood flow, they cannot measure the WSS, which is believed to play a pertinent role in the development of human vascular diseases as demonstrated in many studies [62][63][64][65] . Without the quantitative physiological parameters of the flow and blood-vessel interaction, it remains challenging to fully assess the true risk of vascular lesions to optimize treatment decisions. In this project, we aim to integrate patient-specific ICFD and WSS analysis into diagnostic MRA of vascular diseases focusing on carotid arteries. To demonstrate the reliability and applicability of the VLBM-based ICFD for realistic flow systems, we present one study case for the accuracy of the computed hemodynamics including velocity and WSS in human vertebral and carotid arteries following the flowchart in Fig. 1.
The image data were acquired from the scanning of 3-D time-of-flight (TOF) and peripheral gated (PG) phase-contrast (PC) MRA on a 3T clinical MRI scanner (Siemens Magnetom Prisma) using a 64-channel head and neck coil. As shown in Fig. 9a, the TOF-MRA images with a high resolution of 0.33 2 × 0.5mm 3 are to anatomically extract the artery system. The TOF-MRA images are on DICOM data. In general, the image segmentation can be done by different software such as Materialise Mimics and 3D Slicer and exported in an STL format to follow the illustrated flowchart in Fig. 1. Since we have developed a method to do image segmentation from DICOM format and calculate the volumetric parameter P( x) and the normal direction vector N w ( x b ) of each boundary cell 29,49 , the STL format is not necessary. We segmented the left carotid (with a bifurcation) and vertebral (tube-like) arteries, shown in Fig. 9a from the 3-D TOF-MRA images. The PC-MRA images in Fig. 9b, provide the 4-D velocity profiles on the cross-sections from the bottom up with a spatial resolution, 0.44 2 × 5mm 3 resulting in 10 slices in the vertical direction, and a time resolution of 25ms resulting in 20 time points in one cardiac cycle. We use the in-house Matlab code to extract the 4-D velocity vector fields from the PC-MRA images for both vertebral and carotid arteries based on the method we developed 66 . Figure 10 shows the velocity magnitude contours on Slice 1 of the vertebral artery at 7 representative time points in one cardiac cycle. The 4th time point corresponds to the peak systole of the heart pumping. We simulated the blood flow in both carotid and vertebral arteries. The dimensions of the two arteries are found in Fig. 9a. The TOF-MRA resolution results in grid sizes of 106 2 × 416 for the carotid artery and 60 2 × 416 for the vertebral artery. The time resolution is 1 ms via linear interpolation of the PC-MRA data in Fig. 9b. The blood flow is driven by the 3-D velocity profiles at the inlet (Slice 1) and outlet(s) (Slice 10). We impose the non-equilibrium extrapolation scheme 61 to introduce the velocity boundary conditions at the inlet and outlet(s). The density and kinematic viscosity of blood are 1.06 × 10 3 kg/m 3 and 3.3 × 10 −6 m/s 2 , respectively. The relaxation time τ = 0.505.
We first show the comparisons of computed vs. PC-MRA measured velocity magnitude contours in Fig. 11 (a) vertebral (Slice 4) and (b) carotid (Slice 5) arteries, to demonstrate the accuracy of the computation. On each panel, the top and bottom rows correspond to the peak systole and the end of diastole in a heart beating. Overall, the computed velocity agrees with the MRA measured ones reasonably well. The MRA measured velocity is by nature noisy, appearing at both systole (large velocity) and end of diastole (small velocity.) Fig. 12 shows the 3-D (a) velocity contours with streamlines, (b) pressure contours, and (c) WSS distribution of the carotid artery at the peak systole. While the large velocity and WSS deviate to the internal carotid artery (left), the large pressure deviates to the external artery (right). This is due to the asymmetry between the internal and external carotid arteries. Such a reveal of the 4-D in vivo hemodynamics in live human arteries is meaningful to the patient-specific diagnostics and therapeutics for the advance of precision medicine. Further postprocessing of the computed 4-D velocity, pressure, and WSS will produce more medical relevant results, which will be included in a separate paper in near future.

Summary and discussion
We have presented a unique computational method to quantify 4-D wall stresses of image-based pulsatile flows using the VLBM. With the introduction of a volumetric parameter to distinguish solid, liquid, and boundary cells, the VLBM is well-suited for dealing with complicated flow domains. The bounce-back boundary condition has been formulated in the streaming operation, thus avoiding complications of arbitrary boundaries. In this work, we have generalized the equation for the streaming operation from its original formation. The original streaming equation, Eq. (5), only accounts for streaming from a boundary cell to a fluid cell, which becomes problematic when the streaming is from a boundary cell to another boundary cell. We generalized Eq. (5), by introducing a parameter G, see Eq. (7), to distinguish two streaming scenarios and demonstrated the accuracy improvement.
Starting from the provided image information, generally, in STL data format, we first construct the signed distance field based on the concept of level-set methods and directly calculate the volumetric parameter, P( x) , and the local wall normality, N w ( x) . The P( x) and N w ( x) , together with inlet/outlet boundary conditions and initial conditions are directly fed to the VLBM, resulting in a seamless connection between the image processing and CFD. In the VLBM, the strain rate tensor can be calculated from the non-equilibrium particle population; therefore, the wall stresses can be obtained en-route during the VLBM implementation. These features, together with GPU parallel computing, dramatically reduce resource requirements for ICFD and ease the computational burden for pulsatile flows. We applied the computational method to solve the benchmark pulsatile flows in an image-based pipe. We first simulated a Womersley (laminar) flow with Womersley number α and the base flow Reynolds number Re 0 of 6.89 and 9.4, respectively. The computational results were in good agreement with the analytical solutions, with average relative errors around 1.52% and 3.98% for mean downstream velocity and mean WSS, respectively. We found that the VLBM en-route calculation of WSS is more accurate than the 3-point centered FDMs. We further simulated turbulent flows by varying Re 0 from 535 to 4375 via direct numerical simulation. The computational results of velocity and WSS are again in good agreement with the experimental measurements. We have also presented one application study on the non-invasive quantification of 4-D hemodynamics, including velocity vector and pressure fields and the WSS distribution in MRA-based human carotid and vertebral arteries. The computed velocity agrees reasonably well with the MAR-measured velocity. Three more applications are currently active, including (1) quantification of hemodynamic WSS on image-based human choroidal vasculature endothelium, (2) Drag force on particles in sand flow, and (3) effects of waste streams on ion exchange kinetics in the porous structure. Application (1) is the first-ever computational quantification of WSS in the human choriocapillaris. The findings have been recently submitted for publication 67 . Applications (3) and (4) are in the code adaption stage and results are expected in near future.
Although static boundaries have been the focus in this work, the method can be easily generalized for moving boundaries as the next step. The original VLBM is derived for moving boundaries, and we have herein indicated how this might be applied. The formulation of volumetric LBEs consists of three parts: (1) collision taking into account the momentum exchange between the moving boundary and the flow; (2) streaming accompanying a volumetric bounce-back procedure in boundary cells; and (3) boundary-induced volumetric fluid migration moving the residual fluid particles into the flow domain when the boundary passes over a boundary cell toward a solid cell, and mass conservation is guaranteed in VLBM (no need to solve a pressure Poisson equation as in NS formulations). If time-varying image information is available, the presented computational method can be generalized to solve pulsatile flow problems in image-based moving boundaries.

Data availability
The datasets generated in the current study are available from the corresponding authors upon appropriate request for nonprofit use.