0090-6964/98 $10.50 1 .00 Copyright © 1998 Biomedical Engineering Society
Annals of Biomedical Engineering, Vol. 26, pp. 975–987, 1998 Printed in the USA. All rights reserved.
Finite Element Modeling of Three-Dimensional Pulsatile Flow in the Abdominal Aorta: Relevance to Atherosclerosis CHARLES A. TAYLOR,*† THOMAS J. R. HUGHES,† and CHRISTOPHER K. ZARINS* *Department of Surgery, and †Department of Mechanical Engineering, Stanford University, Stanford, CA (Received 12 November 1997; accepted for publication 28 July 1998)
observed to be devoid of plaque. Ku et al.15 noted strong positive correlations between intimal thickening and the inverse of the maximum wall shear stress, the inverse of the mean wall shear stress, and oscillations in shear characterized by an oscillatory shear index. In the aorta, it is observed that atherosclerotic disease develops first in the abdominal aorta, and is much more common in the abdominal aorta ~i.e., below the diaphragm! than the segment of the aorta above the diaphragm, the thoracic aorta. Roberts et al.34 noted that the abdominal aorta contained by far the most severe aortic atherosclerosis with the greatest involvement occurring below the celiac artery. Glagov et al.6 noted that abdominal aortic atherosclerosis was greater than thoracic aortic atherosclerosis in the majority of cases for males and females, normotensives and hypertensives. Cornhill et al.4 in examining fatty streaks in young subjects noted greater involvement along the lateral and posterior walls of the abdominal aorta. Friedman et al.7 noted increased intimal thickness in regions of low wall shear stress along the lateral walls of the abdominal aortic bifurcation. Moore et al.26 measured intimal thickness in the distal abdominal aorta in subjects with minimal atherosclerotic disease and noted a positive correlation between oscillatory shear and intimal thickness and a negative correlation between mean shear stress and intimal thickening. Clearly, to examine the relationship between vascular disease and hemodynamic conditions detailed quantitative data on flow conditions in the abdominal aorta are required. The pulsatile flow in the abdominal aorta is particularly complex under normal resting conditions as a result of the multiple branches at the level of the diaphragm which deliver blood to the organs in the abdomen. As in the case of carotid artery flow, idealized models representing abdominal aortic anatomy and physiologic conditions have been the basis for most of the investigations into abdominal aortic hemodynamics. Ku et al.16 describe flow visualization techniques applied to study steady flow in an abdominal aorta model under resting,
Abstract—The infrarenal abdominal aorta is particularly prone to atherosclerotic plaque formation while the thoracic aorta is relatively resistant. Localized differences in hemodynamic conditions, including differences in velocity profiles, wall shear stress, and recirculation zones have been implicated in the differential localization of disease in the infrarenal aorta. A comprehensive computational framework was developed, utilizing a stabilized, time accurate, finite element method, to solve the equations governing blood flow in a model of a normal human abdominal aorta under simulated rest, pulsatile, flow conditions. Flow patterns and wall shear stress were computed. A recirculation zone was observed to form along the posterior wall of the infrarenal aorta. Low time-averaged wall shear stress and high shear stress temporal oscillations, as measured by an oscillatory shear index, were present in this location, along the posterior wall opposite the superior mesenteric artery and along the anterior wall between the superior and inferior mesenteric arteries. These regions were noted to coincide with a high probability-of-occurrence of sudanophilic lesions as reported by Cornhill et al. ~Monogr. Atheroscler. 15:13–19, 1990!. This numerical investigation provides detailed quantitative data on hemodynamic conditions in the abdominal aorta heretofore lacking in the study of the localization of atherosclerotic disease. © 1998 Biomedical Engineering Society. @S0090-6964~98!01406-4# Keywords—Hemodynamics, Plaque localization, Computational methods, Parallel computing, Shear stress.
INTRODUCTION Hemodynamic factors are thought to be responsible for the localization of vascular disease in areas of complex flow in the coronary, carotid, abdominal, and femoral arteries. These complex flow regions often occur due to branching, bifurcations, and curvature of the arteries. Zarins et al.42 noted that in the carotid artery, atherosclerotic lesions localize along the outer wall of the carotid sinus region where wall shear stress is low. Conversely, the inner wall of the carotid sinus, an area of rapid, laminar and axial flow, and high wall shear stress, was Address correspondence to Charles A. Taylor, Assistant Professor, Division of Vascular Surgery, Stanford University Medical Center, Suite H3600, Stanford, CA 94305-5450. Electronic mail:
[email protected]
975
976
TAYLOR, HUGHES, and ZARINS
postprandial, and exercise conditions. Moore et al.22 describe pulsatile flow visualization in an idealized abdominal aorta model under resting, postprandial, and exercise conditions. Pedersen et al.28 describe quantitative flow measurements made at three sites in an abdominal aorta flow model using laser Doppler anemometry. Differing physiologic conditions were examined including simulated rest and exercise. Moore and Ku23 describe quantitative flow measurements, using magnetic resonance imaging, for the pulsatile flow in an idealized abdominal aorta model under simulated resting conditions. Noninvasive imaging techniques have also been applied in recent years to measure blood flow in the abdominal aorta in vivo. Maier et al.21 used Doppler ultrasound and magnetic resonance imaging techniques to measure blood flow at a single location between the renal arteries and aortic bifurcation in the abdominal aorta of nine volunteers. Mostbeck et al.27 used MR imaging to measure blood flow patterns in nine volunteers at four locations in the abdominal aorta: 1 cm proximal to the celiac artery, between the superior mesenteric and renal arteries, 1 cm below the renal arteries and 1 cm above the aortic bifurcation. Moore et al.25 measured blood velocity profiles in the abdominal aorta in vivo using magnetic resonance imaging. As in the case of the in vitro model flow studies, current noninvasive imaging technology is limited to the extraction of velocity data at a relatively small number of axial locations. Complete temporal and spatial field variations of velocity and shear stress in the abdominal aorta have not been obtained. In recent years, computational techniques have been used increasingly by researchers seeking to understand vascular hemodynamics. These methods can augment the data provided by in vitro and in vivo methods by enabling a complete characterization of hemodynamic conditions under precisely controlled conditions. Application of these methods to flow in the carotid bifurcation and bypass grafts has provided significant information on vascular hemodynamics. Perktold et al.30 used a finite element method to simulate the pulsatile flow of a Newtonian fluid in a model of a carotid artery bifurcation using a rigid walled approximation. Detailed results on the velocity, pressure, and wall shear stress was presented. Lei et al.18 describe the hemodynamic conditions in a model of a rabbit aortoceliac junction and postulate a role for the wall shear stress gradient in atherogenesis. Numerical methods are well suited to the investigation of phenomena difficult to describe using in vitro techniques including wall compliance, mass transport, particle residence time, and geometric variations. In an investigation of the effect of wall compliance on pulsatile flow in the carotid artery bifurcation, Perktold and Rappitsch32 describe a weakly coupled fluid–structure interaction finite element method for solving for blood flow and vessel
mechanics. Steinman and Ethier37 investigated the effect of wall distensibility in a two-dimensional end-to-side anastomosis. Rappitsch and Perktold33 describe the transport of albumin in a model of a stenosis. Kunov et al.17 proposed a new method for describing particle residence time. Perktold et al.31 examine the effect of bifurcation angle on hemodynamic conditions in carotid artery bifurcation models. Lei et al.19 describe the application of computational methods to the design of end-to-side anastomoses. In contrast to the relatively large number of studies of pulsatile flow in models of the carotid bifurcation and end-to-side anastomosis, there have been few numerical studies of flow in the abdominal aorta. Taylor et al.38 describe the flow in an abdominal aorta model under simulated resting and exercise steady flow conditions. It was noted that a region of flow recirculation and low wall shear stress develop along the posterior wall of the infrarenal abdominal aorta under simulated resting conditions and disappears under simulated moderate and vigorous exercise conditions. Taylor et al.39 describe, qualitatively, the pulsatile flow in a model of an abdominal aorta under simulated resting and moderate exercise conditions. As in the steady flow case, a flow recirculation region was noted at mid-diastole in the infrarenal abdominal aorta under simulated resting conditions. No attempt was made in that study to quantify the hemodynamic conditions under resting or exercise conditions. The computational method described herein was used to quantitatively characterize the hemodynamic conditions under simulated resting pulsatile flow conditions in an idealized model of an abdominal aorta. To date, no numerical results have been published detailing the hemodynamic conditions in the abdominal aorta under pulsatile flow conditions. In particular, although velocity has been measured at a few locations in the abdominal aorta using in vivo and in vitro methods, the spatial variations of wall shear stress have not been reported heretofore. This article details the mean wall shear stress and shear stress oscillations in the abdominal aorta and compares these findings with the observed patterns of sudanophilic lesions as reported by Cornhill et al.4 METHOD A geometric solid model of an idealized abdominal aorta, shown in Fig. 1, was constructed using a custom software system, ‘‘The Stanford Virtual Vascular Laboratory’’ developed to aid in the solution of blood flow problems.39,40 This system was developed using a knowledge based, object oriented programming language and combines geometric solid modeling, automatic finite element mesh generation, multiphysics finite element methods, and scientific visualization capabilities.43–46 The anatomic dimensions of the idealized abdominal
Finite Element Modeling of Flow in the Abdominal Aorta
977
FIGURE 1. Abdominal aorta model with branches identified. The model is not symmetric about the midsagittal plane. Note the tapering of the aorta from the diaphragm to the aortic bifurcation.
aorta model were obtained primarily from Moore et al.22 In the model utilized for the present investigations, the aorta tapers from a circular cross section with a diameter of 2.54 cm at the diaphragm to a circular cross section with diameter of 1.72 cm at the inferior mesenteric artery and then tapers uniformly to an elliptical cross section with a major axis of 1.53 cm and minor axis of 1.3 cm at the aortic bifurcation. It should be noted that the model constructed is not symmetric about the midsagittal plane, but rather includes the feature that the left renal artery is located inferior to the right renal artery. A finite element mesh with 268,563 tetrahedral elements and 58,151 nodes was generated using an automatic mesh generator.36 The surface of the finite element mesh above the aortic bifurcation is displayed in Fig. 2. Note the
FIGURE 2. Closeup of abdominal aorta finite element mesh above the aortic bifurcation. Note the greater mesh refinement used for the infrarenal aorta and the branch arteries.
refinement of the mesh in the region of the aorta below the renal arteries and along the branch vessels. Under resting conditions, approximately 70% of the blood that enters the abdominal aorta is extracted by the celiac, superior mesenteric, and renal arteries. The majority of the remaining 30% flows down the infrarenal segment through the bifurcation into the legs. In contrast to steady flow studies, in vitro and computational studies of pulsatile blood flow require the specification of flow waveforms as well as mean flow rates. The flow rates as a function of time used in the present study for the inflow and branch vessels are shown in Fig. 3. The flow rate waveforms shown in Fig. 3 for the renal artery and iliac artery are those for each of the right and left renal and iliac arteries, respectively. The suprarenal and infrarenal flow rate time functions and the celiac, superior mesenteric, renal, and inferior mesenteric mean flow rates were obtained from Moore and Ku.23 The flow rate waveforms in the celiac, superior mesenteric, renal, and inferior mesenteric arteries were computed to conserve mass and yield the mean flow rates given by Moore et al.23 Note that the abdominal aorta inflow and renal artery outflows are always positive, whereas the iliac flow waveform exhibits the triphasic character of in vivo measurements of infrarenal aortic flow. Based on the volume flow waveforms shown in Fig. 3, pulsatile flow velocity boundary conditions, derived from the Womersley theory, were prescribed for the inflow boundary and all outflow boundaries excluding the left and right iliac outflow boundaries, where zero pressure boundary conditions were prescribed.41 It should be noted that the
978
TAYLOR, HUGHES, and ZARINS
FIGURE 3. Flow rate time functions for resting conditions. Note the triphasic nature of the flow waveforms and the different ordinate scales for each plot.
Womersley boundary conditions are time dependent, axisymmetric velocity profiles at the inlet and outlet boundaries. Computational modeling of blood flow requires solving, in the general case, three dimensional, transient flow equations in deforming blood vessels. The appropriate framework for problems of this type is the arbitrary
Lagrangian–Eulerian ~ALE! description of continuous media in which the fluid and solid domains are allowed to move to follow the distensible vessels and deforming fluid domain.10 The vessel diameter change during the cardiac cycle is observed to be approximately 5%–10% in most of the major arteries, and, as a first approximation, the vessel
Finite Element Modeling of Flow in the Abdominal Aorta
walls are often treated as being rigid. In addition, in diseased vessels which are often the subject of interest, the arteries are even less compliant and wall motion is further reduced. The assumption of zero wall motion was utilized for the computations presented herein. Under this assumption, the ALE description of incompressible flow in a deforming fluid domain reduces to the Eulerian description of a fixed spatial domain. The strong form of the problem governing incompressible, Newtonian fluid flow in a fixed domain consists of the Navier–Stokes equations and suitable initial and boundary conditions. Direct, analytic solutions of these equations are not available for complex domains and numerical methods must be used. The finite element method has been the most widely used numerical method for solving the equations governing blood flow. The finite element method employed in the present investigation is based on the theory of stabilized finite element methods developed by Hughes and colleagues.2,11,12 The method is described in detail by Taylor et al. as it is implemented in the commercial finite element program, SPECTRUM™.39 The essential features of stabilized methods in the context of incompressible flows are the simultaneous stabilization of the advection operator and the circumvention of the Babuska– Brezzi inf–sup condition restricting the use of many convenient interpolations, including the linear velocity, linear pressure interpolations used in the present work. The basic idea is to augment the Galerkin finite element formulation with a least-squares form of the residual, including appropriate stabilization parameters. These stabilization parameters are designed so that the method achieves exact solutions in the case of one-dimensional model problems involving, for example, the steady advection-diffusion equation. A semidiscrete formulation with a second-order accurate time-stepping algorithm is utilized resulting in a nonlinear algebraic problem in each time step. This nonlinear problem is linearized and the resulting linear systems of equations are iteratively solved using conjugate gradient and matrix-free GMRES solvers to reduce memory requirements.13,35 As the flow computations were performed on a parallel computer, the computational domain was divided into 16 subdomains using a graph partitioning method.14 The nonlinear evolution equations were solved for the velocity and pressure fields over three cardiac cycles with 200 time steps per cardiac cycle on an IBM SP-2 parallel computer. Analyses run with a greater number of time steps showed no observable differences in the solution. The method employed has been shown via analytical and experimental validation studies to yield accurate solutions for pulsatile flow problems.39 In addition to the velocity field, quantitative values of the magnitude of the surface tractions are also of interest.
979
The traction vector can be computed from the stress tensor, s, and surface normal vector, n, by t5 s n and then the surface traction vector, ts , defined as the tangential component of the traction vector, can be computed from ts 5t2(t–n)n. We can define the mean shear stress, t mean , a scalar quantity, as the magnitude of the time-averaged surface traction vector as
t mean5
UE U 1 T
T
0
ts dt ,
~1!
and define t mag , another scalar quantity, as the timeaveraged magnitude of the surface traction vector as
t mag5
1 T
E
T
0
u ts u dt.
~2!
Following He and Ku,8 we define the oscillatory shear index ~OSI! as OSI5
S
D
1 t mean 12 . 2 t mag
~3!
Note that in the special case of the steady, uniaxial flow of a Newtonian fluid in a circular cylinder, the magnitude of the surface traction reduces to the wall shear stress, t w 54 m Q/ p r, 3 where m is the viscosity, Q is the flow rate, and r is the lumen radius. However, due to the complex flow fields in the vascular system, the surface traction is a more general, and useful, measure of surface forces as it is a vector quantity.
RESULTS A contour slice of velocity magnitude and a vector slice of the velocity field is displayed in Fig. 4 under simulated resting conditions at three times in the cardiac cycle along the midsagittal plane of the aorta. Note that a vortex develops along the posterior wall of the aorta in late systole and is present throughout diastole. The posterior wall of the aorta is a region with consistently low flow velocities relative to the anterior wall. Also note the orientation of the vector field towards the anterior wall, observed even at peak systole, as a result of the celiac and superior mesenteric outflows. Figure 5 displays the mean surface traction vectors along the anterior wall of the abdominal aorta. It is observed that the direction of the mean surface traction vectors is primarily in the superior to inferior direction along most of the surface of the aorta with the exception of the neighborhoods of the branch vessels. It is of note
980
TAYLOR, HUGHES, and ZARINS
FIGURE 5. Mean surface traction vectors along the anterior wall of abdominal aorta. Note the circumferential orientation of the mean surface traction vectors near the branch vessels.
FIGURE 4. Midplane slice of the abdominal aorta model displaying contours of velocity magnitude and the velocity vector field at „a… peak systole, „b… end systole, and „c… middiastole. Note the low flow velocities and flow recirculation along the posterior wall of the infrarenal aorta.
that the direction of the mean surface traction vectors is from posterior to anterior in the neighborhood of the ostia of the celiac and superior mesenteric arteries but from anterior to posterior at the renal ostia. Figure 6 displays the mean surface traction vectors along the posterior wall of the abdominal aorta. It is noted that the direction of the surface traction vectors is again primarily from superior to inferior with the exception of the surface traction vectors on the aorta wall in the neighborhood of the branch vessels. It is noted that for the two sites of relatively low mean wall shear stress along the posterior wall, namely, the opposite of the celiac and superior mesenteric arteries and distal to the renal arteries, that the direction of the mean surface traction vectors is predominantly circumferential. Furthermore, opposite the celiac and superior mesenteric arteries the mean surface traction vectors is from posterior to anterior whereas the opposite is true distal to the renal arteries where the
mean surface traction vectors are from anterior to posterior. It is also of interest that the mean surface traction vectors form concentric rings around the ostia of the renal vessels as is apparent for the left renal artery in Fig. 6. Figure 7 displays the mean surface traction vectors at the aortic bifurcation. Note that for the iliac vessel wall opposite the aortic bifurcation flow divider, the direction of the mean surface traction vectors is from anterior to posterior along the anterolateral wall and from posterior to anterior along the posterolateral wall. Contour plots of mean wall shear stress, t mean , are displayed in Fig. 8 from anterior and posterior views. Note that under resting conditions a region of low mean wall shear stress develops below the renal artery branch and that, furthermore, shear stresses are high on the inside wall of the aortic bifurcation and low on the lateral walls.
FIGURE 6. Mean surface traction vectors along the posterior wall of abdominal aorta. Note the circumferential orientation of the mean surface traction vectors along the posterior wall of the aorta in the neighborhood of the renal arteries.
Finite Element Modeling of Flow in the Abdominal Aorta
FIGURE 7. Mean surface traction vectors at the aortic bifurcation. Note the orientation of the mean surface traction vectors along the lateral wall of the right iliac artery in the neighborhood of the aortic bifurcation.
981
Mean wall shear stress, t mean , is plotted in Fig. 9 as a function of arc length along the anterior wall. Note that the gaps in the plots of shear stress correspond to the location of the branch vessels. Observe that mean shear stress increases significantly in the neighborhood of the celiac, superior mesenteric and inferior mesenteric vessels, and that the level of shear stress in the distal infrarenal aorta exceeds the shear stress at the level of the diaphragm. Mean wall shear stress is plotted in Fig. 10 as a function of arc length along the posterior wall. The mean shear stress approaches zero opposite the celiac and superior mesenteric vessels and approximately 2 cm distal to the left renal artery. The mean shear stress increases in the distal infrarenal abdominal aorta to a level exceeding that of the diaphragm.
FIGURE 8. Mean shear stress contours along „a… the anterior wall and „b… the posterior wall. Note the regions of low shear stress opposite the celiac and superior mesenteric arteries, in the infrarenal abdominal aorta and along the lateral wall of the aortic bifurcation.
982
TAYLOR, HUGHES, and ZARINS
FIGURE 9. Mean wall shear stress along the anterior wall from the diaphragm to the aortic bifurcation. Note the low mean shear stress in the infrarenal aorta between the superior and inferior mesenteric arteries „arc length of ;9 cm….
FIGURE 11. Oscillatory shear index „OSI… along the anterior wall from the diaphragm to the aortic bifurcation. Note the elevated OSI in the infrarenal aorta between the superior and inferior mesenteric arteries „arc length of ;9 cm….
The OSI is plotted as a function of arc length in Figs. 11 and 12 along the anterior and posterior walls, respectively. Examination of the OSI along the anterior wall reveals a sharp increase from zero in the suprarenal aorta to approximately 0.3 immediately distal to the left renal
artery and then a gradual decline to a value of approximately 0.05 at the level of the aortic bifurcation. The OSI along the posterior wall increases from zero in the suprarenal aorta to approximately 0.4 opposite the superior mesenteric artery, decreases to less than 0.1 at the level of the renal arteries, increases to a maximum value
FIGURE 10. Mean wall shear stress along the posterior wall from the diaphragm to the aortic bifurcation. Note the low mean shear stress opposite the celiac and superior mesenteric arteries „arc length of ;5.5 cm… and in the infrarenal abdominal aorta „arc length of ;9 cm….
FIGURE 12. Oscillatory shear index „OSI… along the posterior wall from the diaphragm to the aortic bifurcation. Note the elevated OSI opposite the celiac and superior mesenteric arteries „arc length of ;5.5 cm… and in the infrarenal abdominal aorta „arc length of ;9 cm….
Finite Element Modeling of Flow in the Abdominal Aorta
of 0.47 distal to the renal arteries, and then decreases to a value of approximately 0.1 at the level of the aortic bifurcation. DISCUSSION Numerical solution of pulsatile flow in the abdominal aorta reveals complex hemodynamic patterns under simulated resting conditions. The present approach to characterize abdominal aorta hemodynamics enables the extraction of complete spatial and temporal variations in field quantities, including velocity and shear stress, for correlation with observed patterns of vascular disease. The significance of the findings regarding velocity field and shear stresses is subsequently described, followed by an examination of the assumptions employed in the present investigation. Two features of the velocity field in the infrarenal abdominal aorta are worth noting. First, it is observed that the velocity field in the infrarenal abdominal aorta is directed preferentially towards the anterior wall. Second, a flow recirculation region appears along the posterior wall of the abdominal aorta in late systole and is present throughout diastole. This recirculation region extends from the level of the renal arteries to the inferior mesenteric artery. Moore et al.,22 using an in vitro model on which the present numerical study was based, noted that, under resting conditions, vortices appeared at the level of the renal arteries and propagated through the infrarenal aorta. Pedersen et al.,28 using laser Doppler anemometry, observed an orientation of the velocity profile in the posterior to anterior direction with the velocities being highest at the anterior wall. Moore and Ku,23 using magnetic resonance imaging, noted extensive flow reversal along the posterior wall and an orientation of the velocity profile in the posterior to anterior direction in the infrarenal portion of an in vitro model. Complex velocity profiles were also observed at the aortic bifurcation with flow reversal and oscillation most prominent along the lateral walls proximal to the bifurcation. Maier et al.,21 using Doppler ultrasound and magnetic resonance imaging techniques in humans, again noted greater anterior flow velocity as compared to posterior velocity during early diastole and sometimes during peak systole. Mostbeck et al.,27 using MR imaging of human subjects, observed that in the infrarenal aorta the mean duration of retrograde flow was greatest along the posterior wall of the aorta. Greater flow velocity along the anterior wall relative to that along the posterior wall and a flow recirculation zone along the posterior wall of the infrarenal abdominal aorta was observed in the computed velocity field. This finding is consistent with prior in vitro and in vivo studies of flow in the abdominal aorta. It should be noted, however, that the extent of the flow recirculation region
983
under resting conditions has not been reported previously due to the fact that velocity measurements in the infrarenal aorta have been restricted to a single cross section. The spatial distribution of the direction of the mean surface traction vectors provides qualitative information on the forces acting on the luminal surface. For most of the abdominal aorta the direction of the mean shear stress vectors is primarily in the superior to inferior direction, but four exceptions are noteworthy. First, in the neighborhood of the ostia of the celiac and superior mesenteric arteries the mean shear stress vectors are radially oriented towards the ostia. Second, in the neighborhood of the renal arteries the shear stress vectors form concentric rings around the renal ostia. Third, for the two sites of relatively low mean wall shear stress along the posterior wall, opposite the celiac and superior mesenteric arteries and distal to the renal arteries, the mean surface shear stress vectors are circumferentially oriented. Finally, the direction of the mean surface traction vectors is from anterior to posterior along the anterolateral wall and from posterior to anterior along the posterolateral wall for the iliac vessel wall opposite the aortic bifurcation flow divider. It should be noted that experimental methods which only measure the axial component of velocity do not permit observation of these nonaxial shear stress vectors. Examination of the spatial distribution of the magnitude of the mean shear stress vector reveals several new findings. It is noted that under resting conditions regions of low mean wall shear stress develop along the lateral and posterior walls opposite the celiac and superior mesenteric arteries and below the renal arteries. It is also noted that mean shear stress increases in the distal abdominal aorta from a level above the inferior mesenteric artery to the aortic bifurcation. This is presumably due to the tapering of the aorta in this region. Although shear stresses were reported by Moore et al.26 for two locations in the infrarenal abdominal aorta, the spatial distribution of shear stress magnitude has not been quantified previously. Finally, as noted by Friedman et al.7 and more recently by Moore et al.26 mean shear stresses are high on the inside wall of the aortic bifurcation and low on the lateral walls. The plots of shear stress and oscillatory shear index along the anterior and posterior walls provide quantitative data on the spatial variation of shear down the length of the abdominal aorta. The mean shear stress in the abdominal aorta at the level of the diaphragm was calculated to be approximately 1.2 dynes/cm2. Moore et al.26 measured the shear stress at this location to be 1.360.6 dynes/cm2 in an in vitro study with approximately the same flow conditions and model dimensions as those used in this study. Along the anterior wall the ostia of the celiac, superior mesenteric, and inferior mesenteric arteries are sites of relatively high mean shear
984
TAYLOR, HUGHES, and ZARINS
stress and low oscillatory shear index. As discussed subsequently, this is believed to be due to an inadequate modeling assumption, incorporated into the idealized model, whereby the junctures between the branch vessels and the aorta are not smoothly tapered. The anterior wall of the infrarenal aorta between the superior mesenteric and inferior mesenteric arteries is a site of low shear stress and high oscillatory shear index. Along the posterior wall, the mean shear stress is low and the oscillatory shear index high opposite the superior mesenteric artery as well as approximately 2 cm distal to the left renal artery in the region noted for flow recirculation. The distal abdominal aorta below the inferior mesenteric artery and proximal to the bifurcation is a site of relatively high mean shear stress and moderate values of the oscillatory shear index. As noted previously, vascular disease is more prominent in the abdominal aorta as compared with the thoracic aorta. The localization of early atherosclerotic lesions was examined by Cornhill et al.4 by staining with Sudan IV and using image processing techniques to identify sudanophilic lesions in the cadaveric aortas of young accident victims. It was noted that the highest probability-of-occurrence of sudanophilia were associated with the inflow tracts of the celiac, superior mesenteric, renal, and inferior mesenteric ostia. In the present study, along the aortic wall, these ostia were not observed to be sites of low mean shear stress, yet regions of low mean shear stress were noted within the branch vessels. Again, note that the model does not incorporate the tapering of the branch vessels in the immediate vicinity of the ostia. Cornhill et al.4 also noted a region of high probability-of-occurrence of sudanophilia along the posterolateral wall of the infrarenal aorta and along the anterior wall between the superior and inferior mesenteric arteries. Cornhill et al.4 state: ‘‘A major portion of areas of high probability of sudanophilia, namely the extensive region along the dorso-lateral aortic surface from the ductus scar to the aorto-iliac bifurcation and the ventral surface of the abdominal aorta midway between the renal and inferior mesenteric ostia, does not occur in regions expected to experience unusual hemodynamic stresses.’’ The results presented herein clearly demonstrate otherwise. Namely, these regions experience low mean wall shear stress and high oscillatory shear index. Finally, Cornhill et al.4 show areas of high probability-ofoccurrence of sudanophilia along the posterior wall at the level of the celiac and superior mesenteric arteries. As noted previously, this is also a region of low mean wall shear stress and high oscillatory shear index. Clearly,
direct correlative studies need to be performed between hemodynamic conditions and plaque localization in the abdominal aorta, but the results presented herein demonstrate that low mean wall shear stress and high oscillatory shear index occur in regions noted to have a high probability-of-occurrence of sudanophilic lesions. The major assumptions employed in this investigation involve the anatomic dimensions of the model, the flow conditions, a rigid wall approximation, and Newtonian viscosity. The utilization of these assumptions, as for prior in vitro investigations7,23–26,28 is believed to constitute a reasonable first approximation to the actual hemodynamic conditions in the abdominal aortas of normal, healthy, humans at rest. The anatomic dimensions employed in the present study were based on those reported by Moore et al.22 and include a tapering of the abdominal aorta and lumbar curvature. Potential limitations of this representation include a relatively abrupt change in dimensions at the branch ostia and the lack of curvature of the branch vessels. Because of a lack of available anatomic data, the tapering of the branch vessels in the neighborhood of the ostia on the aortic wall was not incorporated in the geometric model. This lack of branch vessel taper is expected to significantly alter flow conditions in the neighborhood of the ostia and could explain the relatively high shear stresses at these locations. The effect of this modeling assumption on shear stresses in the neighborhood of the ostia needs to be investigated. Recent work by Caro et al.3 argues that branch vessel tapering as well as the nonplanar curvature of branch vessels beyond the aortic ostia may have a significant influence on hemodynamic conditions. Curvature of the branch vessels was also not considered in the present study, but the importance of including this feature of actual human anatomy in aortic flow studies merits further investigation. The physiologic conditions employed in the present study were based on those reported by Moore et al.22 and Moore and Ku,23 and were chosen to represent resting abdominal aortic flow conditions, yet represent the lower end of resting flow rates. Examination of greater resting flow rates, while maintaining resting flow distribution, merits further study. The physiologic conditions differ from those reported by Moore et al.22 and Moore and Ku23 in that the numerical method utilized herein requires the specification of the time-varying outflow for the celiac, superior mesenteric, renal, and inferior mesenteric vessels. In contrast, only mean flow was determined and reported by Moore et al.22 and Moore and Ku23 for these branch vessels. Three features of abdominal aortic flow waveforms, consistently observed in normal subjects under resting conditions, were used in the present study. First, the flow rate is observed to be positive throughout the cardiac cycle at the entrance to the abdominal aorta at the level
Finite Element Modeling of Flow in the Abdominal Aorta
of the diaphragm. It was noted by Holenstein and Ku9 that 32 of 35 young, normal subjects had no reverse flow as measured by Doppler ultrasound approximately 2–4 cm proximal to the superior mesenteric artery. Bogren and Buonocore1 used MR cine velocity measurements to quantify separately antegrade and retrograde flow and noted low retrograde and positive net flow in the distal descending thoracic aorta. It is important to note that, due to the nature of pulsatile flow, some retrograde flow can appear near the wall during early diastole even in the presence of a positive net flow. In a further study of abdominal aortic flow, Moore et al.,25 using MR velocimetry, noted only a slight flow reversal at end systole in the supraceliac aorta in four healthy volunteers. The second feature of abdominal aortic flow observed in normal subjects under resting conditions is the triphasic flow waveform in the infrarenal aorta. This has been noted in numerous studies employing doppler ultrasound and MRI methods for velocity quantification. Holenstein and Ku9 noted that 34 of 35 young, normal subjects had reverse flow at end systole as measured by Doppler ultrasound. Bogren and Buonocore1 and Moore et al.25 also noted reverse flow at end systole using MR cine velocity measurements in normal, healthy subjects. The third feature of abdominal aortic flow observed in normal subjects under resting conditions is a positive net renal artery flow throughout the cardiac cycle and relatively high diastolic flow. Holenstein and Ku9 noted positive left renal artery flow and an absence of reverse flow using Doppler ultrasound. Bogren and Buonocore1 measured blood flow in the abdominal aorta above and below the renal arteries and noted a total renal artery flow of approximately 1 L/min. Pelc et al.29 used PC-MRI techniques to directly measure renal blood flow in a normal healthy subject and observed flow rates in excess of 0.6 L/min throughout the cardiac cycle. The determination of appropriate boundary conditions and flow waveforms for investigations of abdominal aorta hemodynamic conditions merits further attention. In addition to the specification of flow waveforms, assumptions are required for the variation in velocity at the boundaries where flow is specified. The velocity profiles at the inlet to the abdominal aorta at the level of the diaphragm, and the celiac, superior mesenteric, renal, and inferior mesenteric artery outflows were specified using Womersley theory for pulsatile flow, resulting in an axisymmetric velocity profile. Moore et al.25 confirmed the validity of the assumption of axisymmetric flow in the supraceliac abdominal aorta. This condition was also utilized in the model flow studies reported by Pedersen et al.28 and Moore and Ku.23 Regarding the outflow velocity profiles, the length of the branch vessels was chosen to minimize the effect of the outlet boundary condition on the flow in the aorta. It is likely that the assumption of planar branch vessels in the vicinity of the
985
ostia on the aorta has a greater effect on the flow conditions in the aorta than the axisymmetric outlet velocity profile, but this was not examined in the present study. The effect of aortic compliance was not considered in the present investigation. Prior investigations of flow in deformable models have shown that, in general, incorporating the effect of compliance alters the magnitude of shear stress, but does not change the locations of low shear and high shear regions. Duncan et al.5 examined flow in a deformable model of an aortic bifurcation and showed that at the outer wall of the aortic bifurcation the shear rates were reduced in the deformable model as compared to a rigid model. The radial deformation measured in this model flow study was on the average of 2%–4%. In an investigation of the effect of wall compliance on pulsatile flow in the carotid artery bifurcation, Perktold and Rappitsch32 describe a weakly coupled fluid-structure interaction finite element method for solving for blood flow and vessel deformation. They concluded that the wall shear stress magnitude decreases by approximately 25% in the distensible model as compared to the rigid model, yet the overall effect on the velocity field was relatively minor. Steinman and Ethier,37 based on a study of a compliant two-dimensional end-to-side anastomosis, concluded that the effects of wall distensibility on the flow field are less pronounced as compared to changes in anatomic dimensions or physiologic conditions. It should be noted that although neglecting wall distensibility may not have a major effect on the primary flow field, the incorporation of wall mechanics in these studies is important for other reasons including the description of the stress environment within the vessel walls and the interaction between deformability and mass transport phenomena. In summary, the use of rigid models for flow studies in idealized models of the human abdominal aorta should be viewed as a first approximation. A Newtonian constitutive model for viscosity was employed in the present investigation. It is generally accepted that this is a reasonable first approximation to the behavior of blood flow in large arteries. Perktold et al.31 examined non-Newtonian viscosity models for simulating pulsatile flow in carotid artery bifurcation models. They concluded that the shear stress magnitudes predicted using non-Newtonian viscosity models resulted in differences on the order of 10% as compared with Newtonian models. Future studies of abdominal aorta hemodynamic conditions should assess the effects of non-Newtonian rheologic models. The accuracy of the finite element method for solving the governing equations was not examined quantitatively in the present study due to a lack of published data. In other models Taylor et al.40 examined the accuracy of this method employed for blood flow computations. Numerical results were found to be in excellent agreement
TAYLOR, HUGHES, and ZARINS
986
with Womersley theory and with laser Doppler anemometry velocity data obtained by Loth20 for steady and pulsatile flow in a model of an end-to-side anastomosis. It was not possible to examine further refinements in the finite element mesh employed due to the limitation of computational resources. However, the stabilized finite element method employed has been shown to exhibit excellent coarse grid accuracy.40 It is not expected that further refinements in the mesh will result in substantial quantitative differences in the computed solution. CONCLUSIONS The characterization of the temporal and spatial variations of the velocity field in an idealized model of the abdominal aorta enables new insights into the wall shear stresses acting on the luminal surface under resting conditions and should help resolve the question of which hemodynamic factors correlate with the disease localization patterns observed in cadavers and in vivo. Regions of low mean wall shear stress and high oscillatory shear index were noted along the anterior wall between the superior mesenteric and inferior mesenteric arteries and along the posterior wall opposite the superior mesenteric artery and approximately 2 cm distal to the renal arteries. These locations were also characterized by Cornhill et al.4 as having high probability-of-occurrence of sudanophilic lesions. The numerical studies described herein provide further impetus to examine abdominal aorta hemodynamics in vivo using magnetic resonance imaging techniques, and to assess the validity of the assumptions made regarding vascular anatomy and physiologic conditions, blood rheology, and vessel mechanics. ACKNOWLEDGMENTS The authors gratefully acknowledge the assistance of Jeff Buell, Mary Draney, David Parker, and Arthur Raefsky. Software was provided by Centric Engineering Systems, Inc., TechnoSoft, Inc., XOX Inc., and Rensselaer Polytechnic Institute’s Scientific Computing Research Center. Hewlett-Packard provided a model 755/125 workstation used in this research.
REFERENCES 1
Bogren, H. G., and M. H. Buonocore. Blood flow measurements in the aorta and major arteries with MR velocity mapping. J. Magn. Reson. Imaging 4:119–130, 1994. 2 Brooks, A. N., and T. J. R. Hughes. Streamline upwind Petrov–Galerkin formulation for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. Comput. Methods Appl. Mech. Eng. 32:199–259, 1981. 3 Caro, C. G., D. J. Doorly, M. Tarnawski, K. T. Scott, Q.
Long, and C. L. Dumoulin. Non-planar curvature and branching of arteries and non-planar type flow. Proc. R. Soc. London, Ser. B 452:185–197, 1996. 4 Cornhill, J. F., E. E. Herderick, and H. C. Stary. Topography of human aortic sudanophilic lesions. Monogr. Atheroscler. 15:13–19, 1990 5 Duncan, D. D., C. B. Bargeron, S. E. Borchardt, O. J. Deters, S. A. Gearhart, F. F. Mark, and M. H. Friedman. The effect of compliance on wall shear in casts of a human aortic bifurcation. J. Biomech. Eng. 112:183–188, 1990. 6 Glagov, S., D. A. Rowley, and R. I. Kohut. Atherosclerosis of human aorta and its coronary and renal arteries. Arch. Pathol. 72:82–95, 1961. 7 Friedman, M. H., G. M. Hutchins, and C. B. Bargeron. Correlation between intimal thickness and fluid shear in human arteries. Atherosclerosis (Berlin) 39:425–436, 1981. 8 He, X., and D. N. Ku. Pulsatile flow in the human left coronary artery bifurcation: Average conditions. J. Biomech. Eng. 118:74–82, 1996. 9 Holenstein, R., and D. N. Ku. Reverse flow in the major infrarenal vessels: A capacitive phenomenon. Biorheology 25:835–842, 1988. 10 Hughes, T. J. R., W. K. Liu, and T. K. Zimmermann. Lagrangian–Eulerian finite element formulation for incompressible viscous flows. Comput. Methods Appl. Mech. Eng. 29:329–349, 1981. 11 Hughes, T. J. R., L. P. Franca, and M. Balestra. A new finite element method for computational fluid dynamics: V. Circumventing the Babuska-Brezzi condition: A stable Petrov– Galerkin formulation of the Stokes problem accommodating equal order interpolations. Comput. Methods Appl. Mech. Eng. 59:85–99, 1986. 12 Hughes, T. J. R., L. P. Franca, and G. M. Hulbert. A new finite element method for computational fluid dynamics: VIII. The Galerkin/least squares method for advective-diffusive equations. Comput. Methods Appl. Mech. Eng. 73:173–189, 1989. 13 Johan, Z., T. J. R. Hughes, and F. Shakib. A globally convergent matrix-free algorithm for implicit time-marching schemes arising in finite element analysis in fluids. Comput. Methods Appl. Mech. Eng. 87:281–304, 1991. 14 Karypis, G., and V. Kumar. A fast and high quality multilevel scheme for partitioning irregular graphs. SIAM J. Sci. Comput. ~to appear!. 15 Ku, D. N., D. P. Giddens, C. K. Zarins, and S. Glagov. Pulsatile flow and atherosclerosis in the human carotid bifurcation. Positive correlation between plaque location and low oscillating shear stress. Arteriosclerosis (Dallas) 5:293–302, 1985. 16 Ku, D. N., S. Glagov, J. E. Moore, and C. K. Zarins. Flow patterns in the abdominal aorta under simulated postprandial and exercise conditions: an experimental study. J. Vasc. Surg. 9:309–316, 1989. 17 Kunov, M. J., D. A. Steinman, and C. R. Ethier. Particle volumetric residence time calculations in arterial geometries. J. Biomech. Eng. 118:158–164, 1996. 18 Lei, M., C. Kleinstreuer, and G. A. Truskey. Numerical investigations and prediction of atherogenic sites in branching arteries. J. Biomech. Eng. 117:350–357, 1995. 19 Lei, M., J. P. Archie, and C. Kleinstreuer. Computational design of a bypass graft that minimizes wall shear stress gradients in the region of the distal anastomosis. J. Vasc. Surg. 25:637–646, 1997. 20 Loth, F. Velocity and Wall Shear Stress Measurements inside a Vascular Graft Model under Steady and Pulsatile Flow
Finite Element Modeling of Flow in the Abdominal Aorta Conditions. Atlanta, GA: Georgia Institute of Technology, PhD Dissertation, 1993. 21 Maier, S. E., D. Meier, P. Boesiger, U. T. Moser, and A. Vieli. Human abdominal aorta: Comparative measurements of blood flow with MR imaging and multigated Doppler US. Radiology 171:487–492, 1989. 22 Moore, J. E., D. N. Ku, C. K. Zarins, and S. Glagov. Pulsatile flow visualization in the abdominal aorta under differing physiologic conditions: Implications for increased susceptibility to atherosclerosis. J. Biomech. Eng. 114:391–397, 1992. 23 Moore, J. E., and D. N. Ku. Pulsatile velocity measurements in a model of the human abdominal aorta under resting conditions. J. Biomech. Eng. 116:337–346, 1994. 24 Moore, J. E., and D. N. Ku. Pulsatile velocity measurements in a model of the human abdominal aorta under simulated exercise and postprandial conditions. J. Biomech. Eng. 116:107–111, 1994. 25 Moore, J. E., S. E. Maier, D. N. Ku, and P. Boesiger. Hemodynamics in the abdominal aorta: A comparison of in vitro and in vivo measurements. J. Appl. Physiol. 76:1520– 1527, 1994. 26 Moore, J. E., C. Xu, S. Glagov, C. K. Zarins, and D. N. Ku. Fluid wall shear stress measurements in a model of the human abdominal aorta: Oscillatory behavior and relationship to atherosclerosis. Atherosclerosis (Berlin) 110:225–240, 1994. 27 Mostbeck, G. H., M. C. Dulce, G. R. Caputo, E. Proctor, and C. B. Higgins. Flow pattern analysis in the abdominal aorta with velocity-encoded cine MR imaging. J. Magn. Reson. Imaging 3:617–623, 1993. 28 Pedersen, E. M., H. W. Sung, A. C. Burlson, and A. P. Yoganathan. Two-dimensional velocity measurements in a pulsatile flow model of the normal abdominal aorta simulating different hemodynamic conditions. J. Biomech. 26:1237– 1247, 1993. 29 Pelc, N. J., F. G. Sommer, K. C. P. Li, T. J. Brosnan, R. J. Herfkens, and D. R. Enzmann. Quantitative magnetic resonance flow imaging. Magn. Reson. Q. 10:125–147, 1994. 30 Perktold, K., M. Resch, and R. Peter. Three-dimensional numerical analysis of pulsatile flow and wall shear stress in the carotid artery bifurcation. J. Biomech. 24:409–420, 1991. 31 Perktold, K., R. O. Peter, M. Resch, and G. Langs. Pulsatile non-Newtonian flow in three-dimensional carotid bifurcation models: A numerical study of flow phenomena under different bifurcation angles. J. Biomed. Eng. 13:507–515, 1991.
32
987
Perktold, K., and G. Rappitsch. Computer simulation of local blood flow and vessel mechanics in a compliant carotid artery bifurcation model. J. Biomech. 28:845–856, 1995. 33 Rappitsch, G., and K. Perktold. Pulsatile albumin transport in large arteries: A numerical simulation study. J. Biomech. Eng. 118:511–518, 1996. 34 Roberts, J. C., C. Moses, and R. H. Wilkins. Autopsy studies in atherosclerosis: I. Distribution and severity of atherosclerosis in patients dying without morphologic evidence of atherosclerotic catastrophe. Circulation 20:511–519, 1959. 35 Saad, Y., and M. H. Shultz. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM (Soc. Ind. Appl. Math.) J. Sci. Stat. Comput. 7:856– 869, 1986. 36 Shephard, M. S., and M. K. Georges. Automatic threedimensional mesh generation by the finite octree technique. Int. J. Numer. Methods Eng. 32:709–749, 1991. 37 Steinman, D. A., and C. R. Ethier. The effect of wall distensibility on flow in a two-dimensional end-to-side anastomosis. J. Biomech. Eng. 116:294–301, 1994. 38 Taylor, C. A., B. I. Tropea, T. J. R. Hughes, and C. K. Zarins. Effect of graded exercise on aortic wall shear stress. Surg. Forum XLVI:331–334, 1995. 39 Taylor, C. A., T. J. R. Hughes, and C. K. Zarins. Computational investigations in vascular disease. Comput. Phys. 10:224–232, 1996. 40 Taylor, C. A., T. J. R. Hughes, and C. K. Zarins. Finite element modeling of blood flow in arteries. Comput. Methods Appl. Mech. Eng. 158:155–196, 1998. 41 Womersley, J. R. Method for the calculation of velocity, rate of flow, and viscous drag in arteries when the pressure gradient is known. J. Physiol. (London) 127:553–63, 1955. 42 Zarins, C. K., D. P. Giddens, B. K. Bharadvaj, V. S. Sottiurai, R. F. Mabon, and S. Glagov. Carotid bifurcation atherosclerosis: Quantitative correlation of plaque localization with flow velocity profiles and wall shear stress. Circ. Res. 53:502–514, 1983. 43 SHAPES™ Reference Manual, XOX, Inc., 1994. 44 Adaptive Modeling Language™ Theory Manual, Technosoft, Inc., 1995. 45 Spectrum™ Theory Manual, Centric Engineering Systems, Inc., 1994. 46 Spectrum™ Visualizer Reference Manual, Centric Engineering Systems, Inc., 1994.