1 Introduction

PULSATILEflow through constricted vessels has been widely investigated, primarily to better understand the fluid mechanics of arterial stenoses. Owing to the difficulties involved in making accurate measurements of velocities and shear rates in vivo, there is a strong reliance on model studies employing both experimental and numerical techniques. The arterial blood flow is a critical determinant, for example, in mass transport within the blood and between the blood and the vascular wall (CARO et al., 1971), structural damage by hydrodynamic forces (PATEL and VAISHNAV, 1980) and endothelial-mediated control of vascular tone (LANGILLEand O'DONNELL,1986; BASSENGE and Buss~, 1988). Furthermore, in the presence of stenoses, hydrodynamic forces are involved in plaque fissure (GLAGOV et al., 1988; Ku et al., 1990), post-stenotic dilatation (OJHAet al., 1990) and thrombosis and damage of the arterial wall in the post-stenotic region. Consequently, there is a need to better understand arterial blood flow under physiological and pathophysiological conditions. This would also allow more precise clinical detection of mild to moderate stenoses using techniques such as Doppler ultrasound (RITTGERS and SHU, 1990) and magnetic resonance imaging (Ku et al., 1990). In general, for flow through constricted vessels, a high-velocity jet and the associated flow separation are normally seen. With upstream laminar flow, the velocity field in the post-stenotic regions can become transitional or even turbulent depending on the flow conditions and Correspondence should be addressed to Dr Dvinsky at address 1. First received 30 June and in final form 22 December 1992 9 IFMBE: 1994

138

severity of the stenoses (AHMED and GIODENS, 1984; OJHA et al., 1989; 1990; LEIBER and GIDDENS, 1990). We previously (OJHA et al., 1989) investigated the effects of the severity and shape of stenoses using the photochromic tracer technique. With axisymmetric stenoses of 65 and 7 5 % area reductions, they clearly demonstrated that transition to turbulence occurred through the rolling up of the high shear layer between the separated flow and the jet. In addition, they provided detailed information on the spatial and temporal variation of the wall shear stress by measuring its instantaneous values. With milder stenoses, the formation of coherent flow structures such as vortices was seen during the deceleration phase of the flow cycle and in the vicinity of the reattachment point. In assessing the effect of the shape of the stenosis, the main thrust had been on the analysis of axisymmetric stenosis, even though arterial stenoses are normally asymmetric. The few studies that did consider the effects of asymmetry clearly showed that such a factor played a key role in establishing the poststenotic velocity field (YOUNG and TSAI, 1973; OJHA et al., 1989, 1990). The main objective of this study was to numerically simulate pulsatile flow through a mild asymmetric stenosis. The shape of the stenosis and the flow conditions were similar to those used previously (OJHA et al., 1989) in experimental work using the photochromic tracer technique. This was done to obtain a validation of the numerical code in this traditional flow case and to complement the previous experimental results. 2 Methods

In this study, we used a computer code HEMO to simulate blood flow. H E M O was specifically designed for

Medical & Biological Engineering & Computing

March 1994

Fig. 1 Computer reconstruction of the vessel geometry; the plane intersects the vessel at Z = 0"7

calculating pulsatile flows in vascular geometries and provides all necessary tools to accomplish this efficiently and effectively. The code has a user interface for input of the vessel geometry and for specification of the problem (boundary conditions, physical properties etc.), the governing equations solver, and finally a post-processor for graphic and alphanumeric display and hard-copy output of the calculated results. H E M O is based on a new solution algorithm and a new equation discretisation technique developed previously (DvINSKV and DUKOWlCZ, 1993; DUKOWICZ and DVINSKY, 1992). These techniques are briefly outlined below. The algorithm used in H E M O is based on a high-order approximate factorisation of the incompressible NavierStokes equations (DuKowlcz and DVINSKY, 1992). The new method has several important advantages over the standard first-order Chorin method (CHOR~N, 1968) for computing time-dependent flows. As the new algorithm re.quires about the same computational work to compute one time step as the Chorin method, considerable computing time savings usually result. In addition, the new method does not require any 'numerical' boundary conditions, in contrast to most previous methods, which need a normal pressure gradient at the boundaries in addition to the velocity vector (CHoRIN, 1989). As only the velocity vector is usually known (and required by the Navier-Stokes equations) at the boundary, specification of boundary conditions for the pressure presents a considerable complication. The essence of the Dukowicz-Dvinsky approach is as follows. First, the Navier-Stokes equations are discretised both in space and time. During this process the differential operators are transformed into the difference (matrix) operators incorporating the boundary conditions. The resulting difference problem is then augmented by a term which permits the problem to be written in a factored form. The magnitude of the added term should not be greater than the order of the truncation error of the selected time-integration scheme. The modified system is then written in a split form. No intermediate boundary conditions are required for this splitting. Finite-difference/finite-volume discretisations of the Navier-Stokes equations with the velocity and pressure calculated at the same locations produce the matrix operator for the pressure, which has a null space. This property is very undesirable, as it often results in solutions which exhibit chequerboard-type oscillations in the pressure field. This null space can be eliminated by using a staggered grid, where the velocity vector components and the pressure are calculated on a system of interlaced grids (HARLOW and WELCH, 1965). The staggered grid approach, however, works well on Cartesian (PESKIN, 1977) and on orthogonal curvilinear grids only (HAMAKIOTESand BERGER,1988). As is known from differential geometry (DVlNSKY, 1991), curvilinear orthogonal grids can be constructed only for a limited class of three-dimensional domains and as such are not very useful for modelling flow in blood vessels. To extend the staggered grid concept to the general curvilinear non-orthogonal grids required for real vascular Medical & Biological Engineering & Computing

Fig. 2

Comparison of the centre-line velocities with the experiment; the analytical solution for an infinite circular pipe (Womersley) is shown for reference

geometries, the governing equations have to be written in contravariant components. This considerably complicates the equations, and hence takes much longer to compute (DvINSKY, 1987). By modifying the projection operator, we previously (DvINSKYand DUKOWICZ, 1993) formulated a discretisation scheme which uses the Cartesian form of the governing equations on non-staggered, non-orthogonal curvilinear grids, without the penalty of the oscillating pressure fields. This new scheme was used to discretise the governing equations in H E M O in this study. We have used H E M O to calculate problems with up to 125000 grid cells on our Sun Sparcstation 1+ with 40 Mbytes of memory. Such grid resolution is usually sufficient to calculate most laminar and even some transitional flows in a model stenosis geometry. The calculation usually takes between four and five hours per flow cycle. Although the Sparcstation uses a virtualmemory architecture, which in principle allows it to solve problems requiring more than 40 Mbytes of memory, in practice this limit presents a real barrier as the computer performance deteriorates dramatically once the physical memory limit is exceeded. Three-dimensional flow simulations typically create a large number of field values on a set of discrete space/time points. Thus, a 125000 cell flow simulation produces 600000 data points per time step. To process such an amount of information, the visualisation tools are indispensable. In addition to facilitating the data processing, the flow visualisation methods play an important role in elucidating the physics of flow, In this study, we used one of the simplest techniques for flow visualisation, called particle tracking. The essence of this technique is as follows. First, we take a computed set of data for a given time. Secondly, we display this set using either two-dimensional slices or full three-dimensional view of the velocity vector field. Using a mouse which controls the cross hairs, we then can 'inject' the particles at any point in the field. The three-dimensional trajectory of the particle is found by integrating dr

---~-lY

ds

where r is the position vector, s is the parameter, and v = ~(r), i.e. time is fixed. For two-dimensional slices, the

March 1994

139

Fig. 4

Velocity pattern at the Z = 0"7 plane durin9 maximum

flow Fig. 3

Velocity distribution at the Z = 0 plane durin 9 maximum flow; the solid lines show the trajectories of the particles injected at various points

trajectory is calculated using the velocity vectors tangential to the selected viewplane. For completeness, we briefly describe our previous experimental investigation (OJHA et al., 1989; 1990) used in this study. The experimental data were obtained using a photochromic tracer method. The method makes use of the properties of a normally colourless indicator in the test fluid. Upon irradiation with UV light, the indicator undergoes a reversible photochemical reaction resulting in a colour change. The trace produced by a pulsed ultra-violet laser is photographed at a selected time after its formation, thereby enabling the average velocity distribution within the time interval to be obtained. The multiple-trace photochromic method is a relatively new technique, especially for quantitative flow analysis. Therefore, we think it is appropriate to provide a brief summary of its advantages and limitations. A more complete discussion can be found elsewhere (OJHA et al., 1989). The method permits the acquisition of several velocity profiles simultaneously, which considerably facilitates the description of the flow field. The method is very efficient in finding the location of the reattachment point by bracketing its position. Additional advantages of the photochromic tracer method arise from its non-invasive nature and from the ease in generating the trace at any location in the duct and at any instant of the cycle. 3 Results ..

The computer code H E M O , specifically developed for calculating pulsatile flows in vascular geometries, has been used to investigate a sinusoidal flow through a 38% asymmetrically constricted tube. The tube geometry was a model of a mild arterial stenosis, and the flow parameters approximately corresponded to those of a medium-sized human artery with low peripheral resistance. Specifically, the mean and the modulation Reynolds numbers were 575 and 360, respectively, and the frequency parameter (Womersley number) was 7-5. These flow conditions and the model geometry were similar to those used previously (OJHA et al., 1989) in an experimental study using a photochromic tracer method. A computer reconstruction of the vessel geometry is shown in Fig. 1. 140

Fig. 5

View from the top o f the particle trajectories injected during maximum flow

To ensure reliable results, we tried several grid models ranging from 50 000 to 125 000 grid cells. The 125 000 cells is the largest grid size which fits into the memory of our 40 Mbyte Sun Sparcstation 1 + . The calculation time ranged from two (50000 cells) to five (125000 cells) hours per flow cycle. Apparently, 50000 cell resolution was sufficient to obtain grid-independent results especially during the acceleration of the flow. The grid independence was less satisfactory downstream of stenosis during the deceleration phase when the flow became unstable. In comparing the calculated results with our previous data (OsHA et al., 1989; 1990), we found good agreement with the quantitative results reported. These include centre-line velocities, particularly during the acceleration phase (Fig. 2). The wall shear stress behaviour was in qualitative agreement with the data. The same tendency for a better match during the acceleration phase was again noted. Both numerical and laboratory experiments revealed very similar flow patterns. Complex vortical structures were clearly seen downstream of the stenosis during the deceleration phase of the flow cycle using the photochromic tracer visualisation (OmA et al., 1989). However, with the numerical approach, we were able to obtain detailed information on the characteristics of the complex flow field. A secondary flow pattern consisting of a double helix was clearly seen at the inlet and within the throat of the stenosis during peak flow (Fig. 3). This secondary flow resulted from the centrifugal forces generated b y the sharp change in flow dii'ection due to the stenosis. As the helical vortices were convected downstream, their m o m e n t u m was gradually transferred to another pair of counter-rotating helical vortices. As this second pair o f vortices resided mainly within the separation zone (Fig. 4), they travelled in the upstream direction.

Medical & Biological Engineering & Computing

March 1994

0.0050 y,m

0.0050 f y,m l

0.0000

0.0000 I .

.

.

.

.

-0,0050 0.0000 0.0050 0,0100 z,m a

0.0150

-0.0050

0.0000

0,0050 z,m

0.0100

0.0150

a

0.0050

0.0050 0.0025

0.0025

y,m

y,m

0.0000

0.0000 -0.0025

~

~

~

-0,0050

'

~

~

-

0.0100

~'1 -

-,

." 7 .... " , , ,

0,0150

,

.

-0.0025

0.0000

z,m

b Fig. 6

(a) Particle trajectories in the X = 0 plane during peak flow; (b) Expanded view of flow downstream from the constriction 0,0025

y,m

~

o.oooo -

'

S-I

"~

~//~ / ) ~ ; . " . : : 1 ",

, , \ \ \

\

\\\ ,..

t11I

/

/

I'/i/I//,.]

/

. . . .

I

-0.0025 -0.0025

Fig. 7

b Fig. 8

(a) Particle trajectories in the ) = 0 plane during minimum flow," (b) Expanded view of flow around constriction

