H Gassot

Simulation developments for Spoke superconducting cavities design

In the frame of the HIPPI JRA, which proposed R&D activities meant to foster the development of a common European technology base for the construction of high intensity pulsed linear accelerator, IPN Orsay is involved in the development of the Spoke type superconducting cavities. The 352 MHz 2 gaps Spoke cavity has been proposed by IPN Orsay for beta 0.15 sections, and the 352 MHz 4 gaps Spoke cavity has been proposed by FJZ Jülich for beta 0.65 sections. The important simulation investigations have been realized at IPN Orsay, as the numerical codes adapted for superconducting cavity electromagnetic and mechanical behaviors' studies are different from FJZ laboratory, the simulation results supply confirmations about 4 gaps Spoke design made by FJZ Jülich.

Introduction

Superconducting (SC) RF cavities have much larger efficiency, accelerating gradient and bore aperture than normal conducting (NC) structures. This technology is then expected to be advantageous in a linear accelerator in terms of power consumption, construction cost and beam loss. Although this conclusion is well accepted for the high energy part of the accelerator, this is not the case at low energy, mostly because of the short distance required between transverse focusing magnets, which reduce the energy gain per real estate meter. So for a closed beta, several alternative solutions have been proposed and optimized. All the cavities developed in the HIPPI project are designed in the pulsed mode, the sensitivity to the Lorenz forces is especially critical. It is therefore of high importance to improve the simulation tools which play a really important role in the prediction of Lorentz forces detuning of superconducting cavities. It is thus of high importance to reduce as much as possible the effects of the Lorenz forces detuning for the cavities. The use of several elements in the design such as stiffening rings ribs or tuner to enhance the stiffness of the cavities require optimization and simulations controls. Finally, these studies are helpful to determine the lowest energy at which low beta superconducting cavities could safely and economically operate.

SIMULATION PLATFORM AT IPN ORSAY

1-1. MULTI-PHYSIC BEHAVIOUR

Superconducting RF cavities are formed of niobium walls in which a high electromagnetic field is established in order to provide a longitudinal electric field at the axis of the beam pipe. One of the principal advantages of this structure is its capacity to obtain a high accelerating gradient in a limited distance, minimizing the linear accelerator length. But as a consequence, the level of the electromagnetic field at the inner surface of the niobium wall can be very high; the induced Lorentz forces can become critical for the stability of the resonator. In fact, the Lorentz forces deform the cavity wall's structure, and consequently, the resonance frequency is shifted in a few milliseconds while the electromagnetic field establish instantaneously in the cavity, this phenomena can drive the RF cavity out of stability. For this reason, the prediction and the control of the Lorentz forces detuning is fundamental in the design of the superconducting cavity.

The most significant difficulty concerning the Lorentz forces detuning studies is the result of multi-disciplines behaviours involved and coupled. First, a preferment CAD tool is necessary for the design in all full 3D problems. Starting with the CAD computer design, the electromagnetic optimization must be performed in the vacuum space contained inside the cavity wall, simultaneously, the structure analysis should be performed on the cavity wall modelled with shell and/or volume elements, so the simulations are completed considering the wide aspects to cover and their interconnections.

The major work developed at IPN Orsay consists in creating a platform for the simulation tools, which offers an additional dimension to the traditional specific codes for each specific domain (Catia for CAD design, Cast3m for mechanic calculations and Soprano/Opera for electromagnetic simulations), makes the coupled simulations instantaneously and gives a global vision of the multi-physic problem. The choice of each specific simulation tool is based on the existent codes at IPN Orsay, this guarantees the reliability of some well know tools confirmed in past experiences and reduces the financial investment. Finally, the simulation results can cross-confirm the design studies realized with commercial codes like ANSYS and others.

1-2. PLATE-FORM'S CONFIGURATION

The most important investigation consists in developing a platform which allows the instantaneous linkages between CAD, mechanical simulations and RF simulations. The difficulty is to assume the compatibility of different models used for different physics problems: for electromagnetic, the simulation model is constructed on the cavity inner volume while the mechanical analysis model is constructed on the cavity wall. The mechanical simulation depends on the results of the electromagnetic simulation which determinates the electromagnetic field distribution -and thence forces -on the inner surface of the cavity wall. The crucial point is to transfer these simulation results obtained on the vacuum volume model to the cavity wall model.

This situation determinate the global strategy to make the coupling between mechanical simulations and the electromagnetic simulations: to avoid all problems in terms of the field interpolation or any lost of the numerical precision due to the field transfer between one code to another, the two models used in two different simulations codes should have the same modelling in the common surface which is at same time the inner surface for the cavity wall and the envelop of the vacuum space.

As an additional aspect, we can consider two kinds of model in structure analysis: the cavity wall could be modelled with a shell model under hypothesis of the thick wall or with a volume model. In the case where the simulation concerns only the cavity wall without any additional stiffener, the shell model is simpler to use than the volume model. But in case the stiffeners are necessary, the volume model is more adapted. The development of the global platform is based on an adequate solution which integrates multi-physic tools with all specific constrains; three groups of interface programmes have been developed:

• Interface between CAD code Catia and electromagnetic code Opera • Interface between CAD code Catia and mechanic code Cast3m

• Interface between electromagnetic code Opera and mechanic code Cast3m First, the cavity's design is drawn by Catia. In addition, Catia can make the structures modelling by finite elements. The vacuum space can be meshed by tetrahedrons. The advantage is the easiness of the modelling directly connected to the design: no problem for interpretation of the geometry; time saving. Opera 3D is a software for electromagnetic design. One of its modules Soprano is a specialized analysis module for finding the frequency domain fields, or resonant modes (Eigen values) of a structure. In addition to its own modeler, Opera 3D has a Pre Processor which is able to read external modeling written in IDEAS-Universal format.

The interfaces programmes assure the transfer of all geometry data and the connectivity of elements. More, the interfaces programmes create all surface elements coherent from the tetrahedrons volume modelling for the specific electromagnetic boundary conditions definition: prefect reflexion wall… Another group of the interfaces programmes consists to elaborate the compatible model for mechanical simulations with Cast3m code. The essential matter is to assure the coupling of the electromagnetic analysis with the structure analysis. Here, the interfaces programs offer two possibilities: in case there is no stiffener on the niobium wall, the interfaces programmes built the shell triangular elements model for the cavity all from the Catia's tetrahedrons volume for the vacuum space; in case the stiffeners are used, the cavity wall is modelling at first by the Catia modeller with tetrahedrons volume elements, transferred through the developed programmes to Cast3m, then the inner surface elements are shorted by the interfaces programmes to form the envelop of the vacuum space, which is modelling by Cast3m modeller for electromagnetic simulations. These models were imported to Opera 3D via the third group of the interfaces programmes.

The last group of the interfaces programmes assure the linkage between Cast3m and Opera 3D. The principal functions are, at one side, to transfer the electromagnetic distribution at the inner surface of the cavity wall from Opera to Cast3m, and at the other side, to deliver the deformed cavity geometry due to Lorentz forces from Cast3m to Opera. The basic support is the common meshing of the inner surface of the cavity wall and envelop of the vacuum space. So, if the cavity has not any stiffener, the vacuum space modelling is elaborated by Catia and imported to Opera 3D, the surface elements of his envelop as well as the electromagnetic field after Opera 3D simulations are shorted by interfaces programmes and imported to Cast3m. If the cavity is equipped the stiffeners, the cavity wall modelling is elaborated by Catia at first, the surface elements of the inner surface elements are shorted by interfaces programmes and imported to Cast3m, the Cast3m modeller make the vacuum space modelling from the inner surface elements and this modelling is exported to Opera 3D for electromagnetic simulations.

The interfaces programmes share a set of libraries in which some basic data have been stocked from Catia modelling. So the program is based on objects architecture.

MULTI-GAPS SUPERCONDUCTING SPOKE CAVITY

The 352 MHz β=0.48 multi-Spoke type cavity has been designed by the FZJ Jülich laboratory [START_REF] Zaplatin | Design of a triple-Spoke cavity for 352 MHz and beta=0[END_REF]. The design takes account at first of the electromagnetic optimizations, in order to reach the highest accelerating efficiency [START_REF] Zaplatin | Design of a triple-Spoke cavity for 352 MHz and beta=0[END_REF]. The triple-spoke is made of niobium and defined as a multi-gaps H-cavity. Some parameters of the cavity are listed in table 1. The figure 1 shows the first HIPPI triple-Spoke prototype with a stiffening system realized by FJZ Jülich. The cavity wall is made up of a cylindrical outer conductor. Three spoke tubes, rotated by 90 degrees from one to the next, are picked on the cylinder. The end cell of the cavity is conical, to minimize the magnetic field on the spoke surface.

Parameter

Figure1: Triple-Spoke cavity's design with Catia

The stiffeners on the prototype consist in five niobium ribs (10mm x 20mm) welded on the cylindrical body, two rings (18mm x 40mm) welded at the extremities of the cylinder and a niobium ring (18mm x 35mm) at each end cup.

2-1. MECHANICAL BEHAVIORS

The first mechanical issue concerns the deformation of the cavity wall under vacuum conditions. For the simulation with Cast3m, the niobium structure is modeled by tetrahedrons finite elements. The developed interface programs allow carrying out the geometrical analysis and rearrangement of the Catia meshing into a specific Cast3m meshing, figure 2.

Fig. 2: triple-Spoke cavity wall's meshing

The atmospheric pressure applies on the external wall and deforms the cavity. Two important parameters should be evaluated, the maximum von Mises stress and the resonance frequency shift due to the deformation.

Fig. 3: Von Mises stress distribution, maximum < 18 MPa

The parameters concerning the stiffeners and the results of simulations are listed in table 2 The stress distribution is shown in figure 3. The stress peak is about 18 MPa, largely lower than the elastic limit of the niobium (near 50 MPa at room temperature). The mechanical dynamics simulations have also been performed to find the structure's eigenmodes. In pulsed regime, these eigenmodes could be excited by external noise which is harmful. In general, the first mode should be controlled to be sufficiently high to avoid all risk of excitation due to environment noise. The most significant conclusion from simulations concerns the essential role of the way of fixing. The first eigenmode of the cavity is much higher in case where the cavity is fixed by external stiffening rings on the cylinder compared to the case where the cavity is fixed by his beam pipe, table 3 The simulations show that it is very important to block the cavity at the external wall of the cylinder body in order to prevent the rotation mode around the axis. The first prototype is hold by the rings on the cylinder body, in this way the cavity is largely stabilized, the risks to excite the eigenmodes are minimized.

2-2. COUPLED MECHANICAL-ELECTROMAGNETIC SIMULATIONS

