Meshing strategy for bifurcation arteries in the context of blood flow simulation accuracy

The study presents a mesh dependency study for a carotid artery bifurcation geometry of a real-life specimen. The results of time-averaged velocity profiles at artery control surfaces and wall shear stresses are compared between a set of structured and unstructured meshes, with varying non–dimensional boundary layer first element thickness (y+) parameter. A set of four meshes in total is considered: a full–hexagonal structured mesh, an unstructured tetrahedral mesh with prism inflation layer, both created for y+=1 and y+=30. Apart from numerical results, overall mesh creation work time, overall analysisstability are compared with the mesh quality results: cell non–orthogonality, cell skew and aspect ratio. Numerical results are validated against results of real–life CT examination performed in Poznań Medical University.


Introduction
Application of numerical methods in the field of biomechanics of flows allows distinguishing the mechanical reasons for the occurrence of some pathological changes such as progressive stenosis. Physical quantities used to describe the phenomena occurring in fluid dynamics can also be considered in medical aspects to describe the causes of some specific group of vascular diseases, such as atherosclerosis or cerebral aneurysm [1,2].
When numerical modelling this type of flow, the key aspect in the detection of boundary layer detachment is the turbulence model and the quality of the mesh: especially its resolution and non-dimensional boundary layer first element thickness (y+) parameter.
A literature study conducted for this research [3][4][5] shows that computational grids created for the analysis of blood flow in an artery are, typically, unstructured triangular and tetrahedral meshes prepared for inviscid computations The results of the numerical analysis of blood flow ⇤ e-mail: natalia.gr.lewandowska@doctorate.put.poznan.pl conducted on such meshes are flawed. Low quality and no inflation layer grids do not provide the necessary resolution to model the detachment of the boundary layer in the flow. This contributes to errors in calculation of the values of shear stress on the walls, as well as the velocity distribution and mass flow in the individual branches of the modeled artery.
The research presented in this article shows how the simplifications used in the generation of computational grids a↵ect the accuracy of calculations. Four meshes with di↵erent element types, resolutions, y+ parameter values are given. The comparative analysis will include accuracy in the analysis of wall stresses, detachment of boundary layer and velocity values at outlets. Validation of numerical calculations with results from Doppler examination allows for assessing the accuracy of calculations.

Problem description
The object of study is the specimen for the carotid artery. It is a system of three interconnected arteries. The common carotid artery (CCA) supplies blood from the aorta to an area called the bulb. In the bulb, the CCA is branching into an internal carotid artery (ICA) and external carotid artery (ECA). The ICA supplies blood to brain tissue. The ECA provides the blood circulation of the facial muscles. The plaque most often accumulates in the bulb, on the side of the ICA [6]. The deposition itself is caused by disturbed blood flow in this area. As a result of the increasing crosssectional area of the artery, a positive pressure gradient occurs, which in some cases causes separation of the boundary layer and formation of vortices. The resulting vortex causes the suction of solid particles in the blood, such as atherosclerotic plaque, erythrocytes, etc. Particles begin to adhere to the walls of the artery causing local stenosis.
A metric used for assessing the risk of plaque deposition with CFD methods is Wall Shear Stress (WSS) calculated in the first element at the wall boundary condition. In order to obtain the WSS figure coherent with the Doppler ultrasound examination data, a relevant resolution of the computational grid is required.

Geometry description
The geometry was created based on the computed tomography (CT) imaging of the patient's head. Patient was suspected of progressive stenosis of the artery. The diameters of the CCA, ICA and ECA are 9.6 mm, 8.3 mm and 6.25 mm respectively. The maximum angle of bifurcation between ICA and ECA in the bulb is 45 degrees. Total length of the CA geometry is around 200 mm. Figure 1 shows the outline of the geometry created in ICEM, based on the formatted point data obtained from the CT scan. The CCA branch is the lower left branch, the ECA is right-upper and the ICA is the right-lower branch.

