Turbulent Blood Flow in a Ventricular Assist Device

Application Challenge AC7-03 © copyright ERCOFTAC 2021

CFD Simulations

Overview of CFD Simulations

Various large-eddy simulations (LES) and unsteady Reynolds-averaged Navier-Stokes (URANS) computations were carried out using the commercial flow solver ANSYS CFX. All simulations were performed at the nominal operation (design) point of the VAD. In total, five LES computations on different grid sizes were conducted for the verification of the LES results. The simulation at the finest grid was used as the reference case for the comparison with URANS. For the URANS cases, an extended grid convergence study was performed using seven URANS grids to analyze the influence of the spatial discretization on the main assessment parameters. Additionally, URANS computations with different turbulence models were performed on the finest grid for the comparison with the LES results. The used URANS turbulence models were: a $k$ - $\omega$ model, a $k$ - $\omega$ -SST model, and a $\omega$ -based Reynolds stress model.

Computational Domain

The whole VAD was considered in the numerical analysis. A sketch of the computational domain can be seen in Fig. X. Inflow and outflow cannulas were included in the computational domain. The inlet and outlet of the cannulas were placed sufficiently far away (four and seven impeller diameters respectively) from the pump in order to minimze the influences of the boundary conditions on the results.

Block-structured grids with hexahedral-elements were created using ANSYS ICEM CFD. Since URANS and LES have different requirements for grid resolution and quality, two different, final meshes were created: a.) The final mesh for the LES computation has a size of $102M$ elements. The grid was built according to literature recommendations by Fröhlich and Menter for wall-resolving LES methods. Attention was paid that the near-wall grid fitted the upper limits of $\Delta _{x}^{+}\leq 40$ for the grid with in flow direction and $\Delta _{z}^{+}\leq 15$ for the grid width in spanwise direction. Furthermore, the first wall-normal node had a maximal dimensionless wall distance of $y_{1}^{+}\leq 1$ and the grid growth factor was $r_{g}=1.05$ . This meshing strategy has been implemented throughout the whole domain. Grid angles were larger than 23 $^{\circ }$ and the volume change smaller than 5 with apect ratios smaller than 40 at the pumps wall. The aspect ratios were further reduced with increasing wall distance, so that the values range from 1 to 6 in the core flow region.

Fig.3.1 Computational grid of the LES. (a) Surface grid at the blades and the hub. (b) Volume grid in an axial cut-plane at the leading edge of the impeller blade.

The final URANS grid is coarser with a total size of $20M$ elements. Mesh quality criteria were kept within the ranges of the ANSYS CFX guidelines with maximum aspect ratios around 100, volume changes smaller than 6, and a minimum angle of $20^{\circ }$ . The mesh near the rotor wall has an area-averaged and maximal wall distance of $y_{1}^{+}=0.7$ and $2.1$ . All other pump parts have $y_{1}^{+}$ -values of one or smaller. Care was taken that the nearest wall layer contains more than 10 cells with a growth rate of $r_{g}=1.2$ . An advanced mesh convergence study of Eça and Hoekstra was performed to check whether the mesh size is appropriate to reflect the fluid mechanical and hemodynamical parameter. Therefore, the final URANS mesh was coarsened to 6 coarser grids with a mesh coarsening factor between 1.06 and 1.15.

Solution Strategy

In order to guarantee a properly working numerical scheme, several physical constraints were defined for the simulations. Blood was modeled as a Newtonian fluid with a dynamic viscosity of $\mu =3.5~mPas$ and a density of $\rho =1050~kg/m^{3}$ . This is a valid simplification because blood viscosity is independent of stress when shear rates are larger than $100~1/s$ and this condition is fulfilled in the pump.

The impeller rotated at a constant nominal frequency of 7,900 rpm. Transient rotor-stator coupling was used with a time step equal to 0.36° and 3° rotation of the impeller for LES and URANS, respectively. The resulting RMS Courant numbers were 0.6 for the Large-eddy and 3.4 for the Reynolds-averaged simulation. A total of 33 revolutions were calculated. Convergence was reached, when RMS residuals were at $10^{-5}$ and all monitored values were in a statistically steady condition. Time averaging was done for at least 13 revolutions for statistical analysis of the results.