Multi-disciplinary behaviours are involved concerning the Lorentz forces detuning simulations. In general, the first estimations are based on the cavity without any specific stiffening systems, but it's judicious to plan from the beginning the studies about the stiffeners design. Usually the stiffeners are added at the external surface of the cavity wall whose geometry becomes more complicated. But electromagnetic simulations are based on the vacuum space inside the inner surface of cavity wall, this surface is the only part common to the mechanical model and electromagnetic model. To limit errors in the coupling of the mechanical-electromagnetic simulations, the strategy is to use the same mesh at the inner surface of the cavity wall for both models. With this in mind, the modelling starts with the mechanical model of the bulk cavity. The interface programs have been developed with the possibility to extract the surface zones from volume elements. This way, the envelop of the vacuum space for electromagnetic simulations can be extracted from the mechanical model, then the vacuum space model is build from this envelop by Cast3m modeller.

The electromagnetic fields distribution has been simulated using Opera3d/Soprano on this vacuum space model. As the Lorentz forces depend from the electromagnetic field distribution at the inner surface of the cavity wall, only this part of the electromagnetic field distribution is transferred back to Cast3m via the developed interfaces programs. Then the radiation pressure and the mechanical deformations on the cavity wall modeling are calculated with Cast3m.

The Lorentz forces depend on the electromagnetic field on the inner surface of the cavity wall and deform the cavity wall. The radiation pressure which represents the Lorentz forces distribution acting on the cavity inner wall is expressed as the formula below:

) ( 2 0 2 0 4 1 E H P ε μ - =
From the electromagnetic fields distribution on the inner surface of the cavity wall, the radiation pressure has been calculated by Cast3m.

In the figure 4, different models of triple Spoke cavity are presented. The figure 5 shows on the left side the concentration of the radiation pressure, in middle its vector representation and at the right side the deformation of the cavity wall.

2-3. CAVITY STIFFENING

The simulations based on the design of the 1 st prototype show the efficiency of the stiffeners in reducing Lorentz forces detuning. Several other methods of stiffening have been proposed by FJZ Jülich, and discussed with IPN Orsay, for example the use of ribs at the end cup, figure 6 left side, or of copper coating at the end cup, figure 6 right side.

Fig. 6: Cavity stiffening by ribs at end-cups (left side) or by copper coating (right side)

Even if all options aren't to be realized, an interesting observation from simulations results is that all proposed systems produce quite similar results concerning the Lorentz forces reduction, the order of scale for the Lorentz forces factor's reduction using different stiffeners is compared in table 4.

Design

Δf (Hz/(MV/m)²)

without stiffener -5,3

Prototype's stiffeners The 2 gaps Spoke cavity prototype, figure 8, is built with RRR 250 3 mm thick niobium sheet. The cavity has been stiffened by semi-tube in front of each end cup. The helium vessel is made with stainless steel, it is connected to the cavity by stainless steel bellows.

Simulations are performed at IPN Orsay using the platform of calculation involved previously,

3-1. MECHANICAL BEHAVIORS

The cavity's mechanical behaviours under external 1 bar pressure have been studied. Two extreme situations have been taken into account: the case where the cavity is totally constrained and the case where the cavity is free.

The distribution of Von Mises stress at the cavity wall is shown in figure 9, the pick of the stress level is reduced thanks to the stiffener on the end cups. The maximum level of the stress doesn't exceed 34 MPa.

Fig. 8 : Von Mises stress distribution on the cavity wall: left cavity free; right: cavity fixed

The vacuum load at room temperature is performed on the prototypes, the displacements and the frequency shift have been measured. In these tests, the cavity is free, so the measured values are compared to the simulations results in the same conditions. Simulations (Nb=3mm) Measurements (Nb < 3mm) Max displacement (mm) 0.2 0.17 Max Von Mises stress (MPa) 34 NA Frequency shift (kHz) -140 -440 Table 6: mechanical parameters under vacuum The prototype has been subjected to chemical treatment, its thickness isn't uniform and somewhere the thickness is less than 2 mm, while in the simulation the cavity wall is supposed to be uniform with 3mm everywhere.

3-2. LORENTZ FORCES DETUNING SIMULATIONS

With the same procedure used for multigaps Spoke cavity, the modelling starts with the mechanical model constructed on the cavity wall. The electromagnetic fields distribution has been simulated using Opera3d/Soprano on the vacuum space modelling. The distributions of electromagnetic fields at inner surface of the cavity wall are presented in figure 9. The prototypes of these two kinds of Spoke cavities haven't the same thickness, the choice of the thickness is also suggested by the cost. But, it seems that the 3 mm 2 gaps Spoke prototype is much more sensible to the pressure on the cavity's wall than the 4 mm 4 gaps Spoke prototype, the frequency shift due to the external pressure is four times higher for 3mm 2 gaps cavity than for the 4 gaps 4 mm cavity. The same remark could be observed concerning Lorentz forces detuning factors. In parallel to this, the tuning sensitivity is four times higher for the 2 gaps prototype.

Another important remark concerns the mechanical eigenmodes. It seems if the cavity is fixed at the outside outline of the cylinder body, as shown in this report for the 4 gaps Spoke prototype, the first eigenmode is very high, it means the cavity is very stable versus external vibration excitations. In case the cavity is fixed on the beam pipe, as shown here for the 2 gaps prototype, the lower modes appear. The way of the cavity's fixation influences significantly the level of the first eigenmode.

Conclusion

IPN Orsay is involved in the studies and developments of the Spoke type superconducting cavities for HIPPI project since 2004. Important efforts have been put on numerical simulations in order to give the theoretical prediction for some fundamental parameters such as Lorentz forces detuning factor. Multi-physics analysis has been performed thanks to a platform of simulation tools interconnected by interface programs. Some simulation results have been compared with the experimental results realized on the first prototypes.

It is important to underline the fruitful collaboration between FZJ Jülich laboratory and IPN Orsay on the 4 gaps multi-Spoke 352MHz beta 0.48 cavity's studies. The results presented in this document confirm substantially the simulations made at FZJ simulating the expected performance with other simulations tools.

Fig. 4 :

 4 Fig.4: left: cavity wall model (Cast3m); right: vacuum space model (Opera3d)

  Cu on end cups and Nb ribs on body -0,9

Fig. 7 :

 7 Fig. 7: stiffened 2 gaps Spoke cavity with its helium vessel

Fig. 9 :

 9 Fig. 9: electric field (left) and magnetic field (right) distributions at inner cavity surfaceThe simulation results are listed in the table 7.

Table 1 :

 1 some parameters of triple-Spoke cavity

		Values	Units
	Frequency	352	MHz
	β=v/c	0,48
	Number of gaps	4
	R aperture	2,5	cm
	βλ/2	20,44	cm
	Rcav	21,7	cm
	G	93	Ω
	Epk/Eacc	4.65
	Bpk/Eacc	10.97 mT/MV/m

  .

	Cavity wall thickness	4 mm
	End cup stiffener	18 x 35 mm x mm
	rings	18 x 40 mm x mm
	ribs	10 x 20 mm x mm
	Maximum displacement	0.015 mm

Maximum von Mises stress 18 MPa Δf /ΔP -9.24 Hz/mbar Table 2: triple Spoke cavity under external pressure

  .

	mode	Cavity fixed by beam pipe	Cavity fixed by rings on cylinder body
	1	49 Hz	578 Hz
	2	144 Hz	669 Hz
	3	145 Hz	686 Hz
	4	160 Hz	828 Hz
	5	161 Hz	843 Hz
	6	202 Hz	887 Hz
		Table 3: triple-Spoke's eigenmodes

Table 4

 4 The two gaps Spoke cavity has been studied by IPN Orsay for the low beta section (beta 0.15). Some cavity RF parameters are listed in table 5.

	: reducing Lorentz forces detuning with different systems
	3. TWO-GAPS SUPERCONDUCTING SPOKE CAVITY
		2 gaps Spoke cavity
	Number of gaps	2
	Frequency [MHz]	352
	Beta	0.15
	Bpk/Eacc [mT/(MV/m)]	7.95 / 11.94
	Epk/Eacc	3.24
	G [Ohm]	89
	Lacc [mm]	128
	Optimal beta	0.2

Table 5 :

 5 0.15 two gaps Spoke cavity RF parameters

table 7 .

 7

		2gaps beta 0.15 Spoke prototype (3mm)	4 gaps beta 0.65 Spoke prototype (4mm) with stiffener
	Cavity stiffness	24 [kN/mm]	26 [kN/mm]
	Tuning sensitivity	964 [kHz/mm]	226 [kHz/mm]
	k L with stiffeners, fixed ends	-20 [Hz/(MV/m)²]	-1.4 [Hz/(MV/m)²]
	Pressure sensitivity, fixed ends	41 [Hz/mbar]	9.24 [Hz/mbar]
	Maximum stress	34 MPa	18 MPa
	Table 9: comparison of simulation results for two Spoke cavities

Acknowledgements

We acknowledge the support of the European Community-Research Infrastructure Activity under the FP6 "Structuring the European Research Area" programme (CARE, contract number RII3-CT-2003-506395).

The author would like to thank E. Zaplatin, R. Tölle and M. Pap of FJZ Jülich laboratory for their collaboration on the 4 gaps Spoke cavity design and fabrication, G. Olry, S. Blivet and S. Rousselot of IPN Orsay for their collaborations, also S. Chel of CEA for the coordination on WP3 of HIPPI project.

3-3. COLD TESTS MEASUREMENTS

Several series of tests were carried out on the 2 gaps Spoke prototypes at IPN Orsay in a vertical cryostat at 4.2K. The tests conditions are quite similar to the free ends conditions since the cavity has been put inside a vertical cryostat. The simulation results provide the possible range between the totally free cavity at one side and the completely constrained cavity at the other side.

After simulations, the predicted Lorentz forces detuning factor is estimated be -20 Hz/(MV/m) 2 in case where the cavity is totally constrained and to be -72 Hz/(MV/m) 2 in case where the cavity is free.

The measured Lorentz forces factors are situated between -55 Hz/(MV/m) 2 to -47 Hz/(MV/m) 2 , good agreements between simulations results and experimental data have been observed.

3-4. DYNAMIC ANALYSIS

OBSERVATIONS

Even if the 2 gaps Spoke cavity and the 4 gaps Spoke cavity are designed for the same resonance frequency, 352 MHz, but two different beta (0.15 and 0.48 respectively), the similar studies in mechanical behaviors and about Lorentz forces detuning suggest some observation on the results for the same parameter.


B Fiorina

R Vicquelin

P Auzillon

N Darabiha

O Gicquel

D Veynante

A filtered tabulated chemistry model for LES of premixed combustion

Keywords: Large Eddy Simulation, Turbulent premixed combustion, Tabulated chemistry

A new modeling strategy called F-TACLES (Filtered Tabulated Chemistry for Large Eddy Simulation) is developed to introduce tabulated chemistry methods in Large Eddy Simulation (LES) of turbulent premixed combustion. The objective is to recover the correct laminar flame propagation speed of the filtered flame front when subgrid scale turbulence vanishes as LES should tend toward Direct Numerical Simulation (DNS). The filtered flame structure is mapped using 1-D filtered laminar premixed flames. Closure of the filtered progress variable and the energy balance equations are carefully addressed in a fully compressible formulation. The methodology is first applied to 1-D filtered laminar flames, showing the ability of the model to recover the laminar flame speed and the correct chemical structure when the flame wrinkling is completely resolved. The model is then extended to turbulent combustion regimes by including subgrid scale wrinkling effects in the flame front propagation. Finally, preliminary tests of LES in a 3-D turbulent premixed flame are performed.

Introduction

Flame ignition and extinction or pollutant predictions are crucial issues in LES of premixed combustion and are strongly influenced by chemical effects. Unfortunately, despite the rapid increase in computational power, performing turbulent simulations of industrial configurations including detailed chemical mechanisms will still remain out of reach for a long time. A commonly-used approach to address fluid/chemistry interactions at a reduced computational cost consists in tabulating the chemistry as a function of a reduced set of variables. Some techniques, such as Intrinsic Low Dimensional Manifold (ILDM) developed by Mass & Pope [START_REF] Maas | Simplifying chemical kinetics: Intrinsic lowdimensional manifolds in composition space[END_REF], are based on a direct mathematical analysis of the dynamic behavior of the chemical system response. Alternative approaches are Flame Prolongation of ILDM (FPI) [START_REF] Gicquel | Laminar premixed hydrogen / air counterflow flame simulations using flame prolongation of ildm with differential diffusion[END_REF][START_REF] Fiorina | Modelling non-adiabatic partially-premixed flames using flame prolongation of ildm[END_REF] or Flamelet Generated Manifold (FGM) [START_REF] Van Oijen | Modelling of complex premixed burner systems by using flamelet-generated manifolds[END_REF]. Both techniques assume that the chemical flame structure can be described in a reduced phase subspace from elementary combustion configurations. For instance, the chemical subspace of a turbulent premixed flame can be approximated from a collection of 1-D laminar flames. In such simple situations, all thermo-chemical quantities are related to a single progress variable.

To couple tabulated chemistry with turbulent combustion, mean quantities can be estimated with presumed probability density functions. This approach, that does not require prohibitive resources, has been developed for Reynolds Averaged Navier-Stokes (RANS) computations in the past [START_REF] Vervisch | Three facets of turbulent combustion modelling: Dns of premixed flame, les of lifted nonpremixed v-flame and rans of jet-flame[END_REF][START_REF] Fiorina | Premixed turbulent combustion modelling using tabulated chemistry and pdf[END_REF]. Unfortunately, the extension of RANS turbulent combustion models to LES is not straightforward. Indeed, the primary recurrent problem is that the flame thickness is typically thinner than the LES grid size. As the progress variable source term is very stiff, the flame front cannot be directly resolved on practical LES grid meshes, leading to numerical issues. To overcome this difficulty, dedicated models have been developed under simplified chemistry assumptions. A solution to propagate a flame on a coarse grid is to artificially thicken the flame front by modifying the diffusion coefficient and pre-exponential constant [START_REF] Butler | A numerical method for two-dimensional unsteady reacting flows[END_REF][START_REF] Colin | A thickened flame model for large eddy simulations of turbulent premixed combustion[END_REF]. Following a different strategy and under simplified chemistry assumptions, Boger et al. [START_REF] Boger | Direct numerical simulation analysis of flame surface density concept for large eddy simulation of turbulent premixed combustion[END_REF] and more recently Duwig et al. [START_REF] Duwig | Study of a filtered flamelet formulation for large eddy simulation of premixed turbulent flames[END_REF] have introduced a filter larger than the mesh size to resolve the filtered flame structure. An opposite alternative is to solve a large scalar field where a given iso-surface is identified to the instantaneous flame front position. In such technique, called G-equation model, the inner layer is tracked using a level-set technique. Initially developed in a RANS context [START_REF] Peters | Turbulent Combustion[END_REF], the G-equation has been reformulated for LES [START_REF] Menon | Large eddy simulations of combustion instability in an axisymetric ramjet combustor[END_REF][START_REF] Chakravarthy | Subgrid modeling of turbulent premixed flames in the flamelet regime[END_REF][START_REF] Pitsch | A consistent level set formulation for large-eddy simulation of premixed turbulent combustion[END_REF]. However as level-set techniques provide information only on the thin reaction zone position and not on the filtered flame structure, the coupling with the flow equations is challenging. In particular the knowledge of the temperature field is required for taking into account heat expansion. As recently proposed by Moureau et al. under simplified chemistry assumption [START_REF] Moureau | A level set formulation for premixed combustion les considering the turbulent flame structure, combust[END_REF], a solution is to solve an additional progress variable equation to ensure a consistent coupling with a LES flow solver.

The FPI-PCM (Presumed Conditional Moment) model [START_REF] Galpin | Large-eddy simulation of a fuel-lean premixed turbulent swirl-burner[END_REF], developed to introduce tabulated chemistry effects in LES, combines presumed Probability Density Functions (PDF) and FPI tables to describe the chemical reaction rate of the filtered progress variable accounting for interactions between turbulence and chemistry at the subgrid scale level. However, as will be shown further, this formulation does not guarantee a proper prediction of regimes where the subgrid scale flame wrinkling vanishes. This regime, observed when the subgrid fluctuations are lower than the laminar flame speed, is encountered in practical LES of premixed combustion [START_REF] Pitsch | Large-eddy simulation of turbulent combustion[END_REF][START_REF] Moureau | A level set formulation for premixed combustion les considering the turbulent flame structure, combust[END_REF]. Additionally LES should tend toward DNS when the filter size becomes lower than the Kolmogorov scale. Hawkes & Cant [START_REF] Hawkes | Physical and numerical realizability requirements for flame surface denisty approaches and reynolds averaged simulation of premixed turbulent combustion[END_REF] extensively discussed realizability in premixed combustion LES.

In the present work, it is first demonstrated that the β-PDF formalism applied in the context of premixed combustion LES does not guarantee a proper description of a filtered laminar flame front. Therefore an alternative is proposed to include tabulated chemistry in LES approach ensuring the correct propagation speed of the filtered laminar flame front. The resolved flame structure is mapped from 1-D filtered laminar premixed flames. The idea of tabulating filtered quantities has already been introduced [START_REF] Vreman | Subgrid scale modeling in large eddy simulation of turbulent combustion using premixed flamelet chemistry[END_REF] but unresolved convective and diffusive terms where neglected. As it will be demonstrated further, these assumptions do not allow a proper description of the filtered flame structure and propagation. Here, closure of filtered flow and progress variable equations are first carefully addressed in regimes where the flame wrinkling is fully resolved. One-dimensional computations are performed to investigate the capability of the proposed model to reproduce the correct propagation speed and the filtered flame structure. The model is then extended to turbulent combustion regimes taking into account subgrid scale flame wrinkling. Finally, simulations of a turbulent swirled premixed flame are performed and compared to experimental data.

Coupling tabulated chemistry and LES: filtered equations

Low-dimensional trajectories in composition space are identified in FPI framework from the knowledge of the complex chemical structure of 1-D laminar flames [START_REF] Gicquel | Laminar premixed hydrogen / air counterflow flame simulations using flame prolongation of ildm with differential diffusion[END_REF]. For premixed combustion systems, a 1-D freely propagating flame is first computed using detailed chemical schemes. Thermodynamical and chemical quantities are then tabulated as a function of a unique monotonic progress variable c related to temperature or to a combination of chemical species, where c = 0 corresponds to fresh gases and c = 1 to fully burnt gases. The chemical database is then coupled to the flow field by adding the progress variable balance equation to the Navier-Stokes equations. The progress variable reaction rate and heat release are extracted from the chemical database. For LES, under unity Lewis numbers assumption, these equations are filtered leading to the following system :

∂ ρ ∂t + ∇ • (ρ u) = 0 (1)
∂ ρ u ∂t + ∇ • (ρ u u) = -∇P + ∇ • τ -∇ • (ρ uu -ρ u u) (2)
∂ ρ c ∂t + ∇ • (ρ u c) = ∇ • ρD∇c -∇ • (ρ uc -ρ u c) + ρ ωc (3)
∂ ρ E ∂t + ∇ • (ρ u E) = -∇ • P uδ + ∇ • (τ u) -∇ • ρ uE -ρ u E +∇ • ρD∇h s + ρ ωE (4)
P = ρ r T ( 5
)
where ρ is the density, u the velocity vector, P the pressure, δ the unit tensor, τ the laminar viscous tensor, E = H -P/ρ with H the total nonchemical enthalpy, h s the sensible enthalpy, D is the diffusivity, ωc and ωE , respectively, the progress variable and energy source terms. r = R/W where R is the ideal gas constant and W the mean molecular weight. The overbar denotes the spatial filtering operation,

φ(x) = F (x -x ′ )φ(x ′ )dx ′ , (6)
where φ represents reactive flow variables and velocity components and F the filtering function. The tilde operator denotes the density-weighted filtering defined by ρ φ = ρφ.

The subgrid scale terms, -∇ • (ρ uu -ρ u u) and -∇ • (ρ uϕ -ρ u ϕ), where ϕ denotes c or E quantities, the pressure term P u, as well as the filtered laminar diffusion terms ρD∇ϕ and the filtered source terms ωϕ , require closure models. The model constraints are both to ensure a correct flame propagation and to recover the chemical structure of the filtered flame under two possible situations : (1) the flame wrinkling is fully resolved at the LES filter size, and (2) wrinkling occurs at the subgrid scale and affects the filtered flame speed.

Different strategies exist to model the filtered progress variable reaction rate ωc . An approach that does not require extensive CPU resources is to presume the shape of progress variable PDF, generally by a β function. This formalism has been applied to LES of turbulent premixed combustion [START_REF] Galpin | Large-eddy simulation of a fuel-lean premixed turbulent swirl-burner[END_REF] but, to our knowledge, the ability of the method to reproduce the propagation speed of filtered flame front has not yet been investigated. In the following section the influence of the PDF shape on the filtered flame properties is discussed when the flame wrinkling is resolved at the LES filter scale i.e. when the subgrid scale flame front remains laminar and planar. The use of a β function is found to introduce errors in the filtered flame front propagation speed. A new modeling alternative based on the tabulation of filtered premixed flame elements is then proposed to correct this drawback.

A priori testing of presumed β-PDF formalism in the laminar regime

An unstretched 1-D filtered laminar premixed flame is considered in this section. If no wrinkling occurs at the subgrid scale, the propagation speed S ∆ of the filtered flame front is identical to the laminar flame speed S 0 l . The following relation then needs to be satisfied:

ρ 0 S ∆ = +∞ -∞ ρ ωc (x)dx = +∞ -∞ ρ ωc (x)dx = ρ 0 S 0 l (7)
where ρ 0 is the fresh gases density and x is the spatial dimension.

The ability of presumed β-PDF to satisfy this property is investigated by conducting a priori tests on a 1-D stoichiometric freely propagating laminar premixed propane/air flame computed with PREMIX [START_REF] Kee | A fortran program for modelling steady laminar one-dimensional premixed flames[END_REF] using a modified version of the GRI 3.0 mechanism [21]. The progress variable c is plotted as a function of the spatial coordinate x in Fig. 1(a). The laminar flame thickness, defined by δ l = 1/ max(|dc/dx|) is approximately equal to 0.4 mm. Introducing P , the mass weighted PDF defined by ρ P = ρP , the progress variable filtered reaction rate reads:

ωc (x) = 1 0 ωc (c) P (x, c)dc (8)
Assuming that c follows a β distribution [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]:

P (x, c) = c ac-1 (1 -c) bc-1 1 0 c ac-1 (1 -c) bc-1 dc (9)
where parameters a c and b c are determined from c and the segregation factor

S c = ( cc -c c) / ( c(1 -c)): a c = c 1 S c -1 ; b c = a c 1 c -1 (10)
The knowledge of the first and second moment of the progress variable provides the filtered reaction rate ωc = ωc ( c, S c ). For the configuration considered here, c and S c profiles across the filtered laminar flame front are computed by applying a 1-D Gaussian filter F of size ∆ defined by:

F (x) = 6 π∆ 2 1/2 e -6x 2 ∆ 2 (11)
on the detailed chemistry laminar flame solution.

Favre-filtered progress variable and the segregation factor are shown in Fig. 1(a) for a filter size of ∆ = 20δ l . According to Eq. 9, the presumed β-PDF , P (x, c), is deduced from these two quantities. The reaction rate ωc across the filtered flame front is then estimated from Eq. 8. Finally, the integration of the filtered reaction rate according to Eq. 7 gives an a priori estimation of the filtered flame front propagation speed S ∆ . The ratio S ∆ /S 0 l (square symbols) is plotted as a function of the ratio ∆/δ l in Fig. 1(b). When ∆/δ l < 1 the effect of the β-PDF on the flame structure is moderate and the propagation speed is correctly reproduced. However when the filter size is larger than the flame front, as in LES practical situations, the propagation speed of the filtered progress variable is largely over-estimated by the presumed β function up to a factor of 2.5. In fact, the β-PDF is not relevant when subgrid scale wrinkling is resolved.

A solution to propagate a flame front at the correct speed is to artificially thicken the reaction zone. In the Thickened Flame model for LES (TFLES) [START_REF] Butler | A numerical method for two-dimensional unsteady reacting flows[END_REF][START_REF] Colin | A thickened flame model for large eddy simulations of turbulent premixed combustion[END_REF], both reaction rate and diffusion fluxes are affected in order to ensure a correct propagation of the flame front. However the structure of the thickened flame front does not correspond to the filtered flame front.

An alternative to presumed PDF formalism and TFLES is to directly employ a normalized filter function F (x) to estimate the filtered reaction rate. Then the filtered reaction rate reads:

ωc (x) = 1 ρ +∞ -∞ ρ(x ′ ) ωc (x ′ )F (x -x ′ )dx ′ , (12)
Since by definition, F (x) satisfies +∞ -∞ F (x)dx = 1, Eq. 7 is then always satisfied:

ρ 0 S ∆ = +∞ -∞ ρ ωc (x)dx (13) = +∞ -∞ +∞ -∞ ρ(x ′ ) ωc (x ′ )F (x -x ′ )dx ′ dx (14) = +∞ -∞ ρ(x ′ ) ωc (x ′ ) +∞ -∞ F (x -x ′ )dx dx ′ (15) = +∞ -∞ ρ(x ′ ) ωc (x ′ )dx ′ (16) = ρ 0 S 0 l ( 17
)
This property is verified in Fig. 1(b) where the propagation speed S ∆ of the filtered flame, is a priori computed using Eqs. 11 and 12.

By taking advantages of this property, a model is proposed in Section 4 to ensure the correct propagation of filtered laminar flame front. The closure of unknown terms is carefully addressed and the model is tested on 1-D filtered flame configurations. This approach is extended to turbulent regimes where subgrid flame wrinkling occurs at the subgrid scale level in Section 5 by the introduction of the subgrid flame wrinkling factor.

Filtered laminar premixed flames modeling

Modeling

The flame structure in the direction n normal to the flame front is assumed identical to the structure of a planar 1-D freely propagating premixed laminar flame. A detailed chemical mechanism with N s species is considered. From this reference flame structure and using the filter operators introduced in Section 2, the a priori filtered flame structure is determined. For instance, for a given filter size ∆, any filtered fluxes or filtered thermo-chemical quantities of a planar filtered laminar flames are known.

As an example, a 1-D laminar stoichiometric premixed propane/air flame is computed taking into account detailed chemistry effects. The PREMIX [START_REF] Kee | A fortran program for modelling steady laminar one-dimensional premixed flames[END_REF] solver is combined with a modified version of the GRI 3.0 mechanism [21] involving N s = 70 species and 463 elementary reactions. The filtered operator given by Eq. ( 11) is then applied to the 1-D laminar flame solution. Figure 2 shows all budget terms of the c balance equation in a steady 1-D laminar premixed flame remapped in the c coordinate system for different values of ∆. Using these results, a modeling procedure based on the tabulation of the filtered 1-D laminar flame structure is proposed in the following sections. The closure of each unclosed terms identified in Eqs. 2 to 4 is first carefully addressed in the situation where no flame wrinkling occurs at the subgrid scale level.

Filtered chemical reaction rates

The filtered source terms for c and the energy equations are directly given by the filtered database:

ωϕ = ω * ϕ [ c, ∆]. ( 18
)
where ϕ denotes c or E quantities and the * superscript denotes quantities issued from a 1-D unstretched laminar premixed flame. The notation φ[ c, ∆] means that the variable φ is tabulated in a 2-D look-up table with coordinates c and ∆. Figure 2 shows that the filter operator affects dramatically both the amplitude and the shape of ωc (triangles symbols) profiles.

Filtered laminar diffusion terms ∇ • (ρD∇c) and ∇ • (ρD∇h) These terms are usually neglected or approximated as [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]:

∇ • ρD∇ϕ ≈ ∇ • (ρD ∇ ϕ) . ( 19
)
This approximation is very rough and may introduce large errors. Indeed, in Fig. 2

∇ • (ρD ∇c) = -∇ • (ρD |∇c|n) (20) = -∇ • ρD ∂c * ∂x * n (21) (22)
By introducing a corrective factor α c ( c), one can write:

∇ • (ρD ∇c) = ∇ • (α c [ c, ∆] ρD ∇ c) . (23)
The normal to the flame front n = -∇ c/|∇ c| points into the fresh reactants.

The correction factor α c ( c) is defined as:

α c [ c, ∆] = ρD ∂c * ∂x * ρD ∂ c * ∂x * . ( 24
)
The quantity α c [ c, ∆] is estimated from the 1-D filtered flame solution and is tabulated as a function of c for a given value of filter size ∆.

Similarly, the energy-filtered laminar diffusion term is written as:

∇ • (ρD ∇h s ) = ∇ • α E ([ c, ∆]) ρD ∇ h s , (25)
where the correction factor α E [ c, ∆] is defined as:

α E [ c, ∆] = ρD ∂h * s ∂x * ρD ∂ h * s ∂x * . ( 26
)
The correction factors α c [ c, ∆] and α E [ c, ∆] are plotted in Fig. 3 for different values of filter size ∆. For small values of ∆, as α c [ c, ∆] remains constant and close to 1, effects on the laminar diffusion fluxes modeling will be negligible. However, the profiles present strong variations in terms of c when the filter size ∆ is larger than δ l .

Unresolved convection terms

-∇ • (ρ uϕ -ρ u ϕ)
The displacement speed s d , measuring the flame front local speed relative to the flow, i.e. the difference between the absolute flow velocity u and the absolute flame front speed w, is first introduced:

u = w + s d (27)
The filtered flame front speed w remains constant across the flame brush ( w = w = w), therefore after replacing the flow velocity by relation 27, the subgrid scale convection term then reads: In a 1-D laminar premixed flame the laminar flame speed S 0 l and the fresh gas mixture density ρ 0 are related to the displacement speed through the following relation:

-∇ • (ρ uϕ -ρ u ϕ) = -∇ • (ρ s d ϕ -ρ s d ϕ) (28)
ρ 0 S 0 l = ρs * d ( 29
)
Therefore, under the assumption that the flame remains planar at the subgrid scale level, the unresolved convection terms are directly estimated from the reference laminar 1-D detailed chemistry premixed flame:

-∇ • (ρ uϕ -ρ u ϕ) = - ∂ ∂x * ρ s * d ϕ * -ρ s * d ϕ * (30) = -ρ 0 S 0 l ∂ϕ * ∂x * - ∂ ϕ * ∂x * . ( 31
) = Ω ϕ [ c, ∆] (32)
The 2 for different values of filter size ∆ (squares). For ∆ < δ l , unresolved convective fluxes are very small compared to other fluxes. However, when ∆ ≥ δ l , these fluxes become important and are counter-gradient type. Note that this result is in agreement with recent experiments [START_REF] Pfadler | Direct evaluation of the subgrid scale scalar flux in turbulent premixed flames with conditioned dual-plane stereo piv[END_REF]. The quantity Ω ϕ [ c, ∆], estimated from the 1-D filtered flame solution, is then tabulated as a function of c and ∆. In practice, as the unresolved convective terms are modeled as a source term, only the sum Σ ϕ [ c, ∆] = Ω ϕ [ c, ∆] + ωϕ [ c, ∆] is stored in the filtered database where φ denotes c or E quantities.

term Ω c [ c, ∆] = -ρ 0 S 0 L ∂ ∂x * (c -c) is plotted in Fig.

Pressure term

In a similar way, the pressure term in the energy equation (Eq. 4) is written as:

-∇ • (P uδ) = -∇ • (P uδ) -∇ • (P uδ) -∇ • (P uδ) (33) = -∇ • (P uδ) -∇ • (ρ rT uδ) -∇ • (ρ rT uδ) (34) = -∇ • (P u δ) + Ω p [ c, ∆] (35)
with

Ω p [ c, ∆] = -ρ 0 S 0 l ∂(rT * ) ∂x * - ∂( rT * ) ∂x * . ( 36
)

Momentum equations

Unclosed terms in the filtered momentum equations may be modeled following the same approach. The subgrid scale convection term is written as:

-∇ • (ρ uu -ρ u u) = ∂ ∂x * ρ s * d s * d -ρ s d * s d * n (37) = ρ 0 S 0 l ∂s d * ∂x * - ∂ s d * ∂x * n (38) = Ω u [ c, ∆]n (39)
The strain tensor is expressed by:

∇ • τ = ∇ • (α u [ c, ∆] τ ) with α u [ c, ∆] = τ * τ * . ( 40
)
where τ is defined as:

τ = µ ∇ u + (∇ u) T - 2 3 (∇ • u)δ (41)
However, as shown further, the influence of these terms is moderate and can be neglected.

Summary of the model equations

The momentum, the progress variable and the energy equations are modeled as:

∂ ρ u ∂t + ∇ • (ρ u u) = -∇P + ∇ • (α u [ c, ∆] τ ) + Ω u ( c)n (42) ∂ ρ c ∂t + ∇ • (ρ u c) = ∇ • (α c [ c, ∆] ρD ∇ c) + Σ c [ c, ∆] (43)
∂ ρ E ∂t + ∇ • (ρ u E) = -∇ • P u δ + Ω p [ c, ∆] + ∇ • ( τ u) +∇ • α E [ c, ∆] ρD ∇ h s + Σ E [ c, ∆] (44)
These equations are implemented in the compressible LES code AVBP [25]. The third-order finite element scheme TTGC [START_REF] Colin | Development of high-order taylor-galerkin schemes for unsteady calculations[END_REF] is used. Boundary conditions are prescribed using Navier-Stokes Characteristic Boundary Conditions [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF].

The sum of filtered chemical reactions rates and the subgrid scales fluxes Σ ϕ = Ω ϕ + ρ ωϕ and the diffusion fluxes correction factors α ϕ are estimated after filtering a 1-D laminar stoichiometric premixed propane/air flame. These quantities are stored in a look-up table as a function of c and ∆.

1-D laminar premixed flame simulations

Filtered steady 1-D laminar flames are computed to verify the ability of the present model to reproduce both the correct flame front propagation speed and the filtered flame structure. Computations are performed on uniform meshes with a grid spacing of ∆ x . A parametric study is conducted for different filter sizes relative to the laminar flame thickness. For each case, a reference solution is obtained by filtering the 1-D laminar premixed flame detailed chemistry solution. The simulations are initialized with the reference solution and the overall physical time for each run is t run = 50 δ e c /S 0 l , where δ e c = 1/ max(| ∂e c ∂x |) is an estimation of the filtered flame thickness. A comparison between the numerical solutions on uniform mesh (solid lines) and the reference solution (dashed line) with δ e c /∆ x = 50 and for different values of ∆/δ l is first shown in Fig. 4. The predicted filtered progress 19), i.e., α ϕ = 1. This rough assumption leads to an under-prediction by a factor of 3 of the flame front propagation speed.

