Haemodynamics in a realistic human left heart: Multi-cycle large eddy simulation

1 Introduction Today, phase-contrast magnetic resonance imaging (PC-MRI) is able to provide 3D in-vivo images of the haemodynamics in heart chambers in a non-invasive way [1]. Nevertheless, accuracy of the velocity field is questionable in disturbed flow [2]. In addition, flow quantities such as wall shear stresses and pressure, which constitute relevant clinical information [3], cannot be measured directly with the PC-MRI method. Computational fluid dynamics (CFD) constitutes a possible, but challenging alternative for patient-specific non-invasive flow characterisation. CFD provides data with higher spatio-temporal resolution than actual in-vivo data produced by PC-MRI and gives the opportunity to broadly examine data sometimes inaccessible to medical imaging. The interactions between blood flow and the endocardium have been studied thanks to fluid-structure interaction (FSI) models [4]. However, the required mechanical properties of the heart muscle are virtually unknown and interactions of the heart with its environment have been neglected so far. An alternative consists in using dynamic morphological medical images (4D) for the computation of the heart deformations. In this view, endocardium movements are extracted from morphological medical images and prescribed as boundary conditions of the fluid problem, thus avoiding solving the FSI problem [5]. The present work describes a CFD framework based on this principle and its application to a left heart. 2 CFD framework The YALES2 solver [6] used in this study is a fully explicit 3D CFD in-house code solving the unsteady, Navier-Stokes equations on unstructured grids. This code was widely validated for engineering applications [7]. It is based on a spatial fourth-order finite volume approximation. Time discretisation is performed using an explicit fourth-order Runge-Kutta scheme. Due to the transitional nature of the flow, large eddy simulations are performed using advanced subgrid scale models able to handle wall bounded flows in complex geometries [8][9]. The mesh motion is handled by an arbitrary Lagrangian Eulerian method. The numerical domain is extracted from one particular medical image referred to as native image. The region of interest is segmented by thresholding method using an image processing software (ScanIP;Simpleware Ltd., Exeter, UK). Then, the extracted geometry is imported in a commercial mesher (Gambit, ANSYS) to generate the native numerical domain. Fields of deformations between the native image and the others images taken at different times in the cardiac cycle are computed by an image registration algorithm. These deformations are then applied to the native mesh, producing a set of successive mesh phases matching the physiological cardiac images at different times. Each phase is characterised by a mesh with a constant connectivity, and lasts until the next mesh phase is reached. The computation is stopped and the solution is interpolated to the mesh of the next phase. The computation carries on until the next phase and the same process is applied again. Surface movements are interpolated between each phase to prescribe boundary movements to the numerical domain for every time needed for the simulation. This principle was successfully applied to large vessels as aortas [5]. In this study, it allows us to handle the cardiac chambers and their connected vessels. However, other geometrical elements such as valves are out of reach with this method. As valves are expected to have an impact on the flow structure, they must be accounted for in the simulation. In this study, both the aortic (AV) and mitral (MV) valves are handled by an immersed boundary method (IBM). The AV which has a moderate impact on the ventricular flow is modelled as a planar region being alternatively permeable and impermeable depending on the phase in the cardiac cycle. The MV is modelled by a more realistic geometry based on measurements (position, open surface, leaflets lengths) on the medical images. 3 Simulation A CT exam from a patient is used to apply the numerical procedure to a left heart: 10 CT-scan images are available along the cardiac cycle, which lasts 1000 ms. From the 10 images, the numerical domain and its deformations along the cardiac cycle are computed using the procedure described in the former section. The domain includes the left atrium (LA), left ventricle (LV), the aortic root (AO) and the pulmonary veins. The spatial resolution is imposed to be close to 0.8 mm in all three spatial directions along the cycle, which yields grids of approximately three-million tetrahedral elements. The simulation time step is fixed by a CFL condition, which corresponds to a time step simulation of order of 10−4 s. Blood is modelled as an incompressible Newtonian fluid. A full video of the cardiac cycle can be seen here: youtu.be/1ze6ZxrSDHw. Two periods can be distinguished during the cardiac cycle: the diastolic period (the AV is closed, the MV open), lasting 640 ms for this heart and the systolic one (the AV is open, the MV closed). The IBM source terms mimicking the valves are adjusted to be in phase with the LV volume variations. As either the AV or the MV is closed during the cardiac cycle, the flow waveform imposed to the four inlet conditions of the computational domain can be calculated based on the mass conservation principle. The resulting inlet velocity varies in time during the cycle and is periodic. The inlet Reynolds number for each pulmonary vein varies from 0 to 600, approximately. During diastole, the LV volume increases and blood passes from the LA to the LV, forming a strong jet through the MV while an intense swirling flow can be seen in the atrium. The maximum Reynolds number for the jet is approximately 5000. The jet enters the LV and is deflected along the lateral wall cavity. At the end of diastole the MV closes and the AV opens: systole starts. The LV contracts and its volume decreases by 40%, resulting in a rapid ejection of the blood to the AO. Due to the transitional nature of this complex cyclic flow, the velocity pattern changes from one cycle to another. 20 cardiac cycles are simulated and phase-averages are gathered over the last 15 cycles. During the late diastole, the instantaneous velocity and phase-averaged velocity are different. While a highly disturbed flow is seen for the instantaneous view in the LV, a large recirculation cell and two others smaller recirculation areas are visible in the phase-averaged view. Interestingly, phase-averaged velocity flow fields are very similar to recent PC-MRI measurement of the flow in a left heart [1]. In the full paper, detailed description of this complex flow will be provided, focusing notably on cycle-to-cycle variations in the velocity field. References [1] Markl, M. et al., Comprehensive 4D velocity mapping of the heart and great vessels by cardiovascular magnetic resonance, J. of Cardiovasc. Magnetic Res., Vol. 13, No. 1, 2011, pp. 7. [2] Hollnagel, D. I. et al., Comparative velocity investigations in cerebral arteries and aneurysms: 3D phase-contrast MR angiography, laser Doppler velocimetry and computational fluid dynamics, NMR in Biomed., Vol. 22, No. 8, 2009, pp. 795-808. [3] Chien, S., Mechanotransduction and endothelial cell homeostasis: the wisdom of the cell, American J. of Physiology - Heart and Circulatory Physiology, Vol. 292, No. 3, 2007, pp. H1209-H1224. [4] Tang, D. et al., Patient-Specific MRI-Based 3D FSI RV/LV/Patch Models for Pulmonary Valve Replacement Surgery and Patch Optimization, J. of Biomech. Eng., Vol. 130, No. 4, 2008, pp. 041010. [5] Midulla, M. et al., Haemodynamic imaging of thoracic stent-grafts by computational fluid dynamics (CFD): presentation of a patient-specific method combining magnetic resonance imaging and numerical simulations, Eu. Radiology, Vol. 22, 2012, pp. 2094-2102. [6] Moureau, V. et al., Design of a massively parallel CFD code for complex geometries, Comptes Rendus Mécanique, Vol. 339, No. 2 - 3, 2011, pp. 141-148. [7] Moureau, V. et al., From Large-Eddy Simulation to Direct Numerical Simulation of a lean premixed swirl flame: Filtered laminar flame-PDF modeling, Combustion and Flame,Vol. 158, No. 7, 2011, pp. 1340-1357. [8] Germano, M. et al., A dynamic subgrid-scale eddy viscosity model, Physics of Fluids, Vol. 3, No. 7,1991, pp. 1760-1765. [9] Nicoud, F. et al., Using singular values to build a subgrid-scale model for large eddy simulations, Physics of Fluids, Vol. 23, No. 8, 2011, pp. 085106