Flow in a Ventricular Assist Device - Pump Performance & Blood Damage Prediction

Application Challenge AC7-03 © copyright ERCOFTAC 2022

Description

Introduction

Ventricular Assist Devices (VADs) are implanted in patients with severe heart failure. Today, nearly all VADs are designed as turbomachines, since they have a higher power density than pulsatile pumps. Therefore they can be implanted within the human body. Compared to Total Artificial Hearts (TAHs), the VADs do not replace the heart, but assist a weak heart by creating the pressure needed for supplying sufficient blood flow in the circulatory system. A few examples of VADs and TAHs currently implanted in patients can be seen in Figure 1.1.

By using Computational Fluid Dynamics (CFD), VADs must be designed and optimised in such a way that they produce a physiological pressure increase sufficient to supply the body with enough blood flow. Furthermore, they must be designed so that the blood passing the VAD is not damaged due to non-physiological flow conditions (e.g., high shear stresses, stagnation areas and high turbulent kinetic energies (TKE)) in the device.

The flow simulation in a VAD is challenging, since the inflow is laminar and all turbulence is produced within the pump and decays shortly after the pump outlet. In addition, complex interactions of secondary flows occur in the turbomachine, which must be taken into account in the flow simulation, since they influence the acting stresses in the flow of the VAD [2].

In this respect, the aim of this study is to investigate the suitability of the URANS method for the flow computation in an axial VAD. Here, both, fluid mechanical parameters, such as the pump characteristics, as well as hemodynamic parameters, such as the hemolysis index MIH, are investigated. The stress fields of the URANS will be compared with a highly turbulence-resolving large-eddy simulation, which represents the reference case for comparison. Furthermore, the influence of the grid resolution in the URANS computations will also be investigated based on an extended grid study.

Fig.1.1 Total Artificial Hearts (TAHs) and Ventricular Assist Devices (VADs), which are currently in clinical use. One TAH, one pulsatile VAD and three turbo pumps are shown.

Relevance to Industrial Sector

The flow computation in a Ventricular Assist Device (VAD) is an important procedure for the VAD design and optimization in the pre-clinical evaluation. The aims of these CFD studies are, on one hand, to guarantee that the VAD offers the physiologically relevant pressure increase at the chosen operating points to sufficiently support the VAD patient. On the other hand, hemodynamical parameters must be evaluated in these studies. Here, it is important that the CFD indicates relevant regions for potential blood damage or thrombi formation so that these regions can be minimized in the optimization procedure. Additionally, CFD is important for comparing different designs in order to find the pump with the highest hemocompatibility (lowest blood damage).

When the VAD designer is able to find a good VAD design by CFD, the amount of in-vitro (experimental test of pump performance or hemolysis (red blood cell damage)) and in-vivo testing (animal trials) might be reduced.

Flow Domain Geometry

The geometry of the VAD is available at the link: https://unibox.uni-rostock.de/getlink/fi8AE4mYS4kxu8ZY51oxDh/

The investigated axial VAD is illustrated within the computational domain in Fig. 1.2. It consist of a spinning rotor with two blade channels; a stationary, inlet guide vane with five blades; and a stationary, three-bladed outlet guide vane. While the transfer of angular momentum via the rotor leads to an increase in pressure and velocity (swirl formation), the swirl is converted back into static pressure by the diffuser effect in the outlet guide vane. The inlet guide vane is intended for flow manipulation of the rotor inflow. In the present VAD, a counter-rotation is generated, which theoretically leads to a flatter pressure head curve.

