Turbulent Blood Flow in a Ventricular Assist Device

Application Challenge AC7-03 © copyright ERCOFTAC 2021

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 such that they produce a physiological pressure increase sufficient to supply the body with enough blood flow. Furthermore, they must be designed such that the blood passing the VAD is not damaged due to non-physiological flow conditions (high shear stresses, stagnation areas, high turbulent kinetic energy (TKE) regions, ...) 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 turbomachinery, 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 methods 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 haemolysis index MIH, are investigated. The stress fields of the URANS' will be compared with a highly turbulence-resolving large-eddy simulations, 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 the one side, to guarantee that the VAD offers an physiological relevant pressure increase at the chosen operation points to sufficiently support the VAD patient. On the other side, haemodynamical parameter must be evaluated in these studies. Here, it is important that the CFD reflects relevant regions for potential blood damage or thrombi formation, so that these regions can be minimised in the optimization procedure. Additionally, CFD is important for VAD studies in order to compare different designs to find the pump with the highest haemocompatibility (lowest blood damage).

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

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 (scalar) shear stress $\tau _{s}$

Modified index of hemolysis $MIH$

Volumetric analysis of stress thresholds $I_{\tau _{s}\geq threshold}$

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

Pressure head

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

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

with $p_{tot}$ as the time-averaged total pressure, which is massflow-averaged at the outlet and inlet.

Hydraulic efficiency

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

$Q$ denotes the flow rate and $n$ the rotational speed. The blade torque $M$ is determined at the rotating surfaces of the impeller by accounting 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}={\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}\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

Equivalent shear stress

The equivalent shear stress $\tau _{s}$ is the parameter, by which numerical blood damage is assessed in a VADs flow simulation. The basis for this parameter derives from the shear stress tensor $\tau _{ij}$ , which is built with the velocity gradients in a flow field:

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

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

