Computational and Mathematical Methods in Medicine
Volume 2016 (2016), Article ID 7412926, 16 pages
http://dx.doi.org/10.1155/2016/7412926

Research Article

Effects of Reynolds and Womersley Numbers on the Hemodynamics of Intracranial Aneurysms

Hafez Asgharzadeh and Iman Borazjani

Department of Mechanical and Aerospace Engineering, University at Buffalo, The State University of New York, Buffalo, NY, USA

Received 3 May 2016; Revised 29 July 2016; Accepted 10 August 2016

Academic Editor: Giancarlo Ferrigno

Copyright © 2016 Hafez Asgharzadeh and Iman Borazjani. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The effects of Reynolds and Womersley numbers on the hemodynamics of two simplified intracranial aneurysms (IAs), that is, sidewall and bifurcation IAs, and a patient-specific IA are investigated using computational fluid dynamics. For this purpose, we carried out three numerical experiments for each IA with various Reynolds ( to ) and Womersley ( to ) numbers. Although the dominant flow feature, which is the vortex ring formation, is similar for all test cases here, the propagation of the vortex ring is controlled by both and in both simplified IAs (bifurcation and sidewall) and the patient-specific IA. The location of the vortex ring in all tested IAs is shown to be proportional to which is in agreement with empirical formulations for the location of a vortex ring in a tank. In sidewall IAs, the oscillatory shear index is shown to increase with and because the vortex reached the distal wall later in the cycle (higher resident time). However, this trend was not observed in the bifurcation IA because the stresses were dominated by particle trapping structures, which were absent at low in contrast to higher .

1. Introduction

The rupture of intracranial aneurysms (IAs) is highly associated with mortality and morbidity [1]. Hemodynamics has a significant role in the growth and rupture of IAs [2–4]. Among the hemodynamic factors, vortical structures determine the complexity and stability of the flow pattern in an IA dome, which plays an important role in the rupture of IAs [3, 5–7]. Computational fluid dynamics (CFD) holds an important position in the investigation of hemodynamic factors in aneurysms because of its higher resolution near the walls relative to experimental methods such as laser Doppler velocimetry, particle image velocimetry, and magnetic resonance imaging, which is required to compute hemodynamic factors such as shear stress correctly [8, 9]. Many investigations have been carried out on the hemodynamics of IAs using experimental methods [10–12] and CFD [7, 13, 14].

Reynolds and Womersley numbers (explained in Section 3) are the only two nondimensional parameters required for full dynamic similarity in pulsatile internal flows [15]. Therefore, the investigation of their effects on the hemodynamics of IAs is strongly required. However, contradictory conclusions have been made in the literature on this topic. In fact, Jiang and Strother [16] concluded that increase of Womersley number can significantly increase the complexity of the flow pattern and vortex structures in two patient-specific intracranial aneurysms based on their CFD simulations. In contrast, Le et al. [17] stated that Womersley number does not affect the flow feature and structures based on their CFD simulations of an IA from a rabbit. Furthermore, Gopalakrishnan et al. [18] stated that while Womersley number does not change the vortex mode, but high Womersley number is associated with weak vortex rings in their simulation on abdominal aortic aneurysms. Bouillot et al. [10] concluded that has a negligible effect on the flow structures of an idealized sidewall intracranial aneurysm for a steady inflow according to their PIV measurements. A similar conclusion was made by Le et al. [17] and Cebral et al. [19] based on their CFD simulations of cerebral aneurysms in rabbits and humans, respectively. In contrast, Gopalakrishnan et al. [18] stated that, by increasing , the strength of the main vortex structure increases.

Previous works [7, 19] have shown the significant effect of aneurysm geometry on the aneurysm hemodynamics, which is not investigated further here. Our aim is to compare the effects of the Reynolds and Womersley numbers on the hemodynamics of both sidewall and bifurcations IAs by keeping other parameters, for example, inlet flow waveform and geometry, constant. Furthermore, their effect on the hemodynamics is formulated using dimensional analysis and compared with the literature.

2. Governing Equations and the Numerical Method

In this section, Einstein’s tensor notation, where repeated indices imply summation, is used unless otherwise indicated (, , ). The governing equations are the , unsteady incompressible Navier-Stokes equations for a Newtonian fluid in curvilinear coordinates as follows [20]:

where , , and are the contravariant velocity, curvilinear coordinate, and Cartesian coordinates components, respectively. and are the nondimensional pressure and time, respectively. is the Jacobian of the geometric transformation, , and are the metrics of the transformation. is the nondimensional Cartesian velocity, is the contravariant metric tensor, , and is the Reynolds number of the flow based on characteristic length and velocity scales.

We use curvilinear/immersed boundary (CURVIB) and overset grid methods, which are extensively described and validated [20]. A sharp-interface immersed boundary method is used to handle the , arbitrary complex boundaries (IA geometry in this study) inside the flow domain [21]. The nodes that are outside the flow domain are blanked out and do not affect the solution. These nodes are identified using an efficient ray-tracing algorithm [22]. The boundary conditions are reconstructed on the fluid nodes in the immediate vicinity of the immersed boundary along the normal direction to the boundary [21]. The method has been shown to be second-order accurate for a variety of flows [21, 23]. The overset grid approach is implemented to reduce wasted nodes in a domain, which are blanked out by the immersed boundary method [20]. In this approach, a complex flow domain is divided into several arbitrary subgrids with overlaps. To solve the governing equations at each subgrid, boundary conditions at the interfaces are constructed by interpolation from host subgrid. The details of the overset-CURVIB method can be found in [20]. The method has been validated against experimental and benchmark solutions [20, 24] and has been applied to a variety of problems such as cardiovascular flows [17, 25–28], aquatic swimming [23, 29], and rheology [30]. Furthermore, we have validated our method for flows inside an immersed body by comparing our results with the measurements of the pulsatile flow through a bend in Appendix. As shown in Appendix, the computational results are in excellent agreement with our previous simulations using body-fitted grids [20] and the experimental results [31].

