## Abstract

We investigate the flow dynamics and oscillatory behaviour of wall shear stress (WSS) vectors in intracranial aneurysms using high resolution numerical simulations. We analyse three representative patient-specific internal carotid arteries laden with aneurysms of different characteristics: (i) a wide-necked saccular aneurysm, (ii) a narrower-necked saccular aneurysm, and (iii) a case with two adjacent saccular aneurysms. Our simulations show that the pulsatile flow in aneurysms can be subject to a hydrodynamic instability during the decelerating systolic phase resulting in a high-frequency oscillation in the range of 20–50 Hz, even when the blood flow rate in the parent vessel is as low as 150 and 250 ml min^{−1} for cases (iii) and (i), respectively. The flow returns to its original laminar pulsatile state near the end of diastole. When the aneurysmal flow becomes unstable, both the magnitude and the directions of WSS vectors fluctuate at the aforementioned high frequencies. In particular, the WSS vectors around the flow impingement region exhibit significant spatio-temporal changes in direction as well as in magnitude.

## 1. Introduction

Image-based computational fluid dynamics (CFD) simulations have been used to correlate wall shear stress (WSS), pressure and other indicators with the formation and growth of an aneurysm (Tateshima *et al.* 2003; Cebral *et al.* 2005*a*; Torii *et al.* 2006; Chatziprodromou *et al.* 2007; Sforza *et al.* 2009). In particular, impingement zones and areas of high WSS have been correlated with aneurysm formation, growth in size, or rupture (Cebral *et al.* 2009). It has been argued that blebs, ‘blisters on the aneurysm wall’, form in high WSS regions where the flow impinges into aneurysms (Cebral *et al.* 2007). Some authors (Boussel *et al.* 2008; Jou *et al.* 2008) claim that aneurysms tend to grow at an area where the endothelial layer is subject to abnormally low WSS. Hence, the ability to identify these critical locations in aneurysms, i.e. abnormally high or low WSS regions and impingement locations, would have significant implications for treatment.

Flow impingement regions inside aneurysms are reported to have a particular pattern of WSS/pressure distribution. The apex of a flow impingement area, i.e. the local region of flow stagnation, is characterized by low WSS and higher static pressure, and is surrounded by a band of higher WSS and WSS gradient (Szymanski *et al.* 2008). An experimental study of an impingement area (Meng *et al.* 2007) showed that destructive remodelling involving loss of smooth muscle cells occurs at the region of high WSS and WSS gradient where the flow accelerates. Also, a more sensitive response of endothelial cells to turbulent oscillatory shear rather than to laminar and steady shear has been demonstrated experimentally (Moore *et al.* 1994; Malek *et al.* 1999; Meng *et al.* 2007). Therefore, it seems important to understand the dynamic oscillatory behaviour of WSS vectors in the impingement (stagnation) locations. To this end, we investigate systematically the oscillatory behaviour of WSS vectors inside an aneurysm and the flow instability due to the presence of an aneurysm at the posterior communicating artery (PCoA) origin.

Intraoperative Doppler recording of Steiger & Reulen (1986) showed flow fluctuations in the aneurysms of six patients out of 12. Spectra calculated from the parent vessel's recordings did not show peaks except from the base pulsatile flow while spectra of the data from the aneurysms displayed several small peaks between 5 and 20 Hz. Ferguson (1970) recorded bruits from the sacs of 10 out of 17 intracranial aneurysms during surgery. The murmurs from aneurysms reported by Simkins & Stehbens (1974) were associated with the fluctuating flow (Steiger *et al.* 1987). Sekhar & Wasserman (1984) also reported sounds from saccular aneurysms in eight of 11 patients and later they measured sound from experimental aneurysms in the common carotid bifurcation of dogs (Sekhar *et al.* 1991). Cebral *et al.* (2005*b*) performed simulations to show that ruptured aneurysms tend to have unstable (changing direction of inflow jet) aneurysmal flow unlike unruptured cases, but no further details of unstable flow are provided.^{1}

A lateral aneurysm seems to be geometrically analogous to a cavity in a channel—a prototype flow which has been studied extensively. For the simplest cavity geometries, a streamline departing from the surface on the upstream edge of a rectangular cavity reattaches at or near the downstream edge of the cavity as described by Batchelor (1956). Self-sustained oscillations including the cavity flow were reviewed in Rockwell & Naudascher (1979). Cavity flow and its stability were studied numerically by Ghaddar *et al.* (1986). In aneurysms in cerebral arteries, owing to the effects of complicated geometry, we expect to find distinct regions of inflow and outflow at the mouths of wide-necked aneurysms. Such pattern of inflow and outflow may create inflection (change of curvature) of the velocity profile, which is the classic Rayleigh condition for instability (Drazin & Reid 1981). This, in fact, was first hypothesized by Lighthill (1975). Hence, we expect the instability to develop, probably only transiently because of the pulsatile nature of the blood flow, manifesting itself through velocity fluctuations of a distinctly higher frequency than the cardiac pulse.

Stability of pulsatile flow has been studied experimentally and numerically (Nerem *et al.* 1972; Einav & Sokolov 1993; Straatman *et al.* 2002; Blackburn & Sherwin 2007). Experimental study of flow in a straight pipe by Einav & Sokolov (1993) showed that pulsatile flows become ‘turbulent’ during the deceleration phase of oscillation at Reynolds number around 2000. Straatman *et al.* (2002) performed a linear stability analysis for plane pulsatile Poiseuille flow showing that pulsation may trigger instabilities. Nerem *et al.* (1972) measured velocity in the aorta of dogs and reported transition to turbulence at the end of systole.

