Christophe Chnafa
email: christophe.chnafa@univ-montp2.fr

Simon Mendez

Franck Nicoud

Ao = Aorta

Rspv = Right

Image-based large-eddy simulation in a realistic left heart

Keywords: Patient-specific, Hemodynamics, Heart, Turbulence, Image registration, Arbitrary Lagrangian-Eulerian ALE = Arbitrary Lagrangian-Eulerian, PC-MRI = Phase-Contrast Magnetic Resonance Imaging, CFD = Computational Fluid Dynamics, LV = Left Ventricle, EFSI Tomography, LA = Left Atrium

A numerical framework allowing insight in fluid dynamics inside patientspecific human hearts is presented. The heart cavities and their wall dynamics are extracted from medical images, with the help of a non-linear image registration algorithm, in order to obtain a patient-specific moving numerical domain. Flow equations are written on a conformal moving computational domain, using an Arbitrary Lagrangian-Eulerian framework. Resulting equations are solved numerically with a fourth-order finite-volume technique.

Application of this framework to compute a patient-specific left heart flow is presented as well. The blood flow is characterized by its transitional nature, resulting in a complex cyclic flow. Flow dynamics is analysed in order to reveal the main fluid phenomena and to obtain insights into the physiological patterns commonly detected. It is demonstrated that the flow is neither laminar nor fully turbulent, thus justifying a posteriori the use of Large Eddy Simulation.

PC-MRI studies have contributed to the understanding of hemodynamic features [START_REF] Kim | Left ventricular blood flow patterns in normal subjects: A quantitative analysis by three-dimensional magnetic resonance velocity mapping[END_REF][START_REF] Kilner | Asymmetric redirection of flow through the heart[END_REF][START_REF] Fyrenius | Three dimensional flow in the human left atrium[END_REF][START_REF] Bolger | Transit of blood flow through the human left ventricle mapped by cardiovascular magnetic resonance[END_REF][START_REF] Markl | Comprehensive 4D velocity mapping of the heart and great vessels by cardiovascular magnetic resonance[END_REF][START_REF] Carlsson | Quantification of left and right ventricular kinetic energy using four-dimensional intracardiac magnetic resonance imaging flow measurements[END_REF][START_REF] Eriksson | Four-dimensional blood flow-specific markers of LV dysfunction in dilated cardiomyopathy[END_REF]. Although very comprehensive, the PC-MRI velocity mapping is not real-time. Hence, beat-to-beat variations in the flow cannot be recorded (the k-space is filled over many cardiac cycles). Moreover, PC-MRI suffers from a relatively low spatio-temporal resolution, precluding the observation of small-scale and fast time-varying flow features [START_REF] Zhao | Comparative study of magnetic resonance imaging and image-based computational fluid dynamics for quantification of pulsatile flow in a carotid bifurcation phantom[END_REF][START_REF] Hollnagel | Comparative velocity investigations in cerebral arteries and aneurysms: 3D phasecontrast MR angiography, laser doppler velocimetry and computational fluid dynamics[END_REF].

Echocardiography techniques [START_REF] Kim | Development and validation of echo PIV[END_REF][START_REF] Garcia | Two-dimensional intraventricular flow mapping by digital processing conventional color-doppler echocardiography images[END_REF], with higher spatio-temporal resolution make an alternative to PC-MRI. However, echocardiography only gives access to velocity components directed towards or away from the ultrasonic beam, while one would want to measure the full 3D flow vectors. Nevertheless, investigations have been conducted on normal and abnormal hearts leading to potential hemodynamics-based biomarkers for cardiac health. [START_REF] Gharib | Optimal vortex formation as an index of cardiac health[END_REF][START_REF] Hong | Characterization and quantification of vortex flow in the human left ventricle by contrast echocardiography using vector particle image velocimetry[END_REF][START_REF] Cimino | In vivo analysis of intraventricular fluid dynamics in healthy hearts[END_REF].

In order to obtain more information about the heart hemodynamics, in vitro investigations have been performed in fully controlled experiments [START_REF] Kheradvar | On mitral valve dynamics and its connection to early diastolic flow[END_REF][START_REF] Querzoli | Effect of the prosthetic mitral valve on vortex dynamics and turbulence of the left ventricular flow[END_REF][START_REF] Tanné | Assessment of left heart and pulmonary circulation flow dynamics by a new pulsed mock circulatory system[END_REF][START_REF] Falahatpisheh | High-speed particle image velocimetry to assess cardiac fluid dynamics in vitro: From performance to validation[END_REF]. Blood patterns in heart chambers replications have been studied thanks to particle image velocimetry in healthy and abnormal configurations.

In addition to these studies, computational fluid dynamics (CFD) has been more and more used to predict blood flow in the heart over the last decade. In silico replications of heart chambers, mainly the left ventricle (LV), have been considered. Simulations in idealised ventricles [START_REF] Baccani | Fluid dynamics of the left ventricular filling in dilated cardiomyopathy[END_REF][START_REF] Pedrizzetti | Nature optimizes the swirling flow in the human left ventricle[END_REF][START_REF] Domenichini | Three-dimensional filling flow into a model left ventricle[END_REF] or in more realistic geometry [START_REF] Zheng | Computational modeling and analysis of intracardiac flows in simple models of the left ventricle[END_REF] have been performed. As in vitro experiments, such fundamental CFD studies are particularly useful to isolate and elucidate the effect of well-controlled parameters on the blood flow. However, inherent simplifications raise the question of the relevance of their conclusions for individual clinical cases. In this context, numerical simulation using a combination of computational methods and medical imaging techniques for determining vascular geometry appears to be a relevant strategy. CFD starts to be a mature technique for arterial flows [START_REF] Gerbeau | Fluid-structure interaction in blood flows on geometries based on medical imaging[END_REF][START_REF] Taylor | Patient-specific modeling of cardiovascular mechanics[END_REF][START_REF] Sforza | Hemodynamics of cerebral aneurysms[END_REF][START_REF] Les | Quantification of hemodynamics in abdominal aortic aneurysms during rest and exercise using magnetic resonance imaging and computational fluid dynamics[END_REF], but its application to heart haemodynamics faces additional challenges:

• the geometry of the blood domain is complex and it undergoes large deformations,

• opening and closing valves make the topology of the domain change over the cardiac cycle,

• the flow is the result of a complex electrical-fluid-structure interaction problem,

• the flow regime is most probably transitional between laminar and turbulent and varies over the cardiac cycle.

Two main different strategies have been developed to obtain simulation of the blood flow in realistic heart geometries. The most natural one is to extract the heart geometry at one chosen moment in the heart cycle and to solve an electrical-fluid-structure interaction (EFSI) problem [START_REF] Watanabe | Multiphysics simulation of left ventricular filling dynamics using fluid-structure interaction finite element method[END_REF][START_REF] Cheng | Fluid-structure coupled CFD simulation of the left ventricular flow during filling phase[END_REF][START_REF] Tang | Patient-specific MRI-based 3D FSI RV/LV/patch models for pulmonary valve replacement surgery and patch optimization[END_REF][START_REF] Einstein | Fluid-structure interactions of the mitral valve and left heart: comprehensive strategies, past, present and future[END_REF][START_REF] Nordsletten | Coupling multi-physics models to cardiac mechanics[END_REF][START_REF] Le | On the three-dimensional vortical structure of early diastolic flow in a patient-specific left ventricle[END_REF].