We have assumed rigid walls similar to previous simulations [2, 14, 17, 19] because the displacement of aneurysm’s wall is typically small and the flow patterns of small distensible and rigid models in the carotid region are very similar [32]. In addition, we have assumed Newtonian fluid in our simulations because the non-Newtonian effects are negligible in larger (>500 μm) arteries [33] and previous simulations of both Newtonian and non-Newtonian fluids have shown similar flow patterns [19].

3. Description of Simulated Test Cases

Numerical simulations have been carried out on two simplified geometries (sidewall and bifurcation), which can be considered as simplified models of IA, and a patient-specific IA geometry. Figure 1 shows the configuration of the simplified models for (a) sidewall and (b) bifurcation IAs. The inlet and outlet of the geometries are constructed from a pipe with the diameter . The aneurysm dome is modeled by a cone-like shape with elliptical base and locus, whose radiuses are , , and .

(a)

(b)

Figure 1: Schematic illustrating of IA simplified models (a) sidewall and (b) bifurcation IA, where ellipse radiuses are , , and ). denotes the location of the vortex ring from the proximal wall and inlet/outlet cross section in the sidewall and bifurcation IA, respectively.

Figure 2(a) shows the immersed body and overset grids layout for sidewall type IAs and Figure 2(b) shows the immersed body and a Cartesian grid layout for the patient-specific IA geometry. The inlet and outlet(s) in simplified models are meshed by body-fitted curvilinear grids; the circular base is meshed by grid nodes and the axial (flow) direction is meshed by grid nodes. The geometry of the aneurysm is placed as an immersed boundary onto the uniform grid, and all subdomains (Figure 2(a)) are solved simultaneously [20]. The size of the domain that contains the immersed body is , , and and , , and for sidewall and bifurcation simplified IAs, respectively. The uniform grid is meshed by and nodes for sidewall and bifurcation simplified IAs, respectively. Such resolution was found to be enough to obtain grid independent solutions based on our extensive grid refinement studies. Furthermore, the results obtained using a similar resolution for the pulsatile flow through a bend show good agreement with experimental results (Appendix). The size of the domain that contains the immersed body of a patient-specific geometry is , , and , where ( m) is the hydraulic diameter of the cross section of the inlet vessel and Cartesian grid layout. The uniform grid for the patient-specific geometry is meshed by nodes, which has a similar grid resolution with the simplified IAs.

(a)

(b)

Figure 2: Immersed body and overset grids layout for (a) sidewall IA simplified model and (b) patient-specific IA geometry. The immersed body is meshed with triangular elements.

The Neumann boundary condition is applied at the outlet boundaries, and the volumetric flow rate of each outlet is calculated based on the principle of optimal work [36]. The inlet boundary condition is a uniform (plug) flow in space but pulsatile in time with a waveform as shown in Figure 3 with pulsatility index (where is the difference between the peak systolic and minimum diastolic velocities). The pulsatile velocity waveform is from a patient with a cerebral aneurysm [34]. The numerical simulations were carried out for three cycles to obtain quasi-steady results, which are not affected by the initial condition.

Figure 3: The inlet velocity waveform during one heartbeat. Adapted from [34].

The Reynolds number is defined as , where is the average bulk inlet velocity and  m and  m are the inlet diameter for simplified IAs and the patient-specific IA, respectively, and   is the blood kinematic viscosity (calculated based on and , where and are the blood density and dynamic viscosity, resp. [7]). The Womersley number is defined as , where is the period of the waveform calculated based on the heart rate. The simulations on the sidewall and bifurcation IAs are carried out with two Womersley and two Reynolds numbers. The Womersley number is modified by changing the heart rate and Reynolds number is modified by altering the inlet bulk velocity in the physiologically relevant range [37]. The nondimensional time-step () is calculated by dividing the nondimensional into equal time instants for all simulated cases. Table 1 shows the simulated cases, their specifications, and denoted names. It should be noted that the calculated and in Table 1 are within the physiological range in IAs ( and [7, 14, 38]).

Table 1: Specifications of simulations. H, high; L, low; , diameter; , heart beat period; Re, Reynolds number; Wo, Womersley number; , the average bulk inlet velocity; An, aneurysm number; , the nondimensional time-step.

Simulation
IA type
(m)
Heart beats per minute
(s)
Wo
(m/s)
Re
An

ReLWoL
Simplified sidewall
0.005
75
0.8
7.71
0.1
151.51
1.43
0.002

ReHWoL
Simplified sidewall
0.005
75
0.8
7.71
0.25
378.79
1.43
0.005

ReHWoH
Simplified sidewall
0.005
125
0.48
9.96
0.25
378.79
1.43
0.003

ReLWoL
Simplified bifurcation
0.005
75
0.8
7.71
0.1
151.51
2.86
0.002

ReHWoL
Simplified bifurcation
0.005
75
0.8
7.71
0.25
378.79
2.86
0.005

ReHWoH
Simplified bifurcation
0.005
125
0.48
9.96
0.25
378.79
2.86
0.003

ReLWoL
Patient-specific
0.0048
75
0.8
7.4
0.1
145.45
2.4
0.002

ReHWoL
Patient-specific
0.0048
75
0.8
7.4
0.25
363.64
2.4
0.005

ReHWoH
Patient-specific
0.0048
125
0.48
9.56
0.25
363.64
2.4
0.003

(aneurysm number) in Table 1 is the ratio of the transport to vortex formation time scales in IAs [17]. The transport time scale in sidewall IAs is the time it takes for the parent artery flow to transport a fluid particle across the IA neck and the transport time scale in bifurcation IAs is the time it takes for the inlet artery flow to transport a fluid particle to the outlets. can be defined as follows:

where is the aneurysm neck and the outlet diameter in sidewall and bifurcation IAs, respectively, defined based on the transport time scale. and is determined based on the transport time scale definition for sidewall and bifurcation IAs, respectively. If vortex formation time scale is smaller, that is, , the vortex is formed before it is advected out of the dome area. In contrast, if the flow is dominated by a stationary shear layer from the proximal wall to the distal wall through the IA neck.

4. Results and Discussion

In this section, the effects of Reynolds and Womersley numbers on the flow patterns and vortical structures of IAs are investigated. Because of the substantial difference between vortical structures in sidewall and bifurcation IAs, they are investigated in different sections, that is, Sections 4.1 and 4.2 for sidewall and bifurcation IAs, respectively. In Section 4.3, whether observed effects of Reynolds and Womersley numbers (on flow patterns and vortical structures) in simplified IAs can be seen in a patient-specific (complex) IA is explored. In Section 4.4, the location of the vortex ring in simplified IAs and the patient-specific IA is investigated using the nondimensionalization method. Finally, the effect of and on the oscillatory shear index and wall shear stress of all test cases is investigated in Section 4.5. The three-dimensional vortical structures are visualized by the isosurfaces of -criteria [35].

4.1. Sidewall IA Simplified Models

In this section, the effect of Womersley and Reynolds numbers on the flow patterns and vortical structures of simplified sidewall IAs is investigated. Figure 4 shows the time evolution of the nondimensional out-of-plane vorticity and in-plane velocity vectors for simplified sidewall IAs with different and . In Figure 4, the dominant flow patterns are similar for different and . In all cases a vortex ring, denoted as R1, starts to form from the proximal wall in the early acceleration phase (Figure 4(a)). Before the inlet velocity reaches the peak systole, the vortex ring is convected across the IA neck (Figures 4(b)–4(d)). Eventually it hits the distal wall and diffuses just after the peak systole and in the deceleration phase (Figures 4(e) and 4(f)). A similar behavior can be observed in Figure 5, which shows the time evolution of the three-dimensional vortical structures visualized by the isosurface of -criteria for simplified sidewall IAs with different and . In all cases the vortex ring forms at early acceleration phase (Figure 5(a)) and is convected across the IA neck (Figures 5(b)–5(d)) and hits the distal wall before the peak systole phase. Subsequently, it rolls up to form a recirculation vortical structure, denoted as VS1, which stays until the end of the cycle. The size of vortex ring for all cases is approximately similar (Figure 5). The for all the simplified sidewall cases is equal to , which is consistent with the previous work [11, 17, 26] stating that when the flow is dominated by the vortex formation.

Figure 4: Evolution of the nondimensional out-of-plane vorticity and in-plane velocity vectors for simplified sidewall IAs with different and ; (upper row) and ; (middle row) and ; and (lower row) and at various time instants during the cycle: (a) , (b) , (c) , (d) , (e) , and (f) . R1 depicts the vortex ring in simplified sidewall IA and the green cross shows the location of the maximum -criteria [35] on the midplane of IAs.

Figure 5: Evolution and topology of the three-dimensional vortical structures for sidewall IAs with different and ; (upper row) and ; (middle row) and ; and (lower row) and at various time instants during the cycle: (a) , (b) , (c) , (d) , (e) , and (f) .

The only difference in the vortex ring formation and evolution is the location of the vortex ring for different and (Figures 4 and 5). To further investigate the effect of and , the location of the center of the vortex ring at different time instants in the cycle is identified by the maximum -criteria on the midplane of simplified sidewall IAs and is shown by the green crosses in Figure 4. As it can be observed by comparing different rows of Figures 4(a)–4(d), the location of vortex ring is different at specific time instants. The effect of on the convection of the vortex ring can be observed by comparing middle and upper rows of Figures 4(a)–4(d) at specific time instants; that is, the convection of the vortex ring corresponding to (middle row) is higher than (upper row) with similar because the distance of vortex ring from proximal wall is higher for (middle row) at similar in the cycle. The effect of on the convection of the vortex ring can be observed by comparing middle and lower rows of Figures 4(a)–4(d) at specific time instants; that is, the convection of the vortex ring corresponding to (middle row) is higher than (lower row) with similar . This behavior can be observed by the location of the nondimensional out-of-plane vorticity as well. The convection of the vortex ring in simplified sidewall IAs and its effect on hemodynamic stresses is discussed in detail in Sections 4.4 and 4.5, respectively.

4.2. Bifurcation IA Simplified Models

In this section, the effect of Womersley and Reynolds numbers on the flow patterns and vortical structures of simplified bifurcation IAs is investigated. Figure 6 shows the time evolution of the nondimensional out-of-plane vorticity and in-plane velocity vectors for simplified bifurcation IAs with different and . From Figure 6, the dominant flow patterns, that is, vortex formation, are similar in all simulations of simplified bifurcation IAs. This is consistent with our previous simulations with in terms of vortex formation. The vortex ring, denoted as R2, starts to form at inlet/outlet cross sections in the early acceleration phase (Figure 6(a)). Subsequently, the vortex ring grows in size, while it moves to the outlets (Figures 6(b)–6(d)). Finally in the peak systole, the complicated vortical structures are observed (Figure 6(e)).

Figure 6: Evolution of the nondimensional out-of-plane vorticity and in-plane velocity vectors for simplified bifurcation IAs with different and ; (upper row) and ; (middle row) and ; and (lower row) and at various time instants during the cycle: (a) , (b) , (c) , (d) , (e) , and (f) . R2 depicts the vortex ring in simplified bifurcation IAs and the green cross shows the location of the maximum -criteria on the midplane of IAs.