CFD simulations in the carotid stenosis (Fischer *et al.* 2007; Lee *et al.* 2008; Grinberg *et al.* 2009) were performed and intermittent turbulence during the decelerating systole was reported. However, the current work seems to be the first numerical study to report flow instability in aneurysms. Specifically we consider patient-specific geometries showing the spectra of fluctuating velocity inside aneurysms as well as time traces of WSS vectors on the aneurysm walls under pulsatile flow with volumetric flow rate (VFR) of 150–400 ml min^{−1} in the internal carotid artery (ICA). Through a series of simulations, we report that aneurysmal flows may become unstable during the deceleration phase. We also aim to understand the distribution of WSS and its directional changes in aneurysms during the instability period.

In the following, we first present details of the geometric models and the simulation parameters in §2 and subsequently results of simulations in §3. We discuss findings and limitations of this study in §4 followed by conclusions in §5. In appendix A, we present sensitivity studies, specifically: (i) a systematic resolution study and (ii) the effect of outflow boundary conditions on the velocity fluctuations.

## 2. Methods

In this study, we analyse three representative vessel geometries harbouring aneurysms with different characteristics selected from many patients we considered. For this retrospective review, institutional review board approval was obtained. Specifically, aneurysms in the supraclinoid ICA segment were considered; side (lateral) aneurysms are the most common and hence terminal aneurysms were not considered. Based on patients treated endovascularly, Murayana *et al.* (2003) reported that approximately 40.5 per cent are small aneurysms (less than 10 mm maximum size) with small neck (less than 4 mm), 28.7 per cent are small with wide neck (more than 4 mm), 17.4 per cent are large (more than 10 mm) and 6 per cent giant (more than 25 mm). In general, approximately 15–20% of patients will have multiple aneurysms. According to statistics, the examples chosen in this study are representative of three types of aneurysms occurring in the supraclinoid ICA: wide-necked, small-necked and multiple aneurysms (see figure 1).

The first (patient A) has a very wide-necked aneurysm, resembling a fusiform aneurysm, in the supraclinoid ICA. The second (patient B) has the PCoA harbouring a relatively narrow-necked aneurysm. The third (patient C) developed two wide-necked aneurysms in the supraclinoid ICA. One of them is located at the root of the PCoA and the other one is downstream within a short distance. The aneurysm of patient B and the upstream aneurysm of patient C have the PCoA branches coming out of the aneurysm wall. Hence, all chosen aneurysms are located between the ophthalmic triangle and the ICA bifurcation into the anterior cerebral artery (ACA) and the middle cerebral artery (MCA) and approximately at the level of the eyes, behind the eyes and towards the middle of the head.

Images from computed tomography angiography or three-dimensional digital subtraction angiography were anonymized to remove any patient-specific data, and then exported in a DICOM format. The image resolutions and more detailed information of these data are explained in table 1. The vessel geometries were manually extracted using in-house Matlab R2006b (The Mathworks, Natick, MA) tools, which enable a user to pick a contour of vessel wall on a two-dimensional DICOM image or to obtain an isosurface in a three-dimensional data set. The contour lines and the image intensity level for isosurfaces are decided on a trial and error basis with the consideration of maximum gradient of image intensity, *I*(*x*_{i}, *y*_{i}) or *I*(*x*_{i}, *y*_{i}, *z*_{i}). Unstructured surface and volume meshes were generated using Gridgen version 15.12 (Pointwise Inc., Fort Worth, TX), a mesh generation software. The faces of tetrahedral elements at the vessel walls are projected onto curved surfaces generated by an interpolation algorithm (SPHERIGON), which requires the position of the vertices and the corresponding normal vectors. For details of the algorithm, we refer the reader to Magnenat-Thalmann & Thalmann (2005).

Based on sensitivity tests of inlet boundary conditions (Oshima *et al.* 2005; Castro *et al.* 2006; Moyale *et al.* 2006; Baek *et al.* 2007), the inlet boundary conditions are imposed on the upstream of the lacerum ICA segment, which is shown in figure 1 (centre), to reduce the effect of a simplified boundary condition. Accordingly, long parent feeding vessels, which are not extended pipes but patient-specific ICA segments, are included in the computational domain to guarantee the development of the proper *secondary flow* upstream of aneurysms located in the supraclinoid segment. Small vessels less than 1 mm were not included because their effects are negligible (Cebral *et al.* 2005*a*) and all outlets were assumed to be efferent. The vessel walls were assumed to be rigid due to prohibitive computational cost of compliant vessel walls. The effect of compliant walls is discussed in §4.

The natural boundary condition is ∂**u**/∂*n* = 0 for velocity at outlets. Constant pressure has been used for the pressure boundary condition at the outlets, since the pressure measurement at the outlet vessels or the relation between the flow rate and pressure are not available in most cases. In cases with multiple outlets, however, constant pressure condition should be avoided because flow division is solely determined by the resistance of the downstream vessel segment in the computational domain. In this study, however, flow instability in aneurysms does not seem to be affected even by substantial flow rate changes in the downstream bifurcation; see appendix B for sensitivity studies with different outlet flow rates. Hence, constant pressure boundary condition is employed in most simulations of this study. ‘RC boundary condition’, a variant of the Windkessel model, is used for the outflow sensitivity test with patients B and C geometries. Specifically, the flow rates at the outlets such as MCA, ACA and PCoA are changed from those in constant pressure condition. In this type of boundary condition, the pressure is related to the flow rates through the differential equation *P*(*t*) + *RC* d*P*(*t*)/d*t* = *Q*(*t*)*R*, where *P*(*t*), *R*, *C* and *Q*(*t*) are the pressure at the boundary, the resistance, the capacitance and the flow rate, respectively (Grinberg & Karniadakis 2008).