In this approach, patient-specific data are needed [START_REF] Smith | euheart: personalized and integrated cardiac care using patientspecific cardiovascular modelling[END_REF][START_REF] Vadakkumpadan | Image-based estimation of ventricular fiber orientations for personalized modeling of cardiac electrophysiology[END_REF]. What is the exact rheology of the myocardial muscle? What is the load produced by the heart environment? How to reproduce the mechano-electric coupling in the heart muscle? All these questions make such an approach extremely challenging. Another strategy consists in using realistic heart wall movements extracted from cine MRI or Computed Tomography (CT) scan data. Heart movement is not computed, but prescribed from the patient-specific medical images, which can be acquired using standard clinical imaging procedures. Such a computational approach, where the geometry and the movements are extracted from images will be referred to as image-based computational fluid dynamic (IB-CFD). Different research teams have developed IB-CFD methods for heart flows, more specifically to study the left ventricle alone [START_REF] Saber | Progress towards patient-specific computational flow modeling of the left heart via combination of magnetic resonance imaging with computational fluid dynamics[END_REF][START_REF] Merrifield | Com-54 bined CFD/MRI analysis of left ventricular flow[END_REF][START_REF] Long | Subject-specific computational simulation of left ventricular flow based on magnetic resonance imaging[END_REF][START_REF] Schenkel | MRIbased CFD analysis of flow in a human left ventricle: Methodology and application to a healthy heart[END_REF]. Recently, more advanced work has been published, using a full heart model obtained from CT images [START_REF] Mihalef | Patient-specific modelling of whole heart anatomy, dynamics and haemodynamics from four-dimensional cardiac CT images[END_REF] or a heart model fed from MR images [START_REF] Doenst | Fluid-dynamic modeling of the human left ventricle: Methodology and application to surgical ventricular reconstruction[END_REF]. The feasibility of cardiac IB-CFD has been shown, but the flow results notably suffered from limited spatial resolution or partial geometries (LV only in the majority of the cases). Furthermore, there has been remarkably little focus on the presence of turbulence in the heart, except in experimental works [START_REF] Querzoli | Effect of the prosthetic mitral valve on vortex dynamics and turbulence of the left ventricular flow[END_REF][START_REF] Falahatpisheh | High-speed particle image velocimetry to assess cardiac fluid dynamics in vitro: From performance to validation[END_REF].

In the present paper, an image-based CFD method is presented. As in the aforementioned IB-CFD works [START_REF] Saber | Progress towards patient-specific computational flow modeling of the left heart via combination of magnetic resonance imaging with computational fluid dynamics[END_REF][START_REF] Merrifield | Com-54 bined CFD/MRI analysis of left ventricular flow[END_REF][START_REF] Long | Subject-specific computational simulation of left ventricular flow based on magnetic resonance imaging[END_REF][START_REF] Schenkel | MRIbased CFD analysis of flow in a human left ventricle: Methodology and application to a healthy heart[END_REF], medical images are used to generate a moving patient-specific domain, in which the blood flow equations are solved. The geometry movements are generated from a 4D sequence (MRI or CT scan images) treated by an appropriate image registration algorithm [START_REF] Moreno | Non-linear transformation field to build moving meshes for patient specific blood flow simulations[END_REF][START_REF] Nicoud | Computational hemodynamics in moving geometries without solving the fluid-structure interaction problem[END_REF]. This approach has been used before to compute blood flow in aortas [START_REF] Midulla | Haemodynamic imaging of thoracic stent-grafts by computational fluid dynamics (CFD): presentation of a patient-specific method combining magnetic resonance imaging and numerical simulations[END_REF]. It is further developed for application to the left heart flow, notably by introducing valve modeling. In order to demonstrate the ability of the method to compute flow in the heart, its application to a full patient left heart using 4D CT scan images is presented. Extensive description of the flow obtained is provided and the main flow characteristics usually reported in the literature are emphasised. First numerical insight into turbulence in the left heart is reported as well.

The numerical method is detailed in § 2. The characteristics of the applied case is presented in § 3 and the flow field obtained is described in § 4.

Concluding remarks are given in § 5.

Methodology

In this section, the image-based computational fluid dynamics approach is detailed. First, the fluid problem resolution is detailed in an ALE framework.

Next, § 2.2 presents how the moving computational domain is obtained from the medical images. Specific valve modeling is needed and is the object of § 2.3. Finally the method to obtain the inflow and outflow boundary conditions is described.

Fluid problem

Governing equations

Blood can be modelled as an incompressible fluid, but red blood cells induce a complex rheological behaviour [START_REF] Formaggia | Cardiovascular Mathematics: Modeling and simulation of the circulatory system[END_REF]. However, for high stress levels and in large vessels, non-Newtonian effects are usually neglected and blood is usually modelled as an incompressible Newtonian fluid in numerical simulations [START_REF] Taylor | Patient-specific modeling of cardiovascular mechanics[END_REF][START_REF] Mihalef | Patient-specific modelling of whole heart anatomy, dynamics and haemodynamics from four-dimensional cardiac CT images[END_REF]. Taking the incompressible flow assumption into account and assuming blood as a Newtonian fluid, the fluid motion is described by the Navier-Stokes (NS) equations. These equations are solved on the moving

blood domain Ω f (t) ⊂ R 3 bounded by ∂Ω f (t). The boundary ∂Ω f (t) is such that ∂Ω f (t) = ∂Ω i f (t) ∪ ∂Ω w f (t) ∪ ∂Ω o f (t) and ∂Ω i f (t) ∩ ∂Ω w f (t) ∩ ∂Ω o f (t) = ∅ where ∂Ω i f (t)
∂u f ∂t + (u f • ∇)u f = - 1 ρ ∇p + ν∇ 2 u f + f, ∇ • u f = 0,      on Ω f (t) (1)
where u f is the fluid velocity field, p is the pressure field, ν the kinematic viscosity, ρ the density and f a volumetric force. The corresponding initial and boundary conditions are,

u f (x, 0) = u 0 f (x) on Ω f (0), (2)
u f (x, t)| x∈∂Ω w f (t) = u s (x, t) on ∂Ω w f (t), (3)
u f (x, t)| x∈∂Ω i f (t) = -U in (x, t)n o (x) on ∂Ω i f (t), (4)
where U in (x, t) is the inlet velocity profile imposed as a Dirichlet condition, n o the outward normal at the inlet faces, and u s is the endocardium surface velocity field imposed as a Dirichlet condition as well. A convective outlet boundary condition is imposed on ∂Ω o f (t) as,

∂u f (x, t) ∂t + U conv ∂u f (x, t) ∂n = 0, (5)
where n is the outward normal at the outlet patch and U conv the convective velocity. The uniform convective velocity U conv is imposed in such a way to meet the global mass conservation over Ω f (t). The surface velocity u s is not computed but extracted from the medical images and applied as boundary conditions for the fluid problem (see section 2.2.2).

Time advancement scheme

The time advancement scheme is an explicit low-storage four-step Runge-Kutta scheme [START_REF] Williamson | Low-storage runge-kutta schemes[END_REF] recast in an ALE formalism and coupled with the Chorin's projection correction method [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF] for the pressure term. The grid is displaced during the prediction step only.

The starting point for deriving the time-advancement scheme is the integral form of the NS equations on a node-centred control volume ω(t), its boundary S t moving with the mesh velocity u g . The equations integrated in time between t n =t 0 + n∆t (∆t being the time step size) and t n+1 are classically [START_REF] Donea | Arbitrary Lagrangian-Eulerian methods[END_REF]:

t n+1 t n d dt ω(t) u dωdt + t n+1 t n ω(t) ∇ • ((u f -u g )u f ) dωdt = RHS, (6)
where the RHS containing the viscous fluxes and the pressure gradient is omitted in the following. The four sub-steps of the time advancement are computed as:

u 0 f = u n f , u i f = u n f Ω n Ω i + α i ∆t Ω i ω(t) ∇ • (u i-1 f -u n+1 g )u i-1 f dω for i = 1, ..., 4
u n+1 ⋆ = u 4 f , (7)
where

u n+1 ⋆ =u ⋆ (t n+1 ) is the predicted (non-solenoidal) velocity field, α i is a coefficient as α i = [1/4, 1/3, 1/2, 1] and Ω n (resp. Ω i ) the cell volume at time t n (resp. t i =t n + α i ∆t).
The presented numerical method must satisfy a discrete version of the socalled Geometric Conservation Law (GCL) [START_REF] Farhat | The discrete geometric conservation law and the nonlinear stability of ALE schemes for the solution of flow problems on moving grids[END_REF][START_REF] Boffi | Stability and geometric conservation laws for ale formulations[END_REF]. Referring to the classical interpretation of the GCL, the numerical method has to preserve the state of a constant flow. Considering a constant flow c, a discrete GCL arises from the presented numerical scheme for each sub-step as:

Ω i -Ω n = -α i ∆t ω(t) ∇ • u n+1 g dω for i = 1, ..., 4. (8)
Grid nodal velocity u g is constant and the metrics evolve linearly during the computational time step. In order to satisfy exactly the Eq. ( 8), the integration is computed at the midpoint configuration for each sub-step. Thus, the numerical scheme [START_REF] Eriksson | Four-dimensional blood flow-specific markers of LV dysfunction in dilated cardiomyopathy[END_REF] satisfies the Discrete GCL for each sub-step.

