- Email: [email protected]
A numerical simulation of mechanical heart valve closure fluid dynamics
Journal of Biomechanics 35 (2002) 881–892
A numerical simulation of mechanical heart valve closure fluid dynamics Yong G. Laia, Krishnan B. Chandrana,b,*, Jack Lemmonc b
a IIHR-Hydroscience and Engineering, The University of Iowa, Iowa City, IA 52242-1527, USA Department of Biomedical Engineering, College of Engineering, The University of Iowa, 1402 SC, Iowa City, IA, USA c Medtronic Heart Valves Division, Minneapolis, MN, USA
Accepted 5 March 2002
Abstract A computational fluid dynamics model for the analysis of the bileaflet mechanical heart valve closure process is presented. The objective of the study is to demonstrate the ability of the numerical model to simulate the leaflet motion during the closing phase in order to investigate the closure fluid dynamics and to evaluate the effect of alterations in the leaflet tip geometry. The model has been applied to six different combinations of the leaflet tip geometry and the gap width between the leaflet tip and the housing. The results show that the negative pressure quickly develops on the atrial side of the leaflet tip. The pressure becomes more negative as the leaflet closure progresses and the lowest pressure is reached before the leaflet comes to a stop in the closed position. The flow dynamics at the instant of valve closure is strongly dependent on the leaflet velocity during the closing phase. Decrease of the tip velocity by a factor of three in the last four degrees of leaflet motion leads to a 50% reduction in the negative pressure magnitude. r 2002 Elsevier Science Ltd. All rights reserved. Keywords: Negative pressure transients; Cavitation initiation; Computational flow simulation; Leaflet closing velocity; Clearance flow dynamics
1. Introduction Structural failure of a model of bi-leaflet mechanical heart valves (MHV) implanted in patients has been reported previously in the literature (Klepetko et al., 1989; Kumar et al., 1989; Dimitri and Williams, 1990). Failure due to cavitation type of damage was reported as the mechanism that resulted in the structural failure (Kafesjian et al., 1994). These reports motivated several laboratories to study the dynamics of mechanical valve closure (Garrison et al., 1994; Graf et al., 1992, 1994; Zapanta et al., 1994; Wu et al., 1994; Chandran et al., 1994; Lee et al., 1994; Guo et al., 1994). These studies demonstrated that the fluid dynamics in the vicinity of the valve leaflet and the housing at the instant of valve closure may lead to large negative pressure transients, which may induce cavitation initiation and subsequent *Corresponding author. Department of Biomedical Engineering, College of Engineering, University of Iowa, 1402 SC, Iowa City, IA 52242-1527, USA. Tel.: +1-319-335-5640; fax: +1-319-335-5631. E-mail address: [email protected] (K.B. Chandran).
collapse. When the vapor bubbles collapse close to a wall, the wall boundary is impacted with a very high pressure that may result in stresses beyond the elastic limit. Bubble collapse may also result in hemolysis and platelet activation leading to thrombosis. Mechanisms attributed to cavitation include: water hammer, squeeze flow, tension waves, valve rebound, and vortex formation (Kini et al., 2000). For example, squeeze flow, between the occluder and the valve stop, was postulated as one of the causes for cavitation (Bluestein et al., 1994; Makhijani et al., 1994). In addition, presence of large negative pressure transients on the atrial side and the corresponding positive pressure transients on the ventricular side have also been demonstrated in the in vitro studies with the valve closure in the mitral position (Chandran et al., 1994). The large pressure gradients, present for a short duration during the closing phase, can induce large velocity flow through the clearance region between the leaflet tip and the valve housing. The wall shear stresses in this region can be relatively high during the closing phase and result in platelet activation inducing thrombus initiation. Thrombus deposition and associated breaking away of thrombus to form emboli
0021-9290/02/$ - see front matter r 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 1 - 9 2 9 0 ( 0 2 ) 0 0 0 5 6 - 8
882
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
that lodge in small vessels distally is also a major problem with implanted MHV. Experimental studies on the valve closing dynamics are relatively expensive and only limited amount of data can be obtained because of practical limitations in acquiring data very close to the leaflets and valve housing with the leaflets in motion. A powerful alternative is the use of computational fluid dynamic (CFD) simulation, which is capable of obtaining very detailed flow dynamics at desired locations close to the valve. Once the results of such simulations are validated with experimental data, CFD analysis can provide insightful information on the detailed flow dynamics and flow-induced stresses in the vicinity of the leaflets during valve closure. Parametric studies with geometrical changes of the leaflets and the housing as well as inflow and outflow conditions can be carried out relatively inexpensively. Effect of several design changes can be analyzed in detail with such simulations used for selecting an optimal design before prototypes are developed for further in vitro and in vivo analysis. The literature on computational studies is quite limited possibly due to the difficulties associated with the combined presence of complex MHV geometry, flow unsteadiness, and leaflet motion in simulating valve function. Several studies have been reported on the computation of wall shear stresses in the clearance region with the leaflets of the mechanical valve in the closed position (Haddad, 1990; Reif, 1991; Lee and Chandran, 1995). These studies neglected the effect of leaflet motion and also the large pressure gradients generated due to the closing dynamics demonstrated in the in vitro studies. Bluestein et al. (1994) and Makhijani et al. (1994) simulated the motion of the leaflets during the last few degrees of motion during valve closure to study the fluid being squeezed between the leaflet and the seat stop and the results demonstrated fluid velocities of about 29 m/s. Only recently, simulations coupling the fluid flow with the leaflet motion from the fully open to the closed position have been reported (e.g., Hsu et al., 1999; Aluri and Chandran, 2001) using the prescribed valve motion in a simplified flow domain. In this study, a rigorously developed computational method, capable of incorporating the leaflet motion, is employed to simulate the flow dynamics within an MHV during the period of leaflet motion from the full open position to complete closure. This study represents a first step towards a comprehensive predictive simulation of the MHV operation during the entire cardiac cycle. The major objective of the study is to demonstrate the viability of the numerical model. In particular, the mesh re-generation procedure is developed to provide a simple yet accurate mesh representation of the flow field. Results are presented on the details of the flow dynamics in the vicinity of the leaflets during valve closure with variations in the leaflet tip geometry.
2. Numerical method The numerical method used for the present study is based on the comprehensive CFD code, U2RANS, developed by the authors. The code is based on unstructured grid data structure using arbitrarily shaped elements and has been validated and applied to many fluid flows as described and demonstrated by Lai (1997, 2000). Detailed technical information of the code has been included in those publications and is not repeated here. One of the major difficulties in simulating flows within an MHV is the ability to model the leaflet motion. The bi-leaflet valve leaflets rotate by an angle of about 641 from the fully open position to valve closure. Such a large displacement presents a major challenge for the numerical simulation. In the past, two numerical methods have been very popular for the moving body flows. In the Eulerian method, the computational mesh is fixed in space while the motion of the boundary is incorporated through ‘‘activation’’ or ‘‘deactivation’’ of the computational mesh. Such an approach, e.g., was used by Amsden et al. (1992) to calculate the unsteady scavenging flow within a two-stroke engine. It was also used by Alipour et al. (1996) to simulate the laryngeal flow in a glottal model. In the arbitrary Lagrangian– Eulerian (ALE) method, the computational mesh is allowed to move in space in an arbitrary way to conform to the body motion. The ALE method, was used by Demirdzic and Peric (1990) to calculate a two-dimensional (2D) channel flow with moving indentation, and was extended by Lai and Przekwas (1994) to 3D flows for a number of applications. The two approaches have their own advantages and disadvantages. In general, the Eulerian method has the advantage of using a fixed mesh so that no new grid has to be re-generated at each time step. However, the Eulerian method requires identification of activated or deactivated portion of the mesh; moreover, re-interpolation of flow variables is needed when computational mesh is activated or deactivated to implement the boundary motion. The re-interpolation procedure is quite expensive and often causes degradation of accuracy and conservation property. On the other hand, ALE method allows the mesh to move freely to conform to the boundary motion. Particularly, the moving boundary conditions can be incorporated accurately with the ALE method. The major disadvantage of the ALE method is the need to incorporate a mechanism to regenerate the mesh to conform to the leaflet motion, which may become quite difficult for some applications. In this study, an efficient mesh regeneration procedure is developed, which also has the benefit of maintaining high-resolution mesh when the leaflet is near and at the closed position. A detailed description of the ALE method was presented by Lai and Przekwas (1994),
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
883
where several validation cases have also been included. This simulation model was also used for several numerical studies for the MHVs in the past (e.g., Makhijani et al., 1994,1997). Recently, validation studies were also conducted using the ALE method with the U2RANS algorithm and results were reported by Lai and Chandran (2000). The numerical model used in this study has been extensively validated. The details of the moving mesh generation procedure and its application to study the MHV closure dynamics is briefly described below.
3. Mechanical heart valve closure simulation 3.1. Simulation description Typical bi-leaflet valves have two semi-circular leaflets and a circular housing. In this study, the leaflet closure fluid dynamics is simulated by assuming that the flow is 2D, i.e., only the geometry and flow on the symmetry plane of the valve are simulated. The 2D assumption can be justified based on the following reasons. The aim of the study is to demonstrate the ability of the numerical model to simulate the full motion of the leaflet from the open position to closure. Representation of the actual semi-circular geometry of the leaflet is less critical for the same. Moreover, the 2D approximation is a reasonable assumption to represent the flow on the symmetry plane of the valve. The flow dynamics near the clearance region between the leaflet tip and the valve housing is of major interest for cavitation studies and this region is influenced the most by the leaflet tip velocity. Considering that the linear velocity of the leaflet tip is the largest on the symmetry plane and the radius of the leaflet is relatively large compared to the clearance gap, a 2D simulation should yield reasonable results. Aluri and Chandran (2001) compared the results of a 2D and a corresponding 3D simulation with the leaflets held stationary in the fully closed position and the results on wall shear stresses in the clearance region differed by o10%. The valve closure fluid dynamics is influenced by a number of parameters including the leaflet tip geometry and the gap width between the tip of the leaflet and the valve housing. To examine such geometrical effects, a total of six configurations are simulated in this study representing two different leaflet tip angles and three clearance gaps as illustrated in Fig. 1. The two leaflet tip angles are 01 and 21, and the three gaps used are 0.000 , 0.00100 , and 0.00200 , respectively. For ease of presentation, each configuration will be identified by its leaflet tip angle and the gap between the tip and the housing. For example, case T2-2 refers to the configuration with 21 tip angle and 0.00200 gap, while T0-0 represents 01 tip angle and 000 gap. It is noted that only half of the valve is
Fig. 1. Detailed schematic of the leaflet tip geometry.
simulated taking advantage of the symmetry of the bileaflet valve geometry. The unsteady dynamic closure simulation is performed with the moving grid ALE method so that the leaflet motion is fully incorporated. The simulation starts from the fully open position of the leaflet and continues in time until the leaflet moves to its fully closed position. The pressure at the ventricular side (the outflow side of the valve in the mitral position) is ramped up from 0–120 mm Hg at a rate of 2000 mm Hg/ s (Chandran et al., 1994). The pressure on the atrial side of the valve is held constant at 0 mm Hg during the simulation period. The leaflet motion, i.e., the rotation angle versus time data, is taken from the experimental measurement performed by the Medtronic engineers for the simulated MHV and it is displayed in Fig. 2. This same leaflet motion data is used for all dynamic simulation cases except where noted. In the previous simulation by Aluri and Chandran (2001), leaflet velocity was curve-fitted based on the data of Bluestein et al. (1994). It is believed that the present approach is more accurate as leaflet motion is unique to each design of the heart valve and dependent on the local environment. It is noted that the use of this measured leaflet motion data can be used for even the 2D numerical model as the 2D model intends to simulate the flow near the symmetry plane of the actual 3D geometry of the valve.
884
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
Fig. 2. Ventricular pressure data and measured leaflet rotation angle data.
3.2. Mesh generation Current model uses the ALE method so that the mesh is deforming with time as leaflet is rotating from its full open position (01) to valve closure (641). The ALE method has many advantages but it relies on the regeneration of the continuously deforming mesh capable of maintaining good mesh quality. In this study, an efficient and accurate mesh re-generation procedure is used and is briefly described below. The flow near the leaflet closure is of major interest though the simulation starts from full open position. Therefore, the mesh at the leaflet closure is generated first with more grid points clustered within and near the clearance region between the leaflet tip and valve housing. A general view of the mesh is shown in Fig. 3a while a closed-up view of the mesh in the clearance region is shown in Fig. 3b. The use of such a high-quality mesh near and at the leaflet closure guarantees the accuracy of the numerical results, particularly for the closure dynamics. Such a mesh has not been possible with many previous studies either due to the limitation of the moving grid method used or the mesh-deforming algorithm adopted. Subsequently, the meshes at the full open, as well as several other selected leaflet locations, are generated. The quality and the grid point clustering can be controlled for each mesh based on the mesh configuration and using the elliptic grid smoothers. Such mesh controls are much easier for the mesh at a fixed leaflet position. The beginning, closing, and the intermediate meshes are saved in separate files once they are generated and are used later for generating arbitrary mesh. Finally, the mesh re-generation procedure can be implemented as described below: (a) At each new time step, the leaflet rotation angle is determined based on the measured curve as shown in Fig. 2. (b) Two neighboring angles are found, on which there are meshes already generated and saved in files as
discussed before. One angle is smaller and another is larger so that an interpolation of mesh can be carried out next. (c) An initial mesh at the current time is obtained through interpolation from the two neighboring meshes. (d) The grid points on the leaflet are ‘‘rotated’’ (or, displaced) to the angle obtained in step (a) above, as interpolated angle may not be exact, in general. (e) Finally, an elliptic solver is applied to the entire mesh using the mesh in step (c) as the initial guess and the displacement at the leaflet grid points in step (d) as the boundary conditions. The elliptic solver is to solve the Laplacian equation, r2 d ¼ 0; where d is the displacement of each grid point. A ‘‘final’’ mesh is obtained after the solution of the Laplacian equation using the successive overrelaxation (SOR) method. It should be noted that the ALE method accepts any deformation of the mesh. Therefore, only a few SOR iterations are needed to solve the Laplacian equation in the above procedure. Particularly, the present method prevents too-skewed (or bad) mesh from happening as the pre-generated meshes at selected leaflet locations are recovered and used. To see the results of the above procedure, two dynamic meshes corresponding to two leaflet positions are displayed in Fig. 3c and d. The meshes for other leaflet geometry and gap combinations are similar and, therefore, they are not displayed. Finally, it should be mentioned that the current mesh consists of eight blocks with a total of 9830 grid points and 9550 cells. The study of Aluri and Chandran (2001) indicated that a mesh around 2048 grid points achieved the grid independence for a similar MHV geometry. During the course of present study, a simulation with a finer mesh indicated that the mesh density used by the present study is fine enough to be regarded as grid independent. 3.3. Simulation conditions The fluid is assumed to be incompressible, laminar, and Newtonian with the density of 1056 kg/m3 and viscosity of 0.0035 kg/m s, which is representative of human blood properties at 371C. The ventricular pressure at the inlet was raised from 0 to 120 mm Hg at a pressure rise rate described in Fig. 2 and the atrial pressure at the outlet was maintained at 0 mm Hg. The ventricular and atrial pressures prescribed above are sufficient to be used as the boundary conditions at open boundaries. At the symmetry boundary, normal velocity component is set to zero and all other quantities are extrapolated assuming zero normal gradient. The remaining boundaries are solid walls, and the standard no-slip condition is used. Note that, the leaflet rotation
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
885
Fig. 3. Computational mesh for closure dynamic simulation: (a) Mesh for Case T0-2 at full closed position; (b) Mesh close-up view in the clearance region for Case T2-2; (c) Mesh for Case T0-2 at fully open position (01); (d) Mesh for Case T0-2 with the leaflet at 27.21.
characteristic is obtained from the experimental measurement and the leaflet is fully open (01) at t ¼ 0 and fully closed (641) at about t¼ 33:156 ms:
angle relative to the full open position. The computation advances in time with a time step ranging from 1 to 10 ms until the leaflet comes to a complete stop at 641. A total of six configurations are simulated and the results are discussed below.
4. Results 4.1. Detailed results for T0-1 and T2-1 cases Simulation is carried out with the leaflet moving from the full open position to valve closure. Computation starts from a quiescent flow condition at time 0 when the ventricular pressure is 0 mm Hg and the leaflet is at 01
This subsection presents time-series results for the two leaflet tip designs at the gap of 0.00100 . Comparisons with different gap sizes are discussed in the next section.
886
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
In Fig. 4, the predicted velocity magnitude contours and flow streamlines are displayed for four different leaflet positions for both T0-1 and T2-1 cases. It is noted that the maximum velocity occurs near the leaflet tip and increases with the leaflet rotation. The highest velocity occurs at 62.831 leaflet position with the value of 13.3 m/s. This is compared to the maximum leaflet tip velocity of about 2.1 m/s at 62.831. The highest velocity is mostly located immediately downstream to the leaflet housing. One important observation of the results is that
the leaflet geometry change from 01 to 21 does not induce significant changes in the flow field. The same conclusion holds for other gap widths as well. Fig. 5 displays the predicted pressure fields for the five leaflet positions for T0-1 and T2-1 cases. The lowest pressure computed during the simulation is indicated at the bottom of each figure. It is observed that the negative pressure spot quickly develops, and the magnitude of the negative pressure continues to increase with the leaflet motion. The lowest pressure of
Fig. 4. Predicted velocity contours and flow streamlines at various times (leaflet position) during valve closure in the dynamic simulation. (a) Case T0-1; (b) Case T2-1.
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
887
Fig. 5. Predicted pressure distribution at various times (leaflet position) during valve closure in the dynamic simulation: (a) Case T0-1; (b) Case T2-1.
–68,960 Pa (517 mm Hg) is at the leaflet position of 62.81. It is found that the lowest pressure spot is not necessarily located at the leaflet tip surface. Instead, it is predicted to be behind the leaflet and close to the leaflethousing gap. The leaflet tip pressure is plotted against time (or leaflet position) in Fig. 6. It is seen that the lowest pressure at the tip reaches about 42,000 Pa (317 mm Hg). Overall, it is shown that the results do not change noticeably between the two leaflet tip geometry designs. Simulations have also been carried out for the T0-0, T0-2, T2-0, and T2-2 cases. An examination of the results shows that the overall flow pattern and pressure transient characteristics are quite similar to the results of T0-1 and T2-1 cases and those results are not presented separately.
4.2. Detailed comparison near leaflet closure This section compares results just before the valve closure for various gap widths and leaflet tip geometry, and focus on the region near the gap between the leaflet tip and the housing. Particularly, results at 62.831 leaflet position (1.171 before full closure) are displayed and compared in Figs. 7 and 8. Note that results for T0-1 and T2-1 are not displayed, as they are very similar to those of 0.000 and 0.00200 gaps. These results clearly show that the changes due to six geometric configurations are not appreciable, as discussed earlier. At the 62.831 leaflet position, the lowest pressure spot is behind and away from the atrial side of the leaflet tip. A closer examination of the results indicates that the pressure in the clearance gap between the leaflet tip and the
888
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
Fig. 6. Variation of the leaflet tip pressure during closure.
Fig. 8. Predicted pressure distribution at 62.831 leaflet position (1.171 before full closure).
Fig. 7. Predicted velocity distribution (m/s) at 62.831 leaflet position (1.171 before full closure).
Fig. 9. Pressure distribution along the gap between the leaflet tip and the housing at 62.831 leaflet position (1.171 before full closure).
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
housing is different with different configurations. In order to gain quantitative information, Fig. 9 displays the pressure distribution along the clearance gap between the leaflet tip and the housing at the 62.831 leaflet position (1.171 before full closure) for all six configurations. The pressure gradient is reduced by a change of the leaflet tip angle from 01 to 21 and also reduced by increasing the gap from 000 to 0.00200 . Reduced pressure gradient in the gap may be helpful to the local flow dynamics and indicates a potential benefit of using a 21 tip design. 4.3. Effect of changing leaflet motion The effect of changing the closing velocity of the leaflet on the flow dynamics was also evaluated in this simulation. The original leaflet closing velocity data is
889
designated as ‘‘normal’’ while a slightly modified velocity, named ‘‘modified’’, is used to study the effect. The only modification is that the time taken for the leaflet to travel from 601 to 641 is increased from 0.343 to 1.0 ms. This modification translates into a slowdown of the leaflet rotational speed by almost a factor of 3 during the last 41 of rotational motion before full closure. The computation is carried out for the 21 leaflet tip design and 0.000 gap. Fig. 10 shows the comparison of the calculated results between the Normal and Modified speed of the leaflet at the instance of 62.831 leaflet position (1.171 before full closure). It is observed that a slowdown of the leaflet reduces the maximum positive and negative pressure, as well as the velocity, dramatically. For the case simulated, the maximum positive and negative pressures are reduced by more than 50%. The
Fig. 10. Comparison of results with ‘‘normal’’ and ‘‘modified’’ leaflet closing velocity at 62.831 leaflet position (1.171 before full closure): (a) velocity at normal speed (Vmax ¼ 13:29 m=s); (b) velocity at modified speed (Vmax ¼ 9:01 m=s); (c) pressure at normal speed (Pmax ¼ 31; 814 Pa; Pmin ¼ 68; 555 Pa); (d) pressure at modified speed (Pmax ¼ 14; 279 Pa; Pmin ¼ 30; 075 Pa).
890
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
pressure distribution along the clearance gap between the leaflet tip and the housing is also displayed in Fig. 9 with the symbol ‘‘M’’. It is seen that both the pressure magnitude (positive or negative) and the pressure gradient is dramatically reduced. The above results indicate that the leaflet tip velocity is an important parameter that influences the closure flow dynamics significantly.
5. Discussion and conclusions This study has demonstrated that the CFD simulation can be used to study the complex valve closure fluid dynamics in detail and can be a tool for valve design before machining prototypes for detailed testing. The study concentrates on the simulation of the leaflet closure process in order to evaluate the effect of the design changes of the leaflet geometry on the flow dynamics. Seven dynamic simulations have been carried out with six different combinations of the leaflet tip geometry and the gap between the tip and the housing. Following conclusions can be drawn from this study: (1) The simulation of valve closure dynamics with the leaflet motion incorporated showed that the negative pressure quickly develops on the atrial side of the leaflet tip near the region between the leaflet tip and the housing. As the leaflet closure progresses, the pressure gets more negative with time and the lowest pressure is reached just before the leaflet comes to a complete stop. For the 0.00100 gap case, the lowest pressure is about 68,960 Pa (or 517 mm Hg) at 62.831 leaflet position. The leaflet tip angle change from 01 to 21 does not induce significant changes in the flow field. This finding is also observed for other gap widths. The results show that the vortex dynamics is responsible for the creation of local low-pressure spot for the MHV studied, confirming previous speculation that vortex formation is one of the mechanisms for cavitation. (2) Comparison of six geometric configurations, representing two leaflet tip angles (01 and 21) and three gap clearances between the leaflet tip and the housing (000 , 0.00100 , and 0.00200 ), indicates that the overall closure flow dynamics does not change appreciably between the configurations. The same is true for the maximum positive and negative pressures which are not significantly altered for the six configurations simulated. This finding is contrary to many previous suggestions that local tip geometry may have a large effect on the local closure dynamics and cavitation formation. However, the pressure gradient along the clearance gap between the leaflet tip and the housing is reduced appreciably when the leaflet tip angle is changed
from 01 to 21 and the gap is increased from 000 to 0.00200 . (3) The simulation also demonstrated that the leaflet closing speed has a large influence on the closure flow dynamics. For example, reducing the leaflet tip velocity in the last 41 of motion by a factor of 3, or lengthening the time taken to travel the last 41 of leaflet motion from 0.343 to 1.0 ms, leads to a reduction of both the maximum positive and negative pressures by more than 50%. The results on the presence of negative pressure transients on the atrial side of the leaflets and its dependence on the leaflet closing velocity agree with the in vitro results presented from our laboratory (Chandran et al., 1994; Chandran and Aluri, 1997). The valve leaflet closure is controlled by the rate of ventricular pressure rise dp/dt; similar to that measured in vivo. In an in vitro experimental set up, we varied the dp/dt during closure from 750 to 3000 mm Hg/s and the velocity of the leaflet at the instant of valve closure increased from 2.07 to 3.3 m/s. Correspondingly, the measured negative pressure transients also increased in intensity and duration (Chandran and Aluri, 1997). The dp/dt employed in these studies were the same as measured in vivo where mechanical valves were implanted and negative pressure transients were measured near the valve in the atrial chamber (Chandran et al., 1998). By pharmacological intervention, we were able to vary the dp/dt to a range of magnitudes and demonstrate that the same affected the pressure transients. The agreement between the current simulation and the in vitro results serve as further validation of the results of the simulation.
6. Limitations of the study There are a number of limitations in the simulations reported that must be considered in applying the results from this study. Firstly, the simulation was performed with a 2D model in order to concentrate on the model development. Even though a 2D simulation can be anticipated to reasonably predict the flow in the clearance region around the peripheral zone of the valve, the current simulation did not model the hinge mechanism that guides the motion of the leaflet. Aluri and Chandran (2001) compared the results from a 2D simulation with a corresponding 3D simulation in a quasi-steady case without the leaflet motion. The maximum wall shear stresses in the clearance region between the two cases differed by o10%. Secondly, the motion of the leaflet in the current study was specified from experimental measurements of velocity of the leaflet during valve closure. Such a specification enabled us to avoid the complex problem of leaflet–fluid
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
interaction during the valve closure. Hence in the specification of the boundary conditions, the phase difference between the ventricular pressure rise and the leaflet motion was not modeled. A more realistic simulation will be to model the leaflet motion including the solid fluid interaction rather than specifying the leaflet velocity from in vitro measurements. For the fluid–solid interaction simulation, the governing equations of motion for the leaflet about the hinge must be developed and iteratively solved with the fluid equations with the mesh updated corresponding to the induced motion of the leaflets. Finally, the current simulation was stopped at the instant of valve closure and hence, the interaction of the leaflet with valve housing and subsequent rebound has not been considered in this simulation as well. Kini et al. (2000) employed a particle image velocimetry technique to analyze the flow field near the major radius of the leaflet during valve closure. Their results suggest that the negative pressure drop exist in this region similar to the present simulation. In addition, their study also suggests that additional pressure drop is created due to the rebound of the rigid leaflet after closure. They postulate that the combination of the initial pressure drop and the rebound effect will drive the local pressure below the vapor pressure and hence generate vaporous cavitation. Hence, it is important to analyze the leaflet rebound effect as well in order to study the mechanism of cavitation with mechanical valve closure. Even though our simulations predict large negative pressures on the atrial side of the leaflets in the clearance region between the housing and the tip of the leaflets that can initiate cavitation, we did not attempt to analyze for cavitation initiation. The gap width between the leaflet tip and the housing dimensions are comparable to those of red blood cells. Since our interest was only to delineate the details of the pressure and velocity distribution in this region with valve closure, we used a Newtonian fluid simulation. Two-phase flow simulation will be required if the effect of fluid-induced stresses on the formed elements of blood, as well as initiation of vapor bubbles in the region of negative pressure transients are to be analyzed.
References Alipour, F., Fan, C., Scherer, R.C., 1996. A numerical simulation of laryngeal flow in a forced-oscillation glottal model. Journal of Comp. Speech & Language 10, 75–93. Aluri, S., Chandran, K.B., 2001. Numerical simulation of mechanical mitral heart valve closure. Annals of Biomedical Engineering 29, 665–676. Amsden, A.A., O’Rourke, P.J., Butler, T.D., 1992. Comparisons of computed and measured three-dimensional velocity fields in a motored two-stroke engine. SAE Paper-920418.
891
Bluestein, D.W., Einav, S., Hwang, N.H.C., 1994. A squeeze flow phenomenon at the closing of a bi-leaflet mechanical heart valve prosthesis. Journal of Biomechanics 27, 1369–1378. Chandran, K.B., Aluri, S., 1997. Mechanical valve closing dynamics: relationship between velocity of closing, pressure transients, and cavitation initiation. Annals of Biomedical Engineering 25, 926–938. Chandran, K.B., Lee, C.S., Chen, L.D., 1994. Pressure field in the vicinity of mechanical valve occludes at the instant of valve closure: correlation with cavitation initiation. Journal of Heart Valve Disease 3 (Suppl.I), S65–S76. Chandran, K.B., Dexter, E.U., Aluri, S., Richenbacher, W.E., 1998. Negative pressure transients with mechanical heart valve closure: correlation between in vitro and in vivo results. Annals of Biomedical Engineering 26, 546–556. Demirdzic, I., Peric, M., 1990. Finite volume method for prediction of fluid flow in arbitrarily shaped domains with moving boundaries. International Journal of Numerical Methods in Fluids 10, 771–790. Dimitri, W.R., Williams, B.T., 1990. Fracture of the duromedics mitral valve housing with leaflet escape. Journal of Cardiovascular Surgery 31, 41–46. Garrison, L.A., Lamson, T.C., Deutsch, S., Geselowitz, D.B., Gaumond, R.P., Tarbell, J.M., 1994. An in-vitro investigation of prosthetic heart valve cavitation in blood. Journal of Heart Valve Disease 3 (Suppl. I), S8–S24. Graf, T., Reul, H., Dietz, W., Wilmes, R., Rau, G., 1992. Cavitation at mechanical heart valves under simulated physiological conditions. Journal of Heart Valve Disease 1, 131–141. Graf, T., Reul, H., Detelfs, C., Wilmes, R., Rau, G., 1994. Causes and formation of cavitation in mechanical heart valves. Journal of Heart Valve Disease 3 (Suppl.I), S49–S64. Guo, G.X., Adlparvar, P., Howanec, M., Roy, J., Kafesjian, R., Kingsbury, C., 1994. Effect of structural compliance on cavitation threshold measurement of mechanical heart valves. Journal of Heart Valve Disease 3 (Suppl. I), S77–S84. Haddad, Y.A.M., 1990. A comparative study of the shear stresses induced in the leakage back flow produced by four types of heart valve prostheses. Journal of Engineering in Medicine 204, 111–114. Hsu, A.T., Yun, J.X., Hwang, J.X., 1999. Application of an unstructured grid algorithm to artificial heart valve simulations. ASAIO Journal 45, 581–585. Kafesjian, R., Howanec, M., Ward, G.D., Diep, L., Wagstaff, L.S., Rhee, R., 1994. Cavitation damage of pyrolytic carbon in mechanical heart valves. Journal of Heart Valve Disease 3, S2–S7. Kini, V., Bachmann, C., Fontaine, A., Deutsch, S., Tarbell, J.M., 2000. Flow visualization in mechanical heart valves: occluder rebound and cavitation potential. Annals of Biomedical Engineering 28, 431–441. Klepetko, W., Moritz, A., Mlzoch, G., Schurawitzki, H., Domanis, E., Wolner, E.J., 1989. Leaflet fracture in Edwards–Duromedics bileaflet valves. Journal of Thoracic Cardiovascular Surgery 97, 90–94. Kumar, N., Balasundaram, S., Rickard, M., Al Halees, N., Duran, C.M., 1989. Leaflet embolization from Duromedics valves: a report of two cases. Journal of Thoracic Cardiovascular Surgery 97, 95–97. Lai, Y.G., 1997. An unstructured grid method for a pressure-based flow and heat transfer solver. Numerical Heat Transfer, Part B 32, 267–281. Lai, Y.G., 2000. An unstructured grid arbitrarily shaped element method for fluid flow simulation. AIAA Journal 38 (9), 2246–2252. Lai, Y.G., Chandran, K.B., 2000. Two-dimensional computational simulation of mechanical heart valve closure fluid dynamics. Iowa Institute of Hydraulic Research, Report No. 285, Iowa City, IA. Lai, Y.G., Przekwas, A.J., 1994. A finite-volume method for fluid flow simulations with moving boundaries. Computational Fluid Dynamics 2, 19–40.
892
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
Lee, C.S., Chandran, K.B., 1995. Numerical simulation of instantaneous backflow through central clearance of bi-leaflet mechanical heart valve at closure: shear stresses and pressure fields. Medical and Biological Engineering and Computing 33 (3), 257–263. Lee, C.S., Chandran, K.B., Chen, L.D., 1994. Cavitation dynamics of mechanical heart valve prostheses. Artificial Organs 18, 758–767. Makhijani, V.B., Yang, H.Q., Singhal, A.K., Hwang, N.H.C., 1994. Computational analysis of MHV cavitation: effects of leaflet squeezing and rebound. Journal of Heart Valve Disease 3 (Suppl. 1), S35–S48. Makhijani, V.B., Yang, H.Q., Dionne, P.J., Thubrikar, M.J., 1997. Three-dimensional coupled fluid-structure simulation of pericardial
bioprosthetic aortic valve functioning. ASAIO Journal 43 (5), M387–M392. Reif, T.H., 1991. A numerical analysis of backflow between the leaflets of St Jude medical cardiac valve prostheses. Journal of Biomechanics 24 (8), 733–741. Wu, Z.J., Wang, Y., Hwang, N.H.C., 1994. Occluder closing behavior: a key factor in mechanical heart valve cavitation. Journal of Heart Valve Disease 3 (Suppl. I), S25–S34. Zapanta, C.M., Liszka Jr., E.G., Lamson, T.C., Stinebring, D.R., Deutsch, S., Geselowiz, D.B., Tarbell, J.B., 1994. A method for real-time in vitro observation of cavitation on prosthetic heart valves. ASME Journal of Biomechanical Engineering 116, 460–468.
A numerical simulation of mechanical heart valve closure fluid dynamics Yong G. Laia, Krishnan B. Chandrana,b,*, Jack Lemmonc b
a IIHR-Hydroscience and Engineering, The University of Iowa, Iowa City, IA 52242-1527, USA Department of Biomedical Engineering, College of Engineering, The University of Iowa, 1402 SC, Iowa City, IA, USA c Medtronic Heart Valves Division, Minneapolis, MN, USA
Accepted 5 March 2002
Abstract A computational fluid dynamics model for the analysis of the bileaflet mechanical heart valve closure process is presented. The objective of the study is to demonstrate the ability of the numerical model to simulate the leaflet motion during the closing phase in order to investigate the closure fluid dynamics and to evaluate the effect of alterations in the leaflet tip geometry. The model has been applied to six different combinations of the leaflet tip geometry and the gap width between the leaflet tip and the housing. The results show that the negative pressure quickly develops on the atrial side of the leaflet tip. The pressure becomes more negative as the leaflet closure progresses and the lowest pressure is reached before the leaflet comes to a stop in the closed position. The flow dynamics at the instant of valve closure is strongly dependent on the leaflet velocity during the closing phase. Decrease of the tip velocity by a factor of three in the last four degrees of leaflet motion leads to a 50% reduction in the negative pressure magnitude. r 2002 Elsevier Science Ltd. All rights reserved. Keywords: Negative pressure transients; Cavitation initiation; Computational flow simulation; Leaflet closing velocity; Clearance flow dynamics
1. Introduction Structural failure of a model of bi-leaflet mechanical heart valves (MHV) implanted in patients has been reported previously in the literature (Klepetko et al., 1989; Kumar et al., 1989; Dimitri and Williams, 1990). Failure due to cavitation type of damage was reported as the mechanism that resulted in the structural failure (Kafesjian et al., 1994). These reports motivated several laboratories to study the dynamics of mechanical valve closure (Garrison et al., 1994; Graf et al., 1992, 1994; Zapanta et al., 1994; Wu et al., 1994; Chandran et al., 1994; Lee et al., 1994; Guo et al., 1994). These studies demonstrated that the fluid dynamics in the vicinity of the valve leaflet and the housing at the instant of valve closure may lead to large negative pressure transients, which may induce cavitation initiation and subsequent *Corresponding author. Department of Biomedical Engineering, College of Engineering, University of Iowa, 1402 SC, Iowa City, IA 52242-1527, USA. Tel.: +1-319-335-5640; fax: +1-319-335-5631. E-mail address: [email protected] (K.B. Chandran).
collapse. When the vapor bubbles collapse close to a wall, the wall boundary is impacted with a very high pressure that may result in stresses beyond the elastic limit. Bubble collapse may also result in hemolysis and platelet activation leading to thrombosis. Mechanisms attributed to cavitation include: water hammer, squeeze flow, tension waves, valve rebound, and vortex formation (Kini et al., 2000). For example, squeeze flow, between the occluder and the valve stop, was postulated as one of the causes for cavitation (Bluestein et al., 1994; Makhijani et al., 1994). In addition, presence of large negative pressure transients on the atrial side and the corresponding positive pressure transients on the ventricular side have also been demonstrated in the in vitro studies with the valve closure in the mitral position (Chandran et al., 1994). The large pressure gradients, present for a short duration during the closing phase, can induce large velocity flow through the clearance region between the leaflet tip and the valve housing. The wall shear stresses in this region can be relatively high during the closing phase and result in platelet activation inducing thrombus initiation. Thrombus deposition and associated breaking away of thrombus to form emboli
0021-9290/02/$ - see front matter r 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 1 - 9 2 9 0 ( 0 2 ) 0 0 0 5 6 - 8
882
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
that lodge in small vessels distally is also a major problem with implanted MHV. Experimental studies on the valve closing dynamics are relatively expensive and only limited amount of data can be obtained because of practical limitations in acquiring data very close to the leaflets and valve housing with the leaflets in motion. A powerful alternative is the use of computational fluid dynamic (CFD) simulation, which is capable of obtaining very detailed flow dynamics at desired locations close to the valve. Once the results of such simulations are validated with experimental data, CFD analysis can provide insightful information on the detailed flow dynamics and flow-induced stresses in the vicinity of the leaflets during valve closure. Parametric studies with geometrical changes of the leaflets and the housing as well as inflow and outflow conditions can be carried out relatively inexpensively. Effect of several design changes can be analyzed in detail with such simulations used for selecting an optimal design before prototypes are developed for further in vitro and in vivo analysis. The literature on computational studies is quite limited possibly due to the difficulties associated with the combined presence of complex MHV geometry, flow unsteadiness, and leaflet motion in simulating valve function. Several studies have been reported on the computation of wall shear stresses in the clearance region with the leaflets of the mechanical valve in the closed position (Haddad, 1990; Reif, 1991; Lee and Chandran, 1995). These studies neglected the effect of leaflet motion and also the large pressure gradients generated due to the closing dynamics demonstrated in the in vitro studies. Bluestein et al. (1994) and Makhijani et al. (1994) simulated the motion of the leaflets during the last few degrees of motion during valve closure to study the fluid being squeezed between the leaflet and the seat stop and the results demonstrated fluid velocities of about 29 m/s. Only recently, simulations coupling the fluid flow with the leaflet motion from the fully open to the closed position have been reported (e.g., Hsu et al., 1999; Aluri and Chandran, 2001) using the prescribed valve motion in a simplified flow domain. In this study, a rigorously developed computational method, capable of incorporating the leaflet motion, is employed to simulate the flow dynamics within an MHV during the period of leaflet motion from the full open position to complete closure. This study represents a first step towards a comprehensive predictive simulation of the MHV operation during the entire cardiac cycle. The major objective of the study is to demonstrate the viability of the numerical model. In particular, the mesh re-generation procedure is developed to provide a simple yet accurate mesh representation of the flow field. Results are presented on the details of the flow dynamics in the vicinity of the leaflets during valve closure with variations in the leaflet tip geometry.
2. Numerical method The numerical method used for the present study is based on the comprehensive CFD code, U2RANS, developed by the authors. The code is based on unstructured grid data structure using arbitrarily shaped elements and has been validated and applied to many fluid flows as described and demonstrated by Lai (1997, 2000). Detailed technical information of the code has been included in those publications and is not repeated here. One of the major difficulties in simulating flows within an MHV is the ability to model the leaflet motion. The bi-leaflet valve leaflets rotate by an angle of about 641 from the fully open position to valve closure. Such a large displacement presents a major challenge for the numerical simulation. In the past, two numerical methods have been very popular for the moving body flows. In the Eulerian method, the computational mesh is fixed in space while the motion of the boundary is incorporated through ‘‘activation’’ or ‘‘deactivation’’ of the computational mesh. Such an approach, e.g., was used by Amsden et al. (1992) to calculate the unsteady scavenging flow within a two-stroke engine. It was also used by Alipour et al. (1996) to simulate the laryngeal flow in a glottal model. In the arbitrary Lagrangian– Eulerian (ALE) method, the computational mesh is allowed to move in space in an arbitrary way to conform to the body motion. The ALE method, was used by Demirdzic and Peric (1990) to calculate a two-dimensional (2D) channel flow with moving indentation, and was extended by Lai and Przekwas (1994) to 3D flows for a number of applications. The two approaches have their own advantages and disadvantages. In general, the Eulerian method has the advantage of using a fixed mesh so that no new grid has to be re-generated at each time step. However, the Eulerian method requires identification of activated or deactivated portion of the mesh; moreover, re-interpolation of flow variables is needed when computational mesh is activated or deactivated to implement the boundary motion. The re-interpolation procedure is quite expensive and often causes degradation of accuracy and conservation property. On the other hand, ALE method allows the mesh to move freely to conform to the boundary motion. Particularly, the moving boundary conditions can be incorporated accurately with the ALE method. The major disadvantage of the ALE method is the need to incorporate a mechanism to regenerate the mesh to conform to the leaflet motion, which may become quite difficult for some applications. In this study, an efficient mesh regeneration procedure is developed, which also has the benefit of maintaining high-resolution mesh when the leaflet is near and at the closed position. A detailed description of the ALE method was presented by Lai and Przekwas (1994),
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
883
where several validation cases have also been included. This simulation model was also used for several numerical studies for the MHVs in the past (e.g., Makhijani et al., 1994,1997). Recently, validation studies were also conducted using the ALE method with the U2RANS algorithm and results were reported by Lai and Chandran (2000). The numerical model used in this study has been extensively validated. The details of the moving mesh generation procedure and its application to study the MHV closure dynamics is briefly described below.
3. Mechanical heart valve closure simulation 3.1. Simulation description Typical bi-leaflet valves have two semi-circular leaflets and a circular housing. In this study, the leaflet closure fluid dynamics is simulated by assuming that the flow is 2D, i.e., only the geometry and flow on the symmetry plane of the valve are simulated. The 2D assumption can be justified based on the following reasons. The aim of the study is to demonstrate the ability of the numerical model to simulate the full motion of the leaflet from the open position to closure. Representation of the actual semi-circular geometry of the leaflet is less critical for the same. Moreover, the 2D approximation is a reasonable assumption to represent the flow on the symmetry plane of the valve. The flow dynamics near the clearance region between the leaflet tip and the valve housing is of major interest for cavitation studies and this region is influenced the most by the leaflet tip velocity. Considering that the linear velocity of the leaflet tip is the largest on the symmetry plane and the radius of the leaflet is relatively large compared to the clearance gap, a 2D simulation should yield reasonable results. Aluri and Chandran (2001) compared the results of a 2D and a corresponding 3D simulation with the leaflets held stationary in the fully closed position and the results on wall shear stresses in the clearance region differed by o10%. The valve closure fluid dynamics is influenced by a number of parameters including the leaflet tip geometry and the gap width between the tip of the leaflet and the valve housing. To examine such geometrical effects, a total of six configurations are simulated in this study representing two different leaflet tip angles and three clearance gaps as illustrated in Fig. 1. The two leaflet tip angles are 01 and 21, and the three gaps used are 0.000 , 0.00100 , and 0.00200 , respectively. For ease of presentation, each configuration will be identified by its leaflet tip angle and the gap between the tip and the housing. For example, case T2-2 refers to the configuration with 21 tip angle and 0.00200 gap, while T0-0 represents 01 tip angle and 000 gap. It is noted that only half of the valve is
Fig. 1. Detailed schematic of the leaflet tip geometry.
simulated taking advantage of the symmetry of the bileaflet valve geometry. The unsteady dynamic closure simulation is performed with the moving grid ALE method so that the leaflet motion is fully incorporated. The simulation starts from the fully open position of the leaflet and continues in time until the leaflet moves to its fully closed position. The pressure at the ventricular side (the outflow side of the valve in the mitral position) is ramped up from 0–120 mm Hg at a rate of 2000 mm Hg/ s (Chandran et al., 1994). The pressure on the atrial side of the valve is held constant at 0 mm Hg during the simulation period. The leaflet motion, i.e., the rotation angle versus time data, is taken from the experimental measurement performed by the Medtronic engineers for the simulated MHV and it is displayed in Fig. 2. This same leaflet motion data is used for all dynamic simulation cases except where noted. In the previous simulation by Aluri and Chandran (2001), leaflet velocity was curve-fitted based on the data of Bluestein et al. (1994). It is believed that the present approach is more accurate as leaflet motion is unique to each design of the heart valve and dependent on the local environment. It is noted that the use of this measured leaflet motion data can be used for even the 2D numerical model as the 2D model intends to simulate the flow near the symmetry plane of the actual 3D geometry of the valve.
884
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
Fig. 2. Ventricular pressure data and measured leaflet rotation angle data.
3.2. Mesh generation Current model uses the ALE method so that the mesh is deforming with time as leaflet is rotating from its full open position (01) to valve closure (641). The ALE method has many advantages but it relies on the regeneration of the continuously deforming mesh capable of maintaining good mesh quality. In this study, an efficient and accurate mesh re-generation procedure is used and is briefly described below. The flow near the leaflet closure is of major interest though the simulation starts from full open position. Therefore, the mesh at the leaflet closure is generated first with more grid points clustered within and near the clearance region between the leaflet tip and valve housing. A general view of the mesh is shown in Fig. 3a while a closed-up view of the mesh in the clearance region is shown in Fig. 3b. The use of such a high-quality mesh near and at the leaflet closure guarantees the accuracy of the numerical results, particularly for the closure dynamics. Such a mesh has not been possible with many previous studies either due to the limitation of the moving grid method used or the mesh-deforming algorithm adopted. Subsequently, the meshes at the full open, as well as several other selected leaflet locations, are generated. The quality and the grid point clustering can be controlled for each mesh based on the mesh configuration and using the elliptic grid smoothers. Such mesh controls are much easier for the mesh at a fixed leaflet position. The beginning, closing, and the intermediate meshes are saved in separate files once they are generated and are used later for generating arbitrary mesh. Finally, the mesh re-generation procedure can be implemented as described below: (a) At each new time step, the leaflet rotation angle is determined based on the measured curve as shown in Fig. 2. (b) Two neighboring angles are found, on which there are meshes already generated and saved in files as
discussed before. One angle is smaller and another is larger so that an interpolation of mesh can be carried out next. (c) An initial mesh at the current time is obtained through interpolation from the two neighboring meshes. (d) The grid points on the leaflet are ‘‘rotated’’ (or, displaced) to the angle obtained in step (a) above, as interpolated angle may not be exact, in general. (e) Finally, an elliptic solver is applied to the entire mesh using the mesh in step (c) as the initial guess and the displacement at the leaflet grid points in step (d) as the boundary conditions. The elliptic solver is to solve the Laplacian equation, r2 d ¼ 0; where d is the displacement of each grid point. A ‘‘final’’ mesh is obtained after the solution of the Laplacian equation using the successive overrelaxation (SOR) method. It should be noted that the ALE method accepts any deformation of the mesh. Therefore, only a few SOR iterations are needed to solve the Laplacian equation in the above procedure. Particularly, the present method prevents too-skewed (or bad) mesh from happening as the pre-generated meshes at selected leaflet locations are recovered and used. To see the results of the above procedure, two dynamic meshes corresponding to two leaflet positions are displayed in Fig. 3c and d. The meshes for other leaflet geometry and gap combinations are similar and, therefore, they are not displayed. Finally, it should be mentioned that the current mesh consists of eight blocks with a total of 9830 grid points and 9550 cells. The study of Aluri and Chandran (2001) indicated that a mesh around 2048 grid points achieved the grid independence for a similar MHV geometry. During the course of present study, a simulation with a finer mesh indicated that the mesh density used by the present study is fine enough to be regarded as grid independent. 3.3. Simulation conditions The fluid is assumed to be incompressible, laminar, and Newtonian with the density of 1056 kg/m3 and viscosity of 0.0035 kg/m s, which is representative of human blood properties at 371C. The ventricular pressure at the inlet was raised from 0 to 120 mm Hg at a pressure rise rate described in Fig. 2 and the atrial pressure at the outlet was maintained at 0 mm Hg. The ventricular and atrial pressures prescribed above are sufficient to be used as the boundary conditions at open boundaries. At the symmetry boundary, normal velocity component is set to zero and all other quantities are extrapolated assuming zero normal gradient. The remaining boundaries are solid walls, and the standard no-slip condition is used. Note that, the leaflet rotation
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
885
Fig. 3. Computational mesh for closure dynamic simulation: (a) Mesh for Case T0-2 at full closed position; (b) Mesh close-up view in the clearance region for Case T2-2; (c) Mesh for Case T0-2 at fully open position (01); (d) Mesh for Case T0-2 with the leaflet at 27.21.
characteristic is obtained from the experimental measurement and the leaflet is fully open (01) at t ¼ 0 and fully closed (641) at about t¼ 33:156 ms:
angle relative to the full open position. The computation advances in time with a time step ranging from 1 to 10 ms until the leaflet comes to a complete stop at 641. A total of six configurations are simulated and the results are discussed below.
4. Results 4.1. Detailed results for T0-1 and T2-1 cases Simulation is carried out with the leaflet moving from the full open position to valve closure. Computation starts from a quiescent flow condition at time 0 when the ventricular pressure is 0 mm Hg and the leaflet is at 01
This subsection presents time-series results for the two leaflet tip designs at the gap of 0.00100 . Comparisons with different gap sizes are discussed in the next section.
886
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
In Fig. 4, the predicted velocity magnitude contours and flow streamlines are displayed for four different leaflet positions for both T0-1 and T2-1 cases. It is noted that the maximum velocity occurs near the leaflet tip and increases with the leaflet rotation. The highest velocity occurs at 62.831 leaflet position with the value of 13.3 m/s. This is compared to the maximum leaflet tip velocity of about 2.1 m/s at 62.831. The highest velocity is mostly located immediately downstream to the leaflet housing. One important observation of the results is that
the leaflet geometry change from 01 to 21 does not induce significant changes in the flow field. The same conclusion holds for other gap widths as well. Fig. 5 displays the predicted pressure fields for the five leaflet positions for T0-1 and T2-1 cases. The lowest pressure computed during the simulation is indicated at the bottom of each figure. It is observed that the negative pressure spot quickly develops, and the magnitude of the negative pressure continues to increase with the leaflet motion. The lowest pressure of
Fig. 4. Predicted velocity contours and flow streamlines at various times (leaflet position) during valve closure in the dynamic simulation. (a) Case T0-1; (b) Case T2-1.
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
887
Fig. 5. Predicted pressure distribution at various times (leaflet position) during valve closure in the dynamic simulation: (a) Case T0-1; (b) Case T2-1.
–68,960 Pa (517 mm Hg) is at the leaflet position of 62.81. It is found that the lowest pressure spot is not necessarily located at the leaflet tip surface. Instead, it is predicted to be behind the leaflet and close to the leaflethousing gap. The leaflet tip pressure is plotted against time (or leaflet position) in Fig. 6. It is seen that the lowest pressure at the tip reaches about 42,000 Pa (317 mm Hg). Overall, it is shown that the results do not change noticeably between the two leaflet tip geometry designs. Simulations have also been carried out for the T0-0, T0-2, T2-0, and T2-2 cases. An examination of the results shows that the overall flow pattern and pressure transient characteristics are quite similar to the results of T0-1 and T2-1 cases and those results are not presented separately.
4.2. Detailed comparison near leaflet closure This section compares results just before the valve closure for various gap widths and leaflet tip geometry, and focus on the region near the gap between the leaflet tip and the housing. Particularly, results at 62.831 leaflet position (1.171 before full closure) are displayed and compared in Figs. 7 and 8. Note that results for T0-1 and T2-1 are not displayed, as they are very similar to those of 0.000 and 0.00200 gaps. These results clearly show that the changes due to six geometric configurations are not appreciable, as discussed earlier. At the 62.831 leaflet position, the lowest pressure spot is behind and away from the atrial side of the leaflet tip. A closer examination of the results indicates that the pressure in the clearance gap between the leaflet tip and the
888
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
Fig. 6. Variation of the leaflet tip pressure during closure.
Fig. 8. Predicted pressure distribution at 62.831 leaflet position (1.171 before full closure).
Fig. 7. Predicted velocity distribution (m/s) at 62.831 leaflet position (1.171 before full closure).
Fig. 9. Pressure distribution along the gap between the leaflet tip and the housing at 62.831 leaflet position (1.171 before full closure).
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
housing is different with different configurations. In order to gain quantitative information, Fig. 9 displays the pressure distribution along the clearance gap between the leaflet tip and the housing at the 62.831 leaflet position (1.171 before full closure) for all six configurations. The pressure gradient is reduced by a change of the leaflet tip angle from 01 to 21 and also reduced by increasing the gap from 000 to 0.00200 . Reduced pressure gradient in the gap may be helpful to the local flow dynamics and indicates a potential benefit of using a 21 tip design. 4.3. Effect of changing leaflet motion The effect of changing the closing velocity of the leaflet on the flow dynamics was also evaluated in this simulation. The original leaflet closing velocity data is
889
designated as ‘‘normal’’ while a slightly modified velocity, named ‘‘modified’’, is used to study the effect. The only modification is that the time taken for the leaflet to travel from 601 to 641 is increased from 0.343 to 1.0 ms. This modification translates into a slowdown of the leaflet rotational speed by almost a factor of 3 during the last 41 of rotational motion before full closure. The computation is carried out for the 21 leaflet tip design and 0.000 gap. Fig. 10 shows the comparison of the calculated results between the Normal and Modified speed of the leaflet at the instance of 62.831 leaflet position (1.171 before full closure). It is observed that a slowdown of the leaflet reduces the maximum positive and negative pressure, as well as the velocity, dramatically. For the case simulated, the maximum positive and negative pressures are reduced by more than 50%. The
Fig. 10. Comparison of results with ‘‘normal’’ and ‘‘modified’’ leaflet closing velocity at 62.831 leaflet position (1.171 before full closure): (a) velocity at normal speed (Vmax ¼ 13:29 m=s); (b) velocity at modified speed (Vmax ¼ 9:01 m=s); (c) pressure at normal speed (Pmax ¼ 31; 814 Pa; Pmin ¼ 68; 555 Pa); (d) pressure at modified speed (Pmax ¼ 14; 279 Pa; Pmin ¼ 30; 075 Pa).
890
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
pressure distribution along the clearance gap between the leaflet tip and the housing is also displayed in Fig. 9 with the symbol ‘‘M’’. It is seen that both the pressure magnitude (positive or negative) and the pressure gradient is dramatically reduced. The above results indicate that the leaflet tip velocity is an important parameter that influences the closure flow dynamics significantly.
5. Discussion and conclusions This study has demonstrated that the CFD simulation can be used to study the complex valve closure fluid dynamics in detail and can be a tool for valve design before machining prototypes for detailed testing. The study concentrates on the simulation of the leaflet closure process in order to evaluate the effect of the design changes of the leaflet geometry on the flow dynamics. Seven dynamic simulations have been carried out with six different combinations of the leaflet tip geometry and the gap between the tip and the housing. Following conclusions can be drawn from this study: (1) The simulation of valve closure dynamics with the leaflet motion incorporated showed that the negative pressure quickly develops on the atrial side of the leaflet tip near the region between the leaflet tip and the housing. As the leaflet closure progresses, the pressure gets more negative with time and the lowest pressure is reached just before the leaflet comes to a complete stop. For the 0.00100 gap case, the lowest pressure is about 68,960 Pa (or 517 mm Hg) at 62.831 leaflet position. The leaflet tip angle change from 01 to 21 does not induce significant changes in the flow field. This finding is also observed for other gap widths. The results show that the vortex dynamics is responsible for the creation of local low-pressure spot for the MHV studied, confirming previous speculation that vortex formation is one of the mechanisms for cavitation. (2) Comparison of six geometric configurations, representing two leaflet tip angles (01 and 21) and three gap clearances between the leaflet tip and the housing (000 , 0.00100 , and 0.00200 ), indicates that the overall closure flow dynamics does not change appreciably between the configurations. The same is true for the maximum positive and negative pressures which are not significantly altered for the six configurations simulated. This finding is contrary to many previous suggestions that local tip geometry may have a large effect on the local closure dynamics and cavitation formation. However, the pressure gradient along the clearance gap between the leaflet tip and the housing is reduced appreciably when the leaflet tip angle is changed
from 01 to 21 and the gap is increased from 000 to 0.00200 . (3) The simulation also demonstrated that the leaflet closing speed has a large influence on the closure flow dynamics. For example, reducing the leaflet tip velocity in the last 41 of motion by a factor of 3, or lengthening the time taken to travel the last 41 of leaflet motion from 0.343 to 1.0 ms, leads to a reduction of both the maximum positive and negative pressures by more than 50%. The results on the presence of negative pressure transients on the atrial side of the leaflets and its dependence on the leaflet closing velocity agree with the in vitro results presented from our laboratory (Chandran et al., 1994; Chandran and Aluri, 1997). The valve leaflet closure is controlled by the rate of ventricular pressure rise dp/dt; similar to that measured in vivo. In an in vitro experimental set up, we varied the dp/dt during closure from 750 to 3000 mm Hg/s and the velocity of the leaflet at the instant of valve closure increased from 2.07 to 3.3 m/s. Correspondingly, the measured negative pressure transients also increased in intensity and duration (Chandran and Aluri, 1997). The dp/dt employed in these studies were the same as measured in vivo where mechanical valves were implanted and negative pressure transients were measured near the valve in the atrial chamber (Chandran et al., 1998). By pharmacological intervention, we were able to vary the dp/dt to a range of magnitudes and demonstrate that the same affected the pressure transients. The agreement between the current simulation and the in vitro results serve as further validation of the results of the simulation.
6. Limitations of the study There are a number of limitations in the simulations reported that must be considered in applying the results from this study. Firstly, the simulation was performed with a 2D model in order to concentrate on the model development. Even though a 2D simulation can be anticipated to reasonably predict the flow in the clearance region around the peripheral zone of the valve, the current simulation did not model the hinge mechanism that guides the motion of the leaflet. Aluri and Chandran (2001) compared the results from a 2D simulation with a corresponding 3D simulation in a quasi-steady case without the leaflet motion. The maximum wall shear stresses in the clearance region between the two cases differed by o10%. Secondly, the motion of the leaflet in the current study was specified from experimental measurements of velocity of the leaflet during valve closure. Such a specification enabled us to avoid the complex problem of leaflet–fluid
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
interaction during the valve closure. Hence in the specification of the boundary conditions, the phase difference between the ventricular pressure rise and the leaflet motion was not modeled. A more realistic simulation will be to model the leaflet motion including the solid fluid interaction rather than specifying the leaflet velocity from in vitro measurements. For the fluid–solid interaction simulation, the governing equations of motion for the leaflet about the hinge must be developed and iteratively solved with the fluid equations with the mesh updated corresponding to the induced motion of the leaflets. Finally, the current simulation was stopped at the instant of valve closure and hence, the interaction of the leaflet with valve housing and subsequent rebound has not been considered in this simulation as well. Kini et al. (2000) employed a particle image velocimetry technique to analyze the flow field near the major radius of the leaflet during valve closure. Their results suggest that the negative pressure drop exist in this region similar to the present simulation. In addition, their study also suggests that additional pressure drop is created due to the rebound of the rigid leaflet after closure. They postulate that the combination of the initial pressure drop and the rebound effect will drive the local pressure below the vapor pressure and hence generate vaporous cavitation. Hence, it is important to analyze the leaflet rebound effect as well in order to study the mechanism of cavitation with mechanical valve closure. Even though our simulations predict large negative pressures on the atrial side of the leaflets in the clearance region between the housing and the tip of the leaflets that can initiate cavitation, we did not attempt to analyze for cavitation initiation. The gap width between the leaflet tip and the housing dimensions are comparable to those of red blood cells. Since our interest was only to delineate the details of the pressure and velocity distribution in this region with valve closure, we used a Newtonian fluid simulation. Two-phase flow simulation will be required if the effect of fluid-induced stresses on the formed elements of blood, as well as initiation of vapor bubbles in the region of negative pressure transients are to be analyzed.
References Alipour, F., Fan, C., Scherer, R.C., 1996. A numerical simulation of laryngeal flow in a forced-oscillation glottal model. Journal of Comp. Speech & Language 10, 75–93. Aluri, S., Chandran, K.B., 2001. Numerical simulation of mechanical mitral heart valve closure. Annals of Biomedical Engineering 29, 665–676. Amsden, A.A., O’Rourke, P.J., Butler, T.D., 1992. Comparisons of computed and measured three-dimensional velocity fields in a motored two-stroke engine. SAE Paper-920418.
891
Bluestein, D.W., Einav, S., Hwang, N.H.C., 1994. A squeeze flow phenomenon at the closing of a bi-leaflet mechanical heart valve prosthesis. Journal of Biomechanics 27, 1369–1378. Chandran, K.B., Aluri, S., 1997. Mechanical valve closing dynamics: relationship between velocity of closing, pressure transients, and cavitation initiation. Annals of Biomedical Engineering 25, 926–938. Chandran, K.B., Lee, C.S., Chen, L.D., 1994. Pressure field in the vicinity of mechanical valve occludes at the instant of valve closure: correlation with cavitation initiation. Journal of Heart Valve Disease 3 (Suppl.I), S65–S76. Chandran, K.B., Dexter, E.U., Aluri, S., Richenbacher, W.E., 1998. Negative pressure transients with mechanical heart valve closure: correlation between in vitro and in vivo results. Annals of Biomedical Engineering 26, 546–556. Demirdzic, I., Peric, M., 1990. Finite volume method for prediction of fluid flow in arbitrarily shaped domains with moving boundaries. International Journal of Numerical Methods in Fluids 10, 771–790. Dimitri, W.R., Williams, B.T., 1990. Fracture of the duromedics mitral valve housing with leaflet escape. Journal of Cardiovascular Surgery 31, 41–46. Garrison, L.A., Lamson, T.C., Deutsch, S., Geselowitz, D.B., Gaumond, R.P., Tarbell, J.M., 1994. An in-vitro investigation of prosthetic heart valve cavitation in blood. Journal of Heart Valve Disease 3 (Suppl. I), S8–S24. Graf, T., Reul, H., Dietz, W., Wilmes, R., Rau, G., 1992. Cavitation at mechanical heart valves under simulated physiological conditions. Journal of Heart Valve Disease 1, 131–141. Graf, T., Reul, H., Detelfs, C., Wilmes, R., Rau, G., 1994. Causes and formation of cavitation in mechanical heart valves. Journal of Heart Valve Disease 3 (Suppl.I), S49–S64. Guo, G.X., Adlparvar, P., Howanec, M., Roy, J., Kafesjian, R., Kingsbury, C., 1994. Effect of structural compliance on cavitation threshold measurement of mechanical heart valves. Journal of Heart Valve Disease 3 (Suppl. I), S77–S84. Haddad, Y.A.M., 1990. A comparative study of the shear stresses induced in the leakage back flow produced by four types of heart valve prostheses. Journal of Engineering in Medicine 204, 111–114. Hsu, A.T., Yun, J.X., Hwang, J.X., 1999. Application of an unstructured grid algorithm to artificial heart valve simulations. ASAIO Journal 45, 581–585. Kafesjian, R., Howanec, M., Ward, G.D., Diep, L., Wagstaff, L.S., Rhee, R., 1994. Cavitation damage of pyrolytic carbon in mechanical heart valves. Journal of Heart Valve Disease 3, S2–S7. Kini, V., Bachmann, C., Fontaine, A., Deutsch, S., Tarbell, J.M., 2000. Flow visualization in mechanical heart valves: occluder rebound and cavitation potential. Annals of Biomedical Engineering 28, 431–441. Klepetko, W., Moritz, A., Mlzoch, G., Schurawitzki, H., Domanis, E., Wolner, E.J., 1989. Leaflet fracture in Edwards–Duromedics bileaflet valves. Journal of Thoracic Cardiovascular Surgery 97, 90–94. Kumar, N., Balasundaram, S., Rickard, M., Al Halees, N., Duran, C.M., 1989. Leaflet embolization from Duromedics valves: a report of two cases. Journal of Thoracic Cardiovascular Surgery 97, 95–97. Lai, Y.G., 1997. An unstructured grid method for a pressure-based flow and heat transfer solver. Numerical Heat Transfer, Part B 32, 267–281. Lai, Y.G., 2000. An unstructured grid arbitrarily shaped element method for fluid flow simulation. AIAA Journal 38 (9), 2246–2252. Lai, Y.G., Chandran, K.B., 2000. Two-dimensional computational simulation of mechanical heart valve closure fluid dynamics. Iowa Institute of Hydraulic Research, Report No. 285, Iowa City, IA. Lai, Y.G., Przekwas, A.J., 1994. A finite-volume method for fluid flow simulations with moving boundaries. Computational Fluid Dynamics 2, 19–40.
892
Y.G. Lai et al. / Journal of Biomechanics 35 (2002) 881–892
Lee, C.S., Chandran, K.B., 1995. Numerical simulation of instantaneous backflow through central clearance of bi-leaflet mechanical heart valve at closure: shear stresses and pressure fields. Medical and Biological Engineering and Computing 33 (3), 257–263. Lee, C.S., Chandran, K.B., Chen, L.D., 1994. Cavitation dynamics of mechanical heart valve prostheses. Artificial Organs 18, 758–767. Makhijani, V.B., Yang, H.Q., Singhal, A.K., Hwang, N.H.C., 1994. Computational analysis of MHV cavitation: effects of leaflet squeezing and rebound. Journal of Heart Valve Disease 3 (Suppl. 1), S35–S48. Makhijani, V.B., Yang, H.Q., Dionne, P.J., Thubrikar, M.J., 1997. Three-dimensional coupled fluid-structure simulation of pericardial
bioprosthetic aortic valve functioning. ASAIO Journal 43 (5), M387–M392. Reif, T.H., 1991. A numerical analysis of backflow between the leaflets of St Jude medical cardiac valve prostheses. Journal of Biomechanics 24 (8), 733–741. Wu, Z.J., Wang, Y., Hwang, N.H.C., 1994. Occluder closing behavior: a key factor in mechanical heart valve cavitation. Journal of Heart Valve Disease 3 (Suppl. I), S25–S34. Zapanta, C.M., Liszka Jr., E.G., Lamson, T.C., Stinebring, D.R., Deutsch, S., Geselowiz, D.B., Tarbell, J.B., 1994. A method for real-time in vitro observation of cavitation on prosthetic heart valves. ASME Journal of Biomechanical Engineering 116, 460–468.