Boundary Conditions

A constant, nominal flow rate of ܳFailed to parse (syntax error): {\displaystyle Q= 4.5 ݈l/min } was given at the outlet in both simulations. All walls were assumed to be hydraulically smooth and the no-slip condition applied. Total pressure was set to zero at the inlet. Since the Reynolds number is small at the inlet cannula (ܴ݁ $Re\approx 2,300$ ) and no perturbations are expected upstream, no turbulent inflow condition were specified.

Application of Physical Models

For LES, the governing equations were spatially and temporarly discretized by a bounded CDS scheme and a 2nd order Euler backward scheme. The dynamic Smagorinsky sub-grid scale model was used for closure of the LES equations with a bounded Smagorinsky parameter between $C_{s}=0.0-0.3$ .

For URANS, we used a high-resolution scheme and a 2nd order Euler backward scheme for spatial and temporal discretization, respectively. Three different turbulence models were applied as turbulence models for URANS: a $k$ - $\omega$ model, a $k$ - $\omega$ -SST model, and a $\omega$ -based Reynolds stress model. We expect that streamline curvature and transitional flow effects will influence the physics within the pump. Therefore, the curvature correction model and a Γ-Θ-transition model were applied for the $k$ - $\omega$ and the SST model.

Numerical Accuracy

Numerical accuracy studies were performed with various verification method for the computations, which has been already published in different papers of the author (see XXX). Two verification methods, one for LES and one for URANS, will be presented below.

Extended grid convergency study for URANS

The discretization error is a consequence of the approximation made in finite volume methods due to the transformation of the governing equations into a system of algebraic equations, which are discretely solved on a computational grid. An examination of the result’s discretization error is highly relevant and can be done for RANS methods by grid convergence studies. In this AC, the discretization error is estimated for URANS by means of an advanced grid convergence study, which determines an uncertainty interval as an indicator for the discretization error. A summary of the procedure is given below. For a detailed derivation, please see the study of Eça and Hoekstra. Generally, the discretization error can be estimated by:

\epsilon _{\phi }{\widetilde {=}}\phi _{i}-\phi _{0}=\alpha h_{i}^{p}\qquad \qquad \qquad \qquad (11)

$\phi _{i}$ is an arbritary flow quantity, $\phi _{0}$ is the estimate of the exact result, $\alpha$ is a constant, Failed to parse (unknown function "\inc"): {\displaystyle h_i=(\inc_x \inc_y \inc_z)^{1/3} } is the cell size, and $p$ is the observed order of convergence. Values of $\alpha$ , $p$ and $\phi _{0}$ are needed for the error estimation. These values can be determined by using at least four simulations on different grid sizes and a least-squares error estimation using the minimum of the following equation:

S(\phi _{0},\alpha ,p)={\sqrt {\sum _{i=1}^{N}(\phi _{i}-(\phi _{0}+\alpha h_{i}^{p}))^{2}}}\qquad \qquad \qquad \qquad (12)

The resulting parameter $(\phi _{0}+\alpha h_{i}^{p})$ characterizes a polynomial function, which is called the fit. This fit approximates the development of the computed CFD results. The deviation between the CFD result $\phi _{i}$ and the fit is specified by the fit’s standard deviation $\sigma$ .

The discretization error of the CFD result is indicated by an uncertainty $U_{\phi }$ , which is a 95% confidence interval, wherein the exact solution of the URANS equations is expected. The determination of $U_{\phi }$ depends on several fit and convergence properties, which are explained in detail in the study of Eça and Hoekstra. An example of a “good” error estimation is given by Eq. (13), which uses the fit’s data and an additional safety factor $F_{s}$ .

The resulting uncertainty interval $U_{\phi }$ is an indicator for the deviation between the computed result $\phi _{i}$ on the (typical) cell size $h_{i}$ and the estimated, exact URANS solution $\phi _{0}$ . Furthermore, it is possible to examine whether a computational grid is sufficiently fine or still too coarse to indicate a grid-independent and reliable solution for the analyzed parameter of interest.