Karl Töpperwien
email: karl.topperwien@centralesupelec.fr

Félix Collin-Bastiani

Eleonore Riber

Bénédicte Cuenot

Guillaume Vignat

Kevin Prieur

Daniel Durox

Sébastien Candel

Ronan Vicquelin

Karl T Öpperwien

Élix Collin-Bastiani

B Én Édicte Cuenot

S Ébastien Candel

Large-Eddy Simulation of Flame Dynamics During the Ignition of a Swirling Injector Unit and Comparison With Experiments

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

A safe and reliable ignition of combustion chambers under a wide range of operating conditions is of vital importance for aeronautical applications and, more generally, for gas turbines. This can be a challenging task within the constraints set by emission legislation, lean overall equivalence ratios, strong local inhomogeneities and complex ge- * Address all correspondence to this author. ometries. Aero-engine combustors are operated with liquid fuels which may cause additional difficulties for high altitude reignition when combustor temperatures and pressure are low and flow velocities are reduced, thus leading to poor fuel atomization and mixing. Research is therefore needed not only to better understand the phenomena which lead to successful ignition but also to avoid ignition failure and hardware damage that may be caused by inadequate management of the ignition transient.

Ignition is usually initiated by external forcing, e.g., a spark, although concepts for auto-igniting combustors exist in sequential or multi-staged stationary gas turbines [START_REF] Schulz | Autoignition flame dynamics in sequential combustors[END_REF][START_REF] Ebi | Ignition of a sequential combustor: Evidence of flame propagation in the autoignitable mixture[END_REF]. The classification into the following three phases of forced ignition has been well established in the literature [START_REF] Lefebvre | Gas turbine combustion: alternative fuels and emissions[END_REF][START_REF] Mastorakos | Ignition of turbulent nonpremixed flames[END_REF]: [START_REF] Schulz | Autoignition flame dynamics in sequential combustors[END_REF] Kernel generation, (2) Flame growth and stabilization and (3) Burner-to-burner flame propagation, called light-round ignition. Extensive research has therefore been dedicated to the characterization of a successful ignition process, including detailed investigations of spark characteristics such as the minimum spark energy [START_REF] Ballal | The influence of flow parameters on minimum ignition energy and quenching distance[END_REF], ignition probability and impact of turbulence [START_REF] Ahmed | Spark ignition of turbulent nonpremixed bluff-body flames[END_REF][START_REF] Esclapez | Ignition probability of a partially premixed burner using les[END_REF]. Moreover, to establish a successful flame growth and stabilization, the spark region must offer close to stoichiometric mixtures and flow velocities which do not avert an upstream propagation of the flame [START_REF] Cordier | Laser-induced spark ignition of premixed confined swirled flames[END_REF][START_REF] Sitte | Simulations and experiments on the ignition probability in turbulent premixed bluff-body flames[END_REF]. When fuel is injected as a droplet spray, ignition and propagation may notably differ from the processes corresponding to purely gaseous flames. Concerning the ignition of spray flames, numerical studies performed in [START_REF] Neophytou | The internal structure of igniting turbulent sprays as revealed by complex chemistry DNS[END_REF] highlight differences in the initial kernel structure as well as in the different possible modes of ignition. Flame propagation was found to depend on the droplet size distribution and the amount of prevaporized fuel [START_REF] Ballal | Flame propagation in heterogeneous mixtures of fuel droplets, fuel vapor and air[END_REF]. This is confirmed for instance in [12] in which distinct propagation modes are identified. Moreover, an inverse correlation of the burning velocity with the Sauter Mean Diameter (SMD) is also demonstrated by Neophytou 1 Töpperwien et al.

and coworkers [START_REF] Neophytou | Spark ignition of turbulent recirculating nonpremixed gas and spray flames: A model for predicting ignition probability[END_REF] in a bluff-body configuration. This effect was primarily attributed to the increased time needed for the evaporation of large droplets. Experimental data of light-round ignition in laboratoryscale annular combustors has become available more recently for premixed gaseous [START_REF] Bourgouin | Ignition dynamics of an annular combustor equipped with multiple swirling injectors[END_REF][START_REF] Machover | Experimental investigation on spark ignition of annular premixed combustors[END_REF] and non-premixed [START_REF] Machover | Spark ignition of annular non-premixed combustors[END_REF] systems. Burner spacing and overall equivalence ratio were identified to govern the propagation speed as well as the flow field resulting from the volumetric expansion across the flame. Numerical simulations of premixed [START_REF] Philip | Large Eddy Simulations of the ignition sequence of an annular multiple-injector combustor[END_REF] and nonpremixed [START_REF] Machover | Numerical Investigation of the Stochastic Behavior of Light-Round in Annular Non-Premixed Combustors[END_REF] annular configurations have proven, among others, that Large Eddy Simulations (LES) suitably retrieve experimental light-round trends. In the case of liquid fuel injection as a spray, numerical simulations have been successfully validated against experimental data becoming available recently for swirled annular spray-combustors using an Euler-Euler approach [START_REF] Boileau | LES of an ignition sequence in a gas turbine engine[END_REF][START_REF] Lancien | Large Eddy Simulation of Light-Round in an Annular Combustor With Liquid Spray Injection and Comparison With Experiments[END_REF] or Euler-Lagrange [START_REF]Modeling and Large Eddy Simulation of Two-Phase Ignition in Gas Turbines[END_REF] approach for the liquid phase.

A recent work by Prieur et al. [START_REF] Prieur | Flame and spray dynamics during the light-round process in an annular system equipped with multiple swirl spray injectors[END_REF] offers experimental insights into the dynamic behavior of spray flames during the ignition transient which are less well documented in the previous literature. These recent data are of particular interest as the flame is found to flashback under certain conditions and might damage the injection unit eventually. Flame flashback has been known to occur, for example, due to thermoacoustic instabilities [START_REF] Keller | Mechanism of instabilities in turbulent combustion leading to flashback[END_REF][START_REF] Lieuwen | Fuel Flexibility Influences on Premixed Combustor Blowout, Flashback, Autoignition and Instability[END_REF][START_REF] Schulz | Thermoacoustic instability in a sequential combustor: Large eddy simulation and experiments[END_REF] or transition between combustion regimes [START_REF] Sommerer | Large eddy simulation and experimental study of flashback and blow-off in a lean partially premixed swirled burner[END_REF]. However, experimental measurements reported in [START_REF] Prieur | Flame and spray dynamics during the light-round process in an annular system equipped with multiple swirl spray injectors[END_REF] suggest that flashback events are also possible during the combustor ignition sequence.

The present investigation relies on experimental observations made in a single-injector configuration, the SICCA chamber, operated at the EM2C laboratory. High-speed imaging revealed that the flame penetrates the injector and is initially present in this unit. After a delay, the lower flame brush is quenched, initiating a flame lift-off before stabilization at a distance from the injector. During the same period, the flame is also observed to change its shape.

The objective of the present work is to complement the experimental data and gain further insight into the mechanisms that drive the dynamic behavior which has not been covered in previous numerical studies and cannot be readily measured experimentally. This is achieved by making use of LES.

The outline of this work is as follows: the setup of the SICCA combustor is briefly presented, followed by an overview of the employed models and the numerical setup. Results of the LES are analyzed subsequently and compared with available experimental data.

EXPERIMENTAL CONFIGURATION

The SICCA-combustor (Fig. 1) is a swirled, singleinjector spray-combustor operated at ambient pressure. It is designed to facilitate the investigation of ignition phenomena occurring in aeronautical annular combustion chambers, allowing to characterize the flow aerodynamics, spray structure and flame geometry. The combustor features a pressurized simplex atomizer which feeds liquid fuel at 6 bar into the swirler (measured swirl number S = 0.68) where it mixes with air injected further upstream at the air inlets. The injector end piece has a conical shape with an outlet diameter of d = 8 mm. The fuel injector is mounted 6 mm in recess with respect to the injector outlet. In this configuration, a hollowcone spray is released into the combustion chamber which is confined by a cylindrical quartz tube of 70 mm inner diameter, 150 mm height and 5 mm thickness. Under stationary conditions, a flame stabilizes at a distance from the swirler outlet. For ignition, a conventional spark plug is used which is positioned close to the combustor back-plane and eccentric with respect to the combustor symmetry axis.

Various spray diagnostics have been performed as well as high-speed imaging of the ignition sequence which are detailed in [START_REF] Prieur | Flame and spray dynamics during the light-round process in an annular system equipped with multiple swirl spray injectors[END_REF]. Wall temperature measurements are available at the nominal, steady-state operating point of the combustor which are imposed in the simulations.

MODELING AND NUMERICAL SETUP

In the present two-phase flow approach, the governing equations for each phase are integrated by the explicit flow solver AVBP (www.cerfacs.fr/avbp7x) [START_REF] Colin | Development of High-Order Taylor-Galerkin Schemes for LES[END_REF] on an unstructured mesh. The most important numerical details and physical models are summarized in the following.

Gas phase

The filtered, compressible Navier-Stokes equations are solved for the gas phase. A two-step Taylor-Galerkin (TTGC, [START_REF] Colin | Development of High-Order Taylor-Galerkin Schemes for LES[END_REF]) scheme (third order accuracy in space and time) is used to solve the conservation equations numerically. Sub-grid scale contributions are computed following the classical eddy viscosity assumption. The turbulent eddy viscosity is evaluated according to the SIGMA-model [START_REF] Nicoud | Using singular values to build a subgrid-scale model for large eddy simulations[END_REF]. The turbulent species diffusivity and the turbulent heat conduction coefficient are determined from the turbulent Schmidt and Prandtl numbers (both equal to 0.6).