For the pulse rate, 80 beats per minute were used in all simulations. At the inlet boundary, the Womersley velocity profile, which is the exact solution of pulsatile flow in a straight rigid pipe for a given flow rate, was imposed at the extended inlet. The flow rates either in the upstream or the downstream branches are not available. Hence, the VFR profile *Q*(*t*) at the inlet was taken from the measurements of Marks *et al.* (1992), and then approximated with *a*_{0} + ∑{*a*_{k} cos(2*π* *kt*/*T*) + *b*_{k} sin(2*π* *kt*/*T*)}, where *T* is the period of the cardiac cycle. Fourier coefficients were used to calculate the Womersley profile during simulations. The constant mode *a*_{0} represents the average VFR, and the pairs of *a*_{k} and *b*_{k} (*k* = 1, … , 7) normalized by constant mode *a*_{0} are (−0.152001, 0.129013), (−0.111619, −0.031435), (0.043304, −0.086106), (0.028871, 0.028263), (0.002098, 0.010177), (−0.027237, 0.012160) and (−0.000557, −0.026303). Patient-specific VFR may have higher frequency components, but including more terms showed negligible change in the flow-rate profile because of their small magnitudes compared to the average VFR. The pulsatility index in the VFR profile, i.e. the ratio of peak VFR to average VFR, is around 1.435 in the lower side of the range 1.66 ± 0.16 reported by Ford *et al.* (2005). The mean flow rate into the ICA ranges from 150 to 350 ml min^{−1}, which is in a physiologically reasonable range (Ford *et al.* 2005). Marshall *et al.* (2004) and Deane & Markus (1997) also reported that the averaged ICA VFR is 248.4 and 255 ± 42 ml min^{−1}, respectively.

With the Newtonian assumption of blood, the fluid motion is governed solely by the incompressible Navier–Stokes equations with the no-slip boundary condition on the fixed wall. The possible effect of non-Newtonian property of blood on the result is discussed in §4. The incompressible Navier–Stokes equations were solved with the parallel code NEKTAR (Karniadakis & Sherwin 2005), which implements the spectral/hp element method. NEKTAR employs two different formulations: (i) the Galerkin projection and (ii) the discontinuous Galerkin projection. The first one is similar to the standard finite element method (FEM) whereas the second one is similar to the finite volume method (FVM). Correspondingly, a semi-orthogonal and fully orthogonal Jacobi polynomial basis is used in the form of tensor products of order *p*. For *p* = 1 and 0 we recover the standard FEM and FVM, but the multi-resolution capability of NEKTAR allows us to resolve local flow features without changing the mesh but simply increasing the polynomial order *p*. A typical range for *p* is 5–10 to balance accuracy and computational efficiency. In terms of the mesh discretization, polymorphic elements consisting of hexahedra, tetrahedra, prisms and even pyramids can be employed. The time discretization is based on second- and third-order stiffly stable semi-implicit schemes. The high accuracy obtained with high spatial/temporal accuracy and a small time step *Δ**t* is very important in resolving flow instabilities. For more details of the method and implementation, we refer the reader to Karniadakis & Sherwin (2005). NEKTAR has been validated for many different bioflows (e.g. Sherwin *et al.* 2000). For the simulations in the present study, NEKTAR with the Galerkin projection is employed. We used the second-order semi-implicit scheme and Jacobi polynomial basis of order *p* = 4–8 on tetrahedral meshes. The polynomial order is increased when the mean flow rate increases and especially when the flow instability appears in order to confirm that the flow instabilities are not due to under-resolution. For case CL, systematic *hp*-refinement tests are performed, and convergence of WSS and velocity in *L*_{2}, and frequency spectra are documented in appendix A.

