NUMERICAL INVESTIGATION OF BLOOD FLOW FEATURES IN INTRACRANIAL SACCULAR ANEURYSMS

This study aims to provide insight about how the hemodynamic factors change with artery curvature for a developing aneurysm during a cardiac cycle. The aneurysm is investigated in terms of the vortical structure and the shear stress along the curved artery wall for three developing stages (initial, intermediate and terminal stages), for three instances of a cardiac cycle (diastole end, systole peak and diastole start) and for three different vascular geometries. The stream function vorticity formulation is used with Newtonian constitutive relation. During the systole peak instance for all aneurysm stages, the central vortex squeezes the streamlines towards the distal neck of the aneurysm leading to maximum wall shear stress in the vicinity of the distal wall of the aneurysm. The radius of curvature of the artery and inertial forces increased the wall shear stress along the aneurysm wall. The wall shear stress changes direction and concentrates in the vicinity of the distal neck for all artery geometries. Secondary vortices are observed in the terminal stage during diastole end and diastole start instances for the straight arteries and lead to shear stress fluctuations along the wall. The observations of this study are discussed together with the relevant clinical and numerical literature.


INTRODUCTION
An aneurysm occurs when part of an artery wall weakens, allowing it to widen abnormally or lump out. The blood flow continues inside of this abnormal region and the rupture of the weakened wall causes high rate of mortality [1]. To investigate the hemodynamics of flow inside aneurysms, several studies have been conducted in the context of numerical, experimental and clinical research [2][3][4][5][6][7]. In a glass model of intracranial saccular aneurysm, Steiger and Ruben [2] observed flow instabilities at the Reynolds number 300 during the deceleration of the flow. Meng et al [8] showed that the curved-vessel aneurysm model has fundamentally different hemodynamics compared to straight vessel-aneurysm. The flows in the curved and the straight vessel geometries correspond to 'inertia driven flow' and 'shear driven flow' respectively. Gonzales et al [9] investigated the blood flow inside intracranial aneurysms where they emphasize the rapid change in the flow of a cardiac cycle results in rapid changes in wall shear stress and pressure. This initiates the aneurysm growth and rupture at the cavity neck. Valencia and Solis [10] investigated numerically the terminal saccular aneurysm in basilar artery. They modeled the artery wall using the elastic solid and Mooney-Rivlin hyper elastic models and modeled the blood using Newtonian constitutive relation. It is reported that the effective wall shear stress and deformation are observed at systole.
Clinical experiments and observations are time-consuming and expensive, and can be quite limited due to ethical concerns. Numerical modeling of the blood flow is an ideal way to investigate the efficiency of recommended treatment (such as the introduction of intravenous stent) and to understand the flow mechanism within the aneurysm. The numerical models in the literature either used the flow domain obtained from a three dimensional magnetic resonance imaging of the artery [4][5], or simplified the artery geometry to a constant diameter tube as shown in Figure 1 [6][7][8].
In this study pulsatile flow in an intracranial aneurysm is numerically investigated under the combined effect of the aneurysm size, the curvature of the artery and the flow rates. The aim of this study is to evaluate the flow field and the wall shear stress distribution along the saccular cavity wall during a cardiac cycle and use these findings to have a better understanding of the mechanism leading to the rupture of the vessel.

