Stochastic characterisation of unstable lean hydrogen–air annular premixed flames

This study is devoted to the characterisation of lean premixed hydrogen–air laminar flames stabilised on an annular burner. The flames studied are chosen to lie close to the lean flammability limit, spanning a range of Damköhler and Reynolds numbers, and correspond to situations where instabilities arise. The experimental characterisation is achieved by examining the statistical moments of the flame hydroxyl radical (OH*) chemiluminescence and schlieren. The standard deviation and the skewness, that are needed to fully characterise the measured probability density functions, are shown to exhibit a non-monotonic behaviour with the distance from the burner surface, first increasing, then decreasing downstream the flame tip. An exponentially modified Gaussian distribution model is proposed to describe the skewed schlieren fluctuations probability density function. This model is closed by developing relations for the probability density function parameters based on the Damköhler and Reynolds numbers, and then used to describe the OH* stochastic behaviour. The relation between the distribution parameters and the transport equation of a reactive scalar gradient is addressed also.


Introduction
Since it can be synthesised entirely from renewable sources, the hydrogen energy carrier has the advantage of not producing greenhouse carbon compounds or soot particles during its combustion with air.In comparison with near-stoichiometric mixtures, the use of fuel-lean mixtures in systems destined for power generation is interesting because the associated flame temperatures impose lesser thermal stresses on combustion chambers walls, of the lower fuel consumption, and of the smaller nitrogen oxides emissions.However, the hydrogen-air combustion process can be intrinsically unstable in fuel-lean situations, which should be anticipated to ensure a safe and reliable operation of new practical systems which could be designed to operate using ultra-lean mixtures.Therefore, this work investigates the behaviour of unstable hydrogen-air flames, close to the lean flammability limit, stabilised on a newly developed annular burner.
The stable behaviour of premixed laminar flames of hydrogen and air is well known and has been reported in numerous studies, as can be seen in a recent review of the literature [1].This stable behaviour is observed when the Lewis number controlling the combustion is sufficiently close to one, which is the case for near-stoichiometric mixtures, or when the intrinsically unstable modes are unable to develop.When the concentration of hydrogen in the mixture is much lower than the stoichiometry, thermo-diffusive instabilities may arise, which characteristics have been shown to be controlled by the Lewis number of the limiting reactant [2][3][4].These thermo-diffusive instabilities lead to flame surface wrinkling that eventually evolves to a cellular structure in unconstrained flames.
The presence of such instabilities in lean hydrogen mixtures has been observed in several flame configurations, notably in spherical flames propagation studies [5][6][7][8], and in flames stabilised in planar burners [9][10][11][12].In particular, the first of these configurations has enabled us to identify the effects of the mixture equivalence ratio, temperature and pressure on the growth of such instabilities.More specifically, the cellular structure that results from the unstable behaviour, and manifests as flame surface wrinkling, was found to arise earlier during the spherical flame expansion as the equivalence ratio decreases [8].However, the inherent transient nature of the spherical flames makes it difficult to obtain large samples of the observed phenomena, and thus to perform stochastic analysis.On the contrary, on studies performed on planar burner configurations the long time available for the observations permitted to characterise the dynamic response of the flame front instabilities.This characterisation evidenced, in particular, the effect of the mixture hydrogen mass fraction on the flame front inherent wrinkling and on the corresponding dynamic modes energy content [12].The heat transfer to the burner surface in this planar configuration seems to preclude flame stabilisation in conditions close to the lean flammability limit, though.
Considering the two-dimensional, freely propagating flame configuration, recent numerical studies have addressed both the early, linear, unstable flame development, and also the long term-intrinsic instability behaviour [13][14][15].More precisely, a stochastic analysis of the flame consumption speed and characteristic thickness evidenced the corresponding large variations along the unstable flame front [13].This canonical configuration has enabled to identify and characterise the origin of the flame instabilities.However, even if these studies provide valuable information on the thermo-diffusive instability influence on the flame front, these computed quantities are not readily measurable.Furthermore, practical combustors may offer limits on the flame surface growth, due to a finite premixed flame height, for instance.As a consequence, the instabilities development could be constrained along the flame.
In order to describe the unstable combustion phenomena of lean premixed flames, models of high physico-chemical fidelity are often deployed, in particular those that can be solved by computational fluid dynamics techniques of the balance equations using detailed chemical kinetics and molecular transport descriptions [16][17][18].However, these models are still associated with computation times and costs that are too high to be directly usable in real-time observation and control of industrial processes.To predict and control these unstable behaviours in real time, it is necessary to develop models capable of anticipating them from sensor data that may be installed in the practical systems [19,20].More specifically, it is interesting to combine the detailed results from this type of high-fidelity model with those from measurements that can be deployed in an industrial combustion process environment, thus eventually enabling the construction of reducedorder models [21][22][23][24][25].In particular, stochastic models, based on transported [26,27] or prescribed [28,29] probability density functions of the combustion process properties, are often used to describe hydrocarbon-air combustion in situations where inherent flame instabilities are not relevant.However, the application of these techniques to inherently unstable premixed hydrogen-air combustion could lead to modifications of these models.This is why the knowledge of the different stochastic moments of the flame properties probability density functions is required.
The objective of this study is to study unstable lean premixed hydrogen-air flames stabilised in an annular burner, characterised by a large mixture injection perimeter-to-slit ratio.This configuration has been chosen to enable the stabilisation of flames in situations where circumferential unstable modes develop, but longitudinal ones are constrained by the flame height.More specifically, the aim is to propose a stochastic model of the measured flow-field properties in situations where the flow and the chemical time scalesand thus the Reynolds and Damköhler numbers -exhibit large variations, which is the case when the mixtures are close to the lean flammability limit.Furthermore, the goal is to examine the stochastic properties of schlieren and hydroxyl radical (OH * ) chemiluminescence, thus enabling to model the observed skewed probability density functions, and to propose closure relations for the corresponding moments.Note that OH * chemiluminescence could be used in practical devices to infer hydrogen-air flame properties, in particular those associated with the unstable behaviour.
On what follows, the experimental approach is first described.Then, the original experimental results obtained for a range of Damköhler and Reynolds numbers are discussed.The proposed stochastic model for the skewed probability density functions of the schlieren and OH * chemiluminescence is then developed, together with closure assumptions that describe the measured data.