2

Töpperwien et al.

Dispersed phase

In this work, the liquid droplets are tracked in the Lagrangian framework. The governing equations for droplet motion are given as:

dx p,i dt = u p,i (1)
du p,i dt = f p,i m p (2)
with x p , u p and m p denoting the droplet position, velocity and mass, respectively. The two-way coupling between the polydisperse liquid phase and the gaseous phase causing a drag-force f p is taken into account by the Schiller-Naumann drag force model [START_REF] Schiller | A drag coefficient correlation[END_REF]. Phase interactions associated with evaporation are coupled in a similar manner via exchange terms in their respective mass and energy balances. Fuel addition to the gaseous phase and evaporative cooling are effectively represented by these terms. Unlike the conservation equations for the gaseous phase, a two-step Runge-Kutta scheme is used to integrate the dispersed phase equations which are coupled at every iteration to the gaseous solver.

Evaporation model

Due to droplet evaporation, additional equations are required to account for heat and mass transfer between the phases:

dm p dt = ṁp ( 3
)
dT p dt = - 1 m p c p,l φ c p ( 4
)
introducing the evaporation rate ṁp , the specific heat of the liquid c p,l and the conductive heat-flux φ c p . The well-known Spalding-model [START_REF] Spalding | The combustion of liquid fuels[END_REF] provides expressions for ṁp and heat fluxes in the gas phase φ c g which read as:

ṁp = -πd p Sh ρ g D g ln (1 + B M ) (5)
φ c g = -πd p Nu λ (T ∞ -T p ) ln (1 + B T ) B T (6)
where d p , Sh, ρ g , D g , Nu and λ denote the droplet diameter, the Sherwood number, the density of the gaseous mixture, the diffusion coefficient of the gas phase, the Nusselt number and the thermal conductivity in the gas phase. The vapor film which forms around an evaporating droplet and modifies the mass and thermal transport is accounted for by the Spalding numbers of diffusive mass and thermal transport, B M and B T , evaluated according to the method of Abramzon and Sirignano [START_REF] Abramzon | Droplet vaporisation model for spray combustion calculations[END_REF].

For the evaluation of B M under the assumption of thermodynamic equilibrium, the saturation pressure at the droplet surface is required which can be obtained from the Clausius-Clapeyron formula (pre-tabulated in AVBP). The Sherwood and Nusselt numbers follow from a correlation proposed by Frössling [32] which are functions of the droplet Reynolds number Re p , the Prandtl number Pr f ilm and the Schmidt number Sc f ilm in the vapor film. Its thermal conductivity follows from λ = µc p,g /Pr f ilm . These properties are commonly evaluated at a reference state, which can be interpolated between far-field and droplet surface conditions by the "third-rule" [START_REF] Hubbard | Droplet Evaporation: Effects of Transients and Variable Properties[END_REF].

Y f ilm f uel [-] 0.0 0.5 1.0 T norm [-] 0.7 1.2 Pr f ilm [-]
A balance equation at the phase interface yields a relation for the conductive heat fluxes of each phase, given as:

φ c p = -φ c g + ṁp L v (T p ) (7)
For the evaluation of this equation, the latent heat of vaporization at the droplet temperature is needed which is read from lookup tables. A critical task consists in the computation of the film properties D g , µ and λ which depend on Pr f ilm and Sc f ilm . This has been outlined in detail very recently by Sacomano Filho and co-workers [START_REF] Sacomano Filho | A novel strategy to accurately represent the carrier gas properties of droplets evaporating in a combustion environment[END_REF]. Considerable variations of these dimensionless numbers can be expected during the evaporation of a droplet. Reference data for droplet evaporation of n-heptane (droplet diameter d p = 15 µm, droplet temperature T p = 300 K) in pure, quiescent air computed with the EM2C in-house code AGATH are plotted in Fig. 2 (left). Each curve corresponds to a different surrounding air temperature to which the droplet is exposed. An increase of Pr f ilm can be observed with increasing air temperature while the evaporation time (indicated by markers) is decreasing. For Sc f ilm (not shown here), the trend with increasing surrounding air temperature is inverted. In addition, Pr f ilm appears to change over time in each case (the same holds for Sc f ilm ). Further analysis of the droplet film reveals that the film composition and film temperature are modified during evaporation, hence the observed effect on Pr f ilm and Sc f ilm .

3

Töpperwien et al.

To assess the variations of the film dimensionless numbers, generic droplet films have been systematically calculated in Fig. 2 (right). This map visualizes the variation of Pr f ilm due to changes of the film composition, represented by the film fuel mixture fraction (Y f ilm f uel ) and film temperature, plotted in a normalized form as T norm (normalization between T min = 280 K and T max = 2400 K). The film Prandtl number ranges between 0.7 and 1.2 within the given boundaries. Thus, an accurate description of droplet evaporation requires the computation of the droplet film properties which can be computationally expensive.

This has led to a method derived in [START_REF] Hannebique | Large eddy simulation of reactive two-phase flow in an aeronautical multipoint burner[END_REF] in which Pr f ilm and Sc f ilm are replaced by constants deduced from preliminary detailed 0D evaporation simulations. Very good results in terms of predicted evaporation time can be obtained in those cases in which the film properties can be expected to have little variation only, e.g. in stationary combustion.

During ignition, however, substantial variations in the droplet film can occur due to temperature and mixture fraction gradients which are encountered along the droplet trajectories and this may cause important discrepancies in the predicted droplet evaporation time. Therefore, an extended approach is proposed in this work. Instead of imposing constant values, Prandtl and Schmidt numbers are directly evaluated in AVBP as a function of the film temperature T f ilm and the film fuel mass fraction Y f ilm f uel . In this method, Pr f ilm and Sc f ilm of generic droplet films are fitted with a third-order, bivariate polynomial function f Y f ilm f uel , T norm , normalization of T f ilm using a least-square fit. Note that evaluating the film properties directly from detailed transport properties computation is also a solution but this would significantly increase the computational cost of the LES.

The corresponding fit coefficients are eventually implemented in AVBP and allow to use variable Prandtl and Schmidt numbers in the evaporation model as a function of the governing film properties. For efficient and numerically accurate evaluation of the polynomial function, Horner's rule [START_REF] Horner | A New Method of Solving Numerical Equations of All Orders by Continuous Approximation[END_REF] is applied.

To assess the accuracy of this approach, the error in droplet evaporation time t evap is plotted as a function of the equivalence ratio φ ∞ and temperature T ∞ of the (far-field) gas phase in which a droplet of n-heptane is evaporating. Very good results can be obtained for the droplet evaporation time in pure air (max. error 3 %, see Fig. 3 (left)) and a still acceptable error of max. 8 % is made in burnt gases at high temperature (Fig. 3 (right)).

Fuel injection model

For the injection of liquid fuel, the phenomenological FIM-UR approach (Fuel Injection Method by Upstream Reconstruction, [START_REF] Sanjosé | Fuel injection model for Euler-Euler and Euler-Lagrange large-eddy simulations of an evaporating spray inside an aeronautical combustor[END_REF]) is employed. The injected droplet size distribution is represented by a Rosin-Rammler distribution parameterized by a shape parameter q = 1.3 and a diameter d SMD = 18 µm to fit experimental data reported in [START_REF] Prieur | Combustion dynamics of an annular combustor with multiple spray injectors[END_REF]. The mean half-angle θ = 32 • is chosen to best fit experimen-1000 2000

T ∞ [K] 0.6 0.8 1.0 φ ∞ [-]
1000 2000 [START_REF] Paulhiac | Analysis of the spray flame structure in a lab-scale burner using large eddy simulation and discrete particle simulation[END_REF] is used in this work. Reaction rates following the Arrheniuslaw are corrected according to the methodology referred to as Pre-Exponential Adjustment (PEA) [START_REF] Franzelli | A two-step chemical scheme for kerosene-air premixed flames[END_REF]. A validation against detailed schemes [START_REF] Smallbone | Experimental and modeling study of laminar flame speed and non-premixed counterflow ignition of n-heptane[END_REF] has been performed in [START_REF] Paulhiac | Analysis of the spray flame structure in a lab-scale burner using large eddy simulation and discrete particle simulation[END_REF] proving a reasonably good prediction of the laminar burning velocity at ambient pressure and a wide range of equivalence ratios.

T ∞ [K] 0.6 0.8 1.0 φ ∞ [-] 4
The interaction of the flame with turbulence on a subgrid-scale level is controlled by the Dynamic Thickened Flame Model (TF-LES, [START_REF] Légier | Dynamically thickened flame LES model for premixed and non-premixed turbulent combustion[END_REF]). While this allows for a thin flame front to be resolved on a grid that is coarser than the flame thickness, an efficiency function E [START_REF] Charlette | A power-law wrinkling model for LES of premixed turbulent combustion: Part I -non-dynamic formulation and initial tests[END_REF] with constant β = 0.5 is required to compensate for the decreased sensitivity of the thickened flame to subgrid-scale eddies. Note that due to the presence of a liquid phase, thickening is also applied on drag and evaporation [START_REF] Paulhiac | Analysis of the spray flame structure in a lab-scale burner using large eddy simulation and discrete particle simulation[END_REF]. A flame sensor based on a lagrangian flame-front tracking approach as proposed in [START_REF] Rochette | A generic and self-adapting method for flame detection and thickening in the thickened flame model[END_REF] is employed.

Numerical domain and boundary conditions