Boundary conditions, input data
Inlet boundary condition is implemented at CCA and is a velocity inlet with a time varying velocity magnitude. The velocity function in time domain is obtained from results of carotid ultrasound Doppler examination (USG). The velocity profile given in the USG is approximated with a group of spline functions. The values of coefficients of interpolating polynomials are developed based on a suitable group of equation boundary conditions ensuring continuity of functions at the contact points of these functions and achieving values in characteristic points.The function consists of three 3rd degree polynomial function. The vaules of function coefficients and timesteps are presented in table 1.
The static pressure function at the outlet arteries is described using the 2-parameter Windkessel model [7][8][9]. The blood is modelled as a non-Newtonian fluid. The variable viscosity of the fluid is implemented in the form of the Carreau model with the value coefficients obtained by [10].

Spatial discretization
Four computational grids are prepared for the described study: • structured, hexagonal, multiblock, y+ = 1 • structured, hexagonal, multiblock, y+ = 30 • unstructured, tetrahedral with prism inflation layer, y+ = 1 • unstructured, tetrahedral with prism inflation layer, y+ = 30 Grids were created in ANSYS ICEM CFD meshing software. Structured, multiblock grids utilize a nontrivial Y-O topology. Figure 3 presents a detail of block structure in the area of the artery bulb. Following mesh quality metrics are used to assess the mesh quality: orthogonal quality, skew and maximum aspect ratio. Orthogonal quality is defined as following: For each face, the cosines of the angle between face normal vector (Ai) and the vector from the cell centroid to the centroid of the adjacent cell (ci), and between Ai and the vector from the cell centroid to the centroid of the face (fi) are calculated. The smallest calculated cosine value is the orthogonality of the cell. The worst cells will have an orthogonal quality close to 0 while the best cells will have an orthogonal quality close to 1 [11]. The skewness is defined di↵erently for volume and surface elements. In all cases it is normalized so that 1 is ideal and 0 is the worst possible. For a hexahedral element, skewness is defined as the normalized worst angle between each of the 6 face normal vectors and the vector defined by the centroid of the hexahedron and the center of the face. For tri elements, skewness is defined as the ratio between the area of the element and the area of an equilateral triangle having the same circumcircle.
Aspect Ratio (AR) is computed as the ratio of the maximum value to the minimum value of any of the following distances: the normal distances between the cell centroid and face centroids computed as a dot product of the distance vector and the face normal, and the distances between the cell centroid and nodes. Table 2 presents mean and worst quality parameters of interest for set of meshes. Mean values of quality parameters are similar for all four of the grids. Hexagonal grids have better overall orthogonal quality, however the tetrahedral grids with prism boundary layer (tri, pl) have less skewed elements and better aspect ratio. Both tri and hex meshes for y+1 wall spacing were created in an approximately the same work time of 14 working hours. Meshes with y+30 wall spacing were derived from the y+1 grids by changing the element spacing parameters at wall boundary conditions.

Solvers and setup
The calculations were carried out using the Ansys Fluent software. A pressure-based solver, with Coupled pressurevelocity coupling scheme, was used. A k-omega SST turbulence model with low-Re corrections was implemented [12]. Total flow time is 0.6 s, which comes from patients heart rate of 100BPM during the CT examination. The time step is 1e-04s, and was suggested by [13] as a largest possible time step for obtaining relevant calculation results. A metric used to compare the velocity profiles is a modification of standard deviation of discrete random variable, where the mean value of the measured velocity at given outlet is replaced velocity measured at given time (eq. 1). Smaller value denotes smaller deviation from the measured velocity profile.
Values of the sigma RMS for calculated meshes are presented in Table 3. It is found that mesh no. 3, the structured hexagonal grid with y+=1 has the smallest sigma RMS value. Because the experimental values of the Wall Shear Stress in the given artery are unavailable, it is decided that results obtained from CFD calculations on the hexagonal y+=1 grid are considered as baseline.  the WSS distribution in the bulb area for the particular stages of the cardiac cycle: time of achievement of PSV (t = 0.04 s) and transition time between the systolic and diastolic phase (t = 0.14 s).During the peak systolic phase, WSS values are of highest magnitude, being dependent from the velocity magnitude. Distribution in time 0.04s is similar (Fig. 4) throughout all cases: the most significant concentration of high value WSS occurs at the outflow from the bulb. It is related to the decreasing cross-sectional area: it's causes local acceleration of the flow. Its impact on values of shear stresses.

