A pulsatile 3D flow relevant to thoracic hemodynamics: CFD - 4D 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

Left: Digital geometry (Dimensions are given with respect to the center of the collateral), Right: Mesh at the inlet

Solution Strategy

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

$\nabla \cdot \mathbf {u} =0$

${\frac {\partial \mathbf {u} }{\partial t}}+\mathbf {u} \cdot \nabla \mathbf {u} =-{\frac {1}{\rho }}\nabla p+\nu \nabla ^{2}\mathbf {u}$

where $\mathbf {u} ,\ p,\ \rho =1020\ kg/m^{3}$ and $\nu =4.02\times 10^{-6}m^{2}/s$ are the velocity field, pressure field, density and 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. This 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 zero velocity was imposed at the solid boundaries.

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. To do so a mask was applied in order to impose a linear evolution of the near-wall through-plane velocity with the constraint of zero-velocity at the wall. Thereby for the pixels located at position $r\in [R_{lim},R]$ (where $R$ is the radius of the circular cross-sectional inlet surface and $R_{lim}=0.95R$ an user-defined inner radius), the masked through-plane velocity component $w$ was set equal to:

$w(r)={\frac {w(R_{lim})}{R-R_{lim}}}(R-r)$

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. Thereby an appropriate way seem 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 instantaneous velocity vector and $N=30$ the total number of cycles of period $T=1s$ . As it 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, both modalities have to be expressed on the same grid to allow quantitative comparison. A first 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 coarse 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 fields. The steps of the downsampling process (Fig. 6) were the following: first the MRI grid was subsampled with 5 voxels in each direction (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 downsampling process

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.