Figure 4 shows the numerical domain, a vertical meshcut of the swirler and the lower part of the combustion chamber as well as a cut-plane perpendicular to the swirler. This mesh has already been used for validation purposes in a nonreacting as well as in a reactive, steady-state operation [START_REF]Modeling and Large Eddy Simulation of Two-Phase Ignition in Gas Turbines[END_REF]. The mesh covers the entire air plenum, the swirler, the com-4

Töpperwien et al. bustion chamber and a large cylinder representing the surrounding atmosphere (not shown in Fig. 4). Typical mesh sizes in the flame region are of the order of 0.4 mm. For inlets and outlets, the Navier-Stokes Characteristic Boundary Conditions (NSCBC, [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF]) are applied. Mass flow rates are imposed at the air inlet ( ṁair = 1.94 g.s -1 ) and the fuel inlet ( ṁ f uel = 0.111 g.s -1 ) yielding a global equivalence ratio of φ glob = 0.856. All inlet temperatures are set to T = 300 K.

No-slip walls are assumed for the gas phase in the entire domain with the only exception being the swirler walls for which a wall-law is used. Liquid phase boundary conditions are set to allow elastic rebound in the combustion chamber. In the injector, droplet-wall interactions are considered to be predominant and require a film-type boundary condition.

Previous works on the same combustor have highlighted the importance of an accurate representation of the thermal boundary conditions in terms of flame dynamics and flame shapes [START_REF] Lancien | Etude numérique de l'allumage diphasique de foyers annulaires multi-brûleurs[END_REF]. Therefore, the experimentally measured wall temperature profile is imposed on the inner combustor walls, while all other walls upstream of the fuel injector are adiabatic. The water-cooled combustor backplane is kept at a temperature of T = 410 K. These conditions can be considered to emulate relight conditions immediately after flameout.

Initial conditions and ignition procedure

The two-phase flow field has been converged by first running a non-reactive simulation to reach the targeted global equivalence ratio in the ignition zone of the chamber. In contrast to [START_REF]Modeling and Large Eddy Simulation of Two-Phase Ignition in Gas Turbines[END_REF], the wall temperature profile is imposed ab initio to ensure the convergence of the total heat flux across the chamber wall and backplane, i.e. steady-state conditions. This initial solution is eventually ignited using a 1-step energy deposition model (ED) as developed in [START_REF] Lacaze | Large eddy simulation of spark ignition in a turbulent methane jet[END_REF] to approximate the energy profile that a real spark would deliver, without simulating the initial plasma-phase. A volumetric source term is temporarily added to the energy equation featuring a Gaussian distribution in space and time. The center of the spherical energy source is positioned close to the com-

IGNITION IN RELIGHT CONDITIONS

Three ignition sequences have been recorded with a high-speed camera [START_REF] Prieur | Flame and spray dynamics during the light-round process in an annular system equipped with multiple swirl spray injectors[END_REF] at the same nominal conditions out of which one is used as a reference. Figure 5 (top) shows the experimental images which were averaged over 8.3 ms centered around the respective time instant. Two distinct flame configurations can be identified: during the first two instances, the flame is attached to the injector outlet which can be referred to as shape A (adopting the terminology in [START_REF] Prieur | Flame and spray dynamics during the light-round process in an annular system equipped with multiple swirl spray injectors[END_REF]). In the last instance, a flame detached from the swirler can be observed assuming shape B. Figure 5 (bottom) illustrates the simulated ignition sequence as a line-of-sight integration of the heat release rate inside the combustion chamber and the swirler nozzle. The same temporal averaging has been performed on the numerical data as in the experimental case. Similar flame configurations are found numerically supporting the argument of a dynamic flame behavior during ignition. The numerical images clearly reveal a flashback during the first two instances (shape A) as the lower flame brush is anchored inside the swirler nozzle which is particularly difficult to record experimentally without direct optical access. Conversely, in the following instance, the lower flame brush appears to be entirely extinguished. Therefore, complementary to the definition of the two flame shapes, the initial period will be defined as flashback phase (flame shape A) and the second period as lift-off phase (flame shape B) in the remainder of this work. Note that flashback and liftoff phases are considered to appear within the flame propagation phase (2) of the general definition of forced ignition mentioned in the introduction [START_REF] Lefebvre | Gas turbine combustion: alternative fuels and emissions[END_REF]. Each phase will be analyzed separately in the following sections. [START_REF] Ballal | The influence of flow parameters on minimum ignition energy and quenching distance[END_REF] Töpperwien et al.

Flashback phase

The volume integrated heat release rate inside the combustion chamber sampled during the ignition sequence is visualized in Fig. 6 (top, line with bullet markers). A short and intense peak of the heat release rate of Q ≈ 100 kW can be observed during the consumption of the fresh gases inside the combustion chamber and a fast decay to the nominal thermal power of the combustor at 5 kW. The rapid expansion of burnt gases induces a pressure perturbation that is propagating through the combustor and, more importantly, upstream through the swirler into the air plenum as well. In Fig. 6 (top) the pressure signal fluctuation P is monitored at the swirler inlet and the center of the combustion chamber. As the peak duration of these two signals overlap almost perfectly during the increase of the heat release rate, the effect of the burnt gas expansion on the pressure field can be considered to be immediate inside the entire combustor.

This relation between pressure and heat release is important for the analysis of mass flow rates which are recorded at the injector outlet and upstream in the air plenum in Fig. 6 (bottom). In fact, the pressure perturbation during the burnt gas expansion modulates the mass flow rates upstream of the combustion chamber, that is, at the injector exit (fresh gas mass flow rate) and in the air plenum (pure air mass flow rate): during sustained high pressure fluctuation levels, mass flow rates decrease substantially (cf. light gray areas in Fig. 6). Due to the temporary decrease of the mass flow rates the lower flame brush is able to propagate in the upstream direction and penetrates into the injector which is shown in detail in Fig. 7 (top) in the central cut-plane of the combustor. However, the lower flame brush detaches from the injector immediately after that moment, since the mass flow rates overshoot their respective average values at t = 3.9 ms (peak immediately after the gray area). Upon the subsequent drop of the mass flow rates due to another pressure peak at t = 4.6 ms (marked by a cross in Fig. 6) the flame is able to re-enter the injector up to the atomizer and find a stable configuration inside this unit. This process is illustrated in Fig. 7 (bottom).

In order to understand why this flame configuration remains stable despite further fluctuations of pressure and mass flow rates, it is instructive to examine the flow field and temperature distribution in the Internal Recirculation Zone (IRZ) as shown in Fig. 8 corresponding to the central cut-plane of the combustor. The flame front is visualized as black contour lines of heat release and the IRZ is surrounded by gray contour lines of zero axial velocity emerging from the swirler exit. The left image corresponds to a snapshot at t = 3.3 ms (sustained peak pressure), the right image to the reattached flame configuration at t = 4.6 ms discussed earlier. In the second image, the IRZ entrains burnt gases at higher temperature compared to the preceding snapshot. This leads to high temperature combustion products being fed into the lower flame brush facilitating flame stabilization inside the injector.

With this result in mind, the question might be raised, whether the energy profile introduced in the simulation, the ignition method itself or even the presence of a liquid phase might intensify the pressure excursion in an inappropriate manner leading to a flame flashback which would not have appeared otherwise. While experimental data of a purely gaseous configuration prove that the flashback process also accompanies ignition independently of the presence of the liquid phase (see Appendix of [START_REF] Prieur | Flame and spray dynamics during the light-round process in an annular system equipped with multiple swirl spray injectors[END_REF]), additional data can also be provided to assess the impact of the ignition method. Instead of using the energy deposition model for ignition, a simplified approach can be used by locally replacing the fresh gas mixture at the position of the igniter by burnt gases. In order to compare this method to the energy deposition approach, a separate simulation has been carried out from the exact same initial conditions. Ignition is achieved by introducing a sphere of burnt gases with the same dimensions (d = 11 mm) at the same position as in the energy deposition case to preserve consistency. In this manner, a spherical heat release profile is initialized which corresponds to a laminar, adiabatic flame at the local composition. The resulting simulated evolution of the heat release rate profile as well as the corresponding pressure signal at the center of the combustion chamber is plotted in Fig. 9. The corresponding profiles of the energy deposition method (already shown in Fig. 6) are added to facilitate the comparison. As expected, the volume integrated heat release rate profiles are quite similar in both cases and induce similar pressure signals, although with a steeper initial transient in the case of the energy deposition model. Nevertheless, the same dynamical behavior (flashback and subsequent lift-off) can be observed (not shown here) even with a simplified ignition approach. Thus, the major conclusion that can be drawn from the previous discussion is that the flame flashback occurs independently from the kernel initialization method. In the context of this work, the main purpose of the initialization method is to deliver an appropriate flame kernel which does not require an explicit modeling of the very first flame initiation after the plasmaphase to reproduce the pressure evolution in the combustion chamber which eventually causes the flashback.

Lift-off phase

After the flame stabilization during the flashback phase, this configuration prevails for a total of 14 ms before transitioning into shape B. This transition is shown in the central cut-plane of the combustor in Fig. 10. Starting from the anchored position (Fig. 10 (left)), the flame splits into an upper and lower flame branch (Fig. 10 (center)). The lower branch shrinks and completely vanishes as the transition into shape B takes place. The flame is finally detached from the swirler (Fig. 10 (right)).