An important information for premixed combustion LES is the minimal number of grid points required to capture the filtered flame front without introducing numerical artifacts. The filtered flame front propagation speed is plotted as a function of the mesh resolution ∆ x in Fig. 5(b). The flame speed is recovered with a good approximation for δ e c /∆ x ≥ 5. Below this limit, numerical errors become important and the filtered flame front does not propagate at the correct speed. Then, for numerical reasons, the filter should be at least 5 times larger than the mesh size. Note that even approaches based on level-set transport that use sophisticated numerical methods to track the flame front position also require to filter the flame front at a scale larger than the mesh size in order to resolve density gradients [START_REF] Moureau | A level set formulation for premixed combustion les considering the turbulent flame structure, combust[END_REF]. Finally, a simulation has been performed without considering the filtering effect on the momentum equations (Eq. 42), i.e., with α u = 1 and Ω u = 0 and is compared with the complete model solution in Fig. 6. For both simulations, density as well as velocity profiles match perfectly. In fact, the induced differences are transfered to the pressure that becomes a macro-pressure. As this macro-pressure remains very close to the static pressure, effects on the thermodynamic state are very limited. Then, in order to simplify the model implementation in 3-D configurations, the contribution corresponding to the filtering of a laminar flame in the momentum equation will be neglected.