At peak systole more segments of vortical structures in the dome can be observed for the higher (middle and upper rows of Figure 6(e)) in comparison to the lower (the lower row of Figure 6(e)). This is the consequence of the vortex ring interactions with other vortical structures, which cannot be seen in two-dimensional vorticity on the midplane. In order to show the three-dimensional vortical structures, Figure 7 plots the time evolution of the three-dimensional vortical structures visualized by the isosurfaces of the -criteria for simplified bifurcation IAs with different and . Other than the vortex ring, another vortical structure can be seen at the inlet/outlet cross section close to the IA neck, denoted as VS2, for (middle and upper rows in Figures 5(a)–5(d)). VS2 is observed by Vigolo et al. [39] for the first time and later by Chen et al. [40] in -junction pipes and named “particle trapping vortical structures.” Particle trapping vortical structures are distinguished from other structures by the vorticity magnitude in normal direction to the midplane of models (-direction) in Figure 7, because particle trapping structures have almost zero vorticity in the normal direction to the midplane of the models. Because of the existence of an IA in the current study, the particle trapping is distorted in comparison to that in the -junction pipe [39].

Figure 7: Evolution and topology of the three-dimensional vortical structures colored by vorticity- magnitude for simplified bifurcation IAs with different and ; (upper row) and ; (middle row) and ; and (lower row) and at various time instants during the cycle: (a) , (b) , (c) , (d) , (e) , and (f) .

One effect of in simplified bifurcation IAs is that particle trapping vortical structures do not form at the acceleration phase of (the lower row of Figure 7). This behavior is also observed by Vigolo et al. [39] who reported that particle trapping vortical structures do not form for . As a consequence, more complex and highly asymmetric flow patterns are observed for (upper and middle row of Figure 7(e)) in contrast to because of the vortex ring interaction with the particle trapping vortical structure.

The green crosses in Figure 6 show the location of the center of the vortex ring by the maximum -criteria on the midplane of simplified bifurcation IAs at different time instants in the cycle. The effect of on the convection of the vortex ring can be observed by comparing middle and upper rows of Figures 6(b)–6(d) at specific time instants; that is, the convection of the vortex ring corresponding to (middle row) is higher than (upper row) at similar and . The effect of on the convection of the vortex ring can be observed by comparing middle and lower rows of Figures 6(b)–6(d) at specific time instants; that is, the convection of vortex ring corresponding to (middle row) is higher than (lower row) with similar at the same . This behavior is similar to what was observed in sidewall IAs in the previous section (Figure 4). The convection of the vortex ring in simplified bifurcation IAs and its effect on hemodynamic stresses is discussed in detail in Sections 4.4 and 4.5, respectively.

4.3. Patient-Specific IA

In this section, the effect of Womersley and Reynolds numbers on the flow patterns and vortical structures of a patient-specific IA, which is of sidewall type, is investigated. Figure 8 shows the time evolution of the nondimensional out-of-plane vorticity and in-plane velocity vectors for patient-specific IAs with different and . The plane is selected the way that vortical structures in the dome can be observed clearly (the plane is shown in Figure 2(b)). From Figure 8, the dominant flow patterns are similar to sidewall IAs (Figure 4). In all cases with different and a vortex ring, denoted as R1, starts to form from the proximal wall in the early acceleration phase (Figure 8(a)). The vortex ring evolves as it is convected across the IA neck (Figures 8(b)–8(d)). Finally it hits the distal wall and diffuses after the peak systole and in the deceleration phase (Figures 8(e) and 8(f)). A similar behavior can be observed from Figure 9 for different and , which shows the time evolution of the three-dimensional vortical structures visualized by the isosurface of -criteria for patient-specific IAs. In all cases the vortex ring forms at early acceleration phase (Figure 9(a)) and is convected across the IA neck (Figures 9(b)–9(d)) and hits the distal wall before the peak systole phase. for all cases is equal to , which is consistent with the previous work [11, 17, 26] stating that when the flow is dominated by the vortex formation. The only effect of and on the vortex ring formation and evolution is the location of the vortex ring at specific time (Figures 8 and 9).

Figure 8: Evolution of the nondimensional out-of-plane vorticity and in-plane velocity vectors for patient-specific IAs with different and ; (upper row) and ; (middle row) and ; and (lower row) and at various time instants during the cycle: (a) , (b) , (c) , (d) , (e) , and (f) . R1 depicts the vortex ring in patient-specific IAs and the green cross shows the location of the maximum -criteria on the midplane of IAs.

Figure 9: Evolution and topology of the three-dimensional vortical structures for patient-specific IAs with different and ; (upper row) and ; (middle row) and ; and (lower row) and at various time instants during the cycle: (a) , (b) , (c) , (d) , (e) , and (f) .

By comparing the flow structures in patient-specific IAs (Figure 9) with simplified sidewall IAs (Figure 5), it can be observed that more complex and broken vortical structures are formed in the patient-specific IA in comparison to the simplified sidewall IA, while the dominant flow pattern (vortex ring formation and evolution) is similar. The vortex ring in patient-specific IA is deformed because of the complexity of its geometry.

To further investigate the effect of and , the location of the center of the vortex ring at different time instants in the cycle is identified by the maximum -criteria on the plane (shown in Figure 2(b)) of patient-specific IAs and is shown by the green crosses in Figure 8. The similar behavior with simplified models (Figures 4 and 6) can be observed in the patient-specific geometry (Figure 8); that is, the convection of the vortex ring corresponding to (middle row) is higher than (upper row) with similar because the distance of vortex ring from proximal wall is higher for (middle row) at similar in the cycle. The convection of the vortex ring corresponding to (middle row) is higher than (lower row) with similar . The convection of the vortex ring in the patient-specific IAs and its effect of hemodynamic stresses are discussed in detail in Sections 4.4 and 4.5, respectively.

4.4. The Effect of and on the Location of the Vortex Ring

In this section, the effect of and on the convection of the vortex ring in simplified IAs (sidewall and bifurcation) and a patient-specific IA is investigated, and a formula for the location of the center of the vortex ring is derived using dimensional analysis. Figure 10 plots the distance of the maximum -criteria (location of the center of the vortex ring) on the midplane of simplified IAs (sidewall and bifurcation) and the plane of patient-specific IAs at different and . These distances () are calculated from the proximal wall for sidewall (in -direction shown in Figure 1(a)) and also from the proximal wall for patient-specific IAs (in -direction shown in upper row of Figure 8(f)) and from the inlet/outlet cross section for bifurcation IAs (in -direction shown in Figure 1(b)). is normalized by the diameter of the inlet parent artery () in all test cases. It can be observed that the vortex ring formation starts in and for simplified sidewall and bifurcation IAs, respectively, regardless of and . The location of the vortex ring is different in various and at similar instants in the cycle. Based on Figure 10 and previous discussion in Figures 4, 6, and 8, high and low are associated with faster convection of the vortex ring since the vortex ring corresponding to high and low reaches the distal wall sooner. This behavior can be explained by dimensional analysis. The dimensionless distance of the vortex ring from the proximal wall can be defined as

where and are the vortex convection velocity and time, respectively. is the distance of the vortex ring location from the proximal wall in simplified sidewall IAs (), also the proximal wall in patient-specific IAs (), and the distance of the vortex ring location from the inlet/outlet cross section in bifurcation IAs () as denoted in Figure 1. Here, we simply approximate the vortex ring convection velocity by assuming that it is equal to the bulk flow velocity in the parent artery. Applying the dimensional analysis in (3) leads to

where and . and definitions can be applied in (4), which results in

Based on (5), decreasing and increasing result in decreasing distance of the vortex ring from the proximal wall at specific instants (), which is in agreement with Figure 10. It is noted that () is the dimensionless velocity, which is similar for all simulations at the specific dimensionless time () regardless of and (as can be seen from Figure 3). The effectiveness of and on the normalized distance of the vortex ring from the proximal wall can be compared using (5) as follows:

The trend observed in (6) agrees well with that of Figure 10, meaning that and at specific instants (). The trend of the location of the vortex ring for various and agrees well with the formulation developed from experiments on thin core rings generated by a piston gun in water [41]. By substituting simulation specifications of current study on the formulations expressed in [41] and can be reached which shows promising agreement with the simple formulation obtained in this study (6), considering that in [41] the vortex forms in the tank instead of an IA neck.

(a)

(b)

(c)

Figure 10: The distance of the maximum -criteria [35] on the midplane for different and for (a) sidewall IAs, where the vortex ring distance is measured from the proximal wall in the -direction; (b) bifurcation IAs, where the vortex ring distance is measured from the inlet/outlet cross section in -direction; (c) patient-specific IAs, where the vortex ring distance is measured from the proximal wall in the -direction at various time instants during a part of cycle that vortex ring forms and breaks down. This distance is normalized by the IA width (), distance between proximal and distal wall (), and outlet diameter () in simplified sidewall, patient-specific, and simplified bifurcation IAs, respectively. Values for ReHWoH, ReHWoL, and ReLWoL can be found in Table 1 for different geometries.

4.5. The Effect of and on the Oscillatory Shear Index and Wall Shear Stress

The oscillatory shear index () and wall shear stress () are important parameters in determining aneurysm rupture risk [2, 42, 43]. In this section, we investigate the effect of and on and in terms of the dynamic behavior of the vortex ring. Figures 11 and 12 show the cycle-averaged distribution of oscillatory shear index and normalized wall shear stress, respectively, with different and (Table 1) for (a) simplified bifurcation, (b) simplified sidewall, and (c) patient-specific IAs. For the purpose of comparison, WSS of the IA dome is normalized by that of the inlet parent artery similar to previous work [2]. To quantify the effect of and on the hemodynamic stresses, normalized and are averaged over the dome surface. The results are plotted in Figure 13, which shows oscillatory shear index and normalized wall shear stress for different and on simplified bifurcation, simplified sidewall, and patient-specific IAs.

Figure 11: Cycle-averaged distribution of oscillatory shear index with different and for (a) simplified bifurcation, (b) simplified sidewall, and (c) patient-specific IAs. The oscillatory shear index is calculated based on the cycle-averaged OSI. Values for ReHWoH, ReHWoL, and ReLWoL can be found in Table 1 for different geometries.

Figure 12: Cycle-averaged distribution of normalized wall shear stress with different and for (a) simplified bifurcation, (b) simplified sidewall, and (c) patient-specific IAs. is normalized by that of the inlet parent artery. Values for ReHWoH, ReHWoL, and ReLWoL can be found in Table 1 for different geometries.

(a)

(b)

Figure 13: (a) Oscillatory shear index and (b) normalized wall shear stress for different and at simplified bifurcation, simplified sidewall, and patient-specific IAs. and are averaged on the dome surface. The resulting is normalized by that of the parent artery. Values for ReHWoH, ReHWoL, and ReLWoL can be found in Table 1 for different geometries.

It can be observed in Figure 13(a) that in the simplified sidewall IA for low (lower row of Figure 11) is higher than that of high (middle row of Figure 11) with similar . In addition, for high (upper row of Figure 11) is higher than that of low (middle row of Figure 11) with similar in the simplified sidewall IA. A similar trend can be seen in the patient-specific IA for ; that is, for low (lower row of Figure 11) is higher than that of high (middle row of Figure 11) with similar and for high (upper row of Figure 11) is higher than that of low (middle row of Figure 11) with similar . This is due to the fact that low and high increase the vortex residence time (the time that vortex ring hits the wall and breaks down minus the time that vortex ring starts to form, normalized by the period of the cycle) in the dome of sidewall IAs (Figures 10(a) and 10(c)), which increases the flow disturbance in the cycle and, consequently, results in higher .

