## Unveiling the bi-stable character of stealthy hydrogen-air flames

Rubén Palomeque-Santiago, 1 Alba Domínguez-González, 2 Daniel Martínez-Ruiz, 2 Mariano Rubio-Rubio, 1 Eduardo Fernández-Tarrazo, 1 and Mario Sánchez-Sanz ∗ 1

1) Dpto. de Ing.

Térmica y de Fluidos, Universidad Carlos III de Madrid, 28911, Leganés, Madrid,

España

2)

ETSIAE., Universidad Politécnica de Madrid, Plaza del Cardenal Cisneros 3, 28040, Madrid,

España

(*Electronic mail: mssanz@ing.uc3m.es)

(Dated: 5 August 2024)

Ultra-lean hydrogen-air flames propagating in narrow gaps, under the influence of cold walls and high preferential diffusion, can form two distinct isolated structures. They exhibit either circular or double-cell shapes and propagate at different speeds, with the latter roughly doubling in size and traveling speed to the former. Hydrogen mass diffusivity, convective effects and conductive heat losses are the physical mechanisms that explain the alterations in morphology and propagation speed. In previous experiments, Veiga et al. Phys. Rev. Lett. 124, 174501 (2020) found these clearly distinguished flame structures for different combinations of equivalence ratio, channel gap and the effect of gravity on the dynamics of upwards and downwards propagating flames in a vertical chamber. Present observations in horizontal channels show the simultaneous appearance of these two stable structures, which arise under identical experimental conditions and conform the first evidence that multiplicity of stable solutions coexists in real devices. To explain the observations, we performed numerical simulations using the simplified model, which show that symmetry-breaking details during the ignition transient explain the concurrent emergence of the two stable configurations. This discovery urges the need to implement additional engineering tools to account for the possibility of formation and propagation of isolated flame kernels at different speeds in hydrogen-fueled systems.

## I. INTRODUCTION

Molecular hydrogen (H2) is recognized as a crucial energy carrier for achieving Net Zero greenhouse emission goals in the coming decades. The anticipated increase in its utilization in the power generation and transport sectors has raised legitimate concerns about associated risks and threats 1,2 . Understanding the initiation and progression of chemical reactions in gaseous mixtures is vital for advancing safe hydrogenpowered devices. Among the possible scenarios, we consider here the leaking of hydrogen from a container that, after mixing with the surrounding air, can form a reactive mixture. If this mixture comes into contact with an energy source, it might ignite to form a potentially dangerous traveling premixed flame. The narrow channel configuration addressed here is a confined space that serves as representation of fuel leakage accidents in caged devices. The growing use of hydrogen in the future and its singular properties as a fuel (wide flammability limit, low ignition energy and almost invisible flame) will certainly increase the risks during its utilization.

steady two-dimensional counterpart had been previously considered nonexistent in theoretical analyses 5 , and they have only been found as drifting flame cells in two-dimensional simulations 6 .

The thermo-chemical characterization of this type of flames is typically addressed from the amenable and simplified planar flame (1D) theory 3 , a useful starting point but insufficient to explain experimental observations. These deviations from the theoretical framework are particularly pronounced as the H2 concentration in air decreases. In fact, the planar flame construct does not exist for the ultra-lean conditions studied here, when preferential diffusion effects take part in enabling the survival of curved reactive fronts that display locally richer and hotter conditions.

Stable multidimensional isolated flame structures, such as three-dimensional flame balls, have been predicted theoretically and observed experimentally 4 . On the other hand, their

However, Veiga-López et al. 7 found that various types of nearly two-dimensional stable isolated-flame configurations can be found under confinement due to geometrical constraints and additional viscous and conductive processes. They found multiple quasi two-dimensional propagating fronts in lean hydrogen-air mixtures in a vertical HeleShaw cell. Depending on the controlling parameters, channel width and mixture composition, the flames remained stable propagating at constant speed or divided forming a fractal-like path. These isolated flame fronts, similar to flame balls and drifting flames, are characterized by their small size and low emissivity, which makes them very elusive. Among the atlas of behaviors, only two stable cases were detected. Based on the equivalence ratio of the mixture and the flame propagation direction (upwards or downwards) either circular or doublecell isolated traveling flames were produced, but never observed simultaneously. Circular flames burn slowly and were observed during upward propagation. Double-cell structures travel faster and were only observed during downward propagation. It was thus concluded that buoyancy effects would determine the prevalence of one type over the other by either accelerating or decelerating the front owing to the placement of hot products. This result confirmed the conclusions of previous numerical studies that included buoyancy effects 8 .

An exhaustive parametric study by Domínguez-González et al. 9 detailed the occurrence of numerically stable solutions in the nearly 2D limit with heat losses. Therein, it was shown that both configurations, circular and double, can arise in a common region of the parametric space subject to different

FIG. 1. Schematic of the horizontal chamber and experimental system setup used for image acquisition. The dimensions of the gap are L × W × h , with L = 950 mm, W = 220 mm and 2 &lt; h &lt; 6 mm.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00893v1-figure-6dd87467-93dc-4a5b-b0e5-bb20d93a6f8b.jpg)

initial conditions. It is remarkable that given a prescribed mixture and vessel properties, two stable flames, of different sizes and propagation speeds, can be formed.

To extend these previous studies 7,9 , we experimentally address the formation of both isolated flame structures in real devices. To do so, we designed and constructed a new horizontal experimental set-up to eliminate the alignment of buoyancy effects with the propagation direction. This new experimental campaign shows that the two stable flame structures can coexist under identical geometric and mixture conditions. Such simultaneous coexistence of two stable canonical configurations, which has not been observed before, constitutes the core of this manuscript and represents a unique case in the literature, to the best of the authors' knowledge.

sions to the constrained length h , perpendicular to the propagation direction. The narrow gap between the plates results in a reduced combustion volume, creating a competition between heat generation at the flame and conductive heat losses through the walls. These cases are of wide interest in combustion science 10,11 , being typically easier to control, measure and analyze than 3D experiments of premixed flames 12 in absence of wall effects.