This axial VAD has been designed at the Institute of Turbomachinery at the University of Rostock. The design was inspired by axial VADs, which are currently in clinical use. The design principles are briefly explained: After chosing the inner ( $D_{1}=11~mm$ ) and outer ( $D_{2}=16~mm$ ) diameter, meridian lines were placed between hub and casing to set the blades angles for a chosen nominal operation point $(Q=4.5~l/min,~n=7900~r/min)$ . Afterward, the wrap angle of the blades was adjusted to obtain an even blade progression. After that, the blade thickness and gap height (distance between the casing and the rotor's blade tip; $120~\mu m$ ) were included. The outlet guide vanes were set up so that the swirl is reduced as much as possible. Finally, the coupled rotor and stator were adjusted using URANS simulations, until the pump achieved a pressure head of $H=(70-80)~mmHg$ and a hydraulic effiency of $\eta _{i}\approx 33\%$ at the nominal operating point.

This VAD was designed for reseach purposes only. The aim of the design was not to build a "perfect" implantable VAD with an optimal flow behaviour, but rather a pump, in which the same flow pattern can occur as in real implanted VADs. In these pumps, the inflow angle cannot be optimal during the entire operation time due to varying input from the remaining heart activity. Thus, non-optimal inflow angles and flow paths were delibaretely accepted at the nominal operation point.

Fig.1.2 VAD geometry in a part of the computational domain (inflow and outflow cannulas are much longer in the computational model). The VAD's casing (sketched as black lines) is the radial boundary of the computational domain).

Design or Assessment Parameters

The main assessment parameters for this study are:

Pressure increase via the VAD (pressure head) $H$

Hydraulic efficiency of the pump $\eta _{i}$

equivalent stress $\tau _{s}$ or effective stress $\tau _{eff}$

Modified index of hemolysis $MIH$

Volumes, in which certain stress thresholds $I_{\tau _{s}\geq threshold}$ are exceeded (volumentric analysis)

The rationale behind these parameters and details of how to calculate each of them are described below.

Pressure head

The pressure increase via a VAD is typically defined in millimeters of mercury $[mmHg]$ :

H={\frac {\sum _{outlet}\langle p_{tot}\rangle \cdot {\dot {m}}_{cell}}{\sum _{outlet}{\dot {m}}_{cell}}}-{\frac {\sum _{inlet}\langle p_{tot}\rangle \cdot {\dot {m}}_{cell}}{\sum _{inlet}{\dot {m}}_{cell}}}\qquad \qquad \qquad \qquad (1)

with $\langle p_{tot}\rangle$ as the time-averaged total pressure and ${\dot {m}}_{cell}$ as the mass flow in each cell. These variables are massflow-averaged by taking the sum over the cells at the outlet and inlet.

Inner efficiency

The inner efficiency balances the fluid power that the VAD realizes in the form of a pressure increase $H$ and fluid transport $Q$ , compared to the power required to make the VAD's impeller rotate with a certain torque $M$ and rotational speed $n$ .

\eta _{i}={\frac {H\cdot Q}{M\cdot 2\pi n}}\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad (2)

The blade torque $M$ is determined at the rotating surfaces of the impeller by accounting for the surface pressures $p_{sf}$ and surface shear stresses $\tau _{sf}$ (e.g. the impeller rotates around the z-axis, S denotes the impeller surface):

M_{z}=\langle {\big [}\int _{S}{r_{s}\times (n_{s}p_{sf}dS)}+\int _{S}{r_{s}\times (t_{s}\tau _{sf}dS)}{\big ]}\cdot e_{z}\rangle \qquad \qquad \qquad (3)

with:

$t_{s}$ - unit parallel vector
$n_{s}$ - unit vector parallel to $dS$
$e_{z}$ - unit vector parallel to the z-axis

Two types of inner efficiencies will be evaluated in section Evaluation. Firstly, the inner efficiency of the whole pump (index: $p$ ) and secondly, the inner efficiency of the impeller alone (index: $i$ ). Both types were determined respectively with the pressure build-up at the corresponding domain boundaries.

Equivalent and effective shear stress

The equivalent shear stress $\tau _{s}$ is the parameter, by which blood damage is predicted in blood simulations in complex flow geometries. The basis for this parameter derives from the shear stress tensor $\tau _{ij}$ , which contains the instantaneous spatial gradients of the velocity $u_{i}$ :

\tau _{ij}=\mu ({\frac {\partial u_{i}}{\partial x_{j}}}+{\frac {\partial u_{j}}{\partial x_{i}}})\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad (4)

For blood damage prediction in complex, three-dimensional flows, the stress tensor $\tau _{ij}$ is reduced to a scalar representation - the equivalent shear stress $\tau _{s}$ . Most numerical VAD simulations use a formulation based on the second invariant of the rate-of-strain tensor $S_{ij}$ :

\tau _{s}=\mu {\sqrt {2S_{ij}S_{ij}}}

, with

S_{ij}={\frac {1}{2}}({\frac {\partial u_{i}}{\partial x_{j}}}+{\frac {\partial u_{j}}{\partial x_{i}}})\qquad \qquad \qquad \qquad \qquad (5)

In URANS and LES methods, the instantaneous velocities are decomposed into a resolved component and an unresolved component. Accordingly, the instantaneous strain rate $S_{ij}$ is decomposed into a resolved variable ${\overline {S}}_{ij}$ and an unresolved component $S_{ij}^{'}$ :

S_{ij}={\overline {S}}_{ij}+S_{ij}^{'},\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad (6)

In URANS with unsteady mean flow, ${\overline {S}}_{ij}$ is an ensemble average (but can be considered as time-average, when the averaging time is long compared to the time scales of turbulent motions). In LES, ${\overline {S}}_{ij}$ is the instantaneous, resolved value (with the sub-grid scale fluctuations filtered out) and $S_{ij}^{'}$ are the subgrid-scale (SGS) fluctuations. This decomposition leads to the following expression for the equivalent stresses:

\tau _{s}=\mu {\sqrt {2({\overline {S}}_{ij}{\overline {S}}_{ij}+2{\overline {S}}_{ij}S_{ij}'+S_{ij}'S_{ij}')}}\qquad \qquad \qquad \qquad \qquad \qquad (7)

Taking the square of $\tau _{s}$ and applying the averaging (URANS) or filtering (LES) leads to an instantaneous, effective stress $\tau _{eff}$ , which, when introducing the kinematic viscosity $\nu =\mu /\rho$ , can be related to the dissipation of kinetic energy:

\tau _{eff}={\sqrt {\overline {\tau _{s}^{2}}}}=\mu {\sqrt {2{\overline {S}}_{ij}{\overline {S}}_{ij}+2{\overline {S_{ij}'S_{ij}'}}}}=\mu {\sqrt {{\frac {1}{\nu }}\underbrace {2\nu {\overline {S}}_{ij}{\overline {S}}_{ij}} _{\text{direct dissipation}}+{\frac {1}{\nu }}\underbrace {2\nu {\overline {S_{ij}'S_{ij}'}}} _{\text{turbulent dissipation}}}}(8)

Here, the first part is the direct viscous dissipation, which can be calculated from the resolved instantaneous strain rates ${\overline {S}}_{ij}$ . On the other hand, the second part is the turbulent dissipation $\varepsilon _{turb}$ as defined in equation (8) - in URANS it is the entire turbulent dissipation and in LES the dissipation by the sugrid-scale motion. This second part needs to be obtained from a model and is generally the dominant contributor to the dissipation in a fully turbulent flow.

For URANS, $\varepsilon _{turb}$ is given by the turbulence model (e.g., $\varepsilon =k\cdot \omega$ in the applied $k$ - $\omega$ -SST model). For our LES it is given by[6]:

\varepsilon _{turb}=\nu _{t,SGS}||{{\overline {S}}_{ij}}^{2}||,\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad (9)

where $\nu _{t,SGS}$ is the subgrid-scale eddy viscosity determined by the Smagorinsky SGS model.

As a result, the effective stress $\tau _{eff}$ , which is used for URANS and LES, is obtained from:

\tau _{eff}=\mu {\sqrt {2{\overline {S}}_{ij}{\overline {S}}_{ij}+(1/\nu )\varepsilon _{turb}}}\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad (10)

Please note that $\tau _{eff}$ in Eq. (10) is still an instantaneous value (as shown in Fig. 1.3). A time-averaged, effective stress $\langle \tau _{eff}\rangle$ is defined as the last step and this variable is used for the blood damage evaluation in section Evaluation and shown there in Fig. 5.2.

Modified hemolysis index

The modified index of hemoylsis $MIH$ is a widely used parameter to numerically assess hemolysis in VADs [3-6]. Hemolysis denotes the rupture of red blood cells. The ruptured cells release their hemoglobin into the blood plasma, which can lead to serious complications, like anemia or organ damage. A large number of numerical hemolysis prediction models are based on power laws. These models establish a direct relationship between hemolysis $MIH$ , an equivalent or effective stress $\tau _{eff}$ and the exposure time $t$ . Experimentally determined constants $(\alpha ;\beta ;C)$ are used to combine the variables to an MIH-value:

MIH=H_{l}\cdot 10^{6}=C\cdot \langle \tau _{eff}\rangle ^{\alpha }\cdot t^{\beta }\cdot 10^{6}.\qquad \qquad \qquad \qquad \qquad \qquad \qquad (11)

The hemolysis value $H_{l}(\tau _{eff},t)$ and the exposure time $t$ are needed to compute the MIH-value in the equation above. For numerical design and optimization studies, the MIH-value is tranformed into a formulation, which can be directly calculated from the computed, time-averaged flow field. This results in Eq. (11), which is also used in the present ERCOFTAC KB Wiki entry. The exact derivations of Eq. (10) to Eq. (11) are not given here, but can be found in Ref. [7].

{MIH}={\frac {1}{Q}}\int _{V}{(C^{1/\beta }\langle \tau _{eff}\rangle ^{\alpha /\beta }})dV)^{\beta }\cdot 10^{6}\qquad \qquad \qquad \qquad \qquad \qquad \qquad (12)

Volumetric analysis of stress thresholds

In this blood damage prediction model, which is widely used for design purposes, volumes $I_{\tau _{eff}\geq threshold}$ are calculated in the computational domain, which exceed certain stress thresholds. Three blood damage types can be analysed. For hemolysis, a stress threshold of $\langle \tau _{eff}\rangle \geq 150~Pa$ is defined. Additionally, this model addresses the degradation of von-Willebrandt proteins. Internal bleeding could happen, when this type of protein damage occurs. A stress threshold of $\langle \tau _{eff}\rangle \geq 9~Pa$ is given for the protein damage. Also, a stress threshold for the activation of thrombocytes is defined with $\langle \tau _{eff}\rangle \geq 50~Pa$ . The activation of thrombocytes could lead to a formation of a white thrombus, which can lead to a thromboembolic event or a pump thrombosis.

Flow Physics and Fluid Dynamics Data

The VAD was analysed at the nominal load ( $Q=4.5~l/min$ ) and partial load ( $Q=2.5~l/min$ ) operating point at a rotational speed of $n=7900~1/min$ . These conditions result in a Reynolds number of $Re_{p}=u_{2}\cdot D_{2}/\nu =3\cdot 10^{4}$ (at both operation points under investigation). Here, $u_{2}$ denotes the circumferential velocity of the blades near the casing (with a diameter of $D_{2}$ ). This pump Reynolds number $Re_{p}$ is comparable to the range of Reynolds numbers found in other blood pumps [10].

At the inlet of the VAD, a laminar inflow can be assumed, since the pipe Reynolds number based on the hydraulic diameter $D_{h}$ of the inlet cannula is small $Re_{cannula}=c_{cannula}\cdot D_{cannula}/\nu \approx 2300$ . Hence, all (coherent and turbulent) fluctuations and the associated stresses in the VAD's flow will be produced within the pump. The instantaneous velocities and instantaneous effective stresses in the rotor domain are shown in an animated Figure 1.3 in a cylindrical cut at 0.75 of the casing's radius $R_{2}$ . As displayed in that figure, an unsteady flow develops already at the inlet of the rotor channel. This is triggered by the fluctuating gap vortex in this inlet region, which comes as a result of the superposition of the main flow and a secondary flow from the rotor gap (region between rotor tip and casing). Increased effective stresses are observable in this gap vortex, which have a significant impact on the numerical blood damage prediction (analyzed in section Evaluation). Additionally, further secondary flows - which are typical for turbomachinery flows - such as the passage or horseshoe vortex, interact with each other in the VAD's blade channels, resulting in increased instanteous effective stresses. A detailed analysis of sources of effective stresses in the VAD can be found in Refs. [2] and [27].

Fig.1.3 Instantaneous flow in the rotor domain computed by LES (please click twice on this figure for an animated view). Visualised on a cylindrical cut at 0.75 of the casing's radius. Left: Velocity in the relative frame of reference. Right: Effective stresses (Eq. (10)). Please note that