VIEWS: 12 PAGES: 5 POSTED ON: 3/27/2011
7th Annual Meeting of the CFD Society of Canada, Halifax NS, May 1999. NUMERICAL MODELLING OF THE DISEASED CAROTID ARTERY BIFURCATION USING A PHYSIOLOGICALLY RELEVANT GEOMETRIC MODEL Ronald Gin+ and Anthony G. Straatman++ Department of Mechanical & Materials Engineering The University of Western Ontario London, Ontario, Canada N6A 5B9 David A. Steinman++ The John P. Robarts Research Institute The Department of Medical Biophysics The University of Western Ontario London, Ontario, Canada N6A 5B9 INTRODUCTION The intention of the present work is to carry out a systematic grid-independence study on a model of Atherosclerosis is a disease that is the carotid artery bifurcation for the purpose of characterized by the hardening and thickening of the calculating steady and pulsatile flows. Grid- arterial walls due to the formation of plaque. As this independence is based on predictions of the wall shear disease progresses, the formation of plaque reduces the stress in the normal carotid artery. The normal carotid arterial passage area creating uncharacteristic blood artery is chosen because experimental evidence flow patterns. As a result, this restriction, if severe indicates that no turbulence exists in the flow field. enough, can cause individuals to suffer cardiac arrest or The grid-converged steady flow predictions will then stroke. Also, small pieces of plaque can become be validated by comparisons to Magnetic Resonance dislodged and travel downstream restricting blood flow Imaging (MRI) data. The information obtained from into smaller vessels, such as the arterioles or capillaries, this study will be utilized in a continuing study to again, causing cardiac arrest or stroke. The plaque simulate the diseased carotid artery models created by deposition is focal and usually occurs in sharp Smith et al [7]. curvatures and bifurcations in the cardiovascular system such as the carotid artery. GEOMETRIC MODELS Much research has been performed on the biological aspect of this phenomena and it has lead to A recent study by Smith et al [7] presented a the hypothesis that the flow disturbances will play a family of physiologically relevant geometries that major influence in the localization of the disease [1,2]. describe several stages of lumen narrowing of the Considerable research has been carried out to analyze carotid artery bifurcation. These geometries were the hemodynamic flow patterns in the carotid artery generated by averaging digitized bi-plane angiograms [3], however, a majority of these simulations have been of symptomatic patients with different magnitudes of performed on idealized models of the normal carotid the stenosed artery (see Figure 1). The first geometry, bifurcation. The results obtained from these studies model (a), describes the disease free carotid artery and have enhanced the understanding of the hemodynamic is very similar to the geometry published by Bharadvaj patterns and their association with the development of et al [3]. The main difference between the two atherosclerosis [4,5,6]. However, vital information can geometries is that Smith’s model incorporates a more be determined from analyzing diseased arteries, which physiological realistic description of the downstream may aid in diagnosing the risk factors of stroke. Also, angle between the internal carotid artery (ICA) and the it may assist in distinguishing the appropriate external carotid artery (ECA) branches. The ICA and procedure for the treatment of patients. the ECA branches were bent inwards to be parallel with the common carotid artery (CCA). This implies that + Graduate Student ++ Assistant Professor 7th Annual Meeting of the CFD Society of Canada, Halifax NS, May 1999. where U i is the velocity in the xi directions and ρ and µ are the density and dynamic viscosity, respectively. The commercial Computational Fluid Dynamics (CFD) package, CFX-Tfc [8], was used to impose the boundary conditions and solve the discrete equation set. CFX-Tfc utilizes an element–based, finite volume technique and operates with hybrid, unstructure grids. Internal parameters were adjusted such that the discretized equations were representative of a 3- dimensional, steady, incompressible and Newtonian fluid flow condition. Upon simulation, the most accurate, general second-order convention scheme was used to minimize the effects of false diffusion. The second-order convection scheme is based on upwinding along a streamline with the addition of Numerical Advection Correction (NAC). Modifications were made to parts of the CFX- Tfc source code so that profile conditions could be imposed at the inlet and outlet flow boundaries. With Figure 1. Wireframe drawings of the normal and stenosed this type of boundary condition implementation, carotid bifurcation generated by Smith et al. [7] (a) Modified parabolic boundary profiles could be specified for Bharadvaj normal geometry; concentric geometries: (b) 30%, (c) simulating steady computations. Also, for the 50%, (d) 60% ,(e) 70%, and (f) 80% stenosis, and eccentric geometries: (g) 30%, (h) 50%, (i) 60%, (j) 70%, and (k) 80% continued study, time periodic waveforms of a stenosis. characteristic blood pulse can now be easily imposed. 1200 after the CCA divides, the two branches will proceed vertically at a physiologically realistic inclination, 1000 instead of as continuously diverging branches, as modelled by Bharadvaj et al [3]. The remaining 800 geometries published by Smith et al [7] characterize the carotid artery at varying stages of lumen narrowing, Re ranging from 30% to 80% stenosis, with either 600 symmetric or eccentric stenoses. The symmetric stenosis represents the growth of the disease such that 400 the arterial lumen is positioned in the middle of the artery. The eccentric stenosis describes the case where 200 the disease is focused to one side of the artery, more commonly on the outer wall, within the sinus. 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t/T NUMERICAL SOLVER AND BOUNDARY CONDITIONS Figure 2. Average Reynolds number of the time periodic blood pulse through the carotid artery. To resolve the laminar flow in the carotid artery bifurcation models, the conservation of mass and The geometry that was used in the grid- conservation of momentum equations must be solved. independence study was the normal carotid artery, These equations are given in strongly conservative model (a) in Figure 1. This geometric model was form as: chosen for two reasons. First, on the basics of experimental evidence, it is certain that the normal ∂ρ ∂ (ρU j ) carotid flow is not turbulent and second, MRI data + =0 exists for this case so that the grid converged results (of ∂t ∂x j velocity) can validated. The grid will also be ∂ ( ρU i ) ∂ (ρU j U i ) ∂P ∂ ∂U i ∂U j µ applicable for mildly stenosed cases where no + =− + + ∂t ∂x j ∂x i ∂x j ∂x j ∂x i turbulence exists on the basics that the grid will become condensed in the vicinity of the constriction. 7th Annual Meeting of the CFD Society of Canada, Halifax NS, May 1999. To ensure that the resulting grid will be satisfactory for et al. [9], showed that wall shear stress independence is both steady and pulsatile simulations, the grid- a more accurate convergence criteria (but much more independence is to be based on steady flow calculations difficult to obtain). at the average Reynolds number observed in a characteristic, pulsatile blood flow (see Figure 2). RESULTS Simulations will also be computed at the peak Reynolds number. On the basis of both clinical and Simulations were computed for Reynolds numerical studies, the expected average and peak numbers of Re=250 and 1100 with varying mesh Reynolds number for pulsatile flow in the normal densities, as described in Table 1. The densities range carotid artery is approximately Re=250 and 1100, from approximately 50K to 800K, first order respectively. The peak Reynolds number is observed tetrahedral elements. Comparisons were made between during the average peak flow in the systolic part of the the converged results of each set of subsequent grids blood pulse. (one twice as dense as the other), however, particular In simulating the steady computations, a attention is directed at the 200K, 400K and 800K parabolic profile was imposed on the inlet of the CCA. meshes. Velocity and wall shear stress grid The ECA was set to have a parabolic profile for independence were compared by computed errors. Re=250 and a mass specified boundary condition for From these errors the suggested element sizes are Re=1100. The mass flow split between the ECA and reported for specified areas within the domain of the ICA was set to have 44% of the inlet mass flow leaving carotid artery. Also, comparisons at Re=250 are made the ECA. The ICA boundary condition was with experimental MRI data. implemented using a pressure specified opening which would act as a pressure reference within the carotid Number of elements Number of Nodes Element Size/D artery domain. The Reynolds number was set based on 50K 50757 11458 0.26 the inlet radius of 1 cm. To implement Re=250, the 100K 99746 20986 0.20 kinematic viscosity was chosen to be 8.0*10-3 cm2/s, 200K 400K 199872 399789 40158 77583 0.16 0.12 therefore, producing an average inlet velocity of 1 800K 802402 152579 0.10 cm/s. For simulations at Re=1100, the equivalent average inlet velocity was employed with the Table 1: Grid density information appropriate kinematic viscosity. The initial grid was designed to contain For the Re=250, the results of the velocity approximately 50,000 tetrahedral elements and profiles and the wall shear stress contours can be seen subsequent grids were obtained by sequentially in Figure 3. The velocity profiles, Figure 3a, are doubling the number of elements until the desired plotted in several positions along the symmetry plane accuracy was obtained. In all the models, elements and display the profiles for the grid densities containing were generated based on a specified element size; the 200K, 400K and 800K elements. The profiles are height of a symmetric or perfect tetrahedral. Due to the plotted with values that are scaled by 0.5 in relation to difficulty in generating smooth transitions from course the geometry. The relative errors between each to fine elements, in the grid generator, all elements subsequent grid density, for both velocity and wall were generated based on the same element size. This shear stress, are seen in Table 2. The table reports increased the grid density in areas where fewer errors for five distinct regions; the CCA, the bifurcation elements would have been sufficient, such as areas on ICA side, the bifurcation on ECA side, the ICA and downstream in the CCA and perhaps, downstream in the ECA. The errors are computed for each specified the two branches. However, the benefit of such region by comparing several points within the 3- simulations, will aid in determining the required dimensional geometry. The critical region is found to element size, or the number of elements, that would be in vicinity of the bifurcation where the largest accurately model the area in the vicinity of the 050k-100k 100k-200k 200k-400k 400k-800k bifurcation. The accuracy was assessed by comparing Velocity WSS Velocity WSS Velocity WSS Velocity WSS CCA Max. 3.80 8.43 2.14 4.25 1.02 8.28 0.59 6.07 predicted results for the wall shear stresses at several Ave. 2.14 2.16 0.74 1.61 0.30 1.94 0.21 1.30 points in the region of the bifurcation. The desired BIF. Max. 94.51 46.43 47.67 40.66 35.28 26.34 6.51 25.24 accuracy is obtained when the stresses from two ICA Ave. 9.90 11.84 7.81 10.91 4.82 7.91 1.44 6.85 BIF. Max. 16.38 41.37 10.21 29.42 9.96 23.14 5.23 20.46 subsequent grids (one twice as dense as the other) are ECA Ave. 3.51 13.38 1.94 10.31 2.27 8.07 1.48 5.02 predicted to within the prescribed accuracy. The wall ICA Max. 14.61 32.07 7.50 29.94 10.53 21.06 4.77 9.61 shear stress was chosen over the velocity as the grid Ave. 6.34 5.71 3.34 6.00 3.10 3.53 1.41 3.20 ECA Max. 13.81 23.31 13.05 15.91 9.83 9.22 4.27 3.78 independence parameter because it is thought that the Ave. 7.33 7.12 3.99 4.68 2.53 2.14 1.25 1.23 wall shear stress may have the greatest impact in the assessment of risk. Furthermore, recent work by Ethier Table 2: Velocity and wall shear stress convergence errors. 7th Annual Meeting of the CFD Society of Canada, Halifax NS, May 1999. 6.00 8.00 P O N M 10.00 6.00 8.00 L K I 8.00 J 10.00 6.00 4.00 8.00 G H 8.00 2.00 12.00 12.00 E 1.00 15.00 F 4.00 C 0.25 0.25 D 0.50 4.00 1.00 B Solid - 800K 2.00 Solid - 800k Dashed - 400K 4.00 Dotted - 400k Dotted - 200K 1.0 4.00 4.00 A Figure 3a: Velocity profiles Figure 3b: Wall shear stress contours average and maximum errors occur. This region was velocity converged grid, additional inaccuracies used to establish the criterion for grid convergence produce an estimated error of approximately 8%. Upon because of the relative sizes of the errors in this region, incrementing the grid density to 800K, the wall shear and due to the existence of large and complex flow stress error reduces to approximately 6.8% compared to disturbances. Also, since the sinus normally develops the 1.48% determined in the velocity independence the atherosclerotic lesions, this region would act as an test. As shown here and also mentioned by Either et al appropriate criterion for grid convergence. [9], wall shear stress grid independence is a more In comparing the velocity fields of the grid accurate criteria but more difficult to achieve. The densities 200K, 400K and 800K elements, qualitatively contour plot of wall shear stress at the 800k level, the curves coincide with minute deviations. Figure 3b, indicated good agreement, especially in Quantitatively, the errors that resulted from increasing close proximity to the bifurcation. Thus, from the the density from 200K to 400K and from 400K to 800K results presented, grid independence using wall shear was found to be approximately 4.8% and 1.48%, stress should be considered when reporting quantitative respectively. However, recognizing that 4.8% is an results for wall shear stress. Caution should be advised adequate level of accuracy, the 400K mesh was when reporting velocity grid independence results for considered grid independent based on velocity. In these values. considering wall shear stress grid independence for the 7th Annual Meeting of the CFD Society of Canada, Halifax NS, May 1999. each region. The values chosen are suggested values Element size and may deviate slightly depending on the relative Region Element size for CCA CCA 1.0 accuracy desired for each region. When implementing Bifucation 0.5 the chosen related element sizes, consideration should ICA 0.6 ECA 0.8 be given to the effects of the upstream and downstream results that may influence convergence in the Table 3: Suggested element size for the carotid artery bifurcation region. REFERENCES 1. Zarins, C.F., Giddens, D.P., Bharadvaj, B.K., Sottiurai, V.S., Mabon, R.F., Glagov, S. Carotid Bifurcation Atherosclerosis. Circ Res 53:502-514, 1983. 2. Glagov, S., Zarins, C., Giddens, D.P., Ku, D.N. Hemodynamics and atherosclerosis: Insights and perspectives gained from studies of human arteries. Arch Pathol Lab Med 112:1018-1031, 1988. 3. Bharadvaj, B.K., Mabon, R.F., Giddens, D.P. Steady flow in a model of the human carotid bifurcation. Part 1-Flow visualization. J. Biomechanics 15:349-362, 1982. 4. Ku, D. Blood Flow in Arteries. Annu. Rev. Fluid Mech. 29:399-434, 1997. 5. Xu, X.Y., Collins, M.W. A review of the numerical analysis of the blood flow in arterial bifurcations. Proc Instn Mech Engrs. Part H Journal of Engineering in Medicine 204:205-216, 1990. 6. Gessner, F.B. Hemodynamic Theories of Atherogenesis. Circ Res 33:259-266, 1973. 7. Smith, R.F., Rutt, B.K., Fox, A.J., Rankin, R.N., Holdsworth, D.W. Geometric characterization of stenosed human carotid arteries. Acad Radiol 3:898-911, 1996. 8. Advanced Scientific Computing Ltd., AEA Technology. CFX-Tfc. MRI data 9. Ethier, C.R., Steinman, D.A., Couch G., Ohja, M. 800K Validation of complex separation patterns in an end-to-side anastomosis model. 1998 Advances in Figure 4: Comparison between MRI data and the computational Bioengineering 39:77-78, 1998. data with the 800K grid Along with the number of elements and nodes, found in Table 1, the values of element sizes are recorded. These values are non-dimensionalized with the inlet diameter (D=2), and as mentioned earlier, they represent the height of a symmetric tetrahedral. With the collaboration of the element sizes and the errors presented in Table 2, approximations for the appropriate element size can be generated for the above mentioned regions. The suggested values are summarized in Table 3. The values are non- dimensionalized with the CCA element size and are chosen based on the errors of the wall shear stress convergence, such that the errors are similar within