zone seen previously. Fig. 8 shows the particle trajectories in the midplane of the vessel during minimum flow. Now the flow forms two recirculating zones both upstream and downstream from the constriction, with the downstream zone extending all the way through the outlet. 4 Discussion

~

" -

0.0050 z,m

0.0000 x,m

0,0025

Velocity pattern at the Z = 0 plane during minimum flow

To understand the interaction of these vortices, massless particles were 'injected' into the numerically generated flow field. This simple numerical flow visualisation procedure calculates the particle trajectories under the assumption that the flow field is time-independent, as explained earlier. Fig. 5 shows a three-dimensional view of such trajectories during maximum flow, and Fig. 6 shows the particle trajectories in the midplane, X --- 0, of the vessel. In the acceleration phase, the flow field was symmetric with respect to the symmetry plane of the model (Fig. 4). During the deceleration phase, however, the flow became asymmetric, as seen in Fig. 7, captured during the minimum flow. Such a bifurcation of the flow from a symmetric to an asymmetric pattern is usually an indication of a flow instability that, in turn, appeared to lead to a transitional type flow in the downstream region, as reported previously (OJHA et al., 1989). In addition, an appearance of new discrete frequencies in the wall shear stress was observed at some locations within the separated zone. This phenomenon is usually defined as a subharmonic transition, and it may be related to the non-cyclic fluctuation (with respect to the flow) of the wall shear stress around the reattachment Medical & Biological Engineering & Computing

