Numerical analysis of aneurysm using pulsatile blood flow through a locally expanded vessel

NUMERICAL ANALYSIS OF ANEURYSM USING PULSATILE BLOOD

FLOW THROUGH A LOCALLY EXPANDED VESSEL

H. Niazmand, A. Sepehr 1 and P. B. Shahabi

Mechanical Engineering Dept., Ferdowsi University of Mashhad P.O. Box No. 91775-1111, Mashhad, Iran

1

e-mail: al_se61@stu-mail.um.ac.ir

Key words: Aneurysm, Pulsatile Flow, Numerical Analysis, non-Newtonian, Waveform Abstract. A numerical analysis based on a finite volume approach is employed for a 2-D axisymmetric, incompressible, laminar and transient flow in order to simulate the pulsatile blood flow in aneurysm. The unsteady nature of the blood is studied for three Womersley numbers and for three rheological models including Newtonian, Power law and Quemada. A physiological flow input waveform is presented in terms of Fourier series obtained by curve fitting to the experimental data. Introducing the hemodynamic wall parameters including wall shear stress (WSS), wall shear stress gradient (WSSG) and oscillatory shear index (OSI), results are presented for the variation of these parameters for all rheological models at different Womersley numbers. In addition, the effects of the rheological models and the Womersley number on the formation, growth, and shedding of the vortexes in the aneurysm are studied. The results implicates that the Womersley number has a significant effect on the flow field, while the effect of rheological model on the flow pattern is negligible. Besides, the reattachment and separation points are recognized based on the OSI and WSS diagrams.

1 INTRODUCTION

The precise study of the pulsatile blood flow through the vascular system is complex. The difficulty pertains to the complicated nature of blood as a suspension of red and white blood cells, platelets, and lipid globules dispersed in a colloidal solution of proteins. So far, many experimental relations have been presented for non-Newtonian behavior of blood based on the apparent viscosity in literature1. In the present study, three rheological models including Newtonian, Power Law and Quemada are investigated.

Fukushima et al.2 studied the steady and unsteady aneurysm flow. They showed that even though the steady flow models provide the researchers with some understanding on the flow through arteries, more practical models need to take into account the pulsatile nature of the blood. Consequently, finding an appropriate waveform as a flow input pulse is essential for the simulation of blood flow through a vessel. Sinusoidal function for flow input waveform

1 _{Corresponding author}

has been commonly used in previous studies3,4,5,6,7,8,9. Tu and Deville3 used a flow input waveform as a function of axial velocity in their numerical study based on the experimental data provided by Siouffi et al.5,6. Long et al.7 employed a time-varying volume flow rate function as the flow input waveform. In the numerical study of Buchanan et al.4 a transient Reynolds number was considered for the waveform. All above mentioned waveforms were applied to the pulsatile blood flow through a stenosed vessel. In the present study, an input waveform is developed based on the curve fitting to the experimental data.

Since the vascular diseases such as arteriosclerosis and aneurysm make interesting fluid dynamics and hemdynamics effects, the researchers are stimulated to study the blood flow through an infected vessel. However, most of the available studies are focused on the blood flow through the stenosed vessels and there is limited information about the blood flow in an expanded vessel. Ishikawa et al.9 studied the vortex enhancement in locally stenosed tubes as well as expanded tubes. Present work is concentrated on the pulsatile blood flow through a locally expanded vessel called Aneurysm considering the hemodynamics parameters variations for different rheological models at various Womersley numbers. An aneurysm is a local abnormal dilation of a blood vessel wall that presents a peril to human health from rupture, clotting or degeneration of the arterial wall. The combinations of the biological factors and biophysical forces have increased the complexity of the aneurysm development, growth and rupture10.

