- Email: [email protected]
Computer prediction of thrombogenic sites for a tilting-disc prosthetic heart valve
COMPUTERS
AND
BIOMEDICAL
Computer
Division
RESEARCH,
l&165-182
(1977)
Prediction of Thrombogenic Tilting-Disc Prosthetic Heart
Sites Valve*
ANTHONY
D. Au AND HARVEY S. GREENFIELD
of ArtiJicial
Organs, Department of Surgery, Salt Lake City, Utah 84112
University
for a
of Utah,
Received March 22,1976 In an attempt to identify the hemodynamic factors that contribute to thrombus growth on tilting-disc prosthetic heart valves, a mathematical model was developed to simulate blood flow through the Bjork-Shiley-type aortic prosthesis. The Navier-Stokes equations which describe the flow were solved numerically on a digital computer, and the solutions were displayed on a cathode-ray tube in pictorial form. Results included a detailed pressure variation across the flow regime surrounding the tilted disc. Such a plot provides guidelines for the placement of pressure probes for laboratory simulations. Velocity vector displays indicated certain flow stagnation regions; shear stress and normal stress distributions were made available; simulated steady flow at a relatively low Reynolds number showed interesting comparisons to thrombogenic sites in necropsy findings. Based upon such simulations, it is hypothesized that possible red blood cell damage due to high fluid stress and subsequent deposition of these degenerated cells at certain sites could be the cause of the noted thrombi. The totality of the displays can permit one to obtain proper blood flow patterns necessary for optimal prosthesis design.
Although prosthetic valves have improved during the past few years, it appears that the ultimate designhas not been attained. New models are still being presented to the clinician. Thromboembolism is present and postoperative complications are not rare. An artificial heart valve design study, however, is not an easy one. In such a design problem the resistanceto flow should be minimized, stasisand separation regions must be avoided, as well as unnecessary turbulence, and regurgitation is undesirable. Turbulence may contribute to energy loss, vascular wall damage, and hemolysis. The blood flow field should be disturbed as litt!e as possible, and the flow occluder should respond well to small pressure differences. The design of the flow occluder is important, since area ratio of occluder to orifice and the occluder shape with its accompanying restrictive cage are necessary parameters of study. The flow characteristics of such occluders and their relative performances are usually evaluated in mock circulation equipment; however, it is known that the clinical * This research was supported by the National Institutes of Health, Grant HL-12202 and by Program Project Grant H3-POl-HL-13738-0281. 165 Copyright 0 1977 by Academic Press, Inc. All rights of reproduction Printed in Great Britain
in any form
reserved.
ISSN
00104809
166
AU
AND
GREENFIELD
significance may not be seen in such a pseudo environment. It seems desirable therefore to complement experimental studies with properly formulated theoretical analyses. The application of computer graphics to analyses of medical phenomena permits the solution of involved mathematical equations within the analyses plus a true man-machine interaction within the investigation. Our past computer-formed
FIG. 1. The Bjork-Shiley
tilting disc prosthetic heart valve (without cloth sewing rim).
investigations have resulted in simulations of blood movement through ball and common disc-type cardiac valve prostheses(l-3), and observations have beenmade as related to atherogenesis(4-6). The present discussioncenters about the type of prosthetic heart valve which contains a disc occluder that pivots eccentrically. Such a valve analysis was briefly treated in conjunction with ventricular-aortic flow simulation for the caseof an artifician heart (7). The analysis will be treated here in detail for the tilting-disc situation alone. The tilting-disc prosthesis, with its central flow through the major orifice and low-pressure gradient features, has drawn considerable clinical emphasis. Of the three basic designs-the Wada-Cutter, the Lillehei-Kaster, and the Bjork-Shiley-
PROSTHETIC
HEART
VALVE
SIMULATIONS
167
the latter, and most widely used, was chosen for a computer simulation model. The first design was rejected because of functional irregularities with resulting valvular insufficiency (a), and the second because of its relatively short clinical history. Also, the Bjork-Shiley prosthesis, as shown in Fig. 1, is presently used within Dr. Willem Kolff’s laboratory (authors’ affiliation) as part of the totally implantable, artificial heart studies.
FIG.
2. A necropsy specimen showing heavy thrombus formation on a Bjork-Shiley prosthesis
(IO). Although the dynamics of the Bjork-Shiley valve appear favorable, thromboembolic complications, such as seen in Fig. 2, have been reported (9, IO). Experience within our laboratory with insertion of artificial hearts into calves, has shown similar results (see Fig. 3). The computer-generated simulations presented here depict flow patterns about the heart valve of interest, the stasis regions, the pressure drop, and the stress distributions; suggesting correlation between such parameters and thrombotic conditions.
168
AU
AND
GREENFIELD
FIG. 3. Thrombus formation on a Bjork-Shiley
after placement in a calf (authors’ laboratory) (7).
THEORY
The computer simulation procedure involves solving fluid dynamic equations together with the prescribed boundary conditions which specify the geometry and flow conditions of the physical system. Since an analytical solution to the fluid dynamic equations is possible only for very limited situations, finite difference approximations are generally applied to discretize the solution space into mesh
PROSTHETIC
HEART
VALVE
SIMULATIONS
169
points and to replace the governing partial differential equations with a large but simplified set of algebraic equations. With the solution known only at discrete points, specification of the boundary conditions becomes difficult when the physical boundaries do not coincide with the mesh points. At the full open position the tilting-disc valve presented such difficulty. To reduce the complexity of the problem while maintaining a reasonable approximation to the physical system, a twodimensional model based upon the plane through the axis of symmetry was first adopted. Although the ventricular-aortic section in which the tilting-disc valve resides is cylindrical in nature, a rectangular coordinate system had to be used, since the boundary specification of the tilting-disc valve in a two-dimensional cylindrical coordinate system is meaningless. This simplified model can be mathematically treated as an infinitely long tilted plate situated in an infinitely long channel with varying cross section. During the major portion of the systolic cycle, the aortic valve stays at the fully open position to allow the blood to be discharged from the ventricle into the aorta. For the purpose of correlating blood flow characteristics with thrombus formation, it was decided to characterize the major part of the systolic cycle as steady flow. In this simplified simulation model, it was further assumed that the non-Newtonian contribution of the blood was insignificant and that the blood exhibited Newtonian behavior. Based on the foregoing assumptions, the fundamental equation for two-dimensional incompressible flow of a Newtonian fluid with no body forces and constant properties reduced to the two momentum equations (Navier-Stokes) and the continuity equation
av av "=-3!+L-(~+~), Z+"aX+viY~
ri i
1’1
PI The equations were written in the dimensionless primitive variables of velocity components, u and u, and pressure, p. The variables were rendered dimensionless with respect to the average velocity U and the diameter D of the aorta. The Reynolds number, Re (a similarity parameter that includes the blood viscosity), is defined as pUD/p, where p is the density and p is the viscosity of the blood. The asymmetric property of the tilting-disc geometry made it difficult to determine a stream function boundary value on the disc in the vorticity-stream function approach, and for this reason the marker and cell (MAC) method (II) of solution was applied. In this method the solution space is subdivided into rectangular cells, and the primitive variables are assigned to the cell in the manner shown in the inset
170
AU AND GREENFIELD
of Fig. 4. The geometry of the simulation model represented the ventricular-aortic section of the heart. As seenin Fig. 4, a scaleddrawing of the flow area wasoverlayed with cells for computational purposes. To easecertain difficulties inherent in the numerical technique employed, the boundaries of the blood flow channel were made to lie along the edgesof the computational cells. This necessitatedthe shaping of the curve of the ventricle wall (dash line) to be performed by a mathematical transformation of a straight line in the computational space. Previous analysis without the curved ventricle wall and sewing ring produced different flow patterns (7). To
~SEWING
RING
------
-CD-
VENTRICLE
FIG.
4.
The computational
1 AORTA
model of a tilting-disc heart valve in the ventricular-aortic
region.
accommodate simulation of the fully open disc in the aortic position at a tilted angle of 60”, a rectangular meshin the aortic section was realized, equal to half of the grid spacing upstream. Such an action wasjudiciously chosen to provide sufficient data points for describing the flow variation, and to permit operation within available computer memory with desired computation efficiency. Along boundary lines no-slip boundary treatment was employed, requiring zero fluid velocity along solid boundaries. The finite difference representation of the u-component momentum equation could be written as U?+’ I++.j = YY+;-:. j-
(st/sx>
(Pi+l.
j -Pi,
j),
Pal
PROSTHETIC
HEART
VALVE
171
SIMULATIONS
where (U’)l+1,j $+;.
j
=
Uy++,
j
+
6t
-
(U’):,
j
_
(uv)l++,
J++
-
CUUX+‘,t,
j-+
-
6X
SY
1 Ul++,j - 2Uf++,j + G-t, j + UY+‘+t. J+I - 2uY++,j + G+;O.j-1
+iie
-.
(W2
i
@YY
)I
WI
Similarly, the v-component momentum equation became II+1 ui. j++
=
tl.
J+)
-
(6f/6Y)
(Pi,
.I+1
Pal
-Pi.jh
where
(U”);, j+l
- Co’):,j _ CuvX++.i++ - (uv)Y-+,J++ 6X SY
-
2Vl,
j++
@YY
+
V;, j-+
+
Vy+l,
j+*
-
2Vls
j++
+
VT-l,
j++
Vbl
@x)2
The subscripts i andj indicate the x and y coordinates, respectively, and the superscript n, the time index. Linear interpolation was usedto determine the flow variables not stored in the computer but required in Eqs. [3b] and [4b]; thus
The pressurefield pi, j was to be determined such that the velocity components tPl and v”+l obtained from Eqs. [3a] and [4a] satisfied the continuity equation (12) [71 The pressureand advance-time velocity were then solved simultaneously in an iterative procedure. Upon introducing a relaxation parameter 6r, the I + 1 iterate of the pressurefield was obtained from the velocity field, pi;; =p;, i - sr(Dy.‘i’)‘.
PI
The new pressurefield was substituted into Eqs. [3a] and [4a] to obtain the I + 1 iterate of the velocity components. The new velocity components were then used in Eq. [7] to compute the divergence. This iterative processcontinued until the changes in the iterates pf, j, (uf,‘:, j)’ and (v;,‘&)’ fell below the convergence criteria. To reduce computation time, the advection and viscous terms stored in q and 5 were updated only at the beginning of a time cycle. When a steady-state solution was
172
AU
AND
GREENFIELD
reached, the shearing stress and normal stress values were determined velocity components according to Eqs. [9] and [lo], respectively :
r,, = (2/Re) (&+3x).
from the
UOI
By definition, the stream function values were then computed with Eq. [ 1I] and the vorticity values with Eq. [12] :
0 = (avjax) - (au/i+).
[121
A fully developed velocity profile for channel blood flow was assumed at the inflow boundary. The pressure value at the inlet to the aortic region was determined upon the basis of the inflow velocity value according to Eqs. [7] and [8]. In a confined flow, there are an infinite number of solutions for the pressure field, differing only by an additive constant. It was therefore necessary to establish a reference value in the pressure field; thus, zero pressure value was specified at the outflow boundary. The outflow velocity components were assigned the corresponding velocity values immediately upstream. This approximation is valid only if the local variation in the velocity field is small. Again, for a particular blood flow velocity, this condition is satisfied only if the outflow boundary is far enough downstream from the disc that causes velocity profile deviation from a Poiseuille-flow profile. For a given number of mesh points downstream from the disc, this outflow boundary treatment, in essence, limits the range of Reynolds numbers for which the solution is applicable. METHOD
The solution of the finite difference equations requires a computer with sufficient memory capacity and an efficient algorithm. Since the accuracy of the approximation is dependent on the number of data points chosen to represent the solution space, the numerical scheme generally produces a vast amount of data. In this study the approximation to the flow in the ventricular-aortic root region is represented by approximately 2000 data values for each of the three variables. The analysis of data becomes a very difficult and time consuming task. With the employment of computer graphics, however, the data can be presented in pictorial form almost instantly. Not only is the data analysis simplified, but the graphics representation provides a close resemblance to the physical system under investigation. The computer equipment employed consists of a PDP-10 time-sharing computer, magnetic disk units, an IMLAC graphics display unit with console keyboard, a Tektronix storage display unit, and a teletype. Essential parameters, such as the Reynolds number, the time steps, and the convergence criteria, are fed into the
PROSTHETIC HEARTVALVE
SIMULATIONS
173
computer via the teletype at the initiation of a simulation run. As the computer proceeds to solve the finite difference equations iteratively, the change in the basic variables is monitored for convergence. The changes for all three variables between time steps are printed on the teletype to provide a permanent record. Since a typical simulation may require hours of computation within the time-sharing environment, it is necessary to store the intermediate solution onto disk for continued computation at a later time. This also serves to protect against total loss of data due to system failure. The intermediate solution is copied onto disk upon specific interrupt signals generated from the teletype or at the completion of the prescribed number of iterations. The stored data can be processed by the graphics routines and the results displayed on the screen for analysis. The velocity vector plot is generally used when examining the intermediate results, because such a display can be generated with the least amount of processing time. The vector plot has been found to be very helpful in pointing out program errors during the debugging stage and in signaling numerical instability during simulation. In this study the solution was very sensitive to the change in time step and the relaxation parameter. The solution is considered to have reached a steady state when the maximum change in velocity value between time steps is less than 1 ‘A. A different computer program is initiated to generate stream function, vorticity, and stress solutions from the three primary variables. These available data are then displayed on the screen in pictorial form. The complexity of the display necessitates the employment of a storage display unit. A camera mounted on the display unit records the images on film for future reproduction purposes. RESULTS
The conversion of the numerical values into picture forms was automated within the computer program and results were placed into contour map sequences. The contour map is a common form of data representation in field studies. A family of constant value lines can, in general, provide adequate information about the solution space. However, depending on the contour level selected, significant information can be concealed within a contour map; for example, additional contour levels have to be specified in the stream function display to bring out the eddy formation distal to the sewing ring. Each of the other contour maps shown herein is drawn with 10 equally spaced contour levels which encompass the range of the function values. The positive and negative values are indistinguishable when contours cannot be conveniently labeled, as in the case of computer-generated displays. To supplement the contour map, a program was designed to generate a surface representation of the solution space. The surface shown in the isometric plot represents the function values over the entire region of interest. For clarity, the hidden portion of the surface is eliminated. This very often results in the omission of negative function values. With the surface plot contour values displayed on the
174
AU
AND
GREENFIELD
screen can be assigned to the contour lines. The combination of the contour map and the surface plot can provide sufficient information about the solution to allow detailed analysis. The pressure drop across a heart valve is an important measure of the heart valve performance. In general, lower pressure drop is preferred, as high pressure drop
REYNOLDS
NO.
= 100.0
MAX = CONTOUR
-7.63975
0.68765 9.01544 FIG.
-5.97423
2.35336 10.68096
12.3465 LEUELS
-4.30671 4.01906
-2.64319 5.69440
IIIN
=
-9.3093 -0.97767 7.34992
5. Computer-formed pressure field about the tilting disc within the aortic root region.
across the heart valve means more energy is required by the heart to pump the blood The high pressure gradient across the valve also means higher blood flow velocity, which in turn causes high shear rate. Figure 5 represents the pressure distribution for the simulation model under study. The surface plot provides the detailed pressure variation across the tilting-disc valve which is not generally obtainable from mock circulation experiments. In most experimental studies, the pressure drop is evaluated with the pressure values measured at two arbitrarily chosen points
into circulation.
PROSTHETIC
HEART
VALVE
175
SIMULATIONS
across the valve. It is obvious, from Fig. 5, that correct pressure drops can only be calculated if the measuring devices are placed at appropriate locations. The resulting display, obtained from the numerical simulation can provide guidelines for placement of the pressure measuring probes in the laboratory simulation of a more accurate heart valve model. The pressure data given in Fig. 5 provide a valuable
-------_-__--______-__________________ REYNOLDS
NO.
= 100.0
-0.14648 2.32215
0.34725 2.81587
FIG. 6. Computer-formed the aortic root region.
horizontal
MAX
=
4.7908
CONTOUR LEUELS 0.84097 3.30960
tllN
1.33470 3.60333
=
-0.6402
1.82842 4.29705
velocity component of flow about the tilting disc within
measure of the effectiveness of the tilting-disc heart valve when comparison is made with other heart valve models studied under similar conditions. The two velocity components, u and v, are displayed in Figs. 6 and 7, respectively. High horizontal velocity is indicated near the two ends of the disc as a result of the narrowing of the flow path between the tilted disc and the channel wall. As the fluid passes through the minor orifice it is directed by the sewing ring to move vertically
176
AU AND GREENFIELD
toward the underside of the disc. The change in flow direction is indicated in Fig. 7 by the peak vertical velocity component near the sewing ring. Although sufficient information is contained in the two velocity component displays to determine the flow about this particular tilting-disc configuration, the vector
REYNOLDS
-2.459111 0.23479 FIG. 7. Computer-formed aortic root region.
NO.
= ~00.0
-1.92032 0.77357
IIAX
=
2.9297
CONTOUR LEUELS -1.39154 1.31235
-0.64276 1.95112
IIIN
=
-2.9979
-0.30399 2.39990
vertical velocity component of flow about the tilting disc within the
sum of the velocity components provides a better visual perception of the flow variation. In a velocity vector display the length and direction of the vector at a point represent the magnitude and direction of the velocity. The short vectors shown in Fig. 8 indicate a stagnation region adjacent to the channel wall downstream of the disc. Due to the small-scale factor used to draw the vectors, the eddy formation behind the sewing ring is not as clearly represented as in the streamline plot.
AU AND GREENFIELD
178
A fourth-order Ralston-Runge-Kutta numerical integration technique was employed to evaluate the stream function values with Eq. [l I]. With the velocity defined asthe gradient of the stream function, the closely spacedstreamlinesshown in Fig. 9 represent high-velocity regions. A closed streamline contour representsan eddy formation. The contribution of eddiesin blood tlow to the thrombus formation will be discussed.
---------__ REYNOLDS
-95.49675
47.10165 FIG.
10. Computer-formed
NO.
.._ ._____---___-------= 100.0
-66.97707
75.62133
HAX
= 189.7001
CONTOUR LEVELS -36.45739
104.14101
,,IN
-9.93771
132.66069
--=-124.0164 1a.59197
161.10037
vorticity field about the tilting disc within the aortic root region.
The vorticity function may be considered asa measureof fluid rotation. The signs associatedwith the vorticity values indicate the direction of rotation. In Fig. 10 the peak vorticity values are shown near the solid boundaries. In general, the vorticity contour maps give insight into the interaction of the flow with the boundaries. The shearstressdistribution and the normal stressdistribution within the fluid are presentedin Figs. 11and 12, respectively. The signson the shearstressvalues indicate
PROSTHETIC
REYNOLDS
NO.
= 100.0
HEART
NAX CONTOW?
-1.29109 0.04525 FIG.
11. Computer-formed
-1.02392 0.91252
VALVE
1.3916
=
179
SIMULATIONS
nlN
=
-1.5504
LEVELS
-0.75955 0.57979
-0.40529 0.94709
-0.2aao2 1.11433
shear stress field about the tilting disc within the aortic root region.
the direction of shear, while the signson the normal stressvalues distinguish between the tensive and compressive stresses.These two figures clearly show that the fluid passingthrough the small aperture of the valve experiencesthe most severestress. DISCUSSION
If hemodynamic theories of thrombogenesisare establishedon valid grounds, the results of the blood flow simulation through the tilting-disc heart valve model, although simplified, should bear evidence to support the theories. In order to attempt such verification, it seemsreasonableto inspect the results shown in Figs. 5 through 12. The pressuredistribution is smooth everywhere except at the minor aperture of the valve. The abrupt changein pressureresults in higher velocity which, in turn, causes higher stress.The results shown in Figs. 6, 7, and 8 are summarized in Fig. 9. The
180
AU
AND
GREENFIELD
point of interest in this stream function display is the eddy formation distal to the sewing ring and the less intense eddy further downstream. The vorticity function shows regions of severe fluid-wall interaction at both ends of the disc. Since both the vorticity and shear stress are functions of velocity gradients, some similarity is expected. The shear stress field indicates high shear stress about the upper vascular
REYNOLDS
-1.56791 -0.26877
FIG. 12. Computer-formed
NO.
= 100.0
-1.30808 -0.00894
flax
=
1.0304
CONTOUR LEUELS -1.04826 0.25089
-0.78843 0.51072
MIN
=
-1.8277
-0.52160 0.77055
normal stress field about the tilting disc within the aortic root region.
wall and the two ends of the disc, with stressconcentration at the small aperture of the valve. The normal stressdistribution also shows a high compressive stressfield below the disc. From the foregoing observations, it can be expected that any results of hemodynamically induced thrombogenesiswill occur on the underside of the tilting disc. A study of Figs. 2 and 3 confirms such expectations. A theory for the mechanics of thrombus formation is not possible, of course, without a full understanding of the chemical and biological effects which undoubtedly play important roles. However,
PROSTHETIC
REARTVALVE
SIMULATIONS
181
based on the results of this study perhaps a hypothesis for thrombus formation, although tenuous, might be stated. As red blood cells are destroyed or damaged, becauseof imposed high stress,and are subsequently entrapped in the stasisregion near the vascular wall, the contents of the red cells trigger the extrinsic clotting reaction. The continuous deposition of degeneratedred blood cellsat that site could lead to a thrombus. If the struts were included in the simulation, thrombus growth from the sewing ring position outward might duplicate that seenin Fig. 2. Information on red blood cell destruction due to fluid stresshas been reported in several publications (13,14), and a detailed discussionconcerning the mechanicsof thrombus growth on surfaceshas been given by Blackshear (15). CONCLUSION
To reduce thrombotic stenosis,a heart valve should be designedto minimize stress and eddy formation in the flow field. The natural heart valve with its unobstructed central flow feature possesses these desirable properties. One designobjective of the tilting-disc valve was to provide a similar central flow feature by reducing the flow obstruction. While the obstruction to the flow is undoubtedly less,the disturbance introduced into the blood stream by the tilting disc may be more than that of the ball or the disc prostheses.The stream function display shown in Fig. 9 indicates that the disc affects the flow up to two aortic diameters downstream from the aortic root, while previous analyses on the ball and disc valves at similar Reynolds numbers have shown only minor flow disturbance distal to the valve (26). Since the asymmetrical property of the tilting-disc valve is the causeof the added flow disturbance, the results of this study suggesta reemphasisof symmetrical flow patterns in future heart valve design. It is believed that blood flow experimentation by computer can be useful in designing an optimal heart valve. To reduce or eliminate eddies, a heart valve may be designedsuch that the flow is directed into areaswhich would otherwise contain or produce large eddies.The intensity of stressconcentrations in the flow field may be reduced by streamlining the heart valve components. The use of computer simulation in such a design processwould complement laboratory experimentation and might prove to be both lesscostly and lesstime consuming. In order to vary the shape of the heart valve components, however, the irregular boundary problems discussed in this report must be overcome. The finite element approach could be the solution to this problem, and is presently a topic of research. ACKNOWLEDGMENTS
Appreciationisgivento the editorsof theJournal of Curdiouascular Surgery (Turino) for permissionto publishFig. 2, andto theeditorsof the Transactions of the American Society of Artificial Internal Organs for permission to publishFig. 3.Thanksalsoio Dr. WillemKolff for his continuing supportof thisproject. 7*
182
AU AND GREENFIELD REFERENCES
I. GREENFIELD, H. Special address: Studies of medical phenomena by computer graphics. Trans. Amer. Sot. Artif. Int. Organs l&607-61 2 (1972). 2. GREENFIELD, H., AND KOLFF, W. The prosthetic heart valve and computer graphics. J. Am. Med. Assoc. 219,69-74 (1972). 3. Au, A., AND GREENFIELD, H. Computer graphics analysis of stresses in blood flow through a prosthetic heart valve. Comput. Biol. Med. 4,279-291 (1975). 4. SANDBERG, L., REEMTSMA, K., AND GREENFIELD, H. Some theoretical aspects of vascular degeneration. Amer. J. Surg. 119, 543-552 (1970). 5. GREENFIELD, H., VICKERS, D., SUTHERLAND, I., KOLFF, W., AND REEMTSMA, K. Moving computer graphics images seen from inside the vascular system. Trans. Amer. Sot. Artif. Int. Organs 17,381-385 (1971). 6. CANNON, T., AND GREENFIELD, H. Adaptability of computer graphics to studies of atherosclerosis pathogenesis. J. Assoc. Sot. Adv. Med. Instrum. 6, 250-255 (1972). 7. GREENFIELD, H., Au, A., KELSEY, S., AND KOLFF, W. Simulation of assumed detriments to prosthetic heart function. Trans. Amer. Sot. Artif. Int. Organs 20,673-679 (1974). 8. HALLMAN, G. L., MESSMER, B. J., ELKADI, A., VON DER EMDE, K., AND COOLEY, D. A. Clinical experience with the Wada-Cutter cardiac valve prosthesis. Ann. Thordc. Surg. 10,9-19 (1970). 9. COKKINOS, D. V., VORIDIS, V., BAKOULAS, B., THEODO~~IOU, A., AND SKALKEAS, G. D. Thrombosis of two high-flow prosthetic valves. J. Thoruc. Curdiovusc. Surg. 62,947-949 (1971). 10. BOZER, A. Y., AND KARAMEHETOGLU, A. Thrombosis encountered with Bjork-Shiley prosthesis. J. Cardiovasc. Surg. (Turino) 13,141-143 (1972). Il. WELCH, J. E., HARLOW, F. H., SHANNON, J. P., AND DALY, B. J. “The MAC Method.” Los Alamos Scientific Laboratory Report LA-3425,I966. 12. VIECELLI, J. A. “Computing Method for Incompressible Flows Bounded by Moving Walls.” Lawrence Radiation Laboratory Report UCRL-728 15, 1970. 13. FORSTROM, R. J., BLACKSHEAR, P. L., JR., KESHAVIAH, P.. AND DORMAN, F. D. Fluid dynamic lysis of red cells. Chem. Eng. Prog. Symp. Ser. 67,69-74 (197; ). 24. NEVARIL, C. G., LUNCH, E. C., ALFREY, C. P., AND HELLUMS, J. D. Erythrocyte damage and destruction induced by shearing stress. J. Lab. Clin. Med. 71,784-790 (1968). 15. BLACKSHEAR, P. L., JR., AND WATTERS, C. Observation of red blood cells hitting solid walls. Chem. Eng. Prog. Symp. Ser. 67,60-68 (1971). 16. GREENFIELD, H., AND DEBRY, R. “An Application of Computer Graphics: Two Concurrent Investigations within the Medical Field.” Technical Report UTEC-CSc-71-I 15. Computer Science Department, University of Utah, 1971.
AND
BIOMEDICAL
Computer
Division
RESEARCH,
l&165-182
(1977)
Prediction of Thrombogenic Tilting-Disc Prosthetic Heart
Sites Valve*
ANTHONY
D. Au AND HARVEY S. GREENFIELD
of ArtiJicial
Organs, Department of Surgery, Salt Lake City, Utah 84112
University
for a
of Utah,
Received March 22,1976 In an attempt to identify the hemodynamic factors that contribute to thrombus growth on tilting-disc prosthetic heart valves, a mathematical model was developed to simulate blood flow through the Bjork-Shiley-type aortic prosthesis. The Navier-Stokes equations which describe the flow were solved numerically on a digital computer, and the solutions were displayed on a cathode-ray tube in pictorial form. Results included a detailed pressure variation across the flow regime surrounding the tilted disc. Such a plot provides guidelines for the placement of pressure probes for laboratory simulations. Velocity vector displays indicated certain flow stagnation regions; shear stress and normal stress distributions were made available; simulated steady flow at a relatively low Reynolds number showed interesting comparisons to thrombogenic sites in necropsy findings. Based upon such simulations, it is hypothesized that possible red blood cell damage due to high fluid stress and subsequent deposition of these degenerated cells at certain sites could be the cause of the noted thrombi. The totality of the displays can permit one to obtain proper blood flow patterns necessary for optimal prosthesis design.
Although prosthetic valves have improved during the past few years, it appears that the ultimate designhas not been attained. New models are still being presented to the clinician. Thromboembolism is present and postoperative complications are not rare. An artificial heart valve design study, however, is not an easy one. In such a design problem the resistanceto flow should be minimized, stasisand separation regions must be avoided, as well as unnecessary turbulence, and regurgitation is undesirable. Turbulence may contribute to energy loss, vascular wall damage, and hemolysis. The blood flow field should be disturbed as litt!e as possible, and the flow occluder should respond well to small pressure differences. The design of the flow occluder is important, since area ratio of occluder to orifice and the occluder shape with its accompanying restrictive cage are necessary parameters of study. The flow characteristics of such occluders and their relative performances are usually evaluated in mock circulation equipment; however, it is known that the clinical * This research was supported by the National Institutes of Health, Grant HL-12202 and by Program Project Grant H3-POl-HL-13738-0281. 165 Copyright 0 1977 by Academic Press, Inc. All rights of reproduction Printed in Great Britain
in any form
reserved.
ISSN
00104809
166
AU
AND
GREENFIELD
significance may not be seen in such a pseudo environment. It seems desirable therefore to complement experimental studies with properly formulated theoretical analyses. The application of computer graphics to analyses of medical phenomena permits the solution of involved mathematical equations within the analyses plus a true man-machine interaction within the investigation. Our past computer-formed
FIG. 1. The Bjork-Shiley
tilting disc prosthetic heart valve (without cloth sewing rim).
investigations have resulted in simulations of blood movement through ball and common disc-type cardiac valve prostheses(l-3), and observations have beenmade as related to atherogenesis(4-6). The present discussioncenters about the type of prosthetic heart valve which contains a disc occluder that pivots eccentrically. Such a valve analysis was briefly treated in conjunction with ventricular-aortic flow simulation for the caseof an artifician heart (7). The analysis will be treated here in detail for the tilting-disc situation alone. The tilting-disc prosthesis, with its central flow through the major orifice and low-pressure gradient features, has drawn considerable clinical emphasis. Of the three basic designs-the Wada-Cutter, the Lillehei-Kaster, and the Bjork-Shiley-
PROSTHETIC
HEART
VALVE
SIMULATIONS
167
the latter, and most widely used, was chosen for a computer simulation model. The first design was rejected because of functional irregularities with resulting valvular insufficiency (a), and the second because of its relatively short clinical history. Also, the Bjork-Shiley prosthesis, as shown in Fig. 1, is presently used within Dr. Willem Kolff’s laboratory (authors’ affiliation) as part of the totally implantable, artificial heart studies.
FIG.
2. A necropsy specimen showing heavy thrombus formation on a Bjork-Shiley prosthesis
(IO). Although the dynamics of the Bjork-Shiley valve appear favorable, thromboembolic complications, such as seen in Fig. 2, have been reported (9, IO). Experience within our laboratory with insertion of artificial hearts into calves, has shown similar results (see Fig. 3). The computer-generated simulations presented here depict flow patterns about the heart valve of interest, the stasis regions, the pressure drop, and the stress distributions; suggesting correlation between such parameters and thrombotic conditions.
168
AU
AND
GREENFIELD
FIG. 3. Thrombus formation on a Bjork-Shiley
after placement in a calf (authors’ laboratory) (7).
THEORY
The computer simulation procedure involves solving fluid dynamic equations together with the prescribed boundary conditions which specify the geometry and flow conditions of the physical system. Since an analytical solution to the fluid dynamic equations is possible only for very limited situations, finite difference approximations are generally applied to discretize the solution space into mesh
PROSTHETIC
HEART
VALVE
SIMULATIONS
169
points and to replace the governing partial differential equations with a large but simplified set of algebraic equations. With the solution known only at discrete points, specification of the boundary conditions becomes difficult when the physical boundaries do not coincide with the mesh points. At the full open position the tilting-disc valve presented such difficulty. To reduce the complexity of the problem while maintaining a reasonable approximation to the physical system, a twodimensional model based upon the plane through the axis of symmetry was first adopted. Although the ventricular-aortic section in which the tilting-disc valve resides is cylindrical in nature, a rectangular coordinate system had to be used, since the boundary specification of the tilting-disc valve in a two-dimensional cylindrical coordinate system is meaningless. This simplified model can be mathematically treated as an infinitely long tilted plate situated in an infinitely long channel with varying cross section. During the major portion of the systolic cycle, the aortic valve stays at the fully open position to allow the blood to be discharged from the ventricle into the aorta. For the purpose of correlating blood flow characteristics with thrombus formation, it was decided to characterize the major part of the systolic cycle as steady flow. In this simplified simulation model, it was further assumed that the non-Newtonian contribution of the blood was insignificant and that the blood exhibited Newtonian behavior. Based on the foregoing assumptions, the fundamental equation for two-dimensional incompressible flow of a Newtonian fluid with no body forces and constant properties reduced to the two momentum equations (Navier-Stokes) and the continuity equation
av av "=-3!+L-(~+~), Z+"aX+viY~
ri i
1’1
PI The equations were written in the dimensionless primitive variables of velocity components, u and u, and pressure, p. The variables were rendered dimensionless with respect to the average velocity U and the diameter D of the aorta. The Reynolds number, Re (a similarity parameter that includes the blood viscosity), is defined as pUD/p, where p is the density and p is the viscosity of the blood. The asymmetric property of the tilting-disc geometry made it difficult to determine a stream function boundary value on the disc in the vorticity-stream function approach, and for this reason the marker and cell (MAC) method (II) of solution was applied. In this method the solution space is subdivided into rectangular cells, and the primitive variables are assigned to the cell in the manner shown in the inset
170
AU AND GREENFIELD
of Fig. 4. The geometry of the simulation model represented the ventricular-aortic section of the heart. As seenin Fig. 4, a scaleddrawing of the flow area wasoverlayed with cells for computational purposes. To easecertain difficulties inherent in the numerical technique employed, the boundaries of the blood flow channel were made to lie along the edgesof the computational cells. This necessitatedthe shaping of the curve of the ventricle wall (dash line) to be performed by a mathematical transformation of a straight line in the computational space. Previous analysis without the curved ventricle wall and sewing ring produced different flow patterns (7). To
~SEWING
RING
------
-CD-
VENTRICLE
FIG.
4.
The computational
1 AORTA
model of a tilting-disc heart valve in the ventricular-aortic
region.
accommodate simulation of the fully open disc in the aortic position at a tilted angle of 60”, a rectangular meshin the aortic section was realized, equal to half of the grid spacing upstream. Such an action wasjudiciously chosen to provide sufficient data points for describing the flow variation, and to permit operation within available computer memory with desired computation efficiency. Along boundary lines no-slip boundary treatment was employed, requiring zero fluid velocity along solid boundaries. The finite difference representation of the u-component momentum equation could be written as U?+’ I++.j = YY+;-:. j-
(st/sx>
(Pi+l.
j -Pi,
j),
Pal
PROSTHETIC
HEART
VALVE
171
SIMULATIONS
where (U’)l+1,j $+;.
j
=
Uy++,
j
+
6t
-
(U’):,
j
_
(uv)l++,
J++
-
CUUX+‘,t,
j-+
-
6X
SY
1 Ul++,j - 2Uf++,j + G-t, j + UY+‘+t. J+I - 2uY++,j + G+;O.j-1
+iie
-.
(W2
i
@YY
)I
WI
Similarly, the v-component momentum equation became II+1 ui. j++
=
tl.
J+)
-
(6f/6Y)
(Pi,
.I+1
Pal
-Pi.jh
where
(U”);, j+l
- Co’):,j _ CuvX++.i++ - (uv)Y-+,J++ 6X SY
-
2Vl,
j++
@YY
+
V;, j-+
+
Vy+l,
j+*
-
2Vls
j++
+
VT-l,
j++
Vbl
@x)2
The subscripts i andj indicate the x and y coordinates, respectively, and the superscript n, the time index. Linear interpolation was usedto determine the flow variables not stored in the computer but required in Eqs. [3b] and [4b]; thus
The pressurefield pi, j was to be determined such that the velocity components tPl and v”+l obtained from Eqs. [3a] and [4a] satisfied the continuity equation (12) [71 The pressureand advance-time velocity were then solved simultaneously in an iterative procedure. Upon introducing a relaxation parameter 6r, the I + 1 iterate of the pressurefield was obtained from the velocity field, pi;; =p;, i - sr(Dy.‘i’)‘.
PI
The new pressurefield was substituted into Eqs. [3a] and [4a] to obtain the I + 1 iterate of the velocity components. The new velocity components were then used in Eq. [7] to compute the divergence. This iterative processcontinued until the changes in the iterates pf, j, (uf,‘:, j)’ and (v;,‘&)’ fell below the convergence criteria. To reduce computation time, the advection and viscous terms stored in q and 5 were updated only at the beginning of a time cycle. When a steady-state solution was
172
AU
AND
GREENFIELD
reached, the shearing stress and normal stress values were determined velocity components according to Eqs. [9] and [lo], respectively :
r,, = (2/Re) (&+3x).
from the
UOI
By definition, the stream function values were then computed with Eq. [ 1I] and the vorticity values with Eq. [12] :
0 = (avjax) - (au/i+).
[121
A fully developed velocity profile for channel blood flow was assumed at the inflow boundary. The pressure value at the inlet to the aortic region was determined upon the basis of the inflow velocity value according to Eqs. [7] and [8]. In a confined flow, there are an infinite number of solutions for the pressure field, differing only by an additive constant. It was therefore necessary to establish a reference value in the pressure field; thus, zero pressure value was specified at the outflow boundary. The outflow velocity components were assigned the corresponding velocity values immediately upstream. This approximation is valid only if the local variation in the velocity field is small. Again, for a particular blood flow velocity, this condition is satisfied only if the outflow boundary is far enough downstream from the disc that causes velocity profile deviation from a Poiseuille-flow profile. For a given number of mesh points downstream from the disc, this outflow boundary treatment, in essence, limits the range of Reynolds numbers for which the solution is applicable. METHOD
The solution of the finite difference equations requires a computer with sufficient memory capacity and an efficient algorithm. Since the accuracy of the approximation is dependent on the number of data points chosen to represent the solution space, the numerical scheme generally produces a vast amount of data. In this study the approximation to the flow in the ventricular-aortic root region is represented by approximately 2000 data values for each of the three variables. The analysis of data becomes a very difficult and time consuming task. With the employment of computer graphics, however, the data can be presented in pictorial form almost instantly. Not only is the data analysis simplified, but the graphics representation provides a close resemblance to the physical system under investigation. The computer equipment employed consists of a PDP-10 time-sharing computer, magnetic disk units, an IMLAC graphics display unit with console keyboard, a Tektronix storage display unit, and a teletype. Essential parameters, such as the Reynolds number, the time steps, and the convergence criteria, are fed into the
PROSTHETIC HEARTVALVE
SIMULATIONS
173
computer via the teletype at the initiation of a simulation run. As the computer proceeds to solve the finite difference equations iteratively, the change in the basic variables is monitored for convergence. The changes for all three variables between time steps are printed on the teletype to provide a permanent record. Since a typical simulation may require hours of computation within the time-sharing environment, it is necessary to store the intermediate solution onto disk for continued computation at a later time. This also serves to protect against total loss of data due to system failure. The intermediate solution is copied onto disk upon specific interrupt signals generated from the teletype or at the completion of the prescribed number of iterations. The stored data can be processed by the graphics routines and the results displayed on the screen for analysis. The velocity vector plot is generally used when examining the intermediate results, because such a display can be generated with the least amount of processing time. The vector plot has been found to be very helpful in pointing out program errors during the debugging stage and in signaling numerical instability during simulation. In this study the solution was very sensitive to the change in time step and the relaxation parameter. The solution is considered to have reached a steady state when the maximum change in velocity value between time steps is less than 1 ‘A. A different computer program is initiated to generate stream function, vorticity, and stress solutions from the three primary variables. These available data are then displayed on the screen in pictorial form. The complexity of the display necessitates the employment of a storage display unit. A camera mounted on the display unit records the images on film for future reproduction purposes. RESULTS
The conversion of the numerical values into picture forms was automated within the computer program and results were placed into contour map sequences. The contour map is a common form of data representation in field studies. A family of constant value lines can, in general, provide adequate information about the solution space. However, depending on the contour level selected, significant information can be concealed within a contour map; for example, additional contour levels have to be specified in the stream function display to bring out the eddy formation distal to the sewing ring. Each of the other contour maps shown herein is drawn with 10 equally spaced contour levels which encompass the range of the function values. The positive and negative values are indistinguishable when contours cannot be conveniently labeled, as in the case of computer-generated displays. To supplement the contour map, a program was designed to generate a surface representation of the solution space. The surface shown in the isometric plot represents the function values over the entire region of interest. For clarity, the hidden portion of the surface is eliminated. This very often results in the omission of negative function values. With the surface plot contour values displayed on the
174
AU
AND
GREENFIELD
screen can be assigned to the contour lines. The combination of the contour map and the surface plot can provide sufficient information about the solution to allow detailed analysis. The pressure drop across a heart valve is an important measure of the heart valve performance. In general, lower pressure drop is preferred, as high pressure drop
REYNOLDS
NO.
= 100.0
MAX = CONTOUR
-7.63975
0.68765 9.01544 FIG.
-5.97423
2.35336 10.68096
12.3465 LEUELS
-4.30671 4.01906
-2.64319 5.69440
IIIN
=
-9.3093 -0.97767 7.34992
5. Computer-formed pressure field about the tilting disc within the aortic root region.
across the heart valve means more energy is required by the heart to pump the blood The high pressure gradient across the valve also means higher blood flow velocity, which in turn causes high shear rate. Figure 5 represents the pressure distribution for the simulation model under study. The surface plot provides the detailed pressure variation across the tilting-disc valve which is not generally obtainable from mock circulation experiments. In most experimental studies, the pressure drop is evaluated with the pressure values measured at two arbitrarily chosen points
into circulation.
PROSTHETIC
HEART
VALVE
175
SIMULATIONS
across the valve. It is obvious, from Fig. 5, that correct pressure drops can only be calculated if the measuring devices are placed at appropriate locations. The resulting display, obtained from the numerical simulation can provide guidelines for placement of the pressure measuring probes in the laboratory simulation of a more accurate heart valve model. The pressure data given in Fig. 5 provide a valuable
-------_-__--______-__________________ REYNOLDS
NO.
= 100.0
-0.14648 2.32215
0.34725 2.81587
FIG. 6. Computer-formed the aortic root region.
horizontal
MAX
=
4.7908
CONTOUR LEUELS 0.84097 3.30960
tllN
1.33470 3.60333
=
-0.6402
1.82842 4.29705
velocity component of flow about the tilting disc within
measure of the effectiveness of the tilting-disc heart valve when comparison is made with other heart valve models studied under similar conditions. The two velocity components, u and v, are displayed in Figs. 6 and 7, respectively. High horizontal velocity is indicated near the two ends of the disc as a result of the narrowing of the flow path between the tilted disc and the channel wall. As the fluid passes through the minor orifice it is directed by the sewing ring to move vertically
176
AU AND GREENFIELD
toward the underside of the disc. The change in flow direction is indicated in Fig. 7 by the peak vertical velocity component near the sewing ring. Although sufficient information is contained in the two velocity component displays to determine the flow about this particular tilting-disc configuration, the vector
REYNOLDS
-2.459111 0.23479 FIG. 7. Computer-formed aortic root region.
NO.
= ~00.0
-1.92032 0.77357
IIAX
=
2.9297
CONTOUR LEUELS -1.39154 1.31235
-0.64276 1.95112
IIIN
=
-2.9979
-0.30399 2.39990
vertical velocity component of flow about the tilting disc within the
sum of the velocity components provides a better visual perception of the flow variation. In a velocity vector display the length and direction of the vector at a point represent the magnitude and direction of the velocity. The short vectors shown in Fig. 8 indicate a stagnation region adjacent to the channel wall downstream of the disc. Due to the small-scale factor used to draw the vectors, the eddy formation behind the sewing ring is not as clearly represented as in the streamline plot.
AU AND GREENFIELD
178
A fourth-order Ralston-Runge-Kutta numerical integration technique was employed to evaluate the stream function values with Eq. [l I]. With the velocity defined asthe gradient of the stream function, the closely spacedstreamlinesshown in Fig. 9 represent high-velocity regions. A closed streamline contour representsan eddy formation. The contribution of eddiesin blood tlow to the thrombus formation will be discussed.
---------__ REYNOLDS
-95.49675
47.10165 FIG.
10. Computer-formed
NO.
.._ ._____---___-------= 100.0
-66.97707
75.62133
HAX
= 189.7001
CONTOUR LEVELS -36.45739
104.14101
,,IN
-9.93771
132.66069
--=-124.0164 1a.59197
161.10037
vorticity field about the tilting disc within the aortic root region.
The vorticity function may be considered asa measureof fluid rotation. The signs associatedwith the vorticity values indicate the direction of rotation. In Fig. 10 the peak vorticity values are shown near the solid boundaries. In general, the vorticity contour maps give insight into the interaction of the flow with the boundaries. The shearstressdistribution and the normal stressdistribution within the fluid are presentedin Figs. 11and 12, respectively. The signson the shearstressvalues indicate
PROSTHETIC
REYNOLDS
NO.
= 100.0
HEART
NAX CONTOW?
-1.29109 0.04525 FIG.
11. Computer-formed
-1.02392 0.91252
VALVE
1.3916
=
179
SIMULATIONS
nlN
=
-1.5504
LEVELS
-0.75955 0.57979
-0.40529 0.94709
-0.2aao2 1.11433
shear stress field about the tilting disc within the aortic root region.
the direction of shear, while the signson the normal stressvalues distinguish between the tensive and compressive stresses.These two figures clearly show that the fluid passingthrough the small aperture of the valve experiencesthe most severestress. DISCUSSION
If hemodynamic theories of thrombogenesisare establishedon valid grounds, the results of the blood flow simulation through the tilting-disc heart valve model, although simplified, should bear evidence to support the theories. In order to attempt such verification, it seemsreasonableto inspect the results shown in Figs. 5 through 12. The pressuredistribution is smooth everywhere except at the minor aperture of the valve. The abrupt changein pressureresults in higher velocity which, in turn, causes higher stress.The results shown in Figs. 6, 7, and 8 are summarized in Fig. 9. The
180
AU
AND
GREENFIELD
point of interest in this stream function display is the eddy formation distal to the sewing ring and the less intense eddy further downstream. The vorticity function shows regions of severe fluid-wall interaction at both ends of the disc. Since both the vorticity and shear stress are functions of velocity gradients, some similarity is expected. The shear stress field indicates high shear stress about the upper vascular
REYNOLDS
-1.56791 -0.26877
FIG. 12. Computer-formed
NO.
= 100.0
-1.30808 -0.00894
flax
=
1.0304
CONTOUR LEUELS -1.04826 0.25089
-0.78843 0.51072
MIN
=
-1.8277
-0.52160 0.77055
normal stress field about the tilting disc within the aortic root region.
wall and the two ends of the disc, with stressconcentration at the small aperture of the valve. The normal stressdistribution also shows a high compressive stressfield below the disc. From the foregoing observations, it can be expected that any results of hemodynamically induced thrombogenesiswill occur on the underside of the tilting disc. A study of Figs. 2 and 3 confirms such expectations. A theory for the mechanics of thrombus formation is not possible, of course, without a full understanding of the chemical and biological effects which undoubtedly play important roles. However,
PROSTHETIC
REARTVALVE
SIMULATIONS
181
based on the results of this study perhaps a hypothesis for thrombus formation, although tenuous, might be stated. As red blood cells are destroyed or damaged, becauseof imposed high stress,and are subsequently entrapped in the stasisregion near the vascular wall, the contents of the red cells trigger the extrinsic clotting reaction. The continuous deposition of degeneratedred blood cellsat that site could lead to a thrombus. If the struts were included in the simulation, thrombus growth from the sewing ring position outward might duplicate that seenin Fig. 2. Information on red blood cell destruction due to fluid stresshas been reported in several publications (13,14), and a detailed discussionconcerning the mechanicsof thrombus growth on surfaceshas been given by Blackshear (15). CONCLUSION
To reduce thrombotic stenosis,a heart valve should be designedto minimize stress and eddy formation in the flow field. The natural heart valve with its unobstructed central flow feature possesses these desirable properties. One designobjective of the tilting-disc valve was to provide a similar central flow feature by reducing the flow obstruction. While the obstruction to the flow is undoubtedly less,the disturbance introduced into the blood stream by the tilting disc may be more than that of the ball or the disc prostheses.The stream function display shown in Fig. 9 indicates that the disc affects the flow up to two aortic diameters downstream from the aortic root, while previous analyses on the ball and disc valves at similar Reynolds numbers have shown only minor flow disturbance distal to the valve (26). Since the asymmetrical property of the tilting-disc valve is the causeof the added flow disturbance, the results of this study suggesta reemphasisof symmetrical flow patterns in future heart valve design. It is believed that blood flow experimentation by computer can be useful in designing an optimal heart valve. To reduce or eliminate eddies, a heart valve may be designedsuch that the flow is directed into areaswhich would otherwise contain or produce large eddies.The intensity of stressconcentrations in the flow field may be reduced by streamlining the heart valve components. The use of computer simulation in such a design processwould complement laboratory experimentation and might prove to be both lesscostly and lesstime consuming. In order to vary the shape of the heart valve components, however, the irregular boundary problems discussed in this report must be overcome. The finite element approach could be the solution to this problem, and is presently a topic of research. ACKNOWLEDGMENTS
Appreciationisgivento the editorsof theJournal of Curdiouascular Surgery (Turino) for permissionto publishFig. 2, andto theeditorsof the Transactions of the American Society of Artificial Internal Organs for permission to publishFig. 3.Thanksalsoio Dr. WillemKolff for his continuing supportof thisproject. 7*
182
AU AND GREENFIELD REFERENCES
I. GREENFIELD, H. Special address: Studies of medical phenomena by computer graphics. Trans. Amer. Sot. Artif. Int. Organs l&607-61 2 (1972). 2. GREENFIELD, H., AND KOLFF, W. The prosthetic heart valve and computer graphics. J. Am. Med. Assoc. 219,69-74 (1972). 3. Au, A., AND GREENFIELD, H. Computer graphics analysis of stresses in blood flow through a prosthetic heart valve. Comput. Biol. Med. 4,279-291 (1975). 4. SANDBERG, L., REEMTSMA, K., AND GREENFIELD, H. Some theoretical aspects of vascular degeneration. Amer. J. Surg. 119, 543-552 (1970). 5. GREENFIELD, H., VICKERS, D., SUTHERLAND, I., KOLFF, W., AND REEMTSMA, K. Moving computer graphics images seen from inside the vascular system. Trans. Amer. Sot. Artif. Int. Organs 17,381-385 (1971). 6. CANNON, T., AND GREENFIELD, H. Adaptability of computer graphics to studies of atherosclerosis pathogenesis. J. Assoc. Sot. Adv. Med. Instrum. 6, 250-255 (1972). 7. GREENFIELD, H., Au, A., KELSEY, S., AND KOLFF, W. Simulation of assumed detriments to prosthetic heart function. Trans. Amer. Sot. Artif. Int. Organs 20,673-679 (1974). 8. HALLMAN, G. L., MESSMER, B. J., ELKADI, A., VON DER EMDE, K., AND COOLEY, D. A. Clinical experience with the Wada-Cutter cardiac valve prosthesis. Ann. Thordc. Surg. 10,9-19 (1970). 9. COKKINOS, D. V., VORIDIS, V., BAKOULAS, B., THEODO~~IOU, A., AND SKALKEAS, G. D. Thrombosis of two high-flow prosthetic valves. J. Thoruc. Curdiovusc. Surg. 62,947-949 (1971). 10. BOZER, A. Y., AND KARAMEHETOGLU, A. Thrombosis encountered with Bjork-Shiley prosthesis. J. Cardiovasc. Surg. (Turino) 13,141-143 (1972). Il. WELCH, J. E., HARLOW, F. H., SHANNON, J. P., AND DALY, B. J. “The MAC Method.” Los Alamos Scientific Laboratory Report LA-3425,I966. 12. VIECELLI, J. A. “Computing Method for Incompressible Flows Bounded by Moving Walls.” Lawrence Radiation Laboratory Report UCRL-728 15, 1970. 13. FORSTROM, R. J., BLACKSHEAR, P. L., JR., KESHAVIAH, P.. AND DORMAN, F. D. Fluid dynamic lysis of red cells. Chem. Eng. Prog. Symp. Ser. 67,69-74 (197; ). 24. NEVARIL, C. G., LUNCH, E. C., ALFREY, C. P., AND HELLUMS, J. D. Erythrocyte damage and destruction induced by shearing stress. J. Lab. Clin. Med. 71,784-790 (1968). 15. BLACKSHEAR, P. L., JR., AND WATTERS, C. Observation of red blood cells hitting solid walls. Chem. Eng. Prog. Symp. Ser. 67,60-68 (1971). 16. GREENFIELD, H., AND DEBRY, R. “An Application of Computer Graphics: Two Concurrent Investigations within the Medical Field.” Technical Report UTEC-CSc-71-I 15. Computer Science Department, University of Utah, 1971.