CFD Simulations AC7-04

Front Page

Description

Test Data

CFD Simulations

Evaluation

Best Practice Advice

A pulsatile 3D flow relevant to thoracic hemodynamics: CFD - 4D Flow MRI comparison

Application Challenge AC7-04   © copyright ERCOFTAC 2021

CFD Simulations

Overview of CFD Simulation

Large Eddy Simulations were carried out using the in-house, massively parallel and multiphysics YALES2BIO solver based on YALES2 [4] developed at CORIA (Rouen, France). YALES2BIO is dedicated to the simulation of blood flows at the macroscopic and microscopic scales. The base is a solver for the incompressible Navier-Stokes equations. The equations are discretised using a finite-volume fourth-order scheme, adapted to unstructured meshes [5,6]. The divergence-free property of the velocity field is ensured thanks to the projection method introduced by Chorin [7]. The velocity field is first advanced in time using a low-storage fourth-order Runge-Kutta scheme [6,8] in a prediction step. This predicted field is then corrected by a pressure gradient, obtained by solving a Poisson equation to calculate pressure. This equation is solved with the Deflated Preconditioned Conjugate Gradient algorithm [9]. YALES2BIO was validated and successfully used in many configurations relevant to cardiovascular biomechanics (see [10] for a list of publications). The boundary conditions applied at the inlet came from the data acquired during the experiment (2D cine PC-MRI).

Simulation Case

Computational Domain

The geometry used in the calculation is the same as the one presented in the test data section. The computational domain has one inlet and one single outlet, so that the flow split at the main tube/collateral junction does not depend on the details of the outlet boundary condition. The mesh (Fig. 5) used to perform the simulation consisted of an unstructured tetrahedral mesh generated with GAMBIT 2.4.6 (ANSYS, Inc., Canonsburg, PA). The mesh includes 3,812,438 cells with an average cell volume of 38 mm³ (representative cell size of 0.7 mm).

Figure 5: Computational domain and numerical mesh

Left: Geometry - Dimensions are given with respect to the center of the collateral, Right: Mesh at the inlet plane

Solution Strategy

The blood, or blood-mimicking fluid, is described by the filtered Navier-Stokes equations (NSE):