Considerable difficulties in making accurate measurements of velocities in wall shear stresses in vivo have caused a strong reliance on model studies for obtaining a better understanding of the underlying fluid mechanics. In particular, the problem of flow through constricted arteries has attracted considerable attention from many researchers (OJHA et al., 1989). Unlike laboratory measurements, computer simulations can predict both the core region flow and the boundary layers equally well. Therefore, they have a great potential for becoming a very useful research tool. The main difficulty involved in numerical simulation is the ability of a numerical model to resolve the spatial and temporal significant features of the flow. For example, in the case of H E M O , which uses a second-order accurate spatial discretisation method, the objective of someone constructing a computational grid is to provide a sufficient number of grid cells so that the pressure and the velocity change slowly over the size of the cell. Clearly, this becomes more and more difficult as the Reynolds number increases and the flow structure becomes finer. Fortunately, most interesting haemodynamic problems have Reynolds numbers of the order of 100-1000 which makes them amenable to direct (i.e. without using turbulence models) simulation. The main objective of this work was to evaluate H E M O as a tool for predicting pulsatile flows in constricted tubes on modern desktop workstations. The flows are usually characterised by complex vortical structures and transitional effects and as such present challenging computational problems. As shown in this paper, we can calculate realistic pulsatile flows in constricted tubes using the Sun Sparcstation 1+ several hours per flow cycle. Currently, there are workstations which are significantly faster and which can handle much larger grids. Furthermore, networks of workstations can produce computational power exceeding that of the most powerful supercomputers at a small