The paper is organized as follows. We first offer a description of the configuration at hand and define the characteristic scales of the problem. The detailed presentation of the experimental setup and data processing methods comes next. Along the paper we identify the parameters that control the onset of bistability in these flames to compare our results with recent numerical studies. Finally, we discuss the experimental and numerical results before presenting our conclusions.

## II. NEARLY TWO-DIMENSIONAL CONFIGURATION

We shall consider in the following a horizontal Hele-Shaw chamber, formed by two parallel plates that are separated a small distance apart h , in the millimetre range, as illustrated in figure 1. This particular configuration serves as a representative model for confined geometries and enables the simplification of the analysis by reducing one of the spatial dimen-

Although extinction and quenching can be anticipated for an extreme reduction of the channel gap ( h ≲ 1 mm), it is possible to find stationary flames in the combined regime of ultra-lean hydrogen-air mixtures with a slightly larger plate separation 7 . Considering that the isolated cell flame is of characteristic size δ i and it propagates with velocity ULi , we can estimate of the heat flux from the flame to the horizontal plates as qk ∼ δ 2 i kg ( Ta -Tu ) / h . The subindex i identifies a circular flame when i = 1 and a double-cell flame when i = 2. The factor δ 2 i represents the effective hot area in contact with the plates, kg is the gas conductivity, Ta is the adiabatic flame temperature and the plates are considered to be a heat sink at room temperature Tu .

On the other hand, the heat released by the flame per unit time can be estimated as qf ∼ ρ ULi cp ( Ta -Tu ) h δ i , where cp and ρ represent the heat capacity and density of the gas, respectively, and h δ i is the approximate flame surface. The comparison between the heat losses to the plates and the heat released by the flame yields the relative production-loss parameter

$$\Delta = \frac { q _ { k } } { q _ { f } } = \frac { \delta _ { i } } { h } \frac { \delta _ { T } } { h },$$

which must remain below unity. In this expression, the ther-

mal thickness of the flame is δ T ∼ DT / ULi , where DT = kg / ( ρ cp ) is the thermal diffusivity of the gaseous mixture. Moreover, the value ULi is used to define δ T instead of the planar flame propagation speed SL , traditionally chosen to theoretically estimate the thermal flame thickness. In the ultralean mixture regime tested in our experiments, always with equivalence ratio below φ &lt; . 0 3, planar flames do not exist 3 and, therefore, SL is not a representative value of the flame propagation velocity.

It will be shown below (figure 5) that the characteristic size of the flame kernel is similar to the gap distance δ i ∼ h . Considering this, it is easy to see from equation (1) that a progressive reduction of the distance h will increases the relative relevance of the conductive heat losses through the plates with respect to heat production. Therefore, it is expected that the corrugated continuous flame fronts will split forming isolated weak reaction kernels as h is reduced. This point has been demonstrated in recent studies, both numerical 8 and experimentally 7 , in vertical narrow channels with H2 concentrations as low as 5% vol.

The characteristic size of the isolated flames δ i is determined by the most unstable wavelength of the thermodiffusive instability 13 and by the preservation of activated chemical kinetics through curvature. The latter, in turn, relies on the preferential diffusion of hydrogen-air mixtures, where a large mass-to-thermal diffusivity ratio prescribes low values of the Lewis number Le = DT / Di = 0 3. . The onset of thermodiffusive instability accelerates the wrinkling process of the flame front, forming regions with positive and negative curvature 8 . Moreover, the aforementioned heat losses induce major inhibition of the reaction at the concave part of the wrinkled front, inducing partial extinction and dividing the reactive front into different kernels that propagate steadily, leaving a significant portion of the mixture unburned. Without the effect of thermodiffusive instabilities, minor heat losses would quench the flame, as was demonstrated by Martínez-Ruiz et al. 8 in their calculations with Le = 1. An example of the two different stable isolated fronts that are formed, either circular (slower propagation at speed UL 1 ) or double-cell (faster propagation at speed UL 2 ), can be seen in the inset of figure 1. These flames benefit from conductive heat losses through the confining plates to enable their steady propagation solutions, in a similar manner to flame balls or drifting-comet flames requiring radiation losses 14 or convective effects 6 , respectively.

An attempt to study this problem using direct threedimensional simulations is presented in the work by Melguizo-Gavilanes et al. 15 . These simulations, that used detailed chemistry and transport, validated the quasi-2D approximation for values ∆ ≪ 1 in the adiabatic case. The calculations encompass both the flow field and the thin flame, where the chemical reaction occurs. Typically, the flow field is hundreds or thousands of times larger than the flame, and the grid in both regions must be sufficiently fine to numerically integrate the equations accurately. Such disparity between the characteristics lengths scales quickly escalates the problem to reach tens of millions of unknowns. Consequently, conducting a parametric study to assess the impact of heat losses and buoyancy becomes prohibitively expensive, requir- ing hundreds of thousands of computational hours 16 , even when employing simplified chemical kinetics models.

The calculation cost can be eased by performing an asymptotic analysis in the narrow-channel limit approximation. This enables an adequate separation of scales valid when the characteristic flow length scale ℓ &gt; δ i and plate size L are larger than the gap separation h . The associated residence time at scales of the order of the structure itself is comparable to the heating time of the horizontal plates t r ∼ δ i / ULi ∼ t h ∼ h w 2 / DTw . Therefore, an effective heat transfer to the plates is achieved during the flame propagation inside the channel. In this limit, a linear variation of temperature across the thickness of the walls can be considered 8 .

All things considered, the 3D conservation equations can be substituted by a quasi-two-dimensional viscosity-controlled formulation with conductive heat losses 8 , which was recently used to delineate the parametric space of multiple stable isolated solutions by Domínguez-González et al. 9 . These results were later expanded by Domínguez-González et al. 9 . In their simulations, they found stable isolated circular and double-cell flames with both positive and negative buoyancy by tweaking the initial conditions, a result that seems to broaden the experimental observations by Veiga-López et al. 7 .

