Bastien Duboc

Guillaume Ribert
email: guillaume.ribert@coria.fr

Pascale Domingo

Hybrid transported-tabulated chemistry for partially premixed combustion

Keywords: Hybrid chemistry, Partially premixed, Edge flame, Methane

The integration of combustion chemistry into a fully compressible numerical solver is presently achieved using the hybrid transported-tabulated chemistry (HTTC). With HTTC, the main species are transported while most minor species are tabulated, which means that differences with a fully transported chemistry (FTC) solver are limited and concern mainly table reading for minor species. The implementation steps of HTTC are given in detail and an optimization of the code is proposed by tabulating the properties of the pure species as well as the reaction rates of the elementary reactions as a function of the temperature to speed up simulations. The original version of HTTC, validated for premixed combustion, has been also extended to partially premixed combustion by adding a prolongation of the lookup table for minor species outside the flammability limits. Two strategies are proposed and evaluated on a methane / air edge flame featuring a very high mixing fraction gradient. The results agree favorably by comparison with a reference flame simulated with a detailed chemistry. As the minor species are no longer transported with the flow using HTTC, the calculation cost is found divided by about 5 compared to the FTC solver.

Introduction

The introduction of detailed kinetic mechanisms into numerical codes for combustion is essential for the prediction of the flame behavior in realistic simulations. It consists in calculating the diffusive properties of chemical species as well as their chemical reaction rate in the flame front in order to correctly predict the macroscopic properties of the flame, such as the flame speed, the temperature and the composition of gases.

Two major difficulties make combustion challenging for the numerical simulation of fluid flow. The first is the large number of variables to be transported associated with the large number of elementary reactions involved in the kinetic mechanism [START_REF] Lu | Toward accommodating realistic fuel chemistry in large-scale computations[END_REF]. The second lies in the difference between the macroscopic spatial and temporal scales of the system to be simulated and those of chemistry. In the typical example of the engines, the characteristic length of a combustion chamber is of the order of ten centimeters, while the time required for the flow to cross the chamber is a few tenths of seconds. These quantities have to be compared to the spatial scales of a flame, which vary between a few microns and a few millimeters, and time scales that are spread over several orders of magnitude and can fall below 10 -11 seconds (Fig. 1). Given these difficulties, Figure 1: Width of mass fraction profiles and chemical time characteristic of minor species, calculated in 1D methane/air [START_REF] Lindstedt | Modeling of the chemical complexities of flames[END_REF] and decane/air [START_REF] Luche | Reduction of large detailed kinetic mechanisms : application to kerosene/air combustion[END_REF] flames.

the treatment of chemistry in DNS/LES codes is still the subject of many studies [START_REF] Pope | Small scales, many species and the manifold challenges of turbulent combustion[END_REF].

As the detailed schemes are generally too heavy for practical simulations, an interesting idea is to develop reduced schemes by identifying the important species using Quasi-Steady State Approximation (QSSA), Partial Equilibrium (PE) or Rate-Controlled Constrained Equilibrium (RCCE) techniques [START_REF] Keck | Rate-controlled constrained-equilibrium theory of chemical reactions in complex systems[END_REF][START_REF] Lam | The CSP method for simplifying kinetics[END_REF][START_REF] Lu | A directed relation graph method for mechanism reduction[END_REF][START_REF] Pepiot-Desjardins | An automatic chemical lumping method for the reduction of large chemical kinetic mechanisms[END_REF][START_REF] Goussis | Quasi steady state and partial equilibrium approximations: their relation and their validity[END_REF][START_REF] Jaouen | Automatic reduction and optimisation of chemistry for turbulent combustion modelling: Impact of the canonical problem[END_REF][START_REF] Koniavitis | A methodology for derivation of RCCE-reduced mechanisms via CSP[END_REF].

The obtained kinetics compare well with the original detailed mechanism, but on a much lower range of validity. In order to keep the complete information of the kinetic scheme in flames computation, an alternative is the use of the tabulation of chemical responses of canonical combustion problems such as laminar premixed or non-premixed flames [START_REF] Bradley | Laminar flamelet modeling of recirculating premixed methane and propane-air combustion[END_REF][START_REF] Maas | Simplifying chemical kinetics: intrinsic low-dimensional manifolds in composition space[END_REF][START_REF] Peters | Turbulent combustion[END_REF][START_REF] Van Oijen | Modelling of premixed laminar flames using flamelet-generated manifolds[END_REF][START_REF] Gicquel | Laminar premixed hydrogen/air counterflow flame simulations using flame prolongation of ildm with differential diffusion[END_REF]. For instance with the FPI technique [START_REF] Gicquel | Laminar premixed hydrogen/air counterflow flame simulations using flame prolongation of ildm with differential diffusion[END_REF], the flame structures of premixed flamelets are computed with the full detailed kinetic scheme and the species mass fraction, chemical source terms, the temperature, etc., are stored in a look-up table as a function of a progress variable (Y c ) and a mixture fraction (Z) before the simulation of the target configuration. For the latter, only Y c and Z need to be transported with the flow, resulting in a significant reduction in the computational cost. The description of chemistry is then linked to a database whose size can be very large and not adapted to the context of high-performance computing. This problem was addressed in [START_REF] Ribert | Tabulation of complex chemistry based on self-similar behavior of laminar premixed flames[END_REF][START_REF] Wang | Self-similar behavior and chemistry tabulation of burnt-gas diluted premixed flamelets including heat-loss[END_REF][START_REF] Ribert | A multi-zone self-similar chemistry tabulation with application to auto-ignition including cool-flames effects[END_REF] for laminar flamelets or ignition phenomena using a self-similarity concept (S2FT). This property has been exploited by Ribert et al. [START_REF] Ribert | Hybrid transported-tabulated strategy to downsize chemistry for numerical simulation of premixed flames[END_REF] for a strategy, called Hybrid Transported-Tabulated Chemistry (HTTC), combining the detailed-chemistry solving for the main species with the tabulation of species with very small concentrations using S2FT.

The objectives of this paper is to describe in detail the introduction of HTTC into numerical compressible solvers for Direct Numerical Simulations (DNS) or Large-Eddy Simulations (LES) and to extend the methodology to partially premixed combustion. The DNS of a methane/air edge flame [START_REF] Buckmaster | Edge-flames[END_REF] serves as a numerical challenge for the validation of HTTC compared to the fully transported chemistry.

Coupling compressible solvers with hybrid transported-tabulated chemistry

In the original work of Ribert et al. on HTTC [START_REF] Ribert | Hybrid transported-tabulated strategy to downsize chemistry for numerical simulation of premixed flames[END_REF], the method was implemented and tested in the steady, isobaric one-dimensional numerical code REGATH [START_REF] Ribert | Counterflow diffusion flames of general fluids: Oxygen/hydrogen mixtures[END_REF][START_REF] Pons | Mass transfer and combustion in transcritical non-premixed counterflows[END_REF] for methane-air flames only. Its implementation into a fully compressible DNS numerical code requires additional developments that are now explained by comparison with a fully transported chemistry (FTC) solver. Classic notations are used for the time (t), spatial coordinates (x i ), density (ρ), pressure (P), temperature (T ) and velocity components (u i ) in the x iaxis. The total number of species present in the detailed kinetic mechanism is denoted N sp , and the total number of chemical reactions is N r . The mass fraction of species k is writen Y k with k ∈ [1, N sp ], and the full

vector of species is Y : Y = (Y 1 , • • • , Y N sp ).

Fully transported chemistry (FTC) solver

Neglecting external forces and energy sources, any fully transported chemistry solver considers the balance equations of momentum, ρu, species density k, ρ k (ρ k = ρY k ) and energy. The total non chemical energy, E, is presently used and corresponds to the sum of the kinetic energy, e c , and sensible energy, e s : E = e s + e c . The diffusion velocity is modeled with the Hirschfelder and Curtiss approximation [START_REF] Hirschfelder | Molecular theory of gases and liquids[END_REF] and a correction velocity [START_REF] Coffee | Transport algorithms for premixed, laminar steady-state flames[END_REF] is used to ensure mass conservation. The reaction rate of species k ( ωk ) is evaluated from Arrhenius rate expressions.