M a r c h 1994

141

fraction of the cost, thereby resulting in the solution of larger models of the o r d e r of a million nodes. This creates the o p p o r t u n i t y to model complex h a e m o d y n a m i c problems both in single and bifurcated vessels. The results shown in this paper have also demonstrated that the c o m p u t e r simulations can be very useful as a complementary tool to experimental investigations, particularly in regions of complex flow patterns as seen in the vicinity of the reattachment point d o w n s t r e a m from a stenosis.

Acknowledgments--A. S. Dvinsky was supported in part by US NIH Grant HL 34884.

References AHMED,S., and GIDDENS,D. P. (1984): 'Pulsatile poststenotic flow studies with laser Doppler anemometry,' J Biomech., 17, p. 695 BASSErqGE, E., and BUSSE, R. (1988) Endothelial modulation of vascular tone, Progr. Cardiovasc. Dis., 30, pp. 349-380 CARO, E. G., F~TZGERALO, J. M., and SCHROTER, R. C. (1971): 'Atheroma and arterial wall shear. Observation correlation and proposal for a shear-dependent mass transfer mechanism for atherogenesis,' Proc. R. Soc. London, Bl17, pp. 109-159 CHORtN, A. J. (1968): 'Numerical solution of the Navier-Stokes equations,' Math. Comput., 22, p. 742 CHORIN,A. J. (1989): 'Computational fluid mechanics' (Academic Press) DuKowlcz, J. K., and DVINSKY, A. S. (1992): 'Approximate factorization as a high order splitting for the implicit incompressible flow equations,' J. Comput. Phys., 102, pp. 336-347 DVtNSKu A. S. (1987): 'FLUENT/BFC: a general purpose flow modeling program for all flow speeds' in TAYLOR,C. et al. (Eds.) 'Numerical methods in laminar and turbulent flow' (Pineridge Press) DVINSKV, A. S. (1991): 'Adaptive grid generation from harmonic maps on Riemannian manifolds,' J. Comput. Phys., 95, p. 450 DVlNSKY,A. S., OJng, M., and PERNG, C.-Y. (1991): 'A transitional pulsatile flow through an asymmetrically constricted tube: a numerical investigation" in VANDERBY, R. (Ed.), 'Advances in bioengineering' (ASME) Vol. 20 DvmsKv, A. S., and DuKowlcz, J. K. (1993): 'Null-space-free methods for the incompressible Navier-Stokes equations on non-staggered curvilinear grids,' Computers & Fluids, 22, (6), pp. 685-696 GLAGOV, S., ZARINS, C., G[DDENS, D. P., and Ku, D. (1988): 'Hemodynamics and atherosclerosis," Arch. pathol. Lab. Med., 112, p. 1019 HAMAKIOTES,C. S., and BERGER, S. A. (1988): 'Fully developed pulsatile flow in a curved pipe,' J. Fluid Mech., 195, p. 23