In this work, we performed new transient gravityfree simulations using the quasi-two-dimensional pseudospectral (Fourier-Fourier) code developed by DomínguezGonzález et al. 17 , to obtain the two stable configurations shown in figure 2. A comparison of the ignition transient properties with the experimental measurements will be shown below in section V. In particular, values of the dimensionless heat release q =( Ta -Tu ) / Tu =[ 2 -5 ] are chosen such that the adiabatic temperature reaches Ta =( + ) 1 q Tu ≃ [ 800 -1600 ] K, a characteristic range of ultra-lean equivalence ratios φ ≃ [ 0 2 . -0 35 . ] under atmospheric conditions 18 . Note that the flame front needs to remain curved to reach superadiabatic temperatures that stay above the crossover value, an effect that is only observed in lowLewis number fuels like hydrogen. This crossover temperature is around 1000 K at atmospheric pressure 3 and is defined as the temperature at which the rates of the branching and recombination steps in the chemical mechanism are equal.

The numerical results presented here make use of dimensionless variables, where τ = t / t r is the dimensionless time referred to t r =( DT / ULi 2 ) and x , y are dimensionless coordinates referred to δ T . The transverse-averaged energy conservation equation

$$\rho \left ( \frac { \partial \theta } { \partial \tau } + u _ { x } \frac { \partial \theta } { \partial x } + u _ { y } \frac { \partial \theta } { \partial y } \right ) = \nabla ^ { 2 } \theta + \Omega - b \theta, \quad ( 2 )$$

can be solved for the dimensionless temperature variable distribution θ = ( T -Tu ) / ( Ta -Tu ) with constant transport coefficients ρ DT = const. The parameter b = 2 h λ w / (( h / δ T ) 2 hw λ u ) , assumed to be of order unity, characterizes the off-plane conductive heat losses through the plates and represents the relative importance of the off-plane conductive heat losses through the plates. This parameter, with clear physical meaning, is obtained using a rigorous asymptotic expansions 8,19 in the limit δ i &lt;&lt; h . It involves,

in its definition, the wall thickness hw , and the walls λ w and burned gas λ u conductivities. The off-plane conductive heat losses term b θ , that appears in the right hand side of equation (2)

The conservation of reactant species reads

$$\rho \left ( \frac { \partial Y } { \partial \tau } + u _ { x } \frac { \partial Y } { \partial x } + u _ { y } \frac { \partial Y } { \partial y } \right ) = \frac { 1 } { L e } \nabla ^ { 2 } Y - \Omega. \quad ( 3 ) \quad ^ { 1 5 } \Big \|$$

Preferential diffusion is included in the model through Le , which, in the case of ultralean hydrogen flames is very well characterized 20 by a value Le = 0 3. . Notice the essential effect of preferential diffusion is the enrichment of the mixture because of curvature effects, which produces superadiabatic temperatures that allow the flame to exist below the usual flammability limit of planar flames 3 φ l ∼ 0 25. . This effect is the primary one and we have not included the Soret effect, which is anticipated to produce relatively minor quantitative corrections 21,22 . The system requires the equation of state ρ ( 1 + q θ ) = 1 and the simplified continuity equation in terms of pressure, ∇ 2 p = -q ( △ θ + Ω -b θ ) . The latter is used to compute the velocity components as a function of the pressure field from the momentum equation in the viscous limit

$$u _ { x } = - ( \partial p / \partial x ) \ \text{ and } \ u _ { y } = - ( \partial p / \partial y ). \quad ( 4 ) \quad \text{botto}$$

The reaction rate is modeled using chemical diffusion, that represents adequately the behavior of lean flames while reducing the computational cost of fully detailed kinetics 23 ,

$$\Omega = \rho ^ { 2 } \beta ^ { 2 } \left ( 1 + q \right ) ^ { 2 } \frac { Y } { 2 L e } \exp \left ( \frac { \beta \left ( \theta - 1 \right ) } { 1 + \frac { q } { 1 + q } \left ( \theta - 1 \right ) } \right ), \quad \left ( 5 \right ) \quad \text{reach} \quad \sin, \quad \text{$$

which is a one-step irreversible Arrhenius model which is well characterized for the cases under study by a moderately large Zeldovich number β = Ea ( Ta -Tu ) / ( RT a 2 ) ≈ 10, being Ea / R ≃ 20 × 10 3 K the characteristic value of global activation temperature for lean hydrogen flames 18 .

The time-evolving simulations lead to two distinct steady propagation states under the same set of numerical parameters, encompassing chemical properties (equivalence ratio φ , activation energy Ea , heat release q ) and boundary conditions (constrained flow and conductive heat losses b ). The steady solution that finally emerges as τ ≫ 1 depends only on the initial conditions set at τ = 0. Given appropriate initial profiles, θ ( x y , ; τ = ) 0 and Y x y ( , ; τ = ) 0 , the transient process is integrated in time with a third-order Runge-Kutta algorithm of time step ∆ τ = 10 -6 , limited by reaction rate terms. The interested reader is referred to Domínguez-González et al. 17 for additional details.

The characteristic final steady-state temperature fields of the circular and double-cell canonical flames are illustrated in figures 2( a ) and 2( b ) respectively. In addition, the reacting front is depicted using white isolines at 15% and 90% of the maximum value of the reaction rate Ω max in each case.

Figures 2( c ) and 2( d ) display the distribution of reactant concentration Y alongside the flow streamlines (orange). These results indicate more intense reaction rates in doublecell flame fronts, resulting in a propagation speed UL 2 approximately double that of the circular counterpart UL 1 . In the

FIG. 2. Simulations of stable isolated fronts for the same parametric set q = 5, b = 1. Steady-state dimensionless temperature θ =( T -Tu ) / ( Ta -Tu ) and reactant concentration Y are shown for a circular flame ( a , c ) and a double-cell flame ( b , d ). Isocontours of reaction rate (white) and flow streamlines (orange) are included in top and bottom panels, respectively.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00893v1-figure-ab7fe5db-1acf-4055-bb3c-3ea9c190e529.jpg)

latter, the lower velocity permits higher diffusion of reactants, sustaining the chemical reaction in the rear flame region. Conversely, higher velocities in two-cell flames impede fuel from reaching the rear of the reacting surface through mass diffusion, leading to an elongated shape with a distinct hydrogendepleted zone behind. Most important, the emergence of each fixed-point structure (within a common region of the parametric space in the ultra-lean limit) depends on the initial conditions set for the numerical simulations. This suggests that the ignition-like transient is relevant to the development of either the weaker or the stronger solution, as it will be discussed in the following experimental measurements.