For the simplified bifurcation IA, in contrast to sidewall IAs, it can be observed in Figure 13(a) that corresponding to is slightly () higher than with similar . corresponding to is lower than in simplified bifurcation IA with similar because of the absence of “particle trapping” vortical structures at acceleration phase of , which decreases the segments of vortical structures (by comparing lower and middle row of Figure 7) and flow disturbance.

It can be observed in Figure 13(b) that and variations change the normalized by only a maximum of among all IAs. Therefore, we conclude that the effect of and variations on the normalized is severalfold smaller than that of .

5. Conclusions

The effect of two key parameters, that is, Reynolds and Womersley number, on the hemodynamics of simplified IAs (sidewall and bifurcation) and patient-specific IAs is investigated using CFD simulations. Based on our results, the dominant flow pattern and the vortex structure, for example, the vortex ring formation, remain similar by modifying from to and from to in both simplified bifurcation and sidewall IAs. Similarly, the vortex ring formation remains similar by modifying from to and from to in patient-specific IAs. However, the location of the vortex ring at different instants in a cycle in both simplified and patient-specific IAs is controlled by and . In addition, the interaction of the vortex ring with other vortical structures depends on in bifurcation IAs.

The location of vortex ring as function of time depends on the combination of and . We found that the high and low are associated with the fast convection of the vortex ring in IAs. Using dimensional analysis, a formulation (Equation (5)) is obtained which can clearly demonstrate the trend of vortex ring distance from the proximal wall in simplified sidewall and patient-specific IAs and from the inlet/outlet cross section of simplified bifurcation IAs. Based on the obtained formula, the location of vortex ring is proportional to which shows good agreement with formulations expressed in [41] for the vortex ring location in a tank.

The highly asymmetric and complicated vortical structure observed in bifurcation IAs at in comparison to are a consequence of the interaction between the vortex ring and the particle trapping vortical structures. More organized vortical structures are observed for since the particle trapping vortical structure does not form in the acceleration phase at low ( according to [39]). Therefore, there is no interaction between the vortex ring and the particle trapping vortical structures in bifurcation IAs with .

We found that variations of and slightly affect the normalized (maximum of ), while is proportional to and in sidewall IAs (both simplified and patient-specific). The observed trend of in these IAs is a consequence of higher residence time of the vortex (because the vortex reaches the distal wall later in the cycle at high and low ), which disturbs the flow for a longer time in each cycle and increases the . This trend is not observed in bifurcation IAs because the ring vortex is convected toward the outlets and does not enter the dome. The shear stress in bifurcation IAs is dominated by the particle trapping structures. In fact, the observed lower () in bifurcation IAs at in comparison to is because of the absence of particle trapping structures at . Increasing from to only slightly increased () at similar in simplified bifurcation IAs.

Based on our results, the hemodynamics of the simplified and the patient-specific sidewall IAs is similar in terms of vortex formation, propagation, normalized , and . The bifurcation IA showed different trend relative to the sidewall ones. Note that the conclusions of this work are based on two idealized geometries and on two idealized and one patient-specific geometries for two Re and Wo. We believe our results are valid for other geometries in the physiological range of and for aneurysms. Nevertheless, this will be tested in a cohort of aneurysms rigorously in the future.

Appendix

Validation: Pulsatile Flow through a 90° Bend

The three-dimensional pulsatile flow through a strongly curved pipe bend is simulated to validate our method, which is used to simulate the test cases in this study, for a pulsatile flow inside an immersed body. The geometry of the test case is shown in Figure 14(a). As shown in Figure 14(a), the radius of the curvature of the bend is three times the pipe diameter (). The outlet is placed after the bend and the inlet is placed before the bend. A gear pump providing a steady flow with the Reynolds number equal to (, where is the bulk velocity, is the pipe diameter, and is the kinematic viscosity) in conjunction with a piston pump generating sinusoidal flow waveform with to was used to generate the pulsatile flow waveform in the experiments [31]. The resulting Womersley number of the flow in the experiments was .

(a)

(b)

Figure 14: (a) The geometry for the bend (b) incoming flow waveform. Solid line: analytical nondimensional bulk velocity waveform prescribed at the inlet of the flow domain obtained by solving the Womersley equation; circles: the measurements of nondimensional inlet bulk velocity from experiments [31]. The immersed body is meshed with triangular elements.

To generate a waveform that matches the experiment’s waveform, the inlet boundary condition is the velocity profile set from the Womersley solution of a fully developed pulsatile flow in a circular pipe as follows [44]:

where is the zero-order Bessel function of the first kind, is the radius of the pipe, is the radial distance from the center of the pipe, is the angular frequency of the flow oscillation, and is the flow viscosity. The constant is selected to be to generate a sinusoidal flow waveform, which matches the experiment’s waveform [31, 44]. The above equation is solved using MATLAB, and the resulting solutions are stored and fed into the solver to specify the time-varying inlet flow (Figure 14(b)). As shown in Figure 14(b), the computed inlet waveform is in reasonably good overall agreement with the experimental inlet . The outlet boundary condition is the Neumann boundary condition. The time period of inflow oscillations, which is nondimensionalized based on and , is . The nondimensional time-step is used for the simulations. The Cartesian domain is discretized by grid nodes. The region around the bend is discretized by uniform grid nodes and the grid is stretched to the boundaries.

The simulated flow in the bend with the pulsatile inlet is compared with the experimental results [31] in Figure 15. The streamwise velocity profiles at the symmetry plane are plotted at five different locations (i.e., , , , , and ) for four different time instants during the cycle (, , , and ). The results in Figure 15 are in excellent agreement with the computational results of the same setup with body-fitted grids [20]. Furthermore, as observed in Figure 15, the computational results are in good overall agreement with the experimental results. The maximum discrepancy is at , where the largest deviation from experimental inlet waveform from Womersley inlet waveform exists (Figure 14(b)).