Methodology
This section first describes the developed burner, the associated experimental methods and stochastic characterisation, followed by the relevant non-dimensional numbers definition.

Experimental methodology
Figure 1 shows a scheme of the the annular slit burner 'ANÔ' developed for this study.Note that this configuration bears some resemblance to one recently used to characterise high-frequency combustion instabilities [30].These instabilities are not found to arise in the present study, though.The slit outer and inner diameters are D o = 57.99mm and D i = 53.45mm, which yields a slit width of h = (D o − D i )/2 = 2.27 mm, and the slit length is L = 30 mm.The measurement uncertainty on these diameters is 0.2%.The central part of the burner is hollow so as to enable the free flow of ambient air.These dimensions have been chosen so as to enable the stabilisation of lean premixed flames close to the lean flammability limit, which is found to occur for a hydrogen-air mixture equivalence ratio close to φ = 0.26.Furthermore, the perimeter-to-slit width ratio, ≈ 75 , should enable unstable circumferential modes to evolve, but the development of unstable longitudinal modes is expected to be constrained by the smaller flame height-to-slit width ratio, which will be seen to lie between 1.5 and 4. The burner is made of stainless steel, and its uncooled surface temperature is monitored with thermocouples at the external and internal surfaces.During the experiments, the maximum measured temperature is always smaller than 350 K.
Premixed hydrogen-air enters the slit radially.Fuel-air mixture is effected in a 1.5 m length tube upstream to the burner.Hydrogen (99.999% purity) and air (compressed, oil and moisture free) flow rates are measured with mass flow controllers (Brooks models SLA5850 and SLA5851), which have a calibration uncertainty of 1%.As a consequence, the relative standard uncertainty associated with the hydrogen-air mixture equivalence ratio is 1.4%, whereas the Reynolds number uncertainty is 2%.The flow-meter ranges, i.e. 8, 20, 50 and 100 nlpm are chosen to cover Reynolds and Damköhler numbers values that lie between 100 and 2000, and 2.10 −3 and 100.10 −3 , respectively.These non-dimensional numbers are defined in Section 2.2.
Flame imaging is affected by combining schlieren and OH * chemiluminescence techniques.The latter uses a Princeton Instruments PI-MAX4:1024 EMB low repetition rate (10 Hz) electron multiplying back-illuminated CCD (emICCD, Gen II Intensifier, model RB-FG with P43 phosphor) camera equipped with an objective (CERCO 2085 94 F/4.1, 75% transmissivity in the UV range) and a 10 nm band-pass filter with 65% FWMH transmissivity centred at 311 nm.Note that the OH (A 2 + − X 2 ) transition from the electronic excited state is noted OH * , and that quantification of the measured results, albeit possible, has not been attempted [31].Since obtaining OH * images that encompass the whole annular flame is of interest, the measurement spatial resolution is 13.5 px/mm.This leads to the slit width being resolved along the direction parallel to the ICCD with 31 px.Due to the low OH * luminosity at the lean equivalence ratios of interest, the camera is operated at a gain of 6500, and a gate of 300 µ s is used.The raw, 1024 × 1024 px 2 , images are filtered using a symmetric 5 × 5 two-dimensional median filtering.Image statistics are obtained from a set of 500 images, which represent an averaging time of 0.5 s.The resulting minimum signal-to-noise ratio is 10, which corresponds to the flame obtained with the leanest mixture studied.An example of a mean OH * chemiluminescence field is shown in the left side of Figure 1.In this figure, both the Bunsen-like flame structure above the slit exit, and part of the annular OH * distribution can be seen.
A dual-lens schlieren arrangement is used to obtain time-resolved flame images.After traversing the flame, collimated light from a diode source is focused by a 100 mm diameter, 750 mm focal length, plano-convex lens to a 45 • knife edge.The deviated light is then focused by 95 mm diameter, 150 mm focal distance, bi-convex lens to the Photron FAST-CAM Mini AX200 CMOS detector.The schlieren images spatial resolution is 10 px/mm, as a consequence the burner annular slit is resolved with 23 px along the direction parallel to the camera sensor.Throughout the experiments, this camera is operated at 6400 frames per second, with a shutter speed of 4 µs.For each flow and mixture condition a set of 1000 schlieren images is obtained.The same median filter used to process the OH * images is applied to the 1024 × 1024 px 2 schlieren images.A typical mean schlieren image can be seen in the right side of Figure 1.The interfaces between the fresh unburned mixture and the combustion products, which lies above the slit, and between those and ambient air can be seen.
The set of instantaneous images is then used to obtain the statistical moments of order 1-3 of the experimental images, which should enable to infer the flame surface statistical behaviour.These moments are defined as: where µ and σ are the expected value and the standard deviation of the random variable Z, which has a density f (z), and μ 3 is the skewness.The skewness, that quantifies the departure of a density from a symmetric distribution, will be shown to be needed to fully describe the measured probability density functions.Note that this moment is scaled with respect to the standard deviation, and that the Gaussian distribution skewness is zero.In the present case where the variable Z represents the N sampled images, z i , Equation ( 1) is rewritten as: The different image processing and statistics algorithms used here are those implemented in Matlab.
Aiming to account for differential illumination and for the average flame structure, each schlieren image is background-corrected by dividing the local intensity by the mean.Thus, if S(x, y, t) is the instantaneous schlieren signal, and μ S (x, y) its mean, the stochastic behaviour of s(x, y, t), defined as S(x, y, t) = (s(x, y, t) + 1)μ S (x, y), ( 5 ) is of interest.Note that s(x, y, t) is a zero-mean quantity at each (x, y) position, and that the coordinate axis origin, given in Figure 1, lie at the burner surface along the symmetry axis.

Non-dimensional numbers
The studied flames are characterised using the Reynolds, Damköhler, Zel'dovich and Lewis numbers.The Reynolds number definition used is based on the mixture bulk velocity, Ū, and burner slit hydraulic diameter, D h , and may be written as: where ∀t is the total flow rate and ν is the mixture kinematic viscosity.The Damköhler number is defined as the ratio of a characteristic fluid time scale upstream to the annular flame, τ f , and the unstretched freely propagating planar flame time scale, τ c .The latter is obtained via Cantera [32] freely propagating one-dimensional adiabatic flame computations, being determined as the ratio between the flame thickness based on the maximum temperature gradient and the unstretched one-dimensional flame speed (δ L,max /S 0 L ).These flame computations used the Alekseev chemical kinetics mechanism [33] supplemented by a model for the OH * chemistry [31].Note that the computed flame speed is known to under-predict the existing experimental data for the lean mixtures of interest [33], which is also the case of other available chemical kinetics mechanisms [1].This under-prediction could be attributed to the inherent inability of the one-dimensional model to reproduce the measured flames cellular structure, which would increase the burning rate [34].Despite these shortcomings, these planar flame scales are used for the purpose of determining non-dimensional numbers.The following fits to the computed data are used: for which R 2 > .999 in the 0.28 < φ < 0.6 range.Within this range, the flame speed increases from 2.6 to 103 cm/s, whereas the flame thickness decreases from 2.4 to 0.2 mm as the chemical time scale decreases from 104 to 0.19 ms.Note that this flame model is used with the aim of providing reference time and length scales only.
The fluid time scale is taken as being the inverse of the flow velocity gradient upstream to the flame that, under the hypothesis of fully developed channel flow at the burner annulus (with h D i ), may be written as As a consequence, the Damköhler number reads as: In order to compare the one-dimensional unstretched laminar premixed flame thickness, δ 0 L , with the characteristic length scale of the flow velocity gradient, δ f , a non-dimensional thickness is used: Note that these non-dimensional quantities are related via: where Pr is the Prandt number, the flame thickness based on the maximum temperature gradient and the thermal flame thickness are δ L,max and δ L,T , respectively.Classically, the flame thermal thickness is δ L,T = α/S 0 L , where α is the mixture thermal diffusion coefficient.If the approximation δ L,max /δ L,T = 2(T b /T b ) 0.7 is used [35], Equation ( 11) may be written as: where the temperatures of the fresh and the burnt gases are T u and T b , respectively.Given this relationship, and since Da and δ * exhibit opposite trends with the hydrogen-air mixture equivalence ratio, the choice is made here to characterise the studied flames using Da and Re.Note that if the thermal and maximum temperature gradient flame thickness are supposed to be identical, the right hand side of Equation ( 12) reduces to 1/3.The hydrogen-air mixture effective Lewis number is computed via a formulation [36] that enables a continuous variation between fuel lean, and fuel-rich limits and uses the Zel'dovich to that end: where Le 2 = 2.33 and Le 1 = 0.33 are the fuel-rich and oxidiser-rich values for hydrogenair mixtures, and H is given by: m and n are the reaction orders with respect to the deficient and abundant reactants, respectively, here it is assumed m = n = 1, for simplicity.The normalised final concentration of the abundant reactant is either A = β(φ − 1)/Le 2 or A = β(φ − 1)/Le 1 , and is the Zel'dovich number.In order to determine the Zel'dovich number, β, a recently proposed methodology is employed [37], that permits to continuously compute β throughout the whole equivalence ratio range.This methodology optimally fits an ad-hoc premixed flame heat release rate expression to the values computed of the above-mentioned detailed chemical kinetics calculations.
In this expression θ = (T − T u )/(T b − T u ) is the temperature based reaction progress variable, and γ = (T b − T u )/T u is the flame thermal expansion, whereas β and β control, respectively, the initial and the final q decrease towards equilibrium as θ → 1.The resulting Zel'dovich number correlation that is used in Equation ( 14) reads: β = 5.432.10 4 exp(−30.22φ) + 6.166 exp(−2.657φ), (16) for which R 2 > .999 in the (0.28, 0.6) equivalence ratio range.The Zel'dovich number is found to decrease from ≈ 10 to ≈ 1.3 in this range, whereas the Lewis number increases from 0.38 to 1.0.Finally, it is worth mentioning that, given the above definitions, the ratio between the mixture bulk velocity at the slit and the laminar flame speed may be written as:

Results and discussion
In this section the studied flow conditions are first specified in terms of the hydrogen and air flow rates and the relevant non-dimensional parameters.The overall flame structure, as unveiled by the schlieren and OH * distributions is then presented, followed by an analysis of the corresponding stochastic properties.A model of the probability density function of these measured properties is then developed, together with the required closure hypothesis.Finally, the relation between this new model and the stochastic properties of the gradient of a reactive scalar is assessed.

Studied flow conditions
In this study are examined the results of 12 different flow conditions, given in Table 1 together with selected laminar one-dimensional unstretched flame properties and the corresponding non-dimensional numbers defined in Section 2.2.These flow and mixture studied have been chosen to span a Reynolds number range of (400, 1000) and a Damköhler number range from 2.10 −3 to 16.10 −3 , which correspond to lean fuel/air mixtures close to the lean flammability limit.In this table a given letter corresponds to a Reynolds number, whereas the integers indicate the Damköhler number.The burner is operated in 1 atm and the combustible mixture fed at 294 K.The values of δ * reported in this table indicate that the laminar flame thickness (δ L,max ) may be either smaller or larger than the slit half-width.Indeed, for the leanest and richest fuel/air mixtures studied, the computed flame thickness is 2.8 and 0.76 mm, respectively, which correspond to δ * of 2.5 and 0.67.Furthermore, the leanest mixture Zel'dovich number is of the order of 14, whereas for the richest mixture β < 4. Note that this leads to a significant change of the computed chemical heat release rate distribution within the premixed flame, not shown here for the sake of brevity.Indeed, for the smaller values of β, the heat is released for intermediate temperature values, whereas for the larger values this heat release occurs closer to the maximum flame temperature.
As shown in Table 1, the effective Lewis number of the studied cases are all significantly smaller than one.Therefore, thermo-diffusive instability is expected to influence the flame development above the annular slit.
These studied flow conditions lie closer to the burner blow-off limit than to the flashback one.Blow-off is considered to occur when the flame anchoring ceases at either of the annular surfaces -internal or external, whereas the flashback is annotated once the flame fully enters the annular space, since partial entrances are not visible.Note that the schlieren images are used to characterise these limit situations.The blow-off limit is nearly independent of the Reynolds number, being observed for φ ≈ 0.27, and thus Da ≈ 1.5.10 −3 .The studied conditions are such that blow-off is not observed, and the flame wrinkling will be shown to be purely combustion driven.For the considered Reynolds number range the critical flashback Damköhler number is of the order of 80.10 −3 .

Annular flame structure
Typical instantaneous, average and background-corrected schlieren images are shown in Figure 2(a ) and 2(b).The corresponding flow conditions, given in Table 1, were chosen to represent the smallest (A1) and largest (D3) Damköhler and Reynolds number values studied.A first assessment of the annular flame structure may be obtained by examining the instantaneous schlieren images, S(x, y, t), shown at the top images.At the left and right parts of each image the Bunsen-like conical flame surface is clearly seen above the burner surface, as is the interface between the burned gases and ambient air that surrounds this flame.The centre part of these images is occupied by collected light intensity fluctuations which could stem from a wrinkled flame surface.In the average images, μ S (x, y), given at the centre of these figures, these wrinkled structures are absent for case D3, thus confirming their instantaneous nature.Concerning case A1, a reminiscence of the wrinkles is seen, one which is associated with the use of a fixed sampling time for all experiments, and to the the slower fluctuating motion of this case, when compared to D3.The flame cone and interface between the burned gases and ambient air can be clearly seen also.In these average images, the flame cone tip seems to be open, which is a classical feature of low Lewis number flames, such as those.The bottom images of these figures display the instantaneous background-corrected signal, s(x, y, t), defined in Equation (5).It is possible to verify that the background correction eliminates the conical flame structure above the slit, as well as the interface between the burned gases and ambient air.Furthermore, these wrinkles seem to be vertically aligned close to the burner surface, displaying a seemingly two-dimensional structure along the flame, but exhibit a more complex structure further downstream, in particular, for the low-Da, high-Re situations.Downstream the flame cone these wrinkles progressively disappear, which suggests that they are combustion-driven.The quantitative analysis of the flame wrinkling length scale is deferred to Section 3.3.
To further characterise the annular flow structure, Figure 3 shows the average OH * signal, μ OH , corresponding to cases A1 (top) and D3 (bottom), given in Table 1, and which schlieren results are depicted in Figure 2(a ) and 2(b), respectively.Note that the colour scales used in these images are different, since the maximum value of μ OH in case D3 is approximately five times larger than that of case A1, which is explained by the higher equivalence ratio of the former mixture.In this figure the flame cone can clearly be seen at  the left and right sides, and the cone tip opening that has been hypothesised to exist in the schlieren averaged images analysis is confirmed.Since these are line-of-sight integrated images the OH * intensity is larger at the part of the cone closer to the burner symmetry axis than at the outer parts.Along the centre part of these figures the OH * signal is fainter, which is also a feature of this line-of-sight measurement.Due to the low signal-tonoise ratio of the smaller equivalence ratio cases, no attempts were made to perform image deconvolution, that would have enabled us to describe the OH * signals along the burner symmetry plane.Further comparing the OH * (Figure 3) and schlieren (Figure 2(a ) and 2(b)) distributions, it is clear that the wrinkling is associated with the combustion process, i.e. to the flame surface.Therefore, the analysis of the schlieren and OH * stochastic behaviour should provide a description of these flame front properties.

Flame stochastic behaviour
It may be expected for the statistical properties of the schlieren and OH * images to change with Da and Re.Accordingly, the standard deviation and the skewness of the backgroundcorrected schlieren images, σ S and μ 3,S , given in Figure 4 as a function of the vertical coordinate, enables to quantify the flame surface fluctuations.In this figure are compared, for a given Damköhler number, the results obtained for each Reynolds number.These results correspond to the averaging, along the direction parallel to the burner surface, x, of a given statistical moment computed for the set of 1000 images, i.e.

μi,S (y)
where i = 2 corresponds to the standard deviation, and i = 3 to the skewness.Since only μi,S (y) are analysed henceforth, the bar is omitted in the remaining of the text.In order to restrict the focus to the central portion of the annular flame, this averaging is performed along x = 40 mm only (400 px).The sample size is thus 40,000 at a given vertical position above the burner surface.
Figure 4 shows that a similar behaviour is observed in all cases, i.e. the standard deviation first increases, then decreases above the burner surface.The non-zero value of the standard deviation at the burner surface (σ S (y = 0) = σ S,0 ≈ 0.09) and initial increase along the vertical position seem to be independent both of Da and of Re.However, the maximum standard deviation value and the corresponding vertical position are found to increase with Re, being less sensitive to Da changes.In particular, for the smaller Da values, when the Reynolds number increases from 400 to 1000, the maximum standard deviation goes from σ S,max ≈ 0.13 to ≈ 0.25.Table 2 gives these maximum values for the different studied cases.Furthermore, this non-monotonic standard deviation behaviour with respect to y suggests that this maximum position could be used as a measure of the flame height, which will be examined below.
In order to expand the statistical analysis of the flame front fluctuations, Figure 4 also depicts the corresponding skewness.This statistical moment also exhibits a non-monotonic behaviour above the burner surface, first increasing, then decreasing with height.Furthermore, similarly to the standard deviation, the maximum skewness value is found to increase with the Reynolds number, but is nearly independent of Da.In particular, the skewness increases from ≈ 0.5, for Re = 400 to ≈ 1.0, when Re = 1000.These non-zero values of the skewness indicate that the probability density function of s should be an asymmetrical distribution, which will be examined below.
To further the analysis, the transversely averaged first and second moments of the OH * images, i.e. μ OH and σ OH , are given in Figure 5.The choice of depicting the average OH * evolution is motivated by the classical use of chemiluminescence to characterise combustion heat release, however no attempts to quantify the concentration has been made [31].In this figure, these moments have been arbitrarily divided by the global mean maximum, that is found to occur for case D3, which thus enables to quantitatively compare the cases among themselves.Similarly to what has been observed above for the schlieren moments, both μ OH and σ OH increase with the Reynolds and the Damköhler numbers, the latter seems however more pronounced here than for σ S .The overall evolution of μ OH and σ OH with the vertical position, y, is quite similar.Indeed, a increase near the burner surface which is steeper than that of σ S is first seen.Further downstream the moments decrease to zero, which is consistent with the tendency of reaching chemical equilibrium downstream the flame, but a different behaviour is found to occur for intermediate y values.Indeed, concerning μ OH , a slow decrease is seen to occur, the extent of which increases with Re but decreases with Da.On the contrary, in this intermediate y range, σ OH slowly increases.This slow increase could be a manifestation of the increased schlieren fluctuations, since, by comparing Figures 4 and 5, the decreasing behaviour is seen to initiate immediately downstream the maximum σ S position.Therefore, the y position at which σ S,max occurs, y σ S , will be used as a measure of flame height.The y σ S values obtained for the different studied cases are also given in Table 2.These values may be fitted by with a y = 0.155 mm for a R 2 = 0.96 with 95% confidence bounds.
Another relevant feature of the wrinkled flame front is the characteristic length scale of these wrinkles.The average length scale of the schlieren fluctuating signal seen in Figure 2(a,b), l S , is given in Figure 6 as a function of the non-dimensional vertical position, y/y σ S .This length scale is defined as where δ(x) is the Dirac delta function, and s b (x, y, t i ) is the binarised schlieren signal.This enables to determine the number of positive and negative s(x, y, t i ) regions for each y.The sum is performed throughout the N images, and x is the width of the considered region.More precisely, this length scale has been determined by the following procedure.First each instantaneous image is binarised, thus enabling to identify regions with opposite s signals to finite regions of zeros or ones.Then, at each vertical position, y, a given binarised image is treated as a square wave, in which pulse widths are determined, then averaged at a given height for all images, thus leading to l S .Given this procedure, Figure 6 shows that the value of l S remains nearly constant for 0 < y/y σ S < 1.5, and increases further downstream.The values of l S within this plateau region are independent of the Reynolds number and decrease with the Damköhler number.Indeed, as Da increases from ≈ 2.10 −3 to ≈ 15.10 −3 , l S goes from 1.3 to 1.0 mm.This suggests that the processes that control l S remain invariant in that y/y σ S range, which should be the case in regions of the flow where s is controlled by the premixed flame, i.e. by the balance between chemical and fluid time scales.Furthermore, an increase of l S associated with viscous damping of the fluctuations was to be expected downstream the flame, i.e. for y/y σ S > 1.5 since, as shown above, σ S decreases in this region.These results show that the initial flame wrinkling is purely combustion driven, and that downstream the flame tip viscous dissipation acts to dampen these wrinkles.This is indeed the case of the flame structure results depicted in Figure 2, where both the instantaneous, S(x, y, t), and the fluctuating, s(x, y, t) schlieren exhibit quite regular structures along the transversal direction close to the burner surface (y σ S < 1).Downstream to the average flame tip (y σ S > 1), these structures are progressively dampened.
To further characterise the fluctuating behaviour of the schlieren images, Figure 7(a) and 7(b) presents the histograms (in blue) of the probability density function of s, P(s), at different heights, for cases A1 and D3, respectively.The modelled PDF, given by the red line, will be discussed in Section 3.4.These heights have been chosen to correspond to y/y σ S values of (0.5, 1.0, 1.5 and 2.0) with the purpose of illustrating the P(s) evolution around the maximum σ S region.These similar normalised heights correspond to different y values, which are indicated above each PDF.Following what has been done above for the stochastic moments, this PDF has been constructed by accumulating the s values in a transverse region spanning −20 ≤ x ≤ 20 mm.Aiming to ease the comparison between each case, that correspond to rather different standard deviation values, the horizontal axis has been normalised by the corresponding maximum value, σ S,max , which are given in Table 2.Both these figures show that P(s) spans to smaller s/σ S,max away from y/y σ S = 1.0, i.e. as the value of σ S decreases upstream or downstream this maximum position.However, these cases correspond to rather different P(s) shapes.Indeed, the case A1 exhibits a nearly symmetrical PDF with respect to s = 0, whereas the PDF of case D3 is A model for this probability density function, P(s), will now be developed.However, the discussion of the OH * probability density function, which also exhibits a significantly skewed shape, is deferred to after the presentation of different P(s) moments models.

Proposed stochastic model
A stochastic model of the schlieren fluctuating signal s(x, y, t) is now developed.This development proceeds by considering first the corresponding probability density function, P(s), then a relation between the parameters of this PDF and finally, an evolution equation for these parameters.Relations between the obtained expressions and those that ensue from the transport of a reactive scalar function are discussed also.

Prescribed probability density function
Given the nature of the observed stochastic process of s(x, y, t), which involves a decay downstream to a region where forcing, F(s), occurs, it is proposed that at each height, y = y/y σ S , the PDF of s, P(s), is described by: In the absence of forcing, i.e.F(s) ≡ 0, it is evident that P(s) is an exponential PDF with decay rate τ P .When forcing is present and τ P → 0, one obtains P(s) → F(s).Here it is assumed that the forcing is a Gaussian PDF with zero mean and standard deviation σ P , i.e.
This choice leads to a closed-form solution of Equation ( 21) as an exponentially modified Gaussian function: where erfc is the complementary error function.This PDF has been used to describe the shape of chromatographic measurements, for instance [38,39].The PDF parameters μ P , τ P and σ P are related to the moments of P(s) following , , where the sample mean, standard deviation and skewness (μ 1 , μ 2 and μ 3 ) are given by Equations ( 1)-( 3).Accordingly, provided these moments are known, the PDF is completely determined.Note that the definition of s trivially leads to μ 1 ≡ 0, thus two moments are left to be specified.
A first obvious choice is to use the measured values of μ 2 and μ 3 depicted in Figure 4 to determine the PDF parameters in Equations (24).The resulting P(s) is given in Figure 7(a ) and 7(b) as the solid red lines.The agreement between the measured PDF and the prescribed exponentially modified function is excellent and applies also to all the studied cases and for y values different than those depicted.This excellent agreement permits to say that this PDF model fully explains the measured fluctuating schlieren data.However, the closure of the PDF transport does require that the second-and third-order PDF moments be specified, which is addressed now.

PDF parameters closure
In order to provide closed expressions for the exponentially modified PDF, parameters σ P and τ P , given in Equation (24), are assumed to be functions of the Reynolds and Damköhler numbers, defined in Equations ( 6) and ( 9), respectively.
A first modelling hypothesis is that τ P describes P(s) relaxation and, thus, could be supposed to be dependent on the Reynolds number only.Since this relaxation should reflect the balance between inertia and viscous forces, it is then assumed that τ P ∝ Re.A second modelling hypothesis is that the PDF forcing, supposed to follow a normal distribution (Equation ( 22)), includes both chemistry and transport effects, such that σ P = σ P (Re, Da) is a function to be determined.
However, before developing these closure functions, it is interesting to note that σ P and τ P exhibit a peculiar relationship.This can indeed be observed in Figure 8, where the evolution of τ P /σ 2 P , as determined by using μ 2,S and μ 3,S in Equation ( 24) is given as a function of the non-dimensional vertical coordinate, y/y σ S .This figure permits to verify that a nearly constant value of τ P /σ 2 P arises, for each (Re, Da), in the range 0.25 < y/y σ S < Figure 8. Ratio between the exponentially modified PDF parameters, τ P /σ 2 P , a function of the height above the burner surface, y/y σ S , for the different studied cases.(Colour online.)2.0, i.e. along the flame region.A rather clear separation arises among the different case studies, which should enhance the corresponding discrimination in the (Re, Da) plane, also.Furthermore, this behaviour suggests that the relation between the parameters values at the particular position where the PDF variance is maximum, i.e. y/y σ S = 1, could be used to represent the same relation elsewhere, thus enabling to express τ P as a function of σ P .
Accordingly, these parameters at the maximum standard deviation position are expressed as functions of (Re, Da), where τ P /σ 2 P | max is the value at y/y σ S = 1, and the parameters a τ , a τ ,σ , b τ ,σ and c τ ,σ , obtained by fitting the measured data given in Table 2 are with 95 % confidence bounds.This first closure for the PDF parameters assumes that a relation exists between τ P and σ P which is a known function of the Reynolds and Damköhler numbers.As a consequence, three parameter exponentially modified Gaussian PDF is reduced to a single parameter function, σ P .Thus, reworking Equations (24), one may express this parameter as a function of μ 2,S , the measured PDF standard deviation: Note that the results of Figure 8 show that the order of magnitude of f τ ,σ is ≈ O (10), which implies that is a very good approximation for most of the y/y σ S range, except when y/y σ S → 0.
A final closure step involves determining an expression for σ P .Accordingly, this parameter normalised by the corresponding maximum and minimum values, i.e. σ P = (σ P − σ P,0 )/(σ P,max − σ P,0 ) is shown in Figure 9 as a function of the non dimensional vertical coordinate, y/y σ S for the different studied cases.It is possible to observe that all measured results collapse in the range 0 ≤ y/y σ S ≤ 1.5, exhibiting discrepancies for larger y/y σ S values.The observed behaviour is that of an initial increase, followed by a saturation and, then, a progressive relaxation.This behaviour could be described by prescribing σ P as a function of the non-dimensional coordinate, y = y/y σ S , such as: where n is a parameter that controls the σ P (y ) relaxation as y increases.This proposed expression for the vertical evolution of the PDF forcing standard deviation corresponds to a variable coefficient ordinary differential equation, i.e.
The relation between this equation and a reactive scalar transport equation will be explored in Section 3.5.
The results of such model for two arbitrary values of n (2 and 4) are also given in Figure 9.These values are found to yield to σ P (y ) that encompass the measured values, thus leading to a full closure if σ P,max and σ P,0 are known.More precisely, for y > 1, the smallest and largest Reynolds number σ P (y ) curves are bounded by Equation ( 29) curves for n = 4 and 2, respectively.This inverse behaviour suggests a simple ad-hoc model for the relationship between the Reynolds number and this decay parameter, where the value C n = 1600 has been used to close this relation.The maximum value, σ P,max , is taken to be the function of (Re, Da) determined by relating Equations ( 25) and ( 26), The comparison between the proposed closure models and the measured PDF is given in Figure 10(a ) and 10(b).The solid red lines correspond to P(s) computed with the experimental values of μ 2,S and μ 3,S , which were shown in Section 3.3 to reproduce the experimental histograms exactly and are thus identified here to the experimental PDF.Two models are considered.The first, Mod1, assumes the relationship between τ P and σ P given by Equation (26), thus using the measured σ P together with Equation (29) to describe P(s).The second model, Mod2, further assumes that σ P (y ) is given by Equation (31), thus requiring the values of y σ S , σ P,max , and σ P,0 also.In these figures one may verify that the agreement between the first model and the experimental PDF is excellent.Such an agreement is also observed for the different cases studied, which leads to the conclusion that the closure given by Equation ( 26) and the knowledge of the standard deviation of the forcing function (σ P ) suffice to determine the PDF.Concerning now the second model, where σ P (y ) is prescribed, the agreement with the experimental PDF is also very good for case A1, and for case D3 close to y = 1.However, for this latter case, discrepancies are evidence both close to the burned surface and further downstream.In particular, the modelled standard deviation is larger than the experimental one, which leads to wider P(s).Nevertheless, the model adequately describes the overall variation of the PDF shape for all the studied experimental conditions, which is not shown here for the sake of brevity.
Given the excellent agreement between the prescribed PDF and the schlieren experimental data, and that a this PDF requires a model for the PDF forcing standard deviation, the relation between this parameter and flow-field transported properties is now examined.

Relation to the reactive scalar transport
A relevant analysis point concerns the relation between the schlieren PDF forcing parameter, σ P , and the standard deviation of a transported reactive scalar gradient.In premixed combustion, the temperature is often described by a reactive scalar (0 where ρ, T, v, D and ω(θ) are the mixture density, temperature, velocity, diffusion coefficient and chemical source term, respectively, whereas the initial and adiabatic flame temperature are T u and T b .This equation is obtained from the energy transport equation under the classical hypothesis of constant thermodynamic pressure, low Mach number, heat transfer due to temperature gradients only, unity Lewis number for all chemical species, and constant ρD.
To mimic the schlieren separation in averaged and instantaneous components (Equation ( 5)), it is assumed that the progress variable gradient may be written as ∇θ = (1 + ϑ)∇μ θ , where ϑ is a scalar perturbation to the average progress variable gradient.
Schlieren images result from the integration along the line-of-sight (direction z here) through the refraction index field, m(x, y, z, t), that is inversely proportional to the temperature, followed by the corresponding gradient component selection by the knife edge on the focal plane [40].The schlieren increment may thus be written as: where use has been made of the temperature gradient decomposition, and of the Taylor series expansion of the progress variable.The comparison of Equations ( 5) and (36) provides an analogy between the schlieren and temperature gradient fluctuations, s and ϑ respectively.Even if this analogy is contingent on the used hypothesis, it enables to identify the processes controlling the fluctuations variations in the y direction, and thus of the schlieren standard deviation.Accordingly, Equation ( 35) may be rewritten as an instantaneous equation for ϑ that, after multiplication by ϑ and averaging, leads to an equation for μ ϑ 2 , the variance of ϑ: where the average source term is noted μ ρ ω, and μ ρ ωϑ is the correlation between the chemical source term and ϑ.This equation shows that the μ ϑ 2 production is controlled by this correlation, whereas the term 1 2 ρD μ ρ ω ∇μ θ provides characteristic length scales.To compare with the model of the preceding section, Equation ( 37) is integrated in the (x, z) directions, where represents this spatial integration.Returning now to the ad-hoc model developed in Section 3.4.2, the combination of Equations ( 30) and ( 32) leads to a model equation for the variance of the schlieren signal which -given the analogy between s and ϑ (Equations ( 5) and ( 36)) -suggests a closure expression to Equation (38), i.e.
As shown in Figure 11, this hypothesis is qualitatively supported by the vertical evolution of the average schlieren signal, at least in the range 0.5 ≤ y ≤ 1.25.The discrepancies that are observed below this range are due to the light beams that are deviated by the flame toward the burner surface, which effectively decreases the incident light level.The large variations among the different cases that are seen in this figure for y > 1.5 arise in the burned gases region, as discussed in Section 3.3 in connection to Figures 4 and  5. Accounting for these variations is thus beyond the possibility of the variance model proposed.
The right hand side of Equation (38) represents the correlation factor between ρ ω and ϑ, which extremum corresponds to the maximum σ S , and that should decay to zero in the burned gases region faster than μ ρ ω does.At this stage, a complete closure for this term should require assumptions regarding the velocity field, which has not been characterised.However, the simpler approach adopted, which is limited to the present burner configuration, is to write that The closure of Equation ( 38) is effected by prescribing f y (y ) ∝ (y ) n , and by recognising that σ 2 S = μ ϑ 2 .Given Equation ( 30) and the PDF model (Equation ( 23)), it is clear that this correlation between the source term and the fluctuating schlieren signal is directly related to the PDF forcing, σ P .Improving the proposed ad-hoc model for this parameter could be done by learning from experimental data, for instance, which will be subject of future work.Finally, since the maximum variance value corresponds to y = 1, Equations ( 30), (34) and (41) imply that

OH * stochastic behaviour
Under the premise that the schlieren signal is closely related to temperature gradients, and that temperature is a controlling property in flames, it is interesting to examine the limitations of the prescribed PDF model for the description of the OH * fluctuations.Indeed, in the limit of very lean flames, it has been shown that all radicals follow a quasi-steady-state approximation [34,41].However, to the best of the author's knowledge, the relation between the schlieren and OH * has not been explored.Accordingly, Figure 12 depicts the computed maximum values of ∇T and ∇T/T as a function of the OH * molar fraction maximum, (X OH * ) max .These values have been obtained considering a steady adiabatic freely propagating planar one-dimensional flame for different equivalence ratio of the fresh unburned gases, φ ∈ (0.3, 0.6) -the smaller equivalence ratios corresponding to smaller X OH * values.These temperature gradient expressions have been chosen because the schlieren can be shown to be proportional to ∇m/m, where m is the refraction index of the medium, that is inversely proportional to the temperatures for gases [40].This figure clearly shows that an increasing relationship exists between (∇T) max , (∇T/T) max and (X OH * ) max .More specifically, (∇T/T) max ∝ (X OH * ) 1/3 max , this 1/3 exponent corresponding to the depicted straight line in this log-log representation.Furthermore, for the equivalence ratios considered in the experimental part of this work, which correspond to X OH * ∈ [5,200] • 10 −14 , (∇T) max ∝ (∇T/T) max is a good approximation.Note, also, that the two order of magnitude variation of OH * molar fraction in the studied range underscores the measurement difficulty associated with this species in these lean hydrogen-air flames.As a consequence, in this limit, it could be expected for the OH * stochastic behaviour to closely follow that of the temperature gradient.To verify this, Figure 13(a ) and 13(b) show the OH * fluctuations probability density function for cases A1 and D3, respectively, at different equally spaced heights.Following what has been done for the schlieren PDF in Section 3.3 to enable the comparison between these cases, the OH * fluctuation is divided by the corresponding standard deviation.In these figures the histograms (blue) correspond to the experimental PDF, whereas the curve (red) is the modelled result obtained by combining the measured standard deviation, σ OH * (y), given in Figure 5, with Equations ( 24)- (26).Both these figures show that the measured OH * PDF is significantly skewed, however the agreement between model and experiments is quite different.Indeed, the first of these figures, which correspond to an equivalence ratio of 0.28, show that this agreement is excellent, whereas in the second -obtained for φ = 0.35 -discrepancies exist.More precisely, given σ OH * (y), the OH * PDF skewness is fairly well described by the model relating the Reynolds and Damköhler numbers to the PDF forcing for the smaller equivalence ratios, whereas modelled skewness values that are larger than those measured arise for higher equivalence ratios.Such a behaviour indicates that, even if the OH * fluctuations PDF seems to be described by Equation (21), the functional dependency of the forcing and relaxation parameters on the relevant non-dimensional numbers should be reexamined.The reasons for such a behaviour change are unclear at present, but are subject of an on-going work, where computational fluid dynamics model results of the studied flames are expected to complement this experimental study.Note that, as the equivalence ratio increases from 0.28 (case A1) to 0.35 (case D3), the chemical pathways that control the combustion process have been shown to change [41].More precisely, the chemical equilibrium hypothesis could hold for the OH * radical in the leanest mixtures, but not for the richer ones.Therefore, in the leaner cases the OH * fluctuations would be controlled by the temperature fluctuations, whereas a more complex dependency on other chemical species should arise for richer mixtures.Furthermore, the characteristic length scales of OH * and ∇T/T could be different as well, and could be at the origin of the observed discrepancies.

Conclusions
Hydrogen-air lean premixed laminar flames stabilised on an annular slit burner have been studied.The experimental conditions spanned over a range of Reynolds and Damköhler numbers close to the lean flammability limit.The studied flames, which were characterised via OH * chemiluminescence and time-resolved schlieren, were shown to exhibit unstable flame surface fluctuations along the direction transversal to the flow-field.The analysis of the transversely averaged statistical moments of these schlieren images showed that both the schlieren variance and the skewness first increased with the height above the burner surface, then decreased downstream the flame tip.A similar behaviour has also been evidenced for the OH * mean and variance.
The probability density function of the schlieren fluctuations was shown to be described by an exponentially modified Gaussian function, which is a skewed distribution that models the considered property decay when subject to a Gaussian forcing.Furthermore, these distribution parameters, which are linked to the measured statistical moments, were found to exhibit a peculiar relation, which is a function of the Reynolds and the Damköhler numbers.This relation, together with a hypothesis on the decay dependence on the Reynolds number, enabled to describe the schlieren fluctuations using the Gaussian forcing parameter only.To completely close the stochastic model an ad-hoc, flow-field dependent, expression for this forcing was proposed.The connection between the corresponding equation and a reactive scalar transport equation enabled to relate the proposed closure to the scalar transport terms.
Finally, the extent to which the OH * chemiluminescence probability density function may be described by the proposed model was examined, showing that the model is adequate for the smaller Reynolds and Damköhler numbers values, but should be modified for larger ones.These results are expected to pave the way towards the development of hybrid reduced-order models, in which the unmeasured flow-field dynamics may be learned from the experimental data.
V. Montassier and A. Claverie and to L. Pirateque Henao and N. Lopes M. B. Junqueira for performing the premixed flame computations.This work pertains to the French government programme 'Investissements d'Avenir' (EUR INTREE, reference ANR-18-EURE-0010).For the purpose of Open Access, a CC-BY public copyright license has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission.

Figure 1 .
Figure 1.Annular slit burner 'ANÔ' scheme, with examples of mean OH * chemiluminescence (left) and schlieren (right) images obtained for Reynolds and Damköhler number values of 1000 and 16.10 −3 .The dashed rectangle (40 × 27.5 mm 2 ) represents the area used for averaging the instantaneous schlieren results, whereas the dash-dot line is the burner symmetry axis.Slit outer and inner diameters: D o = 57.99mm and D i = 53.45mm.(Colour online.)

Figure 3 .
Figure 3. Average OH * images (μ OH (x, y)): case A1 (top) and case D3 (bottom) on a 64 × 23.5 mm 2 field of view above the burner surface.The maximum OH * value of case D3 is 5.2× that of case A1. (Colour online.)

Figure 5 .
Figure 5. Vertical evolution of the OH * mean deviation (μ OH , top) and standard deviation (σ OH , bottom) for the different studied cases.(Colour online.)

Figure 9 .
Figure 9. Exponentially modified PDF parameter σ P a function of the height above the burner surface, y/y σ S , for the different studied cases.Also show are Equation (31) model results for two values of n. (Colour online.)

Figure 10 .
Figure 10.Comparison between experiment and models of the probability density function of the schlieren fluctuations s(x, y, t)/σ S,max at different heights above the burner surface, y/y σ S = (0.5, 1, 1.5, 2), the corresponding y values are given above each image.(a) Case A1 and (b) Case D3. (Colour online.)

Figure 11 .
Figure 11.Mean schlieren signal as a function of the height above the burner surface, y/y σ S , for the different studied cases.(Colour online.)

Figure 12 .
Figure 12.Maximum values of ∇T (circles) and ∇T/T (squares) as a function of (X OH * ) max computed for a steady adiabatic freely propagating planar one-dimensional laminar flame.The straight line corresponds to a 1/3 derivative in this log-log representation.(Colour online.)

Figure 13 .
Figure 13.Comparison between experiment (blue histogram) and model (red curve) of the probability density function of the OH * at different equally spaced heights above the burner surface, the corresponding y values are given above each image.(a) Case A1 and (b) Case D3. (Colour online.)

Table 1 .
Experimental flow conditions, laminar flame properties and non-dimensional numbers.