$\left\{{\begin{array}{lcl}{\displaystyle {\frac {\partial {\overline {\left[\mathbf {u} \right]}}}{\partial t}}+\nabla \cdot \left({\overline {\left[\mathbf {u} \right]}}\otimes {\overline {\left[\mathbf {u} \right]}}\right)=-{\frac {\nabla {\overline {\left[p\right]}}}{\rho }}+\left(\nu +\nu _{SGS}\right)\Delta {\overline {\left[\mathbf {u} \right]}}}\\{\displaystyle \nabla \cdot {\overline {\left[\mathbf {u} \right]}}=0}\end{array}}\right.$

where ${\overline {\left[.\right]}},\mathbf {u} ,p,\rho =1020{\text{ kg/m}}^{3},\nu =4.02\times 10^{-6}{\text{m}}^{2}{\text{/s}}$ and $\nu _{SGS}$ are the filter operator, velocity field, pressure field, density, kinematic viscosity and subgrid-scale kinematic viscosity of the fluid, respectively.

The governing equations are discretized and solved using the YALES2BIO solver. In the present simulation a centred fourth-order numerical scheme with an explicit fourth-order Runge-Kutta time advancement scheme was used. To ensure numerical stability, the time step was computed to keep a CFL number equal to 0.9. The pressure advancement and divergence-free condition were met thanks to a fractional-step method (Kim & Moin, 1985 [11]), a modified version of the Chorin’s algorithm (Chorin, 1968 [7]). The Sigma eddy-viscosity-based LES model (Nicoud et al., 2011 [12]) was employed to take into account the turbulence effects. In this model the subgrid-scale kinematic viscosity is modelled as:

$\nu _{SGS}=\left(C\Delta \right)^{2}{\frac {\sigma _{3}\left(\sigma _{1}-\sigma _{2}\right)\left(\sigma _{2}-\sigma _{3}\right)}{\sigma _{1}^{2}}}$

where $C=1.35$ is the model constant, $\Delta$ is the typical size of the local cell of the mesh and $\sigma _{1}\geq \sigma _{2}\geq \sigma _{3}\geq 0$ are the three singular values of the local velocity gradient tensor [12]. The Sigma model guarantees that no eddy viscosity is applied in canonical laminar flows and is well adapted to wall-bounded flows and transitional flows [13,14].

Boundary Conditions

Zero pressure condition was prescribed at the outlet section while no-slip condition was imposed at the solid boundaries. Over the whole cycle, an averaged $y^{+}=2.5$ was reported. The highest values were found at peak systole with a mean $y^{+}=3.1$ and maximal $y^{+}=8.3$.

To prescribe the velocity field at the inlet, a pixel-based inflow was derived from the 2D cine PC-MRI measurements. Beforehand the velocity field from these measurements were corrected to account for partial volume effects caused by the phantom walls and random noise, in particular for the through-plane velocity component $w$. Thereby, a linear $w$-distribution was assumed between $0$ at the wall and $w(R_{lim})$ at $R_{lim}$, where $R_{lim}$ is an user-defined inner radius, here chosen to be $95\%$ of the radius of the circular cross-sectional inlet surface. Then a 2D bilinear interpolation from this corrected velocity field onto the inlet surface of the CFD domain was performed for each of the 30 cardiac phases. Once the velocity field had been spatially interpolated, its discrete temporal evolution had to be converted into a time-continuous signal. As the flow rate delivered by the pump was periodic, a trigonometric interpolation was used by applying at each node a discrete Fourier transform of the temporally discretized velocity signal. Thereby 15 non-redundant complex Fourier coefficients were extracted and the continuous inlet velocity signal could be achieved by injected back these coefficients into a Fourier series expansion for each node of the inlet surface.

CFD post-processing

In order to enable comparison with the 4D Flow MRI measurements, two main post-processing steps were implemented to be coherent with the MRI acquisition process.

Phase-averaging

First, the CFD was phase-averaged. Indeed during a conventional 4D Flow scan, the k-space (spatial frequencies space) is sequentially filled over hundreds of RF-pulses and magnetic gradients sequences. Thus it does not appear relevant to compare instantaneous CFD velocity fields with MRI measurements. Another argument is that high Reynolds number unsteady flows may lead to cycle-to-cycle fluctuations. Hence an appropriate way seems to be to phase-average the simulated flow velocity field over a sufficiently representative number of cardiac cycles. The phase-averaged velocity ${\bar {\mathbf {u} }}(\mathbf {x} ,t)$ at spatial coordinates $\mathbf {x} =(x,y,z)$ and at time $t$ was defined as:

${\bar {\mathbf {u} }}(\mathbf {x} ,t)={\frac {1}{N}}\sum _{k=0}^{N-1}\mathbf {u} (\mathbf {x} ,t+kT)$

where $\mathbf {u} (\mathbf {x} ,t)$ refers to the instantaneous velocity vector and $N=30$ is the total number of cycles of period $T=1$ s. As will be further expanded on in the numerical accuracy section, the first 10 simulated cycles were removed to cancel the effect of the arbitrary initial condition, and thus $N$ corresponds to the number of cycles computed afterwards.

Downsampling

Then, the results of both methods have to be given on the same grid in order to allow a quantitative comparison. One possibility is to interpolate the MRI velocity fields on the finer CFD grid. Another method, which is the one chosen here, is to downsample the CFD velocity fields onto the coarser MRI grid. The idea of the downsampling method is to localize the CFD nodes within a given MRI voxel and to average the contributions of each CFD node corresponding to this MRI voxel in order to build a low resolution CFD (LR-CFD) velocity field. The steps of the downsampling process (Fig. 6) were the following: first the MRI grid was subsampled with 5 voxels in each direction (5$^{3}$ = 125 subvoxels per voxel). Then the CFD solution was linearly interpolated on the subsampled Cartesian grid. Finally these interpolated velocity fields were averaged over the 125 subvoxels included in each MRI voxel, which lead to the LR-CFD velocities expressed on the MRI grid.

Figure 6: Illustration of the 3D downsampling process as viewed from the coronal mid (xz)-plane. The MRI grid is first subsampled with 5 voxels in each direction (125 subvoxels per voxel). The CFD solution is linearly interpolated on the subsampled grid, which is downsampled back to the MRI grid by averaging the velocity over the 125 subvoxels included in each MRI voxel.

Numerical Accuracy

The numerical accuracy studies were partly conducted on previously obtained data and are detailed in T. Puiseux's PhD thesis [1]. The CFD simulation process was the same as the one described above. The only difference to the data exposed in this document is the shape of the inlet velocity profile. Note however that the peak systolic Reynolds numbers at the inlet are very similar (1830 in this previous study against 1915 here). This allows us to have confidence in these previous analyses, which were conducted for mesh convergence and phase-averaging sensitivity. Turbulence modelling employed for the subgrid scales was analysed for the present data.

Mesh sensitivity analysis

To ensure spatial convergence of the simulations, the effects of the mesh were studied on 3 different meshes: a coarse one with 622 thousand cells (cell size = 1.3 mm), a medium one with 1 284 thousand cells (cell size = 1.0 mm) and a fine one with 3 812 thousand cells (cell size = 0.7 mm). The magnitude of the phase-averaged velocity in the aneurysm resulted in a numerical uncertainty of 1.88% for the latter mesh, with an apparent spatial order of convergence p = 2.0 (see Sigüenza et al., 2016 [15] for further details on the method). Therefore the fine mesh was considered fine enough for the velocity field to be spatially converged and independent of the spatial resolution.

Phase-averaging sensitivity analysis

A sensitivity analysis was carried out to estimate the minimal number of cycles needed for the phase-averaged velocities to converge. First, as it was mentioned in the paragraph about the phase-averaging process, the first 10 cycles were removed to cancel the impact of the initial condition. 60 cycles were then simulated and the 60-cycle phase-averaged velocity was compared to the phase-averaged velocity had the simulation been simulated for $n$ cycles, with $n<60$. The results are presented in Fig. 7, where the metric used for analysis is the Pearson's correlation coefficient $r^{2}$. The $r^{2}$-coefficient measures the linear correlation between two sets of data and is defined as:

$r_{ab}^{2}={\frac {(\sum _{i=1}^{N}(a_{i}-{\overline {a}})(b_{i}-{\overline {b}}))^{2}}{\sum _{i=1}^{N}(a_{i}-{\overline {a}})^{2}\sum _{i=1}^{N}(b_{i}-{\overline {b}})^{2}}}$

where $a$ and $b$ are two datasets, $N$ the sample size, the subscript $i$ the index of the individual sample points and ${\overline {\cdot }}$ the sample mean.

The results show that whereas the main components of the velocity ($u$ and $w$) were well-converged after about 5 cycles, around 30 cycles were required for the transverse component of the velocity ($v$). Thereby, 30 cycles were considered sufficient to reach a reasonable asymptotic value for the velocity field. Furthermore, the reproducibility of the process was tested by extracting and independently phase-averaging two subsequent set of 30 cycles from the 60 cycles simulated. High correlation levels of the Pearson’s correlation coefficient $r^{2}$ were found $(r^{2}(u)=0.995,\ r^{2}(v)=0.930,\ r^{2}(w)=0.999{\text{ and }}r^{2}(\vert \vert \mathbf {u} \vert \vert )=0.997)$, which confirms that the process was independent from the set of cycles selected.

Figure 7: Evolution of the Pearson’s correlations ($r^{2}$) between phase-averaged CFD velocity computed over different numbers of cycles and the 60-cycle phase-averaged CFD velocity. All nodes of the computation domain were included in the metric.

Turbulence resolution analysis

In the LES strategy the subgrid scales are not explicitly resolved but modelled. A general criterion of the turbulence resolution in LES was proposed by Pope (2000, [16]), where LES is considered well-resolved if at least 80% of the turbulent kinetic energy is solved. Formally, this can be written as:

$M({\textbf {x}},t)={\frac {k_{sgs}({\textbf {x}},t)}{k_{total}({\textbf {x}},t)}}={\frac {k_{sgs}({\textbf {x}},t)}{k_{sgs}({\textbf {x}},t)+k_{res}({\textbf {x}},t)}}<0.2$

where $k_{res}({\textbf {x}},t)$ is the turbulent kinetic energy of the resolved eddies and $k_{sgs}({\textbf {x}},t)$ the turbulent kinetic energy of the unresolved subgrid scales. This latter quantity is estimated since it is not explicitly calculated in LES. The following estimate was proposed by Yoshizawa and Horitu (1985, [17]):

$k_{sgs}({\textbf {x}},t)={\frac {\nu _{sgs}^{2}}{(C\Delta )^{2}}}$

where $\nu _{sgs}$ is the subgrid scale viscosity added by the model, $\Delta$ the characteristic mesh size and $C\approx 0.1$ the dimensionless dynamic constant of the model.

As mentioned in the CFD post-processing section, the velocity of interest in this AC is not the instantaneous one, but a phase-averaging over 30 cycles after removal of the first 10 initializing cycles. The fluctuating part of the velocity was defined as:

$\mathbf {u} '(x,t)=\mathbf {u} (x,t)-{\overline {\mathbf {u} }}(x,t)$

and the resolved turbulent kinetic energy was computed as:

$k_{res}(x,t)={\frac {1}{2}}({\overline {u'^{2}}}(x,t)+{\overline {v'^{2}}}(x,t)+{\overline {w'^{2}}}(x,t))$

where ${\overline {\ \cdot \ }}$ stands for the phase-averaging. Likewise, a phase-averaged subgrid scale viscosity ${\overline {\nu _{sgs}}}$ was used to compute $k_{sgs}$.

The resolution of the turbulence was studied for two meshes using the MR data acquisitions presented in the 'Test data' section. The first simulation was performed with the fine mesh previously mentioned of 3.8 million cells (characteristic cell size of 0.7 mm) and the second one with a finer mesh of 27 million cells (characteristic cell size of 0.35 mm). As reported in Fig. 8 , the highest ratios of unresolved turbulent kinetic energy arise around peak systole, where M superior to 20% is found for up to 4.1% and 2.6% of the nodes, respectively.

Figure 8: Temporal evolution along a cycle of the nodes fraction having $M>0.2$.

Concerning the near-wall resolution over time, an averaged $y^{+}=1.3$ was reported for the finer mesh, where again the peak systole led to the highest value with a mean $y^{+}=1.6$ and a maximal $y^{+}=5.6$. For both meshes, the higher $y^{+}$ were observed in the collateral narrowing and in the jet outflowing around systole, as shown in Fig. 9.

Figure 9: $y^{+}$ field at A. peak systole and B. end diastole for the fine (characteristic cell size of 0.7 mm) and the finer (characteristic cell size of 0.35 mm) mesh.

The resolved turbulent kinetic energy $k_{res}$ and the total kinetic energy $k_{total}$ are displayed for the two meshes in Figure 10. For the fine mesh, viscosity is added mainly in the jet, as well as around sharp angles (e.g. the junctions in the collateral). For the finer mesh, almost no viscosity is added.

Figure 10: Resolved ($k_{res}$) and total ($k_{total}$) kinetic energies at A. peak systole and B. end diastole for the fine (characteristic cell size of 0.7 mm) and finer (characteristic cell size of 0.35 mm) meshes.

Finally, the turbulence intensity $TI$ was computed as:

$TI(x,t)={\frac {\sqrt {{\frac {2}{3}}k_{res}(x,t)}}{\vert \vert {\overline {\mathbf {u} }}\vert \vert (x,t)}}$

where $\vert \vert {\overline {\mathbf {u} }}\vert \vert (x,t)$ is the magnitude of the phase-averaged velocity vector.

The results are displayed in Fig. 11. For both meshes, at peak systole the turbulence intensity was concentrated in the aneurysm-like region and at the jet issuing from the collateral. At end diastole, high turbulence intensity was found all over the phantom, especially at the outlet, U-bend and aneurysm.

Figure 11: Turbulence intensity at A. peak systole and B. end diastole for the fine (characteristic cell size of 0.7 mm) and finer (characteristic cell size of 0.35 mm) meshes.

The relative error on the phase-averaged velocity magnitude between the two simulations amounts to 0.9% of the maximum velocity magnitude found for the finer mesh. Hence, the finer mesh was considered fine enough to resolve the turbulence.

---

Contributed by: M. Garreau^a, T. Puiseux^a,b, R. Moreno^c, S. Mendez^a, F. Nicoud^a — 
^aIMAG, University of Montpellier, CNRS UMR 5149, Montpellier, France
^bSpin Up, Toulouse, France
^cI2MC, INSERM/UPS UMR 1297, Toulouse, France
^dALARA Expertise, Strasbourg, France

Front Page

Description

Test Data

CFD Simulations

Evaluation

Best Practice Advice

© copyright ERCOFTAC 2021