At the end of this prediction step, the grid reaches the final position of the time step. Hence, the projection step to calculate pressure is performed over a fixed grid. A Deflated Preconditioned Conjugate Gradient algorithm is used to solve the Poisson equation [START_REF] Malandin | Optimization of the deflated conjugate gradient algorithm for the solving of elliptic equations on massively parallel machines[END_REF] involved in the projection step.

A homogeneous Neumann condition is applied for the pressure calculation and the pressure constant is fixed as the averaged pressure in the numerical domain so that the volume-averaged pressure over Ω f (t) is zero.

This ALE method has been implemented in the YALES2BIO 1 solver [START_REF] Mendez | An unstructured solver for simulations of deformable particles in flows at arbitrary reynolds numbers[END_REF].

This code is an adaptation of the YALES2 2 [START_REF] Moureau | Design of a massively parallel CFD code for complex geometries[END_REF] research solver for biomedical applications. As such it inherits from YALES2 4th-order central scheme in space on unstructured meshes and its capabilities of massively parallel computations of turbulent flows [START_REF] Moureau | From large-eddy simulation to direct numerical simulation of a lean premixed swirl flame: Filtered laminar flame-PDF modeling[END_REF].

Large Eddy Simulation

Direct numerical simulation (DNS) is suitable for low Reynolds number flow: all the large and small scales are resolved in this approach, requiring approach is used. While the large scales of the flow are resolved, the smaller subgrid scales (SGS) are modelled using SGS models, usually eddy-viscositybased models.

In this work, a SGS model able to handle wall bounded flows in complex geometries is used [START_REF] Nicoud | Using singular values to build a subgrid-scale model for large eddy simulations[END_REF] as well as a well established Dynamic Smagorinsky formulation [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF]. Filtering is here performed implicitly, the filter width being equal to the grid spacing. When using the Dynamic Smagorinsky model, coefficient C s obtained by the dynamic procedure was locally averaged and the obtained SGS viscosity satisfies the condition ν + ν SGS > 0.0.

Computational Domain

The left heart

The normal left heart has two cavities: the left atrium (LA) which collects blood from four pulmonary veins and the left ventricle (LV) that receives blood from the LA to propel it into the primary blood system through the aorta (AO). The four pulmonary veins are distinct vessels named by their position: the right superior pulmonary vein (RSPV), the right inferior pulmonary vein (RIPV), the left superior pulmonary vein (LSPV) and finally the left inferior pulmonary vein (LIPV).

Extraction of the heart motion

In order to extract the movement of the endocardium from the medical images, image registration is used. It is a classical image processing technique, which consists in finding a transformation to map two images called to the template image and the reference image. It is commonly applied to cardiac images [START_REF] Makela | A review of cardiac image registration methods[END_REF], notably to obtain clinical information about the myocardial contractile function. Here, a classical image registration problem [START_REF] Modersitzki | Numerical Methods for Image Registration[END_REF] is solved to reconstruct the movement of the frontiers of the computational domain. The algorithm used is based on a voxel similarity measures method [START_REF] Maintz | A survey of medical image registration[END_REF], which directly relies on the images grey levels.

A set of N 3D images I i of a patient heart are taken at different times t i , 0 ≤ t 0 , t 1 , ..., t N -1 < T during the heart cycle of period T . One of the N images is selected as a template image and N -1 transformations ψ i between the template image and the i-th image I i (successively considered as the reference image) are computed. Transformation ψ i is searched so that I 0 (ψ i (x)) = I i (x) for each voxel, where I i (x) stands for the voxel grey-level value at position x in the i-th 3D discrete image. Each transformation ψ i is computed through an optimisation problem reading: given two images I 0 and I i , find a mapping ψ i such that the squared intensity difference d between I 0 (ψ i (x)) and I i (x) is minimized for each voxel. However, the problem is illposed and the solution ψ i is not unique: additional constraints are needed.

A constraint based on the Jacobian matrix of the local voxel deformation J and its singular values is applied on ψ i , thanks to prior knowledge of the sought deformation [START_REF] Ashburner | Incorporating prior knowledge into image registration[END_REF]. A detailed description of this procedure can be found in [START_REF] Ashburner | High-dimensional image registration using symmetric priors[END_REF]. This additional constraint plays the role of a regularization term. The multi-objective problem is dealt with by scalarization, a parameter λ weighting the regularization term (see Eq. 11). A compromise for the value of λ has to be found (here empirically) for each set of images, in order to avoid partial deformations (λ too high) or non-physical deformations due to excessive warping (λ too small).

In addition, in order to avoid local minima and to promote smooth transformations, a two-loop algorithm is used, with evolving weighting of the objective functions (inner loop) and progressively refined image filtering (outer loop). During computation of the deformations, d, the sum of the squared intensity difference of each voxel, is weighted. For the j-th iteration of the inner loop, the weight is defined as the inverse of the residual sum of the squared differences computed at the previous inner iteration j -1 and is denoted by 1/d j-1 (see Eq. 9). Because d j-1 has a high value for the first iterations, more weight is initially given to the regularization term, in order to get smooth deformations. As the algorithm gets close to the final solution, d decreases, giving less weight to the priors and letting the algorithm compute more precise deformations.

In order to increase the likelihood of achieving a good solution, a Gaussian pyramid approach is used [START_REF] Adelson | Pyramid methods in image processing[END_REF]. A number of outer iterations consisting in successive lowpass filtering are performed, using a Gaussian smoothing convolution kernel G σ of width σ k decreasing with k. Finally, the sought transformation ψ j,k i at the inner iteration j and the outer iteration k minimizes the function f j,k defined as:

f j,k 1 (ψ j,k i ) = 1 2d j-1 Ω ([G σ k * I 0 ](ψ j,k i (x)) -[G σ k * I i ](x)) 2 dΩ, (9)
f j,k 2 (ψ j,k i ) = g(J(ψ j,k i )), (10)
f j,k (ψ j,k i ) = f j,k 1 (ψ j,k i ) + λf j,k 2 (ψ j,k i ), (11)
where f j,k 1 (ψ j,k i ) is the term related to the difference in grey levels and f j,k 2 (ψ j,k i ) the regularization term. The function g is computed from the Jacobian matrix singular values and determinant. Technical details can be found in [START_REF] Ashburner | Incorporating prior knowledge into image registration[END_REF].

Application of the patient-specific movement to the computational left heart grid

The template image corresponding to volumetric data I 0 is imported into an image processing software (ScanIP;Simpleware Ltd., Exeter, UK). A preprocessing step is applied on the image before segmentation. The region of interest in the image I 0 is isolated and the image is smoothed thanks to a Gaussian filter to erase noise inherent to the medical images protocols. The segmentation itself is done by a classical thresholding method [START_REF] Pham | Current methods in medical image segmentation[END_REF]. A suitable image intensity range which encompasses the voxel intensities of the region occupied by blood in the heart is selected. Note that the quality of the segmentation relies, at least partially, on the operator's skills and knowledge of the human heart morphology and the quality of the chosen medical image. Image resolution plays a role in the potential domain simplification at this step. Trabeculae, left atrial appendage or papillary muscles can either be kept or neglected depending on the image quality and spatial resolution available. The 3D geometric reconstruction covers all the space occupied by blood in the left heart cavities. The surface of the geometric reconstruction of the heart is triangulated.

For each couple of images (I 0 ,I i ) a 3D deformation field ψ i is computed as described in the previous section. If wanted, the deformation fields can also be computed for couples defined as (I i-1 ,I i ); this approach is mathematically equivalent. Trilinear interpolation from each of these deformation fields to the template surface mesh is done. Thus, a set of N -1 successive surface meshes matching the physiological cardiac images at different times t i is produced as schematized in Fig. 1.

Position and velocity of all surface points are needed at any discrete time of the simulation, not only at the times t 0 , t 1 , ..., t N -1 . Since all the generated surface meshes share the same topology (number and connection between nodes, number of cells), a trigonometric interpolation is used to compute the position and the velocity of each node, surface position and velocity then read:

Figure 1: Mesh deformation procedure applied to a left human heart. The template mesh segmented from the image at time t 0 is deformed thanks to ψ i to obtain the mesh at time t i . This procedure is done for each image in the cardiac cycle in order to obtain the corresponding meshes.