Filtered turbulent premixed flames modeling

In practical LES of turbulent combustion, turbulence may cause flame front wrinkling at the subgrid scale level. Here, a strategy is proposed to extend the previously described model to such situations.

Modeling

Turbulent structures induce flame wrinkling that increases the flame surface area at the subgrid scale. As a consequence the filtered flame front propagates at a subgrid scale turbulent flame speed S t [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF] related to the laminar flame speed through the flame wrinkling factor Ξ = S t /S 0 l . The model developed here ensures that the filtered flame front propagates at the turbulent flame speed S t . The filtered flame thickness is assumed to be only related to the filter size ∆ and is not altered by small-scale eddies.

Then, the filtered progress variable turbulent reaction rate is modeled by:

ωct = Ξ . ω * c [ c, ∆] (45)
and the turbulent diffusion term is expressed as follows:

Ω ct = -(∇ • (ρ uc -ρ u c)) t = Ξ Ω c [ c, ∆] + (Ξ -1) ∇ • (α c [ c, ∆] ρD ∇ c) (46)
The first term on the r.h.s corresponds to the thermal expansion and the second one models the unresolved turbulent fluxes. This formulation corresponds to multiply diffusion and source terms by the flame wrinkling factor in the laminar flame balance equation and then ensures that the unstretched filtered flame front propagates at the turbulent flame speed S t = ΞS 0 l in the normal direction.