The simulations presented thereafter have been performed using the parallel numerical code SiTCom-B [START_REF] Bouheraoua | Large eddy simulation of a supersonic lifted jet flame: Analysis of the turbulent flame base[END_REF][START_REF] Petit | Large-eddy simulation of supercritical fluid injection[END_REF][START_REF] Guven | Large-eddy simulation of supersonic H 2 /O 2 combustion: application to a rocket-like igniter[END_REF] that solves the fully compressible multi-species Navier-Stokes equations. It is based on an explicit Finite Volumes scheme written for cartesian grids. The convective terms are computed resorting to a fourth-order centered skew-symmetric-like scheme [START_REF] Ducros | High-order fluxes for conservative skew-symmetriclike schemes in structured meshes: application to compressible flows[END_REF] and the diffusive terms are computed using a fourthorder centered scheme. The time integration is performed using a fourth order Runge-Kutta method [START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes[END_REF].

The spatial discretization scheme is augmented by a blend of second-and fourth-order artificial dissipation terms [START_REF] Tatsumi | Flux-limited schemes for the compressible Navier-Stokes equations[END_REF][START_REF] Swanson | On central-difference and upwind schemes[END_REF][START_REF] Swanson | On some numerical dissipation schemes[END_REF]; these terms are added in order to suppress spurious oscillations and damp high-frequency modes. A sensor based on pressure and density gradients insures that the artificial dissipation is applied only to zones of interest, i.e. where either strong gradients of density or pressure, which cannot be resolved by the mesh, are detected [START_REF] Petit | Large-eddy simulation of supercritical fluid injection[END_REF]. The treatment of boundary conditions is performed with the method of NSCBC (Navier Stokes Characteristic Boundary Conditions [START_REF] Petit | Framework for real-gas compressible reacting flows with tabulated thermochemistry[END_REF][START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF]).

In the FTC solver of SiTCom-B, a temporal iteration starts by the resolution of the balance equations for ρu, ρY k and ρE. Then, thermodynamics variables and transport properties are updated for the next temporal step:

1. The density ρ is given by ρ =

N sp k=1 ρ k , ρ k = ρY k . 2.
Primitive variables such as the velocity vector U, E and Y are obtained dividing by ρ the conservative variables, ρU, ρE and ρY, respectively.

3. The kinetic energy, e c , is computed from the velocity field u i , and e s is given by: e s = Ee c .

4. The temperature, T , is computed from the sensible energy.

5. The pressure is updated through the equation of state: P = ρRT/W, with R = R u /W. 6. Thermodynamics and transport properties as well as the speed of sound are finally updated.

Hybrid transported-tabulated chemistry solver

Introducing the HTTC approach into a numerical code is simplified if a FTC solver is already present. Indeed, the two solvers share many functionalities such as treatment of boundary conditions or the handling of the chemical scheme. However, with HTTC [START_REF] Ribert | Hybrid transported-tabulated strategy to downsize chemistry for numerical simulation of premixed flames[END_REF], the major species, Y M , defined as Y k | k ∈ M with M = {k | species k is transported} are transported with the flow whereas the mass fraction of minor species

Y m (Y k | k ∈ m with m = {k | species k is tabulated}) are extracted from a look-up table.
Then, the focus is on tasks that differ from the FTC solver. For a better understanding, all the steps described thereafter are summed up in table 1.

• The initialization step is nearly identical for the two solvers (FTC and HTTC). The thermodynamical data of species involved in the kinetic mechanism and the Arrehnius coefficients of the chemical reaction are fetched in a data file and the primitive variables (P, T , U and Y) are read from an initial solution.

The HTTC solver requires an additional step, to read and store in the memory of every MPI process the look-up table that contains the tabulated minor species mass fractions (Y m ). The data Y m is accessible via the parameters or coordinates of the look-up table and represent physical phenomena (progress of reaction, mixing, etc.). By following the classic approach of tabulated chemistry [START_REF] Gicquel | Laminar premixed hydrogen/air counterflow flame simulations using flame prolongation of ildm with differential diffusion[END_REF][START_REF] Van Oijen | Modelling of premixed laminar flames using flamelet-generated manifolds[END_REF], it is common to use the progress variable (Y c ), mixture fraction (Z), enthalpy (h), etc. as table coordinates. Y c is often defined as Y c = Y CO + Y CO 2 for methane/air flames but an optimisation process may be used for a better description of the lookup table [START_REF] Ribert | Hybrid transported-tabulated strategy to downsize chemistry for numerical simulation of premixed flames[END_REF]. Computation of remaining variables (diffusion coefficients, T , P, thermal capacities, gradients, etc.), using the whole set of species, or the transported species only, if HTTC "fast" is used. Building the look-up table is detailed in section 3.

• Once the pre-processing step is performed, the temporal loop starts. At the beginning of each temporal iteration, the value of the time step is calculated for both solvers as the minimum value of the acoustic and r tol = 0.1, ωk being the chemical source term. However for the HTTC solver, the tabulated species are not taken into account in these calculations, since they are not transported. This point is essential, since the chemical time step requirement to solve for the tabulated species (i.e. minor species) would be several orders of magnitude below those for transported major species [START_REF] Maas | Simplifying chemical kinetics: intrinsic low-dimensional manifolds in composition space[END_REF]. This is actually the main interest of HTTC, that is responsible for the largest part of the savings on the computational cost.

• Then, the transport of conservative variables is the same for both solvers. Updating the different quantities of the flow follows the procedure given in section 2.1 for FTC. However, to ensure the mass conservation with the solver HTTC, ρ at the current time step n, is given by:

ρ n = k∈M ρ n k + k∈m ρ n-1 k = k∈M ρ n k + k∈m ρ n-1 Y n-1 k with        k∈m ρ k        n-1        k∈M ρ k        n . (1)
With HTTC, the sum of the mass fractions of the full species set is then slightly different from 1:

N sp k=1 Y k = k∈M Y k + k∈m Y k = 1 + (2)
where is a non-zero real number. Following [START_REF] Duboc | Description of kerosene / air combustion with hybrid transported-tabulated chemistry[END_REF][START_REF] Duboc | Evaluation of chemistry models on methane/air edge flame simulation[END_REF], every species mass fraction of the full set is corrected to ensure mass conservation:

Y cor k = Y k 1 + ε , ∀k ∈ [1, N sp ] (3)
where Y cor k is the corrected value of the mass fraction of the species k.

Knowing ρ, the mass fraction of major species is given by Y M = ρ M /ρ. The progress variable, Y c , used to build the look-up table, as well as the mixture fraction, Z, are calculated using Y. Then, the tabulated mass fractions Y m are read in the table.

• Finally, the local value of the pressure is computed using the perfect gas law P = ρRT , where ρ, R and T mainly depends on the transported species, because the mass fraction of the tabulated species is small [START_REF] Ribert | Hybrid transported-tabulated strategy to downsize chemistry for numerical simulation of premixed flames[END_REF].

Finally, all the remaining variables are computed using the full set of species, or the transported species only, depending whether the standard or the fast version (see Section 2.3) of the HTTC solver is selected.

For FTC, these variables are always built with the full set.

Variable

Standard version Fast version

Density ρ = ρ(Y) ρ = ρ(Y) Mixture fraction Z = Z(Y) Z = Z(Y) Diffusion coefficients µ = µ(Y), λ = λ(Y), D k = D k (Y), k ∈ M µ = µ(Y M ), λ = λ(Y M ), D k = D k (Y M ), k ∈ M Molar mass W = W(Y) W = W(Y M ) Thermal capacities c v = c v (Y), c p = c p (Y) c v = c v (Y M ), c p = c p (Y M ) Temperature T = T (Y) T = T (Y M ) Pressure P = ρ(Y) R W(Y) T (Y) P = ρ(Y) R W(Y M ) T (Y M ) Speed of sound c = c(Y) c = c(Y M ) Correction velocity V c = V c (Y) V c = V c (Y M ) Enthalpy diffusion term in the energy equation N sp k=1 h s,k Y k V k,i k∈M h s,k Y k V k,i
Table 2: Differences in the calculation of several variables with the "standard" and the "fast" version of the HTTC solver of SiTCom-B.

"Fast" formulation of the HTTC solver

A considerable amount of computational time is saved with HTTC, by increasing the global time step allowed by the fully explicit detailed chemistry solver. As a consequence, less iterations are required to reach a given physical time. However, with the "standard" version of the HTTC solver presented above, the cost of an iteration is almost identical to the cost of an iteration of the FTC solver. Although the tabulated species are no longer transported, all species are still needed to solve the conservation equations of the transported species to compute, for instance, the transport coefficients and correction velocity.

A simplified version of the original HTTC approach, called "fast" HTTC, is now presented. In this version, the set M of the transported species is used instead of the full set M + m whenever possible, with the aim of decreasing the duration per iteration, and consequently, the global cost of the simulation. In the fast version, the source terms of the transported species and the heat release rate are naturally still computed with the whole set of species. This is a crucial point of the HTTC method in order to take the full detailed kinetic mechanism into account during the simulation. However, several quantities are computed using only the transported species set, as shown in Table 2. Those quantities are selected carefully though, not to violate the mass conservation in the flow, and not to jeopardize the robustness of the solver. With this constraint in mind, the density ρ and the mixture fraction Z are calculated using the full set of species, as it is done for the standard version. Since the sum of the mass fractions of the tabulated species is very small compared to the transported species (< 1% for methane flames and < 5% for kerosene flames), the contribution of the tabulated species to the other quantities presented in Table 2 (diffusion coefficients, thermal capacities, speed of sound, etc.) is not significant, and consequently the latter can be safely computed using the transported species only.

Minor species therefore have a very low concentration and the largest part of the total mass of the flow is due to a small number of species (reactants and main products). The main species such as O 2 , H 2 O, CO, CO 2 , N 2 , etc. are generally present in all mechanisms, and only minor species will differ from one mechanism to another. In practice with HTTC, the number of transported species does not exceed 15, whatever the size of the mechanism, from a few tens to several hundreds of species. Consequently, with the fast version of HTTC the cost per iteration is expected to stay of the same order of magnitude when the mechanism grows, in spite of an additional cost due to a larger number of chemical reactions, and a larger table. As detailed hereafter, the impact of such extra computations can be limited by using a proper code optimization. Thus, the cost per iteration is expected to moderately increase even if the size of the mechanism largely grows.

A direct comparison between the standard and fast formulations of HTTC is provided in Fig. 2 for a stoichiometric mono-dimensional methane/air flame. A very good agreement is found with the FTC formulation. Lines : results from SiTCom-B using different solvers (black solid : FTC, blue solid : HTTC, blue dashed : fast HTTC).

Code optimization

Numerical simulation of reactive flows relies on extremely powerful and costly computer resources which must be employed in a sensible approach. First, a fine optimization of the most time-consuming parts of the code becomes essential to ensure correct performances. The duration of a calculation is controlled by three parameters: 1/ the minimum time step to ensure the stability of the numerical schemes, i.e. presently the chemical time step for reacting multi-species flows, 2/ the mesh size and 3/ the time taken by the code to perform one iteration. For the latter, Lu and Law [START_REF] Lu | Toward accommodating realistic fuel chemistry in large-scale computations[END_REF] propose to decompose the reaction rate according to the values taken by the temperature exponent and the activation temperature with a partial tabulation of exponential functions. In the present study, the reaction rates of each reaction, which depend only on temperature, are calculated once at the beginning of the simulation, for several temperature values uniformly distributed over a sufficiently wide interval to cover all the values that will be encountered during the simulation. The same approach is used to compute the pure species properties such as viscosities, µ k , which are present in the calculation of the mixture viscosity [START_REF] Wilke | A viscosity equation for gas mixtures[END_REF][START_REF] Bird | Transport phenomena[END_REF],

µ = N sp k=1 X k N sp j=1 X j Φ jk µ k with Φ jk = 1 √ 8 1 + W k W j -1 2         1 + µ k µ j 1 2 W j W k 1 4         2 , (4)
thermal conductivities [START_REF] Mathur | Thermal conductivity of binary, ternary and quaternary mixtures of rare gases[END_REF], λ k , binary diffusion coefficients, D jk or coefficients Φ jk [START_REF] Bird | Transport phenomena[END_REF]. A linear interpolation completes the procedure and avoids the calculation of complex, i.e. CPU consuming, formula coming from kinetic gas theory. In Fig. 3, the profiles of µ k , λ k , D jk and Φ jk associated to Lindstedt's mechanism [START_REF] Lindstedt | Modeling of the chemical complexities of flames[END_REF] (N sp = 29, N r = 141) are found smooth enough when described by a tabulation with a discretization of ∆T transp = 5 K, which is a good compromise between accuracy of the description of the transport properties and access time to the stored data. The amount of memory M transp reserved for this tabulation in double-precision of the transport coefficients is calculated as follows:

M transp = C o × T max -T min ∆T transp × 2N sp + 2N 2 sp , (5)
where T min = 200 K and T max = 3000 K delimit the temperature range used to construct the table. C o = 8 for 8-byte real numbers. µ k and λ k are considered evolving as N sp , and D jk and Φ jk as N 2 sp . Considering the kinetic scheme of Luche et al. [START_REF] Luche | Reduction of large detailed kinetic mechanisms : application to kerosene/air combustion[END_REF] for kerosene/air combustion which contains N sp = 91 chemical species (and 991 elementary reactions), M transp = 75 Mo, which is acceptable in the context of massively parallel simulations. Using the chemistry of Curran et al. [START_REF] Curran | A comprehensive modeling study of iso-octane oxidation[END_REF] (N sp = 857, N r = 3598) for the combustion of iso-octane with air would lead to M transp = 6.6 Go which is not affordable. However, the smoothness of the profiles shown in Fig. 3 suggests that ∆T transp can be quite substantially increased without any loss of accuracy.

In the case of reactions rates, the dependence on temperature is exemplified in Fig. 4. For some reactions (see Fig. 5), the reaction rate can exhibit a steep gradient that requires special care, especially for low temperatures. Indeed, even if these temperature values are not encountered in conventional simulations and the chemical species involved in these reactions are generally absent over these temperature ranges, it is nevertheless necessary to ensure an accurate description of profiles. A temperature discretization of ∆T chem = 1 K is then used, and the amount of memory M chem reserved for this tabulation is given by:

M chem = C o × T max -T min ∆T chem × 2N r . (6)
For the chemistry of Curran et al. [START_REF] Curran | A comprehensive modeling study of iso-octane oxidation[END_REF] (N sp = 857 and N r = 3598) a very low value of M chem = 161 Mo is reached.

Finally, the quantification of the computation time can be expressed through the reduced efficiency, RE, defined as

RE = ∆t w N cell × N RK , (7)
with N cell , the number of mesh cells on each processor (N cell = Total number of cells / Total number of processors (N CPU )), and N RK , the step number of the Runge-Kutta temporal scheme. ∆t w is the wall-clock time, i.e. the time between the beginning and the end of one iteration of the solver. The impact of the code optimization based on the tabulation of the reactions rates and the pure species properties as functions of temperature can be evaluated with RE. This strategy leads to substantial performances (See Tab. 3). The computation of reaction rates is 6 to 8 times faster, while the computation time of the transport coefficients, which depends quadratically on the number of species, is reduced by a factor ranging from 14 to 68. These optimizations allow for a significant reduction of the overall calculation time: a factor of 4 is found for the 7)). The percentage of the total time is given in brackets.

HTTC for partially premixed combustion

The formulation of HTTC [START_REF] Ribert | Hybrid transported-tabulated strategy to downsize chemistry for numerical simulation of premixed flames[END_REF] is general for partially premixed combustion. However, this approach is only valid when laminar premixed flames are defined, i.e. for equivalence ratios comprise in laminar flammability limits. Outside such limits, a special treatment is required as for usual tabulation techniques [START_REF] Gicquel | Laminar premixed hydrogen/air counterflow flame simulations using flame prolongation of ildm with differential diffusion[END_REF][START_REF] Van Oijen | Modelling of premixed laminar flames using flamelet-generated manifolds[END_REF]. Usually, flame databases are generated for a small range of equivalence ratio φ between the flammability limits of premixed flames (i.e. φ ∈ [0.6, 1.4] in this work). For cases where chemical reactions take place on a wider range of φ (typically in triple flames), those databases need to be completed or prolonged. The mass fraction of radical species C 2 and C 2 H 5 extracted from 1D premixed methane/air flames is plotted in a (Y c , φ) space in figure 6. For such species, the mass fractions are not equal to zero at the rich boundary (line defined by φ = 1.4), which confirms that the table needs a prolongation if values of φ out of the flammability range of premixed flames are met in the simulation, especially on the rich side.

Most of the radical species fade to zero before the lean boundary is reached, but some of them do not (e.g.

C 2 H 5 ). A prolongation of the table on the lean side will then be needed as well.

Two procedures to get a prolongation of the table out of the flammability limits are proposed here. For the sake of clarity, table prolongations will only be explained for a methane database, with constant fresh gas temperature and pressure (P = 1 bar and T u = 300 K). ratio, the results given by the 1D flame solver do not make any physical sense. Indeed for φ > 2.5, the flame front propagates at velocities smaller than 1 cm/s (Fig. 7(a)), and is very wide, and the 10 meter-long domains used to generate the database are not even sufficient to reach equilibrium values given by EQUIL [START_REF] Lutz | EQUIL: A CHEMKIN implementation of STANJAN for computing chemical equilibria[END_REF] for the progress variable (Fig. 7(b)).

Prolongation of the

The extended HTTC table on φ ∈ [0.5, 5] with the 1D flame solver is displayed in Fig. 8 for C 2 and C 2 H 5 .

For most of the tabulated species, their mass fraction reaches zero for values of the progress variable smaller than the equilibrium value. Moreover, for φ > 2, these species have a negligible mass fraction, whatever the value of Y c , so the part of the flame database for high values of φ is not used. A few radical species are still present at high equivalence ratios with a large mass fraction (e.g. C 2 H 2 ). For these high values of φ, the equilibrium value of Y c computed with REGATH at the end of the computation domain may not be equal to the one given by EQUIL, so a part of the table is missing for these species. Besides, their mass fractions at the rich φ boundary of the table are not zero, which suggests that they should be transported instead of being tabulated. Consequently, when using this extended table, 7 species A better strategy, based on the scaling rules of S2FT 1 , is now presented. The reduction parameters

(C 2 H 2 , C 2 H 3 , C 2 H 4 , C 2 H 5 , C 2 H 6 , CH 2
Y max k (Z), Y A c,k (Z) and Y B c,k (Z), with Z ∈ [Z L , Z R
] are first computed, for every k ∈ m. They are then extrapolated using Gaussian exponential functions. The parameter Y max is then fitted as follows:

             Y max k (Z) = max Z∈[Z L ,Z R ] (Y max k ) exp B L,k (Z -Z max k ) 2 , ∀Z ∈ [0, Z L ] Y max k (Z) = max Z∈[Z L ,Z R ] (Y max k ) exp B R,k (Z -Z max k ) 2 , ∀Z ∈ [Z R , 1] (8)
where Z max k is defined as:

Y max k (Z max k ) = max Z∈[Z L ,Z R ] (Y max k ) (9)
B R,k =                - log Y max k (Z R )/Y max k (Z max k ) (Z R -Z max k ) 2 if Z max k < Z R B L,k if Z max k = Z R (10)
The parameters Y A c and Y B c have a shape more difficult to predict, but most of the time, they increase with Z, because the interval [0, Y c,eq (Z)] widen with Z. They are simply linearly extrapolated, with a slope calculated between Z L and Z R :

Y A c (Z) = Y A c (Z L ) + Y A c (Z R ) -Y A c (Z L ) Z R -Z L (Z -Z L ), ∀Z ∈ [0, Z L ] ∪ [Z R , 1] (11)
Y B c (Z) = Y B c (Z L ) + Y B c (Z R ) -Y B c (Z L ) Z R -Z L (Z -Z L ), ∀Z ∈ [0, Z L ] ∪ [Z R , 1] (12)
As an example, the extrapolated reduction parameters are displayed in figure 9 for species CH. The reconstructed Y k (Y c , Z) maps using the HTTC table extrapolated with this method are displayed for C 2 and C 2 H 5 in figure 10. In the rich extrapolated zone, the Y k (Y c ) profiles are extended beyond the equilibrium point Y c,eq . The impact of this extrapolation error will be assessed in section 5.

The different table extension methods are compared in figure 11. The Y k (Y c , Z) maps are similar for some species, but feature some discrepancies for a few species.

All the extension methods presented here have drawbacks, but the extension of a 1D premixed flame database out of the flammability range imperatively requires the implementation of models.

Automated tool for the generation of tables

A full HTTC simulation requires five main stages (see Fig. 12):

1 With S2FT [START_REF] Ribert | Tabulation of complex chemistry based on self-similar behavior of laminar premixed flames[END_REF], one reduced self-similar profile of species mass fraction, 1. Reading and pre-processing of the 1D premixed flame database;

Y k /Y max k k ∈ m, is stored versus a reduced progress variable Y + c,k = (Y c,k -Y A c,k )/(Y B c,k -Y A c,k ).
2. Processing of the database to get a "raw" table, without self-similarity nor extension;

3. Computation of the reduced profiles and the reduction parameters, and extension of the table ;   4. Writing of the table with the proper format.

Post-processing of the table, to easily check for unexpected behaviors and errors.

For the sake of clarity, the database is supposed to contain one single point in the P and T u directions.

• During the step 1, user-defined parameters are read: definition of Y c , the discretization of Y c and Z directions, choice of XF or YF extension method.

Species names and mass fractions, pressure and temperature are read in the flame database coming from the 1D flame solver. • During the step 2, the full set of species is split to get the tabulated and the transported species, according to the criteria defined in [START_REF] Ribert | Hybrid transported-tabulated strategy to downsize chemistry for numerical simulation of premixed flames[END_REF]:

                             k ∈ m if ∀S ∈ Ω S ,                    Y k (Y c = 0) < ε Y k (Y c = Y c,eq (S )) < ε Y k (Y c = Y c,eq (S )) < ε Y max k (S ) k ∈ M else (13)
where ε = 10 -8 and ε = 0.01. The shape of the mass fraction profiles of each tabulated species is analyzed, for every S ∈ Ω S . The progress variable Y c is then computed using the coefficients provided by the user. The tool checks that Y c increases monotonically in the physical space, or return an error message. A uniform mesh is set for Y c , between zero and the maximum equilibrium value obtained among all S ∈ Ω S . The Y k (Y c ), k ∈ m profiles are interpolated on this new mesh, which is more convenient for the implementation of the self-similarity and the table extension parts, and gives better performances when the table is read during the simulation, than a table with non-uniform indexing. For the same reason, a uniform mesh is set for the mixture fraction Z, from 0 to 1. For each value of Y c , the Y k (Z), k ∈ m profiles are interpolated on this mesh, for Z included in the flammability range.

• The step 3 is optional. If no table extension method is selected, the tabulated mass fractions are set to zero out of the flammability limits. The arrows in Fig. 12 represent the different available options.

• During step 4, the tabulated mass fractions are written in a file, along with a copy of the user-defined input parameters of the generation tool. If the self-similarity has not been requested, the "raw" profiles Y k (Y c , Z)

are stored, else one single reduced profile plus the reduction parameters are stored for each tabulated species.

By default, the stoichiometric reduced profile with P = P min and T u = T min is stored.

Mass fractions of every species are read for each 1D flame in the database. User options are read.

The full set of species is split in a set of tabulated species and a set of transported species.

The shape of the profiles is analyzed to detect anomalies.

The progress variable Y c is computed for every flame, and the maximum of all the equilibrium values Y c,eq max is found. A uniform mesh Y * c is built. The tabulated mass fractions are interpolated on this new mesh.

A uniform mesh Z * is set. The tabulated mass fractions are interpolated on this mesh, in the flammability range.

The mass fraction profiles are extrapolated out of the flammability area, using the XF approach.

The reduction parameters are computed for every mass fraction profile in the flammability area.

The reduction parameters are computed for every mass fraction profile, for Z ∈ [0, 1].

The reduction parameters are extrapolated out of the flammability area, using the YF approach.

The raw table • The step 5 is intended to provide checking and debugging data, such as the maps Y k (Y c , Z) for the whole set of species, or the reduction parameters for the tabulated species.

Y tab k (Z * , Y * c )
For each tabulated species, all the reduced profiles are plotted together to check that they are actually self-similar.

Reference flame

The simulation of edge flames are performed with the solver FTC.

Simulation setup

The simulation setup is a pure methane slot injector surrounded by a co-flow of air. The width of the slot is D f = 2 mm, and the thickness of the injector wall is 0.5 mm. The simulated area along with the boundary conditions is given in Fig. 13. It consists in a two-dimensional domain, 15 millimeters wide and 18 millimeters high, beginning at the injector outlet. The domain is uniformly meshed with 50 micrometer cells.

The detailed kinetic mechanism of Lindstedt [START_REF] Lindstedt | Modeling of the chemical complexities of flames[END_REF] has been used. Species are transported with variable Lewis Symmetry Axis NSCBC Side

Air coflow NSCBC Inlet

T o = 300 K U y,o (x) = U 0 y,o = 0.4 m/s Methane injector NSCBC Inlet U y,f (x) = U 0 y,f " 1 ✓ 2x D f ◆ 2 # T f = 300 K U 0 y,f = 6.0 m/s NSCBC Outlet P = 1 bar U y = U 0 y,o

Results with the FTC solver

The flame tip stabilizes approximately 2 mm above the burner rim, at a radial location where the velocity of the flow is small enough to allow for the flame stabilization. The flame is depicted by few minor species displayed in Fig. 14. Species like H 2 , which is formed in larger quantities at high equivalence ratios, are produced in rich reactive zones and then diffuse in the burnt gases. However some radical species such as CH, are created and directly consumed in the flame. Those radical species are still present downstream of the premixed flame in thin zones, because of the presence ot the diffusion flame. On the lean side of the flow, chemical reactions featuring species containing hydrogen are predominant. The species OH is mainly found downstream of the premixed flame. It is typically a marker of the diffusion flame in edge flames [START_REF] Guo | A numerical study of laminar methane/air triple flames in two-dimensional mixing layers[END_REF][START_REF] Plessing | An experimental and numerical study of a laminar triple flame[END_REF][START_REF] Seo | An experimental study of the fuel dilution effect on the propagation of methane-air tribrachial flames[END_REF].

Indeed, in the present configuration, the flame exhibits a monobrachial structure [START_REF] Chung | Stabilization, propagation and instability of tribrachial triple flames[END_REF]: the lean and rich premixed zones along with a trailing diffusion flame can be observed in Fig. 15(a) using the flame index of Briones et al. [START_REF] Briones | Effects of H 2 enrichment on the propagation characteristics of CH 4 -air triple flames[END_REF], but the premixed wings are merged with the diffusion flame. The maximum of the heat release rate is found located nearby the stoichiometric line, at the triple point, where the three parts of the flame meet (see Fig. 14 and Fig. 15(a)). The contours of heat release rate plotted in Fig. 15 show that the premixed flame front is slanted in relation to the axial direction of the flow. In [START_REF] Kim | Effect of velocity gradient on propagation speed of tribrachial flames in laminar coflow jets[END_REF], Kim et al. have attributed this phenomenon to the effect of the velocity gradient, which is usually strong in jet flames. As a consequence, the propagation velocity of the flame tip cannot be equal to the axial velocity of the flow at the triple point, which is here assumed to be the intersection between the stoichiometric isoline of Z and the isoline Y c = 0.005. In [START_REF] Kim | Effect of velocity gradient on propagation speed of tribrachial flames in laminar coflow jets[END_REF], the propagation velocity is thus assumed to be equal to the velocity U n normal to the premixed front, calculated at the triple point by correcting the axial velocity with the slanted angle of the premixed front. In this work, the velocity U n is calculated from the velocity vector U at the triple point, using the gradient of Y c to determine the normal direction to the premixed front:

U n = ∇Y c |∇Y c | .U (14)
For the present case, U n = 0.209 m/s at the flame tip, while the 1D laminar flame speed at φ = 1 computed with REGATH gives S 0 L = 0.367 m/s. This smaller propagation speed is due to the strong gradient of mixture fraction at the flame tip (Fig. 15(b)).

The reactive points of the simulated domain that contribute to the heat release, i.e. where the heat release rate is at least 1% of its maximum on the whole domain, are plotted in a (Y c , Z) space in figure 16. The flammability limits are here defined as the range of equivalence ratio where reactive points are found (as in [START_REF] Briones | Effects of H 2 enrichment on the propagation characteristics of CH 4 -air triple flames[END_REF]). Even if this definition of the flammability range is different of what is used in experimental studies, it is clear that the limits are extended in the premixed zone, for an equivalence ratio ranging from 0.15 to 2.73. Finally, reactive points are found beyond the equilibrium line extracted from the 1D premixed flame database, when the progress variable is computed as Y c = Y CO + Y CO 2 . Such a phenomena may be due to the presence of the diffusion tail, and to the diffusion of Y c from the rich side to the lean side of the flow.

Edge flame simulations with HTTC

The results of the HTTC simulations obtained with the two prolongation methods (XF and YF) are now compared to the reference FTC results in order to evaluate the capabilities of HTTC to simulate edge flames, and to select the best prolongation method. The two strategies lead to a similar flame structure (Fig. wall, because of a slightly higher propagation speed (U n = 0.209 m/s with FTC, but 0.215 m/s with XF and 0.225 m/s with YF). A very good agreement with the reference temperature and mass fractions, compared on radial profiles at different heights relatively to the flame tip, is obtained (Fig. 18). The agreement is good both in the premixed zone (at height y 0 +1 mm) and in the burnt gases where the diffusion flame shows up (at heights y 0 + 5 mm and y 0 + 10 mm). The method XF is preferred because of its simplicity of implementation.

The simulations performed with HTTC, whatever the methods XF and YF used, are approximately 5 times faster than with the solver FTC. The impact of the coflow of air is provided in Appendix B.

Impact of the progress variable

Until now, the progress variable was simply based on the sum of the mass fractions of the species CO and since the evolution of Y c is monotonic in the physical space, on the whole range of equivalence ratio covered by the table. It has been shown that this definition is well suited to compute 1D premixed flames, where Y c behaves as expected by increasing monoticaly in the physical space, because the same 1D structures are used in the table.

However, the reference simulation of the edge-flame shows here that the validity of this definition can be questioned. Figure 19(a) shows the progress variable Y c along several isolines of Z in the reference simulation, for lean to rich values of the equivalence ratio. For φ = 0.7, 1.0 and 1.2 it decreases in the burnt gases, downstream of the premixed flame front located at a curvilinear abscissa d ≈ 2 mm, whereas it grows monoticaly in the 1D premixed flame database. This phenomenon can be caused by the diffusion of CO and CO 2 toward the outer air flow, caused by the strong gradients of Z and not taken into account in the HTTC tables. For higher values of Z, the other profiles of Y c are monotonic. In other words, in the reference simulation, for a given point (Y c , Z), the mass fraction of the species can have several values, which is not allowed in the framework of tabulated chemistry. This issue is illustrated for the species CH in figure 19(c), especially for φ = 1.2, which is the equivalence ratio where CH is present in the trailing diffusion flame (Fig. 14(b)). This issue highlights the difference between the 2D flame structure met in the reference simulation, and the independent 1D flame structures stored in the HTTC table.

d (mm) (b) Y c = 1 -Y CH 4 -Y O 2 -Y N 2
Y CH + Y c (-) (d) Y c = 1 -Y CH 4 -Y O 2 -Y N 2
To avoid such behavior, a new progress variable, defined as the difference between the total mass and the mass of reactants, is used:

Y c = 1 -Y CH 4 -Y O 2 -Y N 2 . (15)
In Fig. 19(b) Y c now evolves monoticaly along all the mixture fraction isolines in the reference simulation and no value beyond the chemical equilibrium is found. A similar observation is done for Y CH in Fig. 19(d).

This progress variable has been tested in a simulation performed with a self-similar table prolonged with the method XF. The results (see Appendix C for additional profiles) compare well with the same HTTC simulation performed with Y c = Y CO + Y CO 2 and with the FTC solver, even if the height of the flame tip is lower (y 0 = 0.74 mm) with a higher velocity at y 0 (U n = 0.225 m/s). Finally, with Eq. ( 15) the species CO is better described than with Y c = Y CO + Y CO 2 whatever the height considered above the flame tip, contrary to the species CO 2 .

Impact of the mixture fraction

So far, the mixture fraction Z has been computed in the edge-flames by using a combination of all the local species mass fractions of the flow with Bilger's formula [START_REF] Bilger | On reduced mechanisms for methane-air combustion in nonpremixed flames[END_REF]. However, when building the HTTC table, the 1D flames are indexed and located using the mixture fraction in the fresh gases, denoted Z u , and defined as Z u = Z(Y c = 0). In the HTTC simulations presented in this section, the locally computed value Z is assumed to be equal to Z u , and thus is used to access the table. Such an assumption is valid when unity Lewis numbers are used in 1D premixed flames to simplify the diffusion modeling (Z(Y c ) = Z u , ∀Y c ). It is not the case anymore when the differential diffusion effects are taken into account in simulations: the differential diffusion leads to variations of Z within the flame front (Fig. 20), because every species have a different diffusion velocity. In the edge flames presented in this study, strong mixture fraction gradients are met at the flame tip. In this case, the mixture fraction variations due to differential diffusion are small compared to the gradients of Z in the direction normal to the isolines of Z. Consequently, the impact of the differential diffusion in the direction normal to the premixed flame front, on the computation of the value of Z used to access the table, is expected to be small. To test this assumption, a HTTC simulation has been performed, by transporting Z with a Lewis number equal to 1, instead of computing it using the local mass fractions. In the frame of the tabulated thermochemistry method (FPI [START_REF] Gicquel | Laminar premixed hydrogen/air counterflow flame simulations using flame prolongation of ildm with differential diffusion[END_REF] for ex.), solving a transport equation for Z with a unity Lewis number assumption is a usual work-around to limit the effects of the differential diffusion. A self-similar table with a prolongation based on a flame database extended on φ ∈ [0.5, 5] is used (method XF). The progress variable is defined as

Y c = Y CO + Y CO 2 .
The flame structure, the flame stabilization height and the propagation velocity are almost identical, whatever Z is transported with a unity Lewis number, or locally computed with the species mass fractions (y 0 = 0.74 mm and U n = 0.225 m/s). The temperature and mass fractions are also very similar at every heights in the radial profiles, except for the minor species (OH and H 2 ), which are slightly underestimated when Z is transported (see Appendix C). The outline of a method is proposed in Appendix D, to get the proper value of Z to access the HTTC table, by reconstructing the value of Z u corresponding to the local Z.

Self-similarity of radical species in triple flames

The HTTC simulations presented in this study leads to promising results by making use of tables based on 1D premixed flames. This study also suggests that an even better agreement with the reference simulations could be reached if the diffusion effects, due to the strong gradients of Z, are taken into account in the tables, as in [START_REF] Nguyen | Multidimensional flamelet-generated manifolds for partially premixed combustion[END_REF]. Indeed, the mass fraction contours of some tabulated species spread on a larger range of Z in the actual edge flame than in the 1D lookup table. Another strategy is the tabulation of 2D profiles in the frame of S2FT [START_REF] Ribert | Tabulation of complex chemistry based on self-similar behavior of laminar premixed flames[END_REF][START_REF] Wang | Self-similar behavior and chemistry tabulation of burnt-gas diluted premixed flamelets including heat-loss[END_REF][START_REF] Ribert | A multi-zone self-similar chemistry tabulation with application to auto-ignition including cool-flames effects[END_REF], i.e. using the self-similarity property of hydrocarbon flames. Accordingly, reduced mass fraction profiles Y + (Y + c ) have been extracted from the FTC reference simulation along several isolines of mixture fraction, on a range of equivalence ratio (φ ∈ [0.2, 2.7]). The reduced profiles are plotted in figure 22 for a few radical species (CH, CH 2 , C 2 H 3 and C 2 H 5 ), which must be stored in HTTC tables.

The progress variable is computed as

Y c = 1 -Y CH 4 -Y O 2 -Y N 2 .
For all the tabulated species, the reduced profiles are correctly self-similar, even if both a premixed and a diffusion flame structure are met along each isoline. For a few values of the mixture fraction, the mass fraction does not go back to zero in the burnt gases and the reduced profile is not perfectly superimposed with the other profiles. Those values are located around the stoichiometry, where the diffusion flame takes place downstream of the premixed tip region. The reduced profiles extracted from the 2D edge-flame simulation are in a fair agreement with the stoichiometric reduced profiles coming from the 1D flames, showing that the self-similarity is a very generic property which can be found in various configurations of hydrocarbon flames.

Conclusion

The hybrid transported-tabulated chemistry (HTTC) method has been implemented into the reacting fully explicit solver SiTCom-B for DNS/LES. A detailed description is given as well as a strategy to optimize the computation of reaction rates and species properties that are essential when a detailed kinetic is used for flame computation. The original version of HTTC has been extended to partially premixed combustion by adding a prolongation of the self-similar flame tabulation (S2FT) outside the flammability limits. The two methods, either based on extra flames for the database (prolongation XF) or on the selection of an extrapolation function (gaussian, linear, etc.) for the species mass fractions or S2FT parameters (prolongation YF), offer similar results. However, the method XF has the advantage of being straightforward, because the table is built up without using any specific extrapolation procedure (YF). For partially premixed combustion, HTTC simulations cannot be performed using tables without prolongations, because a part of the heat release is missed, especially in the rich premixed region, and because it generates discontinuities that endanger the stability of the solver.

With HTTC the kinetic scheme is left unchanged and only major species are transported with the flow whereas minor species are tabulated from 1D premixed flames to save CPU time. A methane-air edge flame featuring strong mixture fraction gradients has been simulated with success by comparison with the reference flame computed with the fully transported chemistry (FTC). The shape, the propagation speed and the stabilization height of the flames have been predicted with a correct accuracy, with a computational cost divided by approximately 5 compared to the FTC solver, thanks to the increase of the chemical time step. The temperature, the mass fraction of the main products of the combustion and of some minor species have also been computed with a good agreement with the reference results. The differential diffusion in the direction normal to the isolines of Y c does not seem to play a critical role in the calculation of the value of Z used to access the HTTC table, in the edge flames simulated in this paper, since the method used to compute Z does not affect the results. Moreover, solving a transport equation for Z is against the principle of HTTC, and the unity Lewis number assumption should ideally be avoided. However, in cases where the gradients of Z are smaller (e.g. a fully premixed flame), the differential diffusion may become more effective and may affect the local building of Z. Thus, large errors may show up if Z is computed locally, as observed in 1D flames.

Future work should focus on turbulence-flame interaction with HTTC. Since HTTC maintains a part of tabulated chemistry, one may expect its accuracy to be diminished for flames located in the broken reactions regime [START_REF] Borghi | Turbulent combustion modelling[END_REF]. used, but a simple average using the two left and right neighbors is used:

A i * L,k = A i-2 L,k + A i-1 L,k + A i L,k + A i+1 L,k + A i+2 L,k 5 (A.5)
where A i * L,k is the smoothed value of A L,k at the discrete value of Y c with the index i. The extrapolated table with this first sub-method is plotted in figure A.24, for C 2 and C 2 H 5 . Unlike the XF prolongation method, the choice of the transported species is not modified by the prolongation process. Note that, so far, the tables extrapolated with such method have been built without making use of the self-similarity. Enabling the self-similarity to build up the table led to the failure of the simulation. Indeed, when the flame database is extrapolated in the Z direction, the mass fraction profiles Y k (Y c ) artificially built for φ < 0.6 and φ > 1.4

are not self-similar with the "real" profiles extracted from the flame database generated for φ ∈ [0.6, 1.4].

Since raw tables built without self-similarity are out of the framework of HTTC for memory space reasons, the YF prolongation was not further exploited.

) 1 k

 1 Thermodynamic database reading (thermal capacities, enthalpies, diffusion coefficients, etc) Kinetic mechanism reading Mesh reading and primitive variable fields (P, T, U, Y k ) reading from initial solution Initialization of conservative variables ρE, ρU et ρ k from primitive variable fields -At the beginning of every iteration: time step computation ∆t Computation of the convective time step ∆t C and of the diffusive time steps ∆t D ρE and ∆t D ρU Computation of the diffusive time steps ∆t D ρ k and of the chemical time step ∆t chem using the whole set of species Computation of the diffusive time steps ∆t D ρ k and of the chemical time step ∆t chem of convective and diffusive fluxes of ρE and ρU Computation of the convective and diffusive fluxes of the partial densities ρ k , for the whole set of species Computation of the convective and diffusive fluxes of the partial densities ρ k , for transported species only Computation of the species source terms, using every species and reaction of the full kinetic mechanism Update of the conservative variables During every Runge-Kutta step: computation of all the variables of the flow Computation of the mixture density: ρ = N sp k=1 ρ n k Computation of the mixture density: ρ = k∈M ρ n k + k∈m ρ n-Computation of every mass fraction : Y k = ρ k /ρ, k ∈ [1, N sp ] Computation of transported species mass fractions: Y k = ρ k /ρ, k ∈ mixture fraction using the whole set of species: Z = Z(Y)

(

  ∆t CFL ), diffusive (∆t D ) and chemical (∆t Chem ) time steps, thus ensuring the stability of the explicit Runge-Kutta algorithm. ∆t CFL = CFL × min [∆ i / (|u i | + a)], with a, the local speed of sound, ∆ i , the size of the mesh cell in the direction i, and CFL = 1.0. ∆t D = min [∆t D Φ ] with ∆t D Φ = Fo × min [∆ 2 i / D Φ ]. Φ ≡ ρE, ρU or ρ k , D Φ is a diffusion coefficient and Fo = 0.5. ∆t Chem = min [(a tol + r tol × ρY k ) / | ωk |] with a tol = 1.0 × 10 -12 kg.m 3

Figure 2 :

 2 Figure 2: One-dimensional CH 4 -air flame at φ = 1 (T u = 300 K and P u = 1 bar). Squares : reference computed with REGATH.Lines : results from SiTCom-B using different solvers (black solid : FTC, blue solid : HTTC, blue dashed : fast HTTC).

  Coefficients Φ jk (Eq. (4)). Only the two possible cases are plotted for clarity.

Figure 3 :

 3 Figure 3: Diffusion coefficients of each species k of the Lindstedt kinetic scheme (all species shown) as a function of temperature. For each species, coefficients have been normalized by their maximum value.

Figure 4 :

 4 Figure 4: Forward (k f ) and reverse (k r ) reactions rates of elementary reactions from Lindstedt's mechanism as a function of temperature.

  Zoom on range T ∈ [200, 220] K.

Figure 5 :

 5 Figure 5: Reverse (k r ) reactions rates of reaction CH 2 = C + H 2 from Lindstedt's mechanism as a function of temperature.

4 Y 5 Figure 6 :

 456 Figure 6: Species mass fraction of C 2 and C 2 H 5 in the original 1D methane flame database (φ ∈ [0.6, 1.4]), as a function of the progress variable and the equivalence ratio.

  Equilibrium value from EQUIL ( ) and RE-GATH: 10 km-long () and 1000 km-long ( ) mesh.

Figure 7 :

 7 Figure 7: Characteristics of laminar 1D methane/air flames computed with REGATH.

  c,eq (φ)(b) C 2 H 5

Figure 8 :

 8 Figure 8: Species mass fraction of C 2 and C 2 H 5 in the extended 1D methane/air flame database (φ ∈ [0.5, 5.0]), as a function of the progress variable and the equivalence ratio.

3 . 1 . 2 .

 312 CO, CH 3 ) are added to the set M of the transported species (N M = 20), but without modification of the chemical time step. Extrapolation of the 1D premixed flames The second prolongation method (N M = 20), and denoted YF prolongation, is an extrapolation of the database Y k (Y c , φ), k ∈ m, φ ∈ [0.6, 1.4] out of the flammability range. Then, when φ ∈ [0.6, 1.4], the mass fractions provided by the 1D premixed flame database are still used, but when φ < 0.6 or φ > 1.4, an extrapolation model is used to fill the missing parts of the table. A first attempt is proposed in Appendix A based on the shape of Y k (Z) curves. Although the formulation leads to comparable results than the other two procedures, its use in the edge-flame simulation failed, blowing out the code.

Figure 9 : 5 Figure 10 :

 9510 Figure 9: Extrapolation of the reduction parameters for species CH. Solid lines: values computed using the flame database on the flammability range. Dashed lines: extrapolated values. Black: Y max , blue: Y A c , red: Y B c . Zone 1: lean extrapolated zone. Zone 2: flammability range. Zone 3: rich extrapolated zone.

5 Figure 11 :

 511 Figure 11: Contours of species mass fractions in 1D CH 4 /air flames. Prolongation XF ( ), and YF ( ), ( , not used, see Appendix A).

Figure 13 :

 13 Figure 13: Dimensions and boundary conditions of the simulation domain.

( a )Figure 14 :

 a14 Figure 14: Edge flame simulated with the FTC solver. In green, isolines of mixture fraction corresponding to φ = 0.6, 1.0, 1.4 and 5.0. In pink, isolines of heat release rate: 0.05, 0.2, 0.5, 1.0, 1.5 and 2.0 GW/m 3 .

Figure 15 :

 15 Figure15: Flame index[START_REF] Briones | Effects of H 2 enrichment on the propagation characteristics of CH 4 -air triple flames[END_REF] and norm of mixture fraction gradient. Only cells where the heat release rate is larger than 1% of its maximum over the domain are displayed. In black, isolines of mixture fraction corresponding to φ = 0.6, 1.0, 1.4 and 5.0.In white, isolines of heat release rate: 0.05, 0.2, 0.5, 1.0, 1.5 and 2.0 GW/m 3 .

Figure 16 :

 16 Figure 16: Scatter plot in a (Y c , Z) space of the points of the FTC simulation with a heat release rate larger than 1% of its maximum over the whole domain, colored by the flame index. The red line is the equilibrium extracted from the 1D flame database.

Figure 17 :

 17 Figure 17: Isolines of stoichiometric mixture fraction (green) and of progress variable (Y c = 0.005, black). FTC (-) and HTTC (---) simulations. Prolongation XF (left) and YF (right).

CO 2 .Y m at y 0 + 10 mmFigure 18 :

 21018 Figure 18: Radial profiles of temperature (black), major products (M ≡ CO in green, CO 2 in red and H 2 O×1.2 in blue) and minor species (m ≡ OH in purple and H 2 in orange) at different heights above y 0 for FTC ( ) and HTTC (XF -and YF -• -).

  Y c = Y CO + Y CO 2

  Y c = Y CO + Y CO 2

Figure 19 :

 19 Figure 19: Impact of the definition of the progress variable Y c vs. a curvilinear abscissa d (a,b) and for CH mass fraction (c,d), along mixture fraction isolines in the reference flame: φ = 0.7 ( ), 1.0 ( ), 1.2 ( ) and 1.8 ( ).

Figure 20 :Figure 21 :

 2021 Figure 20: Mixture fraction Z versus the progress variable Y c , in methane flames for several φ ∈ [0.6, 1.4] (P = 1 bar and T u = 300 K), computed with the 1D flame solver REGATH.

5 Figure 22 :

 522 Figure 22: Reduced mass fraction profiles of radical species extracted from the reference simulation (black), plotted along isolines of mixture fraction, corresponding to several value of φ in [0.2, 2.7]. The co-flow velocity is 0.4 m/s. The reduced profiles extracted from the stoichiometric 1D flame (P = 1 bar and T u = 300 K) are plotted in red.

2 Figure A. 23 :

 223 Figure A.23: Extrapolation of Y(Z) profiles.

5 Figure A. 24 :

 524 Figure A.24: Species mass fraction of a few radical species in the extrapolated 1D methane/air flame database, as a function of the progress variable and the equivalence ratio. The Y(Z) profiles are extrapolated in areas where φ < 0.6 and φ > 1.4, using extrapolations of Y(Z) profiles.

Figure B. 27 :Figure B. 28 :Figure C. 29 :Figure D. 30 :Figure D. 31 :

 2728293031 Figure B.26: Isolines of mixture fraction isolines (in green, corresponding from right to left to φ = 0.6, 1.0, 1.4 and 5.0), heat release rate (in pink, 0.05, 0.2, 0.5, 1.0, 1.5 and 2.0 GW/m 3 ) and progress variable Y c = Y CO + Y CO 2 (in black, 0.005 and 0.15), in the edge flames simulated with the FTC solver.

Table 1 :

 1

Changes introduced by the HTTC solver in the simulation process, compared to the FTC solver. The steps shared by both solvers are written in black, while the different steps are colored. A set of conditions in the fresh gases is denoted S , S ≡ (Z, P, T ).

  Lindstedt mechanism for methane, 7 for the Luche et al. mechanism for kerosene, and 43 for the mechanism for iso-octane.

				Mechanism
	Task	Methane [2]	Decane [3]	Iso-octane [41]
		29 species	91 species	857 species
		141 reversible Rns	991 elementary Rns	3598 reversible Rns
	Calculation of diffusion coefficients
	Non-Optimized	11.86	(33.9%) 128.44 (60.9%)	20800.41 (96.7%)
	Optimized	0.81	(8.8%)	4.60	(16.6%)	302.02	(61.2%)
	Ratio non-opt./opt.	14.64		27.92		68.87
	Calculation of thermodynamic variables (T , c v , ...)
	Non-Optimized	4.24	(12.1%) 20.72	(9.8%)	234.80	(1.1%)
	Optimized	1.45	(15.9%) 4.64	(16.7%)	52.22	(10.6%)
	Ratio non-opt./opt.	2.92		4.47		4.50
	Calculation of source terms
	Non-Optimized	14.19	(40.6%) 52.06	(24.7%)	382.26	(1.8%)
	Optimized	2.19	(23.8%) 8.85	(32.0%)	45.40	(9.2%)
	Ratio non-opt./opt.	6.48		5.88		8.42
	Total
	Non-Optimized	34.97		210.64		21511.16
	Optimized	9.13		27.71		493.33
	Ratio non-opt./opt.	3.83		7.60		43.60

Table 3 :

 3

Distribution of the computation time for an iteration of the detailed chemistry solver, expressed with the reduced efficiency RE in µs (Eq. (

  is written in a file.

	1 st step	Data	reading
	2 nd step	Raw table building
	3 rd step	Extrapolation and	self-similarity
	4 th step	Table writing	The reduced table, made up of a single reduced profile Y + k (Y + c ) for each species k and the parameters Y A c,k , Y B c,k and Y max k , is written in a file.
	5 th step	Post-processing	For every species, the maps Y tab k (Z * , Y * c ) are displayed.	For every species, all the reduced profiles are displayed together to check their self-similarity property. The maps Y tab k (Z * , Y * c ) are displayed.

Figure 12: Description of the automated table generation tool, designed for the HTTC solver of SiTCom-B.

Acknowledgments

Bastien Duboc was financially supported by the Region Haute-Normandie. Computations were performed using HPC resources from CRIANN. The pre-processing code to build up the HTTC table is available upon request at https://www.coria-cfd.fr/index.php/HTTC.

Since the profiles Y k (Z) are not symmetrical, two coefficients A L and A R are used on the lean and the rich side, respectively, to get a well matching fitting function:

where

Fitted profiles are displayed in figure A. [START_REF] Pons | Mass transfer and combustion in transcritical non-premixed counterflows[END_REF]. If Y L = 0 or Y R = 0, the extrapolation process is not needed on the lean side or the rich side, respectively, and the fitting function is set to zero. For a few species (e.g. [START_REF] Pons | Mass transfer and combustion in transcritical non-premixed counterflows[END_REF], and for some Y c values, Z max k (Y c ) = Z R and (∂Y k /∂Z)(Z R ) 0, so it is impossible to predict the location and the value of the maximum of such Y k (Z) profiles, since they monoticaly increase along Z. In this case, the exponential fit may not be well adapted, but is still used, to keep the extrapolation process as simple as possible. The fitting coefficients are calculated on the lean side (i.e. for Z < Z L ), and symmetrically duplicated on the rich side.

With this 1D extrapolation approach, Y k (Z) profiles are independently treated for each value of Y c . It means that, for a given value of Z, reconstructed Y k (Y c ) profiles in the extrapolated zone may not be continuous, which is expected to be harmful for the table generator, especially during the self-similarity shown that, for strong jet conditions (velocity of the jet much larger than U coflow ), a change in the velocity of the co-flow has a much more appreciable impact on the stabilization height of the flame than a change in the velocity of the fuel jet. In the simulations presented here, strong jet conditions are encountered, since the velocity of the fuel jet is approximately ten times larger than the velocity of the co-flow of air (U 0 y, f = 6 m/s vs. U coflow = 0.4 or 0.6 m/s). This explains why the stabilization is done higher with U coflow = 0.6 m/s, in spite of a larger propagation speed.

The results obtained with HTTC (XF prolongation, Y c = Y CO + Y CO 2 ) give a flame shape similar to the one obtained with the FTC solver (y 0 = 7.76 mm and U n = 0.292 m/s) but with a shift toward the injector: y 0 = 4.67 mm and U n = 0.303 m/s if Z is locally computed with Bilger's formula [START_REF] Bilger | On reduced mechanisms for methane-air combustion in nonpremixed flames[END_REF] or y 0 = 4.84 mm and U n = 0.304 m/s if Z is transported with the flow. The mixture fraction gradient is still large, so that the influence of the differential diffusion along isolines of the progress variable may be small. However, the flames simulated with HTTC propagate faster and stabilize at a lower height than the reference flame (Fig.

B.27).

The agreement with the reference on the temperature and mass fractions measured in the radial profiles is

Appendix D. Correction of the local value of the mixture fraction

It has been shown in section 5 that, because of the differential diffusion effects, the value of the mixture fraction Z varies accross the flame, and therefore differs from its value Z u in the fresh and the burnt gases (Fig. 20). In the frame of HTTC, the local value of Z used as a control parameter for the table and computed from the species mass fractions (e.g. using Bilger's definition [START_REF] Bilger | On reduced mechanisms for methane-air combustion in nonpremixed flames[END_REF]), shows deviations with the value of Z u that should be actually used to access the table, the latter being defined by Y m (Y c , Z u ). No significant effect of this phenomenon has been observed in the edge flames presented in section 5, because of the strong mixture fraction gradients met in this configuration. However, it creates a feedback loop in 1D premixed flames when Z is computed locally, where the wrong value of Z Z u used to access the table leads to inconsistent values of the mass fractions of the tabulated radical species, which in turn leads to even more erroneous values of Z. Mass conservation is not ensured and wrong results are obtained (Fig. 21). The same issue is also expected to happen in flames with low stratification.

The outline of a method to reconstruct locally the proper value of the unburnt mass fraction Z u to access the table is briefly presented here. For the sake of clarity, the pressure P and the fresh gas temperature T u are here supposed constant, but could be added as dimensions of the table. Since the local value of Z based on Bilger's definition cannot be used directly, it is here proposed to use data obtained from the transported

The variable Y * must be built such as there is a unique value of Z u (Y c , Y * ) for a given set (Y c , Y * ). Such a condition can be met by setting Y * equal to the product of linear combinations of the transported species:

where the β k and γ k are real numbers. Several combinations have been tested. For some of them, the uniqueness is not ensured. For example, with equal (according to a given threshold), and the final value of Z u is used to read the tabulated species mass fractions. The simulation of a methane-air flame showed that around 5 iterations are needed to converge to the proper value of Z u . The method using the self-similarity of the Y * profiles led to results similar to those obtained without self-similarity, but the number of data to store is greatly reduced.