x s (t) = m i=0 [a i cos(2iπ t T ) + b i sin(2iπ t T )], u s (t) = 2iπ T m i=1 [-a i sin(2iπ t T ) + b i cos(2iπ t T )],          on ∂Ω f (t) ( 12
)
where m is the number of Fourier modes (m = (N -1)/2 or N/2 depending on the parity of N) and a i , b i the Fourier coefficients. Thus, needed surface velocity u s (Eq.3) of the patient heart is entirely extracted from the medical images.

Volumetric grid

The template surface is imported in a commercial mesher (Gambit, AN-SYS) to generate a template unstructured tetrahedral mesh. The computed boundary Fourier coefficients of Eq. ( 12) are interpolated in this template numerical domain. The computational mesh boundary now follows the shape of the patient endocardium and is updated in every step of the simulation.

Motion of all internal points in the computational mesh is deduced from this prescribed boundary motion thanks to a harmonic extension of u g onto the numerical domain. At each iteration, nodal velocity u g is calculated through the computation of the following problem [START_REF] Masud | An adaptive mesh rezoning scheme for moving boundary flows and fluid-structure interaction[END_REF],

∇ • (k(x)∇u g (x)) = 0, on Ω f (t) (13)
and,

u g (x)| x∈∂Ω f (t) = u s (x) on ∂Ω f (t), (14)
where k(x) is the displacement diffusion coefficient. This coefficient is determined to preserve a good computational grid quality. The larger elements will distort at a faster rate than the smaller ones -a desirable feature for our application in order to preserve grid quality. The grid quality is monitored during the simulations. When the boundary displacement becomes too large compared to the local cell sizes, the cell quality can be highly deteriorated and some cells can even become degenerated. This can lead to convergence problems or negative cell volumes. In this case, the numerical domain is remeshed and the fluid solution is interpolated on this new discretised domain.

Fourier coefficients of Eq. ( 12) are interpolated as well.

Valves

Since the valves are thin highly-moving structures, their precise movements are not extractable from MRI or CT scan exams. Besides, on the numerical side, incorporation of moving valve leaflets in the grid topology would most probably imply a grid quality degradation, making numerous re-meshing mandatory. Hence, it was chosen to model the valves using an immersed boundary technique.

Valves annuli geometries are reconstructed by visual inspections of the medical images. The annulus geometry is represented as a cloud of points p i , whose nodal coordinates are used to define the valves annulus within the numerical domain. These markers allow the reconstruction of the aortic and mitral annuli.

The shape of the aortic valve is simply approximated by the plane passing through the set of markers and obtained by a least squares method. As the focus of the study is on the flow in the atrium and the ventricle, it is not mandatory to develop a more complex model for the aortic valve, since its morphology has very limited effects on the LA/LV haemodynamics. Note also that the physiological aortic valve open position offers virtually no resistance to the main flow passing through the aorta from the LV. The regression plane modeling the aortic valve is thus made alternatively fully permeable or impermeable depending on the phase in the cardiac cycle.

The mitral valve is represented by a more complex model, since its shape is expected to strongly impact the LV haemodynamics. A mitral plane P α is defined to fit the set of markers p i with a least squares method and the mitral annulus geometric center C α is computed. A plane P β parallel to P α is defined at a distance l(t). When the MV is open, it is assumed that the cross section area seen by blood is elliptical. An ellipse ǫ of axes a(t), b(t)

is defined on P β and its angle in its plane set manually to fit the medical images. The ellipse center C ǫ is not a direct projection of C α on plane P β , an eccentricity e(t) is defined. The eccentricity e(t) is the distance between the projection of C α on plane P β and the ellipse center C ǫ . Leaflets are considered as the surface linking the mitral annulus to the elliptical opening ǫ. Position of the annulus markers are projected on the ellipse. Therefore, the leaflets surface is reconstructed by triangles pieces between the annulus markers and the projected annulus markers on the ellipse (see Fig. 2). To summarize, the quantities needed to feed this model are:

• the mitral annulus markers p i ,

• the average leaflets length l(t) which is the distance between C α and the generated plane P β ,

• the eccentricity e(t) of the ellipse center C ǫ compared to the projection of C α on plane P β ,

• the ellipse axis a(t), b(t) and its angle in the P β plane.

The scheme model and its numerical representation are shown in Fig. • MV open, AV closed: this corresponds to the diastolic phase.

Note that MV and AV happen to be closed at the same time in physiological conditions. This corresponds to the isovolumic phases which last for only a small fraction of the cardiac cycle. Given the poor time resolution reachable by 4D medical imaging systems (IRM and CT scan), the isovolumic phase cannot be described accurately anyway. Assuming that the MV and AV are never closed at the same time is thus an acceptable assumption given the accuracy of the medical data used to feed the CFD solver.

Knowing the MV leaflets position during the heart cycle, their effect on the blood flow is accounted for thanks to an immersed boundary method (IBM) [START_REF] Mohd-Yusof | Combined immersed-boundary/B-spline methods for simulations of flow in complex geometries[END_REF]. For this purpose, the leaflets representations (plane for AV and triangle by parts for the MV) are first given a thickness s so that a few mesh nodes are located within the valves.

Grid nodes are tested to decide whether they are in the leaflets or not. For the AV, the distance from the AV plane is computed. If the node distance is

within the closed interval [-s 2 , s 2
] the node is tagged as belonging to the AV valve. For the MV, each node being at a distance within [-s 2 , s 2 ] of one of the MV triangle are treated. The node is projected on the corresponding triangle and the triangle barycentric coordinates are computed. The node is tagged as belonging to the MV valve if all of the barycentric coordinates are greater than zero. Then, a force f in the NS equations ( 1) is set to impose the fluid velocity to zero within the leaflets. Let us assume that the NS equations are discretized in time as:

u n+1 f -u n f ∆t = RHS n + f n , ( 15
)
where ∆t is the simulation time step, f n the force imposed at time t n and RHS n the right hand side containing the convective, viscous and pressure terms at t n . Imposing the velocity to be zero at time t n+1 requires:

f n = -RHS n - u n f ∆t . ( 16
)
This expression is used at each mesh node lying within one of the valve region. The finite-volume scheme used being node centred and since the force f n is imposed at nodes lying into valve leaflets, a null velocity is imposed in the entire dual cell where the governing equations are integrated. Thus, there is no interpolation of the forcing term and the valves geometry is described in a stepwise way. f n is set to zero anywhere else.

Note that, as the opening/closing of the valve is not resolved, valves switch instantaneously from closed position to open position and vice versa. When open, some small displacements may be seen because of the displacement of the valve annulus, but they are here neglected. As a consequence, the source term used mimics the presence of valves as fixed obstacles.

Inlet and outlet boundary conditions

Q i (t) = dV 1 dt , (17)
where Q i (t) is the inlet flow rate to be imposed at the pulmonary veins V 1 , V 2 and V 3 are respectively the LA, LV and the aortic root volumes. In the second case MV open, AV closed, mass conservation imposes:

Q i (t) = dV 1 dt + dV 2 dt , (18)
Therefore, the flow rates are only determined by the time evolution of the studied heart geometry, which is directly related to the medical images. The function U in j (x, t) is used for the j-th inlet condition (j=1,..,4). Functions can be modified to impose velocity profiles if desired. For uniform inflows, the Dirichlet condition (see Eq. 4) reads:

U in j (t) = Q i (t) A j ζ j , (19)
Where A j is the cross-section area of the j-th pulmonary vein and ζ j the flow distribution between the four pulmonary veins (ζ j =0.25 for all j corresponding to an equipartition of the inflow). The template geometry is extracted at a given time t 0 taken during the middiastole.

Due to limitations in spatial and temporal resolutions, the complex intracardiac geometry cannot be fully reconstructed from medical images. The acquisition frequency of the images forbids precise observation of potential fast heart movement. However, the analysis of the flow features in the next sections suggests that this frequency is sufficient for this application. As shown in Fig. 3, the template numerical domain extracted includes the LA, LV, the aortic root and four pulmonary veins. Left atrial appendum papillary muscles are omitted. The LV has a height of 8.8 cm from the MV to the apex (the lowest extremity of the LV in Fig. 3) and a maximum diameter of order of 5 cm. The LA has a height of 5.5 cm from the MV to the upper pulmonary root and a maximum diameter of order of 4 cm.

A focus on the left atrium and its four pulmonary veins is shown in Fig. 4 (left). Sketch of inflow angles of the four pulmonary veins is also shown, giving an idea of the flow trajectory in the atrium. Table 1 shows measures of the ostia in the atrium. Blood velocity at the LA enters will be influenced by these areas. The boundary Fourier coefficients of Eq. ( 12) are computed from the obtained meshes.

