- Email: [email protected]
A CFD model for the prediction of haemolysis in micro axial left ventricular assist devices
Applied Mathematical Modelling 37 (2013) 4199–4207
Contents lists available at SciVerse ScienceDirect
Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm
A CFD model for the prediction of haemolysis in micro axial left ventricular assist devices D. Carswell a,⇑, D. McBride a, T.N. Croft a, A.K. Slone a, M. Cross a, G. Foster b a b
Civil and Computational Engineering Centre, College of Engineering, Swansea University, Swansea, West Glamorgan, Wales SA2 8PP, UK Calon Cardio-Technology Ltd, Institute of Life Science, Swansea University, Swansea, West Glamorgan, Wales SA2 8PP, UK
a r t i c l e
i n f o
Article history: Received 4 October 2011 Received in revised form 13 July 2012 Accepted 10 September 2012 Available online 23 September 2012 Keywords: Haemolysis CFD Blood pump Axial flow
a b s t r a c t This paper describes the development and application of a computational model based upon Computational Fluid Dynamics (CFD) software simulation technology to predict haemolysis in micro Left Ventricular Assist Devices (lLVAD). A CFD model, capturing the full three dimensional geometry of the device, together with an explicit representation of the rotating machinery based upon a rotating reference frame, is solved transiently. Mixed meshes with the order of a million elements are required to resolve the flow adequately and so to enable solutions in a reasonable time (e.g., 3 h) the model is solved on a high performance parallel cluster. Haemolysis is a measure of damage occurring in the blood and is conceived as accumulating as it passes through parts of the device where it encounters high shear forces. As such, the haemolysis model is based upon tracking the behaviour of particles released at the inlet throughout the flow domain and calculating the damage accumulated by each individual particle as it traverses the device. In order to ensure the model predictions of haemolysis are noise free from a statistically significant perspective then it is demonstrated that the number of particles to be tracked must exceed 20 000 in any simulation experiment. Comparisons with experimental data from a companion paper demonstrate the effectiveness of the CFD simulation embedding the haemolysis model. Ó 2012 Elsevier Inc. All rights reserved.
1. Introduction Human, and mammal, red blood cells (RBCs) exist to transport oxygen to every living cell in the body. RBC structure comprises an elastic membrane that encases a volume of haemoglobin in a bi-concave disc of diameter to the order of 8 lm [1]. The membrane comprises a lipid external layer, internal proteins and a cytoskeleton [2]. The cells may deform in different ways under different conditions, usually dependant on shear stress. In healthy human veins and arteries, the shear stresses experienced may range from 0.1 to 20 Pa, increasing in the smaller capillaries [3]. In high shear stress regions the RBC membrane cannot sustain its deformed shape for long and eventually loses its elasticity resulting in mechanical damage. As membrane damage occurs, haemoglobin may escape from the RBC and into the blood plasma. Haemolysis expresses the amount of haemoglobin that escapes in this manner [4]. A relationship linking shear stress (s) and exposure time (t) to plasma-free haemoglobin (pfHb) first proposed by Blackshear et al. [5] is
pfHb ¼ ts2
⇑ Corresponding author. E-mail address: [email protected] (D. Carswell). 0307-904X/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.apm.2012.09.020
ð1Þ
4200
D. Carswell et al. / Applied Mathematical Modelling 37 (2013) 4199–4207
but only tested in the range 1 < s < 10Pa for 1 h. However, the same relationship was found to hold true for s < 8000Pa and t < 102 s [6]. Hb A more general form expressing the proportion of haemoglobin released DHb is related to s and t by
DHb ¼ Ct a sb : Hb
ð2Þ
The values of the constants a; b and C proposed by Giersiepen et al. [7] were 0.785, 2.416 and 3:62 107 respectively. A slightly different set of parameters were cited by Song et al. [8], these being 0.765, 1.991 and 1:8 106 in their CFD simulations. More recently, a different set of parameters (0.747, 2.004 and 1:21 107 respectively) derived for ovine blood 0 1=a were published by Taskin et al. in their CFD simulation [9]. Taskin proposes using a convected variable Hb ¼ Hb to accumulate haemolysis through a pump such that 0
0
dðDHb Þ @ DHb þ qu ¼ S; dt @x
ð3Þ
where the source term S ¼ qðHbC sb Þ1=a . As boundary conditions, it is assumed Hb ¼ 0 at the inlet and zero flux is specified at all solid walls. Using a convected variable such as this it is possible to quickly identify areas of a pump that are more haemolysing than others - something that is not always possible with simply following pathlines. There still appears to be some debate on the upper bound of the safe shear stress, i.e., the level below which haemolysis does not occur or is insignificant. Previous research suggests this value is within the range of 150–800 Pa [10–13]. Computational Fluid Dynamics (CFD) has long be used to aid designers obtain a greater understanding of, not just heart pump operation, but many other applications. The detailed insight it gives may be difficult or nearly impossible to obtain experimentally. Many researchers have already used CFD as a means of analysing both axial and radial heart pumps. Watanabe et al. [14] analysed a centrifugal heart pump using TASCflow (AEA, Pittsburg, USA) utilising fully structured meshes and comparisons were made with experimental data to establish the correctness of the model. Throckmorton et al. [15] used Ansys CFX to evaluate a pump design prior to manufacture to make sure that it would meet the requirements. They then went onto analyse a candidate pump in more detail [16] to ensure the magnetic bearings would be able to support the force exerted on the rotor due to hydrodynamic forces. Another axial pump was also analysed experimentally and numerically, again using TASCFlow, by Chan et al. [17]. Fluent 6.2 was used by Chua et al. [18] in their computational model making use of sliding meshes but has not been compared to experimental data to the best of the authors’ knowledge. An axial pump containing both an inducer and straightener was analysed by Li and Wang [19] using the RNG k- turbulence model and rotating frames of reference. They concluded that inducers can have quite a large effect on pump performance. The NUMECA commercial CFD software was used to perform computations on the axial pump of Zhang et al. [20] but not compared to experimental data. It is clear from the above literature that compare numerical and experimental results that, given an adequate mesh, commercial CFD-based simulation software are capable of matching experimental data. Using pathlines is the general technique used to capture data required for haemolysis and has been employed by a variety of other researchers [21,22,8] as it is, so far, the only way to generate appropriate data for haemolysis estimation. Essentially, the CFD solver produces a flow field through the blood pump and then mass-less particles are tracked through this flow field to assess the damage they experience as a consequence of the shear stress exerted by the characteristics of the flow at any and every location visited by a specific particle. Of course, the data from a large number of particle tracks has to be integrated to provide statistically significant information. To compare haemolysis predictions with experimental data it is often necessary to convert both datasets to the Normalised Index of Haemolysis (NIH). Experimentally, NIH is determined by
NIH ¼ H V ð1 Hct Þ
100 ; Q T
ð4Þ
where V is the volume of blood in the test loop (L), Hct is the haematocrit, Q is the flow rate (L min1) and T time (min). Values of T and V are not available numerically so a slightly different formulation is used
NIH ¼ 100
DHb ð1 Hct Þ j; Hb
ð5Þ
where j is the haemoglobin content of blood. Typical values of Hct and j are 45% and 150 g L1 [23]. This paper describes the computational models employed by the authors to characterise the flow behaviour of lLVADs and to produce predictions of expected performance curves and haemolysis levels. The outcomes of the models are currently underwriting the work of designers in advancing future pump designs by allowing an in-depth analysis of the blood flow behaviour in current pump geometries. 2. Pump design The lLVADs currently being designed are destined as short term devices. That is, they are designed to produce flow rates of around 2.5 L min1 to work alongside a recovering heart that is still functional. Other pumps are designed on the bridgeto-recovery principle or as more permanent devices in which case flow rates beyond 5 L min1 are required.
D. Carswell et al. / Applied Mathematical Modelling 37 (2013) 4199–4207
4201
Fig. 1. Pump geometry considered in this paper.
In this paper, the lLVAD shown in Fig. 1 is considered measuring 7.5 mm in diameter by 40 mm total length. Three blades are 1.5 mm in height with a clearance of just 20 lm from the case. The main inlet is located immediately upstream of the rotor. A three-vaned rotor serves to provide the pumping force and the stator converts whirl energy to a pressure head. Having the main inlet in the side of the pump means the main flow path through the pump is as short as possible. The motor and hydrodynamic bearings are located at the rear of the pump where a lubrication channel is also present. Using blood as the lubricating agent eliminates the need to use other lubricants which also eliminates the need to keep them separate from the blood, hence the design is simplified. The pump is designed to operate between 15,000 and 30,000 r min1. Pump performance is evaluated at a variety of operating pressures and rotation speeds, often outside the anticipated working range to assess the full range of physiological conditions. In laboratory tests, the pump was connected via lengths of pipe to a reservoir, the outlet tubing was restricted via a clamp to provide a pressure at the outlet. Flow rates and inlet and outlet pressures were recorded experimentally. Experiments were conducted using water. 3. Computational models and methods The 3D incompressible Navier–Stokes equations are used to predict the flow field through the pump. In line with other researchers [14–20] we assume the blood flow to be incompressible. Rotating machinery requires special consideration in the formulation of the equations. All rotating parts are enclosed in a rotating domain modelled using a rotating reference frame. In the rotating frame, a relative velocity ur is related to the absolute velocity u, position vector r and angular velocity vector x by ur ¼ u þ x r. Appropriate substitutions into the convection term of the Navier–Stokes equations must be made in the rotating frame by either solving the relative or absolute velocities. In this application the latter was chosen yielding a new formulation of the governing equations
q
@u þ qur ru þ qðx uÞ ¼ rp þ rðlruÞ þ Su : @t
ð6Þ
Mass is also conserved using the continuity equation (Eq. (7)) but requires no special treatment.
qru ¼ Sm ;
ð7Þ
Using the rotor diameter as a characteristic length (‘), Reynolds numbers,
Re ¼
qv ‘ l
ð8Þ
were computed to be never greater than 5000 and for most of the time and space considerably less, so that while the global flow is just outside the laminar flow regime it is by no means fully turbulent. Blood density was kept constant at 1050 kg m3 and the Quemada viscosity model was used to model blood viscosity:
l¼
l
p 1 k H Þ2 2 Q ct
;
ð9Þ
pffiffiffiffiffi k0 þ k1 cr c pffiffiffiffiffi and cr ¼ : cc 1 þ cr
ð10Þ
ð1
where
kQ ¼
Blood plasma viscosity, lp ¼ 1:23 m Pas and the level of haematocrit, Hct ¼ 48%, for a healthy person. Shear strain rate c is obtained from the flow field and k0 and k1 are model constants that represent the rate of RBC rouleaux formation at zero and infinite shear. When it came to considerations of turbulence, a combination of the following reasons: Most of the flow is at worst just bordering on the turbulent regime and much of it is truly laminar, The viscosity is actually non-Newtonian and so someway outside the scope of the original conception of conventional computational turbulence models, and
4202
D. Carswell et al. / Applied Mathematical Modelling 37 (2013) 4199–4207
The very fine meshes and small time steps enables capture of many of the larger more significant vortices led to the conclusion that there was no real justification for employing a formal turbulence model since all were inadequate in some essential manner. Moreover, since the fine meshes and small time steps captured much of key vortex behaviour it is argued that much of the developing turbulent behaviour is represented in this computational implementation. The equations are solved using the in-house code PHYSICA [24] using the vertex-based finite volume discretisation method [25]. In this scheme flow variables are stored at nodes and the finite element shape functions are used in the computation of face fluxes. For any variable / the face gradients are calculated using the shape function gradients such that
@/ @xi
f
n X @Nj j¼1
@xi
/j ;
ð11Þ
where N j is the shape function associated with node j. Face values /f are computed from the nodal values and there are a number of high-order schemes able to do this. In this work a second-order Mass-Weighted (MAW) upwind scheme [26] was implemented for unstructured 3D meshes based on the work of Tran et al. [27] Solving the Navier–Stokes equations requires the use of a pressure–velocity coupling. For co-located variables, the method of Prakash and Patankar [28] was implemented. As the flow rate generated by the rotor is sought, only pressure boundary conditions could be prescribed at both the inlets and outlet. Both inlets were give 0 Pa pressure and the outlet pressure was varied between 3300 Pa and 200 000 Pa to simulate different physiological conditions. A velocity is present on the rotor, in the inertial reference frame, by virtue of the rotating frame. The outer case and stator were specified as no-slip wall conditions. An important part of the overall computational procedure was the generation of meshes of suitable quality for the CFD model where part of the mesh is rotating. The meshing strategy employed a mix of element types. The nature of the device presented a number of specific meshing challenges. In regions of small gaps it was found that it was easier to control the mesh quality and density by using structured, hexahedral elements rather than tetrahedra. Ansys ICEM CFD 13 (Ansys Inc.) was employed as the mesh generation software. Creation of hexahedral elements made use of ICEM’s blocking capabilities. Mesh element quality was kept above 0.4. The cone of the rotor was separated from the stator by only 0.1 mm. Knowing that this gap had to accommodate an interface between the rotating and stationary domains meant that the mesh in this region was, perhaps, finer than it needed to be. The main flow path was meshed using tetrahedra of size 0.1 mm except in the gap between rotor and case in which case the element size was reduced to provide sufficient resolution of the flow in this gap. A range of meshes were built from coarse (500 000 elements, 0.5 mm element size) to fine (2 000 000 elements, 0.1 mm element size) to assess any mesh dependance issues. A coarse mesh is shown in Fig. 2 where the darker region, to the right of the figure, indicates the elements belonging to the rotating reference frame. The figure has been chopped to show the region of interest as the computational inlets and outlets are positioned further away from the actual inlets and outlets so as to not cut short any flow features that may be present. Results showed a semi-fine mesh (1 000 000 elements, 0.2 mm element size) was adequate for the simulation where the predicted flow rate changed by <1%. The solution was deemed converged when residuals of mass and the three velocity components had fallen by three orders of magnitude. The host CFD code is fully parallelised using a distributed memory model [29] and well adapted to run in parallel for various rotating machinery representations [30]. Computations were run in parallel using MPI as the interface between the different processes using two quad-core 64-bit Intel Core i7 chips. A transient analysis was run over a period of 30 time steps with the interval being 0.005 s which was enough for the simulation to produce a steady flow rate at the outlet. One simulation took approximately 3 h to complete. 3.1. Pathlines Once a converged flow field had been computed, pathlines of mass-less particles were tracked through it. A number of pathlines were released at random points within each element at the inlet. The number of particles released per element
Fig. 2. Coarse mesh of a typical pump. The darker elements show the region of the mesh in the rotating reference frame.
D. Carswell et al. / Applied Mathematical Modelling 37 (2013) 4199–4207
4203
is known as the clustering factor and ranges from 1 to 35 resulting between 900 and 30,000 particles. Each step along pathline i is computed from the previous step cni and the current velocity field u by
cnþ1 ¼ cni þ Dt i
@ui : @xi
ð12Þ
A constant time step Dt ¼ 1 105 s was chosen to ensure accurate pathlines. Element sizes and velocities showed the maximum time step permissible was 5 105 s meaning the time step used was suitable. A greater timestep did not provide sufficient accuracy due to the velocities and cell sizes involved. Smaller time steps did not produce any noticeable difference in the computed paths except that it consumed a lot more processing power. A particle ends when it reaches the outlet. Another undesirable condition that causes particles to end is when they hit a wall. One method would be to reflect the particle using the boundary face normal. Occasionally, particles may become trapped if they encounter recirculation regions, therefore never exiting the pump. Only particles that reached the outlet are considered in the haemolysis computations. 3.2. Haemolysis While computing the pathlines, the shear stress is also output at each time step as it is one of the factors influencing haemolysis. The time step at each point along each pathline was used as the exposure time (dt) and the shear stress (s) was assumed constant throughout this period. A number of inconsistencies were found by relying on Eq. (2) to compute haemolysis [31]. In particular it does not account for load history, that is, two particles reaching the same point in the flow via different pathlines may accumulate different levels of damage in accordance with damage previously experienced. This is due to the relaxation time a stresses RBC membrane has once a stress is removed or reduced. Grigioni et al. [31] proposed
" #a1 n i X X b=a Ca sðtj Þ dt þ Dðt0 Þ sðti Þb=a dt i¼1
ð13Þ
j¼1
as a more discriminating means of computing blood damage. The inner sum accounts for the load history where n is the number of steps per particle and the function D allows an initial damage to be specified. Note that s is a function that varies with position along the pathline. In this application, all particles are assumed to have zero initial damage, therefore, Dðt 0 Þ ¼ 0. Damage is then converted to NIH using Eq. (5). 4. Results and discussion Pump performance characteristics are presented in Fig. 3. At the higher rotation rates, the pressure generated by the rotor blades becomes much greater than the outlet pressure and, hence, the curves begin to converge. The predicted flow rates for 0 Pa begin to diverge from the experimental results as flow rate increases. The reason for this is thought to be extra losses caused by the length of pipe used experimentally which was not modelled computationally. It should be noted that for comparison with experiments, the simulations were run with the properties of water (l ¼ 0:001 Pa s, q ¼ 1000 kg m3) but with those of blood in order to compute haemolysis.
Fig. 3. Pump performance curves. Experimental readings are presented using points and numerical predictions shown with lines. Error bars show experimental uncertainty.
4204
D. Carswell et al. / Applied Mathematical Modelling 37 (2013) 4199–4207
A pressure plot for the device operating at 2000 rpm with an outlet pressure of 3300 Pa is presented in Fig. 4. As expected, low pressure is created on the suction side of the rotors and high pressure, which ultimately drives the flow through the pump, is created on the leading edge of the rotors. Maximal pressure of around 12 000 Pa is present between the rotor leading edge and the case but generally the leading edge pressure is between 2500 and 8000 Pa, i.e., a pressure around 4500 Pa greater than the outlet pressure, leading to a flow rate of 1.35 L min1. It was found that by using either of the k— or k—x turbulence models in the computation of the flow domain did not significantly affect the resulting flow domain with respect to pressure and velocity fields. The only differences are greater shear stresses because of the slight increase in the effective viscosity of the fluid which has been shown to dramatically overestimate haemolysis values, as also observed by Taskin et al. [9]. A total of 861 face elements comprise the inlet to the domain, the origin of the particle tracking procedure. The clustering factor (CF) affects the number of particles that are released from each face. For example, setting CF ¼ 1 equates to 861 particles where setting CF ¼ 20 equates to some 17 000 particles. In this investigation, various values of CF were tested in the range between 1 and 35 to assess the effect it had on haemolysis estimation. A typical set of pathlines, albeit only 100 individual tracks (to facilitate clarity), are shown in Fig. 5. The left of the Figure shows a small recirculation region; particles entering this region went onto exit the pump so the damage accumulated in this region is represented in the final NIH value. The graph shown in Fig. 6 shows the correlation between the number of particles computed and the predicted haemolysis index. Due to the random nature by which particles are positioned, it was to be expected that two different particle runs may not produce exactly the same results as clearly shown by the variation in the two lines. It is clear that, for the particular pump examined here a small number of particles do not accurately capture the haemolysis index of the pump. Analysis of the particles reveals that the resultant haemolysis index is dependant on the number of particles entering regions of high shear stress. Contour plots (Fig. 7) show regions experiencing shear strain rates in the region of 20 000 s1 in the gaps between the rotor and outer casing and between rotor and stator. In the work published by Fraser et al. [32] there was considerable variation between haemolysis estimates using Eq. (13) when compared to a slightly different formation of the damage equation, one estimate being nearly 3 times greater for the same pump. Details of the method used to compute the pathlines with regard to time step and number released are elusive but, on the evidence presented herein, if more pathlines were followed the resulting haemolysis indices may change. Song et al. [8] also present estimations of haemolysis that are lower than those recorded elsewhere, which would lend support to the conclusions from this work that a high threshold number of pathlines must be followed in order to generate statistically significant noise free predictions of the haemolysis index. According to Eq. (13) RBCs can survive the high shear stresses found at the rotor tips in this device so long the residence time is low enough. In this pump, the time taken for a particle to traverse the entirety of the fluid domain was of the order of 10–20 ns. The final predicted NIH value, around 0.017 g/100L, shows a good initial design but needs to be improved before it can be clinically tested. An NIH below 0.01 g/100L is deemed acceptable [33] for long-term implantable devices such as these designed here, although the design criteria is much lower than this. Many, if not all, current pumps exhibit NIH values well below this. Another point to mention should be mesh density. Mesh elements at the inlet are of the order of 0.25 mm in size. A number of particles, depending on the clustering factor, are released at random points within each element. With finer inlet meshes the clustering factor may not need to be as large as those suggested here due to the fact that more elements
Fig. 4. Plots of pressure in the vicinity of the rotor and stator.
D. Carswell et al. / Applied Mathematical Modelling 37 (2013) 4199–4207
4205
Fig. 5. Pathlines in the vicinity of the rotor and stator coloured by velocity magnitude.
Fig. 6. Effect of total particles on haemolysis.
will automatically add more particles into the domain. The evidence from extensive testing of this model is that the key factor in obtaining statistically significant noise free predictions of haemolysis is the total number of particle tracks followed through the domain, and in this context it would appear that 20 000 is adequate. 5. Conclusions In the work reported here a comprehensive computational model based upon advanced CFD simulation technology is developed to enable the prediction of the detailed flow behaviour of and damage to blood flowing through a high speed lLVAD, where the blades are rotating at up to 25 000 r min1. To provide reasonable predictions of the flow rate behaviour the computational model has to incorporate: An accurate transient representation of the rotating section of the pump in a mixed rotating-static mesh code. A very fine mesh, which contains multiple element types to enable the capture of the flow details in the very complex shaped flow domain.
4206
D. Carswell et al. / Applied Mathematical Modelling 37 (2013) 4199–4207
Fig. 7. Contours of shear strain (s1) in the vicinity of the rotor tip.
An efficient parallel cluster implementation to enable very small time steps in tandem with the fine mesh to capture the full vortical nature of the flows throughout the flow domain and provide reasonable predictions of the gross flow rates as well as the fine detail in key areas especially where high shear rates are experienced. The implementation and testing of a sophisticated model for the prediction of haemolysis index based upon the integration of the shear stress-time load experienced by blood cells as they traverse the pump. The combination of the CFD model without an added explicit turbulence model and the haemolysis model appears as a key tool in the accurate prediction of the haemolysis index. Indeed, it is argued here that the case for explicitly including a model of turbulence is not clear, since model tests with a two-equation turbulence model have lead to a significant overestimation of haemolysis levels. It is clearly demonstrated that in order for the combined CFD and haemolysis index model to make reliable and repeatable predictions of the haemolysis index in the class of lLVAD reported here then upwards of 20 000 particle track flowlines are required to ensure sufficiently statistically significant numbers of tracks through the high shear rate regions (> 20 000 s1) are generated. Finally, now that the basic computational modelling framework has been developed and evaluated it is now being exploited to improve the design of a range of lVAD blood pumps through the minimising high shear rate areas and the flow of blood through those areas to minimise haemolysis, where at each stage further validation will be performed to provide increased confidence in its predictive capability. Acknowledgement This work was funded by the Technology Strategy Board (TSB) Grant No. TP/8/ADM/6/I/Q2065G entitled Multi-Physics Modelling of a Cardio Circulatory Support System. The authors would like to acknowledge the help of Kevin Fernquest, Alex O’Malley and Andy Hilton for their collaboration and without whom the project would not be possible. The authors would also like to thank Calon-Cardio Technology Ltd. for their permission to publish this work. References [1] J. Woodcock, Physical properties of blood and their influence on blood-flow measurement, Rep. Progr. Phys. 39 (1976) 65–127. [2] A. Hoffbrand, J. Pettit, P. Moss, Essential Haematology, Blackwell Science, 2004. [3] S.S. Lee, K. Ahn, S.J. Lee, K. Sun, P. Goedhart, M. Hardeman, Shear induced damage of red blood cells monitored by the decrease of their deformability, Korea-Aust. Rheol. J. 16 (3) (2004) 141–146. [4] M. Behbahani, M. Behr, M. Hormes, U. Steinseifer, D. Arora, O. Coronado, M. Pasquali, A review of computational fluid dynamics analysis of blood pumps, Eur. J. Appl. Math. 20 (2009) 363–397. [5] P. Blackshear, F. Dorman, J. Steinbach, Some mechanical effects that influence haemolysis, Trans. Am. Soc. Artif. Int. Organs 11 (1965) 112. [6] L. Goubergrits, K. Affeld, Numerical estimation of blood damage in artificial organs, Artif. Org. 28 (5) (2003) 499–507. [7] M. Giersiepen, L. Wurzinger, R. Optiz, H. Reul, Estimation of shear stress-related blood damage in heart valve prostheses - in vitro comparison of 25 aortic valves, Int. J. Artif. Organs 13 (1990) 300–306. [8] X. Song, A. Throckmorton, H. Wood, J. Antaki, D. Olsen, Computational fluid dynamics prediction of blood damage in a centrifugal pump, Artif. Organs 27 (10) (2003) 938–941. [9] M. Taskin, K. Fraser, T. Zhang, B. Gellman, A. Fleischli, K. Dasse, B. Griffith, Z. Wu, Computational characterization of flow and hemolytic performance of the UltraMag blood pump for circulatory support, Artif. Organs 34 (12) (2010) 1099–1113. [10] L. Leverett, J. Hellums, C. Alfrey, E. Lynch, Red blood cell damage by shear stress, Biophys. J. 12 (1972) 257–273. [11] P. Lu, H. Lai, J. Liu, A reevaluation and discussion on the threshold limit for hemolysis in a turbulent shear flow, J. Biomech. 34 (2001) 1361–1364. [12] R. Paul, J. Apel, S. Klaus, F. Schugner, P. Schwindke, H. Reul, Shear stress related blood damage in laminar couette flow, Artif. Organs 27 (6) (2003) 517– 529.
D. Carswell et al. / Applied Mathematical Modelling 37 (2013) 4199–4207
4207
[13] M. Grigioni, C. Daniele, U. Morbiducci, G. D’Avenio, G. Benedetto, V. Barbaro, The power-law mathematical model for blood damage prediction: Analytical developments and physical inconsistencies, Artif. Organs 28 (5) (2004) 467–475. [14] N. Watanabe, O. Karsak, F. Neudel, T. Kink, J. Apel, T. Fujimoto, H. Reul, S. Takatani, Simulation of the BP-80 blood pump, Artif. Organs 25 (9) (2001) 733–739. [15] A. Throckmorton, A. Untaroiu, P. Allaire, H. Wood, G. Matherne, D. Lim, B. Peeler, D. Olsen, Computational analysis of an axial flow pediatric ventricular assist device, Artif. Organs 28 (10) (2004) 881–891. [16] A. Throckmorton, A. Unaroiu, D. Lim, H. Wood, P. Allaire, Fluid force predictions and experimental measurements for a magnetically levitated pediatric ventricular assist device, Artif. Organs 31 (5) (2007) 359–368. [17] W.-K. Chan, Y.-W. Wong, W. Ong, S.-Y. Koh, V. Chong, Numerical investigation of the effects of the clearance gap between the inducer and impeller of an axial blood pump, Artif. Organs 29 (3) (2005) 250–258. [18] L. Chua, B. Su, T. Lim, T. Zhou, Numerical simulation of an axial blood pump, Artif. Organs 31 (7) (2007) 560–570. [19] Y.-J. Li, F.-J. Wang, Numerical investigation of performance of an axial-flow pump with inducer, J. Hydrodyn. 19 (6) (2007) 705–711. [20] Y. Zhang, S. Xue, X.-M. Gui, H.-S. Sun, H. Zhang, X.-D. Zhu, S.-S. Hu, A novel integrated rotor of axial blood flow pump designed with computational fluid dynamics, Artif. Organs 31 (7) (2007) 580–585. [21] D. Arora, M. Behr, O. Coronado-Matutti, M. Pasquali, Estimation of hemolysis in centrifugal blood pumps using morphology tensor approach, Third MIT Conference on Computational Fluid and Solid Mechanics (2005) 578–582. [22] A. Arvand, M. Hormes, H. Reul, A validatied computational fluid dynamics model to estimate hemolysis in a rotary blood pump, Artif. Organs 29 (7) (2005) 531–540. [23] D. Arora, Computational Hemodynamics: Hemolysis and Viscoelasticity, Ph.D. thesis, School of Engineering, Rice University, 2005. [24] Physica,
Contents lists available at SciVerse ScienceDirect
Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm
A CFD model for the prediction of haemolysis in micro axial left ventricular assist devices D. Carswell a,⇑, D. McBride a, T.N. Croft a, A.K. Slone a, M. Cross a, G. Foster b a b
Civil and Computational Engineering Centre, College of Engineering, Swansea University, Swansea, West Glamorgan, Wales SA2 8PP, UK Calon Cardio-Technology Ltd, Institute of Life Science, Swansea University, Swansea, West Glamorgan, Wales SA2 8PP, UK
a r t i c l e
i n f o
Article history: Received 4 October 2011 Received in revised form 13 July 2012 Accepted 10 September 2012 Available online 23 September 2012 Keywords: Haemolysis CFD Blood pump Axial flow
a b s t r a c t This paper describes the development and application of a computational model based upon Computational Fluid Dynamics (CFD) software simulation technology to predict haemolysis in micro Left Ventricular Assist Devices (lLVAD). A CFD model, capturing the full three dimensional geometry of the device, together with an explicit representation of the rotating machinery based upon a rotating reference frame, is solved transiently. Mixed meshes with the order of a million elements are required to resolve the flow adequately and so to enable solutions in a reasonable time (e.g., 3 h) the model is solved on a high performance parallel cluster. Haemolysis is a measure of damage occurring in the blood and is conceived as accumulating as it passes through parts of the device where it encounters high shear forces. As such, the haemolysis model is based upon tracking the behaviour of particles released at the inlet throughout the flow domain and calculating the damage accumulated by each individual particle as it traverses the device. In order to ensure the model predictions of haemolysis are noise free from a statistically significant perspective then it is demonstrated that the number of particles to be tracked must exceed 20 000 in any simulation experiment. Comparisons with experimental data from a companion paper demonstrate the effectiveness of the CFD simulation embedding the haemolysis model. Ó 2012 Elsevier Inc. All rights reserved.
1. Introduction Human, and mammal, red blood cells (RBCs) exist to transport oxygen to every living cell in the body. RBC structure comprises an elastic membrane that encases a volume of haemoglobin in a bi-concave disc of diameter to the order of 8 lm [1]. The membrane comprises a lipid external layer, internal proteins and a cytoskeleton [2]. The cells may deform in different ways under different conditions, usually dependant on shear stress. In healthy human veins and arteries, the shear stresses experienced may range from 0.1 to 20 Pa, increasing in the smaller capillaries [3]. In high shear stress regions the RBC membrane cannot sustain its deformed shape for long and eventually loses its elasticity resulting in mechanical damage. As membrane damage occurs, haemoglobin may escape from the RBC and into the blood plasma. Haemolysis expresses the amount of haemoglobin that escapes in this manner [4]. A relationship linking shear stress (s) and exposure time (t) to plasma-free haemoglobin (pfHb) first proposed by Blackshear et al. [5] is
pfHb ¼ ts2
⇑ Corresponding author. E-mail address: [email protected] (D. Carswell). 0307-904X/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.apm.2012.09.020
ð1Þ
4200
D. Carswell et al. / Applied Mathematical Modelling 37 (2013) 4199–4207
but only tested in the range 1 < s < 10Pa for 1 h. However, the same relationship was found to hold true for s < 8000Pa and t < 102 s [6]. Hb A more general form expressing the proportion of haemoglobin released DHb is related to s and t by
DHb ¼ Ct a sb : Hb
ð2Þ
The values of the constants a; b and C proposed by Giersiepen et al. [7] were 0.785, 2.416 and 3:62 107 respectively. A slightly different set of parameters were cited by Song et al. [8], these being 0.765, 1.991 and 1:8 106 in their CFD simulations. More recently, a different set of parameters (0.747, 2.004 and 1:21 107 respectively) derived for ovine blood 0 1=a were published by Taskin et al. in their CFD simulation [9]. Taskin proposes using a convected variable Hb ¼ Hb to accumulate haemolysis through a pump such that 0
0
dðDHb Þ @ DHb þ qu ¼ S; dt @x
ð3Þ
where the source term S ¼ qðHbC sb Þ1=a . As boundary conditions, it is assumed Hb ¼ 0 at the inlet and zero flux is specified at all solid walls. Using a convected variable such as this it is possible to quickly identify areas of a pump that are more haemolysing than others - something that is not always possible with simply following pathlines. There still appears to be some debate on the upper bound of the safe shear stress, i.e., the level below which haemolysis does not occur or is insignificant. Previous research suggests this value is within the range of 150–800 Pa [10–13]. Computational Fluid Dynamics (CFD) has long be used to aid designers obtain a greater understanding of, not just heart pump operation, but many other applications. The detailed insight it gives may be difficult or nearly impossible to obtain experimentally. Many researchers have already used CFD as a means of analysing both axial and radial heart pumps. Watanabe et al. [14] analysed a centrifugal heart pump using TASCflow (AEA, Pittsburg, USA) utilising fully structured meshes and comparisons were made with experimental data to establish the correctness of the model. Throckmorton et al. [15] used Ansys CFX to evaluate a pump design prior to manufacture to make sure that it would meet the requirements. They then went onto analyse a candidate pump in more detail [16] to ensure the magnetic bearings would be able to support the force exerted on the rotor due to hydrodynamic forces. Another axial pump was also analysed experimentally and numerically, again using TASCFlow, by Chan et al. [17]. Fluent 6.2 was used by Chua et al. [18] in their computational model making use of sliding meshes but has not been compared to experimental data to the best of the authors’ knowledge. An axial pump containing both an inducer and straightener was analysed by Li and Wang [19] using the RNG k- turbulence model and rotating frames of reference. They concluded that inducers can have quite a large effect on pump performance. The NUMECA commercial CFD software was used to perform computations on the axial pump of Zhang et al. [20] but not compared to experimental data. It is clear from the above literature that compare numerical and experimental results that, given an adequate mesh, commercial CFD-based simulation software are capable of matching experimental data. Using pathlines is the general technique used to capture data required for haemolysis and has been employed by a variety of other researchers [21,22,8] as it is, so far, the only way to generate appropriate data for haemolysis estimation. Essentially, the CFD solver produces a flow field through the blood pump and then mass-less particles are tracked through this flow field to assess the damage they experience as a consequence of the shear stress exerted by the characteristics of the flow at any and every location visited by a specific particle. Of course, the data from a large number of particle tracks has to be integrated to provide statistically significant information. To compare haemolysis predictions with experimental data it is often necessary to convert both datasets to the Normalised Index of Haemolysis (NIH). Experimentally, NIH is determined by
NIH ¼ H V ð1 Hct Þ
100 ; Q T
ð4Þ
where V is the volume of blood in the test loop (L), Hct is the haematocrit, Q is the flow rate (L min1) and T time (min). Values of T and V are not available numerically so a slightly different formulation is used
NIH ¼ 100
DHb ð1 Hct Þ j; Hb
ð5Þ
where j is the haemoglobin content of blood. Typical values of Hct and j are 45% and 150 g L1 [23]. This paper describes the computational models employed by the authors to characterise the flow behaviour of lLVADs and to produce predictions of expected performance curves and haemolysis levels. The outcomes of the models are currently underwriting the work of designers in advancing future pump designs by allowing an in-depth analysis of the blood flow behaviour in current pump geometries. 2. Pump design The lLVADs currently being designed are destined as short term devices. That is, they are designed to produce flow rates of around 2.5 L min1 to work alongside a recovering heart that is still functional. Other pumps are designed on the bridgeto-recovery principle or as more permanent devices in which case flow rates beyond 5 L min1 are required.
D. Carswell et al. / Applied Mathematical Modelling 37 (2013) 4199–4207
4201
Fig. 1. Pump geometry considered in this paper.
In this paper, the lLVAD shown in Fig. 1 is considered measuring 7.5 mm in diameter by 40 mm total length. Three blades are 1.5 mm in height with a clearance of just 20 lm from the case. The main inlet is located immediately upstream of the rotor. A three-vaned rotor serves to provide the pumping force and the stator converts whirl energy to a pressure head. Having the main inlet in the side of the pump means the main flow path through the pump is as short as possible. The motor and hydrodynamic bearings are located at the rear of the pump where a lubrication channel is also present. Using blood as the lubricating agent eliminates the need to use other lubricants which also eliminates the need to keep them separate from the blood, hence the design is simplified. The pump is designed to operate between 15,000 and 30,000 r min1. Pump performance is evaluated at a variety of operating pressures and rotation speeds, often outside the anticipated working range to assess the full range of physiological conditions. In laboratory tests, the pump was connected via lengths of pipe to a reservoir, the outlet tubing was restricted via a clamp to provide a pressure at the outlet. Flow rates and inlet and outlet pressures were recorded experimentally. Experiments were conducted using water. 3. Computational models and methods The 3D incompressible Navier–Stokes equations are used to predict the flow field through the pump. In line with other researchers [14–20] we assume the blood flow to be incompressible. Rotating machinery requires special consideration in the formulation of the equations. All rotating parts are enclosed in a rotating domain modelled using a rotating reference frame. In the rotating frame, a relative velocity ur is related to the absolute velocity u, position vector r and angular velocity vector x by ur ¼ u þ x r. Appropriate substitutions into the convection term of the Navier–Stokes equations must be made in the rotating frame by either solving the relative or absolute velocities. In this application the latter was chosen yielding a new formulation of the governing equations
q
@u þ qur ru þ qðx uÞ ¼ rp þ rðlruÞ þ Su : @t
ð6Þ
Mass is also conserved using the continuity equation (Eq. (7)) but requires no special treatment.
qru ¼ Sm ;
ð7Þ
Using the rotor diameter as a characteristic length (‘), Reynolds numbers,
Re ¼
qv ‘ l
ð8Þ
were computed to be never greater than 5000 and for most of the time and space considerably less, so that while the global flow is just outside the laminar flow regime it is by no means fully turbulent. Blood density was kept constant at 1050 kg m3 and the Quemada viscosity model was used to model blood viscosity:
l¼
l
p 1 k H Þ2 2 Q ct
;
ð9Þ
pffiffiffiffiffi k0 þ k1 cr c pffiffiffiffiffi and cr ¼ : cc 1 þ cr
ð10Þ
ð1
where
kQ ¼
Blood plasma viscosity, lp ¼ 1:23 m Pas and the level of haematocrit, Hct ¼ 48%, for a healthy person. Shear strain rate c is obtained from the flow field and k0 and k1 are model constants that represent the rate of RBC rouleaux formation at zero and infinite shear. When it came to considerations of turbulence, a combination of the following reasons: Most of the flow is at worst just bordering on the turbulent regime and much of it is truly laminar, The viscosity is actually non-Newtonian and so someway outside the scope of the original conception of conventional computational turbulence models, and
4202
D. Carswell et al. / Applied Mathematical Modelling 37 (2013) 4199–4207
The very fine meshes and small time steps enables capture of many of the larger more significant vortices led to the conclusion that there was no real justification for employing a formal turbulence model since all were inadequate in some essential manner. Moreover, since the fine meshes and small time steps captured much of key vortex behaviour it is argued that much of the developing turbulent behaviour is represented in this computational implementation. The equations are solved using the in-house code PHYSICA [24] using the vertex-based finite volume discretisation method [25]. In this scheme flow variables are stored at nodes and the finite element shape functions are used in the computation of face fluxes. For any variable / the face gradients are calculated using the shape function gradients such that
@/ @xi
f
n X @Nj j¼1
@xi
/j ;
ð11Þ
where N j is the shape function associated with node j. Face values /f are computed from the nodal values and there are a number of high-order schemes able to do this. In this work a second-order Mass-Weighted (MAW) upwind scheme [26] was implemented for unstructured 3D meshes based on the work of Tran et al. [27] Solving the Navier–Stokes equations requires the use of a pressure–velocity coupling. For co-located variables, the method of Prakash and Patankar [28] was implemented. As the flow rate generated by the rotor is sought, only pressure boundary conditions could be prescribed at both the inlets and outlet. Both inlets were give 0 Pa pressure and the outlet pressure was varied between 3300 Pa and 200 000 Pa to simulate different physiological conditions. A velocity is present on the rotor, in the inertial reference frame, by virtue of the rotating frame. The outer case and stator were specified as no-slip wall conditions. An important part of the overall computational procedure was the generation of meshes of suitable quality for the CFD model where part of the mesh is rotating. The meshing strategy employed a mix of element types. The nature of the device presented a number of specific meshing challenges. In regions of small gaps it was found that it was easier to control the mesh quality and density by using structured, hexahedral elements rather than tetrahedra. Ansys ICEM CFD 13 (Ansys Inc.) was employed as the mesh generation software. Creation of hexahedral elements made use of ICEM’s blocking capabilities. Mesh element quality was kept above 0.4. The cone of the rotor was separated from the stator by only 0.1 mm. Knowing that this gap had to accommodate an interface between the rotating and stationary domains meant that the mesh in this region was, perhaps, finer than it needed to be. The main flow path was meshed using tetrahedra of size 0.1 mm except in the gap between rotor and case in which case the element size was reduced to provide sufficient resolution of the flow in this gap. A range of meshes were built from coarse (500 000 elements, 0.5 mm element size) to fine (2 000 000 elements, 0.1 mm element size) to assess any mesh dependance issues. A coarse mesh is shown in Fig. 2 where the darker region, to the right of the figure, indicates the elements belonging to the rotating reference frame. The figure has been chopped to show the region of interest as the computational inlets and outlets are positioned further away from the actual inlets and outlets so as to not cut short any flow features that may be present. Results showed a semi-fine mesh (1 000 000 elements, 0.2 mm element size) was adequate for the simulation where the predicted flow rate changed by <1%. The solution was deemed converged when residuals of mass and the three velocity components had fallen by three orders of magnitude. The host CFD code is fully parallelised using a distributed memory model [29] and well adapted to run in parallel for various rotating machinery representations [30]. Computations were run in parallel using MPI as the interface between the different processes using two quad-core 64-bit Intel Core i7 chips. A transient analysis was run over a period of 30 time steps with the interval being 0.005 s which was enough for the simulation to produce a steady flow rate at the outlet. One simulation took approximately 3 h to complete. 3.1. Pathlines Once a converged flow field had been computed, pathlines of mass-less particles were tracked through it. A number of pathlines were released at random points within each element at the inlet. The number of particles released per element
Fig. 2. Coarse mesh of a typical pump. The darker elements show the region of the mesh in the rotating reference frame.
D. Carswell et al. / Applied Mathematical Modelling 37 (2013) 4199–4207
4203
is known as the clustering factor and ranges from 1 to 35 resulting between 900 and 30,000 particles. Each step along pathline i is computed from the previous step cni and the current velocity field u by
cnþ1 ¼ cni þ Dt i
@ui : @xi
ð12Þ
A constant time step Dt ¼ 1 105 s was chosen to ensure accurate pathlines. Element sizes and velocities showed the maximum time step permissible was 5 105 s meaning the time step used was suitable. A greater timestep did not provide sufficient accuracy due to the velocities and cell sizes involved. Smaller time steps did not produce any noticeable difference in the computed paths except that it consumed a lot more processing power. A particle ends when it reaches the outlet. Another undesirable condition that causes particles to end is when they hit a wall. One method would be to reflect the particle using the boundary face normal. Occasionally, particles may become trapped if they encounter recirculation regions, therefore never exiting the pump. Only particles that reached the outlet are considered in the haemolysis computations. 3.2. Haemolysis While computing the pathlines, the shear stress is also output at each time step as it is one of the factors influencing haemolysis. The time step at each point along each pathline was used as the exposure time (dt) and the shear stress (s) was assumed constant throughout this period. A number of inconsistencies were found by relying on Eq. (2) to compute haemolysis [31]. In particular it does not account for load history, that is, two particles reaching the same point in the flow via different pathlines may accumulate different levels of damage in accordance with damage previously experienced. This is due to the relaxation time a stresses RBC membrane has once a stress is removed or reduced. Grigioni et al. [31] proposed
" #a1 n i X X b=a Ca sðtj Þ dt þ Dðt0 Þ sðti Þb=a dt i¼1
ð13Þ
j¼1
as a more discriminating means of computing blood damage. The inner sum accounts for the load history where n is the number of steps per particle and the function D allows an initial damage to be specified. Note that s is a function that varies with position along the pathline. In this application, all particles are assumed to have zero initial damage, therefore, Dðt 0 Þ ¼ 0. Damage is then converted to NIH using Eq. (5). 4. Results and discussion Pump performance characteristics are presented in Fig. 3. At the higher rotation rates, the pressure generated by the rotor blades becomes much greater than the outlet pressure and, hence, the curves begin to converge. The predicted flow rates for 0 Pa begin to diverge from the experimental results as flow rate increases. The reason for this is thought to be extra losses caused by the length of pipe used experimentally which was not modelled computationally. It should be noted that for comparison with experiments, the simulations were run with the properties of water (l ¼ 0:001 Pa s, q ¼ 1000 kg m3) but with those of blood in order to compute haemolysis.
Fig. 3. Pump performance curves. Experimental readings are presented using points and numerical predictions shown with lines. Error bars show experimental uncertainty.
4204
D. Carswell et al. / Applied Mathematical Modelling 37 (2013) 4199–4207
A pressure plot for the device operating at 2000 rpm with an outlet pressure of 3300 Pa is presented in Fig. 4. As expected, low pressure is created on the suction side of the rotors and high pressure, which ultimately drives the flow through the pump, is created on the leading edge of the rotors. Maximal pressure of around 12 000 Pa is present between the rotor leading edge and the case but generally the leading edge pressure is between 2500 and 8000 Pa, i.e., a pressure around 4500 Pa greater than the outlet pressure, leading to a flow rate of 1.35 L min1. It was found that by using either of the k— or k—x turbulence models in the computation of the flow domain did not significantly affect the resulting flow domain with respect to pressure and velocity fields. The only differences are greater shear stresses because of the slight increase in the effective viscosity of the fluid which has been shown to dramatically overestimate haemolysis values, as also observed by Taskin et al. [9]. A total of 861 face elements comprise the inlet to the domain, the origin of the particle tracking procedure. The clustering factor (CF) affects the number of particles that are released from each face. For example, setting CF ¼ 1 equates to 861 particles where setting CF ¼ 20 equates to some 17 000 particles. In this investigation, various values of CF were tested in the range between 1 and 35 to assess the effect it had on haemolysis estimation. A typical set of pathlines, albeit only 100 individual tracks (to facilitate clarity), are shown in Fig. 5. The left of the Figure shows a small recirculation region; particles entering this region went onto exit the pump so the damage accumulated in this region is represented in the final NIH value. The graph shown in Fig. 6 shows the correlation between the number of particles computed and the predicted haemolysis index. Due to the random nature by which particles are positioned, it was to be expected that two different particle runs may not produce exactly the same results as clearly shown by the variation in the two lines. It is clear that, for the particular pump examined here a small number of particles do not accurately capture the haemolysis index of the pump. Analysis of the particles reveals that the resultant haemolysis index is dependant on the number of particles entering regions of high shear stress. Contour plots (Fig. 7) show regions experiencing shear strain rates in the region of 20 000 s1 in the gaps between the rotor and outer casing and between rotor and stator. In the work published by Fraser et al. [32] there was considerable variation between haemolysis estimates using Eq. (13) when compared to a slightly different formation of the damage equation, one estimate being nearly 3 times greater for the same pump. Details of the method used to compute the pathlines with regard to time step and number released are elusive but, on the evidence presented herein, if more pathlines were followed the resulting haemolysis indices may change. Song et al. [8] also present estimations of haemolysis that are lower than those recorded elsewhere, which would lend support to the conclusions from this work that a high threshold number of pathlines must be followed in order to generate statistically significant noise free predictions of the haemolysis index. According to Eq. (13) RBCs can survive the high shear stresses found at the rotor tips in this device so long the residence time is low enough. In this pump, the time taken for a particle to traverse the entirety of the fluid domain was of the order of 10–20 ns. The final predicted NIH value, around 0.017 g/100L, shows a good initial design but needs to be improved before it can be clinically tested. An NIH below 0.01 g/100L is deemed acceptable [33] for long-term implantable devices such as these designed here, although the design criteria is much lower than this. Many, if not all, current pumps exhibit NIH values well below this. Another point to mention should be mesh density. Mesh elements at the inlet are of the order of 0.25 mm in size. A number of particles, depending on the clustering factor, are released at random points within each element. With finer inlet meshes the clustering factor may not need to be as large as those suggested here due to the fact that more elements
Fig. 4. Plots of pressure in the vicinity of the rotor and stator.
D. Carswell et al. / Applied Mathematical Modelling 37 (2013) 4199–4207
4205
Fig. 5. Pathlines in the vicinity of the rotor and stator coloured by velocity magnitude.
Fig. 6. Effect of total particles on haemolysis.
will automatically add more particles into the domain. The evidence from extensive testing of this model is that the key factor in obtaining statistically significant noise free predictions of haemolysis is the total number of particle tracks followed through the domain, and in this context it would appear that 20 000 is adequate. 5. Conclusions In the work reported here a comprehensive computational model based upon advanced CFD simulation technology is developed to enable the prediction of the detailed flow behaviour of and damage to blood flowing through a high speed lLVAD, where the blades are rotating at up to 25 000 r min1. To provide reasonable predictions of the flow rate behaviour the computational model has to incorporate: An accurate transient representation of the rotating section of the pump in a mixed rotating-static mesh code. A very fine mesh, which contains multiple element types to enable the capture of the flow details in the very complex shaped flow domain.
4206
D. Carswell et al. / Applied Mathematical Modelling 37 (2013) 4199–4207
Fig. 7. Contours of shear strain (s1) in the vicinity of the rotor tip.
An efficient parallel cluster implementation to enable very small time steps in tandem with the fine mesh to capture the full vortical nature of the flows throughout the flow domain and provide reasonable predictions of the gross flow rates as well as the fine detail in key areas especially where high shear rates are experienced. The implementation and testing of a sophisticated model for the prediction of haemolysis index based upon the integration of the shear stress-time load experienced by blood cells as they traverse the pump. The combination of the CFD model without an added explicit turbulence model and the haemolysis model appears as a key tool in the accurate prediction of the haemolysis index. Indeed, it is argued here that the case for explicitly including a model of turbulence is not clear, since model tests with a two-equation turbulence model have lead to a significant overestimation of haemolysis levels. It is clearly demonstrated that in order for the combined CFD and haemolysis index model to make reliable and repeatable predictions of the haemolysis index in the class of lLVAD reported here then upwards of 20 000 particle track flowlines are required to ensure sufficiently statistically significant numbers of tracks through the high shear rate regions (> 20 000 s1) are generated. Finally, now that the basic computational modelling framework has been developed and evaluated it is now being exploited to improve the design of a range of lVAD blood pumps through the minimising high shear rate areas and the flow of blood through those areas to minimise haemolysis, where at each stage further validation will be performed to provide increased confidence in its predictive capability. Acknowledgement This work was funded by the Technology Strategy Board (TSB) Grant No. TP/8/ADM/6/I/Q2065G entitled Multi-Physics Modelling of a Cardio Circulatory Support System. The authors would like to acknowledge the help of Kevin Fernquest, Alex O’Malley and Andy Hilton for their collaboration and without whom the project would not be possible. The authors would also like to thank Calon-Cardio Technology Ltd. for their permission to publish this work. References [1] J. Woodcock, Physical properties of blood and their influence on blood-flow measurement, Rep. Progr. Phys. 39 (1976) 65–127. [2] A. Hoffbrand, J. Pettit, P. Moss, Essential Haematology, Blackwell Science, 2004. [3] S.S. Lee, K. Ahn, S.J. Lee, K. Sun, P. Goedhart, M. Hardeman, Shear induced damage of red blood cells monitored by the decrease of their deformability, Korea-Aust. Rheol. J. 16 (3) (2004) 141–146. [4] M. Behbahani, M. Behr, M. Hormes, U. Steinseifer, D. Arora, O. Coronado, M. Pasquali, A review of computational fluid dynamics analysis of blood pumps, Eur. J. Appl. Math. 20 (2009) 363–397. [5] P. Blackshear, F. Dorman, J. Steinbach, Some mechanical effects that influence haemolysis, Trans. Am. Soc. Artif. Int. Organs 11 (1965) 112. [6] L. Goubergrits, K. Affeld, Numerical estimation of blood damage in artificial organs, Artif. Org. 28 (5) (2003) 499–507. [7] M. Giersiepen, L. Wurzinger, R. Optiz, H. Reul, Estimation of shear stress-related blood damage in heart valve prostheses - in vitro comparison of 25 aortic valves, Int. J. Artif. Organs 13 (1990) 300–306. [8] X. Song, A. Throckmorton, H. Wood, J. Antaki, D. Olsen, Computational fluid dynamics prediction of blood damage in a centrifugal pump, Artif. Organs 27 (10) (2003) 938–941. [9] M. Taskin, K. Fraser, T. Zhang, B. Gellman, A. Fleischli, K. Dasse, B. Griffith, Z. Wu, Computational characterization of flow and hemolytic performance of the UltraMag blood pump for circulatory support, Artif. Organs 34 (12) (2010) 1099–1113. [10] L. Leverett, J. Hellums, C. Alfrey, E. Lynch, Red blood cell damage by shear stress, Biophys. J. 12 (1972) 257–273. [11] P. Lu, H. Lai, J. Liu, A reevaluation and discussion on the threshold limit for hemolysis in a turbulent shear flow, J. Biomech. 34 (2001) 1361–1364. [12] R. Paul, J. Apel, S. Klaus, F. Schugner, P. Schwindke, H. Reul, Shear stress related blood damage in laminar couette flow, Artif. Organs 27 (6) (2003) 517– 529.
D. Carswell et al. / Applied Mathematical Modelling 37 (2013) 4199–4207
4207
[13] M. Grigioni, C. Daniele, U. Morbiducci, G. D’Avenio, G. Benedetto, V. Barbaro, The power-law mathematical model for blood damage prediction: Analytical developments and physical inconsistencies, Artif. Organs 28 (5) (2004) 467–475. [14] N. Watanabe, O. Karsak, F. Neudel, T. Kink, J. Apel, T. Fujimoto, H. Reul, S. Takatani, Simulation of the BP-80 blood pump, Artif. Organs 25 (9) (2001) 733–739. [15] A. Throckmorton, A. Untaroiu, P. Allaire, H. Wood, G. Matherne, D. Lim, B. Peeler, D. Olsen, Computational analysis of an axial flow pediatric ventricular assist device, Artif. Organs 28 (10) (2004) 881–891. [16] A. Throckmorton, A. Unaroiu, D. Lim, H. Wood, P. Allaire, Fluid force predictions and experimental measurements for a magnetically levitated pediatric ventricular assist device, Artif. Organs 31 (5) (2007) 359–368. [17] W.-K. Chan, Y.-W. Wong, W. Ong, S.-Y. Koh, V. Chong, Numerical investigation of the effects of the clearance gap between the inducer and impeller of an axial blood pump, Artif. Organs 29 (3) (2005) 250–258. [18] L. Chua, B. Su, T. Lim, T. Zhou, Numerical simulation of an axial blood pump, Artif. Organs 31 (7) (2007) 560–570. [19] Y.-J. Li, F.-J. Wang, Numerical investigation of performance of an axial-flow pump with inducer, J. Hydrodyn. 19 (6) (2007) 705–711. [20] Y. Zhang, S. Xue, X.-M. Gui, H.-S. Sun, H. Zhang, X.-D. Zhu, S.-S. Hu, A novel integrated rotor of axial blood flow pump designed with computational fluid dynamics, Artif. Organs 31 (7) (2007) 580–585. [21] D. Arora, M. Behr, O. Coronado-Matutti, M. Pasquali, Estimation of hemolysis in centrifugal blood pumps using morphology tensor approach, Third MIT Conference on Computational Fluid and Solid Mechanics (2005) 578–582. [22] A. Arvand, M. Hormes, H. Reul, A validatied computational fluid dynamics model to estimate hemolysis in a rotary blood pump, Artif. Organs 29 (7) (2005) 531–540. [23] D. Arora, Computational Hemodynamics: Hemolysis and Viscoelasticity, Ph.D. thesis, School of Engineering, Rice University, 2005. [24] Physica,