Dynamics of premixed hydrogen/air flames in unsteady flow Open Access

Turbulent premixed combustion is controlled by flame–turbulence interactions leading to changes in local and global flame dynamics and consumption rates. The flame surface is stretched by tangential straining and/or curvature caused by the non-uniform, unsteady flow, altering the internal flame structure as well as the local and global burning rate of fuel. At high pressures, the flame becomes thinner, and highly turbulent fluctuations may lead to flame wrinkles with extremely large curvatures and flame stretch. Accordingly, the classic flamelet approach, which employs stationary laminar unstretched flames for representing the turbulent flame structure, is not suited in this case.

The flame dynamics in terms of the correlation between laminar and turbulent flame speed with flow stretch has been extensively studied in the last few decades.^1–10 Most of these works were focused on analyzing single snapshots (e.g., for spherically expanding flames) or time-mean flame structures (e.g., for Bunsen-type flames), which assign the influence of flame stretch to integral properties such as the turbulent length scale L_t or the root mean square (rms) values of the velocity fluctuations $u ′$⁠. Hence, the time history of the unsteady flow and its effect on the flame dynamics have been neglected. However, as shown in Ref. 11 by detailed numerical simulation of oscillatory plane-jet H₂/air flames, the flame displacement is attenuated strongly in an unsteady flow, particularly in the high-frequency range of fluctuations. The flame requires a relaxation time to adapt to the imbalance of the prevailing diffusion–reaction processes caused by the unsteady flow. Similar behavior has been reported in Refs. 8, 9, and 12, applying direct numerical simulation (DNS) of 1D counter-flow oscillating, premixed H₂/air flames, and 3D turbulent methane/air flames.¹ It has been shown in the latter case that motion of the flame surface was weakened at fluctuation frequencies above the inverse of the flame transit time $τ c = δ l 0 / S l 0$ that is defined by the ratio of the unstretched laminar flame thickness $δ l 0$ and unstretched laminar flame speed $S l 0$⁠. Therefore, the time scale of the flow represents an essential parameter for studying the effect of flame–turbulence interactions, which should be considered in addition to L_t and $u ′$⁠. For studying the effect of the flow time scale on flame–turbulence interactions, it is regarded as the oscillation period of the velocity fluctuation and is normalized by the chemical time scale to obtain the Damköhler number $D a = τ t / τ c$⁠.

Although implicitly including an integral turbulent time scale via $τ t = L t / u ′$⁠, present models are not able to predict the observed effect of flow unsteadiness on flame dynamics due to incomplete consideration of the underlying flow history. Instead, most of the available combustion models assume the fast chemistry limit with $D a ≫ 1$⁠, i.e., the instantaneous adaption of the chemical reactions to the unsteady flow. Therefore, the objective of the current work is to explore the effect of unsteady stretch on the local and global flame dynamics, which could be incorporated into general combustion modeling concepts. For that purpose, detailed numerical simulations have been conducted for a variety of premixed Bunsen-type, hydrogen/air flames, being excited with harmonically fluctuating inflow velocities at prescribed frequencies. The oscillating inflow is employed to emulate the unsteady nature of turbulent flow. Compared to the previous 1D studies conducted for 1D purely strained counter-flow flames,⁹ this work assesses quantitatively the local and global flame dynamics by means of the time-mean Markstein number and the global consumption speed. Moreover, the study is conducted for a more general 2D Bunsen-type flame setup, which covers both the tangential straining and curvature effects. In this way, correlations of the local and global flame speed with unsteadiness of the flow can be quantitatively assessed. The results can be used to derive first-order estimations for modeling the time history effect of the flow on premixed flame propagation.

The work is organized as follows: Sec. II presents basic parameters used in this work for characterizing flame dynamics in terms of flame speed and stretch. The numerical setups used for the simulations are shown in Sec. III. The simulation results are discussed in Sec. IV with respect to the local and global flame dynamics. Finally, the main conclusions are summarized in Sec. V.