## III. EXPERIMENTAL SETUP AND PROCESSING

The simplifications in the formulation of the problem introduced in the previous section to reduce the computational cost (quasi-two-dimensional problem, viscous flow, fast conductive heat losses or reduced 1-step chemical model) raise questions about the actual existence of the two types of flame structures in real devices. To investigate their simultaneous appearance, we designed and constructed the horizontal setup shown in figure 1 that eliminates the effect of gravity in the flame propagation direction. It consists of two 25 mm-thick methacrylate parallel plates of dimensions L = 950 mm and W = 220 mm, separated by a gap h ≪ W ∼ L . The gap can be modified in the range 2 ≤ h ≤ 6 mm to promote relative heat losses through the walls and prevent autoinduced acoustic pressure waves from interacting with the reactive front 24 .

The chamber is filled with a H2-air mixture prepared before injection using two mass flow controllers VÖTGLIN GSC

FIG. 3. Automatic structure detection over the Schlieren imaging frames. Flame sizes (green and blue circles) and trajectories (dashed lines) are obtained for each kernel through an in-house tracking algorithm.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00893v1-figure-ab6a0e1c-2a8c-4c82-94ae-5547cf7ba5ea.jpg)

red-y smart series. The chemical heat release and flame thickness is adjusted by carefully tweaking the equivalence ratio φ of the mixture. The uncertainty in the initial mixture composition is determined by propagation of errors. The maximum relative deviation is below 2.3%, which corresponds to ∆ φ = ± 0 005. . In order to further reduce the uncertainty, we repeated the experiments using a second set of mass flow controllers, obtaining the same flame behavior in all cases tested. As reported by the manufacturer in their calibration reports, the highest deviation is below 0.12%, which is well below the calculated uncertainty obtained through propagation of errors. The details of these calculations are included in Appendix A.

Once the chamber volume is filled with the desired reactants, we ignite the mixture at the open end using a three-pin spark plug. The spark-plug is activated using a 12 V battery and a 50 Hz (75 Hz charge cycle) square signal. It should be noted that the electric-arc ignition is not a perfectly reproducible initiation as it involves a stochastic process that slightly changes the details of how the energy is transferred to the gas. To minimize this effect, each experiment was repeated at least three times to confirm a high level of repeatability, consistent with the precision of the controllers (0.5% of full scale).

After ignition, the reaction fronts propagate inside the chamber, which is visually accessible through both plates. The reactive front is visualized through density variations using a Schlieren setup that places both the high-speed camera and a light source perpendicular to the plates, as shown in figure 1. A 30 cm-diameter concave mirror is placed below the experiment to reflect the light coming from a 150 W light source (Fiber-Lite LED Fiber Optic Illuminator) aligned with the camera lens axis. The reflected light beam is intercepted with a razor blade at a distance twice the focal length of the mirror (2 5 m). . Depending on the particular case, the experiment is recorded in the range of 800 and 1000 fps using a high-speed camera (Memrecam HX-3). Characteristic experimental images are shown in figure 3. The long-lasting darkened trails are caused by condensation of water vapor of the product gases on the surface of the plates 25,26 , while the hotter flame front can be detected in the picture as a clearer notch. The overall light-gray area illustrates the unburned mixture region.

The propagating isolated fronts are automatically detected and tracked with an in-house python code using the Open Source Computer Vision Library (OpenCV), that provides contour detection at each frame based on the lighter pixels of the flame. The reactive front is identified by thresholding the images to discriminate between burnt and unburned regions. This algorithm detects each front and identifies every same structure along consecutive time steps, which in turn enables the description of divisions and formation of new elements. Frame-by-frame post-processing consists on a series of algebraic operations (background subtraction, addition, threshold selection and binarization), which eliminates the need to correct brightness, contrast or defects between experiments.

Figure 3 illustrates a sample detection and tracking of an experimental run. The detection of the contour of each individual kernel is used to compute the circular region of minimum radius δ i that completely encloses the reactive front at each frame k . The coordinates ⃗ x k i of the center of the circle are used to delineate the path followed by the reactive front. Displacement s = | ⃗ x k i -⃗ x k -1 i | and velocity U k Li = | ⃗ x k i -⃗ x k -1 i | / ∆ t are calculated using backward differentiation with images separated by 70 frames. The estimated uncertainty of the velocity is calculated as ∆ UL / UL = ∆ s / s + ∆ t / t = 0 02. .

## IV. EXPERIMENTAL RESULTS

The following tests have been instrumental to find the simultaneous propagation of two stable flame structures experimentally, confirming the unique bi-stability predicted by the unsteady numerical simulations (figure 2) and providing additional data to map the behavior of lean hydrogen flames. First, in figure 4 (Multimedia available online) we show a composition of experimental images for increasing values of hydrogen concentration with h = 4 mm. Each panel shows an overlay of two frames, separated a time span ∆ t (case dependent) to show the evolution of the premixed flames.

The left panel depicts a group of five isolated circular flames propagating steadily in nearly straight lines for %vol. H2 = 7 95 (equivalence ratio . φ = 0 205). . The time lapse between frames is ∆ t = 2 s, and renders equivalent propagation distances of every isolated flame.

The first evidence of bi-stability emerges spontaneously for a concentration of %vol. H2 = 7 97 ( . φ = 0 206), shown in the . central panel of figure 4, where we identified both steady circular and double-cell flames simultaneously. They are originated at a single ignition event and travel steadily along

FIG. 4. Two-frame composition for propagation in channels with h = 4 mm. Isolated circular flames with %vol. H2 = 7 95 ( . φ = 0 205), . separated a time ∆ t = 2 s (left). Coexistence of circular and double-cell flames at %vol. H2 = 7 97 ( . φ = 0 206), . ∆ t = 1 s (center). Splitting fronts with %vol. H2 = 8 01 ( . φ = 0 207), . ∆ t = 1 5 s (right) (Multimedia available online). .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00893v1-figure-1aade767-0156-487e-a2b3-ca0b23bf58cc.jpg)