The characteristic delay time before the onset of this process is relatively long and may be linked to the cooling of the burnt gases that circulate in the vicinity of the combustor walls as will be shown in what follows. Between t = 8.7 ms and t = 15.4 ms (cf. light gray area in Fig. 11), the burnt gas temperature in the Outer Recirculation Zone (ORZ) is substantially reduced from 2400 K to around 1500 K due to cooling by the combustor walls. By contrast, the temperature in the IRZ remains nearly constant until the lower flame Fig. 6: Top: volume integral of chamber heat release rate (solid line with bullet markers) and pressure signal sampled at the center of the combustor (dash dotted line) and at the injector outlet (solid line). Bottom: mass flow rates in the simulation at the injector outlet (thick line) and inside the air plenum (thin dashed line).

2.9 ms 2.9 ms 2.9 ms 2.9 ms 2.9 ms 2.9 ms 3 ms 3 ms 3 ms 3 ms 3 ms 3 ms 3.2 ms 3.2 ms 3.2 ms 3.2 ms 3.2 ms 3.2 ms 3.6 ms 3.6 ms 3.6 ms 3.6 ms 3.6 ms 3.6 ms 3.9 ms 3.9 ms 3.9 ms 3.9 ms 3.9 ms 3.9 ms branch splits and quenches resulting in a steep temperature decrease in the IRZ (cf. dark gray area in Fig. 11).

However, a heat loss effect in the IRZ which is suggested here as the driving mechanism for the flame lift-off cannot be directly inferred from the temperature signals because no differentiation can be made locally between fresh and burnt gases. For this reason, a qualitative index is introduced to visualize the heat losses from the burnt gases. This index I HL is defined as a normalized enthalpy, written as:

I HL = h -h l h a -h l ∈ [0, 1] ( 8
)
where h is the local enthalpy value, h a the adiabatic enthalpy at the local composition and h l the enthalpy of the local composition at the lowest temperature in the chamber, which is deliberately set to the lowest wall temperature inside the combustion chamber T = 410 K. When I HL is unity, adiabatic conditions prevail in the local composition, while 12.9 ms 12.9 ms 12.9 ms 13. I HL = 0 refers to the lowest enthalpy that a given composition may reach due to heat transfer to the chamber walls (i.e. highest heat loss). The heat loss index I HL distribution in the central cut-plane through the combustor is shown in Fig. 12 after stabilization of the flame inside the injector (a) and during the transition from shape A to B (b-d). The flame is represented by black contour lines of iso-heat release rate, whereas gray contour lines mark regions of zero axial velocity. Note that the value range of I HL is intentionally clipped between 0 and 1, although the fresh gas mixture before ignition can take values of I HL > 1. This is only due to the fact that the fresh gases are heated by the combustor walls in the absence of combustion, resulting in higher enthalpy levels compared to a reference value which is (deliberately) set to the lowest temperature in the combustor. The same applies for burnt gases immediately after ignition and flame stabilization which is shown in Fig. 12 conditions in the center of the combustion chamber which are not yet impacted by heat losses at this stage, near-wall fluid is constantly cooled by the combustor walls and I HL is close to zero in this region. In the ORZ, flow velocities are comparably low (on the order of 15 m.s -1 ) leading to a long exposure time of burnt gases to the walls. A substantial decrease of enthalpy starting in the ORZ and expanding downstream along the chamber walls can therefore be observed (Fig. 12(b-d)). Fig. 13: Gaseous equivalence ratio φ g sampled inside the injector above the fuel atomizer at x = -0.005 mm relative to the chamber backplane.

The flame is affected in the same manner, since cooled burnt gases are transported back upstream to the flame zone by the ORZ. Mixing again with fresh gases then results in a lower enthalpy of the combustion products leading to I HL < 1 even in the flame zone.

It may also be possible that the flame lift-off is locally enhanced by evaporative cooling of the gas phase close to the atomizer. The temporal evolution of the gaseous equivalence ratio φ g inside the injection unit in Fig. 13 suggests that when the flame anchors inside the injector, evaporation is intensified causing a substantial increase of φ g . Consequently, the gas phase locally incurs additional heat losses (cf. inset in Fig. 12(a)) and this may accelerate the extinction of the lower flame brush. However, evaporative cooling does not seem to have a first order impact on the flame lift-off as proven experimentally in a purely gaseous configuration (cf. [START_REF] Prieur | Flame and spray dynamics during the light-round process in an annular system equipped with multiple swirl spray injectors[END_REF]Appendix]), but may alter the transition time from the flashback to the lift-off phase.

Therefore, the flame lift-off can be explained more effectively by considering wall heat-losses: the IRZ which feeds the lower flame brush and sustains the initial flashback entrains burnt gases of decreasing enthalpy. As soon as the enthalpy is not sufficient to maintain a steady flame, the lower flame branch disintegrates and quenches entirely. The flame then stabilizes with a leading edge detached from the swirler.

A final discussion is dedicated to the impact of the twostep reaction mechanism on the prediction of quenching of the lower flame brush quenching. It is known that reduced kinetic mechanisms are less sensitive to strain rate and that in turn may delay quenching and lift-off. This issue is assessed by comparing the different flame regimes for the employed two-step mechanism and for a detailed mechanism for n-heptane (65 species, 315 reaction [START_REF] Smallbone | Experimental and modeling study of laminar flame speed and non-premixed counterflow ignition of n-heptane[END_REF]).

Data is obtained by performing steady-state 1D simulations of stabilized counterflow flames with the code AGATH. The flow configuration is chosen such that the fuel side represents the fresh gas composition at φ glob while the opposing stream side is formed by the corresponding burnt gas composition. Each 1D simulation is performed for a different set of imposed strain rate a and burnt gas temperature T burnt . In the resulting regime diagram plotted in Fig. 14, quenching occurs in the hatched region for the detailed reaction mechanism whereas the dotted region corresponds to quenching predicted by the two-step reaction mechanism. Conversely, a steady flame can be established in the region marked as burning regime. This comparison may be used to assess the flame lift-off observed in the SICCA-combustor in terms of the employed combustion chemistry. While both schemes show similar trends for flame extinction at low burnt gas temperatures and high strain rates, the quenched regime of the detailed mechanism extends towards higher temperatures and lower strain rates than the two-step mechanism. Thus, as the temperature of the IRZ decreases over time (cf. Fig. 11) the lower flame brush might extinguish at an earlier stage with detailed chemistry compared to the two-step mechanism.

CONCLUSION

The present work is focused on the transient dynamical processes that occur during the ignition of a swirling injector in a well-controlled laboratory-scale configuration. The injector offers an idealized version of practical systems used in aero-engines. The investigation is motivated by prior experimental observations in a multi-burner annular combustion chamber, indicating that flame flashback occurs during ignition followed by a final stabilization as a flame detached from the injection unit. This behavior is also retrieved in a single sector combustor equipped with one injector and characterized experimentally. This single-injector geometry is stud-9

Töpperwien et al.

ied in this work using Large Eddy Simulations employing a polydisperse Euler-Lagrange formalism. The simulation suitably retrieves the two periods identified experimentally. During an initial period, the flame is attached to the chamber bottom with a lower flame brush anchoring inside the swirling injector. An analysis of pressure signals inside the combustion chamber and on the upstream side of the swirler indicates a strong transient perturbation of the pressure field which is related to the rapid expansion of the burnt gases inside the chamber. This pressure perturbation controls the mass flow rates of fresh gases at the injector outlet and even further upstream in the air plenum and leads to a temporary decrease of their respective levels. The flame is entrained by the internal recirculation zone which allows for an upstream propagation and eventually stabilizes inside the injector unit. This transient flashback during ignition highlights the potentially detrimental effect on the injector elements. The final stabilization in the form of a flame detached from swirler is found to be controlled by heat transfer taking place between recirculating burnt gases in the vicinity of the combustor walls. A heat-loss index defined as a normalized enthalpy calculated at different instants during the simulation reveals a continuous enthalpy decrease starting in the outer recirculation zone and gradually impacting the flame zone as the cooled burnt gases mix upstream with the fresh stream of reactants. Consequently, the internal recirculation zone progressively entrains colder combustion products and quenches the lower flame root that is initially located inside the injector. Finally, the flame detaches from the injection unit and takes its steady-state position with a leading edge separated from the swirler.

Fig. 1 :

 1 Fig. 1: Experimental setup of the SICCA combustor. Left: confined setup with a single injector unit. The atomizer is recessed with respect to the injector outlet. Right: unconfined setup for spray diagnostics with the atomizer mounted flush with the injector outlet.

Fig. 2 :

 2 Fig. 2: Left: Temporal evolution of Pr f ilm during evaporation of a droplet of n-heptane in pure air at varying air temperatures. Solid line: T air = 300 K, dashed: 1000 K, dash dotted line: 1500 K, dotted: 2000 K. Right: map of Pr f ilm calculated for generic droplet films.

Fig. 4 :

 4 Fig. 4: Left: numerical domain (without ambient atmosphere). Top right: central mesh-cut of the injector and lower part of the chamber. Bottom right: mesh-cut of the swirler.

Fig. 5 :

 5 Fig.5: Top: Experimental ignition sequence recorded by high-speed imaging (reprinted from[START_REF] Prieur | Flame and spray dynamics during the light-round process in an annular system equipped with multiple swirl spray injectors[END_REF]). Bottom: Simulated ignition sequence: line-of-sight integrated heat release rate.

Fig. 7 :

 7 Fig. 7: Instantaneous heat release rate on the central cutplane of the combustor during the initial sustained high pressure levels (top) and during the second drop of the mass flow rates (bottom).

Fig. 8 :

 8 Fig.8: Instantaneous gas temperature in the central axial plane through the combustion chamber. Black contour lines correspond to a heat release rate of 10 8 Jm -3 s -1 , gray contours to zero axial velocity.