RSPV

Computational mesh and simulation details

A nearly isotropic grid is created from the heart model described in the above section using the commercial software Ansys Gambit which was selected for its ability to generate good-quality tetrahedral mesh, appropriate for finite-volume formulations. 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. Computed Fourier coefficients are interpolated on this model surface. In this study, re-meshing was used at each acquisition instants t i , i = 0, ..., 9 in order to ensure good mesh quality and small numerical errors over the whole cardiac cycle.

Valves are modelled as explained in section (2. The simulation time step is fixed by a CFL condition (CFL=0.9) consistent with the explicit time integration used in the CFD solver; this cor-Element responds to a time step varying from 3.0 × 10 -4 s to 8.0 × 10 -4 s in the present computation.

U max D = 2 A
The flow waveform imposed at the four inlet conditions of the computational domain is calculated based on the mass conservation principle as explained in section (2.4). Having no information on the velocity profile, it is assumed uniform for this simulation. The resulting flow rate varies in time during the cycle and is periodic. The partition of the inflow is not known in this case and detailed information about averaged flowrates repartition between PVs have not been found in the literature. From data presented by Dahl [START_REF] Dahl | Impact of pulmonary venous locations on the intra-atrial flow and the mitral valve plane velocity profile[END_REF], it has been chosen to set the flow rate at each of the pulmonary vein equally distributed (ζ j = 0.25). Note that the resulting blood flow, in particular in the atrium, may be sensitive to this repartition. Figure 5

Q (mL.s -1 ) Q o (mL.s -1 ) Q i (mL.s -1
) Figure 5: Flow rates at the aortic valve (top plot), mitral valve (middle plot) and the total heart inflow (bottom plot) imposed at the pulmonary veins during the heart cycle. Vertical dotted lines mark the limit between the systolic phase (t/T between 0.015 and 0.375) and the diastolic phase. The E wave, L wave and A wave are indicated on the mitral flow.

The pulmonary flow rate is coherent with classical medical data, the flow reversal associated with the atrial contraction being even visible from t/T =0.86 onward. The aortic valve flow rate behaves as expected: it increases during systole until its maximum (Q=320 mL.s -1 ) and then decreases until its shutting at t/T =0.375. The aortic flow rate stays null until the next systolic event. The mitral flow rate is usually composed by two peaks.

The E wave or rapid filling and a second one, the A wave corresponding to the late diastole separated by a phase with almost no heart motions called the diastasis. The flow rate shows three peaks in this case: the first one (t/T =0.51, Q=410 mL.s -1 ) represents the E wave, the last one a weak A wave (t/T =0.98, Q=75 mL.s -1 ). A weak third pic is also present during diastasis phase (t/T =0.78, Q=85 mL.s -1 ). This peak may correspond to the so called L wave [START_REF] Ha | Triphasic mitral inflow velocity with middiastolic filling: clinical implications and associated echocardiographic findings[END_REF]. [START_REF] Haugen | Velocity profiles in mitral blood flow based on three-dimensional freehand colour flow imaging acquired at high frame rate[END_REF]. The Reynolds number of the aortic valve is about 5300. These ranges of Reynolds numbers, the complex moving geometry of the heart and the pulsating nature of the inlet flow indicate that the nature of this complex cyclic flow may be transitional if not turbulent.

Phase-averaged and fluctuating velocity definitions

Using the framework described before, twenty cardiac cycles were simulated and phase-averages were gathered over the last n=15 cycles. The phase average over n cardiac cycles is defined as,

u f (x, t) = 1 n n-1 k=0 u f (x, t + kT ), ( 20
)
where u f (x, t) is the fluid velocity, x denotes the spatial coordinates and t the time. Root mean square (r.m.s.) velocity u rms is defined as,

u rms (x, t) = u 2 f (x, t) -u 2 f (x, t). (21)

LES quality assessment

The Kolmogorov length scale η was assessed using the fluctuating part of the velocity u ′ = u fu f (u f being defined in the section 3.3) over the last 15 heart cycles:

η = ν 3/4 /(2(ν + ν SGS )S ′ ij S ′ ij ) 1/4 (S ′ ij
is the symmetric part is the fluctuating strain rate tensor). Smallest length scales are of order of 2 × 10 -5 m in the domain and the spatial averaged value is of order of 10 -4 m. Such length scale makes the DNS of numerous heart cycles hardly reachable today: more than one billion cells would be needed to resolve all turbulence scales. Proper computation are still manageable thanks to LES, as discussed in the section 2.1.3.

In order to achieve a quality assessment of the presented simulation, the Pope criterion [START_REF] Pope | Ten questions concerning the large-eddy simulation of turbulent flows[END_REF] is used. The evaluation of the LES quality is estimated through the fraction of the turbulent kinetic energy in the resolved scales. This fraction M(x, t) is defined as,

M(x, t) = k tot (x, t) -k sgs (x, t) k tot (x, t) , ( 22
) = 1 - k sgs (x, t) k tot (x, t) , ( 23
)
where k sgs is the SGS kinetic energy, k res the resolved kinetic energy of the velocity fluctuations and k tot is the summation of k res and k sgs . Evaluation of this fraction M(x, t) requires computing the resolved turbulent kinetic energy which is defined as,

k res (x, t) = 1 2 u ′ 1 (x, t) 2 + u ′ 2 (x, t) 2 + u ′ 3 (x, t) 2 , ( 24
)
The SGS kinetic energy k sgs was estimated from the expression [START_REF] Sagaut | Large eddy simulation for incompressible flows[END_REF]:

k sgs (x, t) = ν sgs (x, t) 2 (∆(x, t)C s (x, t)) 2 , ( 25
)
where ν sgs is the SGS turbulent viscosity given by the dynamic Smagorinsky-Lilly model, C s the dimensionless dynamic coefficient of this model and ∆ the filter width (equal to the characteristic grid size in this work). According to the Pope criterion, a good LES should be able to resolve at least 80% of the turbulent kinetic energy. Looking at the phase where the turbulent activity is the highest (t/T =0.65), it was found that computation captures more than 80% of the turbulent kinetic energy in 85% of the numerical domain. The last 15% are mainly located on the atrium surface and in the atrial cavity. Sketch of the inflow flow rate presented in Fig. 5 is reported with a time indication. The mapping relating vorticity magnitude with opacity is linear.

Heart wall is made partially transparent to allow observation of the flow behaviour. At the same time, the LA volume increases and is filled with blood. Four jets are present, three of them colliding head on, while the lower part of the LA remains vorticity-free. Blood coming from the RSPV (on the left side of the images) follows a trajectory tangential to the atrial wall giving rise to a net swirling motion in the atrium as observed in vivo [START_REF] Fyrenius | Three dimensional flow in the human left atrium[END_REF] and in experimental work [START_REF] Tanné | Assessment of left heart and pulmonary circulation flow dynamics by a new pulsed mock circulatory system[END_REF]. High vorticity is visible at the pulmonary veins roots due to their narrowed geometry. At the end of systole, the vorticity begins to decrease in the entire heart, as shown at t/T =0.35 (top center plot in Fig. 6). At the end of the ventricle contraction, the aortic valve closes and the mitral valve opens: ventricle filling starts.

Results and discussion

Global description of the flow

The ventricle diastole starts at t/T =0.375 (not shown): the LV volume increases and blood passes from the LA to the LV, forming a strong jet through the MV. The shear layer between the jet generated and the surrounding quiescent fluid rolls-up and shapes the jet head as a vortex ring [START_REF] Le | On the three-dimensional vortical structure of early diastolic flow in a patient-specific left ventricle[END_REF]. The top right plot in Fig. 6 depicts this vortex ring at t/T =0.45. This mechanism was reported and studied by several authors, notably by Domenichini et al. [START_REF] Domenichini | Three-dimensional filling flow into a model left ventricle[END_REF]. A similar process is visible in the LA where formation of four vortex rings takes place at the four pulmonary veins (only two are visible in the top right image in Fig. 6).

At the E wave peak (bottom left image in Fig. 6, t/T =0.55), three of the LA jets collide head on, as described before, but in a more intense way, as the pulmonary veins flow rate is higher at the E wave peak than at t/T =0.25.

Again, Fig. 6 shows that blood coming from the RSPV follows clearly a trajectory tangential to the atrial wall. As the E wave head jet is getting closer to the wall with an angle of approximately 50 degrees, a shear layer is generated at the lateral heart wall (the right one on plots). The part of the vortex ring that is closer to the lateral wall interacts with it and dissipates.

The other part of the vortex ring is marginally affected and thus remains almost intact, moving towards the apex of the ventricle, as described in the literature [START_REF] Mihalef | Patient-specific modelling of whole heart anatomy, dynamics and haemodynamics from four-dimensional cardiac CT images[END_REF].

As expected for a flow at such a Reynolds number [START_REF] Cheng | Numerical study of a vortex ring impacting a flat wall[END_REF], small-scale vortices are also generated with the vortex-wall interaction at t/T =0.65 (Fig. 6, bottom center plot). At the end of the E wave, numerous small vortices are present in the whole LV.

As inlet flow decreases, the vorticity magnitude decreases as well in the entire heart. The LV volume remains stable during this phase called diastasis. Vorticity magnitude progressively decreases in all the left heart due to dissipation of the vortices by viscosity. This is confirmed by the vorticity levels observed in Fig. 6 at t/T =0.99

(bottom right image). Vorticity is now small everywhere, except near the mitral valve, due to the small A wave. Note that contrary to what is seen during the E wave, blood is not entering the LA during the A wave. Contraction of the LA to finish the LV filling results in a small outward flow through the pulmonary veins resulting in non-zero vorticity values.

Velocity fields

The velocity vector fields of the averaged flow solutions are shown in Fig. 7, for the same instants shown in Fig. 6. The non-dimensional velocity vector fields u f /u a are shown over the slice through the left heart indicated in Fig. 3. The reference velocity u a is computed as

u a = ql s /V s = 0.1 m.s -1
where q is the cardiac output ( q=7.50 × 10 -5 m 3 .s -1 ), V s the end systolic volume (V s =5.55 × 10 -5 m 3 ) and l s is the ventricle length at the end of the systole (l s =7.40 × 10 -2 m). Note that due to the strong velocity variations along the cycle, the vector scale is adapted for the different instants. Again, the first two images of Fig. 7 correspond to systole. At the beginning of systole, high velocity is mainly observed in the ascending aortic root as shown in Fig. 7 left top image. Towards the end of systole at t/T =0.35 (top center image in Fig. 7), just before the opening of the MV, blood entering the LA is already directed towards the LV. In the lower half part of the LA, coexistence of an axial movement toward the MV and a rotational movement is observed.

This helical flow has been reported in-vivo by Kilner et al. [START_REF] Kilner | Asymmetric redirection of flow through the heart[END_REF] and Markl et al. [START_REF] Markl | Comprehensive 4D velocity mapping of the heart and great vessels by cardiovascular magnetic resonance[END_REF]. In order to see this flow feature hardly discernible in Fig. 7, Fig. 8 (left plot) shows the velocity field projected on a plane perpendicular to the In the ventricle, the contraction resulting in blood ejection through the aorta generates a small recirculation under the MV (on the right of the plot).

Another recirculation area is noticeable in the aorta: a normal anatomic feature of the ascending aorta is a dilatation of the vessel just above the AV. During the whole systole, blood recirculation occurs in it [START_REF] Kilner | Helical and retrograde secondary flow patterns in the aortic arch studied by three-directional magnetic resonance velocity mapping[END_REF], which can be observed clearly in the velocity field shown at t/T =0.35. However, the aortic valve leaflets dynamics are not modelled here when the valve is open. Consequently, the flow dynamics observed in that region may change if a more realistic aortic valve model were used.

At the beginning of diastole (top right image in Fig. 7), blood enters the LA and the LA contracts, resulting in the E wave through the MV. In the LA, blood is clearly directed towards the MV [START_REF] Fyrenius | Three dimensional flow in the human left atrium[END_REF]. The E wave vortex ring signature can also be seen in the right top plot in Fig. 7 at t/T =0.45 and its evolution can be followed to the left bottom plot, at t/T =0.55. There, the vortex ring is no more symmetric, as the lateral wall prevents its full development. On the other side, the vortex ring gains strength, as described in the previous section. At the same time, the atrial flow seems separated in two halves again. Blood in the lower half is directed to the MV while in the upper half, the flow is a more chaotic due to the collision of the pulmonary veins blood streams.

At t/T =0.65 (bottom, center plot in Fig. 7), a large recirculating cell is visible in the LV, as described classically in the literature [START_REF] Hong | Characterization and quantification of vortex flow in the human left ventricle by contrast echocardiography using vector particle image velocimetry[END_REF][START_REF] Markl | Comprehensive 4D velocity mapping of the heart and great vessels by cardiovascular magnetic resonance[END_REF][START_REF] Long | Subject-specific computational simulation of left ventricular flow based on magnetic resonance imaging[END_REF]. It Sketch of the inflow flow rate presented in Fig. 5 is reported with a time indication.

is characteristic of the flow in the ventricle after the E wave. Note that the intense upward motion along the septum wall is responsible for the high values of vorticity which can be seen in Fig. 6 at t/T =0.65. Two less intense blood recirculation zones can be detected: one at the apex, which is visible during the whole diastole and an intermittent one between the aortic valve and one of the MV leaflet. These blood recirculations are also described in-silico by Mihalef et al. [START_REF] Mihalef | Patient-specific modelling of whole heart anatomy, dynamics and haemodynamics from four-dimensional cardiac CT images[END_REF], Doenst et al. [START_REF] Doenst | Fluid-dynamic modeling of the human left ventricle: Methodology and application to surgical ventricular reconstruction[END_REF] and Schenkel et al. [START_REF] Schenkel | MRIbased CFD analysis of flow in a human left ventricle: Methodology and application to a healthy heart[END_REF].

Figure 8 (right image) shows the averaged non-dimensional velocity field on another plane at t/T =0.65. As at t/T =0.35, in the lower part of the atrium, the downward axial movement towards the MV is rather low compared to the rotational movement: Fig. 8 shows a net swirling motion, visible not only above the MV (right top plot) but also at the tip of the leaflets (right bottom plot). The swirling motion at the MV is visible during all the remainder of diastole, supporting the fact that the velocity profile at the MV is not only skewed [START_REF] Haugen | Velocity profiles in mitral blood flow based on three-dimensional freehand colour flow imaging acquired at high frame rate[END_REF][START_REF] Dahl | Impact of pulmonary venous locations on the intra-atrial flow and the mitral valve plane velocity profile[END_REF] but has also a non-negligible rotational component. This flow feature can actually be observed at each of the decelerations at the pulmonary veins inflow: a swirling motion in the lowest part of the atrium maintains blood motion in the LA even without net incoming flow through the pulmonary veins. Figure 8 (right image) shows the phase-averaged nondimensional velocity magnitude on the same planes. Velocity fields show that for this heart simulation, a skewed velocity profile is obtained. Both the order of magnitude and the structure of the velocity fields show similarities with the ones presented by Dahl and al [START_REF] Dahl | Impact of pulmonary venous locations on the intra-atrial flow and the mitral valve plane velocity profile[END_REF]. Note however that the velocity field at the MV may depend on the chosen flow distribution at the pulmonary veins.

Between the E wave and the A wave, the recirculating cell core in the LV moves from the ventricle center to the septum wall. During the A wave (bottom, right plot in Fig. 7) occurring at t/T =0.99, the blood flux passing though the MV strengthens the recirculating cell in the LV, as classically reported. Atrial contraction expels blood from the LA, both through the MV, as seen in the lower half part of the LA and through the pulmonary veins, as shown by the upward velocity vectors visible in the upper half part of the LA.

Instantaneous structures and cycle-to-cycle variations 4.3.1. Instantaneous structures

Instantaneous flow features have been ignored in the previous sections. In order to describe them, the Q criterion, first introduced by Hunt et al. [START_REF] Hunt | Eddies, stream, and convergence zones in turbulent flows[END_REF],

is displayed. This criterion uses the second invariant of the velocity gradient tensor,

Q = 1 2 (Ω 2 -S 2 ), (26)
where Ω and S are respectively the antisymmetric and the symmetric components of the velocity gradient tensor. Thus, a positive Q criterion indicates a higher rotation rate in comparison to the strain rate, allowing to detect vortical structures. Images of Q criterion isosurfaces are presented in Fig. 9 at different instants of the heart cycle. The instants are the same as in the previous sections and the isosurface at Q = 30u 2 a /R 2 a is shown. At the beginning of systole (Fig. 9 In the course of the following systole, such structures are convected and smoothly elongated towards the aorta. During this process, a unique vortical structure of length up to 10 R a is formed as shown in Fig. 9 top left and center plots. Pieces of this structure remain in the ventricle even after the aortic valve is closed. Figure 9 shows the remnant structure under the aortic valve even after the beginning of the E wave at t/T =0.55. The LV remains almost free from vortical structures during systole, as the lower half part of the LA. However, the upper half of the LA shows many structures and swirling flow is noticeable thanks to the Q criterion at the pulmonary veins on the top left plot in Fig. 9.

At t/T =0.45 (top right plot in Fig. 9), vortex rings are created at the ostia of the pulmonary veins, due to the unsteady flow entering the LA. Two of the vortices are clearly visible in the figure . The vortex ring at the MV described before is also visible. The MV vortex ring travels about two radius R a in direction of the apex during a time interval of T /10 before it hits the lateral wall. The Q criterion isosurfaces in the left and center bottom images in Fig. 9 show the evolution of the vortical structures in the LV from the initial ring-like shape to the final complex 3D flow, when the ventricle cavity is almost completely filled by vortical structures at t/T =0.65.

During the late ventricle filling, the atrium contracts and generates a small vortex ring (right bottom plot in Fig. 9). However, its intensity is much lower than the one generated during the E wave, and it is hardly visible.

The vortex head only travels about one radius R a within a time span of T /12 which corresponds to a velocity 60% smaller than the vortex ring created by the E wave. Finally, vortical structures dissipate almost everywhere in the heart. It has to be noticed that the swirling flows visible at the pulmonary veins at t/T =0.25 and t/T =0.65 through the Q criterion are also visible at t/T =0.99: the pulmonary veins geometry seems to generates this swirling motion.

Velocity fluctuations

The studied flow configuration is a breeding ground for weak turbulence.

The transitional nature of this cyclic flow due to high Reynolds numbers and unsteady inflows results in cycle-to-cycle variations. Investigation of these variations can be quantified through the non-dimensional kinetic energy of the velocity fluctuations. This fluctuating kinetic energy (FKE) is defined as,

E k = 1 2 (u 2 rms + v 2 rms + w 2 rms ), (27)
where u rms , v rms and w rms are the root mean square values of the velocity fluctuations in the three directions. Note that in the present case, E k measures not only the turbulent activity, but more generally cycle-to-cycle variations. Hence, the choice of naming this quantity FKE instead of turbulent kinetic energy was made.

Figure 10 shows the spatial distribution of the FKE along the heart cycle, over the same cutting plane as in Fig. 7.The beginning of systole is associated with a relaminarization process in the LV. FKE intensity levels across the ventricle during systole (top left and center plots in Fig. 10) stay almost null.

The highest FKE levels are present in the root of the aorta between t/T =0.25 and 0.35. Note however that the FKE levels downstream of the AV may be mispredicted here, due to the absence of AV leaflets modelling when the AV 10) until the jet hits the wall, the more intense FKE takes place in the proximity of the jet. These fluctuations are mainly due to small differences in the location of the shear layer and the vortex ring, where the velocity gradients are high. At the E wave peak at t/T =0.55, depicted in the bottom left plot in Fig. 10, high FKE (about 5 u 2 a ) is observed in the jet edge which becomes increasingly unstable throughout the onset of the MV flow deceleration. Furthermore, high FKE of the same order is visible at the lateral wall.

As shown by the FKE field at t/T =0.55, the impact of the vortex ring on the lateral wall is a zone of high velocity fluctuations. At the same time, the E wave starts to decrease. At t/T =0.65 (bottom center image in Fig. 10), the flow through the mitral valve is approximately null. However, the flow in the ventricle is far from quiescent. The mitral jet impact has made the initial vortex ring shatter into small vortices, that occupy the main part of the ventricle. In accordance with this transition to (weak) turbulence, the FKE field shows high values over almost the whole ventricle. Largest values of FKE (from 6 to 10 u 2 a ) are observed at the core of the recirculating cell. These high values of FKE are due to the differences in the position of recirculating cell center between cycles, as reported in this kind of geometrical configuration [START_REF] Druault | Use of proper orthogonal decomposition for time interpolation from PIV data[END_REF]. The FKE in the LA also indicates that the pulmonary jets impact in the upper part of the LA is associated with intense velocity fluctuations from one cycle to another. FKE is about three time lower in the bottom half of the LA. Such a ratio of FKE between the upper and bottom halves of the LA remains approximately the same during the entire heart cycle.

At the end of diastole at t/T =0.99 (bottom right plot in Fig. 10), the FKE becomes more homogeneous in the heart, albeit with higher values at the center of the ventricle recirculating cell (see also Fig. 7) and at the MV tips due to the A wave vortex ring. As soon as systole begins, a relaminarization process is engaged in the ventricle: blood acceleration reorganizes the blood flow and destroys all the remaining vortices. Comparing the FKE fields at t/T =0.99 and at t/T =0.25 shows that late diastole velocity fluctuations do not persist in systole.

Conclusion and outlook

In this paper, a numerical framework and its application to a patient left heart has been exposed. This numerical framework only requires gated 4D images of a patient heart and global morphological parameters of the mitral valve as input data for the computation. A numerical domain is extracted from one 3D medical image and the heart wall movements are automatically calculated to follow the heart wall displacements recorded in the 4D medical images, thanks to a non-linear image registration algorithm. Then, temporal interpolation of the geometry is used to generate a finite-volume mesh of the left heart at any time over the cardiac cycle.

A computation of the blood flow in a full left heart including the left ventricle, atrium and the aorta has been conducted in this domain. The present numerical method is able to provide results consistent with the cur-rent knowledge in terms of left heart flow. The E wave and the A wave are recovered, as well as the pulmonary veins flow reversal associated with the late atrial contraction. In the left atrium, the flow reported shows a clear swirling motion, that maintains blood in motion even without net blood flow coming from the pulmonary veins. Furthermore, the well-known mitral jet, preceded by an energetic vortex ring, is also observed. The classical large recirculating cell, characteristic of the flow during diastasis, can also be seen in the phased-averaged velocity fields. All these features have been reported several times in the literature, both in numerical and experimental studies and by medical imaging. are features tending to promote turbulence. We show that the late diastole large recirculation cell, which is classically reported (and found here in the phase-averaged field) is most probably perturbed by small vortices.

The calculation reported here has several limitations. First, blood has been modelled as a Newtonian fluid. Even if this approximation is commonly accepted for the heart flow, a non-Newtonian model could be included in the present method to assess the impact of the rheological model on the flow.

Then, the spatio-temporal resolution of the medical images is an important limitation, which imposes temporal interpolation and geometrical simplifications of the heart model. It also makes it difficult to accurately account for short phases as the isovolumic contraction and relaxation. Although ten images per cardiac cycle enable a convincing flow prediction, specific tests should be performed to precisely assess the impact of the temporal resolution of the medical images on the flow field. Another limitation is the unknown flow distribution between the pulmonary veins, which has to be prescribed in the simulations.

In the present method, the main element to improve is certainly the valve modelling. The model used for valves highly depends on the anatomical information that can be extracted from the images. Consistently with their insufficient spatio-temporal resolution, a rough model of the valves was used. This is true from the anatomical and numerical points of view. Better valve models should notably include valve opening and closing, change of aperture area along time. Aortic valve leaflets were notably completely omitted during systole, which most probably impacts the flow at the beginning of the aorta.

In addition, the immersed forcing used in this work could be improved [START_REF] Fadlun | Combined immersedboundary finite-difference methods for three-dimensional complex flow simulations[END_REF], notably by accounting for the exact valve velocity.

Still, the presented method does include important aspects of the physiological heart: the entire ventricle and atrium, the aorta root, the pulmonary veins and the valves. Such a geometrical complexity has rarely been achieved.

As a consequence, rarely or never reported features such as secondary vortex recirculation, vortical structure at the aorta, swirling motion at the MV, vortex rings, swirling motion at the pulmonary veins and velocity fluctuations have been discussed for the first time in details. Furthermore, the presented results underline the fact that over-simplified geometries cannot fully reproduced certain physical events as for example the mitral in-plane swirling motion.

To conclude in a clinical point of view, the presented framework presents significant potential to link patient-specific heart geometry to pathologies.

Interpretation of altered cardiac mechanics through blood flow in large cohort is possible, in particularly in diseases such as hypertrophy and heart failure.

The

  LA and the LV are separated by the mitral valve (MV) which prevents blood from flowing back from the LV to the LA. The MV has two leaflets which are prevented from prolapsing by the chordae tendineae and papillary muscles running from the cusps of the valve leaflets to the side of the LV. The LV is approximately conical in shape and presents a nearly circular outline in its transverse section. The LV and the AO are separated by the aortic valve (AV) positioned at a stenosis between the two cavities, inserted into the walls of the sinuses of Valsalva at the root of the aorta. The aortic valve has three simple leaflets without any attachments that come together, providing mutual support when they are closed. Two periods can be distinguished during the cardiac cycle: the diastolic period (the AV is closed, the MV open) and the systolic one (the AV is open, the MV closed). During diastole, the oxygenated blood comes from the LA to the LV. During systole, the myocardium contracts and blood is ejected from the ventricle. During these two phases, the LA and LV walls have inhomogeneous and rapid movements with large volume variations (up to 60% for the LV).

2 .Figure 2 :

 22 Figure 2: Top left: annulus markers and modelled MV leaflets. Parameters a, b of the ellipse ǫ are indicated. Top right: full MV model included in the numerical domain (the AV is represented as well). Bottom: MV scheme.

  Inlet boundary conditions can be extracted from PC-MRI and interpolated on the computational domain inlets surfaces. However, if only CT scan images are available, a different strategy must be used. As blood is considered incompressible, reasoning based on mass conservation can be made to overcome the lack of inflow and outflow information. As already stated in the last section, only two topological configurations are considered (MV closed, AV open or AV closed, MV open). This assumption allows a consistent definition of the inlet/outlet boundary conditions. Indeed, mass conservation imposes, for the first case (MV closed, AV open):

Figure 3 : 3 .

 33 Figure 3: Template computational domain extracted from a 3D medical image. The same domain is represented for four different points of view and the left ventricle (LV), left atrium (LA) and Aorta (AO) are indicated. Black line passing through the left heart indicates the position of slices used to describe the flow in section 4.

Figure 4 :

 4 Figure 4: Left, left atrium with its four pulmonary veins named. Right, sketch of inflow angles and the four pulmonary veins projected on a plane which can be seen on the left picture. The dashed line indicates that the flow is directed a lot more toward the MV compared to the others.

  3). A close examination of the medical images from the CT scan allowed to set the leaflets length to l = 12 mm for the MV. The open area presented to the blood flow is represented by an ellipse of axis a = 15 mm and b = 8 mm. This area is supposed constant over the time when the MV is open. Eccentricity is fixed to e = 5 mm from the same medical data.

  The fluid density ρ and the fluid kinematic viscosity ν are supposed to be constant: ρ = 1040 kg.m -3 and ν = 4 × 10 -6 m 2 .s -1 . Table 2 reports the maximum velocity magnitude and the maximum Reynolds number for different elements of the heart model: pulmonary veins, mitral valve and aortic valve. Inlet Reynolds number for each pulmonary vein varies from 0 to approximately 2000, based on the vein diameter where the flow is imposed. The Reynolds number at the mitral tips varies from 0 to approximately 5000 (based on the effective mitral mean diameter D = 2R a = 2 √ ab, the area of the open mitral, the kinematic viscosity and the maximum flow rate). The maximum transmitral velocity U max falls into the usual measurements