the chamber burning a small fraction of the total fuel available. The two-frame overlay (time lapse ∆ t = 1 s) shows that double-cell flames move, roughly, two times faster than their circular counterpart. Notice that this double-cell flame appears as a single structure, with a narrow extinction gap marked at their center by the separation of the two parallel condensation trails. This configuration is in clear contrast with the unstable adjacent 2D flame structures found numerically by Grcar 6 , in which two isolated flame cells are separated by a gap of unburned gas of roughly the same size as themselves.

Finally, double-cell flames become unsteady for concentrations above %vol. H2 = 8 00. . The rightmost panel of figure 4 shows the flame front splitting into multiple circular and double-cell flames that form an angle with the primary flame displacement. The average flame propagation speed remains constant, but the size of the reactive front oscillates before it eventually divides into more circular and double-cell flames. The frequency of double-cell flame bursting increases with increasing H2 concentration, leading to the continuous breakup of the reactive front when %vol. H2 &gt; . 9 5 ( φ = 0 25). .

as large as their circular counterparts. Furthermore, the unsteady splitting regime is found for greater hydrogen content than %vol. H2 ≃ 8 00 ( . φ = 0 207). Above this critical value, . the size of the double-cell front undergoes another transition from a steady state to a cyclic oscillation branch in which the leading flame grows and shrinks periodically before splitting to form circular flames seeded perpendicularly to the reactive path, see right panel in figure 4. The amplitude of these oscillations is small, and the limited resolution of the images does not allow for their accurate quantification except right before splitting takes place. This is precisely what we represent using the vertical lines to denote the maximum amplitude of these oscillations in the lower panel of figure 5. Both velocity and size of each type of flame remain constant during the whole concentration range tested in the experiments, including the region associated to double-cell cyclic oscillations.

The collection of experimental results for a channel height of h = 4 mm is summarized in figure 5 as bifurcation diagrams for flame speed ULi and size of the reactive front δ i versus fuel concentration φ . The critical values of φ can be drawn to determine the flame properties in the various regimes identified above. First, at ultra-lean propagation conditions %vol. H2 &lt; 7 94, only circular flames exist, with characteristic diameters . of δ 1 = 1 2 cm and steady propagation velocities that range . between UL 1 = 4 and 6 cm/s. This behavior is that of the left panel in figure 4. Above %vol. H2 = 7 95 ( . φ = 0 205), the . system bifurcates into the multiplicity regime causing the simultaneous appearance of two stable configurations, as shown in figure 4, center. Note that the size and speed of the circular fronts do not involve substantial changes compared to the previous regime. Nevertheless, double-cell fronts exhibit a size δ 2 = 2 5 cm and propagation speed . UL 2 = 11 cm/s twice

Finally, equivalence ratios even greater than %vol. H2 ≃ 8 12, produce enough heat release to avoid partial extinction . of the front and formation of isolated structures. In that case, a continuous corrugated and unstable flame front is formed, which burns the whole volume of available reactive mixture.

Although these changes in behavior and critical limits are obtained consistently for a specific separation of the plates, the relative importance of heat losses was also modified by changing the gap size in the range 2 ≤ h ≤ 6 mm. The behavior of the flame remains qualitatively similar to what has been described above for h = 4 mm, with both circular and/or doublecell flames emerging when the equivalence ratio becomes sufficiently small φ &lt; [ 0 181 0 198 0 210 0 240 0 277 . , . , . , . , . ] for h =[ 6 5 4 3 2 , , , , ] mmrespectively.

In contrast, stable circular and double-cell flames coexist exclusively for h = 3 and h = 4 mm gap sizes. In smaller ( h = 2 mm) or larger gaps ( h = 5 and 6 mm), the coexistence of solutions is restricted to stable circular flames alongside unstable double-cell flames that split periodically. This variation in flame behavior indicates that the heat loss-to-production

FIG. 5. Parametric range of multiple stable solutions for h = 4 mm. Steady propagation velocity ULi (top) and size δ i of the reaction front (bottom) of circular i = 1 (circles) and double-cell i = 2 (triangles) flames. Circular flames are always stable while double-cell flames become unstable when φ &gt; . 0 2065.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00893v1-figure-ba9a5adb-b888-453d-a7bf-6e6cfa96dc6b.jpg)

balance, ∆ = qk / qf , becomes either too small or two large to enable the steady state propagation behavior of both kinds of kernels. In fact, the splitting regime emerges when this balance is broken. The thermal energy source increases with higher hydrogen content in the mixture (greater φ ), elevating the flame and exhaust gas temperatures. When heat losses can no longer balance heat production, the flame grows and its tail starts oscillating, initiating the front-splitting phenomenon shown in figure 4 (right).

Additionally, the propagation of a stable ( φ = 0 206) and . an unstable oscillating ( φ = 0 208) double-cell flame is com-. pared in figure 7 (Multimedia available online) In the unstable case, the flame oscillates at approximately 26.6 Hz before splitting, generating multiple circular and double-cell flames. Following this event, the double-cell flame continues its propagation before resuming flickering at the same frequency. In spite of the oscillations, the propagation speed remains unaltered. Furthermore, different events arise when interacting with the bounding lateral walls or between kernels. An illustrative example is given in figure 6, including a sequence of images for %vol. H2 = 7 98 in which several . circular flames coexist with a double-cell flame that propagates steadily. These pictures illustrate the mutual interaction between the reactive fronts, changing the propagation direc- tion but leaving the modulus of the propagation velocity unchanged. Collisions with the sides of the chamber are also observed in the footage, with the only consequence being a change of propagation direction. Contrarily, the change in the morphology of double-cell flames is seen in the middle panel of figure 6, when one of them approaches the hydrogendeprived trail left by other reactive front. Under this scenario, the double-cell structure becomes a circular flame that continues its steady burning at a lower rate. These observations suggest that the smaller circular flames may exhibit greater stability to perturbations compared to the double flame structure.

## V. IRREGULAR IGNITION TRANSIENT

It has been shown that for particular combinations of heat release and relative heat losses (equivalence ratio and channel height), multiple stable solutions may arise. Slight variations in the concentration field, minor fluctuations in ignition spark energy, or subtle heat transfer imbalances between the flame and the plates could explain why a different number of circular or double-cell flames emerged in each experimental run. As a matter of fact, although both canonical solutions are fea-

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00893v1-figure-ba8399b0-2210-4e68-acdf-bd64f249a6fc.jpg)