Fig. 9 :

 9 Fig.9: Volume integrated heat release rate using the energy deposition model (subscript ED, solid line, bullet markers) and a simplified burnt gases approach (subscript BG, dashed line, diamond markers). The corresponding pressure evolutions at the center of the combustion chamber are represented by dash-dotted lines.

  Jm -3 s -1

Fig. 10 :

 10 Fig. 10: Instantaneous heat release rate in the central axial plane of the combustor during flame lift-off and transition from shape A to B.

Fig. 11 :

 11 Fig. 11: Temporal evolution of temperature recorded in the IRZ and ORZ. The IRZ temperature signal is low-pass filtered (linear phase filter) for clearer identification of the curves.

Fig. 12 :

 12 Fig. 12: Heat-loss index I HL in the central axial plane through the combustion chamber. Black contour lines correspond to a heat release rate of 4 × 10 8 Jm -3 s -1 , gray contours to zero axial velocity.

Fig. 14 :

 14 Fig. 14: Regime diagram of a counterflow flame. Thick solid lines mark the critical point of extinction for a given temperature for each reaction mechanism. Dashed and dotted line continuations indicate conditions where extinction and autoignition critical points collapse.

  Fig.3: A posteriori error assessment of the extended evaporation model as a function of temperature and equivalence ratio of the surrounding gas phase. Left: error in t evap for evaporation of n-heptane in pure air. Right: error in t evap for evaporation in burnt gases.tally measured velocity profiles in non-reacting conditions at 2.5 mm above the chamber backplane.Assuming a dilute spray regime, droplet-droplet interactions are neglected which is verified a posteriori. Secondary atomization at the swirler exit is not yet covered.

	0 Error in t evap [%] 1 2 3	0 2 4 6 8 10 Error in t evap [%]
		8
		6
		6
	Combustion model
	A two-step, six species (C 7 H 16 , CO 2 , CO, H 2 O, O 2 , N 2 )
	global reaction scheme for n-heptane as proposed in

Töpperwien et al.

Acknowledgements

This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 766264. Furthermore, the work was performed using HPC resources from GENCI-CINES (Grant 2018-A0042B10118) as well as from the "Mésocentre" computing center of CentraleSupélec and École Normale Supérieure Paris-Saclay supported by CNRS and Région Île-de-France (http://mesocentre.centralesupelec.fr/).


Yunyun Sun

Florian Vixège

Khuram Faraz

Simon Mendez

Franck Nicoud

Damien Garcia

Olivier Bernard
email: olivier.bernard@creatis.insa-lyon.fr.

Creatis

A pipeline for the generation of synthetic cardiac color Doppler

Keywords: Simulations Ultrasound simulator, Color Doppler, Flow imaging, Echocardiography

Color Doppler imaging is the modality of choice for simultaneous visualization of myocardium and intracavitary flow over a wide scan area. This visualization modality is subject to several sources of error, the main ones being aliasing and clutter. Mitigation of these artifacts is a major concern for better analysis of intracardiac flow. One option to address these issues is through simulations. In this paper, we present a numerical framework for generating clinical-like color Doppler imaging. Synthetic blood vector fields were obtained from a patientspecific computational fluid dynamics CFD model. Realistic texture and clutter artifacts were simulated from real clinical ultrasound cineloops. We simulated several scenarios highlighting the effects of i) flow acceleration, ii) wall clutter, and iii) transmit wavefronts, on Doppler velocities. As a comparison, an "ideal" color Doppler was also simulated, without these harmful effects. This synthetic dataset is made publicly available and can be used to evaluate the quality of Doppler estimation techniques. Besides, this approach can be seen as a first step towards the generation of comprehensive datasets for training neural networks to improve the quality of Doppler imaging.

I. INTRODUCTION

Doppler echocardiography is the clinical imaging modality of choice for flow analysis in the cardiac chambers. Several Doppler techniques can be used to assess velocities in ultrasound imaging: i) Continuous-wave Doppler (CWD) to estimate a wide range of velocities along a predefined ultrasound line; ii) Pulsed-wave Doppler (PWD) to estimate velocities within a small spatial volume; iii) Color Doppler imaging (CDI) to estimate velocities over a large scan area. Regardless of the mode used, Doppler velocities represent projections of the actual velocity vectors along the ultrasound lines, thus providing incomplete one-dimensional information of the blood velocity field. Spectral Doppler (CWD and PWD) is used to locally estimate blood velocity in specific situations, such as in the assessment of aortic stenosis severity [START_REF] Pibarot | Energy loss index in aortic stenosis[END_REF]. The ultrasound wave must be carefully aligned with the flow to obtain Doppler velocities closely related to the actual blood speeds. In contrast, CDI is generally used as a visualization modality to display simultaneously the myocardium and intracavitary flow, although at a relatively low frame rate.

CDI is subject to several sources of errors, the main ones being aliasing and clutter noise. Aliasing occurs whenever the velocities to be measured exceed the Nyquist limit, whose value depends on the pulse repetition frequency (PRF) and the pulse center frequency [START_REF] Pozniak | Spectral and color doppler artifacts[END_REF]. In the classic red-blue Doppler color map, when aliased, color-encoded velocities wrap around such that color information turns from red to blue or vice versa [START_REF] Terslev | Settings and artefacts relevant for doppler ultrasound in large vessel vasculitis[END_REF]. Clutter corresponds to different sources of noise that mainly affect the low-frequency content of the slowtime signals. In echocardiography, clutter is mostly due to the motion of highly reflective tissues that surround the blood cavity. Therefore, blood signals must be separated from clutter before flow estimation.

Several approaches have been proposed to remove aliasing [START_REF] Shahin | Cooperation of fuzzy segmentation operators for correction aliasing phenomenon in 3d color doppler imaging[END_REF], [START_REF] Muth | Unsupervised dealiasing and denoising of color-doppler data[END_REF] and reduce clutter noise [START_REF] Yu | Eigen-based clutter filter design for ultrasound color flow imaging: a review[END_REF], [START_REF] Fadnes | Clutter filtering influence on blood velocity estimation using speckle tracking[END_REF]. These techniques are mostly based on standard signal and image processing, which may lead to inherent difficulties in ambiguous and non trivial situations [START_REF] Nahas | A deep learning approach to resolve aliasing artifacts in ultrasound color flow imaging[END_REF]. A recent study reports a deep learning (DL) solution to address the problem of aliasing removal [START_REF] Nahas | A deep learning approach to resolve aliasing artifacts in ultrasound color flow imaging[END_REF]. The proposed method involved a conventional U-Net architecture [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF] to localize aliased and double-aliased regions based on information derived from I/Q signals (i.e. Doppler frequency, power, and bandwidth, speckle flow angle, and speckle flow speed). Once detected, the aliased regions were corrected using a phase unwrapping algorithm as in [START_REF] Muth | Unsupervised dealiasing and denoising of color-doppler data[END_REF]. Since the input of the U-Net method uses physical information, it was necessary to apply a clutter filter at the beginning of the workflow. These results illustrate the capability of DL to segment regions with single or double aliasing. Combined with classical signal processing tools, this DL approach demonstrated its ability to improve the quality of flow visualization by CDI.

DL methods have revolutionized several fields of application in ultrasound imaging, the main ones being related to classification (e.g. automatic view recognition [START_REF] Madani | Fast and accurate view classification of echocardiograms using deep learning[END_REF]) and segmentation (e.g. automatic extraction of clinical indices [START_REF] Leclerc | Deep learning for segmentation using an open large-scale dataset in 2d echocardiography[END_REF]). In light of these successes, one may wonder why no additional work has been done so far to improve CDI. One of the main reasons is probably the difficulty of accessing referenced datasets, i.e. data whose flows are known and which can be used as a reference for learning algorithms. In this context, ultrasound simulations can play a key role in building such datasets. In this regard, our recent pilot study has shown that DL techniques can learn from synthetic ultrasound sequences to improve a targeted task; in this case, the estimation of rotational movements on in vitro data [START_REF] Evain | A pilot study on convolutional neural networks for motion estimation from ultrasound images[END_REF].

Computational ultrasound imaging has been applied for the analysis of the myocardial motion from B-mode sequences [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF], [START_REF] De Craene | 3d strain assessment in ultrasound (straus): A synthetic comparison of five tracking methodologies[END_REF], [START_REF] Alessandrini | A pipeline for the generation of realistic 3d synthetic echocardiographic sequences: Methodology and open-access database[END_REF], [START_REF] Zhou | A framework for the generation of realistic synthetic cardiac ultrasound and magnetic resonance imaging sequences from the same virtual patients[END_REF], [START_REF] Alessandrini | Realistic vendor-specific synthetic ultrasound data for quality assurance of 2-d speckle tracking echocardiography: Simulation pipeline and open access database[END_REF]. Most of the studies involved an electromechanical model (E/M) to generate groundtruth myocardial deformation fields. The E/M model has the advantage of relating cardiac contraction to its biophysical causes, allowing the inclusion of physiological and pathological conditions [START_REF] Marchesseau | Fast parameter calibration of a cardiac electromechanical model from medical images based on the unscented transform[END_REF]. Prakosa et al. first proposed a warpingbased method where a real B-mode sequence was registered to the E/M simulations [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF]. Although the resulting images showed realistic textures, warping defects degraded the quality of the synthetic images. To reduce those warping artifacts, De Craene et al. [START_REF] De Craene | 3d strain assessment in ultrasound (straus): A synthetic comparison of five tracking methodologies[END_REF] used the ultrasound simulator FieldII to generate ultrasound images. However, the scatterers were too simplistic, yielding unrealistic binary-like synthetic ultrasound images. To solve this issue, Alessandrini et al. combined the E/M model, an ultrasound simulator, and B-mode template images into a unified pipeline [START_REF] Alessandrini | A pipeline for the generation of realistic 3d synthetic echocardiographic sequences: Methodology and open-access database[END_REF]. Their simulation pipeline combined the advantages of both [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF] and [START_REF] Alessandrini | A pipeline for the generation of realistic 3d synthetic echocardiographic sequences: Methodology and open-access database[END_REF], and was able to produce realistic synthetic images. However, the backscattering transition between the myocardium and the surrounding environment was not sufficiently controlled, reducing the degree of realism of the simulated sequences. To solve this issue, the authors of two recent studies ( [START_REF] Zhou | A framework for the generation of realistic synthetic cardiac ultrasound and magnetic resonance imaging sequences from the same virtual patients[END_REF] and [START_REF] Alessandrini | Realistic vendor-specific synthetic ultrasound data for quality assurance of 2-d speckle tracking echocardiography: Simulation pipeline and open access database[END_REF]) introduced different strategies to ensure a smoother transition between the myocardium and the blood pool, leading to realistic B-mode echocardiographic imaging simulations (both 2D and 3D). The use of the electromechanical E/M model remains, however, the Achilles heel if a wide spectrum of simulations is required.