Figure 6 :

 6 Figure 6: Volume rendering of non-dimensional vorticity magnitude of the phase-averaged velocity fields at different times of the simulation. First row, left plot: t/T =0.25. Center plot: t/T =0.35. Right plot: t/T =0.45. Second row, left plot: t/T =0.55. Center plot: t/T =0.65. Right plot: t/T =0.99.

Figure 6

 6 Figure6shows the 3D vorticity magnitude map of the phase-averaged velocity field, nondimensionalised by the period T . Six salient instants of the heart cycle are shown: the ventricular mid-systole (t/T =0.25), the end of the ventricular systole (t/T =0.35), the beginning of the E wave (t/T =0.45), its peak (t/T =0.55), its end (t/T =0.65) and the end of the A wave, just before the beginning of the next ventricular systole (t/T =0.99).In the first row of Fig.6, the first two images correspond to systole. The mitral valve is closed, while the aortic valve is open. The ventricle volume is conspicuously decreasing, causing ejection of the blood into the aorta. The narrowed geometry of the aorta root generates high vorticity through shear layer clearly shown at t/T =0.25, while the ventricle remains free of vorticity.

Figure 7 :

 7 Figure 7: Average non-dimensional velocity field (u f /u a ) projected on a slice of the heart. Velocity vector scale is not constant though the heart cycle and is indicated for each plot. First row, left plot: t/T =0.25. Center plot: t/T =0.35. Right plot: t/T =0.45. Second row, left plot: t/T =0.55. Center plot: t/T =0.65. Right plot: t/T =0.99. Sketch of the inflow flow rate presented in Fig. 5 is reported with a time indication.