As illustrated in Fig. 1, the numerical simulations have been conducted for a Bunsen-type flame. The computational domain is built from a wedge with an opening angle of 3°, assuming axial symmetry for a laminar round-jet flame. Only one computational cell is used in the circumferential direction, so that the simulation domain is essentially two dimensional. Since the simulations are performed in the framework of OpenFOAM, a wedge-shaped domain is required for OpenFOAM to treat the governing equations in cylindrical coordinates. y and x indicate the streamwise and radial coordinates. The dimensions of the boundaries have been set as multiples of $δ t h 0$⁠, as indicated in Fig. 1. The domain length is set to be three times the flame height estimated from $S l 0$ and $u bulk$⁠, assuming a perfect conical flame. An equidistant grid is used for the simulations, which resolves the flame thickness with 20 cells with a grid resolution of $Δ ≈ 20 μ$m.

The simulations have been carried out with an in-house solver developed in the framework of the open-source program OpenFOAM,^14–16 which employs detailed calculations of the chemical reaction rates. The reaction mechanism developed by Li et al.¹⁷ is used for hydrogen/air combustion. The chemical source terms are treated via an operator splitting approach, where the governing equations for batch reactors are integrated in each computational cell over the computational fluid dynamics (CFD) simulation time step. Using the species concentrations in each cell at the beginning and at the end of the integration, a chemical source term that is averaged over the CFD time step can be computed and used as a source term in the governing equations. A more detailed description of this procedure can be found in Refs. 14–16 and 18.

The mixture-averaged transport model is used to consider differential diffusion, with transport coefficients of each species computed from the Chapman–Enskog solution of the kinetic gas theory. Soret and Dufour effects are not considered in this work. The latter one can generally be neglected as the contribution to the energy transport is negligible, even for hydrogen flames. As this paper focuses on very lean and rich hydrogen flames, the contribution of the Soret effect is less pronounced than at near-stoichiometric conditions. For example, the difference in the laminar flame speed with and without the Soret effect at $ϕ = 1$ is about 11%, while the difference at $ϕ = 4$ is 3% and at $ϕ = 0.5$ is 6%. The solver has been validated in previous works.^15,18–29 The time steps are 0.5 $μ$s for the $ϕ = 0.5$ flame and 0.25 $μ$s for the $ϕ = 4$ flame, leading to a maximum convective CFL (Courant–Friedrichs–Lewy) number of 0.2 and a maximum Fourier number of 0.35.

A fully implicit scheme of second order accuracy for the time integration and a fourth order, unbounded interpolation scheme for discretization of the convective and diffusive terms have been used. The numerical solver employed in this work uses a pressure–velocity coupling via the PIMPLE (SIMPLE-PISO) algorithm. While this scheme is fully compressible in the sense of fully coupling density and pressure, it behaves more like a low Mach-number approach with regard to acoustics. Due to the implicit treatment and PIMPLE approach, the simulation time step is based on the convective CFL number, not the acoustic CFL number. Therefore, acoustic waves are generally not resolved since they are not the focus of this work. Instead, the non-reflective boundary at the outlet is used to prevent reflection in the case of nonphysical numerical oscillations. Additionally, a sponge region near the outlet avoids reflection of pressure waves back into the domain.

As the simulation is fully compressible, the oscillating mass flow at the inlet results in a fluctuation of the pressure inside the domain due to momentum conservation. The pressure fluctuation is particularly large in the case of high excitation frequencies. However, as the amplitude of the oscillating mass flow and the length of the computational domain are small in this work, the maximum pressure variations observed in the computational domain are on the order of 100 Pa. Therefore, the effect of pressure on chemical reaction rates is negligible.