A. Main contributions

In comparison with the literature, the main novelties introduced in this paper are:

• the development of a simulation pipeline for the generation of clinical-like duplex ultrasound sequences, i.e. the simulation of both B-mode and Doppler images. • the inclusion of a patient-specific computational fluid dynamics (CFD) model into our simulations to produce ground-truth blood vector fields. • the introduction of controlled clutter noise through simple modeling of myocardial motion. • the assessment of the validity of our pipeline by simulating different synthetic duplex sequences.

II. METHODOLOGY

Fig. 1 provides the flowchart of the color Doppler simulation pipeline. In the following, we call "template" any real (clinical) duplex acquisition. The next sections (Sec. II-A, Sec. II-B, Sec. II-C) describe the content of each block. Clinical apical long-axis 3-chamber B-mode recordings are used to simulate sequences with realistic tissue texture. The first step of the simulation pipeline is to delineate the endocardium on the template cineloop. Two regions of interest (ROIs) are then obtained: the intraventicular cavity, and the left ventricular myocardium. Time-varying surface meshes are generated for both ROIs (block A in Fig. 1). Using these echo-derived meshes, a CFD intracardiac flow is anchored to the intracavity ROI of the template cineloop by spatio-temporal registration. It provides a ground-truth sequence of intracardiac blood motion (block B). This reference motion is then used to displace blood point-scatterers, inside the left ventricle, from one simulated frame to the next (block C). As for the myocardial meshes, they are used to move the myocardial scatterers to produce realistic B-mode images and clutter in the Doppler signals. An ultrasound open-source simulator (SIMUS) is finally run on the blood + myocardium scatterer maps to generate synthetic ultrasound signals. These signals are post-processed with conventional beamforming and Doppler algorithms to construct duplex ultrasound cineloops.

A. Pre-processing 1) Template image sequences: A clinical 2-D ultrasound acquisition was used as a template sequence to set the geometry and provide the texture of the simulated ultrasound images. The 2-D sequences (one complete cardiac cycle between two end-diastoles ED1 and ED2) were acquired from an apical three-chamber view using a GE Vivid e9 ultrasound scanner (GE Healthcare) and a 2.5 MHz phased array. All sequences were acquired in duplex mode to obtain both B-mode and Doppler images during a routine clinical examination. The ultrasound machine settings were adjusted to scan the entire left ventricular cavity and myocardium, from the apex to the base. It returned a pulse repetition frequency (PRF) of ∼7000, and a frame rate of 16 fps. All the subjects were scanned with nearly identical conditions.

2) Manual segmentation: The template sequences were segmented manually to get contours over the cardiac cycle. We defined eight instants evenly spaced between ED1 and ED2. Endocardial and epicardial borders were delineated manually at each instant. Each boundary was then resampled to ensure an even distribution of points along the contours. These points were finally interpolated in time using periodic spline to obtain myocardial and intracavity ROIs over the entire cardiac cycle. Once the endocardial and epicardial contours have been manually segmented, the myocardial and left ventricular meshing procedures are fully automatic.

3) Myocardial meshes: The left side of Fig. 2 illustrates the myocardial meshing. For each myocardial ROI, the base of the left ventricle was defined by the line linking point #1 (pt1) to point #2 (pt2). We defined the apex as the point of the epicardial contour whose distance from the base was greatest. 36 points were then evenly distributed over the epicardial contour: 18 on the lateral wall, and 18 on the anterior wall. Intramyocardial segments perpendicular to the epicardial contour and passing through the 36 epicardial points were then automatically drawn to join the epicardial and endocardial contours. Each intramyocardial segment contained 5 evenly distributed points. This resampling scheme allowed myocardial meshing with 180 points (36 longitudinal × 5 radial) and 280 triangle cells. These myocardial meshes were used in Sec. II-C to compute the myocardial displacements.

4) Left ventricle meshes:

The right side of Fig. 2 illustrates the left ventricular intracavitary meshing. As with the myocardial meshes, the basal and apical points were computed automatically from the intracavity ROIs. 21 points were evenly distributed over the endocardial contour. The basal and apical points were used to define the median axis of the left ventricular ROI cavity (see Fig. 2, right panel, gray curve). Two-thirds of the median axis were sampled uniformly, and the corresponding points (green points along the median axis) were associated with their endocardial counterpart to define linear 4-point segments. The apical region was sampled spherically. This resampling scheme allowed intracavitary meshing with 71 points (15 apical + 56 basal) and 114 triangle cells. These left ventricular intracavitary meshes were used in Sec. II-B to register a ground-truth blood flow m otion o n t he template sequences.

B. Blood flow motion 1) CFD model: To incorporate a reference blood-motion field into our simulations, we used a patient-specific heart flow CFD model developed by Chnafa et al. [START_REF] Chnafa | Image-based large-eddy simulation in a realistic left heart[END_REF]. In this CFD model, the cardiac cavities and the wall dynamics were extracted from 4-D clinical images acquired by computed tomography. The CFD computation was based on an ALE (Arbitrary Lagrangian-Eulerian) method. The ALE method is a way to couple Eulerian and Lagrangian formulations to solve fluid-structure coupling problems. This CFD model provided complex 3-D flow motion in a complete left heart including the left ventricle, the left atrium, and the aorta inlet. More details on the CFD model 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], [START_REF] Chnafa | Image-based simulations show important flow fluctuations in a normal left ventricle: What could be the implications?[END_REF]. From this model, an apical three-chamber view was reproduced by locating the probe at the apex, with the sector enclosing the mitral inlet and the left ventricular outflow tract (left part in Fig. 3), as in [START_REF] Assi | Intraventricular vector flow mapping-a doppler-based regularized problem with automatic model selection[END_REF]. In this study, we limited ourselves to planar simulations. Out-of-plane blood motions (in the elevation direction) of the insonified scan-area were not considered since the ultrasound simulations were twodimensional (see Sec. II-D).

2) Spatio-temporal registration: Fig. 3 illustrates the spatio-temporal registration we computed to anchor the CFD intracardiac flow to the template cineloop. We first used a piece-wise linear transformation to match the end-diastole and end-systole times of the CFD model with those of the template sequence. The orange dots in Fig. 3 were used as temporal landmarks to perform this temporal registration. Time-varying surface meshes for the intraventricular cavity of the CFD model were then generated using the same procedure as that described in Sec. II-A4. Each point of the CFD meshes was finally registered to the corresponding template meshes by interpolating the position of the corresponding cell. This procedure was used to spatio-temporally register one blood flow map per emitted beam on the template cineloop. It modified both the amplitudes and directions of the velocity vectors. The mechanical properties of the flow were therefore not totally conserved. For example, the conservation of mass was affected. Our goal here was to obtain clinical-like Doppler images from an imaging point of view, not from a mechanical point of view.

C. Myocardial motion

Our simulation pipeline includes the motion of the myocardium to generate synthetic clutter in left ventricular color Doppler. The myocardial motion was not simulated with the electromechanical model that we used in [START_REF] Alessandrini | Realistic vendor-specific synthetic ultrasound data for quality assurance of 2-d speckle tracking echocardiography: Simulation pipeline and open access database[END_REF]. For each simulation, a set of points was randomly distributed over the myocardial mesh at end-diastole. Each of these points was then propagated over the full sequence by interpolating the displacements of the corresponding cell. This simple procedure allowed us to compute applicable temporal trajectories of any point belonging to the myocardium.

D. Ultrasound duplex simulation

Our pipeline uses a homemade open-source software called SIMUS from the MUST Matlab ultrasound toolbox (https://www.biomecardio.com/MUST) [START_REF] Shahriari | Meshfree simulations of ultrasound vector flow imaging using smoothed particle hydrodynamics[END_REF]. SIMUS simulates backscattered ultrasonic signals for linear, phased, and convex arrays [START_REF] Garcia | Simus: an open-source simulator for ultrasound imaging. part i: theory & examples[END_REF]. It allows the simulation of ultrasound images from different transmission schemes (i.e. using focused beams, plane waves or diverging waves) by adjusting the transmit delay laws. The latest version of the simulator now integrates 3-D acoustics, which was not the case at the time of this study. Our ultrasound simulations were thus planar. The insonified medium is modeled by monopole point sources that do not interact with each other (weak scattering). From the position and reflection coefficient of each scatterer, SIMUS computes the ultrasound signals that are received back by the probe, as illustrated in block C of Fig. 1. We used SIMUS to compute the B-mode and Doppler cardiac images (duplex sequences). The same probe settings used to acquire the real B-mode sequence were simulated: a 2.5 MHz 64-element cardiac phased array emitted at a PRF of 7000 Hz. The bandwidths at -6dB were respectively 60% and 20% for the B-mode and Doppler sequences, which corresponded to transmit pulses of 2 and 6 wavelengths. Assuming a heart rate of 60 beats per minute, a focused-beam configuration returned sixteen 11-cm deep B-mode frames interleaved with color Doppler images obtained with packet lengths of 8.

1) Realistic B-mode sequence:

