Description AC7-03

Front Page

Description

Test Data

CFD Simulations

Evaluation

Best Practice Advice

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 turbomachinery, since they have a higher power density as pulsatile pumps, and therefore 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 needed pressure to sufficiently supply the blood flow in the circulatory system.

By using Computational Fluid Dynamics (CFD), VADs must be designed and optimised in such a way that they reproduce a physiological pressure increase in order to sufficiently supply the body with enough blood flow. Furthermore, they must be designed in order to guarantee that the blood, which passes 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 CFD simulation in a VAD can be challenging, since the inflow is laminar and all turbulence is produced within the pump and decays shortly after the pump outlet. Furthermore, the pump Reynolds number $Re_{p}={u_{2}D_{2}}/{\nu }$ is small with $Re_{p}=10^{3}-10^{4}$ compared to industrial pumps ($Re_{p}\approx 10^{6}$), and transition might occur.

In this respect, the aim of this study is to investigate the suitability of different URANS methods (with different turbulence models and solver settings) for the flow computation in an axial VAD. Here, both fluid mechanical parameters, such as the pump characteristics and velocity fields, as well as haemodynamic parameters, such as the haemolysis index MIH or stagnation zones, are investigated. The flow fields of the URANS simulations 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 a extended grid study.

Fig.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 desing point 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 AC 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}$

• Evaluation of stagnation areas based on the wall shear stress $\tau _{w}$

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}\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 $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}\langle p_{sf}\rangle dS)}+\int _{S}{r_{s}\times (t_{s}\langle \tau _{sf}\rangle 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}=0.5({\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 in the flow $\epsilon =2\nu S_{ij}S_{ij}$. From this, an effective shear stress can be derived from Eq.5, 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}={\sqrt {\rho \mu (\epsilon _{dir}+\epsilon _{turb}+\epsilon _{mod})}}\qquad \qquad \qquad \qquad \qquad \qquad \qquad (6)$

• 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 ruptures cells release their hemoglobin into the blood plasma, which can lead to serious complications, like anemia or organ damage for a VAD patient. 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$. Three experimentally determined constants $(\alpha ;\beta ;C)$ are used to link the parameters:

$MIH=H\cdot 10^{6}=C\tau _{eff}^{\alpha }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 suspension time $t$, Eq. 7 is transformed for a numerical implementation and 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 ${\underline {H^{1/\beta }}}$

$\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 ${\underline {MIH}}$:

${\underline {MIH}}={\underline {H^{\beta /\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 thresold 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.

• Evaluation of stagnation areas based on the wall shear stress

Stagnant, recirculating or detached flow are unfavourable in a VAD because at low velocities (and low wall shear stresses) the coagulation cascade is activated and a red thrombus can form. For this reason, it is useful for VAD simulations to analyse areas where small wall shear stresses $\tau _{w}\leq 1~Pa$ prevail in order to identify potential areas where a red thrombus could form.

Flow Domain Geometry

The flow field of an axial turbo pump was investigated. This design 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 as follows. 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 CFD, 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.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 bahiour, 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.2 VAD geometry in the computation domain.

Flow Physics and Fluid Dynamics Data

Fig.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 ($Q=4.5~l/min$ ) and partial load ($Q=2.5~l/min$ ) operation point at a rotational speed of $n=7900~r/min$. This results 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. 3 shows relevant regions within the rotor domain, where turbulence is produced due to secondary flow interactions, such as the:

• A = gap vortex
• E = turbulent boundary layer at the blades pressure side
• J = Interaction region of the gap vortex and the passage vortex
• L = cylindrical turbulence production region at the entry of the gap
• M1 = turbulent boundary layer at the shroud
• M2 = streak vortices in the turbulent boundary layer at the shroud
• N = turbulent boundary layer at the hub
• O = horseshoe vortex
• P = free-shear layer behind the rotor blades
• R = flow separation at the outlet guide vane

---

Contributed by: B. Torner — University of Rostock, Germany

Front Page

Description

Test Data

CFD Simulations

Evaluation

Best Practice Advice

© copyright ERCOFTAC 2021