Microsoft Word - 48dallavecchia.docx CHEMICAL ENGINEERING TRANSACTIONS VOL. 49, 2016 A publication of The Italian Association of Chemical Engineering Online at www.aidic.it/cet Guest Editors: Enrico Bardone, Marco Bravi, Tajalli Keshavarz Copyright © 2016, AIDIC Servizi S.r.l., ISBN 978-88-95608-40-2; ISSN 2283-9216 Hemodynamic Design Optimization of a Ventricular Cannula Applying Computational Fluid Dynamics (CFD) Alifer D. Bordonesa, Luis R. Rojas-Solórzano*b aDepartment of Biomedical Engineering. University of Texas at San Antonio, San Antonio, TX, USA. bSchool of Engineering. Nazarbayev University, Astana, Rep. of Kazakhstan. luis.rojas@nu.edu.kz Nowadays, a Ventricular Assist Device (VAD) is a hope of clinic solution to patients with heart failure. This study presents the design of the “trumpet mouth” ventricular cannula (TMVC) connected to a pediatric VAD, optimized to reduce the conditions that generate blood damage. The flow through the TMVC is analysed using Computational Fluid Dynamics (CFD) techniques. The optimization of the TMVC considered the turbulent and byphasic nature of the flow at the systolic and diastolic phases of the cardiac cycle. To assess the hemodynamic performance of the cannula, three geometric parameters were chosen as input or independent variables in the optimization: the cannula tip curvature angle (θ) and the inner and outer curvature radii (R1 and R2). The output parameters or objective functions were the Modified Index of Hemolysis (MIH), the Modified Index of Platelet Lysis (MILpl) and the fluid pressure drop (ΔP). The mathematical optimization passed through an enhanced sampling process of combined scenarios of input parameters and thereafter, the response surface was generated using the Kriging algorithm. The selection of the geometric parameters that minimized MIH, MILpl and ΔP were obtained after the application of a Multi-Objective Genetic Algorithm (MOGA) and the generation of Pareto fronts to obtain the points that met the established conditions of reduction of blood damage. The optimal geometric parameters for the ventricular cannula obtained through the optimization process are: θ (27.29 degrees), R1 (0.83 mm) and R2 (0.20 mm). The final results showed a reduction in MIH for both cardiac phase configurations in the zone of the cannula tip and an overall reduction of 25% and 60% in the MIH and MILpl, respectively, compared to the damage for the cannula with nominal input parameters. 1. Introduction When blood is in contact with a non-biological environment, many complications take place. Thermal and osmotic effects and even the contact with the walls of the container can cause damage to the blood cells. However, the most relevant type of damage is caused by the abnormal exposition of the cells to shear stresses (Salazar et al., 2008). Under certain conditions, these shear stresses can cause the rupture or partial damage of the membrane covering the red blood cells allowing the release of hemoglobin into the plasma in a process known as hemolysis. Shear stresses do not only affect the red blood cells; indeed, it has been proven that under their effect platelets can become activated generating changes in the fibrinogen receptors and increasing the possibility to generate some aggregates to them, making them more likely to stick to the surface of the blood vessels and get activated creating thrombus. Therefore, consequences of the platelet activation can be catastrophic for the patient (Tandon et al., 1997). It is of vital importance to take into consideration the potential damage to the blood cells when designing medical devices. Finding the best design involves the construction of multiple prototypes to measure the dependent variable (objective) in each one of them and determine which one depicts the best performance. To accomplish this process experimentally, a considerable amount of tests must be carried on which degenerates in a very expensive, controversial and time consuming task to obtain results that may still not consider all the spectrum of possibilities before choosing the final design. Mathematical optimization is a tool often used in the design of physical components and devices and can be defined as the technique to determine the “best” solution to certain problem defined under a mathematical formulation subject to constraints. These problems are usually DOI: 10.3303/CET1649101 Please cite this article as: Bordones Gonzalez A., Rojas-Solorzano L., 2016, Hemodynamic design optimization of a ventricular cannula applying computational fluid dynamics (cfd), Chemical Engineering Transactions, 49, 601-606 DOI: 10.3303/CET1649101 601 representations of real physical models (Fletcher, 1987). Computational Fluid Dynamics (CFD) has been used for the last several decades to simulate and analyse the flow of fluid in many diverse applications. The objective of this research is to optimize via CFD the design of the cannula of a pediatric Ventricular Assist Device (VAD) in order to reduce the damage caused to the blood cells as a secondary effect associated to the flow of blood through this device. 2. Methods 2.1 Modelling and Meshing The ventricular cannula of a VAD has the function to guide the flow of blood from the left ventricle into the pump of the device. The ventricular cannula analysed in this investigation is based on the model proposed by Antaki et al., (1995) called the “trumpet mouth” ventricular cannula (TMVC). In his investigation, Antaki et al. showed that the TMVC presented a superior hemodynamic performance than its counterparts. The 3- dimensional model, a sketch of the internal contour and internal dimensions of the ventricular cannula are presented in Fig. 1. The geometric parameters φin and φex represent the internal and external diameters, respectively, H represents the height of the cannula, θ refers to the angle of the cannula flare measured from a perpendicular line to the extension pipe connecting to the pump, the radius R1 represents the curvature from the cannula flare into the pipe and R2, the curvature of the outer surface of the device in contact with the left ventricle. Figure 1: Ventricular cannula model, sketch of the internal contour and geometric dimensions. To optimize the ventricular cannula it was necessary to consider not only the device, but also the volume of the left ventricle since the device is directly attached to it and the interaction between them will affect the fluid dynamics of the circulating blood. Since the cannula under study is for pediatric purposes, it was necessary to calculate the volume of the ventricle appropriately. Graham et al. (1971) presented a series of correlations that permitted calculating the ventricle volume based on physical factors such as infant height, weight and body surface area, among others. The equation considered for the calculation of the ventricular volume in diastolic configuration was the corresponding to a two-year old child. To calculate the left ventricle volume in systolic configuration, the ejection fraction (EF) value was used. The EF is a measurement of how much blood is expelled by the left ventricle at each contraction. For children of age equal or higher to 2 years, the EF is approximately 0.65 or 65% (Graham et al., 1971). The weight and height considered was obtained from data published by the U.S. Department of Health and Human Services (2008). In addition to the ventricle volume for both pumping phases, the contact angle of union between the wall of the left ventricle and the cannula tip boundary was considered. Curtis et al. (1998) showed in their experiments that the angles of contact between the cannula wall and the TMVC are 85 degrees for the diastolic phase and 54 degrees for the systolic phase. The computational model of the cannula-ventricle system was generated using CAD tools followed by mesh generation, all embedded in ANSYS Workbench Platform v13. Tetrahedral elements were used with an adequate densification of elements in the zones of interest as, for example, the contraction region or the zones close to the walls where an inflated layer of smaller prismatic elements was applied to capture flow gradients accurately. For every configuration analysed in this investigation, it was verified that numerical results were independent of the mesh discretization. The verification method used was proposed by Roache (1994) and is called the Grid Convergence Index (GCI). The comparison of the values reported by every mesh was made in four (4) lines of study located in the cannula tip area. The CPU-time required to run each 3D simulation using a 2.5GHz Pentium Dual Core with 4 GB of memory RAM was around 24 hours, suggesting a further simplification to make the optimization process viable in a reasonable time framework using modest computational resources. Hence, the dimensionality of the model was reduced by taking advantage of its 602 axisymmetric condition, and a final lighter model was obtained with a reduction of CPU-time by 70%. The GCI study showed that the use of the medium mesh met the criterion established for the error (maximum 5% of difference between the results obtained with two subsequent density meshes) and it was selected as the final mesh. Results for 2D model were compared with 3D preliminary ones and it was proven that the difference was less than 5% for every point in each velocity profile obtained with the fully three-dimensional model. 2.2 Numerical Modelling and Optimization ANSYS Design Exploration v13.0 was used as the platform to optimize the design of the ventricular cannula. For this study, in every case, a second-order discretization scheme was applied for all the equation terms and the turbulence was modelled applying the Shear Stress Transport model (SST) with gamma-theta transition that was proven to be appropriate to reproduce the blood flow in cannula-type devices (Salazar et al., 2008). A constant blood flow of 1.5 L/min was prescribed as inlet boundary condition for the cannula-ventricle system. The blood was modelled as a multiphase fluid with red blood cells and platelets treated as particles of a disperse phase interpenetrating the plasma, considered as a continuous phase. The properties assigned to each phase are presented in Table 1. The remaining boundary conditions include non-slip condition at the wall and static pressure and derivative of velocity in the normal direction equal to zero at the outlet of the system. Table 1: Material properties for the red blood cells (RBCs), platelets and plasma. Phase φparticle (mm) ρ (kg.m-3) μ (cP) Volumetric fraction Plasma Continuous 1025 1.22 58 RBCs Disperse 5.64 1100 2.00 40 Platelets Disperse 2.00 1040 3.00 2 The blood damage was calculated by applying the multiphase CFD models developed by Mubita (2012) and Moreno (2012) to predict Modified Index of Hemolysis (MIH) and Modified Index of Platelet Lysis (MILpl), respectively. The model by Mubita (2012) estimates the hemolytic damage caused to red blood cells based on the formulations proposed by Giersiepen et al. (1990) and Garon and Farinas (2004) where the damage to the membrane of the red blood cells is presented in Eq. (1) as a function of the shear stress (τ) and the exposure time (t) of the cell to it. HGW is the global concentration of free hemoglobin in the plasma. 785.0416.271062.3 t Hb Hb H GW τ −×= Δ = (1) Similarly, platelet lysis depends also on shear stress (τ) and exposure time (t) to it over the platelets membrane. The methodology proposed by Moreno (2012) to calculate the platelet lysis is based on the model by Goubergrits et al. (2006) shown in Eq. (2). Lpl is the index of platelet lysis and measures the variation of the concentration of LDH (Lactate Dehydrogenase) in the plasma. ( ) t LDH LDH L pl 77.0/075.377.0/18 77.0/1 1066.3 τ−×=      Δ= (2) The optimization process was divided in three main steps: design of experiments, construction of the response surface and application of the optimization algorithm. The variables chosen as input parameters were the curvature angle of the cannula tip (θ), the curvature radius of the cannula pipe (R1) and the external border (R2). The internal and external diameters were not considered because these are values that would not only affect the hemodynamics but also the process of insertion to the ventricle, positioning and connection to the pump, going beyond the scope of this research. As for output parameters or objective functions, the MIH, MILpl and the pressure drop through the cannula were chosen. The pressure drop even when is not an indicator of damage, is an important factor because if the pressure differential is too high, the blood flow will decrease affecting the patient’s health. The range of variation for each input parameter is θ (15 to 55 degrees), R1 (0.8 to 1.2 mm) and R2 (0.2 to 0.3 mm). The design of experiments, first step in the optimization process, is a technique that determines the location of sample points so that the input space is explored in the most efficient way allowing to collect all the information required with a minimum number of points. The efficiency of the location of these points increases the precision of the response surface generated from them. The design of experiments scheme applied in this work was the Optimal Space Filling Design with Central Composite Design equivalent sampling units. To obtain the response surface, the Kriging algorithm was used and the surface goodness of fit was controlled to verify that the surface obtained predicted accurately the behavior of the outlet parameters with respect to the 603 inlet parameters. For this investigation the final optimization process was performed applying two sequential techniques. Firstly, the Screening algorithm, using a sample generator based in a variation of the Hammersley method, and secondly, the MOGA algorithm (Multi-Objective Genetic Algorithm) which is based in the concept of controlled elitism to determine the solution from a set of possible solutions or Pareto optimal set (Cavazzutti, 2013). 3. Results 3.1 Blood damage: hemolysis and platelet lysis Performed simulations showed that the concentration of free hemoglobin in the plasma increases in the curved region of the cannula and the region of contact between the cannula tip and the ventricle wall. These higher concentration values were expected according to Eq. (1), since HGW depends on the variation of the shear stress which is closely related to the velocity gradient in the radial direction. In the curved region there is an increment of the speed and a reordering of the flow caused by the sudden change in the area at the entrance of the cannula; the particles close to the wall reduce their velocity to adjust it to the velocity of the plenum. Likewise for hemolysis, platelet lysis was calculated applying a power law depending of two factors: shear stresses and the exposure time (Eq. (2)). High shear stresses were located in the zones close to the wall where the velocity gradients were larger, increasing the Index of Platelet lysis in those regions. The platelet lysis, as for hemolysis, increased when the velocity of the fluid increased relative to the radial distance, causing a rise of the shear stresses. The high speed zone is located in the entrance of the cannula pipe where the fluid needs to accelerate. In systolic and diastolic configuration, the blood damage was calculated in the zone of the cannula tip. Results are presented in Table 2. Table 2: Blood damage calculated in the cannula tip for the cannula-ventricle system in systolic and diastolic configuration. Blood damage MIH MILpl Systole 1.28 0.23 Diastole 15.1 2.83 When comparing the results for diastolic and systolic conditions, there are several important differences. First, the area affected by high shear stresses was smaller for the former as it corresponds to lower velocity gradients in the radial direction. This responds to variations in the velocity field which are present as a consequence of the different angles between the cannula tip and the ventricle wall. Blood damage is considerably less in systolic condition; in fact, 92% lower MIH and MILpl is reported in this condition compared to diastolic, for the nominal geometry. This is explained by the configuration of the joined section between the ventricle wall and the cannula tip; i.e., in systole, the angle (54 degrees) allowed the fluid to move towards the cannula in a harmonious way and with smaller acceleration rates and, therefore smaller shear stresses, than found in diastole (85 degrees angle) for the nominal configuration (i.e., before the optimization algorithm was carried on). 3.2 Optimization As explained in the methodology section, the optimization process starts with the sampling of scenarios associated to the design of experiments. The design space was filled up with 15 design points and the response surface was built for systole and diastole using the Kriging algorithm. Before the final refined response surface was developed, a “goodness of fit” quality check was performed. The response surfaces obtained for the cannula-ventricle system in diastole configuration are presented in Fig. 2. The response surfaces show that the value of MIH decreases when the values for the angle θ are very low or very high, but in specific combination with the radius R1 and R2. In the case of R1 and θ, MIH is smaller when θ is between 15° – 20° and R1 between 0,8 – 0,9 mm. This is explained by the diastolic geometric configuration with a small angle, mimicking a pipe with a small contraction and thus, promoting a velocity gradient at the entrance of the pipe where the flow accelerates. In the case of R2 and θ, the value of MIH decreases when θ is between 15° – 20°, and again close to 55°, for values of R2 very low, between 0,20 and 0,21 mm. Just as in the previous case, a pronounced angle with a small radius between the cannula tip and the ventricle wall generates a geometry similar to a pipe with a contraction (in this case playing a pipe extension) and this generates less intense shear stresses, and in consequence, smaller MIH. For θ between 15° - 20° the behaviour is similar. The MILpl resulted small for θ between 25°- 45°, R1 between 0,8 - 1,2 mm and R2 between 0,23 - 0,3 mm. 604 Figure 2: Response surfaces obtained for the cannula-ventricle system in diastole configuration. Variation of MIH (a. and b.) and variation of MILpl (c. and d.) with changes in the radius R1, R2 and θ. The levels of MIH found in systole configuration are very different from diastole. MIH increases when θ is low, and also for 15°, corresponding to an almost flat cannula, which in contact with the ventricle wall in an angle of 54° generates a slow reduction of the diameter, causing the fluid flow to undergo through a sudden change of area between the cannula tip and the pipe, producing larger changes in the velocity and the increment of shear stresses. The acceptable values of θ for any R1 and R2 are between 30° and 40° because this angle makes the transition smooth between the ventricle wall and the cannula tip, reducing the likelihood of an increment of the shear stresses and MIH. The MILpl is low for angles between 35° and 50°, for values of R1 larger than 0,85 mm and R2 between 0,2 - 0,25 mm. This indicates a difference regarding MIH, where the value of the radius, especially R2 influences the MILpl. The values of R2 between 0,2 - 0,25 mm allow a smaller discontinuity between the cannula tip and the ventricle wall, reducing the likelihood of generating a zone of high concentration of shear stresses which would increase the platelet lysis. This also reduces the zones of stagnation and recirculation of the fluid, a very important factor in the platelet activation. Since the simultaneous minimization of the MIH, MILpl and pressure drop were the two main objectives of the optimization, weighted objectives were used to obtain a unique linear combination objective function with larger weights on blood MIH and MILpl. The optimization algorithm led to six (6) possible candidates for the final design. The values of MIH, MILpl and the pressure drop in all of them are acceptable and all of them reduced the blood damage calculated in the original geometry. The chosen candidate was the best candidate for the system in diastolic configuration because it is the cardiac cycle phase with the greatest amount of blood damage; approximately 92% more than the values estimated in systole. Figure 3 shows two views of the final design of the ventricular cannula and the optimal values for the geometric (input) parameters obtained through the optimization process for the system. Table 3 shows the blood damage comparison between the final design and the initial cannula parameters. Figure 3: Final optimized design of the TMVC obtained by applying CFD tools and mathematical optimization techniques. Table 3: Comparison of the blood damage calculated for both initial and optimized ventricular cannula design. Ventricular Cannula Diastole Systole Optimized design Initial design Optimized design Initial design MIH 12.2 15.1 0.11 1.28 MILpl 0.97 2.83 0.28 0.23 ΔP [Pa] 284 287 150 154.9 605 The final results show a decrease of the values of MIH for both cardiac cycle configurations in the zone of the cannula tip. The reduction is in the order of 18% for diastole and 92% for systole. Nevertheless, for the MILpl the reduction of damage is 65% in diastole, but it increases 17% in systole. Even with the increase of the value of the Modified Index of Platelet Lysis in systole, the overall decrease of damage is 60% compared to initial design. The total reduction of the MIH is 25%. 4. Conclusions This investigation presents the design of a “trumpet mouth” ventricular cannula (TMVC) optimized for a better hemodynamic performance. The optimal design was obtained using mathematical multi-objective optimization techniques and the application of CFD techniques to estimate blood damage through medical devices. The design candidates obtained in the Screening and MOGA algorithms presented the best set of input parameters that reduce the damage caused by hemolysis and platelet lysis by 25 and 60 percent, respectively, compared with the initial design. Although the angle obtained is not the best option when the system is in systolic configuration, it was considered that the ventricle expansion phase (diastolic) is the most critical, because is in this phase of the cardiac cycle when the greater damage to the blood is produced (92% more damage compared with the systolic phase). Acknowledgments The authors thank the Laboratory of Fluid Mechanics at the Simon Bolivar University for all the collaboration given in the development of this project. References Antaki J., Dennis T., Konishi H., Brown M., Tomczak J., Kerrigan J. Kormos, R., 1995, An Improved Left Ventricular Cannula for Chronic Dynamic Blood Pump Support. Artificial Organs, 19(7), 671-675. Cavazzutti, M., Optimization Methods – From Theory to Scientific Design and Technological Aspects in Mechanics. Springer, Modena, Italy. Curtis A., Wu Z., Kormos R., Griffith B., Antaki J., 1998, Novel Ventricular Apical Cannula: in vitro evaluation using transparent, compliant ventricular casts. ASAIO, 44(5), 691-695. Fletcher R., 1987, Practical Methods of Optimization. John Wiley & Sons LTd, West Sussex, England. Garon A., Farinas M., 2004, Fast Three-dimensional Numerical Hemolysis Approximation. Artificial Organs, 28 (11), 1016-1025. Giersiepen M., Wurzinger L.J., Optiz R., Reul H., 1990, Estimation of shear stress-related blood damage in heart valve prostheses-in vitro comparison of 25 aortic valves. International Journal of Artificial Organs, 13(5), 300-306. Goubergrits, L., 2006, Numerical Modeling of Blood Damage: Current Status, Challenges and Future Prospects. Expert Review of Medical Devices, 3(5). pp. 527-531. Graham T., Jarmakani J., Canent R., Morrow M., 1971, Left Heart Volume Estimation in Infancy and Childhood: Reevaluation of Methodology and Normal Values. Circulation, 43, 895-904. Moreno J., 2012, Estimación numérica de lisis plaquetaria en dispositivos de pequeño diámetro empleando modelos multifásico CFD-Eulerianos. Universidad Simón Bolívar, Sartenejas, Venezuela. Mubita T., 2012, Efecto de la Segregación de las células sanguíneas en la estimación numérica de Hemólisis. Universidad Simón Bolívar, Sartenejas, Venezuela. Roache P., 1994, Perspective: A method for uniform reporting of grid refinement studies. J. of Fluids Engineering, 116, 405-413. Salazar F., Rojas-Solórzano L., Blanco A., 2008, Turbulence modeling in the numerical estimation of hemolysis in hemodialisis cannulae. Revista de la Facultad de Ingenieria U.C.V, 23(4), 93-98. Tandon P., Diamond S., 1997, Hydrodynamic effects and receptor interactions of platelets and their aggregates in linear shear flow. Biophysical Journal, 73(5), 2819-2835. U.S Department of Health and Human Services - Center for Disease Control and Prevention, 2008, Anthropometric Reference Data for Children and Adults: United States, 2003–2006 < http://www.cdc.gov/nchs/data/nhsr/nhsr010.pdf > accessed 08.11.2012 606