Figure 8 :

 8 Figure 8: Phase-averaged non-dimensional velocity (u f /u a ) vectors field projected on planes and velocity magnitude. Planes are indicated on the heart sketches and are viewed from the top of the heart in the apex direction (plane's orientation is indicated). The MV leaflets are visible in the right bottom plots and are coloured in gray. Left: t/T =0.35. Right: t/T =0.65.

  Fig.9).

Figure 9 :

 9 Figure 9: Evolution of Q criterion isosurface (Q = 30u 2 a /R 2 a ) showing structures present in the instantaneous field though the heart cycle. First row, left plot: t/T =0.25. Center plot: t/T =0.35. Right plot: t/T =0.45. Second row, left plot: t/T =0.55. Center plot: t/T =0.65. Right plot: t/T =0.99. Sketch of the inflow flow rate presented in Fig. 5 is reported with a time indication.

Figure 10 :

 10 Figure 10: Nondimensional fluctuating kinetic energy E k /u 2 a . First row, left plot: t/T =0.25. Center plot: t/T =0.35. Right plot: t/T =0.45. Second row, left plot: t/T =0.55. Center plot: t/T =0.65. Right plot: t/T =0.99. Sketch of the inflow flow rate presented in Fig. 5 is reported with a time indication.

  Furthermore, the use of fluid numerical method well adapted to turbulent flows enables the observation of cycle-to-cycle variations in the flow field. Such variations are expected in the present flow, due to the high Reynolds numbers encountered and the unsteadiness of the flow incoming by the pulmonary veins. The present results show that in spite of rigorously identical contraction and boundary conditions, fluid inertia makes the flow differ from one cycle to another. More precisely, cycle-to-cycle variations in the left atrium can be observed in its upper part, where the collision of the jets issuing from the pulmonary veins makes the flow particularly chaotic. Cycleto-cycle variations are high after each peak in the pulmonary vein flow rate, during flow deceleration. In the left ventricle, velocity fluctuations are reported mainly in late diastole. Between the impact of the E wave jet on the lateral wall and the end of diastole, the left ventricle displays high levels of cycle-to-cycle fluctuations. Indeed, both the vortex ring impact and the E wave deceleration occur approximately at the same time, and both

Table 1 :

 1 Geometric characteristics of the ostia. Ostia are oval in shape. The areas reported are obtained by assuming elliptic shapes.

	RIPV LSPV LIPV

π

  Re max t m /T

	AV	0.96 m.s -1	2.20 cm	5300	0.16
	MV	1.08 m.s -1	1.87 cm	5000	0.52
	Pulm. vein 0.79 m.s -1	1.00 cm	2000	0.52

Table 2 :

 2 Main flow parameters describing the simulation. The section-averaged maximum velocity is indicated as U max . The maximum Reynolds number Re max is calculated thanks to mean D the diameter. Time when the Re max is reached is reported as t m /T . For valves, A is the area of the lumen when open.

Acknowledgments

The authors would like to express their gratitude to Dr. R. Moreno, Dept. of Radiology, for the CT scan images from the Rangueil University Hospital, Toulouse, France and to Dr. S. Nottin, MD Dr. D. Coisne, MD Dr. I. Schuster for all the fruitful discussions. Dr. V. Moureau is gratefully acknowledged for giving access to the YALES2 solver.