The specific numerical simulation parameters for all cases from case AL to case CH are listed in table 2. VFR ranges from 150 to 400 ml min^{−1}. We performed three, four and two simulations with different VFRs for patients A, B and C, respectively. Simulations with a fixed time step *Δ**t* = 0.001–0.0005 require 145 000–289 000 steps taking 15–20 h per cardiac cycle on parallel supercomputers utilizing 512–768 CPUs. Variations among cases are ascribed to the fact that computing time largely depends on the number of elements and non-dimensional period *TU*/*L*, where *U*, *T* and *L* are mean velocity, period and characteristic length, respectively.

The data from the third cardiac cycle were used for further analysis. We confirmed that the flow field reaches a stationary state within the first two cardiac cycles by calculating *D*_{V} = |*V*(*t*) − *V*(*t* − *T*)|_{∞}/|*V*(*t*)|_{∞} and *D*_{Q} = |*Q*(*t*) − *Q*(*t* − *T*)|_{∞}/|*Q*(*t*)|_{∞}, where *V*, *Q* and *T* are the components of velocity at a history point, VFR at an outlet, and period of the cardiac cycle, respectively; here |·|_{∞} denotes the maximum value. From the second to the third cycle, *D*_{Q} diminishes from *O*(1) to *O*(10^{−4}) and *D*_{V} from *O*(1) to *O*(10^{−4}) computed at some points located at the ICA and both aneurysms of patient C. To quantify the WSS dynamic fluctuations in the aneurysms as well as in the parent vessels, we compute the oscillating shear index (OSI) which is defined as , where we denote the WSS vector by ; is defined by *μ*(∇**u** + ∇^{T}**u**) · **n**, where *μ* and **n** are the dynamic viscosity and the inward normal vector at the blood vessel wall. The angle changes of WSS vectors on the vessel walls were also calculated, where the angle is defined as the amount of deviation from a reference vector. This reference vector at each point can be any unit vector on the tanget plane. However, in this study we chose a unit vector parallel to the WSS vector at the systolic peak. We also visualize the WSS vectors on local patches where OSI is high.

## 3. Results

Our computations provide a very large set of results, and only a small subset is presented here which demonstrate flow instability and WSS fluctuations. In cases AH, BH and CH listed in table 2, the aneurysmal flow becomes unstable with velocity fluctuations in the frequency range 20–50 Hz. Specifically, the flow in patient C became unstable with strong fluctuations in both velocity and WSS above VFR = 150 ml min^{−1}. For patients A and B, the flow instability became pronounced at higher VFRs, i.e. 250 and 400 ml min^{−1}, respectively.

### 3.1. Velocity profiles with inflection points

The presence of inflectional velocity profiles in arteries and associated instability was first hypothesized by Lighthill (1975). Here we document such velocity profiles in different patient-specific aneurysms. Specifically, complicated velocity profiles along the straight lines in aneurysms are observed and inflection points are found, which indicate the possibility of flow instability according to the Rayleigh inflection criterion (Drazin & Reid 1981). Specifically, figure 2 shows several streamlines and velocity profiles along the line A–B inside the aneurysm at the systolic peak; inflection points are present in the velocity profiles. For patient C, figure 3 shows profiles of the centreline velocity along the three lines inside the upstream aneurysm at two instants, i.e. at the systolic peak *and* during the decelerating systole. These graphs also show several inflection points. Here, the velocity profiles at two different instants look similar except that the overall amplitude changes. However, we note that the velocity profile at line ‘K’ shows increased negative velocity near the fundus of the aneurysm at instant (d); see the negative regions in the black lines labelled ‘K’ in figure 3*c*-1,*d*-1. We also checked the corresponding profiles during the diastole and they are closely similar.

### 3.2. Velocity–time traces

The flow instability develops during the deceleration phase resulting in high frequency fluctuations in velocity–time traces. Hence, velocity and pressure were recorded at up to 40 points upstream and downstream of aneurysms as well as inside the aneurysms to see any manifestation of instability; we call these points ‘history points’ in this paper. From velocity–time traces, their spectra are calculated and presented. Some velocity–time traces are shown in figure 4 from the simulations of cases AH and AX (VFR = 250, 150 ml min^{−1}), where case AX does not show strong fluctuations seen in case AH. Patient B exhibits a very stable flow until VFR increases from 200 up to 350 ml min^{−1} as shown in figure 5. Flow at VFR = 400 ml min^{−1}, however, demonstrates strong fluctuations as seen in the other patients A and C. Figure 6 shows velocity fluctuations in case CL (VFR = 150 ml min^{−1}) of patient C. The fluctuation amplitude becomes large between the systolic peak and the dicrotic peak, and decays towards the end of diastole. Figure 6 also shows similar trends in terms of velocity fluctuations, which start from the decelerating diastole and subsequently the flow returns to a smooth pulsatile state. Velocity traces both in the upstream and the downstream aneurysm show fluctuations. The velocity at ‘pt 1’ in the upstream aneurysm starts to fluctuate earlier than at ‘pt 10’ in the downstream aneurysm. The spectra at history points ‘pt 14’ and ‘pt 15’ are shown in figure 7*a* revealing some peaks at frequency *f*/*f*_{heart} around 16–50. The spectrum of patient B decays fast in figure 7*b* while the frequency components in the spectrum of patient C are strong up to *f*/*f*_{heart} = 100, usually a sign of a longer-developed structure in the velocity towards transition to a distinctly turbulent-like motion, but not yet truly turbulent. The changes in spectra over the VFR change are clearly shown in figure 7*c* with the cases AX, AL and AH of patient A. As the VFR increases from 150 to 250 ml min^{−1}, the decay of spectra becomes slower and peaks become more pronounced.

A history point ‘pt 1’ of patient C, which is located inside the upstream aneurysm, has strong peaks around *f*/*f*_{heart} = 22. The frequency ranges of the three cases AH, CL and CH are in good agreement with the Doppler recordings of Steiger & Reulen (1986) as well as the fact that flow instabilities start during the *deceleration phase* and that the flow returns to smooth pulsatile state either at the end of the systolic deceleration phase or at the end of diastole.

### 3.3. Oscillatory WSS vectors

Here we describe the oscillatory behaviour of WSS vectors and the spatial variation of WSS oscillations. We also report on the movement of the stagnation region inside the aneurysm. The direction changes of WSS vectors are compared with the WSS vectors at the systolic peak.

