CFD Simulations 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

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 and a partial operation 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, an URANS computation with a $k$ - $\omega$ - SST turbulence models were performed on the finest grid for the comparison with the LES results.

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 (SVG with PNG fallback (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\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

Various numerical accuracy studies were performed with different verification methods for the URANS and LES computations. These studies 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}$.

$U_{\phi }(\phi _{i})=F_{s}\cdot \mid \alpha h_{i}^{p}\mid +\sigma +\mid \phi _{i}-\phi _{fit}\mid \qquad \qquad \qquad \qquad (13)$

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.

In total, seven meshes were used for the advanced grid convergence study. These meshes are summarised in Tab. 3.1. Since the mesh influence on the computation of, e.g., stresses is one major aspect of this AC, the results of the grid convergency study will be presented in the evaluation Evalution.

Table 3.1 Names and grid size properties of the seven URANS grids.

Grid	Million grid elements $N$	Typical cell size (m)	Grid refinement parameter	Time and space averaged $y_{1}^{+}$-value over blades and maximum
UR-1	3.3	$3\cdot 2.11^{-4}$	1.06	1.43(3.58)
UR-2	4.0	$3\cdot 1.98^{-4}$	1.15	1.25(3.29)
UR-3	5.9	$3\cdot 1.73^{-4}$	1.06	1.12(3.00)
UR-4	7.0	$3\cdot 1.63^{-4}$	1.11	1.00(2.74)
UR-5	9.6	$3\cdot 1.47^{-4}$	1.12	0.89(2.49)
UR-6	13.6	$3\cdot 1.31^{-4}$	1.12	0.79(2.31)
UR-7	19.4	$3\cdot 1.16^{-4}$	-	0.68(2.14)

Power-Loss-Analysis (PLA) for LES

The power loss analysis is an LES quality assessment method, which was developed by the authors of this study. It can be used to investigate whether the performed simulation is adequately capable of computing the internal power losses due to TKE production and dissipation within the flow. These processes, espacially the dissipation, are connected with the computed stresses (see Eq. (6) in the Description). The method is graphically shown in Figure 3.2. In the PLA, three definitions of the total power losses $P_{tot}$ for the mean (time-averaged) flow field are calculated and afterwards, the results from the three total power losses are compared. The rationale and derivations of the PLA are comprehensively explained in Ref. X. In the following, just the most important informations regarding the PLA will be explained.

Fig.3.2 Sketch of the PLA method.

The first total power loss $P_{tot,1}$ is determined by using the pump characteristics of the VAD. The power loss can be calculated by taking the difference between the power at the coupling $P$ (blade's torque $M_{bl}$ multiplied with angular velocity $2\pi n$) and the increase in hydraulic power $P_{hyd}$ via the VAD, see Equation (14):

$P_{tot,1}=P-P_{hyd}=M_{bl}\cdot 2\pi n-\Delta p_{tot}\cdot Q.\qquad \qquad \qquad \qquad (14)$

The total power loss can also be determined by integrating and summing up all internal power losses within the VAD. Three loss parts must be considered as internal losses for LES. The first internal loss $P_{dir}$ (Equation (15)) is due to direct dissipation $\varepsilon _{dir}$ and describes the transfer of mean flow energy directly into heat.

$P_{dir}=\rho \int _{V}\epsilon _{dir}dV=\mu \int _{V}\langle {\frac {\partial u_{i}}{\partial x_{j}}}\rangle {\bigg (}\langle {\frac {\partial u_{i}}{\partial x_{j}}}\rangle +\langle {\frac {\partial u_{j}}{\partial x_{i}}}\rangle {\bigg )}dV\qquad \qquad \qquad (15)$

The second loss part is the energy transfer from mean flow into turbulent motion and is called turbulence production $TKE_{prod}$. The corresponding loss term $P_{prod}$ is defined in Eq. (16). It is taken from the directly resolved flow field of a flow computation.

$P_{prod}=\rho \int _{V}TKE_{Prod}dV=\rho \int _{V}-\langle u_{i}'u_{j}'\rangle {\bigg \langle }{\frac {\partial u_{i}}{\partial x_{j}}}{\bigg \rangle }dV\qquad \qquad \qquad (16)$

Since the LES does not resolve all of the turbulent motions and losses, the modeled loss contribution from the smallest scales must be included into the PLA. In general, the LES turbulence model must ensure that a proper amount of turbulent dissipation from the smallest scales is provided. Due to the equilibrium between turbulence production and dissipation at the smallest scales, the modeled turbulent losses from the applied Smagorinsky model can be written as in Eq. (17):

$P_{mod}=\rho \int _{V}\epsilon _{mod}dV=\rho \int _{V}\langle \nu _{t}{\frac {\partial u_{i}}{\partial x_{j}}}{\bigg (}{\frac {\partial u_{i}}{\partial x_{j}}}+{\frac {\partial u_{j}}{\partial x_{i}}}{\bigg )}\rangle dV.\qquad \qquad \qquad (17)$

Hence, the second total power loss $P_{tot,2}$ is stated in Eq. (18) by summing up the loss parts from Eq. (15), (16) and (17).

$P_{tot,2}=P_{dir}+P_{prod}+P_{mod}\qquad \qquad \qquad \qquad (18)$

The losses from the resolved turbulent flow part in Eq. (16) can also be expressed by another flow variable, which acts on different turbulent scales. It is the resolved turbulent dissipation rate $\epsilon _{turb}$. This variable expresses the energy loss from the small turbulent motions into heat and is defined in terms of a power loss $P_{diss}$ in Eq. (19):

$P_{diss}=\rho \int _{V}\epsilon _{turb}dV=\mu \int _{V}\langle {\frac {\partial u'_{i}}{\partial x_{j}}}{\bigg (}{\frac {\partial u'_{i}}{\partial x_{j}}}+{\frac {\partial u'_{j}}{\partial x_{i}}}{\bigg )}\rangle dV.\qquad \qquad \qquad \qquad (19)$

The summation of Eq. (15), (17) and (19) leads to a third formulation of a total power loss $P_{tot,3}$ for a turbulent flow field, see Eq. (20):

$P_{tot,3}=P_{dir}+P_{diss}+P_{mod}.\qquad \qquad \qquad \qquad (20)$

The PLA method makes it possible to compare the direct losses in the mean flow (due to $\epsilon _{dir}$) with the turbulent losses due to the turbulent energy cascade ($TKE_{prod}$ and $\epsilon _{turb}$). From a numerical point of view, the PLA can be used as a quality assessment method, since the total power losses from Eq. (14), (18) and (20) should be identical per definition. If this is fulfilled, it can be verified that the power losses in form of direct and turbulent dissipation rates as well as turbulence production were adequately computed by a simulation.

Table 3.2 Results of the PLA analysis for the finest LES grid.

PLA variable & operator	Value
$P_{bl}~-$	$2.231~W$
$P_{hyd}~=$	$0.739~W$
${P_{tot,1}}$	${1.492~W}$
Failed to parse (unknown function "\bm"): {\displaystyle \bm{(P_{tot,2}-P_{tot,1})/P_{tot,1}\cdot 100} }	Failed to parse (unknown function "\bm"): {\displaystyle \bm{-1.5 \%} }
$P_{dir}~+$	$1.159~W~(78.8\%)$
$P_{prod}~+$	$0.286~W~(19.5\%)$
$P_{mod}~=$	$0.025~W~(1.7\%)$
${P_{tot,2}}$	${1.470~W}~(100\%)$
Failed to parse (unknown function "\bm"): {\displaystyle \bm{(P_{tot,2}-P_{tot,1})/P_{tot,1}\cdot 100} }	Failed to parse (SVG with PNG fallback (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \bm{-1.5 \%} }
${(P_{tot,2}-P_{tot,1})/P_{tot,1}\cdot 100}$	${-1.5\%}$
Failed to parse (unknown function "\bm"): {\displaystyle \bm{(P_{tot,2}-P_{tot,1})/P_{tot,1}\cdot 100} }	Failed to parse (SVG with PNG fallback (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \bm{-1.5 \%} }
$P_{dir}~+$	$1.159~W~(82.2\%)$
$P_{diss}~+$	$0.225~W~(16.0\%)$
$P_{mod}~=$	$0.025~W~(1.8\%)$
${P_{tot,3}}$	${1.409~W}~(100\%)$
Failed to parse (unknown function "\bm"): {\displaystyle \bm{(P_{tot,2}-P_{tot,1})/P_{tot,1}\cdot 100} }	Failed to parse (unknown function "\bm"): {\displaystyle \bm{-1.5 \%} }
${(P_{tot,3}-P_{tot,1})/P_{tot,1}\cdot 100}$	${-5.6\%}$

Table 3.2 summarizes the results of the PLA for the turbulent production rates $TKE_{prod}$ as well as for the turbulent dissipation rates $\epsilon _{turb}$ for the finest LES grid. From a numerical point of view, the relative deviation between the total power losses is smaller using the turbulence production $P_{tot,2}$ ($-1.5\%$) as with the turbulent dissipation $P_{tot,3}$ ($-5.6\%$). Furthermore, both loss definitions shows that the majority of the turbulent losses is directly resolved within the flow $(16.0\%-19.5\%$) instead of being modeled $1.7\%-1.8\%$).

The high dependency on the spatial discretization can be analyzed using Fig. 3.3. The subfigures show some PLA results, which were ascertained on different grid sizes. It can be seen from subfigure~a.) that the relative deviations between the total power losses decrease with increasing grid size to $-1.5\%$ and $-5.6\%$, respectively. Furthermore, it is noticeable from subfigure~b.) that the total power loss $P_{tot,1}$ (based on the pump characteristics, Eq. (14)) is relatively constant from a grid size of 10M nodes with a value of $P_{tot,1}\approx 1.50~W$. In addition, the internal loss part $P_{dir}$ from direct dissipation $\epsilon _{dir}$ shows a nearly constant progression with a value of $P_{dir}\approx 1.16~W$ between 30M and 102M node points. Therefore, the reduction between the power losses $P_{tot,2}$, respectively $P_{tot,3}$, and $P_{tot,1}$ in subfigure a.) must be a result of a more accurate computation of the turbulent losses based on the turbulent dissipation and turbulence production rates.

Despite the fact that the finest grid cannot resolve the turbulent dissipation rates in total, the relative deviation to the reference value ($P_{tot,1}$) is relatively small with $\approx 5\%$. Contrariwise, that means that $95\%$ of the dissipative losses can be investigated within the flow field of the simulation. In addition, the simulation on the 102M grid is also able to reproduce nearly all losses from the turbulence production $P_{prod}$. Together with the fact that the majority of the losses are directly resolved instead of being modeled, it can be concluded that the 102M grid simulation is appropriate for the flow investigation in the VAD.

Fig.3.3 PLA results. (A) Relative deviations between $P_{tot,2}$ resp. $P_{tot,3}$ and $P_{tot,1}$. (B) Total power loss $P_{tot,1}$ on different grid sizes. Also (B), loss part due to direct dissipation on different grid sizes. (C) Loss part due to turbulence production rates and dissipation rates on different gird sizes.

---

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

Front Page

Description

Test Data

CFD Simulations

Evaluation

Best Practice Advice

© copyright ERCOFTAC 2021