MATHEMATICAL MODEL AND VALIDATION
Three stages of aneurysm with three different artery geometries (straight, δ=1/6 and δ=1/4) are used for the flow domain. The artery curvature is defined as δ=1/R, where R is the radius of curvature of the artery measured from the neck of the aneurysm, as shown in Figure 1. The blood is modeled as a Newtonian fluid [10][11]. The Reynolds number reads: where ρ is the blood density, μ is the blood viscosity, V is the maximum velocity in the cardiac cycle and D is the length of the aneurysm neck. For human blood flow at 37 o C based on an approximation for blood, the viscosity density ratio is set to μ/ρ=0.027 cm 2 /s [3]. The flow rates are 126ml/min and 202.5ml/min; and the corresponding Reynolds numbers are 350 and 500 respectively where they are for the conditions at systole peak stage during a cardiac cycle. These flow rates correspond to realistic flow parameters in human cerebral arteries. The Dean number (Ɗ) which measures the curvature effects is defined as: where δ is the curvature, and the Dean number is ranged from 70 -125 for the curved arteries considered in this study. These values are in the range valid for human cerebral arteries where the Dean number is between 10 and 200 [8].
The pulsatile flow along the artery is replaced with the average velocity at the aneurysm neck shown in Figure 2. This shows a cardiac cycle of 60 beats per minute as given in the experimental study of Milner et al [5]. The curve fitting parameters for the pulsatile flow are listed in Table 1, where the average velocity along the aneurysm neck is defined with the following 9 th degree polynomial: where T is the non-dimensional time during cardiac cycle.
The mass and momentum conservation equations for incompressible fluid under the absence of body forces are given as, = .  where, p is the pressure, I is the identity tensor, μ is the viscosity and D is the rate of deformation tensor. The length of the aneurysm neck (D) and maximum velocity during a cardiac cycle (V) are used to non-dimensionalize the governing equations. The unsteady vorticity-stream function formulation is adopted in generalized body-fitted orthogonal coordinates. The sample elliptic meshes generated for the flow domain are shown in Figure 3 for three developing stages of a saccular aneurysm. The governing equations are solved using a second order finite difference method. The discretized vorticity transport equation is integrated in time using Runge Kutta 4th order scheme. Successive over relaxation (SOR) with Chebychev acceleration is used to solve the stream function equation. The aneurysm wall is considered as a rigid wall and no-slip boundary condition is applied for velocity components. Thoms's formula is used to evaluate the vorticity value at the aneurysm wall [12]. The CFD code is written and compiled using GFortran 95 and GNU Fortran compiler. The numerical scheme is tested for three different grid densities, 41×41, 61×61 and 81×81 for all aneurysm geometries. Mesh convergence results are shown in Table 2 for intermediate stage during a cardiac cycle for Re=350, δ=1/6, Ɗ=71.44. Based on these results all results presented in this study are performed using 61×61 grid elements with time step 10 -5 . The convergence criterion for SOR with Chebychev acceleration to solve the stream-function equation is of the order of 10 -6 .

Intermediate Stage Wall Shear Stress
Zhang et al. [6] investigated experimentally and numerically the saccular aneurysm along large arteries for laminar non-pulsatile flow at Reynolds number 256.6. Shishir et al. [7] used the same geometry as a test case to validate their numerical results in two dimensional flow domain. In Table 3, the average velocity value inside of the aneurysm as well as the average wall shear stress is compared with the results in the literature. It should be noted that Zhang et al. [6] and Shishir et al. [7] modeled the flow inside a saccular aneurysm using a simplified artery geometry and in the present study the pulsatile flow in the artery is replaced with the average velocity along the aneurysm neck. Although the flow inside the aneurysm is investigated with simplified boundary conditions in two dimensions, it can be observed that the average velocity in the aneurysm cavity and the average shear stress value along the aneurysm wall agree well with existing data in the literature. Bouillot et al. [13][14] investigated the saccular aneurysm experimentally and numerically. The saccular aneurysm geometry is idealized as a sphere and a straight tube. The details of the aneurysm geometry are described in Bouillot et al. [14]. Velocity profiles and streamlines are shown in the symmetry plane of the model at the systole peak and the late diastole for particle imaging velocimetry (PIV) and computational fluid dynamics (CFD) results in Figure 4(a). The equivalent Reynolds number for the aneurysm cavity is 250, where the velocity along the neck is around 0.2 m/s. The streamline patterns during the corresponding cardiac cycle instances evaluated using the present model and the simplified boundary conditions showed agreement with the results in the literature, as shown in Figure 4(b).

RESULTS AND DISCUSSION
The shape and the diameter of the artery have a combined effect on the flow structure within an intracranial aneurysm [6][7]. In this study the diameter of the artery is kept constant and the effects of curvature of the artery are investigated together with the shape and size of the developing aneurysm during one cardiac cycle. The solutions presented in this study correspond to the second cardiac cycle. The details of the idealized geometry of the developing aneurysm are shown in Figure 5, where the center of the aneurysm wall and the corresponding angles are given.  A common observation to all contour lines is that the primary vortex is closer to the distal wall (as shown in Figure 1) during the diastole end (1), it then moves closer to the same corner during the systole peak and finally moves towards the center of the aneurysm during the diastole start. This behavior has also been reported by Bouillot et al. [13] where only straight vessels are investigated in the terminal stage. In all initial stages, the primary vortex approaches the aneurysm neck as the cardiac cycle moves through stages 1 to 3. On the other hand, in the intermediate stage, the primary vortex approaches the aneurysm neck in the systole peak. Secondary vortex formation is observed only in the terminal stage of Case 1 (straight artery, Re=350), for diastole end (close to distal neck) and diastole start (close to the center of the aneurysm cavity), where the secondary vortices for this geometry are mostly observed in the opposite side of the cavity under higher inertial effects [13]. This hemodynamic behavior can be considered as a result of the 'shear driven flow' feature [2].  (3)). On each horizontal axis, the angles are measured from the center of the corresponding geometry and they are provided to locate the outermost line of the symmetry plane of the aneurysm wall (shown in Figure 5). The angle ranges are (2π/3:π/3) for initial stage, (π:0) for intermediate stage and (5π/4,-π/4) for terminal stage. On each vertical axis, the shear stress values are provided. The shear stresses are plotted in the clockwise direction along the aneurysm wall. It is observed that for all four cases, the wall shear stress during the systole peak (2) is higher than the diastole end (1) and diastole start (3)  During the terminal stage, the maximum stress changed direction and the systole peak occurred in the vicinity of the distal corner where a stress concentration is observed. These observations are supported by the previous clinical and numerical studies in the literature [9 and 15]. Although the wall shear stress distributions along the aneurysm wall for the curved artery cases (Cases 2-3-4) have a similar trend, the wall shear stress distributions for the straight artery (Case 1) are different. The rich vortical structure during the terminal stage (shown in Table 4), leads to stress fluctuations along the wall in diastole end and diastole start instances (1 and 3 respectively). In Case 1, the change in the shear stress value along the aneurysm wall during the systole peak instance (2) is more intense compared to other cases.  Tables 12 and 13, summarize the peak shear stress values, and demonstrate the effects of the curvature and the flow rate, respectively. In this study, the artery geometry is defined using the radius of the curvature. In the initial and intermediate stages for all cardiac instances the maximum wall shear stress value increases with the radius of curvature, where it is maximum for the straight artery (Case 1). However for the terminal stage, the maximum wall shear stress values for the curved arteries (Case 2 and Case 3) are higher than that of the straight artery (Case 1). As it can be observed in Table 13, for higher flow rates and the same artery geometry (Case 3 and Case 4), the maximum wall shear stress values are observed to be higher for all stages and cardiac cycle instances.
The recirculation occurs inside the intracranial aneurysms mainly due to a sudden increase of the diameter of the artery [15]. The circulation is also enhanced with high radius of curvature at the artery [16][17]. In Case 1 where the artery is straight, the secondary vortex is observed in diastole end and diastole start.
The wall shear stress cannot be evaluated experimentally or measured using current in vivo techniques and can only be determined from the numerical results [17]. The wall shear stress is evaluated from the gradient of the velocity, where high values indicate a sharper change in velocity. The wall shear stress has an important effect on the generation, progression and finally rupture of the aneurysm. The location of the maximum shear stress also has a significant role in the rupture mechanism. In all cases (Cases 1 -4) the peak values of the wall shear stresses are observed to be closer to the distal neck and the maximum value is observed in systole peak in the vicinity of the distal neck.
The vortical structure within the cavity of the aneurysm is laminar [9]. Although the flow is laminar, the wall shear stresses as shown in this study change rapidly within a cardiac cycle, which may cause vibration along the aneurysm wall and bruits which are significant observations during turbulent hemodynamic flow [3]. These vibrations may be one of the main reasons of the origin, progression and rupture of the intracranial aneurysm [9,15].
The mechanism of the origin, progression and rupture of the intracranial aneurysm can be further revealed with an analysis that considers the hemodynamics associated with the vascular tissue mechanics, which is not considered in this study. In this context, a fluid solid interactive analysis can be carried out by considering a nonrigid wall. Also, flow fluctuations and abnormal stress growth affect the alignment of endothelial cells which has to be considered in remodeling. The mechanobiology and the histopathology of these lesions need to be investigated using new computational models like Fluid Solid Growth Models (FSGM), associated with the hemodynamics [18][19][20][21].
In summary, the artery curvature and the developing aneurysm stages play an important role in hemodynamic patterns. The wall shear stress increases with increasing radius of curvature and inertial forces.

CONCLUDING REMARKS
In this study, the hemodynamic patterns are evaluated and analysed for three artery curvatures and developing aneurysm geometries. The governing equations are solved with a finite difference scheme using body fitted curvilinear coordinates. The pulsatile flow in the artery is replaced with a simplified velocity boundary condition along the aneurysm neck. The key observations from the numerical simulation results are: during a cardiac cycle the developing aneurysm vortex is squeezed towards the aneurysm neck in the direction of distal neck and the wall shear stress reaches its maximum value during the systole peak and increases with both the radius of curvature and the flow rate. The results obtained in this study may provide insight about how the hemodynamic factors change with artery curvature for a developing aneurysm during a cardiac cycle. These biomechanical factors can be used to determine the underlying reasons of origin, progression and rupture of an aneurysm.