Owing to the velocity fluctuations we reported in §3.2, the WSS vectors on the aneurysm wall fluctuate as shown in figures 8 and 9 for patients A and B, respectively. Time traces of the magnitude and the direction of the WSS vectors at regions ‘L’, ‘M’ and ‘N’ on the aneurysm wall of patient A are plotted in figure 8. The regions ‘L’ and ‘M’ have high OSI, while the region ‘N’ has low OSI close to 0. The WSS vectors in high OSI regions, i.e. ‘L’ and ‘M’, show direction changes over 180° with large variation between two points at each region. In region ‘N’, however, the WSS vectors show similar fluctuations during the deceleration systole but with much smaller fluctuation amplitudes. Also, traces in region ‘N’ behave much more uniformly than regions ‘L’ and ‘M’ as shown in figure 8. In the bottom two rows of figure 8, WSS time traces of case AX at VFR 150 ml min^{−1} are shown for comparison of WSS between stable flow and unstable flow. Since the magnitudes are significantly different, it is hard to compare them directly. However, the WSS at higher VFR seem to be the scaled WSS of lower VFR superimposed with high frequency fluctuations.

We plot the time traces of the magnitude and direction of WSS vectors at the red points marked in figure 9*a*. The solid yellow lines indicate the locations of converging or diverging *stagnation lines* at the shown WSS vector field. Hence, figure 9*a* shows that the direction change of vectors along the solid yellow lines is minimal but transverse to these lines is maximal. More strikingly, visualizations at different times reveal that these lines are oscillating in space. The dotted yellow lines mark the boundaries of the spatial movement of the solid yellow lines over the cardiac cycle. Hence, the WSS vectors near these stagnation lines are experiencing significant oscillation over the cardiac cycle as their time traces plotted in figure 9*b* indicate. More specifically, the WSS vectors at points ‘L’ and ‘M’ show oscillations in both direction and magnitude. The reason for a smaller angle deviation near the systolic peak is that the angle changes of WSS vectors are computed with respect to the corresponding WSS vectors at the systolic peak. The WSS magnitude, however, shows larger oscillation at the systolic peak and decelerating systole similar to the velocity fluctuations. During the diastole, the WSS vectors change direction from 0° to ±180° while the magnitude does not change.

However, not all regions on the aneurysm wall are experiencing such a large oscillation in WSS magnitude and direction; figure 10 illustrates that the regions of high oscillation amplitude are very localized. In case CH of patient C, time traces of WSS magnitude and WSS vector angle at two points on the upstream aneurysm are plotted in figure 10*b*. The WSS distribution shown in figure 10*a* is a snapshot at the systolic peak. The WSS magnitude at the point marked with letter ‘M’ changes from 0 to 10 N m^{−2} over the cardiac cycle in phase with the flow-rate profile. High frequency oscillations with magnitude 2 N m^{−2} are superimposed on this slowly varying time trace of the WSS magnitude. The angle deviation, however, shows large oscillations from −160° to 180°. The WSS vector at the point ‘L’ demonstrates the opposite behavior; it changes its direction significantly less (±5°) during the cardiac cycle. The WSS magnitude, however, fluctuates rapidly with amplitude ±5–10 N m^{−2} due to the flow instability; it also follows the flow-rate profile varying from 10 N m^{−2} at the diastole to 25 N m^{−2} at the systolic peak.

Inside the aneurysm of patient B, the stagnation area exhibits a low WSS region surrounded by a band of high WSS as shown in figure 11*b*,*c*. The low WSS region has higher pressure than the pressure in the adjacent walls by 1 mm Hg. This pattern of WSS and pressure distribution agrees well with the experimental observations (Imbesi & Kerber 1999; Meng *et al.* 2007) in which it is pointed out that the high WSS region surrounds the impingement point. This pattern of high/low WSS is also noted from the upstream aneurysm of patient C as shown figure 10*a*. In this plot, the low WSS region at the aneurysm fundus looks larger than that of patient B, but the surrounding high WSS is clearly observed. More dramatic is the spatial change in direction of WSS vectors around this spot; all WSS vectors point outward in radial directions over the cardiac cycle. From this figure, it is clear that the flow impinges near the fundus of the aneurysm forming the stagnation point and creating a high WSS area around the impingement (stagnation) area. It is also interesting that this critical location moves around during the cardiac cycle; this is shown by the dots in figure 11*b*,*c*, which correspond to the lowest WSS spot at instants marked with dots of the corresponding colours. Therefore, a WSS vector at a fixed point near the stagnation area changes its direction ‘temporally’ due to change in the impingement location over the cardiac cycle.

### 3.4. Flow structure inside the aneurysm

Next we examine temporal changes of global flow structures inside the aneurysm such as the inflow/outflow regions and their changes along the depth are examined. We present the flow structure inside the aneurysm of case CH because case CH shows the strongest instability among the three patients. Figure 12 shows the distribution of velocity component normal (*a*–*d*) or parallel (*h*–*l*) to the plane on which data are extracted in case CH, VFR = 200 ml min^{−1}. The planes and time instants for the velocity computation are shown in figure 12*f*-1,2,*g*-1,2. The time interval between the time instants in figure 12*f*-2 is close to the half period of the dominant velocity fluctuation. For comparison, black lines in figure 12*h*–*l* are drawn as a reference, which is the border between the positive and negative regions in figure 12*i*.