The overall strategy for the simulation of the B-mode sequence is illustrated in Fig. 4. Each B-mode frame was generated by transmitting 90 focused beams.

Background scatterer maps: N back scatterers were randomly distributed in the sector of the first template image. The scatterer density was 10 per square wavelength, which corresponded to N back ≈ 250, 000. To avoid flickering effects, the background scatterers were kept motionless. We set the reflection coefficients of the scatterers to generate realistic images. To mimic the tissue echogenicity of the recorded model, the local intensities I m of the actual B-mode images were used to calculate the reflection coefficients RC m of the scatterers, according to the following expression:

RC m = (I m /255) (1/γ) • N (0, 1) (1)
where N (•) is the normal distribution, and γ is a constant for gamma compression (set to 0.3).

Myocardial scatterer maps:

The scatterers of the myocardium were selected on the first simulated frame. The positions of these scatterers were then computed for each B-mode frame of the simulated sequence using the strategy described in Sec. II-C. The reflection coefficients of these scatterers were kept constant to maintain the speckle texture throughout the cardiac cycle.

Blood scatterer maps: To simulate Doppler information, blood scatterers were added in each B-mode frame of the simulated sequence. These scatterers were randomly distributed in the intracavity ROIs with a density of 10 scatterers per square wavelength and reflection coefficient values following a Gaussian distribution. To mimic the large difference in amplitude between blood and the surrounding tissues, we set a ratio of 70 dB between the RC values of the myocardial and blood scatterers RC blood m . The following expression was used:

RC blood m = 1 10 (70/20) • N (0, 1) (2)
Final scatterer maps: The final scatterer maps were obtained by combining the background, myocardial, and blood scatterer maps. The background and myocardial scatterers were combined to ensure a smooth transition at the myocardial borders using the same scheme as described in [START_REF] Alessandrini | Realistic vendor-specific synthetic ultrasound data for quality assurance of 2-d speckle tracking echocardiography: Simulation pipeline and open access database[END_REF].

2) Simulation of the Doppler information: Fig 5 illustrates the scheme we implemented to simulate Doppler images. Blood scatterers were moved between each emitted beam by interpolating their displacements through the template flow maps described in Sec. II-B. Each color Doppler image was generated using 44 focused beams in transmission with packet lengths of 8.

3) Reconstructed synthetic B-mode and Doppler data: The synthetic RF signals generated by SIMUS were demodulated to obtain in-phase/quadrature I/Q signals. The I/Q signals were beamformed using a delay-and-sum to obtain B-mode and color Doppler images [START_REF] Perrot | So you think you can das? a viewpoint on delay-and-sum beamforming[END_REF]. For clutter filtering, we removed the mean of the I/Q signals (= zero-order polynomialregression filter). The Doppler velocities were deduced by using a lag-1 autocorrelator [START_REF] Madiena | Color and vector flow imaging in parallel ultrasound with sub-nyquist sampling[END_REF].

E. Simulated scenarios

Our color Doppler simulation pipeline is flexible enough to generate different types of scenarios. This is interesting in the perspective of generating datasets to feed into DL algorithms to improve the quality of flow visualization and quantification by CDI. To illustrate the relevance of each element of the pipeline, we simulated four scenarios ranging from simplified to more realistic situations for five virtual patients. For each simulation, the generated synthetic sequence and the reference motion field are made available to the community, yielding to an open access dataset composed of 20 benchmarked sequences. This synthetic dataset is accessible at the following link: http://humanheart-project.creatis.insalyon.fr/duplex.html. The simulation of a full sequence composed by 16 B-mode and color Doppler images took around 4 hours using a C version of the SIMUS code running on an Intel 2 Xeon E5@2.9Ghz with 512GB RAM.

1) Scenario #1: color Doppler "snapshot": This scenario assumes that one would be able to obtain an instantaneous color Doppler image (color Doppler "snapshot") 1) in the absence of myocardial motion, and 2) assuming that the flow is stationary (no temporal acceleration) during the ultrasound emissions generated to obtain this image. The number of emissions to simulate a Doppler image was 352 (44 focused beams with a packet length of 8, as illustrated in Fig. 5). This prevented i) the introduction of clutter noise due to the movement of high-intensity tissues and ii) the presence of bias during Doppler estimation caused by blood acceleration. In this snapshot scenario, the estimated Doppler velocities should match the radial velocity components of the CFD-based reference flow.

2) Scenario #2: Doppler simulation of a non-steady-state flow: Blood motion is no longer quasi-stationary: the blood scatterers are relocated between two successive ultrasound transmits. This scenario takes flow acceleration into account, which has side effects on Doppler estimation. The motion of the blood scatterers was computed using the template flow maps described in Sec. II-D2.

3) Scenario #3: additional synthetic clutter noise: In contrast with the other scenarios, the third scenario aims at simulating color Doppler as obtained in a clinical setting. All the blood and tissue scatterers are relocated between two successive ultrasound transmits. In this context, the information carried by the slow motion of the tissues is mixed with that of the blood flow, which led to wall clutter noise. 4) Scenario #4: focused vs. diverging waves in transmission: Each color Doppler image presented in the previous scenarios was generated using 44 focused beams in transmission with packet lengths of 8. This led to an accurate estimate but at the cost of a small temporal frequency (e.g. 16 fps for the simulated sequences). To work at higher imaging rates, recent works proposed to estimate Doppler velocities from diverging waves [START_REF] Faurie | Intracardiac vortex dynamics by high-frame-rate doppler vortography-¡italic¿in vivo¡/italic¿ comparison with vector flow mapping and 4-d flow mri[END_REF], [START_REF] Faurie | Coupling myocardium and vortex dynamics in diverging-wave echocardiography[END_REF]. In these conditions, the whole Doppler sector is insonified with each emitted beam, which drastically increases the frame rate of color flow imaging. However, since this strategy spreads the acoustic energy into a large area, the quality of the Doppler estimate is degraded, especially in situations where high-intensity tissue motion is important. In this context, the goal of this scenario was to investigate the influence of the transmission modes on the quality of estimated velocities. Each sequence was simulated in the "snapshot" conditions with the emission of either focused or diverging beams for the estimation of Doppler information.

III. RESULTS

A. Scenario #1: color Doppler "snapshot"

The second line of Fig. 6 shows Doppler images simulated in the snapshot scenario for different instants during filling, when blood velocity was the highest. From these results, one can appreciate the reliability of the simulated sequence compared to the corresponding Doppler reference given by the radial velocities. Most of the reference flow patterns were simulated and estimated accurately. The small differences are in large part explained by smoothing effects introduced by the Doppler estimator. The simulated color Doppler images also contained aliasing artifacts when the velocities to be measured exceeded the Nyquist limit (1.1 m/s in our case).

B. Towards more realistic scenarios 1) Scenario #2: Doppler simulation of an unsteady flow:

The third line of Fig. 6 shows images simulated in such conditions for the same instants as those of the snapshot scenario. The Doppler outputs were slightly different since this scenario does not assume stationary flows. The effect of flow acceleration is more visible during the filling phase.

2) Scenario #3: additional synthetic clutter noise: Because the information conveyed by the motion of the tissues was mixed with that of the blood flow in this scenario, it was necessary to apply a clutter filter on the synthetic I/Q signals (after beamforming) to recover the blood signals. This was realized using a zero-order polynomial-regression filter. The forth line of Fig. 6 shows images simulated in such conditions. The Doppler outputs were significantly different since this scenario did not provide clutter-free instantaneous maps. The effect of the surrounding myocardial motion during the filling phase is well visible during filling onset. As tissue velocities decreased during the diastolic phase, clutter noise (along the ventricular wall) reduced until it vanished at the end of filling. This is in line with real observations, which illustrates the realism of the simulation pipeline. In addition, Fig. 7 shows the Doppler estimates during the systolic phase in the snapshot scenario compared with the wall-clutter scenario. These results further illustrate the impact of tissue motion on Doppler estimation. In this example, Doppler bias appeared along the blood/tissue interface, especially near the septum where tissue motion was most pronounced.

C. Scenario #4: focused vs. diverging waves in transmission

Fig. 8 displays the results obtained for two instants where the motion of high-intensity tissues was either high (i.e. beginning of diastole) or low (i.e. diastasis). These results show that the quality of Doppler estimation is related to the amount of Fig. 6: Illustration of the ability of our Doppler simulation pipeline to generate three different types of scenarios: 1) a color Doppler "snapshot" situation assuming a stationary flow and the absence of myocardial motion (second row); 2) unsteady flow where the blood scatterers are relocated between two successive ultrasound transmits (third row); 3) wall-clutter with unsteady flow where both the blood and tissue scatterers are relocated between two successive ultrasound transmits (fourth row). tissue motion when diverging waves are used in transmission. Indeed, the clutter phenomenon is more prominent when the motion of high-intensity tissues is important.

D. Evaluation of a Doppler estimation algorithm