WSS analysis
The most noticeable di↵erences between the cases occurs in time t = 0.14s (figure 5). In this stage, the velocity value achieves EDV, it is also the moment when the pulse wave is forming. Location of maximum WSS is similar, but the first two cases understate the value of WSS on the right side of the ICA. The region near the ICA wall, just after the outlet from the ICA is one of the most common locations on plaque deposition [6].
To more accurately assess how the type of mesh applied influences WSS values, a point-to-point comparison analysis was carried out. Figures 6-7 show the distribution of stresses against reference WSS obtained from the case no.3, the hexagonal y+ = 1 grid. Details of the validation process are presented in the "Validation" section. The values for the plots are obtained by applying formula 2 to all nodes on the compared mesh: This approach allows for assessing the WSS o↵set from the validated case. The most noticeable di↵erences in stresses occur for the tri grids: the di↵erence in the lower area of the bulb can be observed and on the ECA artery where the WSS value exceeds the up to seven times the reference value (colormap was limited to +/-300Pa for plot clarity). In the figure 7, underestimated values of WSS in the bulb and artery ICA are present. The best correlation with the hex_01 grid is presented by the hex_30 case: the underestimation trend of the WS values in the ICA is comparable with the tri meshes. However, smaller areas of increased WSS can be observed: their concentration is low and mainly occurs in the lower part of the bulb and at the inlet from the bulb to the ECA.

Velocity distribution
The control surface for velocity distribution plots is the axial section of the artery. Accurate resolution of the separation of the boundary layer is desired in the case of analysis of blood flow in the aspect of plaque deposition [14]. Figure 8 shows the separation of the boundary layer in the ECA artery. This flow separation is not in the interest of this study. Usually the plaque deposition in ECA is not considered, because the facial muscles are supplied with blood by many other arteries and blocking the flow in the ECA artery does not have serious consequences for the patient. The velocity distribution during the PSV is very similar for all cases.  Figure 9 shows the velocity distributions in flow time = 0.14 s. The decrease in velocity and separation on the boundary layer in the bulb is most noticeable during the transition between the systolic and diastolic phase. However, in the case of tetrahedral meshes, separation in the widest part of the bulb was not detected. Failure detection of separation in the case of unstructured meshes may a↵ect the size of vortices forming on the side of the artery ICA. Additionally, in case of tri the velocity in ICA is approximately 0.8 m/s, while for hex grids is equal to 0.6 m/s. That may lead to conclusion that tri meshes can overstate the value of the velocity.

Summary
The most significant di↵erences to the reference WSS appear for tetra meshes. The WSS distribution di↵ers significantly in the lower part of the bulb, where the tri grid overstates these values more than thrice. Additionally, in the case of tetra meshes, there was no detection of an increase in WSS in ICA in the region of frequent deposition of plaque. The hex_30 grid has less significant disproportion with respect to the reference grid. Two locations were observed where the WSS were substantially increased: the bottom of the ICA side and at the outflow from ECA. The tri grids did not detect flow separation in the upper part of the bulb and inflate velocity in the ICA by about 0.2 m/s. The velocity distribution for hex grids is comparable.

Conclusions
Based on the conducted research, it can be noted that while choosing the type of computational grid for modelling blood flow in the arteries, the key flow parameters to be analyzed should be taken into account. In the case of flow analysis in the context of aneurysm formation, WSS are mainly analyzed. The use of unstructured meshes leads to an overestimation of the WSS value and an indication of the locations of potential aneurysm growth that have no physical sense. In the case, when most important is the blood flow distribution (in terms of the tendency of the artery geometry to create local stenosis), the y+ parameter is not so important, the type of the elements used is significant. In this case, much more accurate results can be obtained by using structural, hexagonal meshes. They are more e↵ective in detecting separation areas in the bulb. If both of these parameters are taken into consideration, it is recommended to perform a fully structural mesh with the value of y+ equal to 1. The implementation of good quality mesh, especially for such irregular geometries as the carotid artery, is time-consuming and the final grid has the most elements, but calculations made on such a mesh gives the most accurate results both in terms of WSS and flow patterns.