(a)

(b)

(c)

(d)

Figure 15: Pulsatile flow through a bend. The calculated streamwise velocity profiles at the plane of symmetry are compared with experimental results (circles) [31] at different time instants. ( and are the outer and inner bend radius, resp.), is the streamwise velocity at the symmetry plane, and is the maximum streamwise velocity at the symmetry plane at and .

Competing Interests

The authors declare that there is no conflict of interests regarding the publication of this article.

Acknowledgments

This work was supported by National Institute of Health (NIH) Grant R03EB014860 and the Center of Computational Research (CCR) of University at Buffalo. The authors are grateful to Professor Hui Meng for helpful discussions and providing them with the patient-specific geometry of an intracranial aneurysm.

References

F. H. H. Linn, G. J. E. Rinkel, A. Algra, and J. Van Gijn, “Incidence of subarachnoid hemorrhage: role of region, year, and rate of computed tomography: a meta-analysis,” Stroke, vol. 27, no. 4, pp. 625–629, 1996.

J. Xiang, S. K. Natarajan, M. Tremmel et al., “Hemodynamic-morphologic discriminants for intracranial aneurysm rupture,” Stroke, vol. 42, no. 1, pp. 144–152, 2011.

J. R. Cebral, F. Mut, J. Weir, and C. Putman, “Quantitative characterization of the hemodynamic environment in ruptured and unruptured brain aneurysms,” American Journal of Neuroradiology, vol. 32, no. 1, pp. 145–151, 2011.

G. Lu, L. Huang, X. L. Zhang et al., “Influence of hemodynamic factors on rupture of intracranial aneurysms: patient-specific 3D mirror aneurysms model computational fluid dynamics simulation,” American Journal of Neuroradiology, vol. 32, no. 7, pp. 1255–1261, 2011.

J. R. Cebral, F. Mut, J. Weir, and C. M. Putman, “Association of hemodynamic characteristics and cerebral aneurysm rupture,” American Journal of Neuroradiology, vol. 32, no. 2, pp. 264–270, 2011.

D. M. Sforza, C. M. Putman, and J. R. Cebral, “Hemodynamics of cerebral aneurysms,” Annual Review of Fluid Mechanics, vol. 41, pp. 91–107, 2009.

M. Tremmel, S. Dhar, E. I. Levy, J. Mocco, and H. Meng, “Influence of intracranial aneurysmto-parent vessel size ratio on hemodynamics and implication for rupture: results from a virtual experimental study,” Neurosurgery, vol. 64, no. 4, p. 622, 2009.

M. D. Ford, H. N. Nikolov, J. S. Milner et al., “PIV-measured versus CFD-predicted flow dynamics in anatomically realistic cerebral aneurysm models,” Journal of Biomechanical Engineering, vol. 130, no. 2, Article ID 021015, 2008.

S. Tateshima, Y. Murayama, J. P. Villablanca et al., “In vitro measurement of fluid-induced wall shear stress in unruptured cerebral aneurysms harboring blebs,” Stroke, vol. 34, no. 1, pp. 187–192, 2003.

P. Bouillot, O. Brina, R. Ouared, K. O. Lovblad, V. Mendes Pereira, and M. Farhat, “Multi-time-lag PIV analysis of steady and pulsatile flows in a sidewall aneurysm,” Experiments in Fluids, vol. 55, no. 6, article 1746, pp. 1–11, 2014.

T. B. Le, D. R. Troolin, D. Amatya, E. K. Longmire, and F. Sotiropoulos, “Vortex phenomena in sidewall aneurysm hemodynamics: experiment and numerical simulation,” Annals of Biomedical Engineering, vol. 41, no. 10, pp. 2157–2170, 2013.

T. A. Hope, M. D. Hope, D. D. Purcell et al., “Evaluation of intracranial stenoses and aneurysms with accelerated 4D flow,” Magnetic Resonance Imaging, vol. 28, no. 1, pp. 41–46, 2010.

G. Byrne, F. Mut, and J. Cebral, “Quantifying the large-scale hemodynamics of intracranial aneurysms,” American Journal of Neuroradiology, vol. 35, no. 2, pp. 333–338, 2014.

C. Poelma, P. N. Watton, and Y. Ventikos, “Transitional flow in aneurysms and the computation of haemodynamic parameters,” Journal of The Royal Society Interface, vol. 12, no. 105, 2015.

R. T. Jones, “Blood flow,” Annual Review of Fluid Mechanics, vol. 1, no. 1, pp. 223–244, 1969.

J. Jiang and C. Strother, “Computational fluid dynamics simulations of intracranial aneurysms at varying heart rates: a ‘patient-specific’ study,” Journal of Biomechanical Engineering, vol. 131, no. 9, Article ID 091001, 2009.

T. B. Le, I. Borazjani, and F. Sotiropoulos, “Pulsatile flow effects on the hemodynamics of intracranial aneurysms,” Journal of Biomechanical Engineering, vol. 132, no. 11, Article ID 111009, 2010.

S. S. Gopalakrishnan, B. Pier, and A. Biesheuvel, “Dynamics of pulsatile flow through model abdominal aortic aneurysms,” Journal of Fluid Mechanics, vol. 758, pp. 150–179, 2014.

J. R. Cebral, M. A. Castro, S. Appanaboyina, C. M. Putman, D. Millan, and A. F. Frangi, “Efficient pipeline for image-based patient-specific analysis of cerebral aneurysm hemodynamics: technique and sensitivity,” IEEE Transactions on Medical Imaging, vol. 24, no. 4, pp. 457–467, 2005.

I. Borazjani, L. Ge, T. Le, and F. Sotiropoulos, “A parallel overset-curvilinear-immersed boundary framework for simulating complex 3D incompressible flows,” Computers & Fluids, vol. 77, pp. 76–96, 2013.