Figure 2 presents contours of the calculated heat release rate $q ̇$ and temperature T from the steady-state solutions (f = 0 Hz) for $ϕ = 0.5$ at the top and $ϕ = 4$ at the bottom. The solid lines in the temperature fields on the right denote the flame fronts given by the iso-surfaces of $Y Iso$ (see Table I). For the rich flame with $ϕ = 4 , q ̇$ has a maximum at the flame tip and a minimum at the flame base. The behavior is, however, reversed for the lean flame with $ϕ = 0.5$⁠, where the flame is extinguished at the tip, and $q ̇$ is largest at the flame base due to the high diffusivity of hydrogen. The flame is negatively curved at the flame tip, which results in a defocused diffusion of hydrogen from the fresh to the burnt side. Therefore, the local equivalence ratio $ϕ Loc$ is smallest at the flame tip, which leads to the highest $q ̇$ for the rich flame and the smallest $q ̇$ for the lean flame, as shown in Fig. 3 by the contour-plots of $ϕ Loc$⁠, on the left for the $ϕ = 0.5$ and on the right for the $ϕ = 4$ flame. Due to the high diffusivity of hydrogen, $ϕ Loc$ is even smaller than that at the lower ignition limit for $ϕ = 0.5$⁠, so that, the flame extinguishes at the flame tip. On the contrary, the flame yields the highest positive curvature at the flame base, which results in the largest $ϕ Loc$⁠. The correlation of heat release or burning rate with regard to flame curvature results in a negative Markstein number for the current lean flame (S_l increases with stretch) and a positive Markstein number for the considered rich flames (S_l decreases with stretch).

Figure 4 provides an overview of instantaneous contours of $q ̇$ calculated for the unsteady lean-premixed hydrogen–air flame at different f, where the flame shape deviates from that of the stationary flame shown in Fig. 2. From left to right, four snapshots of $q ̇$ are used to represent one whole oscillation period for each Da, having a time interval of $T / 4$ with $T = 1 / f$⁠. For the low frequency case with f = 125 Hz (see the first row with $D a ≈ 10$ in Fig. 4), the flame is able to adapt promptly to the oscillating inflow and moves periodically as a whole with a large displacement from its time-mean position. As f approaches the natural frequency of the flame (see the second and third rows with Da = 2.5 and 1 in Fig. 4), the flame enters a pulsating mode and becomes wrinkled. With further increased f (see the lowest row with Da = 0.25 in Fig. 4), the flame turns into its steady-state solution without discernible motion. The same behavior with respect to flame response to unsteady flow has been found for the case with $ϕ = 4$⁠. Therefore, the overall flame response to an unsteady flow, in terms of changing its flame surface area through stretch, is attenuated with the increased fluctuation frequency of the flow (see also Fig. 7). This effect will be further elucidated in Sec. IV C.

The flame particle (FP) method proposed in Refs. 13, 30, and 31 has been applied in this work to assess the effect of unsteady flow on the local flame dynamics. The method introduces a number of virtual, massless particles along the flame surface and tracks them during flame propagation. In this way, the time evolution of each individual flame surface element along with variations of flow stretch conditions and burning rates can be recorded, allowing a detailed evaluation of the time history effect of the flow on the flame's local dynamics. Figure 5 plots $S l / S l 0$ against Ka for the $ϕ = 0.5$ (top) and $ϕ = 4$ (bottom) hydrogen/air flames, which are calculated from the FPs using Eqs. (1) and (5). For the “on-the-fly” computation of S_l with Eq. (5), the line integration along the flame normal direction is performed locally over a distance of $1.5 δ t h 0$ on both sides of the flame surface. The black points indicate the data obtained from the steady-state solutions, where $S l / S l 0$ increases with Ka for $ϕ = 0.5$ and vice versa for $ϕ = 4$⁠. Therefore, the Markstein number for the stationary flames is negative for the lean flame and positive for the rich flame considered in this work. The flame is stabilized at the burner orifice and yields a radius of approximately 3 mm (8 $δ t h 0$⁠). Therefore, cellular structures attributed to thermodiffusive instabilities cannot grow for the steady-state lean flame case with $ϕ = 0.5$⁠. For the oscillating lean flames, the corrugations on the flame surface are enhanced due to the thermodiffusive effect.

For the unsteady solutions plotted in Fig. 5, the FPs are seeded at the base of the flame at a time distance of $T / 8$ with $T = 1 / f$⁠. These FPs, therefore, follow different spatial trajectories in time and experience different stretch conditions. This leads to a scattering for the distributions of $S l / S l 0$ vs Ka for the oscillating compared with the steady-state flames, as shown in Fig. 5. The diffusion–reaction balance prevailing the flame's local dynamics is disturbed by the unsteady velocity fluctuations. In this case, the flame requires a certain relaxation time for adjusting its internal structure and dynamics according to the transiently changing stretch, so that a previous state of $S l / S l 0$ is recorded by the FPs for the current Ka.