FIG. 6. Time sequence for the simultaneous propagation of circular and double-cell flames at %vol. H2 = 7 98 ( . φ = 0 206). . Interactions between kernels and with the lateral walls of the chamber can be observed (Multimedia available online).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00893v1-figure-30f79825-b830-45f7-9731-9e8605c6cc1b.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00893v1-figure-a362d9e6-3b56-4be1-bb61-8720412c5182.jpg)

FIG. 7. Comparison between the propagation of a stable ( φ = 0 206) . and an unstable double-cell flame ( φ = 0 208). . (Multimedia available online)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00893v1-figure-e682b19d-3a87-455e-8594-e255e246678f.jpg)

sible under these common parametric-space conditions, some symmetry-breaking effects during ignition must determine the actual distribution and propagation of various types of kernels.

To test this hypothesis, we carried out an additional parametric study with several numerical simulations modifying their initial conditions. The initiation profiles involve a hot spot with radial decay, together with consistent species depletion to avoid overshooting of the reaction rate, which are perturbed with N lobes and random noise

$$\theta ( \tau = 0 ) = \theta _ { 0 } \exp \left ( - \frac { r ^ { 2 } } { 1 0 } \left [ 1 + \sum _ { i = 1 } ^ { m } \frac { \cos ( N _ { i } \theta + 2 \pi \epsilon _ { i } ) } { 2 } \right ] \right ), \ ( 6 )$$

$$Y ( \tau = 0 ) = 1 - \exp \left ( - \frac { r ^ { 2 } } { 1 0 } \left [ 1 + \sum _ { i = 1 } ^ { m } \frac { \cos ( N _ { i } \theta + 2 \pi \varepsilon _ { i } ) } { 2 } \right ] \right ), \\ \text{where } r ^ { 2 } - v ^ { 2 } \rightharpoonup v ^ { 2 } \text{ is the dimension locs radial coordinate } m$$

where r 2 = x 2 + y 2 is the dimensionless radial coordinate, m is the number of different modes included in the initiation profile, Ni is the wave-number and ε i ∈ [ 0 1 , ] is a computergenerated random shift of each mode. In axi-symmetric ignition cases ε = 0 to generate a single mode N .

Figure 8 (Multimedia available online) shows numerical solutions of temperature (top) and species (bottom) under a combination of dimensionless heat release q = 5 and conductive heat-loss coefficient b = 1, selected to ensure bi-stable solutions 9 . Moreover, the same peak temperature θ 0 = 1 8 . is selected in all cases with different initial shape distribution, determined by N = 8 modes (left panels), N = 12 modes (center panels) and a superposition of seven different modes Ni =[ 2 5 10 15 20 25 30 , , , , , , ] with random shift. Different regular initiations at τ ≪ 1 (shaped in black isocontours in figure 8) induce the formation of either double ( N = 8) or circular ( N = 12) flames, owing to the transient growth and splitting dynamics of the propagating front. However, non-regular initiation cases evolve into corrugated shapes that do not show azimuthal periodicity, nor axisymmetric solutions with only one kind of kernels. Therein, two different families of flames arise, of different sizes and which propagate at different velocities. Therefore, we expect that minor irregularities in the ignition process suffices to alter the transient history of the front and give rise to both stable-propagating structures, provided that the balance q and b (or qk / qf ) is adequate for multiplicity.

To complete this analysis, we include below in figure 9 an experimental visualization of the area near the spark plug right

FIG. 8. Numerical simulations for different initiation conditions with N = 8 (left), N = 12 (center) and irregular composition of modes Ni = [ 2 5 10 15 20 25 30 , , , , , , ] (right) distributions. Dimensionless temperature profiles at selected instants are shown (top) together with the reactive species field (bottom) (Multimedia available online).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00893v1-figure-2cc16100-ff4b-4156-9e46-2d5eeeae96fc.jpg)

after the ignition. The circular regions at the bottom of the image are the pressure sensor (blue circular region) and the plastic piece that holds the spark plug (gray region). They are both opaque, preventing transverse flame visualization with the Schlieren setup. Shortly after the spark ignition, the front is divided into multiple flames that evolve differently while propagating radially. In this figure we identify two types of reactive fronts propagating in different direction shortly after ignition. After observing the numerical results shown in figure 8, we hypothesize that the irregularities during the ignition process produced the necessary variations in the flow field to facilitate the emergence of both single and double flames over the same mixture composition φ and chamber height h . The control of the ignition details to produce perfectly symmetrical initiation cannot be addressed with the actual experimental set-up and should be studied with additional techniques and ignition devices.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00893v1-figure-0eac3102-08e3-42af-a524-0320ae7b3e42.jpg)

## VI. CONCLUSIONS

In non-linear systems, the presence of multiple solutions typically arises from the equilibrium of competing physical mechanisms. In combustion theory, solution uniqueness fades away in one-dimensional planar flames due to factors like competing non-linear chemical reactions, molecular diffusion,

FIG. 9. Flame structures generated right after the ignition event. Single and double flames can be observed outside the opaque piece mounting the spark plug (gray region) and the pressure sensor (blue region). The time difference between the two images is ∆ t = 0 5s. .

heat losses, or phase-change processes 27,28 . Fixed-point multiplicity can also be observed in higher-dimensional simulations or experimental combustion setups, often linked to hysteresis loops in dynamical systems, such as thermoacoustic

coupling 29 , flame edges 30 , droplet combustion 31 , and cool flames 32 . In these previous examples, solution multiplicity is unveiled by slowly modifying the controlling parameter. Typically, the experiment or simulation starts below a critical value when the system is at a stable steady state. As the value of this parameter is increased, the systems bifurcates and either starts oscillating or evolves to achieve a different steady state. If the experiment is repeated reducing the value of the parameter from a larger initial value, the systems evolves following a different path in which the unsteady solution is not observed or the steady solution is different. Farther reduction of this parameter to reach the critical value would give a sudden transition to return to the initial fixed point. Throughout all these examples, the evolution of the system is determined by its initial state and the bi-stability of the system is only observed in the hysteresis region. Outside this region, both solutions were never observed simultaneously.