142

HARLOW,F. H., and WELCH,J. E. (1965): Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface,' Phys. Fluids, 8, p. 2182 Ku, D. N., ZEmLER, M. N., and DOWNtNG, J. M. (1990): 'One-dimensional steady inviscid flow through a stenotic collapsible tube," J. Biomech. Eng., 112, pp. 444 450 LANC~LLE,B. L., and O'DoNNELL, F. (1986): 'Reductions in arterial diameter produced by chronic decreases in blood flow are endothelium dependent,' Science, 231, pp. 405-407 LIEBER, B. B., and GIDDENS, D. P. (1990): 'Post-stenotic core flow behavior in pulsatile flow and its effects on wall shear stress,' J. Biomech., 6, p. 597 OJHA, M., COBBOLD,R. S. C., JOHNSTON,K. W., and HUMMEL, R. L. (1989): 'Pulsatile flow through constricted tubes: an experimental investigation using photochromic tracer methods,' J. Fluid Mech., 203, p. 173 OJHA, M., COBBOLD,R. S. C., JOHNSTON,K. W., and HUMMEL, R. L. (1990): 'Detailed visualization of pulsatile flow fields produced by modelled arterial stenoses,' J. Biomed. Eng., 12, pp. 463-469 PaTEL, D. J., and VAISHNAV,R. N. (1980): 'Basic hemodynamics and its role in disease process' (University Park Press, Baltimore, Maryland) PESKm, C. S. (1977): 'Numerical analysis of blood flow in the heart,' J. Comp. Phys., 25, p. 220 RITrGERS, S. E., and SHU, M. C. S. (1990): 'Doppler color-flow images from a stenosed arterial model. Interpretation of flow patterns,' J. Vasc. Surg., 12, p. 51 1 YOUNG, D. F., and TSAr, F. Y. (1973): 'Flow characteristics in models of arterial stenoses,' J. Biomech., 6, pp. 395; 547

Author's b i o g r a p h y A 1977 graduate of St. Petersburg Mining Institute, Dr_ Dvinsky received his PhD in Mechanical Engineering from the University of Houston in 1983. His doctoral research involved computer modelling of microcirculation and transport of oxygen from capillaries to the tissue. From 1984 to 1992, Dr. Dvinsky worked at Creare Inc., Hanover, New Hampshire, where he was awarded and directed a number of grants and contracts which involved research in the areas of mathematical modelling of blood flow, heat and mass transfer, and electromagnetic power deposition in the tissue. Since t990, Dr. Dvinsky has been a member of the research staff at Harvard Medical School. In 1992, he became a president of Daat Research Corp., where he has directed research in hyperthermia devices and planning and stenosis severity assessment systems.

Medical & Biological Engineering & Computing

March 1994