F Nicoud

C Chnafa
email: cchnafa@mie.utoronto.ca

J Siguenza

V Zmijanovic

S Mendez

⋅ V Zmijanovic

⋅ S Mendez

Large-Eddy Simulation of Turbulence in Cardiovascular Flows

A 4th-order accurate, low dissipative flow solver is used to perform Large-Eddy Simulations of three typical hemodynamic situations: the flow through the idealized medical device proposed by the American Food and Drug Administration; the intracardiac flow within an actual human left heart whose morphology and deformations are deduced from medical imaging; the flow downstream of an artificial aortic valve which arises from the blood-leaflets interaction problem. In all the cases, the subgrid scale model designed to handle wall-bounded transitional flows is successfully used and the numerical simulations compare favourably with the experimental data available. These results illustrate the potential of the Large-Eddy Simulation methodology to properly handle blood flows. They also support the idea that turbulence, even if not fully developed, may be present in cardiovascular flows, including under non pathological conditions.

Introduction

Turbulence is a natural phenomenon in fluid mechanics, which gives rise to wide spectrum, chaotic fluctuations in the velocity/pressure fields, even in the absence of external forcing. This phenomenon arises when the inertia forces are significantly larger than the viscous effects, viz. when the Reynolds number R e exceeds a few thousands in steady flow, most probably less for pulsatile flows [START_REF] Varghese | Direct numerical simulation of stenotic flows. Part 2. Pulsatile flow[END_REF]. This number usually reaches a few hundreds in medium-sized arteries and a few thousands in the largest vessels, in the heart itself or within biomedical systems such as total artificial hearts or ventricular assist devices. Thus, one may expect turbulence to be present in the cardiovascular system, at least in some favourable regions in space (e.g. close to stenoses, valves, bifurcations) and during specific periods of the cardiac cycle (when the flow decelerates). Still, until recently turbulence was virtually not studied, at least non-invasively, in the cardiovascular biomechanics community, notably because in vivo 3D and unsteady measurements of the blood velocity are very challenging. Recent progresses in both computational fluid dynamics (CFD) and medical imaging techniques give rise to a renewed interest for the turbulent properties of blood flows. Unstable and/or turbulent flows have been observed in numerical simulations in abdominal aortic aneurysms [START_REF] Les | Quantification of hemodynamics in abdominal aortic aneurysms during rest and exercise using magnetic resonance imaging and computational fluid dynamics[END_REF], carotid siphons [START_REF] Valen-Sendstad | High-resolution computational fluid dynamics detects flow instabilities in the carotid siphon: implications for aneurysm initiation and rupture?[END_REF], cerebral aneurysms [START_REF] Mikhal | Immersed boundary method for pulsatile transitional flow in realistic cerebral aneurysms[END_REF] as well as experimentally in an idealized left ventricle [START_REF] Domenichini | Combined experimental and numerical analysis of the flow structure into the left ventricle[END_REF]. Thanks to an MRI method based on the intravoxel velocity standard deviation, it is now possible to estimate the turbulent kinetic energy in vivo. Significant levels of kinetic energy were found not only in pathological left atria [START_REF] Dyverfeldt | Hemodynamic aspects of mitral regurgitation assessed by generalized phase-contrast MRI[END_REF] but also in normal left ventricles [START_REF] Zajac | Turbulent kinetic energy in normal and myopathic left ventricles[END_REF], supporting the idea that turbulence can also occur in healthy human subjects.

From a numerical point of view, representing turbulence is a challenge due to the wide variety of both temporal and spatial scales that are present at the same time. In the ideal case of a decaying isotropic turbulence, one may show that the largestto-smallest length scales ratio behaves like R 3∕4 e [START_REF] Frisch | Turbulence: The Legacy of A.N. Kolmogorov[END_REF]; as a consequence, the number of grid points necessary to properly represent a turbulent flow, including the whole range of fluctuations, increases as fast as R 9∕4 e . This approach where turbulence is accounted for by just and only solving all the scales produced by the Navier-Stokes equations is called Direct Numerical Simulation and remains most of the times out of reach of the current computing power except in simple cases with moderate Reynolds number. The large variety of scales inherent to any turbulent flow motivated the development of turbulent models in order to complement the Reynolds-Averaged Navier-Stokes (RANS) equations which allow to compute the averaged flow quantities without explicitly representing the fluctuations. Despite decades of very active research and some successes in fully developped turbulent flows, there are still many controversies about the proper way to develop 'good' turbulent models [START_REF] Spalart | Progress in aerospace sciences philosophies and fallacies in turbulence modeling[END_REF] and the idea of a universal formulation able to properly deal with any turbulent flows is widely considered as hopeless. On top of not being predictive enough, RANS models are better suited for fully turbulent flows at very high Reynolds numbers than transitional flows. The ability of Large-Eddy Simulation (LES) to better predict a flow without prior knowledge, including regarding its laminar/turbulent nature, is one of the reasons why this technique is now being used in many applications. It is a method of choice for unsteady 3D flows at moderate Reynolds numbers, especially when transition from laminar to turbulence plays a key role. In this paper, we illustrate how LES can be used to represent three configurations typical of cardiovascular applications:

• an idealized medical device where a stationary flow goes from laminar to turbulent in a simple geometry [START_REF] Hariharan | Multilaboratory particle image velocimetry analysis of the FDA benchmark nozzle model to support validation of computational fluid dynamics simulations[END_REF],

• a realistic human left heart whose geometry and cyclic deformations are deduced from medical imaging and where CFD/LES in a time varying domain is used to generate the associated blood flow [START_REF] Chnafa | Image-based large-eddy simulation in a realistic left heart[END_REF],

franck.nicoud@umontpellier.fr

• the flow past an artificial aortic valve where the Fluid-Structure Interaction problem is solved to predict the leaflets dynamics and associated turbulent generation.

Some numerical and physical issues specific to each of these configurations will be discussed in Sects. 3, 4 and 5 respectively, together with comparison with experimental data when available. Some details about the numerics and the subgrid scale modelling strategy are first given in Sect. 2.

LES Requirements

In LES, the modelling effort focuses on the smallest length scales of fluctuations which are presumably less dependent on the geometry of the flow domain and thus more universal. The largest scales (which can be as small as twice the typical cell size, at least in theory) are explicitly resolved as solutions of the low-pass filtered Navier-Stokes equations [START_REF] Sagaut | Large Eddy Simulation for Incompressible Flows. An Introduction[END_REF]. Thus, the finer the mesh, the larger the range of resolved scales and the smallest the contribution of the subgrid scale model. In other words, the modelling effort in less intense in LES compared to RANS and this is all the more true that the grid is refined.

