Document Sample
					                  7th Annual Meeting of the CFD Society of Canada, Halifax NS, May 1999.


                                      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.


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
geometries published by Smith et al [7] characterize
the carotid artery at varying stages of lumen narrowing,

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

             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
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.



           P                                                      O

          N                                                       M                                                                                     10.00


           L                                                      K

                                                                   I                       8.00

           J                                                                                                                             10.00

                                                                                            4.00                                                 8.00
           H                                                                                           8.00
                                                                                           2.00                                  12.00

                                                             E                             1.00                        15.00

                                                         C                                    0.25
                   D                                                                                   0.50

                                                    B Solid - 800K
                                                                                                                                          Solid - 800k
                                                         Dashed - 400K                                                 4.00
                                                                                                                                          Dotted - 400k
                                                         Dotted - 200K
                    1.0                                                                                       4.00


                          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
                            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.


                                                               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,
                                                               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,
                                                               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

Shared By: