Gaël Boury

Yves D'angelo

On third order density contrast expansion of the evolution equation for wrinkled unsteady premixed flames

The dynamics of flat-on-average wrinkled flame front propagating through gaseous premixtures is considered. Leading the asymptotic expansions in powers of the burnt to unburned fractional density contrast (0 < γ < 1) to third order, an evolution equation (called S3) is obtained for the instantaneous front shapes. It reduces to Sivashinsky's original equation (called S1) as γ -→ 0. It also modifies a previous attempt by Sivashinsky & Clavin (called S2) to improve it. Numerical integrations of the S3 equation reveals that the new quadratic and cubic nonlinearities featured at 3rd order happen to mutually compensate partially one another for realistic γ's, and are negligible at γ 1. As a result, the flame shape and speed solutions to S3 nearly coincide with those of a S1/S2 type of equation, even for a 10-fold density variation (γ = 0.9) and for unsteady situations, provided a single O(1) coefficient a(γ) be adjusted therein, once for all for each γ. The O(γ 2 ) (and small) correction to it mainly originates from a quartic nonlinearity of geometrical origin. The agreement carries over to comparisons with some DNS of 2D steady wrinkled fronts. A phenomenological (yet asymptotically correct at γ 1 and exact in the linear limit) interpolating model equation is finally proposed to try and account for inertia effects associated with fast transients (e.g. acoustics related) while reproducing the above results on steady patterns.

Introduction

For all its connections with applications -e.g. turbulent combustion, internal combustion engines, gas turbines, industrial burners ... -the problem of wrinkled flame propagation through premixed gases is a central one in combustion science, which is by essence non-linear. Involving elliptic (for the velocity potentials) or hyperbolic (for vorticity) equations and a free boundary (the flame itself), it is further hampered by vorticity creation across the latter and by the apparent impossibility to analytically solve the Euler or Navier-Stokes equations for generic non potential flows. Direct numerical simulations (DNS) of reactive Navier-Stokes equations may, in principle, give access to "exact results" (i.e. with no large scale modeling) on flame topology and dynamics. However, since these flames are very thin (∼ 0.1 mm in usual conditions), real-scale DNS need very high spatio-temporal resolution and are extremely demanding in terms of computational resources : they are sometimes involving tens of million of CPU hours performed on hundreds of thousand of processors (see for instance recent communications [START_REF] Chen | Petascale direct numerical simulation and modelling of turbulent combustion[END_REF] and [START_REF] Frouzakis | DNS of turbulent autoignition of hydrogen using a spectral element solver for low Mach number combustion[END_REF]). Moreover, due to the very large amount of data to be post-processed, DNS do not necessarily allow for simple physical analysis.

On the other hand, one can take advantage of this scale separation and consider the flame as an infinitely thin front, separating the fresh (or unburned, referred to with u subscript in the sequel) mixture from the burnt (b subscript) gas. Hence, only the flame surface or curve needs be parameterized -one spatial dimension being removed. Also, many physical parameters can be lumped into few ones. The Evolution Equation Modeling approach, consisting in finding an equation for the flame surface dynamics only -and not solving for the whole reactive flow -does not (nor claims to) replace the full 3D reactive equations. It is usually limited to simple geometrical configurations (plane, cylindrical or spherical on average). To date, no exact evolution equation is available. However, if sufficiently precise equations can be derived or built in different contexts (slow or fast transients, expanding flames, strained flames, gravity effects, acoustics...), they may provide pertinent information on flame dynamics or even be used as building blocks of (larger scale) sub grid scale modeling, e.g. in LES (Large Eddy Simulations) of reactive flows.

In a seminal work [START_REF] Sivashinsky | Nonlinear analysis of hydrodynamic instability in laminar flames. Part 1 : derivation of basic equations[END_REF], Sivashinsky realized that the unburned (ρ u ) to burnt (ρ b < ρ u ) density contrast 0 < γ ≡ (ρ uρ b )/ρ u < 1 may be used as a viable bifurcation parameter, albeit of a special kind since all flames have γ > 0 : if γ 1, the Landau-Darrieus instability mechanism of spontaneous wrinkling is weak, enabling the author of [START_REF] Sivashinsky | Nonlinear analysis of hydrodynamic instability in laminar flames. Part 1 : derivation of basic equations[END_REF] to derive a leading-order, weakly non-linear equation, called S1 here, and also known as Michelson-Sivashinsky equation (not to confuse it with the Kuramoto-Sivashinsky equation) for the instantaneous flame shape. Numerics [START_REF] Michelson | Nonlinear analysis of hydrodynamic instability in laminar flames. Part 2 : numerical experiments[END_REF] then analysis [START_REF] Thual | Application of pole decomposition to an equation governing the dynamics of wrinkled flame fronts[END_REF] revealed that the S1 equation describes wrinkled flames qualitatively well. An attempt to go to next (second) order in the γ expansion [START_REF] Sivashinsky | On the non linear theory of hydrodynamic instability in flames[END_REF] was partly successful. Here, we correct this "S2" equation, then go to third order following the same perturbative approach. The solutions to the evolution equation (called S3) so obtained are studied numerically and in fact exhibit striking resemblance with those of an equation (called S-fit) that has the same structure as S1 (and S2, actually), provided a single γ dependent coefficient, featured in S3, be slightly modified therein. Such a resemblance carries over to comparisons about unsteady fronts and with direct simulations of steady fronts, even for realistic γ's. A phenomenological way of extending all this to fast flameshape transients is finally proposed; it can also account for gravity effects, time-dependent or not.

The paper is organized as follows. The model and propagation are introduced in section 2, while the coordinates and expansion scheme are presented in section 3. The evolution equations S1-S3 are derived in section 4 (the most technical part of it being summarized in the Appendix).

Section 5 compares the solutions to S1-S3 among themselves and with others. A model for fast transients is proposed in section 6. We end up with concluding remarks and open questions.

Model

At current time t, the flame is considered here to be a curve x = F (y, t) in the fixed cartesian frame (x, y) defined in figure 1, and separates two twodimensional incompressible flowfields u(x, y, t) where density ρ is either ρ u (upstream, x < F ) or ρ b = (1γ)ρ u < ρ u (downstream). Euler equations are assumed to govern the velocity u = (u, v) and pressure p on both sides. Rankine-Hugoniot relationships, specialized to vanishingly-small Mach numbers [START_REF] Williams | Combustion Theory[END_REF], are meant to hold across the line x = F (y, t). For simplicity, the Markstein [START_REF] Markstein | Non steady flame propagations[END_REF] local propagation law is postulated

n.(u -D) x=F -= S L (1 - γ 2k n F yy /(1 + F 2 y ) 3/2 ), (1)
where n ≡ (1, -F y )/(1 + F 2 y ) 1/2 is the local unit normal to the front, D ≡ (F t , 0), the subscripts y or t represent partial differentations (e.g. (.) y ≡ ∂(.)/∂y) and S L > 0 is the prescribed propagation speed of a flat flame. The prescribed k n > 0, to be later identified with a neutral wavenumber, is related to the effective Markstein length (L) by Lk n = γ/2, that is the only local length scale of the problem. Once endowed with appropriate lateral boundary conditions, here taken to be L box ≡ 2π/k box -periodicity along the y direction for some k box < k n , the Euler equations and Hugoniot relations are in principle enough to compute u x=F -in terms of a presumed -and smooth enough -F (y, t).

Then, equation (1) explicited as

F t + S L 1 + F 2 y -1 + F y v| x=F -= u| x=F --S L + S L γ 2k n F yy 1 + F 2 y , (2)
should provide an evolution equation for F (y, t) itself. Note that the speed S T at which the front advances on average towards the uniform fresh mixture (x = -∞) is given -through overall continuity argument -by

1 ≤ S T (t)/S L = 1 + F 2 y 1 + F 2 y /2(1 -F 2 y /4 + ...) (3)
because the curvature term in (1) vanishes upon averaging along the transverse coordinate y : the operation denoted . is defined as

f (y) L box ≡ L box 0 f (y)dy. (4)
Thus, S T /S L -1 is the fractional increase (per unit y) in flame length caused by wrinkling. Interestingly, equation (3) holds whenever (1) is replaced by n.(u -D) x=F -= S L (1q y (y)/(1 + F 2 y ) 1/2 ), for any q(y). In particular, equation ( 3) is valid for the steady patterns considered in [START_REF] Bychkov | On the dynamics of a curved deflagration front[END_REF], where q y (y) also accounts for stretch. As evoked later on in section (5.4), relation [START_REF] Sivashinsky | Nonlinear analysis of hydrodynamic instability in laminar flames. Part 1 : derivation of basic equations[END_REF] applies to in even more general situations, including those considered in the DNS of the problem at hand.

Curved coordinates and small γ expansions

At this stage, it is convenient to change the streamwise coordinate from x to X = x -F (y, t), then to introduce the scalings summarized in table I.

Fresh side (X < 0) Burnt side (X > 0)

k n X = ξ u = S L (1 + γ 2 U -) u = S L 1 1-γ + γ 2 U + k n y = η v = S L γ 2 V - v = S L γ 2 V + γ S L k n t = τ p = p u + ρ u S 2 L γ 2 P -p = p b + ρ u S 2 L γ 2 P + k n F = γ φ p b -p u = -ρ u S 2 L γ 1-γ
Table I: The various nondimensional variables (ξ, η, τ ) and unknowns (U ± , V ± , P ± , φ); the latter are O(1) but all vanish for a steady, flat flame.

Such scalings, where (U ± , V ± , P ± , ξ, η, τ) are meant to be O(1) as γ -→ 0, were dictated by the need of balancing F t /S L , u| x=F -/S L -1 (estimated from the exact linearized dynamics [START_REF] Darrieus | Propagation d'un front de flamme[END_REF][START_REF] Landau | On the theory of slow combustion[END_REF] at γ 1, see eq. ( 25)), γF yy /k n , and the leading order nonlinearity 1 + F 2 y -1 ∼ F 2 y /2, all featured in (2). As for the estimates on v and p, they follow from the continuity equation and the momentum balances, respectively. The scaling of X is provided by the potential nature (i.e. by Kelvin's theorem) of the upstream flowfield, which cannot introduce any intrinsic length scale. Symbolically, Table 1 assumes

1 S L F t ∼ F y 2 ∼ γF k n y 2 ∼ γ F y , (5)
to be compared with the evolution equations obtained at various orders in γ.

The Euler equations to be solved for ξ < 0 and ξ > 0 then read (with indices ± suppressed)

U ξ + V η = γφ η V ξ ( 6
)
γU τ + 1 1 -γ + γ 2 U -γ 2 φ τ U ξ + γ 2 V U η -γ 3 φ η V U ξ = - 1 1 -γ P ξ (7) γV τ + 1 1 -γ + γ 2 U -γ 2 φ τ V ξ + γ 2 V V η -γ 3 φ η V V ξ = - 1 1 -γ P η + 1 1 -γ γφ η P ξ (8)
where the underlined quantities (representing measures of ρ u /ρ) need be replaced by 1 on the unburned side (ξ < 0). Next, combining the local propagation law (2) and the Hugoniot relations allows one to relate the various jumps [U] (e.g.

[U] ≡ U + (0 + ) -U -(0 -)), [P ],
[V ] across the front (ξ = 0) to the flame shape φ(η, τ ) itself, as

[U] = 1 γ(1 -γ) ⎛ ⎝ 1 1 + γ 2 φ 2 η ⎛ ⎝ 1 - γ 2 2 φ ηη 1 + γ 2 φ 2 η 3 2 ⎞ ⎠ -1 ⎞ ⎠ = O(γ) ,( 9
) [V ] = - φ η 1 -γ 1 1 + γ 2 φ 2 η ⎛ ⎝ 1 - γ 2 2 φ ηη 1 + γ 2 φ 2 η 3 2 ⎞ ⎠ = O(1) , (10)
[P ] = 1 γ(1 -γ) ⎛ ⎝ 1 - ⎛ ⎝ 1 - γ 2 2 φ ηη 1 + γ 2 φ 2 η 3 2 ⎞ ⎠ 2 ⎞ ⎠ = O(γ) , (11)
In terms of the same variables, the propagation law (2) is

γ 2 φ τ + 1 + γ 2 φ 2 η -1 + γ 3 φ η V | ξ=0 -= γ 2 U| ξ=0 -+ γ 2 φ ηη 1 + γ 2 φ 2 η (12)
All the unknowns (U ± , V ± , P ± , φ) are to be expanded in power series of γ -→ 0, e.g. φ = φ (0) +γφ (1) +... . The O(γ) rhs of (9) could also be expanded as 0 + γ[U] (1)

+ γ 2 [U] (2) + ..

., and similarly for [V ] = O(1), [P ] = O(γ).

However, doing so would render the algebra very tedious when carried out to high orders in γ. The following remarks help simplify it. The Laplace problems encountered at each order (see below) being linear, and because

[V ] [U] ∼ [P ], it is expedient to impose [U (0) ] = [P (0) ] = 0 and [V (0) ] = [V ] at leading order, where [V ] is the complete RHS of (11), instead of [V (0) ].
At order γ, one may impose [V (1) ] = 0, [U (1) ] = [U], [P (1) ] = [P ]; and at orders γ i , i ≥ 2, one is left with [U (i) ] = [V (i) ] = [P (i) ] = 0. This shortens the calculations, and avoids the burden of sorting then gathering the various contributions to U(0 -), V (0 -) issued from expanded forms of [START_REF] Bychkov | On the dynamics of a curved deflagration front[END_REF][START_REF] Darrieus | Propagation d'un front de flamme[END_REF][START_REF] Landau | On the theory of slow combustion[END_REF][START_REF] Frankel | An equation of surface dynamics modeling flame fronts as density discontinuities in potential flows[END_REF]. The latter are only needed when writing (2) or ( 12) at various orders, in a consistent way. The above procedure could be useful if more general propagation laws than (1) were made use of.

S1, S2, S3 equations

A key point

At each order in γ, we shall have to solve the Laplace equation Z ξξ +Z ηη = 0 for some velocity potential Z(ξ, η, τ ), with presumed values [Z ξ ] and [Z η ] for the jumps of its space derivatives between ξ = 0 -and ξ = 0 + . To this end, we first split Z as Z + Z , whereby Z has a zero transverse average ( Z = 0); Z ξ is a function of time only, as a consequence of Z ξξ ≡ 0. Next, we Fourier-transform (ξ ↔ k, Z ↔ Z) the equation for Z to get Z = Z(0 ± ) exp(∓|k|ξ) as acceptable solutions; these are bounded since Z has no k = 0 component ( Z = 0). After the inverse Fourier transform is taken, all this ultimately yields

Z ξ (0 ± ) = 1 2 H ([Z η ], η) ± 1 2 [Z ξ ] + Z ξ (13)
Z η (0 ± ) = ± 1 2 [Z η ] - 1 2 H ([Z ξ ], η) (14)
Z ξξ (0 ± ) = ∓ 1 2 ([Z η ]) η + 1 2 H ([Z ξ ] η , η) = -Z ηη (0 ± ) ( 15
)
Z ξη (0 ± ) = - 1 2 H ([Z η ] η , η) ± 1 2 ([Z ξ ]) η (16)
to be used to express the velocities (and their gradients) at the front in terms of their jumps across the latter. Here, the linear operator H(., η) is the Hilbert transform in η defined by H(exp(ikη), η) = i.sign(k) exp(ikη), with the convention that sign(0) = 0, whereby H(1, η) = 0 and H(f, η) = 0 for any f (η). The determination of Z ξ as a function of τ requires another piece of information, that is here provided by the need to satisfy [START_REF] Sivashinsky | Nonlinear analysis of hydrodynamic instability in laminar flames. Part 1 : derivation of basic equations[END_REF]. The necessity to handle the transverse averages separately is the point missed in [START_REF] Sivashinsky | On the non linear theory of hydrodynamic instability in flames[END_REF] when attempting to improve the S1 equation: the overall mass conservation in the fresh medium could not be met, whereas it should, and the resulting predicted S T /S L was erroneous; see sections 4.3 and 6.

The S1 equation

At leading order, the Euler equations are automatically satisfied by

U (0) = Π (0) ξ , V (0) = Π (0) η and P (0) = Π (0) ξ , provided Π (0) ξξ + Π (0)
ηη = 0; to leading order, the flow is potential almost everywhere (ξ = 0), which motivated Frankel's [START_REF] Frankel | An equation of surface dynamics modeling flame fronts as density discontinuities in potential flows[END_REF] work. Using the expanded jumps ( 7)- [START_REF] Darrieus | Propagation d'un front de flamme[END_REF], the results quoted above in [START_REF] Frankel | An equation of surface dynamics modeling flame fronts as density discontinuities in potential flows[END_REF] [START_REF] Dold | An evolution equation modeling inversion of tulip flame[END_REF] and the expanded form of (2) gives the leading order evolution equation for F = F (0) + ...

1 S L F t + a 1 (γ) 2 F 2 y + 1 -a 1 2 F 2 y = Ω 1 I (F, y) + 1 k n F yy , (17)
with a 1 (γ) = 1, Ω 1 (γ) = γ/2; I(F, y) ≡ H(-F y , y) is the Landau-Darrieus operator satisfying I(1, y) ≡ 0 and I • I ≡ -∂ yy when operating on periodic functions of y. This equation is Sivashinsky's [START_REF] Sivashinsky | Nonlinear analysis of hydrodynamic instability in laminar flames. Part 1 : derivation of basic equations[END_REF] original result, and is henceforth denoted S1. Because a 1 = 1, the "counter-term" (1a 1 ) F 2 y /2 does not contribute to F t at this order, and ( 17) is compatible with (3). Frankel's vorticity-free model [START_REF] Frankel | An equation of surface dynamics modeling flame fronts as density discontinuities in potential flows[END_REF], that assumes potential flow on both sides of the flame and does not make use of [START_REF] Darrieus | Propagation d'un front de flamme[END_REF], clearly gives an equation (called POT1) that coincides with [START_REF] Karlin | The rate of expansion of spherical flames Combustion Theory and Modelling[END_REF], because the actual flow is piecewise potential to leading order in γ. As for [START_REF] Darrieus | Propagation d'un front de flamme[END_REF] itself, no contradiction occurs at leading order, by Bernoulli's law for potential flows.

The S2 equation

At next order, V (1) ≡ Π (1) η -φ (0) η Π (0) ξ , P (1) ≡ -Π (1) ξ -Π (0) τ and U (1) ≡ Π (1) ξ + R (1) ± , with R (1)
-≡ 0 on fresh side (by Kelvin's theorem), provided P (1) ≡ Π (1) φ (0) Π (0) ξ has a zero Laplacian. Ultimately (see Appendix for main steps of the calculation procedure) this produces an equation, denoted S2, for F = F (0) +γF (1) +... . Equation S2 has the same form as S1 (equation ( 17)), provided a 1 and Ω 1 are replaced therein by a 2 ≡ 1+γ/2 and Ω 2 = γ 2 a 2 , respectively. The difference between a 2 and a 1 stems from the term F y v| x=F - featured in the propagation law (2), since v| x=F - γ 2 S L F y to leading order (by [START_REF] Bychkov | On the dynamics of a curved deflagration front[END_REF] and ( 24)). Thus, the ratio Ω/a is unchanged from S1 to S2, implying [START_REF] Thual | Application of pole decomposition to an equation governing the dynamics of wrinkled flame fronts[END_REF] that no correction on F t from what S1 yields exists, since F 2 y can be shown [START_REF] Thual | Application of pole decomposition to an equation governing the dynamics of wrinkled flame fronts[END_REF] [START_REF] Joulin | On a tentative, approximate evolution equation for markedly wrinkled premixed flames[END_REF] to be proportional to the square of this ratio; see also section 5.

The S3 equation

At 3rd order, the analysis proceeds similarly but is too lengthy to be reproduced in full in the main text; as for second order, the main steps are compiled in Appendix B. The raw equation deduced for φ = φ (0) + γφ (1) + γ 2 φ (2) then reads

φ τ + 1 2 φ 2 η - γ 2 8 φ 4 η = 1 2(1-γ) I(φ) + γ 2 4 H φ η φ 2 η + φ ηη + γ 2(1-γ) φ ηη -γ 2 1 -γ 2 φτ -φτ 1-γ -γ 2(1-γ) φ 2 η -φ 2 η + γ 2 4 (φ 2 H (φ ηη ) + H (φ 2 φ ηη ) -2φH (φφ ηη )) η + γ 2 4 H φ (φ ηη -φ τ ) η -(φH (φ ηη -φ τ )) η + γ 2 2 I -1 φ τ τ + 1 2 φ 2 η τ -1 2 φ ηητ + 1 2 φ ηη -γ 2 2 φ ηη φ 2 η + O (γ 3 ) , ( 18
)
up to O(γ 3 ) terms. Before proceeding any further a remark is due. The pseudo-inverse I -1 (., η) of I(., η) featured above is defined as I -1 (exp(ikη), η) = exp(ikη)/|k|, and makes sense only when operating on functions with vanishing transverse average; accounting for the equation for φ (0) , one can check that this condition is met in (18).

First or second order in time?

The equation for φ obtained as outlined above contains quantities in third-order (in γ) corrections. The choice we made here, to be commented later on, is to express those by means of the leading order (equations S1 and S2) results. The S3 equation for F = F (0) + γF (1) + γ 2 F (2) is then again first-order in time, and reads

1 S L F t + a 3 (γ) 2 F 2 y - 1 8 F 4 y + 1 -a 3 (γ) 2 F 2 y = Ω 3 (γ) I (F, y) + 1 k n F yy + γ 4 H F y F 2 y + γ k n F yy - γ 2k n F 2 y F yy + γ 4 F 2 H (F yy ) + H F 2 F yy -2F H (F F yy ) y + γ 4 H F F 2 y + γ k n F yy -γI (F, y) y - γ 4 F H F 2 y + γ k n F yy -γI (F, y) y , ( 19
)
with a 3 (γ) ≡ 1 + γ 2 + 3 4 γ 2 and Ω 3 (γ) = γ 2 a 3 (γ), up to terms of higher order in γ. We did not attempt to simplify the RHS of ( 19) any further, e.g. upon using S1 and S2 and such identities as H(f H(g)) + H(gH(f )) ≡ H(f )H(g)f g (see page 192 of reference [START_REF] Ablowitz | Solitons, nonlinear evolution equations and inverse scattering[END_REF]). To 3rd order in γ, Frankel's [START_REF] Frankel | An equation of surface dynamics modeling flame fronts as density discontinuities in potential flows[END_REF] vorticity free model produces an equation for F (called POT3) similar to ( 19), yet with a 3 and Ω 3 replaced by (1 + γ/2 + γ 2 /2) and γ(1 + γ + γ 2 )/2 respectively, implying a larger Ω 2 /a 2 ratio. Also, the two I(F, y)'s in the last two lines of [START_REF] Orszag | On the elimination of aliasing in finite-differences schemes by filtering high wavenumbers components[END_REF] would not show up. In other words, vorticity modifies the coefficients (a 3 , Ω 3 ), and adds two terms in the RHS, that are of hydrodynamical origin since they have no 1/k n factor. The potential analog (POT3) to S3 yields solution profiles (cf. figure 2) that are in even better agreement with DNS results than what S3 itself produces, even though the contributions involving the two aforementionned I(F, y)'s again almost compensate one another (like in figure 6). This motivated the conjecture that S1 (=POT1), S2, POT2, POT3, S3, ... (and, hopefully... reality) might mainly differ by the effective values of a i and Ω i they involve, since the corresponding solutions all have nearly identical shapes.

Solutions to the S3 equation, comparison with others'

Numerics

The above S3 equation ( 19) was integrated numerically with 2π/k boxperiodic boundary conditions in y. The method employed [START_REF] Joulin | On a tentative, approximate evolution equation for markedly wrinkled premixed flames[END_REF] combines a Fourier pseudo-spectral scheme in space and an exponential time differencing scheme in time (of first-order accuracy, ETDRK1, [START_REF] Kassam | Fourth-Order Time Stepping For Stiff Pdes[END_REF]). This scheme is exact in the linear limit and accounts for dealiasing [START_REF] Orszag | On the elimination of aliasing in finite-differences schemes by filtering high wavenumbers components[END_REF]. The main difficulties were time-stepping, because the high-order non-linearities may vary quickly during fast transients, and the intricate ways non-linearities and Hilbert transforms are nested, which required quite many direct and inverse Fourier transforms to be performed at any one step of the marching procedure. They compare the obtained steady shapes to those given by a model equation (called S-fit), designed as follows. Equation S-fit :

Steady patterns, S-fit equation

1 S L F t + a(γ) 2 F 2 y + 1 -a(γ) 2 F 2 y = Ω(γ) I(F, y) + F yy k n ( 20
)
has the same structure as S1-S2, with the Ω i 's replaced by the positive root Ω(γ) of (2γ)Ω 2 + 2Ω = γ/(1γ) to yield the same linear growth at |k| k n as the exact Landau-Darrieus analysis; note that the previous Ω i were Taylor-expanded versions of Ω(γ); in particular Ω(γ) = Ω 3 (γ) + o(γ 3 ) as γ -→ 0, whereby the replacement of Ω 3 (γ) by Ω(γ) in S3 is "lesser sin" (asymptotically speaking) for γ 1. The a i (γ)'s were replaced (once for all at a fixed γ) by a(γ) a 3 (γ) to make the solutions to S3 and S-fit have the same amplitude of wrinkling for this choice of k box ; the needed correction a 3 (γ)a(γ) is positive, and goes like 3 8 γ 2 if γ 1. All the other nonlinearities featured in S3 were simply discarded in S-fit. These are bolder approximations, principaly justified by their consequences.

Clearly the profiles resulting from S3 and S-fit are very close to one another (fig. 3 and4), and are so for all γ ≤ 0.9 (ρ u /ρ b ≤ 10). The predicted flame speeds also are, by [START_REF] Sivashinsky | Nonlinear analysis of hydrodynamic instability in laminar flames. Part 1 : derivation of basic equations[END_REF]. This fact (see also next sections) could not be given a neat explanation; yet a partial one likely is accessible in the framework of Frankel's model [START_REF] Frankel | An equation of surface dynamics modeling flame fronts as density discontinuities in potential flows[END_REF], since the potential flow analog (called POT3) to S3 also exhibits this puzzling property, and the corresponding evolution equation is available [START_REF] Frankel | An equation of surface dynamics modeling flame fronts as density discontinuities in potential flows[END_REF] for any γ. The figures 5-7 and 6 suggest how this may happen. As shown in figure 5, there is a nearly perfect compensation between cubic nonlinearities; a similar, yet less perfect (cf. the scales) cancellation occurs (for γ = 0.8) between cubic and quadratic nonlinearities (fig. 6) , both of which are negligible if γ 1. Overall, (figure 7), the collection of all extra nonlinearities featured at 3rd order nearly modifies the term a 3 (γ)F 2 y -wherever important -by a nearly constant factor ( a(γ)/a 3 (γ)), the main origin of which comes from the quartic term F 4 y /8 featured in the LHS of S3 (all the other extra terms are collectively small and vanish on transverse average): writing a 3 F 2 y -F 4 y /4 as (a 3 -F 2 y /4)F 2 y indicates how the influence of F 4 y might be qualitatively mimicked upon replacing a 3 by a ≤ a 3 , with a 3a = O(γ 2 ), since F 2 y = O(γ 2 ) ≥ 0; the surprise is that a quantitative , for k box = k n /2 and γ = 0.8 : agreement could also be reached. One then conjectures that the whole formal series involved in (3) will approximately behave like F 2 y , up to some γdependent factor (see equation [START_REF] Boury | Etudes théoriques et numériques de fronts de flammes plissées : dynamiques non-linéaires libres ou bruitées[END_REF] and figure 12 below).

γ 4 H F (F 2 y + γ k n F yy -γI(F )) y (solid), - γ 4 F H F 2 y + γ k n F yy -γI(F ) y (

Unsteady S3 and S-fit equations

The a(γ) coefficient in S-fit being selected as above (for k box = k n /2 and in steady cases), how do the solutions to S3 and S-fit equations compare for unsteady flames and/or with k box 's different from k n /2 ? Figure 8 compares those for γ = 0.8 and for the common initial condition F (y, 0) ∼ cos(k box y), k box = k n /2 (corresponding to the most rapidly growing mode) once both results are plotted in terms of F (u t, t), with u = 8S L , so as to represent a smooth oblique cut of the y-periodic, two-variable array F (y, t). Put in words, the S-fit equation nicely captures most of S3's transients, including the process of wrinkle coarsening through crest mergers (fig. 10), even for realistic γ's. The S-fit equation being a MS equation can be solved analytically [START_REF] Thual | Application of pole decomposition to an equation governing the dynamics of wrinkled flame fronts[END_REF], yielding a simple expression for F 2 y when F t is a constant; see (21) below.

S-fit vs direct simulation of steady flames

The results above made it tempting to compare the steady front-cell shapes issued from the S-fit equation to what direct numerical simulations (DNS) of a full reactive Navier-Stokes problem produce. Figure 11 illustrates such a comparison with the DNS by Kadowaki [START_REF] Kadowaki | The influence of hydrodynamic instability on the structure of cellular flames[END_REF] corresponding to γ = 0.8, a unit Lewis number Le, a one-step, one-reactant Arrhenius rate with a large Zel'dovich number Ze = 10, and k box = k n /2. Provided a(γ) is again adjusted (in a suitable O(γ 2 ) -13% neighbourhood of a 3 (γ)) to reproduce the correct amplitude of wrinkling as given by DNS, the agreement between front shapes is quite impressive; that on S T /S L -1 also is, since equation ( 3) was numerically verified by the DNS to hold to better than 1 % when the front is identified with the thin reaction zone (where ρ/ρ b 1 + 1/Ze). The latter property can actually be demonstrated analytically for Le = 1 and Ze 1 [START_REF] Boury | Etudes théoriques et numériques de fronts de flammes plissées : dynamiques non-linéaires libres ou bruitées[END_REF], since the reaction-zone structure then nearly is the same along its normal n as in steady flat flames (see Appendix A).

Our final comparison deals with the variations of S T /S L -1 with γ and k box /k n . For steady front shapes, the S-fit equation predicts [START_REF] Joulin | On the hydrodynamic stability of flat-burner flames[END_REF] F 2 y /2 to be of the form

F 2 y 2 = g(γ)h(k box /k n ) , g(γ) ≡ Ω 2 (γ) 2a 2 (γ) (21)
with [START_REF] Joulin | On the hydrodynamic stability of flat-burner flames[END_REF]

h k box k n = 4N k box k n 1 -N k box k n ( 22
)
where N is a strictly positive integer, less than or equal to the integer-part, and the symbols, redrawn from [START_REF] Bychkov | Nonlinear equation for a curved stationary flame and the flame velocity[END_REF], are from DNS. The dotted lines correspond [START_REF] Vaynblat | Stability of pole solutions for planar propagating flames : I. Exact eigenvalues and eigenfunctions / II. Properties of eigenvalues and eigenfuctions with implication to flame stability[END_REF] to unstable, two-pole steady solutions to S-fit.

N max , of 1 2 (1 + k n k box ). For N = N max ,
(otherwise, S T ≡ S L since the flat flame is then stable). Figure 12 compares the DNS of Bychkov et al. [START_REF] Bychkov | Nonlinear equation for a curved stationary flame and the flame velocity[END_REF] with what (21)( 22) predict. Apart from two points coinciding with unstable (but long-lived [START_REF] Joulin | On a tentative, approximate evolution equation for markedly wrinkled premixed flames[END_REF], and possibly interpreted as converged DNS patterns) steady solutions to S-fit, the agreement with what (21)( 22) predict again is worthnoticing; in particular, the factoring of the S T /S L -1 values found by DNS into such an expression as

m(γ)h(k box /k n ) is remarkable.
The best fit of a(γ) we could achieve as to make the S-fit prediction on S T /S L -1 (i.e. m(γ) = Ω 2 /2a 2 ) coincide with some available DNS results [START_REF] Kadowaki | The influence of hydrodynamic instability on the structure of cellular flames[END_REF][START_REF] Bychkov | On the dynamics of a curved deflagration front[END_REF] was

a(γ) = 1 + 1 2 γ + 3 8 γ 2 + 4 3 (1 -γ) -1/4 -(1 + γ 4 + 5γ 2 32 ) (23)
It agrees with the required correction a 3a 3 8 γ 2 needed to make the S3 and S-fit predictions coincide at γ 1. Because Ω 2 (γ -→ 1 -) 1/(1γ), the above functional form is also consistent with the expected asymptotic behavior S T /S L ∼ (1γ) -1/2 as γ -→ 1 -. Expressing [START_REF] Ruetsch | Effects of heat release on triple flames[END_REF] that the pressure variations (∼ S 2 T ρ u , by Bernoulli's law) ahead of the flame are mechanically compatible with those (∼ (S L ρ u ) 2 (1/ρ b -1/ρ u )) across the front, one indeed arrives at S T /S L ∼ (ρ u /ρ b ) 1/2 when ρ u /ρ b -→ ∞ (i.e. γ -→ 1 -) and curvature effect are omitted; the fit (23) gives 9/32 as proportionality coefficient.

Experiments and S-fit predictions

Once augmented by a term u F y in the LHS to account for an imposed velocity u S L parallel to the mean flame location, a model equation with the same structure as S-fit, with the same a(γ), accurately describes the then convective growth and saturation of wrinkles along the front of an experimental wedge-shaped (2D) Bunsen flame [START_REF] Searby | Comparison of experiments and a nonlinear model equation for spatially developing flame instability[END_REF].

A model equation for fast transients

As mentionned earlier, we made the rather arbitrary choice of expressing the inertia-related terms

F (0) t + S L F (0) 2 y 2 t
, appearing at 3rd order in the γ-expansions, in terms of the lowest order, S1 results, thereby obtaining the S3 equation displayed in [START_REF] Orszag | On the elimination of aliasing in finite-differences schemes by filtering high wavenumbers components[END_REF], that is first-order in time. It is known, however, that the exact linearized flame dynamics corresponding to [START_REF] Sivashinsky | On the non linear theory of hydrodynamic instability in flames[END_REF][START_REF] Kassam | Fourth-Order Time Stepping For Stiff Pdes[END_REF][START_REF] Williams | Combustion Theory[END_REF][START_REF] Markstein | Non steady flame propagations[END_REF][START_REF] Bychkov | On the dynamics of a curved deflagration front[END_REF][START_REF] Darrieus | Propagation d'un front de flamme[END_REF][START_REF] Landau | On the theory of slow combustion[END_REF] is 2nd order in time [START_REF] Markstein | Non steady flame propagations[END_REF][11] [START_REF] Landau | On the theory of slow combustion[END_REF], small harmonic wrinkles (F ∼ exp(ωt + iky)) being governed by ω = S L |k|Ω(γ, |k|/k n ), with Ω(γ, |k|/k n ) given by either roots Ω ± of

(2 -γ)Ω 2 + 2Ω 1 + γ|k| 2k n = γ 1 -γ 1 - |k| k n ; ( 24
)
the previous Ω(γ) was Ω + (γ, 0). Equivalently, the exact linearized equation for F (y, t) is

(2 -γ) (F tt ) + 2S L I F t + γ 2k n I(F t ) = γ 1 -γ S 2 L I (I (F ) + F yy /k n ) (25)
in which the second argument y of I(., y) was omitted for clarity. The following interpolating model equation, that shares many features with that [START_REF] Dold | An evolution equation modeling inversion of tulip flame[END_REF] proposed to mimic tulip flame inversions, reads

(2 -γ) (F tt -F tt ) + 2S L I F t + b(γ) 2 S L F 2 y + γ 2k n I(F t ) = γ 1 -γ S 2 L I (I (F ) + F yy /k n ) + (g -F tt )γI(F ). ( 26
)
Combined with (3) as to evaluate F tt ≡ F t t it presents several interesting, complementary, features : i) in the linear limit, it reduces to the exact result [START_REF] Dold | An evolution equation modeling inversion of tulip flame[END_REF], and to its generalization [START_REF] Markstein | Non steady flame propagations[END_REF] when a gravity (of intensity g) along the x-axis is included;

ii) for steadily propagating patterns (F t =const.) and zero gravity, its solutions coincide with what a S-fit type of equation would give, provided b(γ) is selected to yield 2b(γ)(1-γ)Ω(γ)/γ ≡ a(γ) (or any accurate empirical fit of DNS);

iii) for γ 1, equation ( 26) resumes the asymptotically correct equation needed to account for weak gravities [START_REF] Boury | Etudes théoriques et numériques de fronts de flammes plissées : dynamiques non-linéaires libres ou bruitées[END_REF] [START_REF] Boury | Mean cell wavelength selection of wrinkled premixed flames in weak gravity fields : spontaneous evolutions[END_REF]. It is still meaningful when g is a function of time, and can then be used as a model to take up the nonlinear developments of parametric instabilities [START_REF] Markstein | Non steady flame propagations[END_REF] and the partial (or total) suppression of Landau-Darrieus instability by oscillating g's.

iv) it is also compatible with the fact that, even for g = 0, F tt -S L F 2 y /2 t + ... = 0 provides the flame with a gravity effect in its own frame; as first pointed out in [START_REF] Gostintsev | Self-similar propagation of a turbulent flame in a mixed gaseous mixture Comb[END_REF], this might well be a key ingredient, but possibly not the only one [START_REF] Boury | Etudes théoriques et numériques de fronts de flammes plissées : dynamiques non-linéaires libres ou bruitées[END_REF], to explain the "self" acceleration of outwardly expanding wrinkled flame fronts [START_REF] Angelo | On model evolution equations for the whole surface of three-dimensional expanding wrinkled premixed flames Combustion Theory & Modeling[END_REF].

Conclusions

The expansion at small density contrasts γ ≡ (ρ u -ρ b )/ρ u 1, pioneered in [START_REF] Sivashinsky | Nonlinear analysis of hydrodynamic instability in laminar flames. Part 1 : derivation of basic equations[END_REF], was extended -using the same perturbative approach -up to the 3rd order in γ and did not encounter inconsistancy nor difficulty, apart from the length of the computations. The results (flame shape and velocity) were found to smoothly approach DNS results on steady fronts as higher powers of γ are included, even for realistic γ's. Even though γ -→ 0 + implies γ ≡ ρ u ρ b -1 -→ 0 + , the two limit processes are not equivalent when the radius and speed of convergence are concerned. To wit, the convergence of power series expansions of

γ 1 + γ = γ -γ 2 + ...(≡ γ) or of log(1 + γ ) =
γγ 2 /2 + ...(≡log(1γ)) = γ + γ 2 /2 + ... are noticeably accelerated upon use of γ > 0 instead of γ > 0 as expansion parameter, and is not limited to γ < 1 any longer : whatever γ > 0, γ is less than 1. This likely contributed to the success of γ-expansions when coping with the flame problem considered here. Yet it is not known at present whether the series in γ of the present flame are convergent or merely formal. Besides, the obtained evolution equation at 3rd order in γ included has solutions that are tantalizingly similar to those of a model equation (S-fit) of the Michelson-SIvashinsky type, with a unique renormalized coefficient (we insist that this renormalization is done once for all at a given γ). All this fosters the conjecture that a simple underlying analytic structure, akin to what pole decompositions [START_REF] Thual | Application of pole decomposition to an equation governing the dynamics of wrinkled flame fronts[END_REF] revealed for Sivashinsky's equation and its second order in γ extension, would still be present or acting at 3rd order, and possibly at all orders (see the DNS results plotted in figures [START_REF] Darrieus | Propagation d'un front de flamme[END_REF][START_REF] Landau | On the theory of slow combustion[END_REF]. or that there would exist a simple, e.g. integral, argument to bound the error made when using S-fit type of equations. A similar remark holds out about Frankel's potential flow model. Contrary to the latter, no simple non-perturbative (cf. [START_REF] Kazakov | On-shell description of stationary flames[END_REF][START_REF] Joulin | Potential-flow models for channelled two-dimensional premixed flames around near-circular obstacles[END_REF]) approach to any non-trivial vortical model of flame propagation is as yet available. Still, due to its simplicity, the S-fit model equation and its 2nd order (in time) generalization could be of some help one get semi-quantitative answers (or, at least, make predictions that can be tested) on flames, including 3D planar [START_REF] Denet | Sivashinsky equation in a rectangular domain[END_REF] or 3D expanding ones [START_REF] Angelo | On model evolution equations for the whole surface of three-dimensional expanding wrinkled premixed flames Combustion Theory & Modeling[END_REF][START_REF] Karlin | The rate of expansion of spherical flames Combustion Theory and Modelling[END_REF]; or in some astrophysical contexts where γ is known to be comparatively small [START_REF] Timmes | The conductive propagation of nuclear flames[END_REF]. At this order, the kinematic relation writes

φ (1) τ + φ (0) η φ (1) η = U (1) (0 -) -φ (0) η V (0) (0 -) + 1 2 φ (1) ηη ,
Substituting directly U (1) (0 -) and V (0) (0 -) by the above expressions, one obtains

φ (1) τ + φ (0) η φ (1) η = 1 2 η [V ] τ - η ] τ + 1 2 ([P ] -[P ] ) (B.6) + 1 2 φ (0) η [V ] -φ (0) η [V ] + 1 2 φ (1) ηη + O (γ) .
One can NOW use the expansion for the jumps

[V ] = - φ η 1 -γ + O(γ 2 ) , η [V ] τ = -φ (0) τ + O(γ) , [P ] = φ (0) ηη + O(γ) ,
Then combining the two kinematic relations at the 2 first orders yields

φ τ - γ 2 φ τ + 1 2 1 + γ 2 φ 2 η - γ 2 φ 2 η = 1 2 1 + γ 2 (I (φ) + φ ηη ) + O(γ 2 ) .
Using the kinematic relation at first order (in order to compute γ φ τ = -γ φ 2 η /2 + O(γ 2 )), this finally yields the S2 equation (in dimensionalized form) :

F t + S L 2 1 + γ 2 F 2 y -S L γ 4 F 2 y = S L γ 2 1 + γ 2 I (F ) + 1 k n F yy . (B.8)

Appendix C. Resolution at third order

At third order in γ, the Euler equations write

U (2) ξ + V (2) η = φ (0) η V (1) ξ + φ (1) η V (0) ξ , (C.1a) U (1) τ + U (2) ξ + P (2) ξ = U (0) τ -U (0) -φ (0) τ U (0) ξ -V (0) ξ U (0) η , (C.1b) V (1) τ + V (2) ξ + P (2) η = V (0) τ -U (0) -φ (0) τ V (0) ξ -V (0) ξ V (0) η + φ (0) η P (1) ξ + φ (1) η P (0) ξ . (C.1c)
The two momentum equations (C.1b) and (C.1c) are verified if 1) τ (η , τ)dη ,

U (2) = Π (2) ξ -R (1) τ ξ + R (2) (η, τ ) , V (2) = Π (2) η -φ (0) η Π (1) ξ -φ (1) η Π (0) ξ + η R (
P (2) = -Π (2) ξ -Π (1) τ + Π (0) τ - 1 2 Π (0) ξ 2 - 1 2 Π (0) η 2 + φ (0) τ Π (0) ξ ,
and the continuity equation (C.1a) shall impose Π

ξξ + Π (2) ηη = φ (0) ηη Π

ξ + φ (1) ηη Π (0)

ξ + 2φ (0) η Π (1)
ξη + 2φ (1) η Π (0)

ξη -φ (0) η 2 Π (0)
ξξ .

But one may write Π (2) = φ (0) Π

(1)

ξ + φ (1) Π (0) ξ -1 2 φ (0) 2 Π (0)
ξξ + P (2) in order to once again get P

ξ ] = φ (1) [V ] η + 1 2

η [V ] τ .
Equation (C.4a) allows one to compute the vorticity R (2) (η, τ ) in the burnt side : At this order, the kinematic relation becomes : After solving equation (C.3), one gets U (2) (0 -), that can be substituted in (C.5) to yield

φ (2) τ + 1 2 φ (1)
φ (2) τ + 1 2 φ (1)
η 2 + φ (0) η φ (2) η -

1 8 φ (0) η 4 = + 1 4 φ (0) 2 I ([V ]) η - 1 2 φ (0) H η [V ] τ + [P ] -φ (0) [V ] η η + 1 2 H φ (0) [V ] τ + 1 2 H φ (0) [P ] η + 1 4 I φ (0) 2 [V ] η - 1 2 I -1 [U] + η [V ] τ + [P ] τ - 1 4 η [V ] τ - η [V ] τ + 1 2 φ (1) η [V ] -φ (1) η [V ] + 1 2 φ (2) ηη - 1 2 φ (0) ηη φ (0) η 2 + O (γ) . (C.6)
Then, computing the expansion at 3rd order (in γ 2 ) for the kinematic relation, one may get One can then NOW introduce the expansion for the jumps, to get

φ τ + 1 γ 2 1 + γ 2 φ 2 η -1 = 1 2 H ([V ]) + γ 2 ([P ] -[P ] ) + γ 2 1 - γ 2 η [V ] τ - η [V ] τ + γ 2 (φ η [V ] -φ η [V ] ) - γ 2 4 φ 2 H ([V ] η ) + H φ 2 [V ] η -2φH (φ[V ] η ) η + γ 2 2
[V ] = - φ η 1 -γ + γ 2 2 φ η (φ η ) 2 + φ ηη + O γ 3 , γ η [V ] τ = - γ 1 -γ φ τ + O γ 3 , γ 2 [U] = - γ 2 2 (φ η ) 2 + φ ηη + O γ 3 , γ 2 [W ] = γ 2 (φ ηη -φ τ ) + O γ 3 , γ 2 ([U] + [W ]) = -γ 2 φ τ + 1 2 (φ η ) 2 - 1 2 φ ηη + O γ 3 .
Replacing directly these expression into (C.7) yields the Evolution Equation

φ τ + 1 2 φ 2 η - γ 2 8 φ 4 η = 1 2(1 -γ) I(φ) + γ 2 4 H φ η φ 2 η + φ ηη + γ 2(1 -γ) φ ηη - γ 2 1 - γ 2 φ τ -φ τ 1 -γ - γ 2(1 -γ) φ 2

Figure 1 :

 1 Figure 1: Two-dimensional configuration : a plane-on-average flame front propagating along a 2-D channel. The coordinates are such that the far-upstream velocity is S L along the x-axis, whatever F (y, t).

Figure 2 :

 2 Figure 2: Comparison of steady half-cell shapes given by various evolution equations (S1, S2, S3 and POT1, POT2, POT3) for k box = k n /2 and γ = 0.8. The symbols are from the direct simulation results of Kadowaki[START_REF] Kadowaki | The influence of hydrodynamic instability on the structure of cellular flames[END_REF], obtained in the same conditions.

Figures 3 -

 3  show half the final F (y, t) -F t t resulting from S3 when k box = k n /2 (the most unstable mode) and for γ = 0.8 (i.e. ρ u = 5ρ b ) or γ = 0.9 (ρ u = 10ρ b ).They compare the obtained steady shapes to those given by a model equation (called S-fit), designed as follows. Equation S-fit :

Figure 3 :Figure 4 :

 34 Figure 3: Steady half-cell shapes corresponding to the S-fit (solid lines) and S3 (symbols) evolution equations for k box = k n /2 and γ = 0.8.

Figure 5 : 4 F 2 Figure 6 :

 5426 Figure 5:Comparison of some cubic terms of the MS3 equation (19) :γ 4 F 2 H(F yy ) y + H F 2 F yy y (solid line), -γ 2 [F H (F F yy )] y (dashed line) and their sum (dotted), for k box = k n /2 and γ = 0.8. Note that only one fourth of the cell length is shown; beyond, the sum is undistinguishable from zero at the scale shown.

Figure 7 :

 7 Figure 7: The quadratic term a 3 F 2 y (solid line) compared to the sum of all nonlinear terms featured in S3 (dotted line); their difference is plotted as the dashed line.

Figure 9

 9 does the same for k box F (y, 0) = 0.2 cos(k box y)+0.1 cos(2k box y-0.2), again with γ = 0.8 but for k box = k n /4 and u = 4S L . At the scale of fig.9, the two histories are undistinguishable; yet fig.10exhibits the differences in trajectories of local front maxima/minima, and shows how small they are.

Figure 8 :

 8 Figure 8: F (u t, t) vs time (t) as predicted from the S3 (solid line) or S-fit (symbols) equations, both for u = 8S L , k box = k n /2, γ = 0.8, and the same harmonic F (y, 0).

Figure 9 :

 9 Figure 9: Same as in figure 8, but for u = 4S L , k box = k n /4 and k box F (y, 0) = 0.2 cos(k box y) + 0.1 cos(2k box y -0.2). The two curves are undistinguishable. Note the crest merger.

20 -

 20

Figure 10 :

 10 Figure 10: Enveloppes of flame crests (top) or troughts (bottom) as predicted by the S3 (solid lines) and S-fit equations (dashed); same conditions as in figure 9.

  corresponding to the only stable steady pattern[START_REF] Vaynblat | Stability of pole solutions for planar propagating flames : I. Exact eigenvalues and eigenfunctions / II. Properties of eigenvalues and eigenfuctions with implication to flame stability[END_REF] h(.) is piecewise parabolic with maxima of 1 if k box < k n

Figure 11 :Figure 12 :

 1112 Figure 11: The predictions of S-fit (solid lines) vs. the DNS results of [20], for γ = 0.8 and k box = k n /2.

Figure A. 13 :

 13 Figure A.13: Sketch of a steady 2D flame in a channel, in a frame attached with its mean location. The dashed line represents the control volume; its rightmost boundary has Y 1, hence ρ ρ b and < T T b

( 2 )

 2 ξξ + P(2) ηη = 0 .(C.3)The jump [V ] has already been taken into account at first order, and the jumps [U] and [P ] at second order; hence the jumps to be considered now are reduced to zero : η ] = φ (0) [V ] τ +φ(0) [P ] η -

U ( 2 ) 2 +

 22 (0 -)φ (0) η V (1) (0 -)φ (1) η V (0) (0 -O (γ) .

γ 2 2 I 2 φ ηη 1 + γ 2 φ 2 η+ O γ 3 ,

 22123 H (φ[W ] η ) -(φH ([W ])) η --1 (([U] + [W ]) τ ) + 1 η [V ] τ + [P ].

  H (φ ηη ) + H φ 2 φ ηη -2φH (φφ ηη ) η (φ ηηφ τ ) η -(φH (φ ηηφ τ )) η >We wish to derive pertinent evolution equations for premixed wrinkled flames >The Michelson-Sivashinsky perturbative approach is extended to 3rd order in gamma >Obtained flame shape and speed quantitavely compare to 1rt order equation S1 >An S1-type model equation is proposed, with only one gamma-dependent coefficient

									η
	+ φ 2 + γ 2 4 γ 2 4 H φ + γ 2 2 I -1 φ τ τ +	1 2	φ 2 η τ -	1 2	φ ηητ +	1 2	φ ηη -	γ 2 2	φ ηη φ 2 η + O γ 3 .
									(C.9)

ηφ 2 *Highlights

Aknowledgements

The authors wish to thank Senior Researcher G. Joulin for sharing unpublished results, and also for many fruitful suggestions and discussions. One of the author (YD) benefitted from discussions with X.-S. Bai, V. Sabelnikov, N. Swaminathan, C.O. Paschereit, A. Gruber, N.E. Haugen, E. Oran, V. Yakhot and other participants of the Nordic Institute for Theoretical Physics (NORDITA) Programme on Turbulent Combustion in 2010.

Appendix A. Deriving S T /S L for steady flame shapes, if Le=1

and Ze 1

In [START_REF] Kadowaki | The influence of hydrodynamic instability on the structure of cellular flames[END_REF][START_REF] Bychkov | On the dynamics of a curved deflagration front[END_REF], the DNS made use of a one-reactant (mass fraction Y ), one-step Arrhenius kinetics of local rate ω(T, Y ), with a large activation temperature T a . The DNS also assumed constant values of heat conductivity λ, specific heat at constant pressure c, density times reactant diffusivity ρDc and a unit Lewis number Le ≡ λ/ρDc. The first integral Combining the energy balance (cu.∇T = λ∇ 2 T + Qω), the mass conservation (∇.(ρu) = 0) and the equation of state (ρ ∼ 1/T ) yields

Upon integration over the control volume of the figure A.13, equation (A.1) produces

where

(.)dl is shorthand for the integration of (.) along the control- Most of the above reasoning can be adapted to show that (3) still approximately holds, yet with a somewhat larger uncertainty, whenever Le -1 = O(1/Ze).

Appendix B. Resolution at second order in γ

At second order, Euler equations write

ξ + P

(1) ξ = 0 , (B.1b)

and the two momentum equations (B.1b) and (B.1c) now verified if there exist Π (0) (ξ, η), Π (1) (ξ, η) and R (1) (η, τ ) such that

τ . (recall that underlined terms only exist in the burned side).

The continuity equation (B.1a) shall impose Π

ξξ + Π (1)

Writing Π (1) = φ (0) Π (0) ξ + P (1) , one can again get a Laplace equation P

η ] = 0. From (B.4c), one deduces [P


Tryfon Antonakakis

Fadi Baïda

Abderrahmane Belkhir

Kirill Cherednichenko

Shane Cooper

Richard Craster

Guillaume Demesy

John Desanto

Gérard Granet

Editor Evgeny Popov

Boris Gralak

Sébastien Guenneau

Daniel Maystre

André Nicolet

Brian Stout

Fréderic Zolla

Benjamin Vial

GratinGs: theory and numeric applications

Editorial Preface

A typical question that almost all of us (the authors' team and other colleagues) has been asked not only once has in general the meaning (although usually being shorter): "What is the best method for modeling of light diffraction by periodic structures?" Unfortunately for the grating codes users, and quite fortunately for the theoreticians and code developers, the answer is quite short, there is no such a bird like the best method.

In the more than 30 years active studies on the subject, I have worked on the theory and numerical applications of several approximate methods, like Rayleigh expansion, coupled-wave theory, beam propagation method, first-order approximations, singular Green's function approximation, effective index medium theory, etc. My conviction is that they are quite useful (otherwise why to exist) for physical understanding, but my heart lies in what is considered as rigorous grating theories. Name 'rigorous' is used in the sense that in establishing the theories, exact vector macroscopic Maxwell equations and boundary conditions are applied without approximations. The approaches become approximate after computer implementation, due to the impossibility to work with infinite number of equations and unknowns, and due to the finite length of the computer word.

Of course, there are always initial approximations and assumptions, like the infinite dimensions of the grating plane, linearity of the optical response, etc. From physical point of view, the main feature of the methods, presented in this book are characterized by the use of optical parameters of different substances as something given by other physical optics theories and the experiment as an ultimate judge.

The necessity to use more than a single rigorous method comes from practice: different optogeometrical structures made of different materials and working in different spectral regions require a variety of methods, because each one is more effective in some cases, and less effective (or failing completely) in others. In addition, each approach is a subject of constant research and development. Grating modeling, grating manufacturing and grating use go hand in hand, and practice provides strong stimuli for the theory development. Vice versa, recent grating technologies and application cannot advance without proper theoretical and numerical support.

When I started my grating studies, the method of coordinate transformations that uses eigenvector technique to integrate the Maxwell equations (sometimes known as the Cmethod) has just been formulated. It worked perfectly for holographic grating whatever the polarization and the grating material, but failed completely for grooves with steep facets. It took more than 15 years to refine its formulation, so that now it can deal with echelles and pyramidal bumps (in the case of two-dimensional periodicity) with slopes up to 87 deg steepness. However, the method is not at all adapted to lamellar gratings. On the other hand, the Fourier modal method (also known as Rigorous coupled-wave approach, RCW) is perfect for such profiles, but its use in the case of arbitrary grating profiles (e.g., sinusoidal or triangular profiles) in case of metallic grating material causes problems when using a staircase approximation of the profile. The differential method does not use this approximation and could deal with arbitrary profiles, but it took more than 20 years to make it working with metallic gratings in TM (p, or S) polarization. And quite ironically, the improvement came from advances in the competing RCW approach.

These methods are relatively easy for programing nowadays, after solving the numerical problems due to growing exponentials and factorization rules of the product of permittivity and electric field, however there are still some persisting problems for highly conduction metals. In addition, neither the differential, nor the Fourier modal methods can deal with infinitely conducting gratings.

Several methods are quite flexible concerning the geometry of the diffracting objects and the grating material. For example, the integral method can treat inverted profiles, rod gratings with arbitrary cross section, finitely or infinitely conducting materials in any polarization, but its programming require deep mathematical understanding of the singularities and integrability of the Green's functions. Other two flexible methods are quite famous and widely used, even in the form of commercially available codes. These are the finite-element method, and the finite-difference time domain method. The flexibility with respect to the geometrical structure, optical index inhomogeneity and anisotropy, etc. has to be paid by the necessity of sophisticated meshing algorithms and very large sparse matrix manipulations.

These few examples represent only the top of the iceberg, and are invoked to illustrate the basic idea that the best method has not been invented, yet. Probably never.

We have tried to gather a team of specialists in rigorous theories of gratings in order to cover as large variety of methods and applications as practically possible. The last such effort dates quite long ago, and it has resulted in the famous Electromagnetic Theory of Gratings (ed. R. Petit, Springer, 1980), a book that has long served the community of researchers and optical engineers, but that is now out of press and requires a lot of update and upgrade, something that we hope to achieve, at least partially with this new book.

Our choice of electronic publishing is determined by the desire to ensure larger free access that is not easily available through printed editions. I want to thank all the contributors to this Edition. Special thanks are due to my colleagues Fréderic Forestier and Boris Gralak for the technical efforts to make the electronic publishing possible.

Marseille, FranceEvgeny Popov December 2012