The localized aneurysm formation strongly implicates the relationship between the hemodynamic forces and arterial geometry in the disease process11. Several authors have performed experimental and statistical analysis on aneurysm size, or aneurysm aspect ratio to evaluate the risk of aneurysm rupture12-17. The statistical analysis depends on the population and cannot be freely applied to the group outside the study population. In an attempt to understand the development, growth and rupture of aneurysm, aneurysm located on straight vessel model is frequently used18_{. This study is focused on three-dimensional axisymmetric} sinusoidal aneurysm geometry in a straight vessel with particular attention to the non-Newtonian behavior and pulsatile nature of blood flow. The sinusoidal geometry of aneurysm is close enough to the clinical observations. Utilizing the numerical analysis, the hemodynamic of aneurysm, with different Womersley numbers, for different rheological models are investigated. The results can provide some insight for the study of the vascular disease, and with a medical investigation can aid the clinicians to decide the suitable treatment for individual patients.

2 GEOMETRICAL PARAMETERS OF ANEURYSM

The geometrical parameters are usually used to determine the category and characteristic of aneurysms for making treatment decisions. Without the clear understanding of how these geometrical parameters influence the onset and the treatment of aneurysm, the precise analysis of the problem is not feasible. These parameters include aneurysm size, aspect ratio and vessel curvature.

should be observed rather than treated19. However, a lot of clinical reports revealed that aneurysm can rupture at smaller than 5mm or remain non-ruptured at greater than 10mm. In the present study, we consider a local expansion for aneurysm geometry, which is corresponding to a 75% of expansion in a blood vessel.

The schematic of local expansion inside a blood vessel is shown in Fig 1. In this study, the following sinusoidal relation is employed for a tubular geometry by radius r0 and the local, smooth, and axisymmetric expansion. This choice is considered because the sinusoidal profile well approximates the real cases.