As shown in Fig. 5, the scattering of $S l / S l 0 ( K a )$ is particularly large when Da approaches unity (⁠ $D a ≈ 1$⁠). In contrast, $S l / S l 0 ( K a )$ resembles the steady-state solutions for Da <1 and $D a ≫ 1$⁠. The results indicate that the largest deviations of the local flame dynamics compared with that of the steady-state solution are due to the unsteady flow condition with $f ≈ 1 / τ c$ (⁠ $D a ≈ 1$⁠) or when the turnaround frequency of the flow approaches the flame's natural frequency. The internal structure of the flame sustained by the intrinsic reaction–diffusion balance yields the smallest resistance against external flow disturbances in this case, so that the flame is excited and falls into resonance with the unsteady flow fluctuations. Accordingly, the internal flame structure and the flame dynamics are altered considerably at $D a ≈ 1$⁠, which leads to enhanced wrinkling of the flame (see also Fig. 4 for the case with Da = 1) and strong scattering of $S l / S l 0$ with respect to Ka.

Figure 6 plots the Markstein number $M a ¯$ calculated from linear fitting of all data pairs of $S l / S l 0$ vs Ka assigned to each excitation frequency shown in Fig. 5. As the linear Markstein correlation is generally valid only for small stretch values with Ka close to 0, the fitting procedure has been performed for the range of $− 1 < K a < 1$⁠. In this way, $M a ¯$ corresponds to an unsteady or “turbulent” Ma averaged over different states of the flame during one oscillation period. The steady-state Markstein numbers Ma⁰ used in the ordinate axis of Fig. 6 for normalizing $M a ¯$ are $M a 0 = − 0.78$ for $ϕ = 0.5$ and $M a 0 = 0.10$ for $ϕ = 4$⁠. $M a ¯ / M a 0$ is less than unity in general under unsteady flow conditions, which indicates a decreased sensitivity of the local flame dynamics to stretch. The smallest values of $M a ¯ / M a 0$ are found for Da close to unity, where the oscillation frequency matches with the intrinsic frequency of the flame. The flame is excited and wrinkled in this case, as shown in the third row of Fig. 4. Consequently, $S l / S l 0$ vs Ka reveals a scattering in Fig. 5. The correlation between $S l / S l 0$ and Ka is partly positive in the range of $− 1 < K a < − 0.7$ for $ϕ = 4$ due to formation of the enclosed envelope structures for $S l / S l 0$ and Ka, so that a significant decrease in $M a ¯ / M a 0$ at Da = 1 can be detected in Fig. 6 for this case.

Figure 7 depicts the temporal evolution of the global consumption speed S_t calculated from integration of the fuel reaction rate over the whole domain, see Eq. (6). In this way, S_t corresponds to the bulk flow rate of fuel prescribed at the inflow, assuming a balance of incoming and consumed fuel mass. The time axis is normalized with the oscillation period $T = 1 / f$ in order to compare S_t calculated from different excitation frequencies or Da. The dashed lines with a sinusoidal function in Fig. 7 show the oscillation curve of the bulk inflow velocity $u bulk$⁠. S_t exhibits a phase delay $Δ α$ with respect to the inflow, as indicated in the lower part of Fig. 7, which strictly increases with the increasing f or decreasing Da. In addition, S_t yields a damped oscillation amplitude $Δ S t$ compared with that of the inflow. In the cases with $D a ≤ 1$⁠, the flame cannot respond to the highly unsteady flow and approaches its stationary solution, so that $Δ S t$ is decreased significantly and approaches 0. Figure 7 reveals that the ability of the flame to follow the flow fluctuation is strongly attenuated under the conditions of $D a ≤ 1$⁠. In this case, the oscillation amplitude of S_t is significantly weakened, but the time periods remain as given by $T = 1 / f$⁠.