The price to pay is that the dynamics of the largest scales, as well as the related macroscopic mixing and turbulent stresses, must be properly represented by the numerics. As a consequence the numerical errors associated to both the spatial and temporal discretizations must be kept as small as possible. Notably, dissipation-free schemes must be preferred so that the (relative) intensity of the different resolved scales is not spoiled. An important feature of LES is that the scheme should be as accurate as possible for all the scales represented on the mesh, the largest as well as the smallest. This means that the order of accuracy of a numerical method is not a good measure of its capability to perform well in the LES; the effective wave number [START_REF] Vichnevetsky | Fourier Analysis of Numerical Approximations of Hyperbolic Equations[END_REF] and kinetic energy conservation [START_REF] Morinishi | Fully conservative higher order finite difference schemes for incompressible flow[END_REF] concepts are better suited to characterize a numerical scheme dedicated to LES. For example, better results may be obtained with a second order centered (dissipation-free at all scales) scheme than with a higher-order upwind biased (dissipative) scheme [START_REF] Mittal | Suitability of upwind-biased finite differance schemes for large eddy simulation of turbulent flows[END_REF][START_REF] Park | Discretization errors in large eddy simulation: on the suitability of centered and upwind-biased compact difference schemes[END_REF]. As a consequence, a RANS solver, usually based on stabilized scheme coupled to implicit time marching methods designed to speed up the convergence towards a steady solution cannot be easily adapted to LES. In the present study, the in-house flow solver YALES2BIO (http://www.math.univ-montp2.fr/~yales2bio) was used to solve the filtered flow equations. Dedicated to the computations of blood flows, this general purpose solver relies on a 4th-order accurate centred (dissipation-free) finite-volume formulation where the projection method [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF] is used to meet the divergence-free condition. The time-stepping is also 4th-order accurate, based on the dissipationfree explicit low-storage Runge-Kutta scheme [START_REF] Williamson | Low-storage Runge-Kutta schemes[END_REF]. At the end, the algorithm is non dissipative and one relies on the physical dissipation, either laminar or issued by the subgrid scale model, to ensure numerical stability. This is made possible because the Reynolds number is low to moderate in biomedical applications so that numerfranck.nicoud@umontpellier.fr ical stabilization is not necessary; this also requires to use high quality 3D meshes. Finally, YALES2BIO inherits its parallel efficiency from the YALES2 package developed by V. Moureau and co-workers [START_REF] Moureau | Design of a massively parallel CFD code for complex geometries[END_REF][START_REF] Malandain | Optimization of the deflated Conjugate Gradient algorithm for the solving of elliptic equations on massively parallel machines[END_REF].

From the modelling point of view, the subgrid scale stress tensor SGS which results from the filtering operation applied to the Navier-Stokes equations is most of the time modelled thanks to the eddy-viscosity assumption (introduced in the RANS concept in the 70s). In the case of an incompressible fluid, this reduces to

SGS = 2 SGS
where is the fluid density, is the strain rate tensor of the resolved scales and SGS is the so-called subgrid scale viscosity. Several models were proposed over the years for this quantity and most of them share the following form:

SGS = (C ) 2 D( ) (1)
where C is the model constant which is usually tuned (either theoretically or numerically) so that the model produces the proper amount of dissipation in the simple case of decaying isotropic turbulence. The length scale denotes the typical size of the local cell of the mesh used to solve the filtered Navier-Stokes equations and D is a differential operator which defines the model and operates on the resolved velocity field . The -model [START_REF] Nicoud | Using singular values to build a subgridscale model for large eddy simulations[END_REF] is used in this study, meaning that:

D(u) = 3 ( 1 -2 )( 2 -3 ) 2 1
(2)

In this expression, 1 ≥ 2 ≥ 3 ≥ 0 are the three singular values of the local velocity gradient tensor and can be efficiently computed [START_REF] Nicoud | Using singular values to build a subgridscale model for large eddy simulations[END_REF]; the model constant is C = 1.35. This model was selected because it meets several useful properties relevant in terms of SGS modelling, although not shared by the other SGS models. On top of vanishing for a variety of canonical laminar flows for which no SGS viscosity is expected, it has the proper cubic behavior in near-wall regions. In [START_REF] Nicoud | Using singular values to build a subgridscale model for large eddy simulations[END_REF], the model was implemented in different numerical solvers and validated in academic cases (homogeneous isotropic turbulence, turbulent channel); its effectiveness in a more demanding configuration was also demonstrated recently by considering the impingement of an unsteady jet with a rigid wall [START_REF] Baya Toda | Assessment of subgrid-scale models with a large-eddy simulation-dedicated experimental database: the pulsatile impinging jet in turbulent cross-flow[END_REF]. In this particular case which is relevant to intra-cardiac hemodynamics, the comparison with experimental data shows that the -model outperforms the well-known Dynamic Smagorinsky model [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF].

In the next three sections, results from the -model as implemented in the YALES2BIO solver are discussed to illustrate how LES can be used to analyse flows relevant to cardiovascular applications.

franck.nicoud@umontpellier.fr

The FDA Medical Device Test Case

The United States Food and Drug Administration (FDA) medical device benchmark was introduced at the beginning of the current decade as an answer to the ubiquitous presence of various numerical methods in biomedicine. The purpose of the proposed benchmark program was to validate the CFD codes on idealized medical device examples. The first FDA test model [START_REF] Hariharan | Multilaboratory particle image velocimetry analysis of the FDA benchmark nozzle model to support validation of computational fluid dynamics simulations[END_REF][START_REF] Stewart | Assessment of CFD performance in simulations of an idealized medical device: results of FDA's first computational interlaboratory study[END_REF] aims at challenging CFD codes on a case whith transition to turbulence, which is frequently encountered in biomedical devices. Accurately predicting turbulence is necessary to properly assess shear stresses, themselves relevant to blood damage prediction in ventricular assist devices (VAD) [START_REF] Marsden | Recent advances in computational methodology for simulation of mechanical circulatory assist devices[END_REF][START_REF] Mendez | Assessing macroscopic models for hemolysis from fully resolved simulations[END_REF], but also in blood flows through stents and pumps, to cite a few. Mispredicting the turbulent nature of blood flows heavily impacts the results and may cause critical faulty conclusions.

The FDA medical device model, as illustrated in Fig. 1, consists of a fixed axisymmetric geometry with a long inlet section, a convergent nozzle, a 10 diameters long throat section and a sudden expansion at the downstream end. This geometry submits the flowing fluid to several events as acceleration through the convergent, pipe type of flow in the throat where maximal Reynolds number (Re th ) is reached and jet formation into the sudden expansion section characterized by deceleration, mixing shear layer generation and possible turbulence transition. Five flow regimes were considered, which correspond to fully laminar case (Re th = 500), transient case (Re th = 2000), laminar inlet-transition to turbulence in sudden expansion (Re th = 3500), laminar inlet-turbulent case (Re th = 5000) and fully turbulent case (Re th = 6500).

The reported computational investigations on this FDA benchmark can be divided into two major groups: blinded studies that were conducted without foreknowledge of the experimental results and the studies conducted after the publication of first benchmark results. While none of the first 28 CFD research groups which used various RANS methods managed to fully reproduce the FDA experimental results [START_REF] Stewart | Assessment of CFD performance in simulations of an idealized medical device: results of FDA's first computational interlaboratory study[END_REF], almost all succeeding CFD studies [START_REF] Passerini | Validation of an open source framework for the simulation of blood flow in rigid and deformable vessels[END_REF][START_REF] Delorme | Large eddy simulation of FDA's idealized medical device[END_REF][START_REF] Bhushan | Laminar, turbulent, and transitional simulations in benchmark cases with cardiovascular device features[END_REF][START_REF] Janiga | Large eddy simulation of the FDA benchmark nozzle for a Reynolds number of 6500[END_REF] reported good agreement with the FDA experiments. This raises a number of questions regarding both the numerical Fig. 1 Geometry specifications of the first FDA ideal medical device franck.nicoud@umontpellier.fr methods for turbulence transitional blood flows and the overall methodology to represent the FDA blood flow case. While the RANS codes [START_REF] Stewart | Assessment of CFD performance in simulations of an idealized medical device: results of FDA's first computational interlaboratory study[END_REF][START_REF] Sotiropoulos | Computational fluid dynamics for medical device design and evaluation: are we there yet?[END_REF] in this case clearly lag behind the LES due to their inability to accurately predict the laminar to turbulent transition, the absence of analysis of the sensitivity to numerical parameters in the reported DNS and LES results [START_REF] Passerini | Validation of an open source framework for the simulation of blood flow in rigid and deformable vessels[END_REF][START_REF] Delorme | Large eddy simulation of FDA's idealized medical device[END_REF][START_REF] Janiga | Large eddy simulation of the FDA benchmark nozzle for a Reynolds number of 6500[END_REF] questions the actual predictive character of these simulations. Indeed, the FDA case proves to be highly sensitive to any physical or numerical specificities as illustrated in the remaining of this section.

Simulations with Perturbation-Free Inlet

Investigation of the effects of numerics on the results was conducted by considering three grid levels as well as a variety of CFL numbers (dimensionless time step). In order not to presume the resulting flow field in any way, unstructured finite volume grids were designed with almost no-local refinement in the domain with only slight stretching towards the upstream and downstream ends. The three grids contain respectively 5 million tetrahedral elements with h = 0.34 mm average cell height (Coarse), 15 million with h = 0.2 mm (Medium) and 50 million with h = 0.14 mm (Fine). Further details on the numerics and grids are available in [START_REF] Zmijanovic | About the numerical robustness of biomedical benchmark cases: Interlaboratory FDA's idealized medical device[END_REF].

Snapshots of instantaneous velocity fields obtained with the medium grid are displayed in Fig. 2 for the transitional (Re th = 3500) case. Except for the CFL number, the same physical and numerical conditions were used in both cases. This figure illustrates that the location of the transition to turbulence is very sensitive to the numerical details, as confirmed by the profiles displayed in Fig. 3. The latter also shows that the upstream part of the flow (x < 0, the sudden expansion starts at x = 0) is robust to the details of the simulation. However, both the streamwise mean velocity and the static pressure strongly depend on the numerics in the downstream region. This means that the axisymmetric configuration considered, when fed with an ideal Poiseuille velocity profile without perturbation, is very sensitive to any physical (boundary conditions) or numerical (mesh stretching, schemes, time step) modification. This makes hopeless any effort towards a predictive simulation of this configuration. A cure to this issue is considered in the next section.

Small Perturbations at the Inlet

When dealing with transitional flows, the perturbations injected at the inlet of the computational domain are of prime importance as already discussed in the literature [START_REF] Wu | Osborne Reynolds pipe flow: direct simulation from laminar through gradual transition to fully developed turbulence[END_REF][START_REF] Avila | The onset of turbulence in pipe flow[END_REF]. In the particular case of the FDA nozzle, the absence of such inlet perturbations leads to an extreme sensitivity of the results which makes the simulations unreliable. Note that this may also lead to simulations which may compare well with the experimental data set but which are nonetheless inadequate since not robust (e.g. the case CFL = 0.6 and coarse grid in Fig. 3a agrees well with the measurements but should not be trusted since the results depend on the time step). As shown in Fig. 4, even the smallest amount of perturbations injected at the inlet (thereafter denoted turbulence injection, TI) greatly improves the situation in the sense that the results are now robust to any tested numerical condition. More importantly, the numerical results are also systematically in good agreement with the experimental data, suggesting that the simulations outcomes are both robust and accurate. Profiles of the root-mean-square (RMS) velocity, Fig. 5a, show that the amount of flow perturbations at inlet is significantly increased by TI, as expected. The small amplitude perturbations are quickly dissipated downstream, being decreased by more than 2 orders of magnitude as shown in Fig. 5b. Figure 5c andd show that the flow in the throat is characterized by an increase of the near wall RMS. Just after the franck.nicoud@umontpellier.fr sudden-expansion, Fig. 5e andf, both RMS profiles obtained with TI are approaching the original no TI cases before merging into laminar flow profile (Re th = 500) or advancing into turbulence breakdown (Re th = 3500).

Figure 6 shows that TI greatly improves the LES prediction as well as its robustness for laminar, mildly turbulent cases and even in the most demanding, transitional case (Re th = 2000). The injected perturbations dissipate without inflicting transition in the laminar case, while in the transitional and turbulent cases they develop and contribute to more accurate and robust results.

franck.nicoud@umontpellier.fr The study thus reveals that the upstream injection of small-scale perturbations into transition sensitive flows is an indispensable tool for any biomedical case where accurate prediction of turbulence is critical.

Intracardiac Turbulence

The hemodynamics of the left heart (LH) conveys information regarding the heart function [START_REF] Kheradvar | On mitral valve dynamics and its connection to early diastolic flow[END_REF][START_REF] Pedrizzetti | On the left ventricular vortex reversal after mitral valve replacement[END_REF]. Observation of LH flow may thus reflect the presence of an existing pathology. In addition, there is substantial evidence that the LH hemodynamics can be responsible for the initiation of ventricular remodelling through mechanosensitive feedbacks modulating cardiomyocytes architecture [START_REF] Davies | Turbulent fluid shear stress induces vascular endothelial cell turnover in vitro[END_REF][START_REF] Olesen | Haemodynamic shear stress activates a K+ current in vascular endothelial cells[END_REF] and thus cardiac function [START_REF] Pasipoularides | Mechanotransduction mechanisms for intraventricular diastolic vortex forces and myocardial deformations: Part 1[END_REF]. Therefore, accurate assessment of the intracardiac flow is of paramount importance to get further comprehension of the role played by the hemodynamics in normal and abnormal LH. Used with caution, image-based computational fluid dynamics (CFD) can retrieve all the scales of the instantaneous flow, hence being able to capture disturbed flows. However, pioneering CFD studies focused only on the large-scale features of the ventricular flow i.e. the jets and the large recirculating cell. Apart from a few studies that mentioned the potential presence of flow instabilities [START_REF] Chnafa | Image-based large-eddy simulation in a realistic left heart[END_REF][START_REF] Domenichini | Three-dimensional filling flow into a model left ventricle[END_REF][START_REF] Chnafa | Image-based simulations show important flow fluctuations in a normal left ventricle: what could be the implications?[END_REF][START_REF] Le | On the three-dimensional vortical structure of early diastolic flow in a patient-specific left ventricle[END_REF], little focus has been given to the disturbed nature of the LV flow, despite its potential importance.

Method

Large-eddy simulations were preformed to compute 35 cardiac cycles using the inhouse YALES2BIO solver in a patient-specific left heart model obtained from CMR images. The set of images used to define the time evolution of the model over the cadiac cycle consists in 20 3D images of spatial resolution 5.0 × 1.1 × 1.1 mm 3 (21 × 256 × 256 voxels). The subject was 26 years old and his cardiac cycles lasted franck.nicoud@umontpellier.fr on average T = 750 ms. As shown in Fig. 7, the model includes the valves, atrium, ventricle, papillary muscles and ascending aorta. The deformations of the heart all along the heart cycle match the CMR exam thanks to an image registration technique. The technical procedure can be found in [START_REF] Chnafa | Image-based large-eddy simulation in a realistic left heart[END_REF][START_REF] Chnafa | Using image-based CFD to investigate the intracardiac turbulence[END_REF]. The flow domain within the heart is discretized using ten million tetrahedral elements. The average edge length is close to 0.55 mm during diastole. The time step, based on a Courant-Friedrichs-Lewy stability number of 0.9, varies from 0.2 ms during the beginning of diastole to 0.5 ms during diastasis.

Results

The LV large-scale flow features i.e. the blood ejection, the two vortex rings for the E and A waves and the recirculating cell, are well retrieved (Fig. 8) in accordance with the numerous observations performed in vivo [START_REF] Le | On the three-dimensional vortical structure of early diastolic flow in a patient-specific left ventricle[END_REF][START_REF] Kilner | Asymmetric redirection of flow through the heart[END_REF], in silico [START_REF] Domenichini | Three-dimensional filling flow into a model left ventricle[END_REF][START_REF] Mihalef | Patient-specific modelling of whole heart anatomy, dynamics and haemodynamics from fourdimensional cardiac CT images[END_REF] and in vitro [START_REF] Falahatpisheh | High-speed particle image velocimetry to assess cardiac fluid dynamics in vitro: From performance to validation[END_REF][START_REF] Querzoli | Effect of the prosthetic mitral valve on vortex dynamics and turbulence of the left ventricular flow[END_REF]. However, our computation reveals also large velocity fluctuations, which are usually not reported in silico but which are in line with the results of previous experimental work that used both simplified ventricle geometries and inflow boundary condition [START_REF] Domenichini | Combined experimental and numerical analysis of the flow structure into the left ventricle[END_REF][START_REF] Querzoli | Effect of the prosthetic mitral valve on vortex dynamics and turbulence of the left ventricular flow[END_REF]. In order to quantify and study the regional distribution of these fluctuations, we computed them as the difference between the instantaneous velocity components u i (i = 1, 2, 3) and the phase-averaged velocity components U i . The fluctuating part of the fluid velocity is then:

u ′ i = u i -U i = u i -⟨u i ⟩
while the corresponding fluctuating kinetic energy (FKE) per unit volume is:

franck.nicoud@umontpellier.fr

FKE = 2 ⟨ u ′ i u ′ i ⟩ ,
where ⟨⟩ denotes phase-averaged values, is the fluid density (here 1040 kg/m 3 ) and the implicit summation rule over repeated indices is used. Figure 9 displays the FKE in the left heart during systole (left figure) and diastasis (right figure). Note that non negligible levels of FKE are visible at these instants. Detailed results of this simulation and its biomedical consequences can be found in [START_REF] Chnafa | Image-based simulations show important flow fluctuations in a normal left ventricle: what could be the implications?[END_REF].

Role of SGS Model

As a common practice when using properly resolved LES, all quantities are computed from the resolved velocity field in this study [START_REF] Celik | Index of resolution quality for large eddy simulations[END_REF][START_REF] Pope | Turbulent Flows[END_REF]. We note that during our computations the subgrid scale viscosity remains low in the LV during the whole the heart cycle, showing that the subgrid scale model dissipates a moderate amount of energy. Figure 10 shows the ratio between the subgrid scale viscosity and the fluid viscosity. As expected, since the flow is laminar, the subgrid scale viscosity is virtually zero during the systole in the LV and has moderate values, but non-null, during the most turbulent part of the cycle. Thus, the -model [START_REF] Nicoud | Using singular values to build a subgridscale model for large eddy simulations[END_REF] mainly introduces dissipation in regions where there are significant turbulent fluctuations and franck.nicoud@umontpellier.fr almost no dissipation in regions where the flow is laminar. This is one of the features of the -model that makes it suitable for computing cardiovascular flows. For this particular simulation, an estimation of the Pope criterion [START_REF] Pope | Ten questions concerning the large-eddy simulation of turbulent flows[END_REF] can be computed franck.nicoud@umontpellier.fr as k sgs ∕k ≈ 3C∕2( ∕ L) 2∕3 [START_REF] Pope | Turbulent Flows[END_REF] where C = 1.5, = 0.55 mm and L = 0.04 m is the characteristic length of the largest structures. For the present LES, this leads to k sgs ∕k ≈ 0.06 which remains three times below the 15-20% threshold usually used to evaluate if a LES is sufficiently resolved [START_REF] Pope | Ten questions concerning the large-eddy simulation of turbulent flows[END_REF].

Discussion

The present results question the common idea of laminar flow in the normal left heart. Imposing the flow to be laminar (as often done in CFD for the LV using commercial software) may be sufficient to retrieve main large-scale hemodynamic features, namely jets, recirculating cell, and ejection. However, transitional or turbulent regimes feature small-scale phenomena that cannot be retrieved under the laminar hypothesis. The typical Reynolds number value (a few thousand) would make the DNS of one single cardiac cycle extremely CPU/memory demanding, if not impossible. For example, Chnafa et al. [START_REF] Chnafa | Image-based large-eddy simulation in a realistic left heart[END_REF] estimated that more than 1 billion nodes would be necessary to represent the Kolmogorov scale in their realistic human heart model. On the contrary, well-resolved large-eddy simulation allows retrieving the main characteristics of the flow and its transient properties at a lower cost. This allows simulating several cycles (here 35) at a tractable numerical cost, thus studying more complex flow features compared to what is usually done based on RANS simulations, including cycle-to-cycle variations.

Cardiac Valves

The aortic heart valve separates the left ventricle from the aorta. It is composed of three thin deformable leaflets that open and close passively during the cardiac cycle, preventing blood from flowing back into the left ventricle, and thus ensuring an unidirectional flow through the cardiovascular system. Aortic valves often degenerate and lead to either insufficiency or stenosis, which can cause the death of the patient if not treated. Usually, medication is not a sufficient treatment, and the aortic valve needs to be replaced by a prosthesis. One of the basic engineering concern for artificial aortic valve design is the hemodynamics. Indeed, the ideal aortic valve design should minimize production of turbulence which is notably known to have effect on thrombus formation (blood coagulation) [START_REF] Stein | Measured turbulence and its effect on thrombus formation[END_REF][START_REF] Yoganathan | Fluid mechanics of heart valves[END_REF]. In Sect. 4, the aortic valve was represented thanks to a simple model fed by morphological data gained from medical imaging [START_REF] Chnafa | Image-based large-eddy simulation in a realistic left heart[END_REF][START_REF] Chnafa | Image-based simulations show important flow fluctuations in a normal left ventricle: what could be the implications?[END_REF][START_REF] Chnafa | Using image-based CFD to investigate the intracardiac turbulence[END_REF]. A more realistic modelling based on a Fluid-Structure Interaction (FSI) formulation is described in this section where an experimental configuration relevant to an artificial aortic valve [START_REF] Pott | Dynamics of artificial aortic valves: a combined experimental and numerical study[END_REF] is also considered for validation.

franck.nicoud@umontpellier.fr

Numerical Method

The FSI model relies on the immersed thick boundary method (ITBM) [START_REF] Sigüenza | Validation of an immersed thick boundary method for simulating fluid-structure interactions of deformable membranes[END_REF] which was adapted from the original immersed boundary method (IBM) introduced by Peskin [START_REF] Peskin | The immersed boundary method[END_REF] to deal with 3D membranes. Two independent meshes are considered to discretize the solid valve and the fluid. The valve leaflets are represented by a moving Lagrangian mesh while the fluid is discretized by a fixed Eulerian unstructured mesh. The different steps of the ITBM are the following:

1. Knowing the displacement ⃖⃖⃖⃖ ⃗ U m of each Lagrangian node, the mechanical force ⃖⃖⃖⃖ ⃗ F m resulting from the membrane deformation is calculated, 2. The mechanical force ⃖⃖⃖⃖ ⃗ F m is regularized on the fluid mesh, giving rise to the volumetric force ⃖ ⃗ f j on each fluid node, 3. The Navier-Stokes equations (forced by the regularized mechanical forces) are solved on the fluid mesh, yielding the velocity of the fluid ⃖⃖ ⃗ v j on each fluid node, 4. The velocity of the membrane ⃖⃖⃖⃖ ⃗ V m at each Lagrangian node is interpolated from the ⃖⃖ ⃗ v j field, enabling to deduce the new position ⃖⃖⃖⃖ ⃗ X m from the position at the pre-

vious timestep ⃖⃖⃖⃖ ⃗ X m = ⃖⃖⃖⃖ ⃗ X m previous + t ⃖⃖⃖⃖ ⃗ V m
The displacement of each Lagrangian node is then updated thanks to

⃖⃖⃖⃖ ⃗ U m = ⃖⃖⃖⃖ ⃗ X m - ⃖⃖⃖⃖ ⃗ X 0
m , where ⃖⃖⃖⃖ ⃗ X 0 m stands for the initial stress-free position (also referred to as the reference position).

The regularization and interpolation of steps 2 and 4 require dealing with Dirac functions that must be properly regularized on the Eulerian fluid mesh [START_REF] Mendez | An unstructured solver for simulations of deformable particles in flows at arbitrary Reynolds numbers[END_REF]. To this respect, the Reproducing Kernel Particle Method [START_REF] Pinelli | Immersed-boundary methods for general finitedifference and finite-volume Navier-Stokes solvers[END_REF] is used to deal with unstructured meshes. The computation of the mechanical force ⃖⃖⃖⃖ ⃗ F m (step 1) is performed by the LMGC90 solid mechanics solver [START_REF] Radjai | Discrete Numerical Modeling of Granular Materials[END_REF], using the classical finite-element method. In the present work, a quasi-incompressible Neo-Hookean material is used to model the valve, defined by its strain energy function:

W = G 2 Ī1 -3 + K 2 (ln J) 2 , ( 3
)
where G and K are the shear and bulk modulus, respectively. Ī1 is the first invariant of the isochoric right Cauchy-Green deformation tensor ̄ , and J is the Jacobian of the transformation.

franck.nicoud@umontpellier.fr

Computational Setup

The computational setup was designed in such a way to represent the experimental test rig studied at the Helmholtz Institute of Aachen [START_REF] Pott | Dynamics of artificial aortic valves: a combined experimental and numerical study[END_REF]. The aortic geometry includes the three sinuses of Valsalva after the calculations by Reul et al. [START_REF] Reul | The geometry of the aortic root in health, at valve disease and after valve replacement[END_REF] (see Fig. 11). The valve is built on a tri-commissure frame positioned within the aortic geometry, and is designed in an almost closed position (see Fig. 11). Both the arterial wall and the tri-commissure frame are defined as rigid bodies whose borders are non-slip wall boundary conditions of the fluid domain. As shown in Fig. 11, a flow rate signal is imposed as a periodic inlet condition either upstream or downstream of the valve, depending if the signal is positive or negative. The whole fluid geometry is represented by approx. 3.3 million tetrahedral elements, and the zone of interest (displayed in Fig. 11) is meshed with a uniform mesh resolution of 0.5 mm. The valve mesh is composed of approx. 1500 quadratic hexahedral elements, with an effective resolution of 0.5 mm.

The characteristics of the flow and the valve are relevant of physiological data and are summarized in Table 1. In view of the Reynolds (R e ) and Womersley (W 0 ) numbers, the viscosity of the Newtonian blood analogue fluid is set in the high-shear limit values [START_REF] Robertson | Rheological models for blood[END_REF].

Results and Discussion

Figure 12 shows the streamwise velocity field downstream of the valve. Four characteristic instants of the cardiac cycle are shown: t = 0.20 s, just before the flow rate reaches its maximal value, referred to as Early Systole (ES); t = 0.26 s, when the flow rate is maximal, referred to as Peak Systole (PS); t = 0.35 s, just after the flow franck.nicoud@umontpellier.fr The main feature of the flow is the propagation of a jet emerging from the valve at ES (see label a), propagating between ES and MS (see labels b and c), and dissipating at LS. This jet is most probably hydro-dynamically unstable, promoting transition to turbulence downstream of the valve. In order to analyse the structure of this jet and highlight the cycle-to-cycle variations of the flow, two velocity profiles are extracted over 10 distinct cycles (displayed in Fig. 12): profile 1 (P1), located immediately behind the valve; profile 2 (P2), located further downstream of the valve. Cycle-to-cycle variations are mainly observed on P2, suggesting that the jet emerging from the valve is mostly laminar and that it becomes turbulent when propagating further downstream. Regarding P2 at PS, MS and LS, it is seen that cycle-to-cycle variations increase as the inlet flow rate decreases, which is consistent with the observations made by Chnafa et al. [START_REF] Chnafa | Image-based large-eddy simulation in a realistic left heart[END_REF][START_REF] Chnafa | Image-based simulations show important flow fluctuations in a normal left ventricle: what could be the implications?[END_REF] in the case of whole heart geometries with simple (without FSI) valve model. This feature should be related to the general trend of all decelerated flows to destabilize. Few cycle-to-cycle variations are nonetheless observed on P1 at ES; these fluctuations are in fact induced by variations in the valve opening process and not related to any turbulent motion. This is illustrated in Fig. 13 which displays some typical shapes of the aortic valve. All these views correspond to the same three instants (the third one corresponding to ES) but were extracted at four different cycles. If the overall opening process is always the same (buckling of the leaflets, either successive or simultaneous), very large fluctuations are also present from one cycle to the other. These variations most probably result from two main ingredients: the sensitivity of the buckling phenomenon to the details of the load on the first hand, the existence of a chaotic pressure load induced by the downstream turbulent fluctuations on the other hand.

Conclusion

Given the flow regimes encountered in macro-circulation, turbulence is most probably present in many practical configurations, either physiological or within biomedical devices. Large-Eddy Simulation has the potential to properly represent the effects of turbulence provided that it is used with appropriate numerics, subgrid scale model and boundary conditions. The numerical solver YALES2BIO (http://www. math.univ-montp2.fr/~yales2bio/) gathers all the required properties and was successfully used in three different situations including the intracardiac flow and the interaction with a realistic aortic valve. Transitional flows are even more demanding and should be addressed very carefully since they may be extremely sensitive to the details of the simulation if no perturbation is introduced in the upstream condition. This is illustrated by considering the simplified medical device proposed by FDA as a validation test case.

Fig. 2 6 Fig. 3

 263 Fig. 2 Instantaneous velocity fields (medium grid) for two different values of the CFL numbers (top CFL = 0.6; bottom CFL = 0.1)-Re th = 3500

Fig. 4 Fig. 5

 45 Fig.[START_REF] Mikhal | Immersed boundary method for pulsatile transitional flow in realistic cerebral aneurysms[END_REF] Mean axial velocity and pressure profiles with small TI at inlet (Medium grid, Re th = 3500)

Fig. 6

 6 Fig. 6 Axial velocity profiles for the laminar, transitional and mildly turbulent cases with and without TI (Medium grid)

Fig. 7

 7 Fig. 7 Human left heart extracted from CMR images. The same domain is shown from four different points of view. The inlets and outlet flow extensions are visible in the left figure. The left ventricle (LV), left atrium (LA), Aorta (AO), Antero Lateral (A-L) and Postero Medial (P-M) papillary muscles are indicated. A black line passing through the LH indicates the position of the slice which will be used to display the results in Figs. 8 and 10

Fig. 8

 8 Fig. 8 Phase-averaged velocity field over a cutting plane through the LH (see Fig. 7 for the plane position). The mitral and the aortic valves are depicted in grey. The figure shows the classically reported structure of the heart cycle: the strong blood ejection at systole (0.15 T), the E waveinduced vortex ring at diastole (0.65 T) and the recirculating cell during diastasis (0.75 T)

Fig. 9 Fig. 10

 910 Fig. 9 3D volume rendering of the FKE in the left heart during systole (left, 0.15 T) and diastasis (right, 0.75 T). The opacity is set to 0% for FKE values below 30 J/m 3 . The mitral and aortic valves are not displayed for more clarity but were accounted for in the simulation

Fig. 11

 11 Fig. 11 Details of the computational setup. The flow is generated using switching inlet and outlet conditions, imposed far from the zone of interest. The assumption of a flat velocity profile is made to impose the inlet/outlet flow rate signal

Fig. 12 Fig. 13

 1213 Fig. 12 Streamwise velocity field downstream of the valve. Two velocity profiles (P1 and P2) are extracted over 10 distinct cycles. Four different instants of the cardiac cycle are depicted (ES, PS, MS and LS)

12

 12 Fig.[START_REF] Vichnevetsky | Fourier Analysis of Numerical Approximations of Hyperbolic Equations[END_REF] Typical shapes of the aortic valve during the opening phase. The last instant (t = 0.20 s) corresponds to ES (see Fig.12). Each column corresponds to a specific cycle

  Fig.[START_REF] Vichnevetsky | Fourier Analysis of Numerical Approximations of Hyperbolic Equations[END_REF] Typical shapes of the aortic valve during the opening phase. The last instant (t = 0.20 s) corresponds to ES (see Fig.12). Each column corresponds to a specific cycle

Table 1

 1

	Characteristics of the flow and the valve set to reproduce physiological flow	Flow parameters Density	= 1100 kg/m 3
	conditions	Dynamic viscosity	= 3.6 × 10 -3 Pa s
		Heart rate	n bpm = 60
		Mean cardiac output	Q mean = 3.48 L/min
		Reynolds number	R e = 1388
		Womersley number	W 0 = 17
		Valve parameters
		Density	s = 1000 kg/m 3
		Shear modulus	G = 2.4 MPa
		Bulk modulus	K = 1.6 MPa
		Thickness of the valve	e = 0.15 mm
		Radius of the valve	R = 12.5 mm

rate reaches its maximal value and begins to decrease, referred to as Mid-Systole (MS); t = 0.55 s, when the flow rate is decreasing, referred to as Late Systole (LS).

franck.nicoud@umontpellier.fr

Acknowledgements

The authors thank ANR and BPI for their supports through the Laboratory of Excellence NUMEV (ANR-10-LABX-20), the FORCE (ANR-11-JS09-0011) and the DAT@DIAG (ISI-I1112018W) project. Dr S. Sonntag and PhD student D. Pott from the Helmholtz Institute of Aachen are gratefully acknowledged for providing the details of their experimental test rig about the aortic valve dynamics. CC also thanks CNRS for funding his thesis. Dr. V. Moureau is gratefully acknowledged for giving access to the YALES2 solver. This work was performed using HPC resources from GENCI-CINES (Grants 2014-, 2015-and 2016-c2015037194) and the HPC@LR Center.

franck.nicoud@umontpellier.fr


F Della Croce

T Garaix

A Grosso

ITERATED LOCAL SEARCH AND VERY LARGE NEIGHBORHOODS FOR THE PARALLEL-MACHINES TOTAL TARDINESS PROBLEM

We present computational results with a heuristic algorithm for the parallel machines total weighted tardiness problem. The algorithm combines generalized pairwise interchange neighborhoods, dynasearch optimization and a new machine-based neighborhood whose size is non-polynomial in the number of machines. The computational results significantly improve over the current state of the art for this problem.

Introduction

In the P m| | j w j T j problem a set of jobs N = {1, 2, . . . , n} are given, with their processing times p j , weights w j and due dates d j . The jobs are to be processed in a schedule S on a set M = {1, 2, . . . , m} of identical parallel machines so that their completion times C j minimize the objective function

T (S) = n j=1 w j T j = n j=1
w j max {C j -d j , 0} .

For m = 1 the problem reduces to the single-machine total tardiness problem, that is well studied and solved in both the exact and heuristic frameworks -we refer to [START_REF] Bigras | Time-Indexed Formulations and the Total Weighted Tardiness Problem[END_REF][START_REF] Pan | On the equivalence of the max-min transportation lower bound and the timeindexed lower bound for single machine scheduling problems[END_REF][START_REF] Rodrigues | Heuristic algorithm for the parallel machine total weighted tardiness scheduling problem[END_REF][START_REF] Rodrigues | Algorithms over Arc-time Indexed Formulations for Single and Parallel Machine Scheduling Problems[END_REF] and [START_REF] Congram | An iterated dynasearch algorithm for the single machine total weighted tardiness problem[END_REF][START_REF] Ergun | Fast neighborhood search for the single machine total weighted tardiness problem[END_REF][START_REF] Grosso | An enhanced dynasearch neighborhood for the single machine total tardiness problems[END_REF] for recent developments.

The literature seems to be fairly limited for the problem with parallel machines; the most recent references are [START_REF] Anghinolfi | Parallel machine total tardiness scheduling with a new hybrid metaheuristic approach[END_REF][START_REF] Bilge | A tabu search algorithm for parallel machine total tardiness problem[END_REF][START_REF] Koulamas | Decomposition and hybrid simulated annealing heuristics for the parallelmachine total tardiness problem[END_REF][START_REF] Rodrigues | Heuristic algorithm for the parallel machine total weighted tardiness scheduling problem[END_REF][START_REF] Rodrigues | Algorithms over Arc-time Indexed Formulations for Single and Parallel Machine Scheduling Problems[END_REF] to the authors' knowledge.

Iterated Local Search (ILS, see [START_REF] Lourenco | Iterated local search, in Handbook of MetaheuRistics[END_REF] for a survey) is a local search framework that can be seen as a tradeoff between the naive multistart and complex metaheuristics. In multistart a local search driven optimization starts several times (often a huge number of times) from randomly generated initial solutions, in order to achieve a wide exploration of the solutions set. In metaheuristics a number of sophisticated devices (genetic crossover, short or long-term memory, etc) are employed in order to escape poor local optima. In ILS the search is simply restarted from a slightly perturbed version of the best-known solution. With this type of restart, the starting point of each local search is not completely random, and the perturbation -called "kick" -aims at projecting the search "not too far" from previously explored local optima, without completely loosing their partially optimized structure.

Very Large Neighborhood Search (VLNS) denotes local search methods that define and explore complex neighborhoods for combinatorial optimization problems; such neighborhoods are characterized by having an exponential number of neighbor solutions -with respect to the problem size -but can be explored in polynomial time by means of exact or heuristic procedures (see [START_REF] Ahuja | A survey of very large-scale neighborhood search techniques[END_REF]).

ILS is often successfully coupled with VLNS, hence moving the complexity of the search from the overall algorithm to the neighborhood exploration. We refer to [START_REF] Congram | An iterated dynasearch algorithm for the single machine total weighted tardiness problem[END_REF][START_REF] Ergun | Fast neighborhood search for the single machine total weighted tardiness problem[END_REF] for a successful application of such a VLNS technique (called dynasearch) to the 1| | w j T j problem.

Rodrigues et al. [START_REF] Rodrigues | Heuristic algorithm for the parallel machine total weighted tardiness scheduling problem[END_REF] proposed a simple and quite effective ILS algorithm for the P m| | j w j T j problem, using a local search based on pairwise interchange operators. That algorithm was tested on a batch of 100 instances with n = 40, 50 and m = 2, 4 derived from a subset of the 1|| w j T j problem instances available in the OR-library1 . Notice that, on that batch, the algorithm was able to detect all but one optimal solutions. This paper aims at defining an improved ILS algorithm for P m| | j w j T j by incorporating VLNS techniques. Particularly, we introduce a dynasearch optimization on each machine in the shop and a new "Very Large" neighborhood whose size is non-polynomial in the number of machines.

We illustrate the basic building blocks of the algorithm and present computational experiments for assessing their effectiveness in Section 2. The complete algorithm is described in Section 3, where also the computational results are discussed. The proposed ILS algorithm outperforms the ILS of [START_REF] Rodrigues | Heuristic algorithm for the parallel machine total weighted tardiness scheduling problem[END_REF] on instances with a number of jobs n ranging from 40 to 300, and a number of machines m ranging from 2 to 20. The advantage of the new algorithm grows on instances with large m thanks to the new neighborhood.

The basic neighborhoods

2.1. Generalized pairwise interchanges. The well-known GPI operators work on a sequence of jobs σ producing a new sequence σ Let σ = αiπjω, with jobs i and j in position k and l respectively. The most common GPI operators are (1) Swap αiπjω → αjπiω (π may be empty);

(2) Forward insertion αiπjω → απjiω;

(3) Backward insertion αiπjω → αjiπω.

(4) Twist αiπjω → αj πiω with π = π reversed. The implementation of such operators is straightforward in single-machine sequencing problems with regular cost functions, since the machine is never idle and the sequence σ is the schedule. The so-called GPI dynasearch neighborhood for singlemachine sequencing problems combines possibly many independent moves of types (1)-(4); two moves are said to be independent if the pairs of positions (k, l) and (p, q) on which they act are non-overlapping, i.e max {k, l} < min {p, q}. In a single machine environment with an additive objective function the contributions of independent moves combine additively, and the best set of independent moves can be worked out by dynamic programming (see [START_REF] Congram | An iterated dynasearch algorithm for the single machine total weighted tardiness problem[END_REF] for details). A GPI dynasearch neighborhood exploration for an n jobs sequence requires O(n 2 ) time with its best implementation (see [START_REF] Ergun | Fast neighborhood search for the single machine total weighted tardiness problem[END_REF]).

In parallel-machines environments, GPI operators can be applied provided that a sequence σ can be converted to a schedule. Rodrigues et al. [START_REF] Rodrigues | Heuristic algorithm for the parallel machine total weighted tardiness scheduling problem[END_REF] proposed a simple yet quite effective ILS algorithm for the P m| | j w j T j problem. The algorithm applies GPI operators -limited to (1)-( 3) in their implementation -on a sequence of jobs; the schedule on parallel machines associated with this sequence σ = (j 1 , j 2 , . . . , j n ) is computed from scratch by means of the most natural dispatching rule: assign the next job in the sequence to the earliest available machine. The neighborhood exploration is performed with a first-improve strategy, and frequent restarts are applied (one kick every five complete descents).

Whereas the basic GPI neighborhood can be easily adapted to the parallel machines environment, this is not the case for the GPI dynasearch neighborhood: since the job starting times are determined by applying the dispatching rule, the contribution of independent moves is no longer purely additive. Rodrigues et al. [START_REF] Rodrigues | Heuristic algorithm for the parallel machine total weighted tardiness scheduling problem[END_REF] do not provide a different notion of independent moves, neither it is easy to see an obvious one.

2.2.

Integrating GPIs on parallel machines and Dynasearch. A possible drawback of the basic GPI neighborhood is that, in a parallel environment, the working sequence on each single machine is poorly optimized, since the machinesequencing criterion is extremely crude. We then investigated the opportunity of adding a single-machine optimization phase through the use of a dynasearch neighborhood. We tested the following algorithms, called 1 and 2 respectively.

Algorithm A1: The GPI iterated local search of [START_REF] Rodrigues | Heuristic algorithm for the parallel machine total weighted tardiness scheduling problem[END_REF] (kindly provided by the authors). Algorithm A2: The same algorithm, where after building a schedule by the dispatching rule, each single machine is optimized by a full descent using the GPI dynasearch neighborhood where moves ( 1)-( 4) are used. We considered a batch of 125 randomly instances with n = 100, m = 4; we refer to Section 3.3 for details on the generation scheme. Two parameters R, T determine the practical difficulty of the instances. We recorded the performances of the algorithms in terms of time spent for reaching the best solution and number of local search descents performed. For both algorithms we allowed one hour of CPU time. Table 1 points to a comparison of the computational costs of the two algorithms in terms of CPU time and number of descents, detailing them by (R, T ) pairs. Out of the 125 instances in the batch, Algorithm A2 delivered better solutions in 37 cases, and worse solutions in 27 (columns labeled "#Bests"). The higher number of better solutions comes at the cost of higher CPU times to be spent in the search. The number of descents required to reach the best solution is always consistently less for Algorithm A2 than for Algorithm A1, but Algorithm A2 -quite expectedly -exhibits in most cases higher CPU times, since every solution undergoes a full dynasearch descent on each machine. Anyway, in the details of the tests we were able to observe that on 18 instances of the batch, Algorithm A2 finds a better solution and requires less CPU time to reach the optimum; this behaviour comes out with dramatic evidence on some classes of instances like R = 0.6, T = 0.6, and the classes with R = 1.0. This test suggests that an effort for keeping highly optimized sequences on the machines can be worth, if a clever search strategy can be developed in order to limit the growth of the CPU time 2.3. A parallel-machines neighborhood. A simple notion of independent moves arises if, instead of using the dispatch rule on a sequence, one works directly on the jobs-to-machines assignment. Assume a schedule is given. We consider, for a given pair of machines m 1 , m 2 , the respective job working sequences σ 1 , σ 2 and the following moves:

(a) extract a job j from σ 1 and insert it into σ 2 ;

(b) extract a job j from σ 2 and insert it into σ 1 ; (c) extract jobs i ∈ σ 1 , j ∈ σ 2 and insert i in σ 2 -{j}, and j in σ 1 -{i}. The insertion position for a job j in a sequence σ is chosen in linear time, so that σ = αjω and T (αjω) is as small as possible. We note that if such moves are executed on disjoint pairs of machines (m 1 , m 2 ), (m 3 , m 4 ), . . . , they are independent moves, in the sense that their contributions to the schedule's cost combine additively. Hence a neighborhood whose size is non-polynomial in the number of machines can be defined as follows, for a given schedule.

(1) For each pair of machines (m 1 , m 2 ), compute the maximum decrease in tardiness ∆ m1,m2 that can be obtained by applying moves (a), (b) and (c) to σ 1 , σ 2 , for all j ∈ {σ 1 } ∪ {σ 2 }.

(2) Build a weighted improvement graph G(M, E) where

E = {m 1 , m 2 } : ∆ m1,m2 > 0 .
(3) A neighbor schedule is generated by taking a matching on G and executing the moves associated to the matching edges. Note that, for each pair of machines:

• no more than O(n) jobs have to be considered for moves (a) and (b);

• no more than O(n 2 ) pairs i, j have to be considered for move (c);

• evaluating each move of type (a), (b), (c) requires O(n) operations for each job (or job pair) -note that the schedule defines the machine to which the job is assigned. Hence evaluating all the possible moves (a), (b), (c) requires O(n 3 ) operations. The best neighbor schedule can be computed by taking a maximum weighted matching in G; the whole process for building G and selecting the matching can be implemented with running time

O(n 3 ) + O(m 3 ) ≈ O(n 3 )
(as m < n in non trivial instances).

Simple testing showed that the solutions provided by the GPI ILS technique of [START_REF] Rodrigues | Heuristic algorithm for the parallel machine total weighted tardiness scheduling problem[END_REF] are often not locally optimal with respect to the parallel machines neighborhood. Particularly, applying one search of the parallel machine neighborhood on the best solution provided by Algorithm A1 after a 1-hour run, we found 11 improvements (over 125 solutions) for the instances in the n = 100, m = 10 batch; the number of improvement rises to 96 for the n = 300, m = 20 batch.

The exploration of the parallel neighborhood is usually fast (less than 0.6 seconds on the n = 300 instances). Hence, we keep the parallel machines neighborhood as a cheap and useful tool whose impact becomes more and more important on instances with many machines.

Combining different neighborhoods

3.1. Neighborhoods and refinements. In view of the experimental observations reported in the previous section, a careful combination of GPI moves, dynasearch optimization and the parallel machines neighborhood can be the key instrument for handling larger instances of the P m| | j w j T j problem. We used three neighborhoods called N 1 , N 2 , N 3 .

• Neighborhood N 1 is the GPI neighborhood of Rodrigues et al. [START_REF] Rodrigues | Heuristic algorithm for the parallel machine total weighted tardiness scheduling problem[END_REF]. Note that in [START_REF] Rodrigues | Heuristic algorithm for the parallel machine total weighted tardiness scheduling problem[END_REF] a somewhat sophisticated rule is used for breaking ties in choosing the best neighbor. We avoided it in favor of a random tie-breaking rule. • Neighborhood N 2 is the parallel-machines neighborhood described in Section 2. • Neighborhood N 3 is a GPI neighborhood where every neighbor schedule generated by a GPI operator is improved by a dynasearch descent applied on each machine schedule. In order to reduce the computational effort, we applied the following refinements to N 1 and N 3 . Following [START_REF] Rodrigues | Heuristic algorithm for the parallel machine total weighted tardiness scheduling problem[END_REF], the neighborhood search proceed by first-improve, hence the first improving neighbor is adopted as new solution. We observed that often profitable moves happen between jobs that appear in relatively "close" positions. Hence the GPI operators are applied in "stages"; each stage applies the GPI operators between jobs in positions i and i + γ mod n, for i = 1, . . . , n, γ fixed.

In exploring N 1 , N 3 , at each successive stage γ is set to 1, 2, . . . , n-1. The improving neighbor is often found at early stages.

Accordingly with [START_REF] Rodrigues | Heuristic algorithm for the parallel machine total weighted tardiness scheduling problem[END_REF], a full exploration of a N 1 neighborhood is accomplished in O(n 3 log m) operations. An exploration of N 2 requires O(n 3 ) operations (see Section 2). Exploring N 3 takes up to O(n 4 ) operations.

3.2. The algorithm. We now describe the complete proposed algorithm. It uses all three neighborhoods N 1 , N 2 , N 3 . Among these neighborhoods, the local search phase for N 1 , N 2 is exploited to a full descent, while for N 3 it is limited to a single neighborhood exploration because of the higher computational cost of such procedure. This is denoted in the pseudocode by the operations called FullDescent and Search.

The search starts from the usual EDD (dispatched) sequence, which is often accepted in literature as a reasonable quick-and-dirty starting point. A full descent of N 1 performs the first optimization in the main loop, then N 2 and one exploration of N 3 are invoked iteratively one after the other, as long as the iteration is profitable. The Kick phase -i.e. perturbation of the best-known solution -consists of a limited number of random swaps, removal and insertions performed among different machines.

The algorithm uses a time-limit as stopping criterion, since we considered it the most simple and tunable one.

Algorithm A3

1: Set S * := S := The EDD dispatched sequence ; 2: Set IterCount := 0; 3: repeat Improve S by executing dynasearch on each machine; until N 3 failed improving S 2 ; end if 29: until Time-limit exceeded Following [START_REF] Rodrigues | Heuristic algorithm for the parallel machine total weighted tardiness scheduling problem[END_REF], a kick is executed each MaxNoImprove non-improving iterations, with MaxNoImprove = 5. Also, as a further refinement, the number of stages in the exploration on N 3 was fixed to a maximum of γ max = 5; note that a higher value of γ max causes more time to be spent in exploring N 3 and, correspondingly, a lower number of kicks executed in the allowed time limit. The values MaxNoImprove = 5, γ max = 5 gave the best results in some preliminary testsonly a modest amount of testing was needed to identify this value, without need for extensive calibration. The value γ max = 5 actually lowers the time spent for exploring N 3 to O(γ max n 3 ).

3.3.

Evaluation of the algorithm. We tested the hybrid ILS algorithm on batches of random instances adapted from the well established literature on tardiness problems in single-machine environments. The single-machine instances are characterized by uniformly distributed random data with processing times p i and weights w i from [START_REF] Ahuja | A survey of very large-scale neighborhood search techniques[END_REF]100], and due dates from a uniform distribution whose bounds are determined by two parameters R, T called due date range and tardiness factor -see for example References [START_REF] Crawuels | Local search heuristics for the single machine total weighted tardiness scheduling problem[END_REF][START_REF] Potts | A branch and bound algorithm for the total weighted tardiness problem[END_REF]. Specifically, the d i values are randomly drawn from

[(1 -T -R/2) n i=1 p i , (1 -T + R/2) n i=1 p i ].
For fixed n, m, we considered five instances for each (R, T ) pair, with R, T ∈ {0.2, 0.4, 0.6, 0.8, 1.0} -25 (R, T ) pairs, and 125 instances. For n = 40, 50, 100 we used times, weights, and due dates from the 375 single-machine OR-library instances. For n = 300 we used the instances from Tanaka et al. [START_REF] Tanaka | An exact algorithm for single-machine scheduling without machine idle time[END_REF]. The number of machines m was set to 2, 4 and 10, and pushed up to 20 for the largest instances. The due dates are adapted to the parallel machines case by scaling the due dates by 1 m (rounding down the obtained values). Tables 2-4 focus on the comparison of Algorithm A1 by Rodrigues et al. [START_REF] Rodrigues | Heuristic algorithm for the parallel machine total weighted tardiness scheduling problem[END_REF] and Algorithm A3. Accordingly with [START_REF] Rodrigues | Heuristic algorithm for the parallel machine total weighted tardiness scheduling problem[END_REF] we allowed at most one hour CPU time to each test and report an aggregated comparison in Table 2 2 ; the results obtained within shorter CPU times are also presented in Tables 3 and4. The total number of instances tested was 1625. For both A1 and A3 the results of a single run are reported; athough the kick phase accounts for some nondeterminism in A3, we did not observe it delivering significantly different tadiness values in different runs, as far as the time limits reported in the tables are allowed.

In Table 2 we report information on the behaviour of the algorithms running with a 3600 seconds time limit. The "dev" columns report the percentage average and maximum deviations of the objective function delivered by A1 and A3 with respect to the best found value -the "best" value is defined as the minimum between the tardiness value obtained by the two algorithms within the allowed 1hour run. The column # best counts the number of instances (out of the 125) where Algorithm A1 or Algorithm A3 delivered a better solution. CPU best reports the average time-to-best for both algorithms and N desc the average number of descent performed. The performances of the two algorithms are basically comparable for "small" instances (say for n = 40, 50); with such limited problem sizes, both A1 and A3 often find an optimal solution.

Tables 3 and4 compare the behaviour of the two algorithms for different values of the time limits enforced. The n = 300 cases are not reported for time limits ≤ 30 secs (Table 3) because in several instances such time limits where not enough to perform a full execution of A3 -this is due to the heavier computational requirements of the dynasearch component of A3. Aside from this limitation, A3 is seen to strongly outperform A1 in terms of solution quality.

A3 becomes apparently the best option for large n (n = 100, 300), and especially for instances with a large number of machines (m = 10, 20). On the latter 2 per-instance results are available as a compressed file at www.di.unito.it/~grosso/solution.tar.7z instances the key factors for the success of A3 are the ability to exploit the parallelmachines "very large" neighborhood, and the powerful dynasearch neighborhood for optimizing each machine sequence. Table 3. Comparison for small time limits (≤ 30 sec).

5 :S 1 :S 2 : 8 :S 3 :

 51283 = FullDescent(S, N 1 ); = FullDescent(S 1 , N 2 ); = Search(S 2 , N 3 );

9 :

 9 if T (S 3 ) < T (S 2 ) then 10:Set S 1 := S 3 ;

13 :if N 2

 132 if T (S 3 ) ≥ T (S * ) then 14: Set IterCount := IterCount + 1; and N 3 failed to improve S 1 then 20: if IterCount > MaxNoImprove then 21: S := Kick(S * );

Table 1 .

 1 Basic

					Algorithm 2				Algorithm 1
		R	T	CPU avg N desc avg #Bests		CPU avg N desc avg #Bests
		0.2 0.2	21.17	3.60	0		1.36	4.40	0
		0.2 0.4	27.65	3.20	0		7.10	27.00	0
		0.2 0.6 1164.99	147.80	1		301.41	1266.40	0
		0.2 0.8 2246.95	252.00	1		540.10	3793.40	0
		0.2 1.0 1690.30	205.20	2		723.50	6689.40	0
		0.4 0.2	121.81	21.00	0		5.43	22.60	0
		0.4 0.4	46.43	4.80	0		25.07	99.20	0
		0.4 0.6	752.06	90.80	2		1167.06	5924.60	1
		0.4 0.8 1662.31	178.60	1		1057.19	8114.60	4
		0.4 1.0 1370.86	165.20	2		976.78	8721.60	3
		0.6 0.2	4.00	0.00	0		0.07	0.00	0
		0.6 0.4	143.21	19.00	0		163.62	614.40	0
		0.6 0.6	414.91	45.00	2		1472.98	8029.20	0
		0.6 0.8 2075.77	223.00	3		2085.35 14456.20	2
		0.6 1.0	901.48	91.40	3		2178.29 18166.40	2
		0.8 0.2	4.20	0.00	0		0.07	0.00	0
		0.8 0.4	546.14	84.40	0		210.41	890.80	1
		0.8 0.6 1808.82	196.40	2		1781.18 10169.60	2
		0.8 0.8 2378.09	244.40	4		1562.70 10504.00	1
		0.8 1.0 3004.82	339.80	3		1149.06	9383.20	1
		1.0 0.2	4.05	0.00	0		0.07	0.00	0
		1.0 0.4	591.58	84.00	0		420.00	2359.80	0
		1.0 0.6 1459.21	147.00	3		1855.63 10910.00	2
		1.0 0.8 1015.41	107.20	5		2073.00 13795.40	0
		1.0 1.0 1879.91	210.40	3		2338.92 17907.40	2
						GPI local search (Algorithm 1) and
		GPI+dynasearch (Algorithm 2). Comparison for n = 100,
		m = 4.
	n m			dev %			# best	CPU best (sec)	N desc
		A3 avg A1 avg A3 max A1 max	A3 A1	A3	A1	A3	A1
	40 2	0.0	0.0	0.0	0.0	1	0	0.0	0.2	464	4819
	40 4	0.0	0.0	0.0	5.5	2	0	0.9	0.1	2174	8072
	40 10	0.0	0.0	0.0	1.5	8	0	0.3	0.2	34780 64940
	50 2	0.0	0.0	0.0	0.0	0	0	0.1	0.1	486	4414
	50 4	0.0	0.0	0.0	0.6	4	1	2.5	1.5	3077	9885
	50 10	0.0	0.1	0.0	8.1	13	0	0.2	29.6	75898 111079
	100 2	0.0	0.0	0.0	0.0	14	0	52.3	35.3	5480 44630
	100 4	0.0	0.0	0.0	0.1	55	4	835.7 1030.2	8958 23856
	100 10	0.0	0.3	0.0	9.0	86	1	1525.6 1139.7	11113 11562
	300 2	0.0	0.0	0.1	0.3	72	5	2800.9 3308.3	183	1192
	300 4	0.0	0.2	0.0	11.3	89	3	3290.4 1148.8	279	628
	300 10	0.0	0.8	0.0	21.3	102	0	2819.9 3022.6	346	287
	300 20	0.0	0.8	0.0	28.3	104	1	1663.1 2149.4	336	207

Table 2 .

 2 Overall results (1 hour time limit).

	n m	dev % 5 s.		dev % 10 s.	dev % 30 s.
		avg	max	avg	max	avg	max
		A3 A1 A3	A1 A3 A1 A3	A1 A3 A1 A3

Table 4 .

 4 3 16.4 76.2 0.3 1.3 16.4 76.2 0.3 1.2 16.4 76.2 0.3 0.7 13.6 18.1 0.3 0.6 13.6 17.9 300 10 0.7 3.2 17.3 72.6 0.6 2.1 17.2 41.9 0.6 1.6 17.2 37.5 0.5 1.5 15.8 37.5 0.4 1.5 15.8 37.5 300 20 0.7 4.3 11.4 96.1 0.6 2.2 11.4 50.4 0.5 2.0 11.4 50.4 0.5 1.9 11.4 50.4 0.4 1.8 10.3 50.4 Comparison for large time limits

http://people.brunel.ac.uk/~mastjjb/jeb/info.html