Summary of the model equations

To summarize, momentum, progress variable and energy equations for this new model called Filtered Tabulated Chemistry for LES (F-TACLES) can be written as follows:

∂ ρ u ∂t + ∇ • (ρ u u) = -∇P + ∇ • τ + ∇ • τ t (47) ∂ ρ c ∂t + ∇ • (ρ u c) = Ξ∇ • (α c [ c, ∆] ρD ∇ c) + ΞΣ c [ c, ∆] (48)
∂ ρ E ∂t + ∇ • (ρ u E) = -∇ • (P u δ) + ΞΩ p [ c, ∆] + ∇ • ( τ u) +Ξ∇ • α E [ c, ∆] ρD ∇ h s + ΞΣ E [ c, ∆]. ( 49
)
Note that here the effect of the flame filter ∆ on the momentum equations is neglected and the subgrid scale turbulent fluxes ∇ • τ t are modeled using the Smagorinsky model. Different alternatives exist to estimate the subgrid flame wrinkling factor that appears in Eqs 48 and 49. It can be either estimated from analytical models [START_REF] Colin | A thickened flame model for large eddy simulations of turbulent premixed combustion[END_REF][START_REF] Pitsch | A consistent level set formulation for large-eddy simulation of premixed turbulent combustion[END_REF][START_REF] Charlette | A power-law flame wrinkling model for les of premixed turbulent combustion, part i: non-dynamic formulation[END_REF][START_REF] Charlette | A power-law flame wrinkling model for les of premixed turbulent combustion, part ii: dynamic formulation[END_REF] or from the solution of a surface density balance equation [START_REF] Hawkes | A flame surface density approach to largeeddy simulation of premixed turbulent combustion[END_REF][START_REF] Richard | Towards large eddy simulation of combustion in spark ignition engines[END_REF].

Large Eddy Simulation of a swirled premixed burner

The proposed method is applied to the simulation of the complex PREC-CINSTA swirled burner experimentally investigated by Meier et al. [START_REF] Meier | Detailed characterization of the dynamics of thermoacoustic pulsations in a lean premixed swirl flame[END_REF]. The geometry, shown in Fig. 7, derives from an aeronautical combustion device. It features a plenum, a swirl-injector and a combustion chamber. Details of the burner geometry and of the measurement can be found in Ref. [START_REF] Meier | Detailed characterization of the dynamics of thermoacoustic pulsations in a lean premixed swirl flame[END_REF]. Different modeling strategies for LES have been used to numerically investigate this configuration : an LES of the combustor using the thickened flame model and a two-step mechanism has been first performed by Roux et al.. [START_REF] Roux | Studies of mean and unsteady flow in a swirled combustor using experiments, acoustic analysis, and large eddy simulations[END_REF]. Moureau et al. [START_REF] Moureau | A ghost-fluid method for large-eddy simulations of premixed combustion in complex geometries[END_REF] used this configuration to validate a new level-set algorithm to track the flame front position. Recently, Galpin et al. [START_REF] Galpin | Large-eddy simulation of a fuel-lean premixed turbulent swirl-burner[END_REF] performed the LES of this lean premixed burner by using a presumed β-PDF to couple a thermo-chemical look-up table with the filtered flow equations.

The operating conditions chosen in the present study correspond to an air mass flow rate of 12.2 g/s and to a methane mass flow rate of 0.6 g/s. In the experiment, air and methane are injected separately in the plenum inlet, however in the present simulation the mixing is assumed to be fast enough to burn a perfect mixing of oxidizer and fuel in the combustion chamber. Thus methane injection is not taken into account and a methane/air mixture characterized by an equivalence ratio of 0.83 is injected at the plenum inlet. These conditions correspond to a stable regime where laser Raman scattering has been performed, allowing comparison between predicted and measured thermo-chemical quantities such as temperature and species mass fractions.

The boundary conditions and the computational geometry have been already described in [START_REF] Roux | Studies of mean and unsteady flow in a swirled combustor using experiments, acoustic analysis, and large eddy simulations[END_REF]. The mesh used to perform the computation is unstructured and made of 12.7 millions elements. The third-order finite element scheme TTGC [START_REF] Colin | Development of high-order taylor-galerkin schemes for unsteady calculations[END_REF] is retained. For building-up the chemical look-up table, a 1-D laminar methane/air flame is first computed for an equivalence ratio equal to 0.83 using the GRI 3.0 mechanism [21]. Then, according to the modeling procedure discussed previously, this laminar flame solution is filtered by the Gaussian function defined by Eq. 11.

Note that, as the mesh considered here is almost uniform in the filtered flame front region, an unique filter width ∆ is considered. In order to ensure a sufficient meshing of the filtered flame front, the filter width has been set to ∆ = 20δ l . The progress variable is defined by c = Y CO 2 /Y eq CO 2 , where Y eq For strongly non-uniform meshes this procedure is not optimized and could lead to overrefined or under-refined flame front regions. Then, an additional coordinate, the filter width, can be easily considered when computing the look-up table.

Following the system of equations 47-49, this new model F-TACLES has been implemented into the compressible LES code AVBP [25]. The subgrid flame wrinkling factor Ξ is estimated from the analytical model developed by Colin et al. [START_REF] Colin | A thickened flame model for large eddy simulations of turbulent premixed combustion[END_REF]. Mean and resolved Root Mean Square (RMS) quantities are computed by time averaging LES solutions over a physical time that correspond to 6 flow-through times based on the fresh gas inlet velocity. Mean temperature and CO 2 mass fractions are plotted on Figs. 8 (top) and 9 (top), respectively. A very good agreement is observed between experimental and numerical profiles, which demonstrated that the correct flame angle and mean flame thickness are reproduced by the model. Because heat losses have not been considered when generating the chemical database and in the numerical simulation, the LES slightly over-estimates the temperature profiles close to the combustion chamber wall, in the outer recirculation zone for x < 20 mm and at a distance larger than 20 mm from the jet axis. Note that heat losses effects on the flame structure can be taken into account with the addition of the enthalpy as a control parameter of the chemistry tabulation [START_REF] Fiorina | Modelling non-adiabatic partially-premixed flames using flame prolongation of ildm[END_REF][START_REF] Fiorina | Premixed turbulent combustion modelling using tabulated chemistry and pdf[END_REF].