⎪ ⎩ ⎪ ⎨ ⎧ > ≤ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + + = = 0 0 0 0 z z if 1 z z if z z cos 1 2 1 r ) z ( r ) z ( R π δ (2.1)

where, r0 is the vessel radius, z0 is the aneurysm half-length and δ is the thickness of expanded vessel as shown in Fig 1.

-4

-2

0

2

4

z coordinate

0

1

2

r

c

oor

dina

te

Fig 1. Geometrical schematic of aneurysm

The characteristic of aneurysm neck and size defines the volume and affect the velocity of flow into and out of the aneurysm, thereby influences the distribution of hemodynamic stresses12. Thus, Ujiie et al.12,13 defined the geometrical parameter, aspect ratio as the ratio of aneurysm depth to aneurysm neck width. They showed that 80% of all aneurysms with aspect ratio of 1.6 or higher were ruptured.

3 FLOW PARAMETERS OF ANEURYSMS

The physiology of aneurysm is closely related to the local hemodynamic stresses. Besides the geometrical factors, the flow parameters are also associated with the magnitude and distribution of hemodynamic stresses. Unsteady nature, nonlinear arterial wall property and non-Newtonian effects are important flow parameters.

Flow instabilities are developed with pulsatile flow because of more important influences of the inertia of the fluid. With unsteady flow in arterial curvature and bifurcation sites, there

δ

2

are regions in which the flow is reversed during part of the cardiac cycle. Hence the wall shear stress (WSS) varies from a large magnitude in one direction to negative values during part of the cardiac cycle4. Therefore, to clarify the physiological hemodynamic condition and to precisely predict the outcome of treatment, pulsatile flow simulation is adopted in this study. In the present study, curve fitted formula in terms of a Fourier series function is obtained based on the experimental data used by Tu and Deville3 as shown in Fig 2. A time dependent spatial average velocity profile according to the equation 3.1 is specified as the flow input waveform.

∑ _⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + = = n i i mean _T t b T t 1 aicos 2 sin 2 U U(t) π π (3.1)

The transient Reynolds number Re(t)=ρU(t)D/µ∞ is based on the transient input velocity and the vessel diameter, where µ∞ is the limiting high shear rate Newtonian viscosity. The mean Reynolds numbers is defined based on the time-averaged mean velocity Remean=ρUmeanD/µ∞. 0 0.5 1 1.5 2 t/T -5 0 5 10 15 20 25 30 v( c m /s ) Souffi et al., 1984. 5 fourier series expressions 7 fourier series expressions 10 fourier series expressions 15 fourier series expressions

0.86 0.88

6.0 6.2 6.4 6.6

Fig 2. the flow input waveform for a given period based on Fourier series fitted to experimental data

and time consuming. Thus, rigid wall model is assumed in this study.

Blood behaves as a non-Newtonian fluid, i.e. the viscosity decreases with increasing shear rate. This non-Newtonian effect becomes more significant in small diameter vessels (d<0.1mm) and low shear rate (τ< 100 s-1)19. In the large vessels, the viscosity of blood is approximately 3.5 centi Poise19.

Two characteristic features determine the certain conditions for the blood as a non-Newtonian fluid (i) yield stress and (ii) apparent viscosity as a strong function of shear rate. In order to take these effects into account, several equations have been presented in literature for blood behavior1,3,5,6_{. In this study, three fluid rheological models, including Newtonian,} Power law, and Quemada are considered. For Newtonian behavior, a linear stress-strain relationship exhibits a constant viscosity for blood. For non-Newtonian models, a relation is used in the form of apparent viscosity.

= =η∞ η constant Newtonian (3.2) 1 n m − = γ

η & Power law (3.3)

where m and n are constant coefficients depending on the type of fluid. Since blood is known as a shear-thinning fluid, thus n<1.

For the Quemada model the apparent viscosity is given with the following correlation4_. 2 0 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + = _∞ γ λ τ η η & Quemada (3.4)

The various coefficients are listed in the Table 1, and the variation of viscosity with respect to shear rate is demonstrated in Fig 3. The Quemada model maximum shear rate,γ& , is max 256.7 s-1, while for the Power law model γ&max =76.8 s

Table 1. Thermophysical and geometrical information

Rheological information for three models in cgs system η=µ∞=0.0309 g/cm.s m=0.1260 g/cm.sn n=0.8004 η∞=0.02654 g/cm.s τ0=0.04360 g/cm.s2 λ=0.02181 s-1 Newtonian model

Power law model Quemada model Geometrical information 5 . 0 =

δ for local expansion of 75% Half length of aneurysm Z0 ( normalized ) =2

10

10

10

10

10

10

10

γ

10

10

η

Fig 3. Viscosity variation vs. shear rate compared with experimental data

4 GOVERNING EQUATIONS

0 v . = ∇ r (4.1)

( )

Since the apparent viscosity depends on the local shear rate, the momentum equation is derived with taking this effect into account. In this equation, τt represents the stress tensor that is expressed as:

( )

[

]

where tr indicates the transpose. The relation of stress tensor is varied corresponding to the rheological model.

The boundary conditions are no slip at the wall, stress free at the outlet and symmetry at the centerline. A time dependent parabolic velocity profile based on the time variation of the Re(t) is specified at the inlet. The input velocity profile is determined numerically for an adequately long straight tube with the same diameter and a mass flow rate associated with the Re=1 and related rheological properties.

5 NUMERICAL MODEL

The 2-D axisymmetric time-dependent Navier-Stokes equations are solved numerically by a finite volume algorithm in a body fitted coordinate system. The SIMPLE algorithm provides the linkage between the velocity and pressure fields. Staggered grid is used to avoid pressure field checker board effect. All calculations are carried out on a mesh size of 300×30 with an expansion factor of 0.65 in both directions to cluster the mesh grid points near the wall and in the expansion region. Computations are performed with different mesh grid points of 350×40, 300×25, 250×40 and it is revealed that the mesh grid points of 300×30 satisfy the mesh-independency of the results.

During the iterative sequence, the residuals of u, v and mass conservation are calculated and the convergence criteria are set to 10-4. In order to eliminate the start-up disturbances of the pulsatile flow, calculations are continued for a minimum 4 periods. It was found that for the rheological and physical data considered here the maximum difference in the results between the fourth and the fifth period was less than 0.5%. The number of time steps varies for different Womersley numbers. For Wo=4.0, greater time step is required due to the larger time span of a period. The number of time steps used in the numerical calculations are 5200, 1600 and 960 for Wo=4, 7.5, and 12.5, respectively.

6 RESULTS AND DISCUSSION

The key parameters influence the flow pattern in the expanded region are the Womersley number, the type of rheological model, the flow input pulse, and the degree of the local expansion.

Results include the details of the flow field for each rheological models and Womersley numbers, as well as hemodynamic wall parameters. The mean Reynolds number based on the mean input velocity is 614.95. All thermophysical properties and geometrical parameters are given in Table 1.

The velocity fields are presented for three rheological models at different Womersley numbers (Figs 5-13). The sequence of the flow field streamlines corresponds to five time levels in one cycle as shown in Fig 4. Time levels are chosen at points to clearly show the formation, growth, expansion, and shedding of the vortexes in the local expansion of the vessel.

A sequence of transient flow patterns are plotted for a Newtonian fluid at Wo=7.5 in Fig 8. The effects of Womersley number and rheological models on the flow pattern will be discussed later. At time level T1, the Reynolds number of the fluid increases and thus the flow accelerates. It is expected that at high enough Reynolds number a recirculating zone forms inside the expanded area. Present calculation shows that a small vortex forms in the expansion before the time level T2. This vortex grows in time and at time level T2, it occupies almost all the expansion region. As Fig 4 shows, flow deceleration starts after the time level T2 and continues to time level T4, therefore, time level T3 in the middle of the deceleration process results in the growth of the vortexes in the expanded zone. T4 time level corresponds to the start of the next acceleration phase. At this time level the vortex moves to the core region of the flow due to the presence of a strong reverse flow near the wall. At the end of this acceleration phase, flow near the wall changes direction again, which leads to the reduction of the vortex strength and causes the vortex to move toward the expansion and disappear later leading to the generation of a free-vortex flow pattern in the expansion region. It is worth mentioning that in the deceleration phase, time level T5, corresponding to the time interval 0.85< t/T<1, a small vortex forms inside the expansion and disappears shortly due to the major reduction in the flow Reynolds number.

t/T V( c m /s ) 0 0.5 1 1.5 -5 0 5 10 15 20 25 30

t

t

t

t

t

Fig 4. Time levels chosen during a cycle

Effect of Womersley number on the flow pattern can be seen in Figs 5,11, where the same flow fields are presented for Wo=4.0 and 12.5. Considering the definition of the Womersley

number ( 1/2

0(2 / T) r

Wo= π υ ), lower Womersley numbers are associated with larger period and vice versa. This means that the flow patterns are developing in larger time span for smaller Womersley numbers. Therefore, there is more time for the formation and growth phases (t2,t3) of the vortexes for Wo=4.0 as compared to Wo=7.5 and can be seen in Fig 5 compared to Fig 8. However, during the deceleration phase, which corresponds to the weakening of the vortex, at Wo=4.0 the vortex is smaller as compared to Wo=7.5 Since there is more time for the reverse flow to reduce the vortex strength. The change in the flow patterns at Wo=12.5 as compared to Wo=7.5 can be explained in the same line.

Effects of the rheological models on the flow pattern are shown in Fig 8,9 and 10, where, streamlines at Wo=7.5 are presented for the power law and Quemada models. Clearly, different rheological models have slight affects on the flow patterns. The same feature is reported for the flow in stenosed vessels23.

Hemodynamics can be defined as the study of blood flow and the physical forces concerned in the circulatory system. Although patients’ genetic, anatomic, histological, and physiological factors have been implicated in the pathogenesis of aneurysm, the hemodynamic forces play significant roles in the development and progression of aneurysm.

dt w U T 1 WSS T 0 2

∫

For indicating the changes due to the shearing forces, and hence ‘aggravating effects’ on the endothelium of arterial wall, the wall shear stress gradient is defined as4_:

∫

To describe the unsteady nature of blood flow through vascular system, the oscillatory shear index is expressed as4.

⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − =

∫

∫

OSI is an effective parameter for identifying the points of time-averaged separation and reattachment and has a value between 0 and 0.5.

In above Equations T, ρ, U, S and D are the period, density, area-averaged velocity of mean Reynolds number, local coordinate along the wall and tube diameter respectively.

The response to an increase wall shear stress is to dilate and then remodel to a larger diameter with the same arterial structure. Conversely, decreased shear stress resulting from lower flow induces vasoconstriction and remodel to decrease in internal vessel radius19. The responses of arteries to the different hemodynamic environment may constitute to the normal adaptation or pathological disease19.

The hemodynamic wall parameters for three Womersley numbers, Wo=12.5, 7.5 and 4.0, are shown in Figs 14-16. Time-averaged separation and reattachment points can be estimated from both WSS and OSI graphs. The first axial location where WSS becomes zero indicates the time-averaged separation point, while the second one corresponds to the reattachment point. These points are also associated with the first and second maxima of the OSI curves, respectively.

Time-averaged reattachment point shows the interaction of viscous and inertial forces and it differs upon the rheological model. Reattachment point for all models slightly decreases with increasing Womersley number. The time-averaged reattachment point is maximum for the Newtonian model and minimum for the Power law model at each Womersley number. The Quemada model is located between the Newtonian and Power law values.

where the velocity profiles are reserved at the wall. The reverse flow is associated with the atheroma development, which in turn is related to the degeneration of the arteries wall.

Comparing axial time averaged variation of OSI in Figs 14-16 indicates that in expansion zone all three cases reaches their maximum values of about 0.5, while outside this region there is a significant difference on OSI values as the Womersley number increase. Apparently, increasing in Womersley number increases the OSI values outside the expansion region. Furthermore, one can find more similar behavior for Newtonian and Quemada models, while Power law behaves differently.

The wall shear stress gradient shows the strength of the vortex, which has been generated in the flow pattern. The strength of the vortex that forms in the expansion region increases with the decrease of the Womersley number (Figs 14-16). Thus, fluid flow with the lower Womersley numbers generates the stronger vortexes. The vortex strength at each Womersley number has no considerable difference between the three different rheological models. This is in contrast to the blood flow in the stenosed vessel23.

7 CONCLUSION

Three rheological models are used in the study of the two dimensional pulsatile blood flows through a locally expanded vessel. A non-sinusoidal periodic function is considered for the flow input waveform based on the experimental data. The effects of Womersley number variations as a key parameter in pulsatile flows are investigated for three rheological models through the wall hemodynamic parameters (WSS, OSI and WSSG). The time-averaged separation and reattachment points are clearly indicated by the OSI and WSS parameters. Furthermore, it is found that the flow pattern around the expansion is strongly affected by the Womersley number, while it is less sensitive to the rheological models.

REFRENCES

[1] Zhang Jian-Bao, Kuang Zhen-Bang, “Study on Blood Constitutive Parameters in Different Blood Constitutive Equations,” Journal of Biomechanics, 33:355-360,2000 [2] Fukushima T., Matsuzawa T., Homma T., “Visualization and finite element analysis of

pulsatile flow in models of the abdominal aortic aneurysm,” Biorheology, 26:109-130, 1989

[3] Tu, C., Deville, M., “Pulsatile Flow of Non-Newtonian Fluids through Arterial Stenoses,” J. Biomechanics, Vol. 29, No. 7, 899-908, 1996

[4] Buchanan, J.R., Kleinstreuer, C., Comer, J.K., “Rheological Effects on Pulsatile Hemodynamics in a Stenosed Tube,” Computers & Fluids, Vol. 29, 695 -724, 2000

[5] Siouffi, M., Deplano, V., Pelissier, R., “Experimental Analysis of Unsteady Flows through a Stenosis,” Journal of Biomechanics, Vol. 31, 11-19, 1998

[7] Long, Q., Xu, X.Y., Ramnarine, K.V., Hoskins, P., “Numerical Investigation of Physiologically Realistic Pulsatile Flow trough Arterial Stenosis,” Journal of Biomechanics, Vol. 34, 1229-1242, 2001

[8] Ishikawa, T., Guimaraes, L.F.R., Oshima, S., Yamane, R., “Effect of Non-Newtonian Property of Blood on Flow through a Stenosed Tube,” Fluid Dynamics Research, Vol. 22, 251-264, 1998

[9] Ishikawa, T., Oshima, S., Yamane, R., “Vortex Enhancement in Blood Flow through Stenosed and Locally Expanded Tubes,” Fluid Dynamics Research, Vol. 26, 35-52, 2000. [10] Kayembe KNT, Sasaharam M, Hazama F. “Cerebral Aneurysms and Variations of Circle

of Willis,” Stroke, 15:846-850, 1984

[11] Burleson A.C., Turitto V.T., “Identification of quantifiable hemodynamic factors in the assessment of cerebral aneurysm behavior,” Thromb Haemost, 76:118-123, 1996

[12] Ujiie H., Tachibana H., Hiramatsu O., et al.,”Effects of size and shape (aspect ratio) on the hemodynamics of saccular aneurysm: a possible index for surgical treatment of intracranial aneurysms,” Neurosurgery, 45:119-130, 1999

[13] Ujiie H., Tamano Y., Sasaki K., Hori T., “Is the aspect ratio a reliable index for predicting the rupture of a saccular aneurysm?,” Neurosurgery, 48(3):495-503, 2001 [14] Kassell, N.F., Torner, J.C., “Size of Intracranial Aneurysms,” Neurosurgery,

12(3):291-297, 1983

[15] Crompton M.R., “Mechanism of growth and rupture in cerebral berry aneurysms,” Br Med J, 5496:1138-42, 1966

[16] Kumar B.V., Naidu K.B., “Hemodynamics in aneurysm, ”Computers and Biomedical Research,” 29:119-139. 1996

[17] Liou T.M., Chang W.C., Liao C.C., “Experimental study of steady and pulsatile flows in cerebral aneurysm model of various sizes at branching site,” J Biomech Eng , 119(3):325-32, 1997

[18] Aenis M., Stancampiano A.P., Wakhloo A.K., et al., “Modeling of Flow in a Straight Stented and Nonstented Side Wall Aneurysm Model,” J. Biomech. Eng , 119:206-213, 1997

[19] Hoi Yie Meng, “Correlation of Hemodynamic Forces and Aneurysm Geometry: Results of a Three-Dimensional Computational Fluid Dynamics,” MSc thesis, Dept. of Mechanical Aerospace Engineering, University of New York at Buffalo, US, 2003

[20] Shahabi, P.B., “Study of pulsatile blood flow in stenosed vessels,” MSc thesis, Dept. of Mechanical Engineering, Ferdowsi University of Mashhad, Iran, 2003

[21] Shahabi, P.B., Modarres Razavi, M.R., “Pulsatile blood flow simulation in stenosed vessels by non-Newtonian fluid,” Proceedings of 8th Fluid Mechanic Conference (in Persian), Tabriz University, Tabriz, Iran, September 2003

[22] Seyedein, S.H., “Simulation of Fluid Flow and Heat Transfer in Impingement Flows of Various Configurations,” MS Thesis, Dept. of chemical Engineering, McGill University, 1993

Fig 13. Flow patterns for Quemada model at Wo=12.5 Z OS I -5 0 5 10 0.1 0.2 0.3 0.4 0.5 Newtonian Power Law Quemada Z WS S -5 0 5 10 -0.02 0 0.02 0.04 0.06 0.08 Newtonian Power Law Quemada Z WSS G -5 0 5 10 0 0.04 0.08 0.12 Newtonian Power Law Quemada

Z OS I -5 0 5 10 0.1 0.2 0.3 0.4 0.5 Newtonian Power Law Quemada Z WS S -5 0 5 10 -0.02 0 0.02 0.04 0.06 0.08 Newtonian Power Law Quemada Z WSS G -5 0 5 10 0 0.04 0.08 0.12 Newtonian Power Law Quemada

Fig 15. Hemodynamic wall parameters for Wo = 7.5

Z OS I -5 0 5 10 0.1 0.2 0.3 0.4 0.5 Newtonian Power Law Quemada Z W SS -5 0 5 10 0 0.02 0.04 0.06 0.08 Newtonian Power Law Quemada Z W SSG -5 0 5 10 0 0.04 0.08 0.12 Newtonian Power Law Quemada