A key distinction between these examples and our findings lies, precisely, on the absence of hysteresis in our results and the simultaneous observation of two stable solutions. The same behaviour is observed regardless of how the equivalence ratio, our controlling parameter, is varied in the experiments. Each experiment is completely independent from the others, and the evolution of the system is not affected by its initial state. In other words, the simultaneous observation of two stable solutions does not depend on how the equivalence ratio is modified.

The critical equivalence ratio φ ≃ 0 205 identifies the bifur-. cation point above which stable circular and double-cell reactive fronts coexist within a single experimental. The doublecell solution losses stability for φ &gt; . 0 207, with the flame oscillating at constant frequency before splitting, while the circular reactive front remains stable. This region extends up to φ ≃ 0 25, . when the homogeneous burning front solution is recovered.

Finally, the numerical predictions and experimental verification of this study indicate that the emergence of the two stable configurations is determined by symmetry-breaking details during ignition transients before reaching the final steady state. The importance of stochasticity during the ignition process to explain the multiplicity of solution found in the experiments has been tested by changing the initial condition in our numerical simulations. The resemblance of the initial stages of the propagation to the experimental results is very good, confirming the relevance of the ignition details to explain the onset of bistability and the subsequent evolution of the reactive front

The differences in morphology and burning speed are the kind of fundamental aspects that need to be understood to ensure a safe deployment of technology for low-emission hydrogen, whose production is expected to reach 20 Mt/year by 2030 33 . We believe our findings will pave the way for further research to enhance our understanding of the physics of H2 flames, ultimately promoting safer and more efficient utilization of this essential fuel.

## ACKNOWLEDGMENTS

This work was funded by the Regional Government of Madrid and by MCIN of Spain with funding from European Union NextGenerationEU (PRTR-C17.I1) and by Agencia Estatal de Investigación under projects TED2021-129446BC41 &amp; TED2021-129446B-C43 (MICINN/FEDER, UE).

The authors wish to thank the technical knowledge and assistance of David Díaz, Israel Pina and Manuel Santos in the design and construction of the experimental setup.

Declaration of Interests. The authors report no conflict of interest.

## DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## Appendix A: Error propagation

To check the accuracy of the measurements, we repeated the experiments using two sets of VÖGTLIN GSC red-y smart series controllers (high-performance and standard controllers), with the highest deviation being below 0.12%, as reported by the manufacturer in the calibration reports. Within the margin of error, the threshold values of the equivalence ratio separating different flame behaviors were identical with both sets of controllers. The uncertainty in the measurement of the equivalence ratio was determined using the propagation of errors theory and it represents the maximum deviation expected in the performed experiments. This value was calculated using the following expression

$$\left | \frac { \Delta \phi } { \phi } \right | = \left | \frac { \Delta \dot { m } _ { F } } { \dot { m } _ { F } } \right | + \left | \frac { \Delta \dot { m } _ { a } } { \dot { m } _ { a } } \right |,$$

with ∆ ˙ mF and ∆ ˙ ma the accuracy of the hydrogen and air mass flow controllers, respectively. The values of ∆ ˙ mF and ∆ ˙ ma are provided by the manufacturer, being ± 0 3% of full scale plus . ± 0 5% of reading in the case of high-performance controllers . and ± 1% of full scale if a standard controller is used.

As seen in the derivation above, the uncertainty depends on the specific mass flow rates used in each experiment. Hence, the values of ˙ mF and ˙ ma in our experiments were carefully selected to maintain the measurements as accurate as possible. In the most unfavorable case φ = 0 2, . we chose ˙ ma = 80 0 . g/min, ˙ mF = 0 477 g/min, . ∆ ˙ ma / ˙ ma = ± 0 0121, . ∆ ˙ mF / ˙ mF = ± 0 0107, leading to . ∆ φ / φ &lt; . 2 3%, or ∆ φ = ± 0 005. .

1 I. Staffell, D. Scamman, A. Abad, P. Balcombe, P. Dodds, P. Ekins, N. Shah, and K. Ward, 'The role of hydrogen and fuel cells in the global energy system,' Energy Environ. Sci. 12 , 463-491 (2019).

2 B. Salvi and K. Subramanian, 'Sustainable development of road transportation sector using hydrogen energy system,' Renew. Sust. Energ. Rev. 51 , 1132 - 1155 (2015).

3 A. Sánchez and F. Williams, 'Recent advances in understanding of flammability characteristics of hydrogen,' Prog. Energy Combust. Sci. 41 , 1-55 (2014).