Our simulation pipeline can also be used to provide benchmark datasets to assess the quality of Doppler estimation methods. We evaluated the effect of wall filtering in the presence of myocardium-based clutter by using two different datasets: one generated from scenario #1, and another from the same scenario but with additional wall clutter. Each dataset was composed by 16 Doppler images for 5 virtual patients. In each case, the estimated Doppler velocities were compared with the radial velocity components of the CFD-based reference flow. Fig. 9 displays the mean correlation plots computed over the five virtual patients. Results show a high agreement between the CFD-based and simulated-derived Doppler velocities for the dataset generated from scenario #1, with a coefficient of determination of 0.92. The small differences can be explained by the smoothing effect of the estimator. The degradation induced by the wall clutter noise can be investigated from the right part of Fig. 9. In particular, a decrease of 19% of the coefficient of determination can be observed.

IV. DISCUSSIONS

A. On the benefits of training DL methods through simulations Deep learning techniques have been successfully applied to echocardiography for several years. These methods have allowed major advances in many specific areas such as view classification and segmentation of anatomical structures. Almost all the best performing techniques correspond to supervised approaches, requiring the deployment of largescale datasets with manual annotations. As far as ultrasound motion estimation is concerned, very few referenced data are available, the main one being an open access dataset proposed for myocardial strain benchmarking [START_REF] Alessandrini | Realistic vendor-specific synthetic ultrasound data for quality assurance of 2-d speckle tracking echocardiography: Simulation pipeline and open access database[END_REF]. A supervised DL network was recently trained on this open-access dataset to automatically estimate myocardial strains [START_REF] Østvik | Myocardial function imaging in echocardiography using deep learning[END_REF]. Although the results were very promising, the lack in quantity and diversity of training data did not allow the authors to obtain significantly better scores than the state-of-the-art. Regarding blood flow estimation, there is currently no baseline dataset. Thus creating large-scale synthetic datasets that mimic hundreds or thousands of patients would be the key to train data-driven DL methods for flow estimation. In this context, we proposed the first simulation pipeline that generates clinical-like duplex ultrasound sequences. Although the size of the derived dataset is not sufficient to train DL models, our framework constitutes the first step toward achieving more complete datasets. A major weakness of our protocol is the need to manually segment the left ventricle and the myocardium throughout the sequence. In the near future, we plan to make this step fully automatic to enhance the scalability of our simulation pipeline, with the goal of creating simulated data that replicates a large number of patients.

B. An innovative simulation pipeline for duplex sequences

The objective of our study was to simulate clinical-like B-mode and color Doppler ultrasound images. By "clinicallike" we mean that these images should resemble those obtained by a cardiologist with a clinical ultrasound machine. For the B-mode images, a set of focused waves scanned the cardiac region, and then a standard delay-and-sum technique was used for image construction. Although the transmission/beamforming scheme may differ from those used by clinical ultrasound scanners, the similarity to clinical images is apparent, particularly with respect to speckle patterns. The realism is also found in the color Doppler simulations. Not only is the similarity observed on the Doppler fields, but also on classical artifact phenomena, such as aliasing and wall clutter.

In this study, we focused only on the generation of clutter due to myocardial motion. We neglected other sources of clutter such as those related to rib reverberations, multiechoes, or twinkling artifacts [START_REF] Nilsson | Artefacts in sonography and doppler[END_REF]. Our goal was not to perform a thorough evaluation but rather to show that our pipeline could generate datasets for the evaluation of Doppler estimation techniques. A thorough investigation is beyond the scope of this paper. The evaluation showed that the datasets can lead to accurate estimates of Doppler velocities in the absence of wall-based clutter. In addition, we quantified the decrease in color-Doppler performance when wall clutter was present. This is a direct consequence of the motion of highly reflective tissues that surround the blood cavity. This shows the difficulty in developing techniques that are robust to the clutter phenomenon. We believe that DL techniques could help to mitigate wall clutter, based on adequate simulations and appropriately designed references.

C. Potential improvements

The proposed pipeline is a relevant advance in the field of cardiac ultrasound simulation. Despite the variety of scenarios already provided, it would be of interest to integrate out-ofplane motions. Indeed, the 3-D nature of blood flow induces some decorrelation of the slow-time signals, which may affect the quality of Doppler estimates. In this study, our goal was to develop the first duplex simulation pipeline by focusing on key properties such as unsteady flow and wall clutter. The 3-D CFD model will allow us, in the future, to study and quantify the impairment caused by 3-D motion. We also evaluated the effectiveness of our simulation pipeline from healthy subjects. It would be interesting in the future to extend our framework to simulate large scale datasets including pathological cases. We will need to incorporate blood and tissue motions that reflect anatomical and functional characteristics of patients with heart disease. It should be noted however that some pathologies (e.g. valvular diseases, septal defect,...) cannot be considered due to the limitations of the CFD model.

D. Perspectives

The pipeline was used to generate an open access dataset of 20 synthetic sequences available at the following link: http://humanheart-project.creatis.insa-lyon.fr/duplex.html. Part of these data can already be used to benchmark the quality of Doppler estimation techniques. Our goal is to progressively enrich this dataset with new simulations to cover an exhaustive and diverse range of blood flow dynamics, from healthy to pathological cases. This will provide richer reference datasets that will help improve the quality of flow visualization and quantification by CDI using DL methods. We are targeting three main applications in the near future. Our framework will first be used to estimate alias-free Doppler velocities from the backscattered ultrasonic signals. We have shown that our pipeline can also generate wall clutter noise. Despite a large literature dedicated to this phenomenon, clutter remains one of the main sources of error in intraventricular flow estimation. Our solution therefore offers new perspectives to develop DL-based clutter filters more efficient than the state-of-theart (such as singular value decomposition filters). Finally, it has been shown that it was possible to decipher 2-D velocity vector fields using color Doppler images [START_REF] Garcia | Two-dimensional intraventricular flow mapping by digital processing conventional color-doppler echocardiography images[END_REF]. This can be used to characterize the vortex that forms in the left heart and investigate the filling capacity of the heart. The underlying method involves several tedious steps, such as segmentation of the left ventricle and aliasing/clutter filtering. We thus plan to investigate the application of DL for the automatic and robust estimation of two-dimensional intraventricular flows.

V. CONCLUSIONS

We developed a novel simulation pipeline for the generation of clinical-like duplex ultrasound sequences. The proposed pipeline combines state-of-the-art solutions in the fields of computational fluid dynamics and computational ultrasound imaging. The CFD model determines the reference value of the blood flow. The ultrasound simulator, in conjunction with scatterer strategies, ensures realistic speckle patterns and Doppler estimates. We evaluated the genericity of our pipeline in four different scenarios. In particular, we showed that our pipeline allows the integration of flow acceleration and wall clutter noise. This led to side effects and made Doppler estimation challenging, as in real cases. The pipeline was used to generate an open access dataset of 20 synthetic sequences. This synthetic dataset can already be used to benchmark the quality of Doppler estimation techniques.

Fig. 1 :

 1 Fig. 1: Proposed pipeline for the simulation of clinical-like duplex ultrasound sequences. (A) A clinical recording works as a template for speckle texture, anatomy definition, and myocardial motion estimation; (B) A CFD model controls the synthetic blood flow motion; (C) An ultrasound simulation environment merging information from the template image sequence and the CFD model accounts for the image formation process. In the simulated sequence, the blood flow is fully controlled by the CFD model while the visual appearance is very similar to the one of a real acquisition. Each block specifies the section number where the description is detailed.

Fig. 2 :

 2 Fig. 2: Meshing of the segmentation masks for the myocardium (left) and left ventricular cavity (right). The proposed schemes allowed myocardial meshing with 180 points (36 longitudinal × 5 radial) and 280 triangle cells, and intracavitary meshing with 71 points (15 apical + 56 basal) and 114 triangle cells.

Fig. 3 :

 3 Fig. 3: Illustration of the spatio-temporal alignment used to register the flow of the CFD model to the template sequence. The color map indicates the amplitude of the velocity along the cross-beam x-direction. The orange markers were used as key instants to perform a piece-wise linear registration between the CFD model and the template sequence. Each point of the CFD meshes was registered to the corresponding template meshes by interpolating the position of the corresponding cell.

Fig. 4 :

 4 Fig. 4: Overall strategy for the simulation of B-mode sequences. Each B-mode frame of the template sequence was used to define a set of scatterers characterized by their positions and reflection coefficients. This information was then provided into the SIMUS ultrasound simulator to generate synthetic images.

Fig. 5 :

 5 Fig. 5: Overall strategy for the simulation of Doppler images. 352 registered flow maps were used between each consecutive B-mode frame. This leads to a different map per emitted beam at the PRF frequency. 44 focused beams with a packet length of 8 were used to estimate each color Doppler image.

Fig. 7 :

 7 Fig. 7: Color Doppler snapshot vs. wall-clutter scenarios simulated during systolic phase. Large tissue motions along the septum generate synthetic clutter that degrades the color Doppler image estimated from the synthetic I/Q signals.

Fig. 8 :

 8 Fig. 8: Focused vs. diverging beams for the estimation of Doppler velocities. While the diverging wave strategy increases the temporal frequency, the spread of the acoustic energy into a large area degrades the quality of the Doppler estimate, especially in situations where the displacements of high-intensity tissues are significantly large.

Fig. 9 :

 9 Fig. 9: CFD-based vs, simulation-derived Doppler velocities (ν D ) for (a) the scenario #1 and (b) the scenario #1 with additional wall clutter noise. Velocity data from sixteen Doppler images for five different virtual patients.

ACKNOWLEDGMENT

This work was performed within the framework of the LABEX PRIMES (ANR-11-LABX-0063) of Universite ´ de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR).