X. Y. Xu Research Fellow.
M. W. Collins Professor and Deputy Director. Thermo-Fluids Engineering Research Centre, City University, London EC1V OHB, U.K.
C. J. H. Jones Lecturer. Department of Cardiology, University of Wales College of Medicine, Heath Park, Cardiff CF4 4XN, U.K.
Flow Studies in Canine Artery Bifurcations Using a Numerical Simulation Method Three-dimensional flows through canine femoral bifurcation models were predicted under physiological flow conditions by solving numerically the time-dependent threedimensional Navier-stokes equations. In the calculations, two models were assumed for the blood, those of (a) a Newtonian fluid, and (b) a non-Newtonian fluid obeying the power law. The blood vessel wall was assumed to be rigid this being the only approximation to the prediction model. The numerical procedure utilized a finite volume approach on a finite element mesh to discretize the equations, and the code used (ASTEC) incorporated the SIMPLE velocity-pressure algorithm in performing the calculations. The predicted velocity prof lies were in good qualitative agreement with the in vivo measurements recently obtained by Jones et al. [1]. The non-Newtonian effects on the bifurcation flow field were also investigated, and no great differences in velocity profiles were observed. This indicated that the nonNewtonian characteristics of the blood might not be an important factor in determining the general flow patterns for these bifurcations, but could have local significance. Current work involves modeling wall distensibility in an empirically valid manner. Predictions accommodating these will permit a true quantitative comparison with experiment.
Introduction The development of atherosclerosis is determined by two pathogenetic factors: the arterial wall and the contents of the vessel. Fluid dynamics has for a long time been recognized as playing an important role in the development of the first pathogenetic concept, due to the regular localization of atherosclerotic lesions in the vicinity of arterial branchings. This has provided a primary motivation for the rise of interest in the study of blood flow through the major mammalian arteries, particularly at the bifurcation sites. Such flows are pulsatile and pass through vessels which are complex both in geometry and in their elastic nature. These factors are all important, not only in determining the detailed characteristics of arterial flows, but also the spatial nature of any fluid dynamics involvement in the disease process. Over the past two decades, considerable progress has been made in both experimental and numerical investigations of the flow near arterial bifurcations in models. Laser-Doppler anemometers have played an important role in obtaining quantitative measurement of velocities in in vitro model studies. Recent technological advances in Doppler ultrasound also provided a powerful tool for the measurement of detailed velocity patterns in vivo. In such a study carried out by Jones et al. [1], velocity profiles at surgically exposed arterial bifurcations Contributed by the Bioengineering Division for publication in the JOURNAL OF BIOMECHANICAL ENGINEERING. Manuscript received by the Bioengineering Division January 31, 1991; revised manuscript received February 28, 1992. Associate Technical Editor: T. J. Pedley.
in anaesthetized dogs were measured at 8.4 ms intervals in 7 different planes by using a 20 MHz 80 channel Doppler ultrasound device. These results probably represent the most detailed measurements of the flow field at intact bifurcations in vivo so far obtained. An alternative approach to gain more knowledge of flow patterns and wall shear stresses at bifurcations is to make use of the newly developed computational fluid dynamics tools and supercomputers, which allow the calculation of threedimensional pulsatile flow in a bifurcation segment. Two-dimensional numerical studies have been extensively carried out (e.g., references [2-5]), while three-dimensional simulations have just started in the last few years. As an initial step, Wille [6] developed a three-dimensional mathematical model of steady flow in a symmetrical bifurcation by using a finite element method. Nevertheless, this method required too long a computation time, thus prohibiting further simulation development. Another computational scheme based on finite difference methods was presented by Dinnar et al. [7]. In an application to a 90 deg T-bifurcation of rectangular crosssection, the suggested scheme proved to be extremely efficient when compared with a finite element approach. The latest three-dimensional predictions have been performed by Rindt [8] for steady flow in a carotid artery bifurcation, Perktold and Peter [9] for pulsatile flow in an arterial T-bifurcation model, and Perktold and Resch [10] for pulsatile flow in human carotid artery bifurcations, by using penalty function and newly developed pressure correction finite element computational
504/Vol. 114, NOVEMBER 1992
Transactions of the ASME
Copyright © 1992 by ASME Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 05/22/2015 Terms of Use: http://asme.org/terms
MODEL ( A )
riODEL
Fig. 1 Schematic illustration and three-dimensional representation of two different canine ilio-femoral bifurcation models
methods respectively with supercomputers.' In all these studies, the blood was assumed to be homogeneous, incompressible and Newtonian, the vessel wall being rigid. The present study attempts to provide comparative predictive data for the above in vivo canine femoral bifurcation results by means of a computer simulation approach. In this study of two of the tested bifurcations three-dimensional flows through these bifurcations were calculated under physiological flow conditions, the non-Newtonian viscosity of blood also being accommodated. For these results the only effect not being treated is the distensibility of the wall. All calculations presented were performed using ASTEC (developed at UKAEA), which has a finite element definition and finite volume solution methods, on a supercomputer CRAY X-MP/28.
a
0.0
MODEL ( A )
Geometrical Models and Grid Generation It is a well accepted fact that there are rather large individual variations in bifurcation geometry. The two bifurcation models adopted in the present study were based on data from the in vivo measurements [1], in which branching angles, B-mode ultrasound images of the cross-sections of the upstream parent and both downstream daughter vessels, and the origin of both vessels were all photographed. The bifurcation geometrical data were measured from photographs. Figure 1 gives a schematic illustration and three-dimensional representation of two different canine ilio-femoral bifurcation models. In both cases, parent and daughter vessels were straight, with the bifurcation lying in one plane which is the plane of symmetry. These assumptions are consistent with the measurements. The two models have similar configurations and differ only in vessel diameters and parent to large daughter angle. MODEL Figure 1 also illustrates the sites at which the comparison Fig. 2 Measured blood velocity waveforms and the time-dependent between predictions and measurements is presented. In terms velocity profiles at site 1 in the canine femoral arteries of the parent vessel diameter, site 1 is two diameters upstream of the flow divider in the parent artery, site 2 is at the level of flow divider in the larger daughter vessel, sites 3 and 4 are one and four diameters downstream respectively, and site 5 is one diameter downstream of the flow divider in the smaller daughter vessel. For convenience of comparison, they are placed at 60° to the vessel axes, since the measurements were taken ' Of course, despite the added computational costs, for a true physiological geometry, a finite element approach has been virtually essential to date. Now, at such positions. however, the latest possibility to emerge involves multi-block finite difference For the application to modeling of an arterial bifurcation methods. In our future work we intend to use this approach with an alternative the grid generation problem is crucial, since the complex code from Harwell Laboratory UK, FLOW3D, Release 3, late 1991. (see refboundary shape has to be presented accurately to ensure that erence [19]). Journal of Biomechanical Engineering
NOVEMBER 1992, Vol. 114 / 505
Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 05/22/2015 Terms of Use: http://asme.org/terms
no unnecessary errors are introduced. Using the semiautomatic mesh generator SOPHIA [11] available with ASTEC, a purpose built routine for handling intersections of vessels was developed for general bifurcation studies. The rationale for this routine was to be able to generate any given bifurcation, with a minimum amount of input data. By using the routine, a generalized three-dimensional arterial bifurcation geometry can be gridded into a mesh consisting of arbitrarily shaped 8node blocks, on which finite volume or finite element based fluid flow solvers can easily be applied.
account at the solid wall, i.e., each velocity component was set equal to zero. At the bifurcation inlet, the time-dependent axial velocity profiles obtained from the in vivo measurements were directly imposed. As shown in Fig. 2, these velocity profiles are relatively symmetric, although in some phases they show a tendency to be skewed towards the larger daughter vessel side.2 It should also be noted that the measured inlet profiles are only axial, and we have assumed negligible secondary velocities, which is an approximation or small error, possibly of most• significance in reverse flow. This issue is discussed later in the paper when predictions of the secondary flow field are presented. These will clarify the quantitative Numerical Models significance of the above assumption. 1 Basic Equations and Boundary Conditions. The probFor the downstream boundary conditions, there are genlem considered herein is a pulsatile flow of a homogeneous, erally three types of treatment at the outlet: (1) to assume fully incompressible, non-Newtonian fluid through a three-dimen- developed flow; (2) to assume a zero condition of surface sional model with a rigid wall. Such a problem can be described traction; and (3) to specify a constant pressure. The first type by the time-dependent three-dimensional Navier-Stokes equa- is more restrictive as it requires the flow rate ratio between the tions. In integral form, the basic momentum and continuity two branches to be known a priori, also the length of the principles for blood of constant density may be written as: specified branch has to be large enough to permit development. The other two types are more realistic for this application, the p^-\udV= - p ( u - u - r f A - \ VPdV+ U e ( V u ) - d A (1) constant pressure outflow condition (3) being used here. For the two bifurcations investigated, zero pressures were specified at both outlets which were placed at five diameters (in terms of parent vessel diameter) downstream from the flow divider. u-rfA = 0 (2) Outflow condition (3) resulted in the flow division ratio changing throughout the cycle. The fraction of parent flow rate where u = («, v, w)Tis a velocity vector; dA and dVrepresent passing through the larger daughter vessel varied from 0.61 to elements of control area and control volume respectively; and 0.67 for model A, and 0.52 to 0.60 for model B. P is pressure. The density of blood is p, and the blood viscosity The calculations were carried out under the corresponding is ne, which is a function of strain rate in the non-Newtonian pulsatile flow conditions for individual canine models. In this case. It is here derived from a power law expression: study the average Reynolds number is defined as n (3) i?e = 4V/7rD, where V is the average flow rate in the vessel T=-m\y\ over a whole cycle, D is the diameter of the vessel, and v the hence the non-Newtonian viscosity \ie can be written as kinematic viscosity of the blood (4x 10~6 m 2 /s). The Womn (4) ersley parameter is defined as a = D/2 sjoi/v, where oi = 2ir / i>,e = v.e{\y\) = m\y\ -' where T and y denote shear stress and shear strain, m and n is the angular frequency. In the flow simulation, the average being constants. For 3-D flows, Eq. (4) can be modified using inflow Reynolds numbers were 92 and 108, and the Womersley parameters were 2.04 and 2.39, for bifurcation models A and the second invariant I2 of the strain rate tensor: B, respectively. These values are in the lower physiological (5) range, and therefore turbulent effects are unlikely.
i^SS^' I
j
where djj (ij = 1 , 2 , 3 ) are the components of strain rate defined as:
(6)
^= 2fe+ ^j 1 fdUi
du,\
the modified Eq. (4) is: lxe
= m\\IW2\n-x
(7)
where
( Yy Tz) dv\
dw
+
(du +
dw\ +
\Tz YX)
(8)
These equations are discretized by employing the finite volume method. The procedure for solving the discrete equations is iterative and based upon the SIMPLE solution algorithm [12]. The pressure corrections are calculated using a preconditioned conjugate gradient algorithm. A vector upwinding scheme is employed, which greatly reduces the false diffusion errors associated with the advection terms. The time differencing is fully implicit, so avoiding Courant stability restriction on the time-step size. The discretized Eqs. (1), (2), and (7) are solved together under appropriate boundary conditions. The no-slip boundary conditions and the rigidness of the bifurcation were taken into
2 Newtonian and Non-Newtonian Models. It has been believed for a long time that blood behaves essentially as a Newtonian fluid during flows through large blood vessels. However, in certain experimental studies [13-15] there is evidence which indicates some non-Newtonian effects even in large arteries. In view of this it is desirable to accommodate some non-Newtonian models of possible blood behavior. Several theories have been postulated to describe such behavior of the blood, for instance, the Casson and Bingham models have found a wide application in the past [16-17]. However, there are still doubts regarding the existence of the yield stress in dynamic situations since it has been observed only under static loading conditions. Therefore, a more appropriate form of the constitutive equation would be a general power law. As shown in Eq. (3), the model utilized in this study is a relatively simple one which depends upon shear rate only. For human blood the constants m and n are taken to be 0.042 Pa-s and 0.61, respectively [13]. A modified power law expression for whole blood which depends upon shear rate, hematocrit, and total protein minus albumin [16] is expected to be adopted in our future study. 3 Wall Distensibility.
Arterial walls are viscoelastic in-
2 An alternative upstream boundary condition comprised a typical average canine velocity waveform supplied us by Dr. Kim Parker, and was allowed to develop in a downstream direction over the relatively straight parent tube upstream of site 1. Although predictions are omitted here for reason of space, they compared substantially less well than those of Fig. 6. Hence it is necessary to use a specific waveform in making detailed bifurcation predictions.
506 / Vol. 114, NOVEMBER 1992 Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 05/22/2015 Terms of Use: http://asme.org/terms
Transactions of the ASME
I10DEL
MODEL
(A)
Fig. 4 Three-dimensional perspective view of axial velocity (vertical axis) versus radial position (r) and phase (time, f) at five different sites (as defined in Fig. 1) in the bifurcation model B
MODEL
CA3
Fig. 3(6)
jO.Bm/s
(Ap)
(Cp) MODEL
(B)
Fig. 3(c) Fig. 3 Axial velocity profiles of a Newtonian fluid in the bifurcation plane at four different phases Ap, Bp, Cp, Dp during a cardiac cycle (see text for explanation)
homogeneous multilayered tissues. They are composed mainly of collagen, elastin fibers and smooth muscle cells. In most of the previous studies, the arterial wall was assumed to be rigid, due to the general belief that its distensibility may be regarded as a "second-order" effect on arterial bifurcation flows and thus negligible. However, some recent experimental studies [13, 15] have revealed the differences between rigid and elastic wall behavior during pulsatile flow, and demonstrated the fact that wall distensibility may play an important role in the understanding of bifurcation flows. The difficulties encountered in the incorporation of the moving wall are twofold: one is to Journal of Biomechanical Engineering
choose a simple and accurate mathematical expression of the mechanical behavior of the arterial wall, i.e., a constitutive equation, and another is to construct an interface between the moving wall and the fluid flow. The former has been partially solved by utilizing the pseudo strain energy functions which were concluded as the most practical approximation of the stress-strain relationship for the arterial wall [18]. The latter requires moving boundary numerical techniques, such as adaptive gridding (see reference [19]). Research is going on in our group to enable the wall behavior to be modeled in a comprehensive manner; however, in the present study a rigid wall is used. 4 Grid Resolution Tests. For the two cases presented, it was assumed that a symmetry plane existed so that simulations were confined to only half of the bifurcation. The computational mesh consists of 5040 elements and 6090 nodes, with 18 elements in the radial direction at each cross-section. The grid chosen is a compromise between accuracy and available computer time and storage. Preliminary calculations were carried out for bifurcation model B under steady flow conditions using three grid distributions, i.e., 12, 18, and 24 elements in the radial direction. Results of these calculations indicated that a grid with 18 elements in the radial direction was sufficient for the present purpose. 3 The average velocity differences were 0.19 percent with the finer grid and 6.17 percent with the coarser grid. With the finer grid (24-18 elements) the maximum axial velocity differed by 1.36 percent and the maximum wall shear stress differed by 2.5 percent. Before application to the current bifurcation study, the computer code ASTEC was tested by making a prediction for fully developed laminar flow in a straight tube with both Newtonian and non-Newtonian fluids. The computer velocity profiles were in good agreement with the theoretical solutions (maximum local velocity error 0.3 percent). The code has also been previously validated for various bifurcation flow cases, i.e., two-dimensional steady and pulsatile flows in a T-junction, three-dimensional steady and pulsatile flows in a T-junction with circular cross-section. The comparison between the computer results and laboratory measurements were very satisfactory and typical presentations are in reference [20]. Results and Discussions Initial flow predictions were performed under Newtonian assumptions for the blood. Figures 3(6) and (c) show the axial velocity profiles in the bifurcation plane at the pulse 3 This grid convergence check confirmed the results of an earlier exercise on a two-dimensional T-junction bifurcation geometry.
NOVEMBER 1992, Vol. 114 / 507
Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 05/22/2015 Terms of Use: http://asme.org/terms
INNER WALL
MODEL B
Fig. 5 Wall shear stress versus time at (a) one diameter above the flow divider in the parent vessel, and (b) 0.2 diameter downstream of the flow divider In the larger daughter vessel on the outer wall, and the inner wall for bifurcation model B
phases of accelerated flow Ap, maximum positive flow Bp, decelerated flow Cp, and maximum reverse flow Dp (as defined in Fig. 3(a)) for the two bifurcation models. Comparison of Fig. 3 (b) and (c) for bifurcation models A and B respectively, shows that the essential features of the flow are the same. Due to space limitation, we present results in the remainder of the paper for model B only. Figure 4 provides a three-dimensional view of the axial velocity at five different sites in the plane of the bifurcation for model B. From these results, it is observed that in the parent' vessel velocity profiles are hardly affected by the flow branching downstream. They are relatively symmetrical with reverse flow occurring during later diastole. In the two daughter vessels, velocity profiles vary with position and phase during the cycle. At the level of the flow divider, velocity profiles are skewed towards the inner wall (at the flow divider side) due to the presence of the apex and the axial momentum of the fluid when it leaves the parent vessel. As the flow moves downstream, the peak axial velocity gradually shifts to the center of the vessel. Near the outer wall of the larger daughter vessel, flow separation does not occur during flow acceleration but is present during later flow deceleration just before flow reversal. This suggests that in bifurcation regions separated flow might not always exist throughout the cycle, its presence possibly depending upon some physical parameters which characterize the flow, such as Reynolds number and flow boundary conditions. This finding is consistent with previous experimental observations [21-22]. Under the current flow conditions, only small differences in velocity profiles at two bifurcations with different branching angles are noticed. Wall shear stress was also calculated for bifurcation model
•o^~4000.9.~
feS^ 4000000000005^
sftoooooooooooooso*^
We*U2!°l
• -»««*"^
-•-a-
•"*»
^onWOnonnrtaa
»'«-«'O"0^«'«'«"«
|c)
(a) AT SITE 1 OF DOG B COMPUTATIONAL RESULTS oa°°«
EXPERIMENTAL RESULTS
508 / Vol. 114, NOVEMBER 1992
(b) at s i t e 2 of dog
AT S I T E 3 OF DOG 1
(d)
AT S I T E
4 OF DOG 1
(e)
AT S I T E
5 OF DOG 1
Fig. 6 Comparison of axial velocity profiles between predictions (—) using a non-Newtonian model and the in vivo measurements (© © o) at five different sites (as defined in Fig. 1) in bifurcation model B. In each case, the right hand side is the outer wall, and the left-hand side is the inner wall (flow divider side).
Transactions of the ASME
Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 05/22/2015 Terms of Use: http://asme.org/terms
1.0 0.2 0.1 0.5 0 (a) AT SITE 1 OF DOG B (SOLID LINE: NEWTONIAN MODEL) (DIAMONDS: NON-NEWTON IAN MODEL)
AT SITE 2 OF DOG I
B when a Newtonian fluid was used. The wall shear stress obtained at one diameter above the flow divider in the parent vessel is presented in Fig. 5(c) as a function of time during the cycle. The peak stress is 2.36 N/m2, occurring during the acceleration phase of the systole, and the time-averaged mean value is 0.28 N/m2. At 0.2 diameter downstream of the flow divider in the large daughter vessel, the wall shear stress behavior is considerably different as shown in Fig. 5(b). On the outer wall the magnitudes of the stress are rather low and the time-averaged mean value is close to zero. On the inner wall the shear stress is much greater reaching the maximum and time-averaged mean values of 5.9 N/m2 and 1.10 N/m2, respectively. The latter is far above that found in the parent vessel. To verify the computational results, comparison between calculations and measurements is necessary. In Fig. 6, predicted axial velocity profiles for the bifurcation model B with non-Newtonian fluid are compared with those of the in vivo measurements at five different sites shown in Fig. 1. The predictions represent our latest model. However, there are several facts which must be taken into account in making the comparison. First of all, the measurements were performed by using the Doppler ultrasound device which measures the velocity component in the direction of beam but calculates the velocity as though it is parallel to the vessel axis. This assumption is likely to give misleading results in such complex flows. Secondly, the arterial wall was treated as rigid in the predictions due to difficulties in coding a moving boundary problem at the current stage. According to the experimental studies published so far, the wall distensibility might have a significant, even if local, effect on arterial bifurcation flows, especially when it is coupled with the non-Newtonian character of the blood. Thirdly, there are uncertainties about the velocity very close to the wall and the position of the moving wall itself in the measurements. Therefore, the comparison between the predictions and the measurements has to be made very carefully, and at this stage it can only be done with diffidence. In Fig. 6, it is observed that during diastole there is a very good agreement between the calculations and the measurements in parent and both daughter vessels at all sites. Good agreement is also obvious in the parent vessel at site 1, and in the larger daughter vessel at site 4, where the velocity profiles are rather symmetric. Some differences are noticed in the daughter vessels at sites 2, 3, and 5, which are near the flow divider, during Journal of Biomechanical Engineering
0.6 0.8 1.0 AT SITE 3 OF DOG B
0 Id) AT SITE 4 OF DOG B (SOLID LINE: NEWTONIAN MODEL) (DIAMONDS: NON-NEWTONIAN MODEL)
0.2
0.1
0.6
0.8
t .0 (x/d)
(e) AT SITE 5 OF DOG !
Fig. 7 Comparison of axial velocity profiles of a Newtonian (—) and a non-Newtonian fluid (o o o) for bifurcation model B at five different sites. In each case, the right-hand side is the outer wall, and the left hand side is the inner wall.
the acceleration phase of systole. The predicted velocity profiles show greater skewing towards the inner wall. We feel this is mainly due to the assumption of parallel flow made in the measurements and the rigidness of the vessel wall adopted in the calculations. In order to investigate quantitatively the non-Newtonian effects on the bifurcation flow field, calculated axial velocity profiles of a Newtonian fluid are compared with those of a non-Newtonian fluid for bifurcation model B at five different sites. As shown in Fig. 7, in general, the velocity profiles appear rather flattened in the non-Newtonian fluid case. This flow NOVEMBER 1992, Vol. 114 / 509
Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 05/22/2015 Terms of Use: http://asme.org/terms
MODEL
(B)
(Bp)
A2
(Ap)
Fig. 8
B
(Cp)
Secondary velocity vector fields for bifurcation model B
behavior has already been demonstrated by others [14]. At the bifurcation region, no large variations in the velocity distribution are observed, which indicates that the non-Newtonian characteristics of the blood may not play an important role in determining the general flow field at these bifurcations for the current flow conditions. However, it is still possible that these characteristics may significantly affect a flow situation which is sensitive to small changes. Such a postulation is in accordance with previous experimental studies [13-15], In these the velocity distributions behind the bifurcation are substantially different in Newtonian and non-Newtonian fluids. In general, then, the quantitative effects of non-Newtonian viscosity appear to depend very much upon the bifurcation geometry and the flow conditions. It is, therefore, valuable to have the option of a non-Newtonian model in the code. Since the in vivo data essentially gives axial flow information, we have concentrated so far on this in the paper. However, the predictions provided a full three-dimensional flow field. In Fig. 8 are presented typical secondary velocity vector fields for bifurcation model B at three different phases Ap, Bp, and Cp as defined. These velocity vectors represent therefore the remainder of the flow field information, and give an indication of any swirling type structures. The results are qualitatively consistent with other published data for bifurcations [8, 20]. It has previously been mentioned that the upstream boundary condition assumed a zero secondary velocity field, due to the lack of measuring this in experiment. Such an approximation is analogous to the accepted procedure for a numerical experiment to test equation stability. For example, in Collins [23] Appendix, even a substantial deliberate error in axial velocity was shown rapidly to diffuse out to a negligible effect, thus demonstrating good stability in the algorithm. Here, the secondary velocities at the near inlet position A1-A2 show that they are relatively small during positive flow. Although during reverse flow they are slightly distorted, these errors are shown
to diminish downstream. However, it would be preferable to include measurements of secondary as well as axial velocities for an upstream boundary condition in future studies. Conclusions To provide predictive data for in vivo canine femoral bifurcation measurements, a numerical investigation was carried out for the tested bifurcations under physiological flow conditions by means of a three-dimensional fluid flow code and a supercomputer. The results obtained for the rigid (the only model approximation) canine femoral bifurcation models with (a) a Newtonian, and (b) a non-Newtonian fluid, are summarized as follows: 1) Under the current flow conditions, i.e., average inflow Reynolds number of 92 and 108, Womersley parameter of 2.04 and 2.39 for bifurcation models A and B, respectively, flow separation was not obvious and existed only for a small period of the whole cycle. 2) Shear stress on the outer wall of the larger daughter vessel was very low and the time-averaged mean value was close to zero. Shear stress on the inner wall of the larger daughter vessel was much greater than either on the outer wall or in the parent vessel. 3) While there was good qualitative and not unreasonable quantitative agreement between the predicted velocity profiles and those of the in vivo measurements, they differed somewhat in detail in the daughter vessels in the branching region, indicating the importance of the secondary flow influences on the flow patterns there. 4) For the two bifurcation models and the current flow conditions, the non-Newtonian characteristic of the blood was not found to be an important factor in determining the general flow field parameters. The numerical results obtained are believed to be important
510 / Vol. 114, NOVEMBER 1992 Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 05/22/2015 Terms of Use: http://asme.org/terms
Transactions of the ASME
in interpreting the measured velocity profiles and acquiring more information on the flow field in bifurcation regions. The present study represents an initial attempt to provide comparable results to in vivo measurements of the flow in arterial bifurcations. Work is continuing to enable wall distensibility effects to be studied. When these are incorporated into the predictions, a detailed quantitative assessment will be made of both experiment and theory.
Acknowledgments The authors are very grateful to Dr. K. H . Parker for his help and discussion, and especially for providing detailed experimental data. This project was supported by a British Heart Foundation Research Grant. Finally, ASTEC is a proprietary code of the United Kingdom Atomic Energy Authority. References 1 Jones, C. J. H., Lever, M. J., Ogasawara, Y., Parker, K. H., Hiramatsu, O., Mito,K.,Tsujioka, K.,andKajiya, F., "Blood Velocity Distributions Within Intact Canine Arterial Bifurcations," American Journal of Physiology: (Heart and Circulatory Physiology, Vol. 31, 1992, pp. H1592-1599. 2 O'Brien, V., and Ehrilich, L. M., "Simulation of Unsteady Flow at Renal Branches," Journal of Biomechanics, Vol. 10, 1977, pp. 623-631. 3 Liepsch, D., and Moravec, S.,' 'Measurement and Calculations of Laminar Flow in a Ninety Degree Bifurcation," Journal of Biomechanics, Vol. 15, 1982, pp. 473-485. 4 Perktold, K., and Hilbert, D., "Numerical Simulation of Pulsatile Flow in a Carotid Bifurcation Model," Journal of Biomechanics, Vol. 8, 1986, pp. 193-199. 5 Rindt, C. C. M., Vosse, F. N. N. V. D., Steenhoven, A. A. V., and Janssen, J. D., "A Numerical and Experimental Analysis of the Flow in a TwoDimensional Model of the Human Carotid Bifurcation," Journal of Biomechanics, Vol. 20, 1987, pp. 499-509. 6 Wille, S., "Numerical Simulation of Steady Flow Inside a Three-Dimensional Aortic Bifurcation Model," ASME JOURNAL OF BIOMECHANICAI ENGINEERING, Vol. 6, 1986, pp. 49-55. 7 Dinnar, U., Enden, G., and Israeli, M., " A Numerical Study of Flow in a Three-Dimensional Bifurcation," Cardiovascular Systems Dynamics Society Meeting, 1988, Canada. 8 Rindt, C. C. M., "Analysis of the Three-Dimensional Flow Field in the
Carotid Artery Bifurcation," PhD Thesis, 1989, Eindhoven University of Technology, The Netherlands. 9 Perktold, K., and Peter, R., "Numerical 3D-Simulation of Pulsatile Wall Shear Stress in an Arterial T-Bifurcation Model," ASME JOURNAL OF BIOMECHANICAL ENGINEERING, Vol. 12, 1990, pp. 2-13.
10 Perktold, K., and Resch, M. "Numerical Flow Studies in Human Carotid Artery Bifurcations: Basic Discussion of the Geometric Factor in Atherogenesis," ASME JOURNAL OF BIOMECHANICAL ENGINEERING, Vol. 12, 1990, pp. 111-
123. 11 Stainsby, R., " A Guide to Using the SOPHIA Mesh Generator," UKAEA (Draft), 1989. 12 Patankar, S. V.] and Spalding, D. B., "A Calculation Procedure for Heat, Mass and Momentum Transfer in Three-Dimensional Parabolic Flows," International Journal of Heat and Mass Transfer, Vol. 15, 1972, pp. 1787-1800. 13 Ku, D. N., and Liepsch, D., "The Effects of Non-Newtonian Viscoelasticity and Wall Elasticity on Flow at a 90° Bifurcation," Biorheology, Vol. 23, 1986, pp. 359-370. 14 Moravec, S., and Liepsch, D., "Flow Investigations in a Model of a ThreeDimensional Human Artery with Newtonian and Non-Newtonian Fluids, Part I , " Biorheology, Vol. 20, 1983, pp. 745-759. 15 Liepsch, D. W., "Flow in Tubes and Arteries—A Comparison," Biorheology, Vol. 23, 1986, pp. 395-433. 16 Walburn, F. J., and Schneck, D. J., " A Constitutive Equation for Whole Human Blood," Biorheology, Vol. 13, 1973, pp. 201-210. 17 Nakamura, M., and Sawada, T., "Numerical Study on the Flow of a NonNewtonian Fluid Through an Axisymmetric Stenosis," ASME JOURNAL OF BIOMECHANICAL ENGINEERING, Vol. 110, 1988, pp. 137-143.
18 Fung, Y. C , Pronek, K., and Patitucci, P., "Pseudoelasticity of Arteries and the Choice of its Mathematical Expression," American Journal of Physiology, Vol. 237, 1979, pp. H620-H631. 19 Xu, X. Y., and Collins, M. W., " A Review of the Numerical Analysis of Blood Flow in Arterial Bifurcations," Proc. Instn. Mech. Engrs., Vol. 204, Part H, Journal of Engineering in Medicine, 1990, pp. 205-216. 20 Collins, M. W., and Xu, X. Y., "A Predictive Scheme for Flow in Arterial Bifurcations: Comparison with Laboratory Measurements," Biomechanics Transport Processes, ed., F. Mosora, C. Caro et al., Plenum Press, New York and London, 1990, pp. 125-133. 21 Richardson, P. D., and Christo, J. L., "Flow Separation Opposite a Side Branch," Biofluid Mechanics, Blood Flow in Large Vessels, ed., D. Liepsch, Springer-Verlag, 1990, pp. 275-284. 22 Cho, Y. I., Back, L. H., and Crawford, D. W., "Experimental Investigation of Branch Flow Ratio, Angle, and Reynolds Number Effects on the Pressure, and Flow Fields in Arterial Branch Models," ASME JOURNAL OF BIOMECHANICAL ENGINEERING, Vol. 107, 1985, pp. 257-267.
23 Collins, M. W., "Finite Difference Analysis for Developing Laminar Flow in Circular Tubes Applied to Forced and Combined Convection," International Journal for Numerical Methods in Engineering, Vol. 15, 1980, pp. 381-404.
(Contents continued) Announcements and Special Notices 460 Transactions Change of Address Form 496
Errata on a Previously Published Paper
553 Calendar of Events
Journal of Biomechanical Engineering
NOVEMBER 1992, Vol. 114 / 511
Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 05/22/2015 Terms of Use: http://asme.org/terms
M. W. Collins Professor and Deputy Director. Thermo-Fluids Engineering Research Centre, City University, London EC1V OHB, U.K.
C. J. H. Jones Lecturer. Department of Cardiology, University of Wales College of Medicine, Heath Park, Cardiff CF4 4XN, U.K.
Flow Studies in Canine Artery Bifurcations Using a Numerical Simulation Method Three-dimensional flows through canine femoral bifurcation models were predicted under physiological flow conditions by solving numerically the time-dependent threedimensional Navier-stokes equations. In the calculations, two models were assumed for the blood, those of (a) a Newtonian fluid, and (b) a non-Newtonian fluid obeying the power law. The blood vessel wall was assumed to be rigid this being the only approximation to the prediction model. The numerical procedure utilized a finite volume approach on a finite element mesh to discretize the equations, and the code used (ASTEC) incorporated the SIMPLE velocity-pressure algorithm in performing the calculations. The predicted velocity prof lies were in good qualitative agreement with the in vivo measurements recently obtained by Jones et al. [1]. The non-Newtonian effects on the bifurcation flow field were also investigated, and no great differences in velocity profiles were observed. This indicated that the nonNewtonian characteristics of the blood might not be an important factor in determining the general flow patterns for these bifurcations, but could have local significance. Current work involves modeling wall distensibility in an empirically valid manner. Predictions accommodating these will permit a true quantitative comparison with experiment.
Introduction The development of atherosclerosis is determined by two pathogenetic factors: the arterial wall and the contents of the vessel. Fluid dynamics has for a long time been recognized as playing an important role in the development of the first pathogenetic concept, due to the regular localization of atherosclerotic lesions in the vicinity of arterial branchings. This has provided a primary motivation for the rise of interest in the study of blood flow through the major mammalian arteries, particularly at the bifurcation sites. Such flows are pulsatile and pass through vessels which are complex both in geometry and in their elastic nature. These factors are all important, not only in determining the detailed characteristics of arterial flows, but also the spatial nature of any fluid dynamics involvement in the disease process. Over the past two decades, considerable progress has been made in both experimental and numerical investigations of the flow near arterial bifurcations in models. Laser-Doppler anemometers have played an important role in obtaining quantitative measurement of velocities in in vitro model studies. Recent technological advances in Doppler ultrasound also provided a powerful tool for the measurement of detailed velocity patterns in vivo. In such a study carried out by Jones et al. [1], velocity profiles at surgically exposed arterial bifurcations Contributed by the Bioengineering Division for publication in the JOURNAL OF BIOMECHANICAL ENGINEERING. Manuscript received by the Bioengineering Division January 31, 1991; revised manuscript received February 28, 1992. Associate Technical Editor: T. J. Pedley.
in anaesthetized dogs were measured at 8.4 ms intervals in 7 different planes by using a 20 MHz 80 channel Doppler ultrasound device. These results probably represent the most detailed measurements of the flow field at intact bifurcations in vivo so far obtained. An alternative approach to gain more knowledge of flow patterns and wall shear stresses at bifurcations is to make use of the newly developed computational fluid dynamics tools and supercomputers, which allow the calculation of threedimensional pulsatile flow in a bifurcation segment. Two-dimensional numerical studies have been extensively carried out (e.g., references [2-5]), while three-dimensional simulations have just started in the last few years. As an initial step, Wille [6] developed a three-dimensional mathematical model of steady flow in a symmetrical bifurcation by using a finite element method. Nevertheless, this method required too long a computation time, thus prohibiting further simulation development. Another computational scheme based on finite difference methods was presented by Dinnar et al. [7]. In an application to a 90 deg T-bifurcation of rectangular crosssection, the suggested scheme proved to be extremely efficient when compared with a finite element approach. The latest three-dimensional predictions have been performed by Rindt [8] for steady flow in a carotid artery bifurcation, Perktold and Peter [9] for pulsatile flow in an arterial T-bifurcation model, and Perktold and Resch [10] for pulsatile flow in human carotid artery bifurcations, by using penalty function and newly developed pressure correction finite element computational
504/Vol. 114, NOVEMBER 1992
Transactions of the ASME
Copyright © 1992 by ASME Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 05/22/2015 Terms of Use: http://asme.org/terms
MODEL ( A )
riODEL
Fig. 1 Schematic illustration and three-dimensional representation of two different canine ilio-femoral bifurcation models
methods respectively with supercomputers.' In all these studies, the blood was assumed to be homogeneous, incompressible and Newtonian, the vessel wall being rigid. The present study attempts to provide comparative predictive data for the above in vivo canine femoral bifurcation results by means of a computer simulation approach. In this study of two of the tested bifurcations three-dimensional flows through these bifurcations were calculated under physiological flow conditions, the non-Newtonian viscosity of blood also being accommodated. For these results the only effect not being treated is the distensibility of the wall. All calculations presented were performed using ASTEC (developed at UKAEA), which has a finite element definition and finite volume solution methods, on a supercomputer CRAY X-MP/28.
a
0.0
MODEL ( A )
Geometrical Models and Grid Generation It is a well accepted fact that there are rather large individual variations in bifurcation geometry. The two bifurcation models adopted in the present study were based on data from the in vivo measurements [1], in which branching angles, B-mode ultrasound images of the cross-sections of the upstream parent and both downstream daughter vessels, and the origin of both vessels were all photographed. The bifurcation geometrical data were measured from photographs. Figure 1 gives a schematic illustration and three-dimensional representation of two different canine ilio-femoral bifurcation models. In both cases, parent and daughter vessels were straight, with the bifurcation lying in one plane which is the plane of symmetry. These assumptions are consistent with the measurements. The two models have similar configurations and differ only in vessel diameters and parent to large daughter angle. MODEL Figure 1 also illustrates the sites at which the comparison Fig. 2 Measured blood velocity waveforms and the time-dependent between predictions and measurements is presented. In terms velocity profiles at site 1 in the canine femoral arteries of the parent vessel diameter, site 1 is two diameters upstream of the flow divider in the parent artery, site 2 is at the level of flow divider in the larger daughter vessel, sites 3 and 4 are one and four diameters downstream respectively, and site 5 is one diameter downstream of the flow divider in the smaller daughter vessel. For convenience of comparison, they are placed at 60° to the vessel axes, since the measurements were taken ' Of course, despite the added computational costs, for a true physiological geometry, a finite element approach has been virtually essential to date. Now, at such positions. however, the latest possibility to emerge involves multi-block finite difference For the application to modeling of an arterial bifurcation methods. In our future work we intend to use this approach with an alternative the grid generation problem is crucial, since the complex code from Harwell Laboratory UK, FLOW3D, Release 3, late 1991. (see refboundary shape has to be presented accurately to ensure that erence [19]). Journal of Biomechanical Engineering
NOVEMBER 1992, Vol. 114 / 505
Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 05/22/2015 Terms of Use: http://asme.org/terms
no unnecessary errors are introduced. Using the semiautomatic mesh generator SOPHIA [11] available with ASTEC, a purpose built routine for handling intersections of vessels was developed for general bifurcation studies. The rationale for this routine was to be able to generate any given bifurcation, with a minimum amount of input data. By using the routine, a generalized three-dimensional arterial bifurcation geometry can be gridded into a mesh consisting of arbitrarily shaped 8node blocks, on which finite volume or finite element based fluid flow solvers can easily be applied.
account at the solid wall, i.e., each velocity component was set equal to zero. At the bifurcation inlet, the time-dependent axial velocity profiles obtained from the in vivo measurements were directly imposed. As shown in Fig. 2, these velocity profiles are relatively symmetric, although in some phases they show a tendency to be skewed towards the larger daughter vessel side.2 It should also be noted that the measured inlet profiles are only axial, and we have assumed negligible secondary velocities, which is an approximation or small error, possibly of most• significance in reverse flow. This issue is discussed later in the paper when predictions of the secondary flow field are presented. These will clarify the quantitative Numerical Models significance of the above assumption. 1 Basic Equations and Boundary Conditions. The probFor the downstream boundary conditions, there are genlem considered herein is a pulsatile flow of a homogeneous, erally three types of treatment at the outlet: (1) to assume fully incompressible, non-Newtonian fluid through a three-dimen- developed flow; (2) to assume a zero condition of surface sional model with a rigid wall. Such a problem can be described traction; and (3) to specify a constant pressure. The first type by the time-dependent three-dimensional Navier-Stokes equa- is more restrictive as it requires the flow rate ratio between the tions. In integral form, the basic momentum and continuity two branches to be known a priori, also the length of the principles for blood of constant density may be written as: specified branch has to be large enough to permit development. The other two types are more realistic for this application, the p^-\udV= - p ( u - u - r f A - \ VPdV+ U e ( V u ) - d A (1) constant pressure outflow condition (3) being used here. For the two bifurcations investigated, zero pressures were specified at both outlets which were placed at five diameters (in terms of parent vessel diameter) downstream from the flow divider. u-rfA = 0 (2) Outflow condition (3) resulted in the flow division ratio changing throughout the cycle. The fraction of parent flow rate where u = («, v, w)Tis a velocity vector; dA and dVrepresent passing through the larger daughter vessel varied from 0.61 to elements of control area and control volume respectively; and 0.67 for model A, and 0.52 to 0.60 for model B. P is pressure. The density of blood is p, and the blood viscosity The calculations were carried out under the corresponding is ne, which is a function of strain rate in the non-Newtonian pulsatile flow conditions for individual canine models. In this case. It is here derived from a power law expression: study the average Reynolds number is defined as n (3) i?e = 4V/7rD, where V is the average flow rate in the vessel T=-m\y\ over a whole cycle, D is the diameter of the vessel, and v the hence the non-Newtonian viscosity \ie can be written as kinematic viscosity of the blood (4x 10~6 m 2 /s). The Womn (4) ersley parameter is defined as a = D/2 sjoi/v, where oi = 2ir / i>,e = v.e{\y\) = m\y\ -' where T and y denote shear stress and shear strain, m and n is the angular frequency. In the flow simulation, the average being constants. For 3-D flows, Eq. (4) can be modified using inflow Reynolds numbers were 92 and 108, and the Womersley parameters were 2.04 and 2.39, for bifurcation models A and the second invariant I2 of the strain rate tensor: B, respectively. These values are in the lower physiological (5) range, and therefore turbulent effects are unlikely.
i^SS^' I
j
where djj (ij = 1 , 2 , 3 ) are the components of strain rate defined as:
(6)
^= 2fe+ ^j 1 fdUi
du,\
the modified Eq. (4) is: lxe
= m\\IW2\n-x
(7)
where
( Yy Tz) dv\
dw
+
(du +
dw\ +
\Tz YX)
(8)
These equations are discretized by employing the finite volume method. The procedure for solving the discrete equations is iterative and based upon the SIMPLE solution algorithm [12]. The pressure corrections are calculated using a preconditioned conjugate gradient algorithm. A vector upwinding scheme is employed, which greatly reduces the false diffusion errors associated with the advection terms. The time differencing is fully implicit, so avoiding Courant stability restriction on the time-step size. The discretized Eqs. (1), (2), and (7) are solved together under appropriate boundary conditions. The no-slip boundary conditions and the rigidness of the bifurcation were taken into
2 Newtonian and Non-Newtonian Models. It has been believed for a long time that blood behaves essentially as a Newtonian fluid during flows through large blood vessels. However, in certain experimental studies [13-15] there is evidence which indicates some non-Newtonian effects even in large arteries. In view of this it is desirable to accommodate some non-Newtonian models of possible blood behavior. Several theories have been postulated to describe such behavior of the blood, for instance, the Casson and Bingham models have found a wide application in the past [16-17]. However, there are still doubts regarding the existence of the yield stress in dynamic situations since it has been observed only under static loading conditions. Therefore, a more appropriate form of the constitutive equation would be a general power law. As shown in Eq. (3), the model utilized in this study is a relatively simple one which depends upon shear rate only. For human blood the constants m and n are taken to be 0.042 Pa-s and 0.61, respectively [13]. A modified power law expression for whole blood which depends upon shear rate, hematocrit, and total protein minus albumin [16] is expected to be adopted in our future study. 3 Wall Distensibility.
Arterial walls are viscoelastic in-
2 An alternative upstream boundary condition comprised a typical average canine velocity waveform supplied us by Dr. Kim Parker, and was allowed to develop in a downstream direction over the relatively straight parent tube upstream of site 1. Although predictions are omitted here for reason of space, they compared substantially less well than those of Fig. 6. Hence it is necessary to use a specific waveform in making detailed bifurcation predictions.
506 / Vol. 114, NOVEMBER 1992 Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 05/22/2015 Terms of Use: http://asme.org/terms
Transactions of the ASME
I10DEL
MODEL
(A)
Fig. 4 Three-dimensional perspective view of axial velocity (vertical axis) versus radial position (r) and phase (time, f) at five different sites (as defined in Fig. 1) in the bifurcation model B
MODEL
CA3
Fig. 3(6)
jO.Bm/s
(Ap)
(Cp) MODEL
(B)
Fig. 3(c) Fig. 3 Axial velocity profiles of a Newtonian fluid in the bifurcation plane at four different phases Ap, Bp, Cp, Dp during a cardiac cycle (see text for explanation)
homogeneous multilayered tissues. They are composed mainly of collagen, elastin fibers and smooth muscle cells. In most of the previous studies, the arterial wall was assumed to be rigid, due to the general belief that its distensibility may be regarded as a "second-order" effect on arterial bifurcation flows and thus negligible. However, some recent experimental studies [13, 15] have revealed the differences between rigid and elastic wall behavior during pulsatile flow, and demonstrated the fact that wall distensibility may play an important role in the understanding of bifurcation flows. The difficulties encountered in the incorporation of the moving wall are twofold: one is to Journal of Biomechanical Engineering
choose a simple and accurate mathematical expression of the mechanical behavior of the arterial wall, i.e., a constitutive equation, and another is to construct an interface between the moving wall and the fluid flow. The former has been partially solved by utilizing the pseudo strain energy functions which were concluded as the most practical approximation of the stress-strain relationship for the arterial wall [18]. The latter requires moving boundary numerical techniques, such as adaptive gridding (see reference [19]). Research is going on in our group to enable the wall behavior to be modeled in a comprehensive manner; however, in the present study a rigid wall is used. 4 Grid Resolution Tests. For the two cases presented, it was assumed that a symmetry plane existed so that simulations were confined to only half of the bifurcation. The computational mesh consists of 5040 elements and 6090 nodes, with 18 elements in the radial direction at each cross-section. The grid chosen is a compromise between accuracy and available computer time and storage. Preliminary calculations were carried out for bifurcation model B under steady flow conditions using three grid distributions, i.e., 12, 18, and 24 elements in the radial direction. Results of these calculations indicated that a grid with 18 elements in the radial direction was sufficient for the present purpose. 3 The average velocity differences were 0.19 percent with the finer grid and 6.17 percent with the coarser grid. With the finer grid (24-18 elements) the maximum axial velocity differed by 1.36 percent and the maximum wall shear stress differed by 2.5 percent. Before application to the current bifurcation study, the computer code ASTEC was tested by making a prediction for fully developed laminar flow in a straight tube with both Newtonian and non-Newtonian fluids. The computer velocity profiles were in good agreement with the theoretical solutions (maximum local velocity error 0.3 percent). The code has also been previously validated for various bifurcation flow cases, i.e., two-dimensional steady and pulsatile flows in a T-junction, three-dimensional steady and pulsatile flows in a T-junction with circular cross-section. The comparison between the computer results and laboratory measurements were very satisfactory and typical presentations are in reference [20]. Results and Discussions Initial flow predictions were performed under Newtonian assumptions for the blood. Figures 3(6) and (c) show the axial velocity profiles in the bifurcation plane at the pulse 3 This grid convergence check confirmed the results of an earlier exercise on a two-dimensional T-junction bifurcation geometry.
NOVEMBER 1992, Vol. 114 / 507
Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 05/22/2015 Terms of Use: http://asme.org/terms
INNER WALL
MODEL B
Fig. 5 Wall shear stress versus time at (a) one diameter above the flow divider in the parent vessel, and (b) 0.2 diameter downstream of the flow divider In the larger daughter vessel on the outer wall, and the inner wall for bifurcation model B
phases of accelerated flow Ap, maximum positive flow Bp, decelerated flow Cp, and maximum reverse flow Dp (as defined in Fig. 3(a)) for the two bifurcation models. Comparison of Fig. 3 (b) and (c) for bifurcation models A and B respectively, shows that the essential features of the flow are the same. Due to space limitation, we present results in the remainder of the paper for model B only. Figure 4 provides a three-dimensional view of the axial velocity at five different sites in the plane of the bifurcation for model B. From these results, it is observed that in the parent' vessel velocity profiles are hardly affected by the flow branching downstream. They are relatively symmetrical with reverse flow occurring during later diastole. In the two daughter vessels, velocity profiles vary with position and phase during the cycle. At the level of the flow divider, velocity profiles are skewed towards the inner wall (at the flow divider side) due to the presence of the apex and the axial momentum of the fluid when it leaves the parent vessel. As the flow moves downstream, the peak axial velocity gradually shifts to the center of the vessel. Near the outer wall of the larger daughter vessel, flow separation does not occur during flow acceleration but is present during later flow deceleration just before flow reversal. This suggests that in bifurcation regions separated flow might not always exist throughout the cycle, its presence possibly depending upon some physical parameters which characterize the flow, such as Reynolds number and flow boundary conditions. This finding is consistent with previous experimental observations [21-22]. Under the current flow conditions, only small differences in velocity profiles at two bifurcations with different branching angles are noticed. Wall shear stress was also calculated for bifurcation model
•o^~4000.9.~
feS^ 4000000000005^
sftoooooooooooooso*^
We*U2!°l
• -»««*"^
-•-a-
•"*»
^onWOnonnrtaa
»'«-«'O"0^«'«'«"«
|c)
(a) AT SITE 1 OF DOG B COMPUTATIONAL RESULTS oa°°«
EXPERIMENTAL RESULTS
508 / Vol. 114, NOVEMBER 1992
(b) at s i t e 2 of dog
AT S I T E 3 OF DOG 1
(d)
AT S I T E
4 OF DOG 1
(e)
AT S I T E
5 OF DOG 1
Fig. 6 Comparison of axial velocity profiles between predictions (—) using a non-Newtonian model and the in vivo measurements (© © o) at five different sites (as defined in Fig. 1) in bifurcation model B. In each case, the right hand side is the outer wall, and the left-hand side is the inner wall (flow divider side).
Transactions of the ASME
Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 05/22/2015 Terms of Use: http://asme.org/terms
1.0 0.2 0.1 0.5 0 (a) AT SITE 1 OF DOG B (SOLID LINE: NEWTONIAN MODEL) (DIAMONDS: NON-NEWTON IAN MODEL)
AT SITE 2 OF DOG I
B when a Newtonian fluid was used. The wall shear stress obtained at one diameter above the flow divider in the parent vessel is presented in Fig. 5(c) as a function of time during the cycle. The peak stress is 2.36 N/m2, occurring during the acceleration phase of the systole, and the time-averaged mean value is 0.28 N/m2. At 0.2 diameter downstream of the flow divider in the large daughter vessel, the wall shear stress behavior is considerably different as shown in Fig. 5(b). On the outer wall the magnitudes of the stress are rather low and the time-averaged mean value is close to zero. On the inner wall the shear stress is much greater reaching the maximum and time-averaged mean values of 5.9 N/m2 and 1.10 N/m2, respectively. The latter is far above that found in the parent vessel. To verify the computational results, comparison between calculations and measurements is necessary. In Fig. 6, predicted axial velocity profiles for the bifurcation model B with non-Newtonian fluid are compared with those of the in vivo measurements at five different sites shown in Fig. 1. The predictions represent our latest model. However, there are several facts which must be taken into account in making the comparison. First of all, the measurements were performed by using the Doppler ultrasound device which measures the velocity component in the direction of beam but calculates the velocity as though it is parallel to the vessel axis. This assumption is likely to give misleading results in such complex flows. Secondly, the arterial wall was treated as rigid in the predictions due to difficulties in coding a moving boundary problem at the current stage. According to the experimental studies published so far, the wall distensibility might have a significant, even if local, effect on arterial bifurcation flows, especially when it is coupled with the non-Newtonian character of the blood. Thirdly, there are uncertainties about the velocity very close to the wall and the position of the moving wall itself in the measurements. Therefore, the comparison between the predictions and the measurements has to be made very carefully, and at this stage it can only be done with diffidence. In Fig. 6, it is observed that during diastole there is a very good agreement between the calculations and the measurements in parent and both daughter vessels at all sites. Good agreement is also obvious in the parent vessel at site 1, and in the larger daughter vessel at site 4, where the velocity profiles are rather symmetric. Some differences are noticed in the daughter vessels at sites 2, 3, and 5, which are near the flow divider, during Journal of Biomechanical Engineering
0.6 0.8 1.0 AT SITE 3 OF DOG B
0 Id) AT SITE 4 OF DOG B (SOLID LINE: NEWTONIAN MODEL) (DIAMONDS: NON-NEWTONIAN MODEL)
0.2
0.1
0.6
0.8
t .0 (x/d)
(e) AT SITE 5 OF DOG !
Fig. 7 Comparison of axial velocity profiles of a Newtonian (—) and a non-Newtonian fluid (o o o) for bifurcation model B at five different sites. In each case, the right-hand side is the outer wall, and the left hand side is the inner wall.
the acceleration phase of systole. The predicted velocity profiles show greater skewing towards the inner wall. We feel this is mainly due to the assumption of parallel flow made in the measurements and the rigidness of the vessel wall adopted in the calculations. In order to investigate quantitatively the non-Newtonian effects on the bifurcation flow field, calculated axial velocity profiles of a Newtonian fluid are compared with those of a non-Newtonian fluid for bifurcation model B at five different sites. As shown in Fig. 7, in general, the velocity profiles appear rather flattened in the non-Newtonian fluid case. This flow NOVEMBER 1992, Vol. 114 / 509
Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 05/22/2015 Terms of Use: http://asme.org/terms
MODEL
(B)
(Bp)
A2
(Ap)
Fig. 8
B
(Cp)
Secondary velocity vector fields for bifurcation model B
behavior has already been demonstrated by others [14]. At the bifurcation region, no large variations in the velocity distribution are observed, which indicates that the non-Newtonian characteristics of the blood may not play an important role in determining the general flow field at these bifurcations for the current flow conditions. However, it is still possible that these characteristics may significantly affect a flow situation which is sensitive to small changes. Such a postulation is in accordance with previous experimental studies [13-15], In these the velocity distributions behind the bifurcation are substantially different in Newtonian and non-Newtonian fluids. In general, then, the quantitative effects of non-Newtonian viscosity appear to depend very much upon the bifurcation geometry and the flow conditions. It is, therefore, valuable to have the option of a non-Newtonian model in the code. Since the in vivo data essentially gives axial flow information, we have concentrated so far on this in the paper. However, the predictions provided a full three-dimensional flow field. In Fig. 8 are presented typical secondary velocity vector fields for bifurcation model B at three different phases Ap, Bp, and Cp as defined. These velocity vectors represent therefore the remainder of the flow field information, and give an indication of any swirling type structures. The results are qualitatively consistent with other published data for bifurcations [8, 20]. It has previously been mentioned that the upstream boundary condition assumed a zero secondary velocity field, due to the lack of measuring this in experiment. Such an approximation is analogous to the accepted procedure for a numerical experiment to test equation stability. For example, in Collins [23] Appendix, even a substantial deliberate error in axial velocity was shown rapidly to diffuse out to a negligible effect, thus demonstrating good stability in the algorithm. Here, the secondary velocities at the near inlet position A1-A2 show that they are relatively small during positive flow. Although during reverse flow they are slightly distorted, these errors are shown
to diminish downstream. However, it would be preferable to include measurements of secondary as well as axial velocities for an upstream boundary condition in future studies. Conclusions To provide predictive data for in vivo canine femoral bifurcation measurements, a numerical investigation was carried out for the tested bifurcations under physiological flow conditions by means of a three-dimensional fluid flow code and a supercomputer. The results obtained for the rigid (the only model approximation) canine femoral bifurcation models with (a) a Newtonian, and (b) a non-Newtonian fluid, are summarized as follows: 1) Under the current flow conditions, i.e., average inflow Reynolds number of 92 and 108, Womersley parameter of 2.04 and 2.39 for bifurcation models A and B, respectively, flow separation was not obvious and existed only for a small period of the whole cycle. 2) Shear stress on the outer wall of the larger daughter vessel was very low and the time-averaged mean value was close to zero. Shear stress on the inner wall of the larger daughter vessel was much greater than either on the outer wall or in the parent vessel. 3) While there was good qualitative and not unreasonable quantitative agreement between the predicted velocity profiles and those of the in vivo measurements, they differed somewhat in detail in the daughter vessels in the branching region, indicating the importance of the secondary flow influences on the flow patterns there. 4) For the two bifurcation models and the current flow conditions, the non-Newtonian characteristic of the blood was not found to be an important factor in determining the general flow field parameters. The numerical results obtained are believed to be important
510 / Vol. 114, NOVEMBER 1992 Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 05/22/2015 Terms of Use: http://asme.org/terms
Transactions of the ASME
in interpreting the measured velocity profiles and acquiring more information on the flow field in bifurcation regions. The present study represents an initial attempt to provide comparable results to in vivo measurements of the flow in arterial bifurcations. Work is continuing to enable wall distensibility effects to be studied. When these are incorporated into the predictions, a detailed quantitative assessment will be made of both experiment and theory.
Acknowledgments The authors are very grateful to Dr. K. H . Parker for his help and discussion, and especially for providing detailed experimental data. This project was supported by a British Heart Foundation Research Grant. Finally, ASTEC is a proprietary code of the United Kingdom Atomic Energy Authority. References 1 Jones, C. J. H., Lever, M. J., Ogasawara, Y., Parker, K. H., Hiramatsu, O., Mito,K.,Tsujioka, K.,andKajiya, F., "Blood Velocity Distributions Within Intact Canine Arterial Bifurcations," American Journal of Physiology: (Heart and Circulatory Physiology, Vol. 31, 1992, pp. H1592-1599. 2 O'Brien, V., and Ehrilich, L. M., "Simulation of Unsteady Flow at Renal Branches," Journal of Biomechanics, Vol. 10, 1977, pp. 623-631. 3 Liepsch, D., and Moravec, S.,' 'Measurement and Calculations of Laminar Flow in a Ninety Degree Bifurcation," Journal of Biomechanics, Vol. 15, 1982, pp. 473-485. 4 Perktold, K., and Hilbert, D., "Numerical Simulation of Pulsatile Flow in a Carotid Bifurcation Model," Journal of Biomechanics, Vol. 8, 1986, pp. 193-199. 5 Rindt, C. C. M., Vosse, F. N. N. V. D., Steenhoven, A. A. V., and Janssen, J. D., "A Numerical and Experimental Analysis of the Flow in a TwoDimensional Model of the Human Carotid Bifurcation," Journal of Biomechanics, Vol. 20, 1987, pp. 499-509. 6 Wille, S., "Numerical Simulation of Steady Flow Inside a Three-Dimensional Aortic Bifurcation Model," ASME JOURNAL OF BIOMECHANICAI ENGINEERING, Vol. 6, 1986, pp. 49-55. 7 Dinnar, U., Enden, G., and Israeli, M., " A Numerical Study of Flow in a Three-Dimensional Bifurcation," Cardiovascular Systems Dynamics Society Meeting, 1988, Canada. 8 Rindt, C. C. M., "Analysis of the Three-Dimensional Flow Field in the
Carotid Artery Bifurcation," PhD Thesis, 1989, Eindhoven University of Technology, The Netherlands. 9 Perktold, K., and Peter, R., "Numerical 3D-Simulation of Pulsatile Wall Shear Stress in an Arterial T-Bifurcation Model," ASME JOURNAL OF BIOMECHANICAL ENGINEERING, Vol. 12, 1990, pp. 2-13.
10 Perktold, K., and Resch, M. "Numerical Flow Studies in Human Carotid Artery Bifurcations: Basic Discussion of the Geometric Factor in Atherogenesis," ASME JOURNAL OF BIOMECHANICAL ENGINEERING, Vol. 12, 1990, pp. 111-
123. 11 Stainsby, R., " A Guide to Using the SOPHIA Mesh Generator," UKAEA (Draft), 1989. 12 Patankar, S. V.] and Spalding, D. B., "A Calculation Procedure for Heat, Mass and Momentum Transfer in Three-Dimensional Parabolic Flows," International Journal of Heat and Mass Transfer, Vol. 15, 1972, pp. 1787-1800. 13 Ku, D. N., and Liepsch, D., "The Effects of Non-Newtonian Viscoelasticity and Wall Elasticity on Flow at a 90° Bifurcation," Biorheology, Vol. 23, 1986, pp. 359-370. 14 Moravec, S., and Liepsch, D., "Flow Investigations in a Model of a ThreeDimensional Human Artery with Newtonian and Non-Newtonian Fluids, Part I , " Biorheology, Vol. 20, 1983, pp. 745-759. 15 Liepsch, D. W., "Flow in Tubes and Arteries—A Comparison," Biorheology, Vol. 23, 1986, pp. 395-433. 16 Walburn, F. J., and Schneck, D. J., " A Constitutive Equation for Whole Human Blood," Biorheology, Vol. 13, 1973, pp. 201-210. 17 Nakamura, M., and Sawada, T., "Numerical Study on the Flow of a NonNewtonian Fluid Through an Axisymmetric Stenosis," ASME JOURNAL OF BIOMECHANICAL ENGINEERING, Vol. 110, 1988, pp. 137-143.
18 Fung, Y. C , Pronek, K., and Patitucci, P., "Pseudoelasticity of Arteries and the Choice of its Mathematical Expression," American Journal of Physiology, Vol. 237, 1979, pp. H620-H631. 19 Xu, X. Y., and Collins, M. W., " A Review of the Numerical Analysis of Blood Flow in Arterial Bifurcations," Proc. Instn. Mech. Engrs., Vol. 204, Part H, Journal of Engineering in Medicine, 1990, pp. 205-216. 20 Collins, M. W., and Xu, X. Y., "A Predictive Scheme for Flow in Arterial Bifurcations: Comparison with Laboratory Measurements," Biomechanics Transport Processes, ed., F. Mosora, C. Caro et al., Plenum Press, New York and London, 1990, pp. 125-133. 21 Richardson, P. D., and Christo, J. L., "Flow Separation Opposite a Side Branch," Biofluid Mechanics, Blood Flow in Large Vessels, ed., D. Liepsch, Springer-Verlag, 1990, pp. 275-284. 22 Cho, Y. I., Back, L. H., and Crawford, D. W., "Experimental Investigation of Branch Flow Ratio, Angle, and Reynolds Number Effects on the Pressure, and Flow Fields in Arterial Branch Models," ASME JOURNAL OF BIOMECHANICAL ENGINEERING, Vol. 107, 1985, pp. 257-267.
23 Collins, M. W., "Finite Difference Analysis for Developing Laminar Flow in Circular Tubes Applied to Forced and Combined Convection," International Journal for Numerical Methods in Engineering, Vol. 15, 1980, pp. 381-404.
(Contents continued) Announcements and Special Notices 460 Transactions Change of Address Form 496
Errata on a Previously Published Paper
553 Calendar of Events
Journal of Biomechanical Engineering
NOVEMBER 1992, Vol. 114 / 511
Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 05/22/2015 Terms of Use: http://asme.org/terms
