J Binmrchann
1976. Vol
9. pp 575-580.
Pergamon Press
Prmtcd m Great Bntam
PULSATILE FLOW THROUGH A BIFURCATION APPLICATIONS TO ARTERIAL DISEASE*
WITH
R. C. FERNANDEZ, K. J. DE Wtrr and M. R. BOTWIN Departments
of Chemical and Civil Engineering, College of Engineering, The University of Toledo.
Toledo, OH 43606. U.S.A. Abstract-Fatty streaks and fibrous plaques are known to be localized in bends, junctions, and branches of the large arteries, suggesting a significant role for fluid mechanics in the genesis of arterial disease. Various mechanisms have been proposed in the literature based on the known occurrence of eddies and high surface traction. For the purpose of helping in the evaluation of these mechanisms. the flow behavior in a bifurcation was chosen for study. The equations of continuity and motion were solved using the marker and cell numerical technique for the two-dimensional, periodic flow of a Newtonian fluid through a horizontal, bifurcating, rigid channel. Both regions of eddying at the outer wall and of high shear stress at the inner wall were shown by the calculations for both steady and pulsatile flows. The eddies disappeared completely during a portion of each cycle for pulsatile flow.
1. INTRODUCTION
Although in the early stages of development morphological generalizations are not easily made about arterial lesions, two striking features stand out. Arterial lesions (1) tend to be localized in bends, junctions and branches of large arteries, and (2) are characterized by the presence of lipids. Mitchell and Schwartz (1965) broadly classify arterial lesions as either fatty streaks or raised plaques. With respect to localization in the aorta, they report that fatty streaks tend to be found in the aortic arch and upper thoracic aorta, while raised plaques are more commonly found in the abdominal aorta. Another interesting feature of localization is that fatty streaks have a tendency to spare the aortic wall around and immediately distal to the ostia of the intercostal arteries, while in direct contrast these same regions are common sites for raised plaques. Calcification is commonly found at these as well as other branch openings. The presence of lipids in lesions cannot explain the localization of arterial disease. In addition their presence does not necessarily imply their participation in atherogenesis. It is possible that they appear only later during growth, but this speculation has not yet been substantiated. The localization of arterial disease could conceivably be explained by structural or functional pecularities of the vessel wall. But Caro et al. (1969) dismiss this explanation and assert that the innominate (brachiocephalic) artery and the left common carotid artery are indistinguishable in structure and are not strikingly different in function, yet one (the carotid) is involved in arterial disease and the other is spared. Understandably, investigators have turned to arterial blood mechanics. A study of the characteristics of flow in the sites of predilection could reveal their common features. Mechanisms could then be proposed and tested. Although some mechanisms will be * Receiced 21 August 1975. I, M 9,9--u
precluded, more than one mechanism can be accommodated. The various mechanisms proposed can be classified into: (a) autoregulatory response to reversals of wall shearing stress; (b) direct damage to the arterial wall; (c) wall starvation; (d) material deposition on the endothelium from the blood; and (e) suppression of mass transfer from the wall to the blood. Reversals of wall shearing stress necessitates pulsatile flow. This fact may explain the involvement of arteries instead of veins. Fry (1968) attributed wall damage to mechanical drag or shear stress that eroded the endothelium. Fry reports that exposure of canine walls to a shear stress of 380dyn/cm’ for as short as one hour causes degradation. He does not quote any value for humans. Texon (1960) blamed Bemoullitype suction that lifted the intima, and Spurlock er al. (1968) have advanced impaction. Starvation, deposition, and suppression can be due to partial isolation of a region, variously called stasis. eddy, separation, etc. from the main stream. Mass transfer between the stasis and the mainstream is minimal but energy transfer is not. Nutrients in the mainstream are prevented from reaching the wall (Dintenfass, 1962), and lipids are kept from leaving the wall (Car0 et at., 1969). The materials trapped in the stasis (Mustard et al., 1963; Fox and Hugh, 1966) may be healthy erythrocytes or thrombocytes which, if imprisoned long enough (life span of an erythrocyte is 90-120 days), may lyse. The released products, along with fibrin and other debris (Merrill, 19663, could interact with the wall leading to its degeneration. Low shear rates in the stasis could also allow agglutination and syneresis of otherwise healthy trapped cells, again leading to unfavorable results for the wall. 2. PREVIOUSWORK ON THE BIFURCATION Attinger (1964) and Crowe (1969) experimentally for three-dimensional flow using
studied flow patterns 575
576
R. C. FERNANDEZ, K. J. DE WILT and M. R. J~OTWIN
birefringence techniques. For pulsatile flow, Attinger observed the on-set of turbulence in the branches at an upstream Reynolds number of 1700 and cavitation at the outer walls of the branches at a Reynolds number of cu. 10,000. At these high Reynolds numbers, disturbances were sucked back into the mother tube during flow deceleration. Crowe found high shear stresses at the inner wall about a diameter downstream from the stagnation point. He also observed separation along the outer wall for an upstream Reynolds number between 20 and 50 depending on the branching angle. More recently, Brech and Bellhouse (1973) studied flow patterns and velocity and shear stress distributions downstream of symmetrical branches using thin film probes. The model branches were rigid and water was used as the fluid medium. The ratio of the cross-sectional area of the combined daughter tubes to the parent tube was 1.12, which compares with the value of 1.0 for arterial branches as measured by Caro et al. (1971). For both steady and pulsatile flow, their results showed high shear on the inside wall of the daughter tube ca 0.2 diameters downstream, but no regions of stasis were observed on the outside wall. They attribute this latter fact to secondary flow characteristics. The upstream Reynolds numbers used were 750 and 1500. Davids and Cheng (1972) applied a finite element numerical analysis to steady flow through a bifurcation, and reported high stresses at the centerline of the crotch and peak stress at the stagnation point. Their analysis was restricted to constant cross-sectional areas (from mother to daughters) and very small branching angles. All of the above experiments and theoretical analyses were on tubes of circular cross-section. The numerical studies by Naff (1969), Lynn er af. (1972), and O’Brien et al. (1973) are restricted to a two-dimensional bifurcation. In these three studies computer solutions were found for the vorticity and stream functions. Both Naff and O’Brien et al. indicate that separation occurs for Reynolds numbers near the physiological range, but Lynn et al. report no separation for steady flow with Reynolds numbers up to 1,000. The Reynolds numbers in these three investigations were based on the width of the inlet channel. Only the work by O’Brien et al. considered pulsatile as well as steady flow. In the present study, the original variables of velocity and pressure are retained in the numerical work. and the flow results are easily presented. 3. TAECOVERNINCEQUATlONSANDNUMERlCAL PROCEDURE The problem considered herein is the periodic flow of a homogeneous, incompressible, Newtonian fluid through a horizontal, bifurcating, rectangular, rigid channel (Fig. 1). The channel was conceived as being deep enough so that, at some horizontal section
Fig. 1. The model bifurcation. between the top and bottom walls, the viscous effects of these walls could be ignored. The llow patterns can then be considered as being planar, or twodimensional. Since the channel is symmetric, we may restrict our analysis to the left half only. This part can be further treated as the coupling of two sections, the mother and the daughter. Each part can then be analyzed separately without forgetting to insure conservation of mass and momentum at their line of union. The total branching angle, K, was taken to be 90”. 3.1 Field equations For a homogeneous, incompressible fluid of constant viscosity, the equations of conservation of mass and momentum (Navier-Stokes) reduce to Continuity:
v*ti=o
at2 &
+V*(tx)=
-+vyo
x (v x iI-)],
(2)
where u’ is the velocity vector, cm/set, p is the pressure, dyn/cm2, p is the fluid density, 1.05 g/en?, v is the fluid kinematic viscosity, 0.0381 cm2/sec, t is physical time, set, and V is the de1 or nabla operator. These equations are next transformed into dimensionless form and applied to the two-dimensional bifurcation. For the daughter portion, the quantities U = ub/v and V = vb/v are defined, where u and v are the axial and transverse velocities, respectively. The quantity b is onahalf the width of the daughter portion, l.OOcm. Also X = x/b, Y = y/b, a = wb2/v, 0 = tv/b2 and P = pb’/pv’, where o, Hertz, is the frequency of the forcing function, and x and y are the axial and transverse coordinates, respectively. Substituting these definitions into equations (1) and (2), the equations for the two-dimensional daughter portion of the bifurcation become Continuity:
au
(ygf+z=o.
av
_ (3)
$7:
Pulsatile flow through a bifurcation
pulse at the exit. In reality, the fluid is compressible and the wall is distensible so that the pulse wave velocity is not infinite, being about several meters per second. However, for the short segment considered here the assumption of infinite wave velocity is quite Momentum-Y: good. In dimensionless form. with I, being the length of (5) the daughter portion as measured from the apex of the bifurcation centerlines. For the mother portion, the same process is folAt X = 0, P,, = P, + P3 sin (2naB) lowed. resulting in the same field equations, including (mother entrance) (9) the use of h as the characteristic length. At X = 1,/b, P,., = P2 + P4 sin (2na0) 3.2 Boundary and initial conditions (daughter exit) Momentum - X:
The vessel walls are assumed to be impervious, i.e. there is no mass flow through them. In addition, they are assumed to be rigid so that transverse velocity there can be set to zero. Finally, it is assumed that there is no slip at this rigid wall so that the axial velocity can also be set to zero. For the mother portion. pressure and axial velocity are symmetric with respect to the centerline. Transverse velocity, however. is antisymmetric with respect to the centerline, and is zero at the centerline. In dimensionless variables the transverse boundary conditions for the mother and daughter portions are 1. mother portion At Y = 0, g
= 0 and I’ = 0. (centerline)
At Y = a/b, U = 0 and V = 0. (mother wall),
(6)
where a is one-half the width of the mother portion, 1.41 cm. 2. daughter portion At Y = - 1, U = 0 and V = 0. (inner wall) (7) At Y = + 1. U = Oand V = O.(outerwall). The axial boundary conditions used were that of fully-developed parabolic flow at the mother entrance and the daughter exit. Regardless of the initial state of the fluid, it will respond to the periodic forcing function until it is in harmony with it. Analogously, the values of the field variables will move toward a steady periodic condition even if initially they are given zero or any other arbitrary values. Hence 3.At0=O,U=OandV=OatallXandY locations for both mother and daughter portions. (8)
so that AP = P,., - Pi, = P, + P6
sin
(2naf3).
(10)
The values used for P5 and P6 were 1210 and 1100, respectively. The value of P5 was determined by performing the calculations for the steady flow case Reynolds number (a = O), where an upstream (Re = Zau,,/v) of 146 was obtained, and the value of P6 was then set to give the flow a reasonable periodic behavior and Reynolds number variation. The dimensionless frequency, a, was set equal to 1.0 for the periodic flows, where instantaneous upstream Reynolds numbers varied from 130 minimum to 173 maximum, the mean value being different from 146 because of the problem nonlinearity. The quantity Id was determined by performing the calculations for various lengths until fully developed flow was obtained at the daughter outlet. It was found that an Id of 6 cm satisfied this fully developed criteria. 3.4 Numerical procedure The above governing equations were solved numerically using the Marker and Cell (MAC) method, which was developed at the Los Alamos Scientific Laboratory and is the subject of a comprehensive report by Welch et al. (1966). Unlike the stream function-vorticity method, the MAC method involves the primary variables U, V, and P, thus enabling the direct application of the physical boundary conditions. In the present study, the MAC method was extended in that it was necessary to set up two grids, one for the mother and one for the daughter, and to couple the two grids by weaving pressure and velocity values along their line of union. Six forcing cycles were needed before steady periodic conditions were achieved. The computer used was the IBM 360/44. Complete details of the numerical work are given by Fernandez (1972).
3.3 Forcing ,funcrion In order to simulate the pulsatile nature of arterial flow, a sinusoidal periodic forcing function was used. Because the fluid is incompressible and the wall is rigid, the pulse wave will propagate at infinite velocity from the flow entrance at the mother portion to the flow exit at the daughter portion. Hence, there is no phase fag between the pulse at the entrance and the
A INTERPRETATION
OF COMPUTER RESULTS
Although the model used is two-dimensional and has a sharp corner, the results seem to be in fair agreement with Crowe’s (1969) experimental results from branched configurations formed by the intersection of axisymmetrical tubes. Both regions of eddying
578
R. C. FERNANDBZ,K. J. DE WITT and M. R. B~TWIN
at the outer wall and of high shear stress at the inner wall were found in our calculations for both steady and pulsatile flows. Since the Reynolds numbers were not high enough, the eddy was not sucked back into the mother channel during flow deceleration. The position of high shear stress was also unchanged by the pulse. No secondary flow, which Brech and Bellhouse (1973) attribute in their experiments as the reason for preventing the appearance of zones of stasiq could be handled in the present study because of the two-dimensional model used. Figure 2 shows the results for steady’ flow at upstream Reynolds number of 146. The flow is from left to right. The arrows are the total velocity vectors, i.e. the vector sum of transverse velocity, V, and axial velocity, U. The heavy curves are the velocity profiles at the locations indicated by the heavy vertical lines, and are approximately the same as the usual velocity profiles drawn for the axial velocity alone. The approximation is not too much in error since the transverse velocities at the locations considered are nearly zero. Upstream, the profile is parabolic. Downstream the profile shows flow reversal at the outer wall near the comer and also steepness indicating high shear rate at the inner wall. The rectangular portion at the outer wall is also shown magnified four times in area. Here only the axial velocity is plotted so that the envelope shown by the dotted line going through points of reversal may be accurately drawn. Farther downstream the profile shows modulation and a trend towards recovery of the new parabolic shape. For pulsatile flow between the instantaneous Reynolds numbers of 130 to 173, the same general features are evident; however, there is the added feature of growth and disappearance of the envelope as the pressure forcing function goes through a cycle. At a Reynolds number of 130 (25” into the pressure pulse), no envelope can be seen as is illustrated in Fig. 3. Figure 4 shows a small envelope starting to appear at a Reynolds number of 142 (68” into the pressure cycle). As the Reynolds number increases the
It-i
/
Fig. 3. (a) Velocity profiles for pulsatile flow for Re = 130 (25” into pressure pulse). (b) Enlarged view outside wall region.
envelope grows both upstream and downstream and also away from the wall into the mainstream. Growth in the downstream direction is thd most extensive and that towards the mainstream is the least extensive. Figure 5 shows the maximum Reynolds number reached of 173 (158’ into the pressure cycle), but the size of the envelope is not maximum then. Instead, as seen in Fig. 6, it reaches a maximum at a Reynolds number of 169 (202” into the pressure cycle). Afterwards, the envelope contracts and eventually disap pears, having been present over the interval from approx. 60” to 335” of the pressure cycle. These results indicate that, for pulsatile flow, the instantaneous Reynolds number is not sticient for specifying the flow, as the pressure pattern must also be specified. Similar qualitative results using different numerical techniques were obtained by O’Brien et al. (1973). Pressures around the stagnation point on the inner wall were on the order of 1.7 times as great as the pressures at the input boundary. This pressure difference could be of importance in evaluating the impaction mechanism if the values required for damage of vessel wall or blood cells themselves were known, since the difference in pressures is a measure of the
It-i Fig. 2. (a) Velocity profiles for steady flow for Re = 146. (b) Enlarged view of the separation region.
--
---L--f
*
Fig. 4. (a) Velocity profiles for pulsating flow for Re = 145 (68” into pressure pulse). (b) Enlarged view of separation region.
579
Pulsatile flow through a bifurcation
---_) /
Fig. 5. (a) Velocity profiles for pulsatile flow for Re = 173 (158” into pressure pulse). (b) Enlarged view of separation region.
square of the velocity of impaction. In the region of stasis or eddying at the outer wall, the pressure variance among neighboring cells was less than 0.15 dyn/cm’. This fact somewhat questions Texon’s (19-60) suction hypothesis unless the endothelium proves to be extremely pliant. A region of high shear stress was found at the inner wall about a daughter half-width downstream from the stagnation point, but the highest shear stress at a Reynolds number of 146 for the steady case was ca. ‘10 dyn/cm2, which is an order of magnitude smaller than either the experimental value of 107 dyn/cm’ reported by Crowe (1969) for a tube Reynolds number of 100, or the extrapolated experimental value for blood of 220 dyn/cm2 by Brech and Bellhouse (1973) for a Reynolds number of 750. For the pulsatile case at a Reynolds number of 173, the highest wall shear stress was 12 dyn/cm’. No value was reported by Crowe for the pulsatile case, but Brech and Bellhouse’s extrapolated blood value is 310 dyn/cm’. These values, expecially those. from the present study, are all lower than the 380dyn/cm2 for wall tearing reported by Fry (1968) for dogs. These
----w
/
Fig. 6. (a) Velocity profiles for pulsatile flow for Re = 169 (202” into pressure pulse). (b) Enlarged view of separation region.
values Seem to suggest that if direct wall damage were to be a plausible mechanism then the pathological values for damage to humans should be smaller than that for canines. The periodic appearance of the eddy at the outer wall suggests a further investigation into the plausibility of the proposed mechanism on autoregulatory response. It also provides a test on the mechanisms based on mass transfer, such as deposition and suppression. If it can be shown that the disappearance of the eddy allows surface renewal, i.e. material previously trapped in the eddy has a chance to escape allowing influx of fresh material to the outer wall, then these mechanisms will have to be revised to incorporate this added phenomena. S. DISCUSSION The flow behavior of a Newtonian fluid in a twodimensional bifurcation has been investigated numerically in order to aid in the evaluation of current, proposed mechanisms for atherogenesis. The key results for pulsatile flow illustrate that a region of separation exists on the outer wall of the daughter channel over a large portion of the pressure cycle, the possible effects of which have been interpreted with regard to these existing atherogenesis theories. In order to completely evaluate all of the presently proposed mechanisms, further work should include the experimental and numerical study of particle trajectories, the determination of the stress magnitudes necessary to damage vessel walls, platelets, erythrocytes, etc., and the determination of the nature and magnitude of the forces between cells that would cause them to agglutinate or undergo sedimentation.
REFERENCES
Attinger, E. 0. (1964) Flow patterns and vascular geometry. In Pufsarile Blood Flow. (Edited by Attinger, E. O.), p. 179. McGraw-Hill, N.Y. Brech. R. and Bellhouse, B. J. (1973) Flow in branching vessels. Cardiouasc. Res. 7, 593400. Caro, C. G., Fitz-Gerald, J. M. and Schroter. R. C. (1969) Arterial wall shear and distribution of early atheroma in man. Nature 223, 1159-1161. Caro, C. G., F&x-Gerald, J. M. and Schroter, R. C. (1971) Atheroma and arterial wall shear: observation, correlation and proposal of a shear dependent mass transfer mechanism for atherogenesis. Proc. Roy. Sot. Lond. B. 177, 109-159. Crowe, Jr, W. J. (1969) Studies of arterial branching in models using flow birefringence. PhD. Thesis, The University of Florida, Gainesville, Florida. Davids, N. and Cheng, R. (1972) Finite element analysis of flow in blood vessels with arbitrary cross-section. Progress Report, The Pennsylvania State University, University Park, Pennsylvania. Dintenfass, L. (1962) Thixotropy of blood and proneness to thrombus formation. Circ. Res. 11, 233-239. Fernandg R. C. (1972) Pulsatile flow through a bifurcation with applications to arterial disease. PhD. Thesis, The University of Toledo, Toledo, OH.
580
R. C. FERNANDEZ, K. J. DE Wtrr and M. R.
Fox. J. and Hugh, A. E. (1966) Localization of atheroma: a theory based on boundary layer separation. Brit. Heart J. 28, 388-399. Fry, D. L. (1968) Acute vascular endothelial changes associated with increased blood velocity gradients. Circ. Res. 22, 165-197. Lynn, N. S., Fox. V. G. and Ross, L. W. (1972) Computation of fluid-dynamical contributions to atherosclerosis at arterial bifurcations. Biorheology 9, 61-66. Merrill. W. W., Gilliland, E. R., Lee, T. S. and Salzman, E. W. (1966) Blood rheology: effect of fibrinogen deduced by addition. Circ. Res. 13, 437-446. Mitchell, J. R. A. and Schwartz, C. J. (1965) The morphology of arterial plaques. In Arterial Disease, p. 18. Blackwell, Oxford. Mustard, J. F., Rowsell, H. C., Murphy, E. A. and Downie, H. G. (1963) Intimal thrombosis in atherosclerosis. In Euolurion of the Atheroscleroric Plaque (Edited by Jones, R. L.), p. ,181. University of Chicago Press, Chicago, IL. Naff, S. A. (1969) Accelerated laminar flow through a twodimensional bifurcation. M. S. Thesis, Carnegie-Mellon University, Pittsburgh, PA. O’Brien, V., Ehrlich, L. W. and Friedman, M. H. (1973) Nonlinear simulation of unsteady arterial flows in a branch. Presented at the American Institute of Chemical Engineers Annual Meeting, Philadelphia, PA. Spurlock, 1. M., Homa, M. and Durfee, R. L. (1968) Atherosclerosis etiology: further evidence of the role of hemodynamic effects. Research Monograph of the Atlantic Research Corporation, Alexandria, VA.
BOIWIN
Texon, M. (1960) The hemodynamic concept of atherosclerosis. Bull. N.Y. Acad. Med. 36, 163-274. Welch, J. E., Harlow, F. H., Shannon, J. P. and Daly, B. J. (1966) The MAC Method-A computing technique for solving viscous, incompressible, transient fluid-flow problems involving free surfaces. Los Alamos Scientific Laboratory. Los Alamos, NM.
NOMENCLATURE axial and transverse coordinates, cm dimensionless axial and transverse coordinates velocity vector, cm/set axial and transverse velocities, cm/set dimensionless axial and transverse velocities half-widths of the mother and daughter channels, cm length of daughter channel, cm pressure, dyn/cm* dimensionless pressure = pb*/pv* fluid density = 1.05 gm/cm3 fluid kinematic viscosity = 0.0381 cm*/sec time, set dimensionless time = tv/b* forcing frequency, Hertz dimensionless forcing frequency = r&*/v Reynolds number = Zau,,/v average axial velocity, cm/set total branching angle = 90”.
1976. Vol
9. pp 575-580.
Pergamon Press
Prmtcd m Great Bntam
PULSATILE FLOW THROUGH A BIFURCATION APPLICATIONS TO ARTERIAL DISEASE*
WITH
R. C. FERNANDEZ, K. J. DE Wtrr and M. R. BOTWIN Departments
of Chemical and Civil Engineering, College of Engineering, The University of Toledo.
Toledo, OH 43606. U.S.A. Abstract-Fatty streaks and fibrous plaques are known to be localized in bends, junctions, and branches of the large arteries, suggesting a significant role for fluid mechanics in the genesis of arterial disease. Various mechanisms have been proposed in the literature based on the known occurrence of eddies and high surface traction. For the purpose of helping in the evaluation of these mechanisms. the flow behavior in a bifurcation was chosen for study. The equations of continuity and motion were solved using the marker and cell numerical technique for the two-dimensional, periodic flow of a Newtonian fluid through a horizontal, bifurcating, rigid channel. Both regions of eddying at the outer wall and of high shear stress at the inner wall were shown by the calculations for both steady and pulsatile flows. The eddies disappeared completely during a portion of each cycle for pulsatile flow.
1. INTRODUCTION
Although in the early stages of development morphological generalizations are not easily made about arterial lesions, two striking features stand out. Arterial lesions (1) tend to be localized in bends, junctions and branches of large arteries, and (2) are characterized by the presence of lipids. Mitchell and Schwartz (1965) broadly classify arterial lesions as either fatty streaks or raised plaques. With respect to localization in the aorta, they report that fatty streaks tend to be found in the aortic arch and upper thoracic aorta, while raised plaques are more commonly found in the abdominal aorta. Another interesting feature of localization is that fatty streaks have a tendency to spare the aortic wall around and immediately distal to the ostia of the intercostal arteries, while in direct contrast these same regions are common sites for raised plaques. Calcification is commonly found at these as well as other branch openings. The presence of lipids in lesions cannot explain the localization of arterial disease. In addition their presence does not necessarily imply their participation in atherogenesis. It is possible that they appear only later during growth, but this speculation has not yet been substantiated. The localization of arterial disease could conceivably be explained by structural or functional pecularities of the vessel wall. But Caro et al. (1969) dismiss this explanation and assert that the innominate (brachiocephalic) artery and the left common carotid artery are indistinguishable in structure and are not strikingly different in function, yet one (the carotid) is involved in arterial disease and the other is spared. Understandably, investigators have turned to arterial blood mechanics. A study of the characteristics of flow in the sites of predilection could reveal their common features. Mechanisms could then be proposed and tested. Although some mechanisms will be * Receiced 21 August 1975. I, M 9,9--u
precluded, more than one mechanism can be accommodated. The various mechanisms proposed can be classified into: (a) autoregulatory response to reversals of wall shearing stress; (b) direct damage to the arterial wall; (c) wall starvation; (d) material deposition on the endothelium from the blood; and (e) suppression of mass transfer from the wall to the blood. Reversals of wall shearing stress necessitates pulsatile flow. This fact may explain the involvement of arteries instead of veins. Fry (1968) attributed wall damage to mechanical drag or shear stress that eroded the endothelium. Fry reports that exposure of canine walls to a shear stress of 380dyn/cm’ for as short as one hour causes degradation. He does not quote any value for humans. Texon (1960) blamed Bemoullitype suction that lifted the intima, and Spurlock er al. (1968) have advanced impaction. Starvation, deposition, and suppression can be due to partial isolation of a region, variously called stasis. eddy, separation, etc. from the main stream. Mass transfer between the stasis and the mainstream is minimal but energy transfer is not. Nutrients in the mainstream are prevented from reaching the wall (Dintenfass, 1962), and lipids are kept from leaving the wall (Car0 et at., 1969). The materials trapped in the stasis (Mustard et al., 1963; Fox and Hugh, 1966) may be healthy erythrocytes or thrombocytes which, if imprisoned long enough (life span of an erythrocyte is 90-120 days), may lyse. The released products, along with fibrin and other debris (Merrill, 19663, could interact with the wall leading to its degeneration. Low shear rates in the stasis could also allow agglutination and syneresis of otherwise healthy trapped cells, again leading to unfavorable results for the wall. 2. PREVIOUSWORK ON THE BIFURCATION Attinger (1964) and Crowe (1969) experimentally for three-dimensional flow using
studied flow patterns 575
576
R. C. FERNANDEZ, K. J. DE WILT and M. R. J~OTWIN
birefringence techniques. For pulsatile flow, Attinger observed the on-set of turbulence in the branches at an upstream Reynolds number of 1700 and cavitation at the outer walls of the branches at a Reynolds number of cu. 10,000. At these high Reynolds numbers, disturbances were sucked back into the mother tube during flow deceleration. Crowe found high shear stresses at the inner wall about a diameter downstream from the stagnation point. He also observed separation along the outer wall for an upstream Reynolds number between 20 and 50 depending on the branching angle. More recently, Brech and Bellhouse (1973) studied flow patterns and velocity and shear stress distributions downstream of symmetrical branches using thin film probes. The model branches were rigid and water was used as the fluid medium. The ratio of the cross-sectional area of the combined daughter tubes to the parent tube was 1.12, which compares with the value of 1.0 for arterial branches as measured by Caro et al. (1971). For both steady and pulsatile flow, their results showed high shear on the inside wall of the daughter tube ca 0.2 diameters downstream, but no regions of stasis were observed on the outside wall. They attribute this latter fact to secondary flow characteristics. The upstream Reynolds numbers used were 750 and 1500. Davids and Cheng (1972) applied a finite element numerical analysis to steady flow through a bifurcation, and reported high stresses at the centerline of the crotch and peak stress at the stagnation point. Their analysis was restricted to constant cross-sectional areas (from mother to daughters) and very small branching angles. All of the above experiments and theoretical analyses were on tubes of circular cross-section. The numerical studies by Naff (1969), Lynn er af. (1972), and O’Brien et al. (1973) are restricted to a two-dimensional bifurcation. In these three studies computer solutions were found for the vorticity and stream functions. Both Naff and O’Brien et al. indicate that separation occurs for Reynolds numbers near the physiological range, but Lynn et al. report no separation for steady flow with Reynolds numbers up to 1,000. The Reynolds numbers in these three investigations were based on the width of the inlet channel. Only the work by O’Brien et al. considered pulsatile as well as steady flow. In the present study, the original variables of velocity and pressure are retained in the numerical work. and the flow results are easily presented. 3. TAECOVERNINCEQUATlONSANDNUMERlCAL PROCEDURE The problem considered herein is the periodic flow of a homogeneous, incompressible, Newtonian fluid through a horizontal, bifurcating, rectangular, rigid channel (Fig. 1). The channel was conceived as being deep enough so that, at some horizontal section
Fig. 1. The model bifurcation. between the top and bottom walls, the viscous effects of these walls could be ignored. The llow patterns can then be considered as being planar, or twodimensional. Since the channel is symmetric, we may restrict our analysis to the left half only. This part can be further treated as the coupling of two sections, the mother and the daughter. Each part can then be analyzed separately without forgetting to insure conservation of mass and momentum at their line of union. The total branching angle, K, was taken to be 90”. 3.1 Field equations For a homogeneous, incompressible fluid of constant viscosity, the equations of conservation of mass and momentum (Navier-Stokes) reduce to Continuity:
v*ti=o
at2 &
+V*(tx)=
-+vyo
x (v x iI-)],
(2)
where u’ is the velocity vector, cm/set, p is the pressure, dyn/cm2, p is the fluid density, 1.05 g/en?, v is the fluid kinematic viscosity, 0.0381 cm2/sec, t is physical time, set, and V is the de1 or nabla operator. These equations are next transformed into dimensionless form and applied to the two-dimensional bifurcation. For the daughter portion, the quantities U = ub/v and V = vb/v are defined, where u and v are the axial and transverse velocities, respectively. The quantity b is onahalf the width of the daughter portion, l.OOcm. Also X = x/b, Y = y/b, a = wb2/v, 0 = tv/b2 and P = pb’/pv’, where o, Hertz, is the frequency of the forcing function, and x and y are the axial and transverse coordinates, respectively. Substituting these definitions into equations (1) and (2), the equations for the two-dimensional daughter portion of the bifurcation become Continuity:
au
(ygf+z=o.
av
_ (3)
$7:
Pulsatile flow through a bifurcation
pulse at the exit. In reality, the fluid is compressible and the wall is distensible so that the pulse wave velocity is not infinite, being about several meters per second. However, for the short segment considered here the assumption of infinite wave velocity is quite Momentum-Y: good. In dimensionless form. with I, being the length of (5) the daughter portion as measured from the apex of the bifurcation centerlines. For the mother portion, the same process is folAt X = 0, P,, = P, + P3 sin (2naB) lowed. resulting in the same field equations, including (mother entrance) (9) the use of h as the characteristic length. At X = 1,/b, P,., = P2 + P4 sin (2na0) 3.2 Boundary and initial conditions (daughter exit) Momentum - X:
The vessel walls are assumed to be impervious, i.e. there is no mass flow through them. In addition, they are assumed to be rigid so that transverse velocity there can be set to zero. Finally, it is assumed that there is no slip at this rigid wall so that the axial velocity can also be set to zero. For the mother portion. pressure and axial velocity are symmetric with respect to the centerline. Transverse velocity, however. is antisymmetric with respect to the centerline, and is zero at the centerline. In dimensionless variables the transverse boundary conditions for the mother and daughter portions are 1. mother portion At Y = 0, g
= 0 and I’ = 0. (centerline)
At Y = a/b, U = 0 and V = 0. (mother wall),
(6)
where a is one-half the width of the mother portion, 1.41 cm. 2. daughter portion At Y = - 1, U = 0 and V = 0. (inner wall) (7) At Y = + 1. U = Oand V = O.(outerwall). The axial boundary conditions used were that of fully-developed parabolic flow at the mother entrance and the daughter exit. Regardless of the initial state of the fluid, it will respond to the periodic forcing function until it is in harmony with it. Analogously, the values of the field variables will move toward a steady periodic condition even if initially they are given zero or any other arbitrary values. Hence 3.At0=O,U=OandV=OatallXandY locations for both mother and daughter portions. (8)
so that AP = P,., - Pi, = P, + P6
sin
(2naf3).
(10)
The values used for P5 and P6 were 1210 and 1100, respectively. The value of P5 was determined by performing the calculations for the steady flow case Reynolds number (a = O), where an upstream (Re = Zau,,/v) of 146 was obtained, and the value of P6 was then set to give the flow a reasonable periodic behavior and Reynolds number variation. The dimensionless frequency, a, was set equal to 1.0 for the periodic flows, where instantaneous upstream Reynolds numbers varied from 130 minimum to 173 maximum, the mean value being different from 146 because of the problem nonlinearity. The quantity Id was determined by performing the calculations for various lengths until fully developed flow was obtained at the daughter outlet. It was found that an Id of 6 cm satisfied this fully developed criteria. 3.4 Numerical procedure The above governing equations were solved numerically using the Marker and Cell (MAC) method, which was developed at the Los Alamos Scientific Laboratory and is the subject of a comprehensive report by Welch et al. (1966). Unlike the stream function-vorticity method, the MAC method involves the primary variables U, V, and P, thus enabling the direct application of the physical boundary conditions. In the present study, the MAC method was extended in that it was necessary to set up two grids, one for the mother and one for the daughter, and to couple the two grids by weaving pressure and velocity values along their line of union. Six forcing cycles were needed before steady periodic conditions were achieved. The computer used was the IBM 360/44. Complete details of the numerical work are given by Fernandez (1972).
3.3 Forcing ,funcrion In order to simulate the pulsatile nature of arterial flow, a sinusoidal periodic forcing function was used. Because the fluid is incompressible and the wall is rigid, the pulse wave will propagate at infinite velocity from the flow entrance at the mother portion to the flow exit at the daughter portion. Hence, there is no phase fag between the pulse at the entrance and the
A INTERPRETATION
OF COMPUTER RESULTS
Although the model used is two-dimensional and has a sharp corner, the results seem to be in fair agreement with Crowe’s (1969) experimental results from branched configurations formed by the intersection of axisymmetrical tubes. Both regions of eddying
578
R. C. FERNANDBZ,K. J. DE WITT and M. R. B~TWIN
at the outer wall and of high shear stress at the inner wall were found in our calculations for both steady and pulsatile flows. Since the Reynolds numbers were not high enough, the eddy was not sucked back into the mother channel during flow deceleration. The position of high shear stress was also unchanged by the pulse. No secondary flow, which Brech and Bellhouse (1973) attribute in their experiments as the reason for preventing the appearance of zones of stasiq could be handled in the present study because of the two-dimensional model used. Figure 2 shows the results for steady’ flow at upstream Reynolds number of 146. The flow is from left to right. The arrows are the total velocity vectors, i.e. the vector sum of transverse velocity, V, and axial velocity, U. The heavy curves are the velocity profiles at the locations indicated by the heavy vertical lines, and are approximately the same as the usual velocity profiles drawn for the axial velocity alone. The approximation is not too much in error since the transverse velocities at the locations considered are nearly zero. Upstream, the profile is parabolic. Downstream the profile shows flow reversal at the outer wall near the comer and also steepness indicating high shear rate at the inner wall. The rectangular portion at the outer wall is also shown magnified four times in area. Here only the axial velocity is plotted so that the envelope shown by the dotted line going through points of reversal may be accurately drawn. Farther downstream the profile shows modulation and a trend towards recovery of the new parabolic shape. For pulsatile flow between the instantaneous Reynolds numbers of 130 to 173, the same general features are evident; however, there is the added feature of growth and disappearance of the envelope as the pressure forcing function goes through a cycle. At a Reynolds number of 130 (25” into the pressure pulse), no envelope can be seen as is illustrated in Fig. 3. Figure 4 shows a small envelope starting to appear at a Reynolds number of 142 (68” into the pressure cycle). As the Reynolds number increases the
It-i
/
Fig. 3. (a) Velocity profiles for pulsatile flow for Re = 130 (25” into pressure pulse). (b) Enlarged view outside wall region.
envelope grows both upstream and downstream and also away from the wall into the mainstream. Growth in the downstream direction is thd most extensive and that towards the mainstream is the least extensive. Figure 5 shows the maximum Reynolds number reached of 173 (158’ into the pressure cycle), but the size of the envelope is not maximum then. Instead, as seen in Fig. 6, it reaches a maximum at a Reynolds number of 169 (202” into the pressure cycle). Afterwards, the envelope contracts and eventually disap pears, having been present over the interval from approx. 60” to 335” of the pressure cycle. These results indicate that, for pulsatile flow, the instantaneous Reynolds number is not sticient for specifying the flow, as the pressure pattern must also be specified. Similar qualitative results using different numerical techniques were obtained by O’Brien et al. (1973). Pressures around the stagnation point on the inner wall were on the order of 1.7 times as great as the pressures at the input boundary. This pressure difference could be of importance in evaluating the impaction mechanism if the values required for damage of vessel wall or blood cells themselves were known, since the difference in pressures is a measure of the
It-i Fig. 2. (a) Velocity profiles for steady flow for Re = 146. (b) Enlarged view of the separation region.
--
---L--f
*
Fig. 4. (a) Velocity profiles for pulsating flow for Re = 145 (68” into pressure pulse). (b) Enlarged view of separation region.
579
Pulsatile flow through a bifurcation
---_) /
Fig. 5. (a) Velocity profiles for pulsatile flow for Re = 173 (158” into pressure pulse). (b) Enlarged view of separation region.
square of the velocity of impaction. In the region of stasis or eddying at the outer wall, the pressure variance among neighboring cells was less than 0.15 dyn/cm’. This fact somewhat questions Texon’s (19-60) suction hypothesis unless the endothelium proves to be extremely pliant. A region of high shear stress was found at the inner wall about a daughter half-width downstream from the stagnation point, but the highest shear stress at a Reynolds number of 146 for the steady case was ca. ‘10 dyn/cm2, which is an order of magnitude smaller than either the experimental value of 107 dyn/cm’ reported by Crowe (1969) for a tube Reynolds number of 100, or the extrapolated experimental value for blood of 220 dyn/cm2 by Brech and Bellhouse (1973) for a Reynolds number of 750. For the pulsatile case at a Reynolds number of 173, the highest wall shear stress was 12 dyn/cm’. No value was reported by Crowe for the pulsatile case, but Brech and Bellhouse’s extrapolated blood value is 310 dyn/cm’. These values, expecially those. from the present study, are all lower than the 380dyn/cm2 for wall tearing reported by Fry (1968) for dogs. These
----w
/
Fig. 6. (a) Velocity profiles for pulsatile flow for Re = 169 (202” into pressure pulse). (b) Enlarged view of separation region.
values Seem to suggest that if direct wall damage were to be a plausible mechanism then the pathological values for damage to humans should be smaller than that for canines. The periodic appearance of the eddy at the outer wall suggests a further investigation into the plausibility of the proposed mechanism on autoregulatory response. It also provides a test on the mechanisms based on mass transfer, such as deposition and suppression. If it can be shown that the disappearance of the eddy allows surface renewal, i.e. material previously trapped in the eddy has a chance to escape allowing influx of fresh material to the outer wall, then these mechanisms will have to be revised to incorporate this added phenomena. S. DISCUSSION The flow behavior of a Newtonian fluid in a twodimensional bifurcation has been investigated numerically in order to aid in the evaluation of current, proposed mechanisms for atherogenesis. The key results for pulsatile flow illustrate that a region of separation exists on the outer wall of the daughter channel over a large portion of the pressure cycle, the possible effects of which have been interpreted with regard to these existing atherogenesis theories. In order to completely evaluate all of the presently proposed mechanisms, further work should include the experimental and numerical study of particle trajectories, the determination of the stress magnitudes necessary to damage vessel walls, platelets, erythrocytes, etc., and the determination of the nature and magnitude of the forces between cells that would cause them to agglutinate or undergo sedimentation.
REFERENCES
Attinger, E. 0. (1964) Flow patterns and vascular geometry. In Pufsarile Blood Flow. (Edited by Attinger, E. O.), p. 179. McGraw-Hill, N.Y. Brech. R. and Bellhouse, B. J. (1973) Flow in branching vessels. Cardiouasc. Res. 7, 593400. Caro, C. G., Fitz-Gerald, J. M. and Schroter. R. C. (1969) Arterial wall shear and distribution of early atheroma in man. Nature 223, 1159-1161. Caro, C. G., F&x-Gerald, J. M. and Schroter, R. C. (1971) Atheroma and arterial wall shear: observation, correlation and proposal of a shear dependent mass transfer mechanism for atherogenesis. Proc. Roy. Sot. Lond. B. 177, 109-159. Crowe, Jr, W. J. (1969) Studies of arterial branching in models using flow birefringence. PhD. Thesis, The University of Florida, Gainesville, Florida. Davids, N. and Cheng, R. (1972) Finite element analysis of flow in blood vessels with arbitrary cross-section. Progress Report, The Pennsylvania State University, University Park, Pennsylvania. Dintenfass, L. (1962) Thixotropy of blood and proneness to thrombus formation. Circ. Res. 11, 233-239. Fernandg R. C. (1972) Pulsatile flow through a bifurcation with applications to arterial disease. PhD. Thesis, The University of Toledo, Toledo, OH.
580
R. C. FERNANDEZ, K. J. DE Wtrr and M. R.
Fox. J. and Hugh, A. E. (1966) Localization of atheroma: a theory based on boundary layer separation. Brit. Heart J. 28, 388-399. Fry, D. L. (1968) Acute vascular endothelial changes associated with increased blood velocity gradients. Circ. Res. 22, 165-197. Lynn, N. S., Fox. V. G. and Ross, L. W. (1972) Computation of fluid-dynamical contributions to atherosclerosis at arterial bifurcations. Biorheology 9, 61-66. Merrill. W. W., Gilliland, E. R., Lee, T. S. and Salzman, E. W. (1966) Blood rheology: effect of fibrinogen deduced by addition. Circ. Res. 13, 437-446. Mitchell, J. R. A. and Schwartz, C. J. (1965) The morphology of arterial plaques. In Arterial Disease, p. 18. Blackwell, Oxford. Mustard, J. F., Rowsell, H. C., Murphy, E. A. and Downie, H. G. (1963) Intimal thrombosis in atherosclerosis. In Euolurion of the Atheroscleroric Plaque (Edited by Jones, R. L.), p. ,181. University of Chicago Press, Chicago, IL. Naff, S. A. (1969) Accelerated laminar flow through a twodimensional bifurcation. M. S. Thesis, Carnegie-Mellon University, Pittsburgh, PA. O’Brien, V., Ehrlich, L. W. and Friedman, M. H. (1973) Nonlinear simulation of unsteady arterial flows in a branch. Presented at the American Institute of Chemical Engineers Annual Meeting, Philadelphia, PA. Spurlock, 1. M., Homa, M. and Durfee, R. L. (1968) Atherosclerosis etiology: further evidence of the role of hemodynamic effects. Research Monograph of the Atlantic Research Corporation, Alexandria, VA.
BOIWIN
Texon, M. (1960) The hemodynamic concept of atherosclerosis. Bull. N.Y. Acad. Med. 36, 163-274. Welch, J. E., Harlow, F. H., Shannon, J. P. and Daly, B. J. (1966) The MAC Method-A computing technique for solving viscous, incompressible, transient fluid-flow problems involving free surfaces. Los Alamos Scientific Laboratory. Los Alamos, NM.
NOMENCLATURE axial and transverse coordinates, cm dimensionless axial and transverse coordinates velocity vector, cm/set axial and transverse velocities, cm/set dimensionless axial and transverse velocities half-widths of the mother and daughter channels, cm length of daughter channel, cm pressure, dyn/cm* dimensionless pressure = pb*/pv* fluid density = 1.05 gm/cm3 fluid kinematic viscosity = 0.0381 cm*/sec time, set dimensionless time = tv/b* forcing frequency, Hertz dimensionless forcing frequency = r&*/v Reynolds number = Zau,,/v average axial velocity, cm/set total branching angle = 90”.