The same behavior is illustrated in the phase diagram in Fig. 8, where S_t (solid curves in Fig. 7) is plotted against $u bulk$ (dashed curves in Fig. 7). The delayed response of global flame dynamics in terms of S_t to the unsteady inflow results in a hysteresis loop with a tilted elliptical shape. At the same time, the oscillation amplitude of S_t is reduced with the decreasing Da or increasing f, which is close to 0 for $D a ≤ 1$ (quasi-horizontal curves in Fig. 8). The results reveal that the overall flame dynamics are sensitive to flow stretch only at low excitation frequencies, whereas its response to the flow is strongly attenuated in the case of high frequencies.

In summary, the flame wrinkling or enlargement of flame surface caused by turbulent flow is mainly attributed to unsteady fluctuations with Da >1. The flame regime “ideally stirred reactor” proposed in the classic Borghi–Peters regime diagram for Da <1 and $R e t > 1$⁠, which is characterized by a considerably thickened flame brush and homogeneously mixed chemical state, is not visible.

Figure 9 depicts $S t / S l 0$ as a function of the total flame surface area $A t / A 0$⁠, where A₀ is the area of the inlet boundary and the black dots represent the steady-state solutions. The solid line indicates the flamelet relation with $S t / S l 0 = A t / A 0$⁠, assuming that the local consumption speed S_l is equal to the unstretched laminar flame speed $S l 0$⁠. Although not shown here, A_t oscillates in a similar way as S_t in time (as shown in Fig. 7) and no distinct phase delay between A_t and S_t has been observed. Therefore, $S t / S l 0$ and $A t / A 0$ yield a quasi-linear correlation in Fig. 9 but deviate from equality. The result reveals that the increase in S_t is due to the increase in A_t and almost not affected by the unsteady nature of the flow.

The current round-jet flame experiences an overall high negative stretch due to curvature in the circumferential direction, which is compensated by the positive stretch component due to tangential straining, see Eq. (1). For the $ϕ = 0.5$ flame, $S t / S l 0$ over $A t / A 0$ in Fig. 9 lies closely together with the reference curve (solid line) with $I ≈ 1$⁠, indicating that the global stretch Ka_t considering both curvature and straining of the whole flame surface is close to 0. In the case with $ϕ = 4$⁠, the bulk flow velocity at the inlet (⁠ $u bulk = 3 S L 0$⁠, see Sec. III and Table I) is much higher than that of the $ϕ = 0.5$ flame. Therefore, the rich flame yields overall a positive Ka_t due to the higher positive tangential strain compared with the negative stretch caused by flame curvature. This leads to I <1 according to Ma >0 for the case with Φ = 4. As a result, the correlation curves of $S t / S l 0$ vs $A t / A 0$ in Fig. 9 lie below the reference line indicated with $S t / S l 0 = A t / A 0$⁠. Note that as the lean flame with $ϕ = 0.5$ is locally extinguished at its tip (see Figs. 2 and 4), the flame surface area evaluated from the iso-surface of $Y Iso$ is overestimated. As a consequence, the distributions of $S t / S l 0$ vs $A t / A 0$ in Fig. 9 for the lean flame should be shifted toward smaller $A t / A 0$⁠, which results in an increased I. As a conclusion, the averaged flame speed over the flame surface area is larger than $S l 0$ in the case of Ma <0 and smaller than $S l 0$ for flames with Ma >0, which is valid for both stationary and transient flames.

According to the results shown in Fig. 10, the flame can follow the flow fluctuations and change its flame surface area only under the conditions of Da >10 with $E > 80 %$⁠; the ability of the flame to respond to velocity fluctuations is strongly weakened for Da <1 with $E < 5 %$⁠. This behavior for the effectiveness of global flame response to unsteady flow is illustrated in Fig. 11. The result suggests an extension of the regime diagram for flame–turbulence interaction with the asymptotic transition from a strongly responsive flame at Da >10 to a weakly responsive flame at Da <1 by a given turbulent Reynolds number. Its influence is particularly important in highly unsteady turbulent flows, where the energy cascade covers a wide range of time scales or Da. In this case, the flamelet correlation assuming a proportional increase in S_t with $u ′$ (⁠ $S t = u ′$⁠) is strongly violated due to Da <1. In addition, even an excessively large stretch rate will not lead to flame quenching at Da <1 because the flame is not able to sense the highly unsteady flow stretch.