Figs. 8 (bottom) and 9 (bottom) show a comparison between resolved LES RMS and measured RMS of the temperature and the CO 2 mass fraction, respectively. As the plotted LES RMS does not include the subgrid scale RMS, conclusions regarding the model performance in terms of flame turbulence interactions are more difficult. However, it is observed that LES RMS remains lower than measured RMS, as expected from theory.

As all thermo-chemical variables are related to c, the post-processing of the filtered progress variable solution with the filtered chemical database allows to access all chemical species. As an example, Fig. 10(a) shows 2-D contours of c used to estimate HCO mass fraction plotted in Fig. 10(b).

Finally, Fig. 11 indicates the flame position in the Pitsch LES regime diagram for turbulent premixed combustion [START_REF] Pitsch | Large-eddy simulation of turbulent combustion[END_REF], where the ratio ∆/δ l is expressed as a function of the Karlovitz number Ka in logarithmic scale. The Karlovitz number is related in LES to the subgrid velocity fluctuations v ′ ∆ and laminar flame scales [START_REF] Pitsch | Large-eddy simulation of turbulent combustion[END_REF]:

Ka 2 = δ l S 0 l 3 ε = v ′ ∆ S 0 l δ l ∆ ( 50
)
where ε is the kinetic energy transfer rate. The subgrid velocity fluctuations are computed as follows:

v ′ ∆ = µ t ρC k ∆ 3/2 (51)
where the turbulent viscosity µ t is estimated from Smagorinsky model. For Ka < 1, combustion takes place in the corrugated flame regime while the thin reaction zone regime is observed when Ka > 1. Computational nodes located in the filtered flame front are considered, i.e. for 0.01 < c < 0.99, and are plotted in the LES diagram (horizontal thick solid black line in Fig. 11). As a a unique filter width ∆ is considered in the present simulation, the scatter plot reduced to the line ∆/δ l = 20. The smallest size of the flame wrinkling is given by the Gibson length [START_REF] Peters | Turbulent Combustion[END_REF]:

∆ l G = v ′ ∆ S 0 l (52)
The substitution of Eq. 52 into Eq. 50 shows that ∆ = l G condition corresponds to ∆/δ l = Ka -2 represented by a line of slope -2 in the LES diagram (Fig. 11). In the corrugated flame regime, when the filter width becomes smaller than the Gibson length, the subgrid velocity fluctuation v ′ ∆ is smaller than the laminar flame speed S 0 l . In such cases, the flame wrinkling is fully resolved at the LES filter scale. At the opposite, on the right side of the l G = ∆ line, subgrid scale wrinkling exists and will impact the filtered flame front propagation speed S ∆ . The node distribution versus the Karlovitz number is plotted in Fig. 12. First, it can be observed that most of the points are located in the corrugated flame regime (Ka < 1). The chemical flame structure remains therefore laminar as assumed in the present model. Secondly, for a substantial area of the flame surface ( about 30 %), the Gibson length l G is larger than the filter width and consequently the flame wrinkling is fully resolved. With future increase of computational power, as meshes will be finer, this trend should be emphasized. It demonstrates the crucial need of ensuring a proper propagation of the laminar flame front when deriving a turbulent combustion model.

Conclusion

A new modeling strategy called Filtered Tabulated Chemistry for LES (F-TACLES) has been developed to introduce tabulated chemistry methods in premixed combustion LES. A filtered 1-D laminar premixed flame is used to build a filtered chemical look-up table. The model performances are demonstrated on 1-D filtered laminar flame computations. Finally the proposed strategy has been applied to perform a 3-D simulation of a swirled turbulent premixed flame. Good agreement between the numerical simulation and the experiments is observed.

Figure 1 :

 1 Figure 1: A priori test of the β-PDF formalism in laminar regime. Left (a): progress variable c (solid line) and filtered progress variable c (bold line) profiles as a function of the spatial coordinate x. Dashed line is the subfilter progress variable segregation factor S c . Right (b): a priori computations of the filtered progress variable propagation speed for different values of filter size. The filtered progress variable reaction rate is modeled by a β-PDF (squares) or by a Gaussian filter (triangles).

Figure 2 :

 2 Figure 2: Budget terms (in kg.m -3 .s -1 ) as a function of c of the filtered progress variable balance equation of a steady 1-D laminar planar filtered premixed flame for different values of filter size ∆ : ∂ ρe u * e c * ∂x * = ∂ ∂x * ρD ∂c * ∂x * -∂ ∂x * ρ u * c * -ρ u * c * + ρ ω * c . -: ∂ ρe u * e c * ∂x * . :

  the exact laminar diffusion fluxes ∂ ∂x * ρD ∂c * ∂x * (filled diamonds) and the approximation by ∂ ∂x * ρD ∂e c * ∂x * (empty diamonds) are shown for different values of the filter size ∆. When the filter size is smaller than the laminar flame thickness δ l , approximation by Eq. 19 remains valid. However as soon as the filter size ∆ becomes larger than δ l , important differences are observed between ∂ ∂x * ρD ∂c * ∂x * and ∂ ∂x * ρD ∂e c * ∂x * . As shown further, these errors impact dramatically the prediction of the propagation speed. In the present work, the filtered diffusion term for the c equation is modeled by:

Figure 3 :

 3 Figure 3: Diffusion correction factor α c (left) and α E (right) as a function of c for different values of ∆. Dashed dotted dotted lines: ∆ = 0.2δ l . Dashed lines: ∆ = δ l . Dashed dotted lines: ∆ = 5δ l . Solid lines: ∆ = 25δ l .

Figure 4 :

 4 Figure 4: Filtered 1-D premixed flame solutions. Filtered progress variable (solid) compared to the reference solution (dashed) for ∆/δ l = 2, 10 and 20.

  (a) δ e c /∆ x = 50 (b) ∆/δ l = 20

Figure 5 :

 5 Figure 5: Predicted flame speed as a function of ∆/δ l (left) and δ e c /∆ x (right). Square symbols are the complete model solution and the triangle symbol is the solution with α ϕ = 1.

Figure 6 :

 6 Figure 6: Filtered 1-D premixed flame solutions. Effects of the flame filter in the momentum equation. Solid: α u = 1 and Ω u = 0. Symbols: α u ( c) and Ω u ( c) from the filtered database.

CO 2 is

 2 the equilibrium CO 2 mass fraction in the fully burnt gases. The filtered quantities required by the model: Σ c [ c, ∆], α c [ c, ∆], Ω p [ c, ∆], Σ E [ c, ∆] and α E [ c, ∆] are then tabulated as a function of c for ∆ = 20δ l .

Figure 7 :

 7 Figure 7: LES of Preccinsta with F-TACLES turbulent combustion model. The computational domain features the plenum, the swirl-injector and the combustion chamber. An instantaneous view of the filtered flame front iso-surface ( c=0.8) is shown.

Figure 8 :Figure 9 :

 89 Figure 8: Mean (top) and RMS (bottom) of temperature, case φ = 0.83. Symbols: measurements. Lines: simulation with F-TACLES. x = 0 matches the swirler exit.

Figure 10 : 2 -

 102 Figure 10: 2-D instantaneous view of c and Y HCO .

Figure 11 :

 11 Figure 11: LES regime diagram for turbulent premixed combustion. The thick solid black line represent the range covered by the Preccinsta flame simulation.

Figure 12 :

 12 Figure 12: Node distribution versus the Karlovitz number. Only nodes located into the filtered flame front have been considered, i.e. for 0.01 < c < 0.99.

Acknowledgments

This work was supported by the ANR-07-CIS7-008-04 Grant of the French Ministry of Research. We are grateful to the CERFACS (Toulouse, France) combustion team for providing the PRECCINSTA burner geometry. The authors warmly acknowledge the support of the 2008 Summer Program of the Center for Turbulence Research (Stanford University -NASA Ames) during which this work was initiated. This work was granted access to the HPC resources of IDRIS under the allocation 2009-i2009020164 made by GENCI (Grand Equipement National de Calcul Intensif)


Olivier Esnault

Matthieu Boileau
email: matthieu.boileau@em2c.ecp.fr

Ronan Vicquelin

Benoît Fiorina

Olivier Gicquel

A Method To Accelerate LES Explicit Solvers Using Local Time-Stepping

In practical flow configurations, a large disparities of geometrical length scales are often encountered. Inside a combustor, for example, the ratio between the diameter of the injection holes and the size of the entire combustion chamber may present several orders of magnitude. When considering an explicit solver for fully compressible Navier-Stokes equations, the global time step is constrained through a CFL-like condition by the size of the smallest cells in the overall computational domain. Local refinement of the injector leads to an inhomogeneous mesh and the former restriction drastically alters the overall solver efficiency. A new local time-stepping (LTS) method is proposed to address this issue. The domain is divided into subgrids composed of cells that have similar sizes. Flow equations are simultaneously advanced on each subgrid which have a local time step adapted to satisfy the local CFL condition. The accuracy of the method has been verified on a simple convection case using a test code. The method has also been implemented in a large eddy simulation (LES) explicit solver and successfully tested for an acoustic wave propagation. It has been finally used in the two-dimensional large eddy simulation of a turbulent jet.

I. Introduction

Thanks to the progress of turbulent combustion modeling and to the increase in computing power, today numerical simulations are now part of the design process of aeronautical turbine combustors. Main challenges when designing a gas turbine combustor are to control flame ignition and extinction as well as combustion instabilities that may cause the destruction of the engine. Predicting these phenomena requires to take into account flow unsteady effects and their interaction with the reaction zone. The Large Eddy Simulation (LES) technique, combined with significant improvements in the field of combustion modeling, have provided researchers with the tools needed to study a wide variety of turbulent reacting flows. Although LES has the capability to accurately predict the key phenomena, some studies remain beyond the capacities of existing CFD tools. For example, simulations of a turbine combustor are usually limited to a single sector. However, to accurately predict ignition, reignition, quenching and instabilities taking place in a combustor, computations of the full chamber are required. Recently, LES of reacting flows in full combustion chambers of gas turbines have been successfully performed with the AVBP code on several massivelly parallel machines. [START_REF] Allsopp | Unfolding the IBM eServer Blue Gene Solution[END_REF][START_REF] Boileau | Toward LES of an ignition sequence in a full helicopter combustor[END_REF][START_REF] Boileau | LES of an ignition sequence in a gas turbine engine[END_REF][START_REF] Wolf | Massively parallel LES of azimuthal thermo-acoustic instabilities in annular gas turbines[END_REF][START_REF] Staffelbach | Large Eddy Simulation of self excited azimuthal modes in annular combustors[END_REF] However, such calculations are still too expensive to be used in industry.