Comparison of normal velocity distribution among figure 12*a*–*d* shows that negative (into the aneurysm) and positive (out of the aneurysm) regions do not change significantly over the cardiac cycle. Inside the aneurysm, the negative (inflow) region gradually moves toward the distal part of the aneurysm along the depth, following the aneurysm wall due to the inertia until it impinges on the distal dome of the aneurysm. Although the flow is unstable as demonstrated in the time traces of velocity at history points in figure 6, flow structures of positive and negative regions on planes over the cardiac cycle seem to change insignificantly. However, a closer look into the flow during the decelerating systole shows more changes as shown in figure 12*h*–*l* marked with black arrows, i.e. a patch of negative region appearing in the positive territory, and the countour lines crossing the black border line; these small scale changes are not significant enough to change the global flow structure. While the velocity fluctuates inside the aneurysm, pressure also fluctuates near the wall with peak-to-peak magnitude of about 100–200 N m^{−2} during the decelerating systole as shown in figure 13*c*. Pressure increases locally near the outer neck of the aneurysm wall, and a region of positive pressure difference *P*_{n} − *P*_{m} is clearly shown in figure 13*a* while the flow rate is decreasing.

## 4. Discussion

Our study shows that flow instability may occur in certain saccular aneurysms but not in all aneurysms. Specifically, the narrow-necked lateral aneurysm of patient B does not show velocity fluctuations seen in other patients until the VFR reaches 400 ml min^{−1} which is beyond the physiologically correct VFR. This flow instability does not seem to be affected by flow rate changes in the bifurcations and small vessels as documented in our sensitivity study in appendix B. Moreover, this study has not considered terminal aneurysms. Based on the fact that wide-necked aneurysms in patients A and C show much stronger velocity fluctuations from the VFR as low as 250 and 150 ml min^{−1}, we anticipate that flows in terminal aneurysms may be more susceptible to instability.

Oscillations reported herein seem to be induced by a hydrodynamic instability that is described by the Rayleigh inflection condition (Drazin & Reid 1981). Hydrodynamic instability and transition into turbulence in arterial pulsatile flow has been studied extensively in the past (Nerem *et al.* 1972; Steiger & Reulen 1986; Einav & Sokolov 1993; Blackburn & Sherwin 2007). More recently, high-resolution simulations were performed for stenosed carotids by Fischer *et al.* (2007), Lee *et al.* (2008) and Grinberg *et al.* (2009). In all cases, turbulence or temporal fluctuations were observed in the decelerating systole or near the end of systole. The observed frequencies, however, are much higher than the frequencies observed in the current work and intraoperative Doppler recordings of Steiger & Reulen (1986). The reasons for the discrepancy seem to be the different flow conditions such as the much higher Reynolds number 1000–2000 of the flow in the aorta or the stenosed carotid. In the numerical study of flow over periodic grooves in a channel by Ghaddar *et al.* (1986), the frequency of the self-sustained oscillation corresponding to the least unstable Tollmien–Schlichting mode was predicted from linear stability analysis. Our peak frequencies of patient C (around 24–58 Hz) seem to agree well with a frequency range (16–24 Hz) found in their study; this frequency range corresponds to frequencies in grooves with non-dimensional separation distance *L* = 5 and non-dimensional groove length *l* = 2.5 with different depth (non-dimensionalized with the half-channel width). The Reynolds numbers in our study based on the flow rate and the kinematic viscosity are 877 and 1271 on average and at the systolic peak, respectively. They are in good agreement with the critical Reynolds number for the onset of instability which is close to *Re* = 975 (Ghaddar *et al.* 1986).

Blood is a non-Newtonian fluid with shear-thinning and yield-stress. The effects of non-Newtonian property on the flow distribution in large arteries have been investigated and some studies show dramatic effects (Gijsen *et al.* 1999*a*,*b*) while others report a small effect limited to a fraction of the cardiac cycle (Perktold *et al.* 1991; Johnston *et al.* 2004). When the shear rates are higher than 50 or 100 s^{−1} (Long *et al.* 2004; Fischer *et al.* 2007), blood can be treated as a Newtonian fluid with constant viscosity. In our simulations, we assumed that the physiological flow rates in the ICA are in the range 200–250 ml min^{−1}, which gives an average velocity in the range of 0.265–0.331 m s^{−1}. This is based on the diameter of the ICA near the supraclinoid artery in our geometry, which is less than 4 mm, and the characteristic shear rate 8*V*_{a}/3*r*, 353–441 (s^{−1}), where *V*_{a} is the average velocity and *r* is the radius of the blood vessel. We confirmed that even inside the aneurysm, the shear rates are well above 100 s^{−1} even at the end of diastole. Our *a posteriori* analysis is shown in figure 14 with three non-Newtonian models and two different parameters per model, which are summarized in table 3. The characteristic shear rate where *D*_{ij} = 1/2(∂*u*_{i}/∂*x*_{j} + ∂*u*_{j}/∂*x*_{i}), *i*, *j* = 1, 2, 3 (Macosko 1994) is above 100 s^{−1} even at the end of diastole of a case with the smallest flow rate 150 ml min^{−1} except near the fundus of the downstream aneurysm and near the centre of upstream ICA as shown in figure 14*a*,*c*. For example, the shear rates at the black dots in figure 14*b*,*d* are 1352 and 1154.8 s^{−1} for patient A and C, respectively. We also calculate the ratio, *μ*/*μ*_{N}, of apparent viscosity (*μ*) using the three non-Newtonian models, i.e. the Carreau, Carreau–Yasuda and Ballyk models, to the Newtonian viscosity (*μ*_{N}) used in our simulations as shown in figure 14*b*,*c*. Except for the centre of the aneurysm of patient A and of the parent vessels of patient C, the shear rates are well above 200 s^{−1} and the apparent viscosity is close to or smaller than that used in the present study. Hence, we expect that the Newtonian assumption is acceptable in our simulations as it will not affect the onset of flow instability.

Our current study shows that the impingement locations inside aneurysms exhibit dynamically changing patterns of high or low shear and high static pressure as well as temporal oscillations in the direction of WSS vectors. Hence, the oscillatory behaviour of WSS vectors at such critical locations inside an aneurysm may be an important factor that should be considered for the pathology of an aneurysm. We have not attempted yet to establish a model which would represent the effects of flow on aneurysm development, or of systemic hypertension. In order to understand such pathological interactions between unsteady flow and aneurysms, the aggregation of red blood cells or adhesion of white blood cells and platelets to the aneurysmal wall has to be simulated in detail. However, the existing non-Newtonian models are inadequate in modelling the discrete nature of the blood rheology inside the aneurysm. Ultimately, we aim to assess the effects of non-Newtonian behaviour and elastic wall from ‘atomistic’ models (i.e. discrete at the scales of active entities such as cell adhesion molecules) and not from continuum-level phenomenological models. Hence, it seems necessary to develop a new multi-scale methodology to simulate the interaction of the unsteady flow with the blood cells and quantify their aggregation, accumulation or adhesion to endothelial wall inside the aneurysmal cavity. To this end, we are developing a hybrid continuum–atomistic formulation on two overlapped domains, with the latter employed in a small domain that includes the aneurysm while the former is employed in a much larger domain to represent accurately the large-scale flow dynamics. The atomistic formulation is based on the dissipative particle dynamics method whereas the continuum method is based on the incompressible Navier–Stokes equations using the spectral element method (Fedosov & Karniadakis 2009).

As regards the *assumption of rigid walls*, it can possibly affect the magnitude of WSS and pressure and/or the spectrum of velocity/WSS fluctuations. From the reports of Torii *et al.* (2007) and Oubel *et al.* (2007), we expect that our simulations with rigid walls may overestimate somewhat the WSS magnitude. However, as Torii *et al.* (2007) demonstrated that the effect of the deformable walls in a lateral aneurysm case is not as significant as in the terminal aneurysm case, we estimate that wall compliance changes only slightly the flow structure and WSS distribution in our case. We have also performed two-dimensional simulations with compliant walls and found that the velocity–time traces change by less than 10 per cent over a 10-fold change in the elastic modulus of the wall. However, we have not performed full three-dimensional simulations and the rigid wall assumption is a limitation of this study requiring further investigation.

## 5. Conclusions

Through a series of highly resolved CFD simulations, aneurysmal flow in wide-necked saccular intracranial aneurysms is shown to exhibit instabilities in the VFR 150–250 ml min^{−1}, which is lower than the mean VFR of 275 ml min^{−1} with standard deviation ±52 ml min^{−1} (Ford *et al.* 2005). Flow in a narrow-necked lateral aneurysm becomes unstable as the VFR in the ICA reaches 400 ml min^{−1}. Specifically, we observed that the pulsatile flow in aneurysms is subject to a hydrodynamic instability during the decelerating systolic phase resulting in a high-frequency oscillation over that period (Steiger & Reulen 1986). The flow returns to its original pulsatile state near the end of diastole. The frequency of velocity fluctuations in our simulations ranges from 20 to 50 Hz, and is in good agreement with the Doppler frequency recordings (Steiger & Reulen 1986). Inside saccular aneurysms, the pattern around the impingement areas shows relatively high static pressure and low WSS region surrounded by a band of higher WSS. At such locations, the WSS vectors are spatially changing significantly in their direction. Owing to flow instability, the WSS vectors fluctuate in magnitude and direction. In stagnation regions, the pressure/WSS patterns also fluctuate and move around spatially. The temporal change in these locations increases the temporal variation of WSS vectors. The oscillatory behaviour of WSS vectors, in combination with high/low WSS magnitude and locally elevated static pressure in stagnation areas, may play an important role in the formation and rupture of an aneurysm.

## Acknowledgements

We would like to acknowledge support by the NSF OCI grants 0845449, 0636336 and 0904288. Simulations were performed on the Cray XT3 (Bigben) at PSC, the Sun constellation linux cluster (Ranger) at TACC and the Cray XT5 (Kraken) at NICS.

## APPENDIX A

#### A.1. Resolution studies *(*‘hp-refinement’ test*)*

We performed a systematic ‘*hp*-refinement’ test under the same conditions used for our simulations. The time-dependent flow at the anatomically correct geometric model is tested at four different polynomial orders and four different meshes at fixed polynomial order *p* = 4. This study is intended to make sure that the oscillation is not due to numerical artefacts caused by under-resolution and to confirm that our NEKTAR code can resolve high frequency oscillation caused by the hydrodynamic instability. Hence, we selected a case in which the velocity oscillates during the decelerating systole.

We used patient C's two-aneurysm model. This geometric model was chosen because our study showed that the flows in the two adjacent aneurysms are prone to become unstable. The type of element used for space discretization is tetrahedral. For *p*-refinement test, the tested polynomial orders are *p* = 4, 6, 8, 10 and the grid points used in the calculations are (*p* + 3)(*p* + 2)^{2} = 252, 576, 1100, 1872 per element, respectively. Since the total number of elements is 89 324, the total grid points per variable for polynomial order *p* = 10 is 167, 214, 528. The effective distance between grid points in this case is about 0.0091–0.0364 mm in the aneurysm, which is 110–440 times smaller than the diameter of the inlet. For *h*-refinement, the number of elements are 81 162, 103 693 and 121 807 with a polynomial order 4 and the reference was taken from the simulation with the largest number of elements 121 807 and a higher order of polynomial ‘*p* = 6’.

#### A.2. h-refinement

We checked convergence in the simulation using a specific metric, the temporal average of spatial *L*_{2} error in WSS, against the results from a simulation with the finest mesh as a reference using equation (A 1). The frequency spectra of pointwise velocity are also compared among the results from different meshes.

The formula for the WSS error calculation is defined as
A 1
where the surface integration is done on the curved surface by numerical quadrature and the reference WSS (*p* = 6, number of elements 121 000) is denoted by .

The frequency spectra are obtained from the time traces at history points close to the points used in *p*-refinement test shown in figure 6 and a point in the upstream parent vessel during the entire third cardiac cycle. The WSS distribution on a patch shown in figure 15*a* is calculated from the velocity field during the entire third cycle.

*Convergence and spectra change*: the *L*_{2} errors of WSS from the *h*-refinement test show linear convergence over the number of elements or the size of element edge as shown in figure 15*b*,*c*. Spectra from different meshes show similar peaks and the changes are insignificant among them as shown in figure 15*d*.

#### A.3. p-refinement

We checked convergence in the simulation using two metrics: the *L*_{2} error in time of pointwise velocity and temporal average of spatial *L*_{2} error in WSS. The errors are compared against the results from a simulation with the highest order polynomial, i.e. *p* = 10, as a reference.

Simulations for error computations start from the systolic peak and continue until the dicrotic peak for the aforementioned reasons. This time interval for error computation is about one-sixth of the cardiac cycle as indicated by the grey box in figure 6.

The history points and the time interval during which temporal averages were taken are shown in figure 6. A patch on which WSS errors are computed is shown in figure 15*a*. The formula for WSS error calculation is equation (A 1), where the surface integration is done on the curved surface by numerical quadrature and the reference WSS (*p* = 10) is denoted by . The velocity error at each history point is taken as
A 2
where the reference velocity (*p* = 10) is denoted by .

We performed a simulation with polynomial order *p* = 6 up to four cardiac cycles. During the last cycle, the check points files were saved every time unit with time step *Δ**t* = 0.00025. Using the check point file at the systolic peak of the fourth cycle, three more simulations with different polynomial orders, i.e. *p* = 4, 8, and 10, were performed with the same time step, which satisfies the CFL condition for all polynomial orders.

*Convergence*: time traces of velocity at two history points are plotted in figure 15*f* showing no discernible difference. The *L*_{2} errors of WSS and velocity from the *p*-refinement test show exponential convergence over the order of polynomial as shown in figure 15*g*,*h*, respectively. Average WSS over the patch is about 10 N m^{−2} and the error changes from 0.2956 N m^{−2} (*p* = 6) to 0.0591 N m^{−2} (*p* = 8). Hence, they are about 3 and 0.6 per cent of the average WSS. For velocity, errors are even smaller: 0.1 per cent (*p* = 6) and 0.01 per cent (*p* = 8). WSS has one order lower accuracy because WSS is calculated through differentiation of field variables.

## Appendix B. Effect of outflow conditions

Figure 5 shows quantitative changes in the flow when flow rates at the outlets change by a factor of 5 in small arteries and by 20 per cent in large outlets of patient B. Specifically, using RC boundary conditions, the flow rates at the PCoA and the ophthalmic artery are reduced by factor of 5, and the increase in MCA flow rate compensates the flow reduction in the small arteries. The velocity–time traces at the upstream and downstream history points are plotted and compared between constant pressure (CP) and RC-type (RC) outflow boundary condition. Velocities at more points are examined as well as the points in figure 5. The zoom-in plots of time traces at pt 11, 15, and 18 show that the frequencies of small fluctuations do not change although the magnitudes of fluctuation change somewhat at some history points. The peak at pt 1 seems to disappear due to the flow rate increase in the MCA. However, the flow rate changes seem not to affect the flow instability.

Figure 16 shows flow rate changes in the MCA/ACA. In this case, we apply two different RC boundary conditions, one of which (RC1) makes the flow rates close to those corresponding to the constant pressure condition. By imposing the other RC outflow boundary (RC2) condition, we increase the VFR at the MCA and reduce the flow at the ACA by 30 per cent. Upstream history points do not show any high frequency oscillations seen in velocity–time traces at history points inside the aneurysms. Velocity–time traces at history points inside the aneurysms and upstream of aneurysms are almost the same among the different boundary conditions. Those downstream of the aneurysms show different magnitudes but similar fluctuations. Hence, we conclude that outlet flow rates do not seem to affect aneurysmal flow instability.

## Footnotes

↵1 In these simulations, about 200 time steps per cardiac cycler were used. Hence, possible temporal fluctuations were suppressed (J. R. Cebral 2009, personal communication).

- Received October 28, 2009.
- Accepted November 24, 2009.

- © 2009 The Royal Society