The current work has reproduced the existing findings from the previous studies by Im and Chen⁹ applied to 1D oscillating counter-flow flames, where the flame speed has been shown to be less sensitive to stretch under unsteady flow conditions. In the current work, a round-jet flame setup is used, which leads to a substantially extended stretch range, considering both effects caused by tangential strain and curvature, and covering both negative and positive stretch. Furthermore, the global flame response in terms of the overall flame surface area is used here to evaluate the global flame dynamics, whereas the concept of flame surface is less meaningful for the planar flame setup. Therefore, the current work extends considerably the available knowledge based on 1D studies.

The obtained results can be incorporated into combustion modeling concepts in order to account for the time history effect on the response of global flame dynamics. For that purpose, a fitting function may be derived to model the correlation of E with Da shown in Fig. 10, which can be used to evaluate a cell-averaged burning rate based on S_t, for instance, in the framework of a RANS (Reynolds-Averaged Navier–Stokes) or LES (Large Eddy Simulation) concept.

Note that Da yields a non-uniform distribution along the flame surface in a general turbulent flame: on one hand, the broadband fluctuations of the flow cause variation of the local turbulent time scale, and on the other hand, the flame thickness and the flame speed can be changed locally by the turbulent flow, which results in a different chemical time scale. Furthermore, the flame–turbulence interactions take place over a wide range of time scales in a real turbulent flame, whereas a fixed forcing frequency is applied to interact with the flame in the current study.

Detailed numerical simulations have been performed for a series of oscillating premixed H₂/air flames with pre-defined frequencies f. The objective is to assess systematically the impact of flow unsteadiness on the local and global flame dynamics. The flame requires a relaxation time to respond to the oscillating flow so that flame speed S_l is less responsive to stretch under unsteady conditions. The statement is valid for both lean and rich premixed flames with oblique flow and strained and curved flame structure. As a result, the magnitude of the Markstein number $| M a ¯ |$ decreases in an unsteady flow. The flame is excited with an enhanced wrinkling, and $| M a ¯ |$ is smallest when f approaches the natural frequency of the flame, i.e., $D a ≈ 1$⁠.

The global consumption speed S_t, which correlates almost linearly with the total flame surface area, exhibits a phase delay of $Δ α$ and a reduced amplitude of $Δ S t$ with respect to the oscillating inflow. $Δ α$ increases with f, whereas the reversed trend is found for $Δ S t$⁠. The damped oscillation of S_t due to flow unsteadiness has been modeled with an efficiency factor E in dependence of Da, which reveals the asymptotic limits with $E → 0$ for $D a → 0$ and $E → 1$ for $D a → ∞$⁠.

In conclusion, the time history effect or unsteadiness of the flow plays an essential role for the correct description of local and global flame dynamics regarding flame–flow interaction. The generally used flamelet assumption is valid only for large Da, having E >0.8 for Da >10 in the current work. This work has focused on premixed hydrogen/air flames excited at fixed single frequencies. Further studies for additional fuel/air mixtures and considering the interactions between different flow time scales are required to justify the general validity of the obtained results.

The authors acknowledge the financial support of the Helmholtz Association of German Research Centers (HGF), within the research field MTET (Materials and Technologies for the Energy Transition), subtopic “Anthropogenic Carbon Cycle” (38.05.01). This work utilizes computing resources from the supercomputers Hawk at HLRS Stuttgart and HoreKa at SCC/KIT. The work leading to this publication was supported by the PRIME programme of the German Academic Exchange Service (DAAD) with funds from the German Federal Ministry of Education and Research (BMBF).

The authors have no conflicts to disclose.

Feichi Zhang: Investigation (lead); Methodology (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review and editing (lead). Thorsten Zirwes: Investigation (lead); Methodology (lead); Validation (lead); Writing – review and editing (supporting). Yiqing Wang: Investigation (equal); Methodology (equal); Software (equal). Zheng Chen: Conceptualization (equal); Investigation (equal); Methodology (equal). Henning Bockhorn: Conceptualization (equal); Methodology (equal); Supervision (equal); Writing – review and editing (supporting). Dimosthenis Trimis: Funding acquisition (equal); Project administration (equal); Supervision (equal). Dieter Stapf: Funding acquisition (equal); Project administration (equal); Supervision (equal).

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