So, although the increasing efficiency of parallel computing, computational time remains an important issue. The CPU time directly relies on restrictions imposed by the numerical schemes stability. When considering an explicit solver for fully compressible Navier-Stokes equations, the time step is constrained through an accoustic CFL-like condition by the smallest cells in the overall computational domain. For combustion chamber simulation, local refinement of the fuel injector leads to an inhomogeneous mesh and the former restriction affects the overall solver efficiency. One way to address this issue is to use local time-stepping (LTS) methods. They allow to define local time steps whose values depend on the corresponding local mesh properties. The aforementioned severe CFL restrictions may then be overcome.

Existing LTS methods mainly group in two classes: hierarchical and non-hierarchical methods. The category of hierarchical methods counts AMR-like strategies and adaptative multiscale schemes. AMR methods were introduced in the 80s by Berger and co-workers [START_REF] Berger | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF][START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF] and are widely used and suitable for parallel solvers using structured meshes. They are based on a hierarchical tree of overlapping grids. Adaptative multiscale schemes [START_REF] Coquel | Local time stepping for implicit-explicit methods on time varying grids[END_REF][START_REF] Müller | Fully adaptive multiscale schemes for conservation laws employing locally varying time stepping[END_REF] are multiresolution techniques in the line of Harten's work. [START_REF] Harten | Multiresolution algorithms for the numerical solution of hyperbolic conservation laws[END_REF] On the contrary, other LTS methods are rather concerned by interfacing neighboring grids that use different time steps. In this category, one can find some domain decomposition methods for treating evolution problems that base the coupling between subdomains on Schwarz waveform relaxation methods. [START_REF] Gander | Optimal Schwarz waveform relaxation for the one dimensional wave equation[END_REF][START_REF] Martin | Méthode de décomposition de domaine et de couplage pour des problèmes d'évolution[END_REF] However, these mathematical studies are not ready for application to non-academic multidimensional problems. Still in this second category, one can find methods [START_REF] Crossley | Time accurate local time stepping for the unsteady shallow water equations[END_REF] in the line of a pioneering work by Osher and Sanders. [START_REF] Osher | Numerical approximations to nonlinear conservation laws with locally varying time and space grids[END_REF] Their work on one-dimensional scalar conservation laws has been extended in [START_REF] Dawson | High resolution schemes for conservation laws with locally varying time steps[END_REF] that show the possibility for their LTS method to reach a second-order in time. Other LTS methods have also been designed for moving mesh methods. [START_REF] Tan | Moving mesh methods with locally varying time steps[END_REF] Most of papers use for tests some academic fluid flow problems, but one also find several publications from the electromagnetic community. [START_REF] Fumeaux | A generalized local time-step scheme for efficient FVTD simulations in strongly inhomogeneous meshes[END_REF] The objective of this work is to develop a new LTS method -called DECCOUP -adapted to explicit and parallel solvers. The principle is to balance the computing load by splitting the computational domain into subgrids which size depends on the local time step. To the authors' knowledge, DECCOUP is the first LTS method suited for massively parallel simulations except AMR-like methods. [START_REF] Berger | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF][START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF] The DECCOUP technique is developed in the scope of the ANR-CIS SIMTUR project. Final objectives are to implement it and to use it in the AVBP code. [START_REF] Avbp | AVBP Code[END_REF] AVBP is an unstructured code that solves the fully compressible LES equations for reacting flows. [START_REF] Moureau | Numerical methods for unsteady compressible multi-component reacting flows on fixed and moving grids[END_REF] Using a cell-vertex formulation, it can work with both structured and unstructured grids which makes it easily applicable to complex geometries. [START_REF] Schonfeld | A cell-vertex approach to local mesh refinement for the 3-D Euler equations[END_REF] Centered spatial schemes and explicit timeadvancement are used to control numerical dissipation. [START_REF] Rudgyard | Development of high-order Taylor-Galerkin schemes for LES[END_REF] Several numerical schemes are available ; among them are a second order finite-volume Lax-Wendroff scheme and a third order finite-element TTGC scheme.

Making an estimate based on the case illustrated on Fig. 1, one can broadly evaluate the speedup resulting from DECCOUP. Consider a mesh involving:

• n small cells from blocks advanced with δt small ,

• n large = R n • n small larger cells from blocks advanced with δt large = R δt • δt small ,

• n buf f er = λR n • n small cells from buffer blocks.

The role of these buffer blocks will be described later in the document. They are intermediate blocks whose cells are calculated twice: first taken as small cells, then taken as large ones. Without any LTS method, the overall load to be handled to reach the physical time t sim in the simulation writes:

Load = (n small + n buf f er + n large ) • t sim δt small (1)
while using DECCOUP it reduces to the following:

Load LT S = (n small + n buf f er ) • t sim δt small + (n buf f er + n large ) • t sim δt large (2)
Thus, assuming that computational time is proportional to the load, the following speedup is expected:

Speedup = Load Load LT S = 1 + (1 + λ)R n (1 + λR n ) + (1 + λ) Rn R δt (3)
Assuming λ ≪ 1, the former expression underlines that values of Speedup ranging from 1 (in the case corresponding to R n = R δt = 1) to δt large /δt small (in the ideal case obtained when R n → 0) can be expected. Intermediate cases are illustrated on the graph from Fig. 2. In practice, this speedup can be obtained on parallel computations for an optimal load balancing. Indeed, Eq. 3 is satisfied if the load is balanced on processors such as:

R δt n small N P,small = n large N P,large = (1 + R δt ) n buf f er N P,buf f er , (4)
where N P,small , N P,buf f er , N P,large are the number of processors handling small, buffer and large cells, respectively. The document is organized as follow. The strategy of DECCOUP is presented in section II. A test configuration and corresponding results are described in section III. Results from the implementation of DECCOUP in the LES code AVBP are presented in section IV. Finally, section V concludes on the present work.

II. The DECCOUP method

In most Navier-Stokes explicit and parallel solvers, a domain decomposition method is used to split the overall computational domain into several blocks (noted B j ). Each block is then affected to one processor. Generally, the initial solver runs over all blocks using a unique global time step, determined through a reduction operation (each processor evaluates the time step value allowed by the CFL criterion over its cells and a reduction operation keeps the most restrictive value). DECCOUP is designed to allow each block B j to keep its proper minimum time step δt j . The method is described using 1D examples for clarity but the generalisation to multidimensional problems is straightforward. Consider a given spatial numerical scheme in an explicit context. To be applied at bordering cells of a block, the corresponding stencil generally requires data from neighbouring blocks at the same time t n . These data can be made available by using ghost cells surrounding each block and are updated from communications between processors at the begining of each new iteration. When using a unique time step for each block, no approximation is needed to fill in ghost cells since all blocks are simultaneously advanced in time. On the contrary, when using LTS, two neighbouring blocks may not use the same time steps to advance the solution. Consequently, the data required to fill in ghost cells may not be available at each discrete time values. Ghost cells values are then approximated between two communications steps (see Fig. 3).

= m j • δt j , ∀j ∈ [1, N B ],
where m j are integers, exists. Times corresponding to the global time step are noted (t n ) n=1,N∆t , where N ∆t is the total number of iterations using ∆t. By construction, t n is common to all block solutions. Intermediate times between t n and t n+1 using the local time step δt j are expressed by t = t n + ν • δt j , ν = 1, m j . The computed B j block solution vector at time t is noted w(B j , t). Ghost cells solution values belonging to B j and shared by B j and B k are noted w(B j ∩ B k , t).

At a given interface between two neighbouring one-dimensional blocks, two cases can be distinguished: either blocks using the same time step value so no LTS handling is required or blocks using different time steps so the DECCOUP strategy applies. In that case, DECCOUP requires a buffer block to be inserted at the LTS interface during a DECCOUP-adapted domain decomposition phase. This preprocessing step is detailed further. To illustrate the method, a sequence of 3 blocks B small , B buf f er and B large is considered. Their corresponding parameters are all noted with the appropriate subscript (see Fig. 4). In the present case, δt small < δt large and m small • δt small = m large • δt large . The buffer block B buf f er is dedicated to advance the solution twice from t n to t n+1 . Grid spacing on this block is designed in such a way that stability criterion is respected for both time step values δt small and δt large . The solution is computed firstly through m small subcycle iterations with δt small = ∆t/m small , that lead to a first solution w(B buf f er , t n + m small • δt small ), and secondly using δt large = ∆t/m large , leading to a second solution w(B buf f er , t n + m large • δt large ). If δt small is used to advance the buffer block, interior block data are available to compute w(B small , t n+1 ) but ghost cells data w(B large ∩ B buf f er , t n + ν • δt large ), ν = 1, m large are missing. They are roughly approximated by:

w(B large ∩ B buf f er , t n + ν • δt large )) = w(B large ∩ B buf f er , t n ), ∀ν = 1, m large .
(5)

Here DECCOUP takes advantage of the explicit nature of the solver. Indeed, given an explicit solver that use a stencil requiring s lef t (resp. s right ) cells on the left (resp. right) of the current point, any perturbation introduced at time t n cannot propagate through the grid of a distance exceeding -per time step δt j -s lef t cells to the right and s right cells to the left a . Consequently, errors introduced by the previous approximation (Eq. 5) propagate over a known number of cells from the B large ∩ B buf f er ghost cells. This number of cells is equal to s lef t • m j at time t n+1 . The corresponding solution w(B buf f er , t n + m small • δt small ) is then polluted by errors near the interface with B large but remains unaffected in the vicinity of the interface with B (see Fig. 5). When using δt large to advance B buf f er in time, enough data are available to obtain w(B large , tn+1) and a partially valid solution on B buf f er at the same time.

When δt large is used to advance B buf f er , the same difficulty occurs at the interface B buf f er ∩ B small (see Fig. 6): w(B large , t n+1 ) is obtained, as well as the polluted solution w(B buf f er , t n + m large • δt large ). Finally, the way of interfacing grids using different time steps whithin DECCOUP consists in using an intermediate grid. Advancing the solution with δt small allows to locally preserve a continuity in the accuracy of the solution through the interface between B small and B buf f er but introduces errors at the other interface, whereas the use of the biggest time step locally preserve a continuity in the accuracy of the solution through the interface between B small and B buf f er . On each cell of B buf f er , two estimations of the solution w(B buf f er , t n+1 ) are available. Both of them are polluted by errors but, if B buf f er is long enough, cells affected by errors in the first computation do not overlap with the polluted ones during the second computation. It is then possible to built w(B buf f er, t n+1 ) as a smooth transition between the two polluted solutions that keeps only the unaffected part of them and filters all wrong values. This is done using the following equation:

w(x i , t n+1 ) = (1 -p(x i )) • w(x i , t n + m small • δt small ) + p(x i ) • w(x i , t n + m large • δt large ) , (6)
where x i ∈ B buf f er and the function p(x) can be an hyperbolic tangent located at the center x 0 of B buf f er p(x) = 1 2 (1 + tanh (A(x -x 0 ))) , so as to keep the not polluted part of each solution.

a This is true if the spatial scheme is explicit and cannot be applied if implicit schemes such as compact schemes [START_REF] Lele | Compact finite difference schemes with spectral-like resolution[END_REF] are used.

III. Evaluation of the method in a test code

In order to evaluate the accuracy of the DECCOUP method, calculations are performed using a test code which solves the one-dimensional convection equation: ∂ρ ∂t + U ∂ρ ∂x = 0 with U = cst > 0. Periodic boundary conditions are used and the temporal integration is explicit. The initial solution is a Gaussian wave in a domain of length L B . This preliminary study focuses only on the accuracy of the underlying schemes and not on the speedup of the method. The overall domain of calculation is splitted into 6 blocks, involving two successive LTS interfaces with the following alternance: B small , B buf f er , B large , B large , B buf f er , B small . Thus, the wave crosses LTS-interfaces twice a period, once from small to large time steps , and a second time from large to small time steps. Here, a value of m large = 1 is considered, leading to ∆t = δt large , m small = 8. The number of iterations is set so as to end the simulation at t sim = k • (L B /U ), i.e. until the wave has traveled k times through the overall domain.

The code developed for testing DECCOUP is used in 8 configurations, depending on the choices for the scheme (HOSRK4/WENO5RK3), for the use of LTS (with/without) and for the nature of the mesh (uniform/non-uniform). Main input parameters are:

• a level L i to define the size of mesh cells,

• the number of subcycles m small ,

• the wave speed U ,

• a CF L number

• the number of flow-through times k.

A set of meshes L 0 to L 4 is used to compare the solution for various levels of refinement (the width of the Gaussian wave being fixed once for all). These meshes are built on a root one (L 0 ), called level 0, and they are defined recursively such that L n+1 is obtained by splitting every cells of L n in two identical cells (see Fig. 7). Level 0 characteristics, namely the cell size δx large on B large (and B buf f er ) as well as cells numbers n small , n large and n buf f er , are hardcoded. The cell size δx small on B small is equal to δx large for the uniform mesh and equal to δx small = δx large /m small in the non-uniform case. For level L i , δx values are divided by 2 i while cells numbers are multiplicated by the same factor, insuring that block lengths are the same for all levels to allow comparisons. In the non-uniform case, B buf f er and B large remains uniform but B small is meshed with cell sizes ranging from δx large /m small to δx large through an arithmetic progression. In terms of spatial resolution, L 0 is such that the solution goes from almost 0 to 1 over 5 cells (and over 10, 20, 40 and 80 cells for L 1 , L 2 , L 3 and L 4 , respectively). Concerning time steps, δt large is set using the CF L and wave speed U (δt large = CF L • δx large /U ). When LTS is used δt small = δt large /m small . If not, δt small is equal to δt large . Two schemes are retained: a high order finite differences scheme 23 having a third order in space even on non-uniform grids coupled to a RK4 time integration method (HOSRK4), as well as a fifth order finite volumes scheme of the WENO type 24 coupled to a TVD-RK3 time integration method (WENO5RK3). This WENO scheme is also able to handle non-uniform meshes. The run parameters are: U = 1, CF L = 0.5, δt large /δt small = 8, where CF L is based on δt small for the uniform mesh cases. A solution with DECCOUP for the convection equation is shown on Fig. 8. After the wave has travelled 10 times through the periodic domain, the DECCOUP solution perfectly matches the theoretical one. For refinement levels L 0 to L 4 , the accuracy of the solution is compared when DECCOUP is used and when it is not. The error is defined as the difference -based on a L ∞ -norm -between the computational solution and the exact one (i.e. the Gaussian wave translated by U • (t -t 0 )), following the equation: Three main results appear on Fig. 9 and 10: i/ DECCOUP do not necessary introduces a loss of accuracy; ii/ the impact of the underlying scheme on performances of DECCOUP is significant; and iii/ the method does not seem to be sensitive to the non-uniformity of the mesh. For HOSRK4 scheme, the accuracy of the solution appears to be either better when DECCOUP is used, in the case of a uniform mesh or unchanged, in the non-uniform case. The performance of DECCOUP also appears to be dependent of the scheme, since results with the WENO scheme are not as accurate as those corresponding to the HOSRK4 scheme. On Fig. 11, scheme orders with and without DECCOUP have been computed for both schemes. The HOSRK4 scheme appears to be very efficient for all considered levels of refinement, showing a better scheme order with DECCOUP than without. As noticed previously, the method is less accurate when used with the WENO scheme. In that case, DECCOUP introduces a slight loss in the scheme order. Figure 13 shows that a very good agreement is obtained between the AVBP solution calculated with DECCOUP and the one calculated without. Another result of interest is the speedup obtained with DEC-COUP. When no LTS method is used, the CPU time required to reach a given physical time t sim on 5 processors is 4.95 times longer than when LTS is used. Actually, the LTS method requires 2 additional processors to handle the duplicated buffer blocks. The CPU cost due to these additional processors is included in the calculation of the speedup using a correction factor of 5/7. The corrected speedup obtained with DECCOUP is thus equal to 3.53 while the theoretical value given by Eq. 3 (using the same correction) is 3.25. DECCOUP has also been tested in the 2D large eddy simulation of a plane air jet. This jet is characterized by an inlet Reynolds number of Re 0 = 30 000, an inlet Mach number of M 0 = 0.1 and a surrounding co-flow at Mach number 0.057M 0 . The subgrid stress tensor is modelled by the filtered Smagorinsky model. [START_REF] Ducros | Large-eddy simulation of transition to turbulence in a boundary layer developing spatially over a flat plate[END_REF] A third order finite-element TTGC scheme 21 using a two-step time integration is applied to an unstructured triangular grid (see Fig. 14a). A passive scalar is injected at the inlet to visualize the turbulent jet mixing. The aim of this test case is to see how DECCOUP may affect of the flow solution in a turbulent configuration. A simple decomposition in three transverse blocks is performed as shown by Fig. 14a. According to the cell size ratio between large and small cells blocks, a ratio of δt large /δt small = 10 is used. The CF L number based on the smallest cell is equal to 0.7. At some time t 0 of the jet simulation, the DECCOUP method is started and performed during t sim = 49 flow times (Fig. 14b) corresponding to 30 000 time iterations of δt small . The scalar field concentration at t 0 + t sim is compared to the solution obtained at the same physical time without using DECCOUP (Fig. 14c). Very small differences in the scalar transport can be seen between the basic and the DECCOUP solutions. These differences may be due to the growth of rounding errors in parallel large eddy simulations, an effect that has been identified by Senoner et al. [START_REF] Senoner | Growth of Rounding Errors and Repetitivity of Large-Eddy Simulations[END_REF] Figure 14b shows that the solution is continuous and regular inside, outside and at the boundaries of the buffer block where the exchange of the flow solution takes place between the small and large time steps regions.

error L∞ (t n ) = max i |ρ n i -ρ(x i , t n )| max i |ρ(x i , t n )| (7)

V. Conclusion

An acceleration method for LES explicit solvers based on local time-stepping has been presented. This method called DECCOUP uses a buffer block to handle interfaces between each LTS zone. DECCOUP has first been evaluated using a test code in a one-dimensional configuration. Results prove that the method has a minor incidence on the accuracy of the solution, whatever the nature of the underlying scheme. DECCOUP has also been implemented in the LES code AVBP. First academic tests shows that DECCOUP keeps the accuracy of the original numerical method while providing a speedup which is in agreement with the theoretical value. A 2D large eddy simulation test case has shown that DECCOUP is also accurate in turbulent configurations.

In practical applications, DECCOUP will be particularly efficient when only a few small cells penalize the overall calculation. Moreover, two requirements have to be fullfilled by the domain decomposition for DECCOUP: i/ the size of buffer blocks has to be large enough to allow to built a valid solution from nonpolluted parts of intermediate solutions, and ii/ a relation between the small-size, large-size and buffer cells have to be verified to ensure an efficient load balancing. This load balancing can be achieved using a multi criteria mesh decomposition algorithm.

Figure 1 .

 1 Figure 1. Cell sizes distribution in a combustion chamber mesh.

Figure 2 .

 2 Figure 2. Speedup isolines in the R δt -Rn plane.

Figure 3 .

 3 Figure 3. Approximations are required when using LTS : blocks Bj-1 and Bj use the same time step δt while Bj+1 uses a larger time step ∆t = 3δt. Values to fill in the ghost cells in Bj at its interface with Bj+1 are missing so they must be approximated.

Figure 4 .

 4 Figure 4. Within DECCOUP, a buffer block is used at LTS interfaces between B small and B large .

Figure 5 .

 5 Figure 5. When using δt small to advance B buf f er in time, enough data are available to obtain w(B small , tn+1) and a partially valid solution on B buf f er at the same time.

Figure 6 .

 6 Figure 6.When using δt large to advance B buf f er in time, enough data are available to obtain w(B large , tn+1) and a partially valid solution on B buf f er at the same time.

Figure 7 .

 7 Figure 7. First refinement levels of a schematic LTS interface, in the uniform case, and the Gaussian wave profile with a fixed-width. Here m small = 3.

Figure 8 .

 8 Figure 8. Comparison between the solution obtained by DECCOUP with the WENO scheme (•) and the initial solution (• • •) after exactly 10 periods. The mesh is uniform (level L1), U = 1, CF L = 0.5 (based on δt small ), δt large /δt small = 8.

Figure 9 .

 9 Figure 9. Error as a function of the refinement level Li for uniform meshes.

Figure 10 .

 10 Figure 10. Error as a function of the refinement level Li for non-uniform meshes.

Figure 11 .

 11 Figure 11. Order of the method as a function of the refinement level Li for HOSRK4 and WENO schemes.

Figure 12 .

 12 Figure 12. Decomposition of the AVBP computational domain for the acoustic wave test case.

Figure 13 .

 13 Figure 13. Comparison between the AVBP solution (pressure field in P a) obtained with (-) and without (•) DECCOUP after 20 000 time iterations.

Figure 14 .

 14 Figure 14. Snapshots of passive scalar fields at initial time t0 (Fig. a.) and current time t0 + t sim where tsim = 49 flow times (corresponding to 30 000 time iterations). The flow at t0 + t sim have been calculated with DECCOUP (Fig. b.) and without DECCOUP (Fig. c.).Figure a. shows the mesh grid (zoom) and a three-block decomposition: a small-step block, a buffer block and a large-step block. The gray scale has been saturated at 42% of the inlet value.

VI. Acknowledgments

This work was founded by the Agence Nationale de la Recherche through the scientific program ANR-CIS SIMTUR. Calculations were performed using HPC resources from GENCI-IDRIS (Grant 2009 -i2009080164). The authors gratefully acknowledge the assistance of G. Staffelbach from CERFACS.