\tau _{s}=\mu {\sqrt {2\cdot S_{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 \qquad \qquad (5)

The expression within the square-root in Eq. (5) can be further connected to a flow variable, the energy dissipation rate $\epsilon =2\nu S_{ij}S_{ij}$ . From this, an effective shear stress can be derived, which connects the effectiv stress $\tau _{eff}$ with the computed dissipation rates in the time-averaged flow, namely the direct dissipation $\epsilon _{dir}=2\nu \langle S_{ij}\rangle \langle S_{ij}\rangle$ , the resolved turbulent dissipation $\epsilon _{turb}=2\nu \langle S_{ij}^{'}S_{ij}^{'}\rangle$ and the modelled turbulent dissipation from the turbulence model $\epsilon _{mod}=2\langle \nu _{t}S_{ij}S_{ij}\rangle$ :

\tau _{eff,1}={\sqrt {\rho \mu (\epsilon _{dir}+\epsilon _{turb}+\epsilon _{mod})}}\qquad \rightarrow (w/~\epsilon _{mod})\qquad \qquad \qquad \qquad \qquad \qquad \qquad (6.1)

The formulation in Eq. (6.1) has been used in recent VAD flow studies. In earlier VAD flow studies, also an equivalent stress (Eq. (6.2)) was used, which ommit any component from the turbulence model $\epsilon _{mod}$ . In this article, the influence of this turbulent term on the blood damage prediction will be discussed and hence, both equivalent stresses will be considered.

\tau _{eff,2}={\sqrt {\rho \mu (\epsilon _{dir}+\epsilon _{turb}+0)}}\qquad \rightarrow (w/o~\epsilon _{mod})\qquad \qquad \qquad \qquad \qquad \qquad \qquad (6.2)

Modified hemolysis index

The modified index of hemoylsis $MIH$ is a widely used parameter to numerically assess hemolysis in VADs. 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 haemolysis prediction models are based on power laws. These models establish a direct relationship between hemolysis $MIH$ , an equivalent 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\cdot 10^{6}=C\cdot \tau _{eff}^{\alpha }\cdot t^{\beta }\cdot 10^{6}\qquad \qquad \qquad \qquad \qquad \qquad \qquad (7)

Because of the non-linear relationship between the hemolysis value $H(\tau _{eff};t)$ and time $t$ , Eq. (7) is transformed for a numerical implementation and then written as a transport equation:

{\frac {\partial H^{1/\beta }}{\partial t}}+u_{i}{\frac {\partial H^{1/\beta }}{\partial x_{j}}}=C^{1/\beta }\tau _{eff}^{\alpha /\beta }\qquad \qquad \qquad \qquad \qquad \qquad \qquad (8)

If it is now assumed that the hemolysis is time-independent and homogeneous in space, Eq. (8) can be further simplified to a (temporally and spatially) averaged hemolysis value:

\int _{V}{(u_{i}{\frac {\partial H^{1/\beta }}{\partial x_{j}}}})dV=\int _{V}{(u_{i}\cdot n_{i}\cdot \partial H^{1/\beta }})dS=\int _{V}{(C^{1/\beta }\tau _{eff}^{\alpha /\beta }})dV\qquad \qquad \qquad \qquad \qquad \qquad \qquad (9)

Assuming that the hemolysis is zero at the VADs inlet, Eq. (9) can be further simplified to geht an averaged MIH-value ${MIH}$ :

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

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 $\tau _{eff}\geq 150~Pa$ is defined. Additionally, this model address the degradation of von-Willebrandt proteins. Internal bleeding could happen, when this type of protein damage occur. A stress threshold of $\tau _{eff}\geq 9~Pa$ is given for the protein damage. Also, a stress threshold for the activation of thrombocytes is defined with $\tau _{eff}\geq 50~Pa$ . The activation of thrombocytes could lead to a formation of a white thrombus, by which a thromboembolic event or a pump thrombosis could happen.

Flow Domain Geometry

The flow field of an axial turbo pump is investigated. The 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 ( $11~mm$ ) and outer ( $16~mm$ ) diameter, meridian lines were placed between hub an 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 ( $120~\mu m$ ) were included. The outlet guide vanes were set 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 operation point.

The axial VAD within the computational domain is illustrated in Fig. 1.2. It consist of a two-bladed rotor, an inlet guide vane with five, slightly bended blades, to apply a counter-swirl, which theoretically leads to an increased pressure head for a specific rotational speed, and a three-bladed outlet guide vane.

This VAD was designed for reseach purposes solely. Hence, a bearing concept for a clinical use was not considered. Furthermore, no axial gaps between the rotor and stators were included in the original design. 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 the computation domain.

Flow Physics and Fluid Dynamics Data

Fig.1.3 Secondary flow interactions within the VAD, which lead to an increase in production of turbulent kinetic energy.

The VAD was analysed at the nominal load ( $Q=4.5~l/min$ ) and partial load ( $Q=2.5~l/min$ ) operation point at a rotational speed of $n=7900~r/min$ . These condition result in a pump Reynolds number of $Re_{p}={\frac {u_{2}\cdot D_{2}}{\nu }}=3\cdot 10^{4}$ . Here, $u_{2}$ denotes the circumferential velocity at the blades outer diameter $D_{2}$ .

The blood was treated as Newtonian, single-phase fluid with a density ( $\rho =1050~kg/m^{3}$ ) and dynamic viscosity ( $\mu =3.5~mPas$ equal to human blood at a hematocrit of $H_{ct}=45~\%$ . Assuming a Newtonian fluid behaviour in a VAD is reasonable, since blood shows an asymptotic viscosity at shear rates greater than $S_{ij}\geq 100~1/s$ , which is fulfilled in nearly all parts of the VAD.

At the inlet of the VAD, a laminar inflow can be assumed, since the pipe Reynolds number based on the hydraulic diameter is small $Re_{pipe}={\frac {u_{bulk}\cdot D_{h}}{\nu }}\approx 2300$ . Hence, all turbulence will be produced within the VAD. This is realised due to secondary flow interactions within the VAD, which lead to bypass transition and to an increase in turbulence kinetic energy. Fig. 1.3 shows relevant regions within the VAD at nominal operation point, where turbulence is produced due to secondary flow interactions, such as the: