- Email: [email protected]
Particle paths and stasis in unsteady flow through a bifurcation
I
Biomechonrc.
L97i.
Vol.
IO. pp. 561-568.
Perpamon
Press
Prmted m Great Brnaln
PARTICLE PATHS AND STASIS IN UNSTEADY FLOW THROUGH A BIFURCATION*?LOUISW. EHRLICH and MORTON H. FRIEDMAN Research Center, Applied Physics Laboratory. The Johns Hopkins University, Laurel. MD 20810, U.S.A. Abstract-In
our study of the role of fluid mechanics in atherogenesis, particle paths and stasis in a two-dimensional symmetric branch were examined, using the results of earlier fluid mechanical simulations of pulsatile blood flow to provide the velocity field. In regions where separation occurs. the
particles trace paths which cannot be predicted by a steady flow calculation, but are complex functions of the fluid velocity field. All particles examined were downstream from their starting point at the end of one pulsatile cycle, implying the absence of “perpetual” stasis. The numerical results suggest that stasis may not be as significant a factor in atherogenesis as may have been thought. at least for the area ratio considered here.
INTRODUCTION Among the sites in the vascular tree at which athero-
sclerotic lesions appear most likely to form, a substantial number might be expected from fluid mechanical considerations to be sites of low shear, separation or stasis. Indeed, all three of these phenomena have been proposed as responsible at least in part for the distribution of early vascular lesions (e.g. Gessner, 1973). In Friedman et al. (1975) we used numerical calculations of pulsatile flow in a Y-branch to examine the shear distribution and occurrence of separation along the walls of the branch; here we examine the “distribution of stasis” in the flow. The use of quotation marks in the preceding sentence emphasizes the fact that stasis has not the precise definition that shear and separation enjoy. It is a Lagrangian property, associated with particle trajectories. and cannot be measured by many techniques conventionally employed to probe the Eulerian flow field. This distinction is particularly important for the unsteady pulsatile flows of principal interest here. In steady flow, fluid particles remain on streamlines; if the streamlines are closed, as in a separation bubble, the particles in principle never leave the region. In contrast, when the flow is unsteady, the velocity field changes with time and the particle does not follow any single streamline. The imprecise definition of stasis and the greater difficulty of measuring particle paths in unsteady flow is mirrored in a paucity of published data on this subject, particularly as compared to the associated Eulerian phenomenon of separation. In the cardiovascular literature, Fox and Hugh (1970) used in V~DO radiography to examine qualitatively the stasis of contrast medium in the human internal carotid artery.
* Received 21 January 1977. t Supported in part by Research Grant HL14207 from the National Heart and Lung Institute.
They demonstrated persistence of medium from the carotid origin along the length of the vessel wall for several pulsatile cycles under conditions of elevated intracranial pressure, and in the proximal 1-3 cm of the artery in patients with signs of cerebrovascular disease. The contrast medium also appeared to remain in recesses in the vessel wall. In an accompanying paper, Hugh and Fox (1970) showed that contrast medium persisted immediately distal to an atheromatous plaque at the origin of the internal carotid artery, and noted that the site of this lesion was coincident with that of contrast medium stasis in the undiseased vessel. The authors proposed that thrombotic changes in the static blood could involve the wall to produce mural thrombi predisposing to atheroma a view which has also been put forward by French (1971). Goldsmith’s (1972) model studies of particle dynamics in flowing blood have shown that rouleaux can form in the low shear rate regions which might reasonably be expected to correspond to static zones. In steady model experiments past a spherical obstruction affixed to the tube wall, he found that red cells were trapped in the vortex behind the sphere for times up to a thousandfold longer than that required for nearby cells in the mainstream to pass this region of flow separation. In the studies to be reported here, the flowing blood will be regarded as homogeneous and the calculated particle trajectories will be those of “fluid particles” rather than formed elements. Our understanding of the interaction of hydrodynamic and inertial forces in fluids whose particle density is comparable to that in blood is simply insufficient to attempt the latter. As noted earlier, emphasis will be placed on fluid particle dynamics in unsready flows, where the streamline representation appropriate to the description of particle trajectories in steady flow cannot be used. Some steady flow calculations were done to validate the computational scheme, to demonstrate the important effect of pulsatility, and to determine
561
562
LOUISW. EHRLKHand
the rates at which fluid particles move along the streamlines which define their path. For both steady and unsteady flows, emphasis will be placed on the distances traveled by the particles whose trajectories are followed, since this quantity can be used to define stasis in a quantitative fashion. The unsteady flow field used to calculate the particle trajectories is two-dimensional, three-dimensional unsteady calculations being too expensive to perform. The precise definition of stasis and the computational techniques developed here for evaluating this property can readily be adapted to those cases where three-dimensional flow fields are available. However, care must be taken in extrapolating the results reported here to the case of blood flow in elastic arteries, where the geometry is less regular and three-dimensional phenomena such as secondary flows are to be expected. Similar caution should be exercised in drawing conclusions regarding the role of stasis in atherogenesis, though the implication of this study is that stasis may not be as significant as some others have thought. MATHEMATICAL
METHODS
Governing flow and trajectory equations The model geometry is shown in Fig. 1. The geometry and the flow are two-dimensional in the plane of the figure. The branch and the flow are symmetrical so that calculations only have to be carried out on one side of the indicated plane of symmetry. The flow is fully-developed and parallel at the inlet to the parent vessel. At the sharp “corner”, the outer wall of the parent becomes the outer wall of the daughter vessel; the daughter’s inner wall originates at the flow divider “tip” on the plane of symmetry. The entire branch angle is 90”. The area ratio of the branch, a, is equal to the spacing between the daughter vessel walls divided by L, the half-width of the parent chan-
MORTON
H. FRIEDMAN
nel. For the present studies, CL= 0.796, which is a typical area ratio for aortic bifurcations (e.g. Caro. Fitz-Gerald and Schroter. 1971). The equations which govern the viscous flow in this region are derived from the normalized twodimensional Navier-Stokes equation for a Newtonian incompressible fluid. $&
V2$ = -R,
(la) (lb)
where q = (2L20/v) “’ , !2(%, y, 7) is nondimensional vorticity, differentiation is indicated by subscripting. 7 = cot, Re, = UL/v (note that this value is half the Reynolds number as it is usually defined), I&, J, 7) is the nondimensional stream function, f = x/L, 8 = y/L, w is angular frequency, v is the kinematic viscosity of the fluid, t is time, and U is the mean inlet flow velocity averaged over y and t. In two dimensions, the vorticity and the stream function are related to the nondimensional velocity components i@, j, 7) and 4% 8. 7) by R = r%ld~ - diijir~, where ii = u/U, t? = ii = a+@, and fi = -a$/&, v/U, and u and v are the components of fluid velocity in the x and y directions, respectively. The method of solution of these equations is described in Ehrlich (1974). Using these definitions, the trajectory equations for a fluid particle, dx/dt = u and dyldt = v, can be written in nondimensional terms: d%/dr = (2Re,,,/&i,
(2a)
dJ/dr = (2Re,/n2)8.
(2b)
and
All unsteady-flow calculations are for an idealized pulsatile flux wave form consisting of a constant Re component of 100 on which is superimposed a sinusoidal flow whose Re amplitude is 100 and whose frequency is such that q = 10 in the parent vessel. Thus, the maximum instantaneous flow (at 7 = 0) is twice the mean, the minimum flow (at 7 = n) is zero, and the values of Re, and ~7are 100 and 10, respectively. The resulting fluid velocity field was the starting point for the trajectory calculations presented here. DeJnition
Fig. 1. Model geometry, showing 3 and JJ coordinates. Because of symmetry, calculations only have to be carried out for the half-branch region outlined heavily. See text for definition of symbols.
- Re,(&.R,- - I,@&) = V’R,
of stasis
As noted in the Introduction, a prime object of this study was to examine the levels of stasis in pulsatile flow through the branch. Stasis is ultimately a Lagrangian property, in that it depends on the distances traveled by individual “fluid particles” in the stream; yet it is normally described in an Eulerian fashion in terms of its distribution within the geometry of interest. To reconcile these aspects of stasis, its definition is here based on the distances traveled by particles passing through selected points within the branch. The basic computed quantity is thus ?@r, Ti, A), which is the vector distance traveled in time A by
Stasis in unsteady flow a particle initially at a location r at time ri. In all calculations to be reported here, r was in the daughter vessel and A was set equal to 211; that is. one pulsatile cycle. 2 (r. Ti, 2n) is then essentially the average velocity of a particle through one cycle. Note that W depends on the distance of r from the vessel wall (2 + 0 as the wall is approached) and the peculiarities of the local flow. To isolate the effect of the latter, which is of greater interest, we consider point r’ far downstream in the branch. where the flow is parallel; in this region. .a(~-‘, 5i, 2n) is independent of ri. Let g,,(i) = n(r’, ri. 2~). We define a normalized distance travelled &‘((I;ri) = B(r, Tit 2x)/w,, (t’), where r and r’ are the same distance from the wall of the branch. W(r, ri) goes to unity across the entire inlet and outlet cross sections. Since particles whose stasis is “large” travel little, the value of relative stasis at a point was defined as (a)-‘, where (2) = (1/2n) jiR 9 dri. The integral was approximated by averaging the values of $? for eight equally-spaced values of ri, starting at TV = 0. Because of the normalization by 8,,, (n)-’ is “relative” in the sense that it exceeds unity if the stasis is greater than that which would obtain in fully-developed parallel flow at that distance from the wall. Velociry field and trajectory
calculations
To determine particle paths, it was necessary to have a complete velocity vector field. The computational scheme described above yielded the stream function only at grid points. Thus, these data were splined as described below, using two dimensional cubic splines, and the splines were then differentiated to obtain velocity vectors. When particles crossed the outlet, their trajectories were continued using the analytic solution for fully developed flow in a planeparallel channel. To determine the velocity vectors within the branch, we considered the following. Since the normal derivative of the stream function JI. = 0 along the wall. @ was replaced near the wall by +(a, 9) = II/, + j2 g(% P),
(3)
where $,,, is the value of I/I on the wall and z is the nondimensional perpendicular distance from the wall to the point (X. 7). Thus $. = I,& = 2Zg + Z2gi, where g; is the derivative of g in the I-direction. On the wall, f = 0 and $, = 0. The function g(x‘, p) was determined from
Applying l’Hospital’s Rule to equation (4) and using equation (3). we have at the wall, g&, j) = -Q&Z, ~)/2. Since the spline routine we used required a rectangular array of grid points, the derived values of g were reflected symmetrically across the wall. The function g(f, 3) was splined and differentiated to yield the required velocities. Since the velocity vectors are obtained from an analytic spline on g, hence effec-
563
tively on $I, the scheme is mass-conserving. The particles were advanced by forward-differencing the trajectory equations, (equation (2)), using AT = 2ni880. This was the same time step as that used in the unsteady fluid mechanical calculations in Ehrlich (1974). To establish the validity of the above technique, the problem of equation (1) was solved in the steady state to yield the constant velocity field. A variety of starting points were chosen and particle paths were computed for A = 2n. Since each particle must stay on a streamline in steady flow, the change in stream function along the trajectory was considered a measure of the error in the technique. In all cases tried, the error was less than 0.019;, i.e. at least 4 figure accuracy, and often more. This level of error did not change when AT was halved. RESULTSAND DISCUSSION Separation
area: choice of starting points
When a pulsatile flow having the parameters given in the preceding section is passed through the branch, a region of transient separation develops distal to the corner, along the outer wall of the daughter. lasting a little over 10% of the cycle (Friedman et al., 1975). Figure 2 shows roughly the greatest extent of this region, at maximum flow. Outside of this region the flow was such that most fluid particles in the daughter either were swept out of the bifurcation in much less than one cycle or. if they started close to the wall, traveled nearly parallel to the wall for a distance similar to the parallel-flow value (9 = 1). The more “interesting” particles were those which started near or in the region in Fig. 2, which we shall term the “sep aration area”. Our attention was thus directed to this part of the branch.
1 1375--
11251.
2.0
2.0625
2.125
2.1875
2.25
Fig. 2. Region of separation at maximum flow (T = 0). Grid lines are for AR = Ajj = l/16, the mesh size used for the fluid mechanical calculations. The position of the point of separation is uncertain, owing to the steep gradient of wall vorticity immediately distal to the corner.
564 Efects
LOIJS W. EHRLICH and MORTON H. FRIEDMAN of starting
times and positions
on %’
The distance a particle travels during a 2n cycle is a function of both the starting position and the time during the cycle at which the position is occupied. Figures 3(a) and 3(b) show the termini of particles leaving two positions equidistant from the wall at various times during the cycle. The terminal locations vary considerably and depend on the starting
location and ri. However, some rough generalizations are possible. If a particle originates distal to the separation area. then it travels more-or-less parallel to the wall. The particle paths terminate downstream near the wall. From Fig. 3(a) it appears that particles starting when the flow is approaching its maximum value travel least.
2.25 (a) 200-
1
175-
outer
y 1.50-
1.25-
1.75 2.75
2.00
2.25
2.50
2.75
3m
3.25
3.50
3.75
’
’ 4.00
3.oO
3.25
3.50
3.75
4.00
’
T7
( b)
2.50
2.25
2.00 Outer wall- daughter ci 1.75
1.25
/ 2.00
2.25
2.50
2.75
I
x
Fig. 3. Termini and trajectories of a set of particles. 0: starting point; + : terminus for indicated value of starting time si; A: terminus in steady flow, A = 2x; ---: separation area. Trajectories are given for q = 0 and q = II. (a). Point initially at (2.227, 1.211), outside separation area. The two trajectories are not well resolved. (b). Point initially at (2.102, 1.086), inside separation area; terminus for q = 7 z/4 is (5.36, 4.27).
565
Stasis in unsteady flow
On the other hand, if the particle starts in the separation area. its final position depends upon whether separation occurs while the particle is there, or not (see Fig. 3b). Particles which start nearer zero flow parallel to the wall and do not travel far. However. near maximum flow. the particle travels in a swirling path and can be swept out of the area into the rapid mainstream flow. and carried a considerable distance. To see some of these details more clearly. the region of Fig. 3 has been rotated 4Y and the direction normal to the wall magnified, giving Figs. 4(a) and
* The term “separation bubble” refers to the separated flow itself. while “separation area” refers to that portion of the branch occupied by the bubble at 5 = 0.
0-l
Outer
i’ 7n 7 X
I 7001 1
4(b). Figure 4(a) shows that particles starting distal to the separation area near maximum flow are driven toward the wall by flow around the separation bubble. where they travel more slowly. Otherwise. the particle drifts toward and parallel to the wall, oscillating when the flow reverses but never returning to its starting point. Figure 4(b) shows the paths and termini of two particles starting at a point within the separation area. If the flow is not separated at r,, and the particle can leave the separation area before separation commences, then its path is similar to those in Fig. 4(a) (Ti = rtin Fig. 4b). If. on the other hand, the particle starts while the separation bubble* is present (TV = 0 in Fig. 4b). then its path is more complicated.
wall- daughter
3n T,
3na b-5
-5-J
X -Ti P5n 4T
0 02 -1
Outer
o-
wall-daughter
$ -Tj x
f O.Ol-
.
0.02I 0.03-
71 = 0 ?x
ble. -). Replot of Fig. 3. rotated 45’. Coordinate normal to wall fi) magnified to show paths more clearly; E-coordinate measures distance from corner C parallel to wall. 0: position of particle at intervals of t = n/4.
566
Louts 1.25
W. EHRLKH and
MORTONH. FRIEDMAN
1
1.21875 1.1875 1.15625
1.0525 1.03125
0.96875 1.96875 2.00
2.0625
2.125
2.1875
2.25
x
Fig. 5. Trajectories of a set of particles starting at r = 0 inside the separation bubble. 0 : starting area. Path 5 is the path fi = 0 point ; 0 :position of particle at intervals of n/4; --:separation
on Figs. 3(b) and 4(b). starting position; all proceed downstream. Since a value of A = 2n has been used, this observation implies the absence of “perpetual” stasis in the pulsatile flow case.
The trajectories of particles starting within the sep aration bubble depend on where in the bubble they start. Figure 5 shows the paths of several particles originating in the separation bubble at r = 0. The paths vary considerably and bear no resemblance to the closed streamlines followed by a particle in a steady separation bubble. Paths 1 and 4 are illustrative of the trajectories of a large fraction of the particles originating in the bubble; the circulation in this region causes the particle to move out towards the mainstream which sweeps it away. Particles starting close to the wall (Paths 2 and 5) move more slowly and do not reach the mainstream before the separation bubble dissipates. The particle which follows Path 3 appears to have originated near the vortex center. As with particles originating outside the separation area, no particle examined here returns to its
Relative stasis in steady and pulsatile flow
In Figs. 3 and 4 are also shown the termini which would have been reached by the sample particles had the flow been steady at Re = 100. These termini yield a relative stasis in steady, flow, 9; *, independent of q. The effect of pulsatility on stasis can be assessed by comparing (9) - ’ and 9; ‘; see Table 1. In general, the average stasis of particles in pulsatile flow is greater than that in steady flow at the same mean flow. Two exceptions to this generalization are points inside the separation area, denoted by a dagger in Table 1; particles starting at these points near or at
Table 1. Relaiive stasis in pulsatile and steady flow 0.022 1 (w>-’ 0.0442
0.133 0.22 1 0.309 0.398 0.486 0.575 0.663 0.751 0.840
1.488 1.675
1.192 1.330
1.490
1.284
1.314 1.241 1.142 1.127 1.091 1.070 1.047
1.233 1.190 1.156 1.127 1.103 1.083 1.067
t 1.114t 1.603 1.464 1.280 1.205 1.156 1.100 1.059 1.042
9;’ 1.034 1.361 I .337 1.276 1.220 1.174 1.138 1.109 1.086 1.067
(a)-’
0.01105 a;’
: 1.080t 2.725 2.347 2.105 1.880 -
1.247 1.048 1.429 1.410 1.416 1.366 1.342
* See Fig. 4 caption for definition of coordinates. t For at least one q. particle traveled into mainstream and was swept outside of the computational region, so % could not be found. $ See text.
Stasis in unsteady flow Table 2. Normalized distance traveled near wall Starting point. r _~ M 0.177 0.177 01177 0.177 0.177 0.177
0.0221 0.0110 0.00552 0.00276 0.00138 0 (wall)*
ri = 0 1.177 0.292 0.156 0.161 0.174 0.189
Normalized distance travelled, & 5i = ?I Ti = 3ni’2 K, = ni? 0.756 0.495 0.288 0.208 0.196 0.189
0.826 0.473 0.287 0.210 0.191 0.189
0.321 0.212 0.182 0.194 0.191 0.189
* Normalized distance computed from wall vorticity
r = 0 are swept into the mainstream, travel far, and bias (W)-’ downward. If no particles starting at points inside the separation area experience such large excursions, then the stasis values of these points are particularly high. In all cases, the relative stasis values exceed unity and approach unity distally as the flow becomes more parallel. The non-parallel flow does in general increase stasis at points in the branch, though individual values of d in pulsatile flow may be greater than unity. At a given distance from the wall, the relative stasis values pass through maxima before leveling off towards unity at the outlet. Limiting stasis
as wall is approached
In the Appendix it is shown that as the starting point approaches the wall, the normalized distance traveled in one cycle should approach a vaiue which depends on only the starting position along the wall. This value is the time-average vorticity (proportional to shearing stress) at the wall position divided by the average parallel flow vorticity. To confirm this result, we selected a wall point and a series of positions, a11 within the separation area, on the perpendicular to the wall at this point. The normalized distances traveled by particles starting on the perpendicular at four times during the cycle are displayed in Table 2. It is clear that as the wall is approached, the distance traveled over a cycle is decreasingly dependent on starting time. Furthermore, the values of W do indeed approach the limiting value at the wall. The consistency of the trajectory calculations with the prediction based on wall vorticity lends credence to our technique for calculating stasis. However, it should be noted that the smallest values of 5 used in this numerical verification correspond to distances from the wall at which the continuum assumption underlying the analysis is no longer realistic for even large arteries.
complex functions of the time dependent velocity field. Further, the average distance traveled by particles in pulsatile flow is generally less than the corresponding distance in steady flow. Since particles. in one cycle, always move downstream from where they started (see Table l), the implication is that no “perpetual” stasis occurs in the region of bifurcations. It would seem that unless the relative stasis in arteries is substantially greater than unity, the role of this phenomenon in atherogenesis is questionable. The highest (Se> -i which we have found is less than three, which does not seem very large. Recalling the limitations noted in the Introduction, we hesitantly conclude that stasis per se may not be a major factor in the localization of atherosclerotic lesions, at least for area ratios near that considered here.
REFERENCES
Care. C. G.. Fitz-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. Land. B 177. 109-l 59. Ehrlich, L. W. (1974) Digital simulation of periodic fluid flow in a bifurcation. Comput. Fluids 2, 237-247. French, J. E. (1971) Atherogenesis and thrombosis. Semin. Hemat. 8, 84-94.
Friedman, M. H., O’Brien, V. and Ehrlich, L. W. (1975) Calculations of pulsatile flow through a branch: Imphcations for the hemodynamics of atherogenesis. Circulation Res. 36, 277-285. Fox, J. A. and Hugh, A. E. (1970) Static zones in the internal carotid artery: correlation with boundary layer separation and stasis in model flows. hit. J. Radial. 43, 37&376.
Gessner. F. B. (1973) Hemodynamic theories of atherogenesis. Circulation Rex 33, 259-266. Goldsmith, H. L. (1972) The flow of model particles and blood cells and its relation to thrombogenesis. Prog. Hemostasis Thromh. 1, 97-139. Hugh, A. E. and Fox. J. A. (1970) The precise localisation of atheroma and its association with stasis at the origin of the internal carotid artery-a radiographic investigation. &it. J. Radio!. 43, 377-383.
CONCLUSIONS
Unlike steady flow, where the particle paths and streamlines coincide, the particle paths in unsteady pulsatile flows are highly dependent on both location and starting time, especially in the separation area. This dependency is evidenced by the paths shown in Fig. 5. In general, the shapes of these paths are
APPENDIXNORMALIZED TRAVEL
NEAR WALL
The particle trajectories in the daughter near the wall are more easily examined in the rotated coordinates i and %, where t is perpendicular distance from the wall and i? is measured parallel to the wall. Let the velocity components in these directions be c(E, f. T) and c(E, 2, T), respectively.
Louts W. EHRLICHand MORTON H. FRIEDMAX
568
Expanding < in a Taylor series at the wall. <(S. 5. T) = <(ii’,0, r) + z &(k, 0, T) + P(z’). By the non-slip boundary condition at the wall, @9,0, t) = 0. Furthermore, &(iG,0, r) in the daughter corresponds to -u+ x, 1, 7) (minus sign because p directed oppositely to 2; see Fig. 1) in the parent. These velocity gradients equal the wall vorticity R,. by equation (lb), since I,& (in the parent) and I(/= (in the daughter) are zero at the wall because u (in the parent) and c (in the daughter) are zero there. The flow is incompressible; hence 5, = -cr. This equation is integrated from c = 0 at i = 0:
S-history is 1
i+K
;-(5)=to
2
(A3)
Substitution of f(7) from equation (A3) into (A2a) yields an expression for dC/dr in terms of only i? and r; integrating from time 7i to ri + 2~.
i(i+,,,T) = -
=-
; _as2,
22 zx
+
=--
1
O(r’) d?
f[0
an, _ 2 -p
(A4)
7) + W3),
(Al)
where i in the integrals is understood to be a dummy variable. Consider that 2 is sufficiently close to zero to permit r and c to be described by the leading non-zero terms in their respective expansions. The trajectory equations [cf. equations (2)3 are then
For starting points suthciently close to the wail, the first term in the brackets in the r-integral will be much greater than the second term. Dropping the second term (the r’-integral), 5 - R, becomes proportional to fe. Thus, if ie is small enough, W remains sufficiently close to 5, so that Q,fi, T) in the r-integral can be approximated by t-U=,, 5): qri + 2x) - wa % Kie
r,+Zll R,(ii’,. r)dr. s r,
(A5)
From Equation (A3). it is seen that neglect of the r’-integral implies 4~) = Z,; that is, the flow becomes parallel (A2a) to the wall as the wall is approached. Hence the left-hand side of Equation (A5) approaches s({i&,. r,), TV,2~). since and all motion is then in the %-direction. The right-hand side of Equation (A5) equals 2x K &(Q,>(F,), where (Q,,)(iQ (A2b) is time-average wall vorticity at the wall coordinate G,. Similarly, in parallel flow far downstream, ~,,((cc, I,}) approaches 2n K io(QJ ,,, where (82,) ,, is time-average where K = 2Re,J$. Denote by the subscript zero the starting coordinates (ii+,,&,)of a particle at starting time TV. wall vorticity in parallel Row. Hence the normalized distance traveled in one cycle as & +O equals (Q,)(&)/ Rearranging and integrating equation (A2b), an expression for the i-coordinate of the particle, in terms of its G-U /,7independent of ~~and a function of only &,. dw -& = K< = Kzn,[w(r),r],
Biomechonrc.
L97i.
Vol.
IO. pp. 561-568.
Perpamon
Press
Prmted m Great Brnaln
PARTICLE PATHS AND STASIS IN UNSTEADY FLOW THROUGH A BIFURCATION*?LOUISW. EHRLICH and MORTON H. FRIEDMAN Research Center, Applied Physics Laboratory. The Johns Hopkins University, Laurel. MD 20810, U.S.A. Abstract-In
our study of the role of fluid mechanics in atherogenesis, particle paths and stasis in a two-dimensional symmetric branch were examined, using the results of earlier fluid mechanical simulations of pulsatile blood flow to provide the velocity field. In regions where separation occurs. the
particles trace paths which cannot be predicted by a steady flow calculation, but are complex functions of the fluid velocity field. All particles examined were downstream from their starting point at the end of one pulsatile cycle, implying the absence of “perpetual” stasis. The numerical results suggest that stasis may not be as significant a factor in atherogenesis as may have been thought. at least for the area ratio considered here.
INTRODUCTION Among the sites in the vascular tree at which athero-
sclerotic lesions appear most likely to form, a substantial number might be expected from fluid mechanical considerations to be sites of low shear, separation or stasis. Indeed, all three of these phenomena have been proposed as responsible at least in part for the distribution of early vascular lesions (e.g. Gessner, 1973). In Friedman et al. (1975) we used numerical calculations of pulsatile flow in a Y-branch to examine the shear distribution and occurrence of separation along the walls of the branch; here we examine the “distribution of stasis” in the flow. The use of quotation marks in the preceding sentence emphasizes the fact that stasis has not the precise definition that shear and separation enjoy. It is a Lagrangian property, associated with particle trajectories. and cannot be measured by many techniques conventionally employed to probe the Eulerian flow field. This distinction is particularly important for the unsteady pulsatile flows of principal interest here. In steady flow, fluid particles remain on streamlines; if the streamlines are closed, as in a separation bubble, the particles in principle never leave the region. In contrast, when the flow is unsteady, the velocity field changes with time and the particle does not follow any single streamline. The imprecise definition of stasis and the greater difficulty of measuring particle paths in unsteady flow is mirrored in a paucity of published data on this subject, particularly as compared to the associated Eulerian phenomenon of separation. In the cardiovascular literature, Fox and Hugh (1970) used in V~DO radiography to examine qualitatively the stasis of contrast medium in the human internal carotid artery.
* Received 21 January 1977. t Supported in part by Research Grant HL14207 from the National Heart and Lung Institute.
They demonstrated persistence of medium from the carotid origin along the length of the vessel wall for several pulsatile cycles under conditions of elevated intracranial pressure, and in the proximal 1-3 cm of the artery in patients with signs of cerebrovascular disease. The contrast medium also appeared to remain in recesses in the vessel wall. In an accompanying paper, Hugh and Fox (1970) showed that contrast medium persisted immediately distal to an atheromatous plaque at the origin of the internal carotid artery, and noted that the site of this lesion was coincident with that of contrast medium stasis in the undiseased vessel. The authors proposed that thrombotic changes in the static blood could involve the wall to produce mural thrombi predisposing to atheroma a view which has also been put forward by French (1971). Goldsmith’s (1972) model studies of particle dynamics in flowing blood have shown that rouleaux can form in the low shear rate regions which might reasonably be expected to correspond to static zones. In steady model experiments past a spherical obstruction affixed to the tube wall, he found that red cells were trapped in the vortex behind the sphere for times up to a thousandfold longer than that required for nearby cells in the mainstream to pass this region of flow separation. In the studies to be reported here, the flowing blood will be regarded as homogeneous and the calculated particle trajectories will be those of “fluid particles” rather than formed elements. Our understanding of the interaction of hydrodynamic and inertial forces in fluids whose particle density is comparable to that in blood is simply insufficient to attempt the latter. As noted earlier, emphasis will be placed on fluid particle dynamics in unsready flows, where the streamline representation appropriate to the description of particle trajectories in steady flow cannot be used. Some steady flow calculations were done to validate the computational scheme, to demonstrate the important effect of pulsatility, and to determine
561
562
LOUISW. EHRLKHand
the rates at which fluid particles move along the streamlines which define their path. For both steady and unsteady flows, emphasis will be placed on the distances traveled by the particles whose trajectories are followed, since this quantity can be used to define stasis in a quantitative fashion. The unsteady flow field used to calculate the particle trajectories is two-dimensional, three-dimensional unsteady calculations being too expensive to perform. The precise definition of stasis and the computational techniques developed here for evaluating this property can readily be adapted to those cases where three-dimensional flow fields are available. However, care must be taken in extrapolating the results reported here to the case of blood flow in elastic arteries, where the geometry is less regular and three-dimensional phenomena such as secondary flows are to be expected. Similar caution should be exercised in drawing conclusions regarding the role of stasis in atherogenesis, though the implication of this study is that stasis may not be as significant as some others have thought. MATHEMATICAL
METHODS
Governing flow and trajectory equations The model geometry is shown in Fig. 1. The geometry and the flow are two-dimensional in the plane of the figure. The branch and the flow are symmetrical so that calculations only have to be carried out on one side of the indicated plane of symmetry. The flow is fully-developed and parallel at the inlet to the parent vessel. At the sharp “corner”, the outer wall of the parent becomes the outer wall of the daughter vessel; the daughter’s inner wall originates at the flow divider “tip” on the plane of symmetry. The entire branch angle is 90”. The area ratio of the branch, a, is equal to the spacing between the daughter vessel walls divided by L, the half-width of the parent chan-
MORTON
H. FRIEDMAN
nel. For the present studies, CL= 0.796, which is a typical area ratio for aortic bifurcations (e.g. Caro. Fitz-Gerald and Schroter. 1971). The equations which govern the viscous flow in this region are derived from the normalized twodimensional Navier-Stokes equation for a Newtonian incompressible fluid. $&
V2$ = -R,
(la) (lb)
where q = (2L20/v) “’ , !2(%, y, 7) is nondimensional vorticity, differentiation is indicated by subscripting. 7 = cot, Re, = UL/v (note that this value is half the Reynolds number as it is usually defined), I&, J, 7) is the nondimensional stream function, f = x/L, 8 = y/L, w is angular frequency, v is the kinematic viscosity of the fluid, t is time, and U is the mean inlet flow velocity averaged over y and t. In two dimensions, the vorticity and the stream function are related to the nondimensional velocity components i@, j, 7) and 4% 8. 7) by R = r%ld~ - diijir~, where ii = u/U, t? = ii = a+@, and fi = -a$/&, v/U, and u and v are the components of fluid velocity in the x and y directions, respectively. The method of solution of these equations is described in Ehrlich (1974). Using these definitions, the trajectory equations for a fluid particle, dx/dt = u and dyldt = v, can be written in nondimensional terms: d%/dr = (2Re,,,/&i,
(2a)
dJ/dr = (2Re,/n2)8.
(2b)
and
All unsteady-flow calculations are for an idealized pulsatile flux wave form consisting of a constant Re component of 100 on which is superimposed a sinusoidal flow whose Re amplitude is 100 and whose frequency is such that q = 10 in the parent vessel. Thus, the maximum instantaneous flow (at 7 = 0) is twice the mean, the minimum flow (at 7 = n) is zero, and the values of Re, and ~7are 100 and 10, respectively. The resulting fluid velocity field was the starting point for the trajectory calculations presented here. DeJnition
Fig. 1. Model geometry, showing 3 and JJ coordinates. Because of symmetry, calculations only have to be carried out for the half-branch region outlined heavily. See text for definition of symbols.
- Re,(&.R,- - I,@&) = V’R,
of stasis
As noted in the Introduction, a prime object of this study was to examine the levels of stasis in pulsatile flow through the branch. Stasis is ultimately a Lagrangian property, in that it depends on the distances traveled by individual “fluid particles” in the stream; yet it is normally described in an Eulerian fashion in terms of its distribution within the geometry of interest. To reconcile these aspects of stasis, its definition is here based on the distances traveled by particles passing through selected points within the branch. The basic computed quantity is thus ?@r, Ti, A), which is the vector distance traveled in time A by
Stasis in unsteady flow a particle initially at a location r at time ri. In all calculations to be reported here, r was in the daughter vessel and A was set equal to 211; that is. one pulsatile cycle. 2 (r. Ti, 2n) is then essentially the average velocity of a particle through one cycle. Note that W depends on the distance of r from the vessel wall (2 + 0 as the wall is approached) and the peculiarities of the local flow. To isolate the effect of the latter, which is of greater interest, we consider point r’ far downstream in the branch. where the flow is parallel; in this region. .a(~-‘, 5i, 2n) is independent of ri. Let g,,(i) = n(r’, ri. 2~). We define a normalized distance travelled &‘((I;ri) = B(r, Tit 2x)/w,, (t’), where r and r’ are the same distance from the wall of the branch. W(r, ri) goes to unity across the entire inlet and outlet cross sections. Since particles whose stasis is “large” travel little, the value of relative stasis at a point was defined as (a)-‘, where (2) = (1/2n) jiR 9 dri. The integral was approximated by averaging the values of $? for eight equally-spaced values of ri, starting at TV = 0. Because of the normalization by 8,,, (n)-’ is “relative” in the sense that it exceeds unity if the stasis is greater than that which would obtain in fully-developed parallel flow at that distance from the wall. Velociry field and trajectory
calculations
To determine particle paths, it was necessary to have a complete velocity vector field. The computational scheme described above yielded the stream function only at grid points. Thus, these data were splined as described below, using two dimensional cubic splines, and the splines were then differentiated to obtain velocity vectors. When particles crossed the outlet, their trajectories were continued using the analytic solution for fully developed flow in a planeparallel channel. To determine the velocity vectors within the branch, we considered the following. Since the normal derivative of the stream function JI. = 0 along the wall. @ was replaced near the wall by +(a, 9) = II/, + j2 g(% P),
(3)
where $,,, is the value of I/I on the wall and z is the nondimensional perpendicular distance from the wall to the point (X. 7). Thus $. = I,& = 2Zg + Z2gi, where g; is the derivative of g in the I-direction. On the wall, f = 0 and $, = 0. The function g(x‘, p) was determined from
Applying l’Hospital’s Rule to equation (4) and using equation (3). we have at the wall, g&, j) = -Q&Z, ~)/2. Since the spline routine we used required a rectangular array of grid points, the derived values of g were reflected symmetrically across the wall. The function g(f, 3) was splined and differentiated to yield the required velocities. Since the velocity vectors are obtained from an analytic spline on g, hence effec-
563
tively on $I, the scheme is mass-conserving. The particles were advanced by forward-differencing the trajectory equations, (equation (2)), using AT = 2ni880. This was the same time step as that used in the unsteady fluid mechanical calculations in Ehrlich (1974). To establish the validity of the above technique, the problem of equation (1) was solved in the steady state to yield the constant velocity field. A variety of starting points were chosen and particle paths were computed for A = 2n. Since each particle must stay on a streamline in steady flow, the change in stream function along the trajectory was considered a measure of the error in the technique. In all cases tried, the error was less than 0.019;, i.e. at least 4 figure accuracy, and often more. This level of error did not change when AT was halved. RESULTSAND DISCUSSION Separation
area: choice of starting points
When a pulsatile flow having the parameters given in the preceding section is passed through the branch, a region of transient separation develops distal to the corner, along the outer wall of the daughter. lasting a little over 10% of the cycle (Friedman et al., 1975). Figure 2 shows roughly the greatest extent of this region, at maximum flow. Outside of this region the flow was such that most fluid particles in the daughter either were swept out of the bifurcation in much less than one cycle or. if they started close to the wall, traveled nearly parallel to the wall for a distance similar to the parallel-flow value (9 = 1). The more “interesting” particles were those which started near or in the region in Fig. 2, which we shall term the “sep aration area”. Our attention was thus directed to this part of the branch.
1 1375--
11251.
2.0
2.0625
2.125
2.1875
2.25
Fig. 2. Region of separation at maximum flow (T = 0). Grid lines are for AR = Ajj = l/16, the mesh size used for the fluid mechanical calculations. The position of the point of separation is uncertain, owing to the steep gradient of wall vorticity immediately distal to the corner.
564 Efects
LOIJS W. EHRLICH and MORTON H. FRIEDMAN of starting
times and positions
on %’
The distance a particle travels during a 2n cycle is a function of both the starting position and the time during the cycle at which the position is occupied. Figures 3(a) and 3(b) show the termini of particles leaving two positions equidistant from the wall at various times during the cycle. The terminal locations vary considerably and depend on the starting
location and ri. However, some rough generalizations are possible. If a particle originates distal to the separation area. then it travels more-or-less parallel to the wall. The particle paths terminate downstream near the wall. From Fig. 3(a) it appears that particles starting when the flow is approaching its maximum value travel least.
2.25 (a) 200-
1
175-
outer
y 1.50-
1.25-
1.75 2.75
2.00
2.25
2.50
2.75
3m
3.25
3.50
3.75
’
’ 4.00
3.oO
3.25
3.50
3.75
4.00
’
T7
( b)
2.50
2.25
2.00 Outer wall- daughter ci 1.75
1.25
/ 2.00
2.25
2.50
2.75
I
x
Fig. 3. Termini and trajectories of a set of particles. 0: starting point; + : terminus for indicated value of starting time si; A: terminus in steady flow, A = 2x; ---: separation area. Trajectories are given for q = 0 and q = II. (a). Point initially at (2.227, 1.211), outside separation area. The two trajectories are not well resolved. (b). Point initially at (2.102, 1.086), inside separation area; terminus for q = 7 z/4 is (5.36, 4.27).
565
Stasis in unsteady flow
On the other hand, if the particle starts in the separation area. its final position depends upon whether separation occurs while the particle is there, or not (see Fig. 3b). Particles which start nearer zero flow parallel to the wall and do not travel far. However. near maximum flow. the particle travels in a swirling path and can be swept out of the area into the rapid mainstream flow. and carried a considerable distance. To see some of these details more clearly. the region of Fig. 3 has been rotated 4Y and the direction normal to the wall magnified, giving Figs. 4(a) and
* The term “separation bubble” refers to the separated flow itself. while “separation area” refers to that portion of the branch occupied by the bubble at 5 = 0.
0-l
Outer
i’ 7n 7 X
I 7001 1
4(b). Figure 4(a) shows that particles starting distal to the separation area near maximum flow are driven toward the wall by flow around the separation bubble. where they travel more slowly. Otherwise. the particle drifts toward and parallel to the wall, oscillating when the flow reverses but never returning to its starting point. Figure 4(b) shows the paths and termini of two particles starting at a point within the separation area. If the flow is not separated at r,, and the particle can leave the separation area before separation commences, then its path is similar to those in Fig. 4(a) (Ti = rtin Fig. 4b). If. on the other hand, the particle starts while the separation bubble* is present (TV = 0 in Fig. 4b). then its path is more complicated.
wall- daughter
3n T,
3na b-5
-5-J
X -Ti P5n 4T
0 02 -1
Outer
o-
wall-daughter
$ -Tj x
f O.Ol-
.
0.02I 0.03-
71 = 0 ?x
ble. -). Replot of Fig. 3. rotated 45’. Coordinate normal to wall fi) magnified to show paths more clearly; E-coordinate measures distance from corner C parallel to wall. 0: position of particle at intervals of t = n/4.
566
Louts 1.25
W. EHRLKH and
MORTONH. FRIEDMAN
1
1.21875 1.1875 1.15625
1.0525 1.03125
0.96875 1.96875 2.00
2.0625
2.125
2.1875
2.25
x
Fig. 5. Trajectories of a set of particles starting at r = 0 inside the separation bubble. 0 : starting area. Path 5 is the path fi = 0 point ; 0 :position of particle at intervals of n/4; --:separation
on Figs. 3(b) and 4(b). starting position; all proceed downstream. Since a value of A = 2n has been used, this observation implies the absence of “perpetual” stasis in the pulsatile flow case.
The trajectories of particles starting within the sep aration bubble depend on where in the bubble they start. Figure 5 shows the paths of several particles originating in the separation bubble at r = 0. The paths vary considerably and bear no resemblance to the closed streamlines followed by a particle in a steady separation bubble. Paths 1 and 4 are illustrative of the trajectories of a large fraction of the particles originating in the bubble; the circulation in this region causes the particle to move out towards the mainstream which sweeps it away. Particles starting close to the wall (Paths 2 and 5) move more slowly and do not reach the mainstream before the separation bubble dissipates. The particle which follows Path 3 appears to have originated near the vortex center. As with particles originating outside the separation area, no particle examined here returns to its
Relative stasis in steady and pulsatile flow
In Figs. 3 and 4 are also shown the termini which would have been reached by the sample particles had the flow been steady at Re = 100. These termini yield a relative stasis in steady, flow, 9; *, independent of q. The effect of pulsatility on stasis can be assessed by comparing (9) - ’ and 9; ‘; see Table 1. In general, the average stasis of particles in pulsatile flow is greater than that in steady flow at the same mean flow. Two exceptions to this generalization are points inside the separation area, denoted by a dagger in Table 1; particles starting at these points near or at
Table 1. Relaiive stasis in pulsatile and steady flow 0.022 1 (w>-’ 0.0442
0.133 0.22 1 0.309 0.398 0.486 0.575 0.663 0.751 0.840
1.488 1.675
1.192 1.330
1.490
1.284
1.314 1.241 1.142 1.127 1.091 1.070 1.047
1.233 1.190 1.156 1.127 1.103 1.083 1.067
t 1.114t 1.603 1.464 1.280 1.205 1.156 1.100 1.059 1.042
9;’ 1.034 1.361 I .337 1.276 1.220 1.174 1.138 1.109 1.086 1.067
(a)-’
0.01105 a;’
: 1.080t 2.725 2.347 2.105 1.880 -
1.247 1.048 1.429 1.410 1.416 1.366 1.342
* See Fig. 4 caption for definition of coordinates. t For at least one q. particle traveled into mainstream and was swept outside of the computational region, so % could not be found. $ See text.
Stasis in unsteady flow Table 2. Normalized distance traveled near wall Starting point. r _~ M 0.177 0.177 01177 0.177 0.177 0.177
0.0221 0.0110 0.00552 0.00276 0.00138 0 (wall)*
ri = 0 1.177 0.292 0.156 0.161 0.174 0.189
Normalized distance travelled, & 5i = ?I Ti = 3ni’2 K, = ni? 0.756 0.495 0.288 0.208 0.196 0.189
0.826 0.473 0.287 0.210 0.191 0.189
0.321 0.212 0.182 0.194 0.191 0.189
* Normalized distance computed from wall vorticity
r = 0 are swept into the mainstream, travel far, and bias (W)-’ downward. If no particles starting at points inside the separation area experience such large excursions, then the stasis values of these points are particularly high. In all cases, the relative stasis values exceed unity and approach unity distally as the flow becomes more parallel. The non-parallel flow does in general increase stasis at points in the branch, though individual values of d in pulsatile flow may be greater than unity. At a given distance from the wall, the relative stasis values pass through maxima before leveling off towards unity at the outlet. Limiting stasis
as wall is approached
In the Appendix it is shown that as the starting point approaches the wall, the normalized distance traveled in one cycle should approach a vaiue which depends on only the starting position along the wall. This value is the time-average vorticity (proportional to shearing stress) at the wall position divided by the average parallel flow vorticity. To confirm this result, we selected a wall point and a series of positions, a11 within the separation area, on the perpendicular to the wall at this point. The normalized distances traveled by particles starting on the perpendicular at four times during the cycle are displayed in Table 2. It is clear that as the wall is approached, the distance traveled over a cycle is decreasingly dependent on starting time. Furthermore, the values of W do indeed approach the limiting value at the wall. The consistency of the trajectory calculations with the prediction based on wall vorticity lends credence to our technique for calculating stasis. However, it should be noted that the smallest values of 5 used in this numerical verification correspond to distances from the wall at which the continuum assumption underlying the analysis is no longer realistic for even large arteries.
complex functions of the time dependent velocity field. Further, the average distance traveled by particles in pulsatile flow is generally less than the corresponding distance in steady flow. Since particles. in one cycle, always move downstream from where they started (see Table l), the implication is that no “perpetual” stasis occurs in the region of bifurcations. It would seem that unless the relative stasis in arteries is substantially greater than unity, the role of this phenomenon in atherogenesis is questionable. The highest (Se> -i which we have found is less than three, which does not seem very large. Recalling the limitations noted in the Introduction, we hesitantly conclude that stasis per se may not be a major factor in the localization of atherosclerotic lesions, at least for area ratios near that considered here.
REFERENCES
Care. C. G.. Fitz-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. Land. B 177. 109-l 59. Ehrlich, L. W. (1974) Digital simulation of periodic fluid flow in a bifurcation. Comput. Fluids 2, 237-247. French, J. E. (1971) Atherogenesis and thrombosis. Semin. Hemat. 8, 84-94.
Friedman, M. H., O’Brien, V. and Ehrlich, L. W. (1975) Calculations of pulsatile flow through a branch: Imphcations for the hemodynamics of atherogenesis. Circulation Res. 36, 277-285. Fox, J. A. and Hugh, A. E. (1970) Static zones in the internal carotid artery: correlation with boundary layer separation and stasis in model flows. hit. J. Radial. 43, 37&376.
Gessner. F. B. (1973) Hemodynamic theories of atherogenesis. Circulation Rex 33, 259-266. Goldsmith, H. L. (1972) The flow of model particles and blood cells and its relation to thrombogenesis. Prog. Hemostasis Thromh. 1, 97-139. Hugh, A. E. and Fox. J. A. (1970) The precise localisation of atheroma and its association with stasis at the origin of the internal carotid artery-a radiographic investigation. &it. J. Radio!. 43, 377-383.
CONCLUSIONS
Unlike steady flow, where the particle paths and streamlines coincide, the particle paths in unsteady pulsatile flows are highly dependent on both location and starting time, especially in the separation area. This dependency is evidenced by the paths shown in Fig. 5. In general, the shapes of these paths are
APPENDIXNORMALIZED TRAVEL
NEAR WALL
The particle trajectories in the daughter near the wall are more easily examined in the rotated coordinates i and %, where t is perpendicular distance from the wall and i? is measured parallel to the wall. Let the velocity components in these directions be c(E, f. T) and c(E, 2, T), respectively.
Louts W. EHRLICHand MORTON H. FRIEDMAX
568
Expanding < in a Taylor series at the wall. <(S. 5. T) = <(ii’,0, r) + z &(k, 0, T) + P(z’). By the non-slip boundary condition at the wall, @9,0, t) = 0. Furthermore, &(iG,0, r) in the daughter corresponds to -u+ x, 1, 7) (minus sign because p directed oppositely to 2; see Fig. 1) in the parent. These velocity gradients equal the wall vorticity R,. by equation (lb), since I,& (in the parent) and I(/= (in the daughter) are zero at the wall because u (in the parent) and c (in the daughter) are zero there. The flow is incompressible; hence 5, = -cr. This equation is integrated from c = 0 at i = 0:
S-history is 1
i+K
;-(5)=to
2
(A3)
Substitution of f(7) from equation (A3) into (A2a) yields an expression for dC/dr in terms of only i? and r; integrating from time 7i to ri + 2~.
i(i+,,,T) = -
=-
; _as2,
22 zx
+
=--
1
O(r’) d?
f[0
an, _ 2 -p
(A4)
7) + W3),
(Al)
where i in the integrals is understood to be a dummy variable. Consider that 2 is sufficiently close to zero to permit r and c to be described by the leading non-zero terms in their respective expansions. The trajectory equations [cf. equations (2)3 are then
For starting points suthciently close to the wail, the first term in the brackets in the r-integral will be much greater than the second term. Dropping the second term (the r’-integral), 5 - R, becomes proportional to fe. Thus, if ie is small enough, W remains sufficiently close to 5, so that Q,fi, T) in the r-integral can be approximated by t-U=,, 5): qri + 2x) - wa % Kie
r,+Zll R,(ii’,. r)dr. s r,
(A5)
From Equation (A3). it is seen that neglect of the r’-integral implies 4~) = Z,; that is, the flow becomes parallel (A2a) to the wall as the wall is approached. Hence the left-hand side of Equation (A5) approaches s({i&,. r,), TV,2~). since and all motion is then in the %-direction. The right-hand side of Equation (A5) equals 2x K &(Q,>(F,), where (Q,,)(iQ (A2b) is time-average wall vorticity at the wall coordinate G,. Similarly, in parallel flow far downstream, ~,,((cc, I,}) approaches 2n K io(QJ ,,, where (82,) ,, is time-average where K = 2Re,J$. Denote by the subscript zero the starting coordinates (ii+,,&,)of a particle at starting time TV. wall vorticity in parallel Row. Hence the normalized distance traveled in one cycle as & +O equals (Q,)(&)/ Rearranging and integrating equation (A2b), an expression for the i-coordinate of the particle, in terms of its G-U /,7independent of ~~and a function of only &,. dw -& = K< = Kzn,[w(r),r],