A. Gilmanov and F. Sotiropoulos, “A hybrid cartesian/immersed boundary method for simulating flows with 3D, geometrically complex, moving bodies,” Journal of Computational Physics, vol. 207, no. 2, pp. 457–492, 2005.

I. Borazjani, L. Ge, and F. Sotiropoulos, “Curvilinear immersed boundary method for simulating fluid structure interaction with complex 3D rigid bodies,” Journal of Computational Physics, vol. 227, no. 16, pp. 7587–7620, 2008.

I. Borazjani and F. Sotiropoulos, “Numerical investigation of the hydrodynamics of carangiform swimming in the transitional and inertial flow regimes,” Journal of Experimental Biology, vol. 211, no. 10, pp. 1541–1558, 2008.

I. Borazjani and F. Sotiropoulos, “Vortex-induced vibrations of two cylinders in tandem arrangement in the proximity-wake interference region,” Journal of Fluid Mechanics, vol. 621, pp. 321–364, 2009.

I. Borazjani, “Fluid–structure interaction, immersed boundary-finite element method simulations of bio-prosthetic heart valves,” Computer Methods in Applied Mechanics and Engineering, vol. 257, pp. 103–116, 2013.

H. Asgharzadeh, I. Borazjani, J. Xiang, and H. Meng, “Vortex generation in two intracranial aneurysms,” in Proceedings of the ASME International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, V003T12A010, American Society of Mechanical Engineers, Buffalo, NY, USA, 2014.

Z. Song and I. Borazjani, “The role of shape and heart rate on the performance of the left ventricle,” Journal of Biomechanical Engineering, vol. 137, no. 11, Article ID 114501, 2015.

I. Borazjani, “A review of fluid-structure interaction simulations of prosthetic heart valves,” Journal of Long-Term Effects of Medical Implants, vol. 25, no. 1-2, pp. 75–93, 2015.

R. G. Bottom II, I. Borazjani, E. L. Blevins, and G. V. Lauder, “Hydrodynamics of swimming in stingrays: numerical simulations and the role of the leading-edge vortex,” Journal of Fluid Mechanics, vol. 788, pp. 407–443, 2016.

M. Daghooghi and I. Borazjani, “The influence of inertia on the rheology of a periodic suspension of neutrally buoyant rigid ellipsoids,” Journal of Fluid Mechanics, vol. 781, pp. 506–549, 2015.

C. C. M. Rindt, A. A. Van Steenhoven, J. D. Janssen, and G. Vossers, “Unsteady entrance flow in a 90° curved tube,” Journal of Fluid Mechanics, vol. 226, pp. 445–474, 1991.

K. Perktold and G. Rappitsch, “Computer simulation of local blood flow and vessel mechanics in a compliant carotid artery bifurcation model,” Journal of Biomechanics, vol. 28, no. 7, pp. 845–856, 1995.

D. N. Ku, “Blood flow in arteries,” Annual Review of Fluid Mechanics, vol. 29, no. 1, pp. 399–434, 1997.

N. M. Naughton, B. D. Plourde, J. R. Stark, S. Hodis, and J. P. Abraham, “Impacts of waveforms on the fluid flow, wall shear stress, and flow distribution in cerebral aneurysms and the development of a universal reduced pressure,” Journal of Biomedical Science and Engineering, vol. 7, no. 1, pp. 7–14, 2014.

J. C. R. Hunt, A. Wray, and P. Moin, “Eddies, stream, and convergence zones in turbulent flows,” Center for Turbulence Research Report CTR-S88, 1988.

S. Oka and M. Nakai, “Optimality principle in vascular bifurcation,” Biorheology, vol. 24, no. 6, pp. 737–751, 1987.

K. Y. Lee, Y. H. Sohn, J. S. Baik, G. W. Kim, and J.-S. Kim, “Arterial pulsatility as an index of cerebral microangiopathy in diabetes,” Stroke, vol. 31, no. 5, pp. 1111–1115, 2000.

Z. Kulcsár, L. Augsburger, P. Reymond et al., “Flow diversion treatment: intra-aneurismal blood flow velocity and WSS reduction are parameters to predict aneurysm thrombosis,” Acta Neurochirurgica, vol. 154, no. 10, pp. 1827–1834, 2012.

D. Vigolo, S. Radl, and H. A. Stone, “Unexpected trapping of particles at a T junction,” Proceedings of the National Academy of Sciences of the United States of America, vol. 111, no. 13, pp. 4770–4775, 2014.

K. K. Chen, C. W. Rowley, and H. A. Stone, “Vortex dynamics in a pipe T-junction: recirculation and sensitivity,” Physics of Fluids, vol. 27, no. 3, Article ID 034107, 2015.

I. S. Sullivan, J. J. Niemela, R. E. Hershberger, D. Bolster, and R. J. Donnelly, “Dynamics of thin vortex rings,” Journal of Fluid Mechanics, vol. 609, pp. 319–347, 2008.

H. Meng, V. M. Tutino, J. Xiang, and A. Siddiqui, “High WSS or low WSS? Complex interactions of hemodynamics with intracranial aneurysm initiation, growth, and rupture: toward a unifying hypothesis,” American Journal of Neuroradiology, vol. 35, no. 7, pp. 1254–1262, 2014.

J. R. Cebral and H. Meng, “Counterpoint: realizing the clinical utility of computational fluid dynamics—closing the gap,” American Journal of Neuroradiology, vol. 33, no. 3, pp. 396–398, 2012.

L. Ge and F. Sotiropoulos, “A numerical method for solving the 3D unsteady incompressible Navier-Stokes equations in curvilinear domains with complex immersed boundaries,” Journal of Computational Physics, vol. 225, no. 2, pp. 1782–1809, 2007.