- 4 P. Ronney, K. Whaling, A. Abbud-Madrid, J. Gatto, and V. Pisowicz, 'Stationary premixed flames in spherical and cylindrical geometries,' AIAA J. 32 , 569-577 (1994).
- 5 B. Matkowsky and G. Sivashinsky, 'Asymptotic derivation of two models in flame theory associated with the constant density approximation.' SIAM J. Appl. Math. 37 , 686-699 (1979).
- 6 J. Grcar, 'A new type of steady and stable, laminar, premixed flame in ultra-lean, hydrogen-air combustion,' Proc. Combust. Inst. 32 , 1011-1018 (2009).
- 7 F. Veiga-López, M. Kuznetsov, D. Martínez-Ruiz, E. Fernández-Tarrazo, J. Grune, and M. Sánchez-Sanz, 'Unexpected propagation of ultra-lean hydrogen flames in narrow gaps,' Phys. Rev. Lett. 124 , 174501 (2020).
- 8 D. Martínez-Ruiz, F. Veiga-López, D. Fernández-Galisteo, V. Kurdyumov, and M. Sánchez-Sanz, 'The role of conductive heat losses on the formation of isolated flame cells in Hele-Shaw chambers,' Combust. Flame 209 , 187199 (2019).
- 9 A. Domínguez-González, D. Martínez-Ruiz, and M. Sánchez-Sanz, 'Stable circular and double-cell lean hydrogen-air premixed flames in quasi two-dimensional channels,' Proc. Combust. Inst. 39 , 1731-1741 (2023).
- 10 Y. Ballossier, P. Boivin, and C. Almarcha, 'Three dimensional shapes of hydrogen-air flames within millimetric hele shaw cells,' International Journal of Hydrogen Energy 60 , 333-341 (2024).
- 11 D. Martínez-Ruiz, F. Veiga-López, and M. Sánchez-Sanz, 'Premixed-flame oscillations in narrow channels,' Phys. Rev. Fluids 4 , 100503 (2019).
- 12 P. Ronney, M. Wu, H. Pearlman, and K. Weiland, 'Experimental study of flame balls in space: Preliminary results from sts-83,' AIAA J. 36 , 13611368 (1998).
- 13 P. Clavin and G. Searby, Combustion waves and fronts in flows: flames, shocks, detonations, ablation fronts and explosion of stars (Cambridge University Press, 2016).
- 14 Y. Zeldovich, Theory of Combustion and Detonation in Gases (USSR Academy of Science, 1944).
- 15 J. Melguizo-Gavilanes, D. Fernández-Galisteo, and A. Dejoan, 'Threedimensional simulations of lean h2-air flames propagating in a narrow gap: n the validity of the quasi-two-dimensional appoximation,' International Conference on Hydrogen Safety 9 (2021).
- 16 T. Lu and C. K. Law, 'Toward accommodating realistic fuel chemistry in large-scale computations,' Prog. Energy Combust. Sci. 35 , 192-215 (2009).
- 17 A. Domínguez-González, M. P. Encinar, and D. Martínez-Ruiz, 'Numerical analysis of the bistable propagation of lean premixed hydrogen flames,' in Journal of Physics: Conference Series , Vol. 2753 (IOP Publishing, 2024) p. 012014.
- 18 D. Fernández-Galisteo, A. Sánchez, A. Liñán, and F. Williams, 'One-step
- reduced kinetics for lean hydrogen-air deflagration,' Combust. Flame 156 , 985-996 (2009).
- 19 D. Fernández-Galisteo, V. N. Kurdyumov, and P. D. Ronney, 'Analysis of premixed flame propagation between two closely-spaced parallel plates,' Combustion and Flame 190 , 133-145 (2018).
- 20 M. D. Smooke, Reduced kinetic mechanisms and asymptotic approximations for methane-air flames: a topical volume (Springer, 1991).
- 21 E. Fernández-Tarrazo, A. Sánchez, A. Liñán, and F. Williams, 'The structure of lean hydrogen-air flame balls,' Proc. Combust. Inst. 33 , 1203-1210 (2011).
- 22 E. Fernández-Tarrazo, R. Gómez-Miguel, and M. Sánchez-Sanz, 'Minimum ignition energy of hydrogen-ammonia blends in air,' Fuel 337 , 127128 (2023).
- 23 J. Chung, X. Zhang, C. Kaplan, and E. Oran, 'The structure of the blue whirl revealed,' Sci. Adv. 6 , eaba0827 (2020).
- 24 F. Veiga-López, D. Martínez-Ruiz, E. Fernández-Tarrazo, and M. SánchezSanz, 'Experimental analysis of oscillatory premixed flames in a HeleShaw cell propagating towards a closed end,' Combust. Flame 201 , 1-11 (2019).
- 25 B. Bregeon, A. Gordon, and F. Williams, 'Near-limit downward propagation of hydrogen and methane flames in oxygen nitrogen mixtures,' Combust. Flame 33 , 33-45 (1978).
- 26 M. Kuznetsov and J. Grune, 'Experiments on combustion regimes for hydrogen/air mixtures in a thin layer geometry,' Int. J. Hydrog. Energy 44 , 8727-8742 (2019).
- 27 P. Clavin, P. Fife, and B. Nicolaenko, 'Multiplicity and related phenomena in competing reaction flames,' SIAM J. Appl. Math. 47 , 296-331 (1987).
- 28 A. Bonnet, 'Non-uniqueness for flame propagation when the lewis number is less than 1,' Eur. J. Appl. Math. 6 , 287-306 (1995).
- 29 S. Etikyala and R. Sujith, 'Change of criticality in a prototypical thermoacoustic system,' Chaos 27 , 023106 (2017).
- 30 R. Daou, J. Daou, and J. Dold, 'Effect of heat-loss on flame-edges in a premixed counterflow,' Combust. Theory Model. 7 , 221 (2003).
- 31 J. Carpio, D. Martínez-Ruiz, A. Liñán, A. Sánchez, and F. Williams, 'Hysteresis in the vaporization-controlled inertial regime of nonpremixed counterflow spray combustion,' Combust. Sci. Technol. 192 , 433-456 (2019).
- 32 C. Reuter, S. Won, and Y. Ju, 'Experimental study of the dynamics and structure of self-sustaining premixed cool flames using a counterflow burner,' Combust. Flame 166 , 125-132 (2016).
- 33 IEA, 'Global hydrogen review 2023,' Tech. Rep. (International Energy Agency, 2023).
Ruben Palomeque-Santiago Alba Dominguez-Gonzalez Daniel Martinez-Ruiz Mariano Rubio-Rubio Eduardo Fernandez-Tarrazo Mario Sanchez-Sanz
Unveiling the bi-stable character of stealthy hydrogen-air flames
Ultra-lean hydrogen-air flames propagating in narrow gaps, under the
influence of cold walls and high preferential diffusion, can form two distinct
isolated structures. They exhibit either circular or double-cell shapes and
propagate at different speeds, with the latter roughly doubling in size and
traveling speed to the former. Hydrogen mass diffusivity, convective effects
and conductive heat losses are the physical mechanisms that explain the
alterations in morphology and propagation speed. In previous experiments, Veiga
et al. Phys. Rev. Lett. 124, 174501 (2020) found these clearly distinguished
flame structures for different combinations of equivalence ratio, channel gap
and the effect of gravity on the dynamics of upwards and downwards propagating
flames in a vertical chamber. Present observations in horizontal channels show
the simultaneous appearance of these two stable structures, which arise under
identical experimental conditions and conform the first evidence that
multiplicity of stable solutions coexists in real devices. To explain the
observations, we performed numerical simulations using the simplified model,
which show that symmetry-breaking details during the ignition transient explain
the concurrent emergence of the two stable configurations. This discovery urges
the need to implement additional engineering tools to account for the
possibility of formation and propagation of isolated flame kernels at different
speeds in hydrogen-fueled systems.