Contributions of flame thickening and local extinctions to burning rate of intensely turbulent premixed flames

The internal structure of premixed flames can be highly influenced by the Karlovitz number, $Ka=\tau_{\mathrm{f}}/\tau_{{\eta}_{\mathrm{K}}}$, where $\tau_{\mathrm{f}}$ and $\tau_{{\eta}_{\mathrm{K}}}$ are the flame and Kolmogorov time scales, respectively [22, 23]. Flames with $Ka\lesssim 1$ usually feature both thin preheat and reaction zones, and the turbulent eddies only tend to wrinkle the flame surfaces [22, 24, 25, 26]. For $Ka\gtrsim 1$, the results reported in the literature for the preheat zone thickness are controversial. On one hand, as predicted by Peters [22], for $1\lesssim Ka\lesssim 100$, small eddies may penetrate into and broaden the preheat zone. This is confirmed through both experimental [27, 28, 9, 12, 29, 30, 31, 11] and Direct Numerical Simulation (DNS) studies [32, 33]. However, several studies [34, 25, 35, 36, 9] suggest that the preheat zone may not broaden for $Ka\gtrsim 1$. The reason for this discrepancy is elaborated in [9, 3]. Driscoll et al. [3] and Skiba et al. [9] argue that, eddies associated with $Ka\gtrsim 1$ may feature sufficient turbulent kinetic energy allowing them to penetrate and broaden the preheat zone. Skiba et al. [9] and Driscoll et al. [3] suggested that, in addition to $Ka\gtrsim 1$, flames with $u^{\prime}\Lambda/(S_{\mathrm{L,0}}\delta_{\mathrm{L}})\gtrsim 180$ may feature broadened preheat zones, where $\Lambda$, $S_{\mathrm{L,0}}$, and $\delta_{\mathrm{L}}$ are the integral length scale, unstretched laminar flame speed, and laminar flame thickness, respectively [3, 9].

Similar to the preheat zone, the results reported in the literature regarding the turbulent flames reaction zone thickness are also controversial. Experimental [9] and DNS studies of [37, 38, 39] show that the reaction zone remains thin upto $Ka\approx 550$ and 1000, respectively. In contrast to the studies of [9, 37, 38, 39], those of [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21] suggest that the reaction zone is in fact broadened compared to the laminar flame counterpart. For example, in a recent investigation [11], the authors show that for flames with $Ka\approx 76$, the reaction zone thickness can increase up to 3.9 times that of the corresponding laminar flame. Mohammadnejad et al. [11] speculated that relatively small size energetic eddies may penetrate into and broaden the reaction zone. It was shown [11] that the probability density function of eddy size distribution and kinetic energy of the eddies can be influenced by the type of the utilized turbulence generating mechanism, and this may be the reason for the discrepancy in the literature regarding the turbulent premixed flames reaction zone thickness.

The internal structure of turbulent premixed flames has implications for estimation of the turbulent burning velocity. This parameter can be estimated globally and locally, which are referred to as the global and local consumption speeds, respectively [2, 3]. These speeds are estimated using [40, 3, 2]

where $S_{\mathrm{T,GC}}$ and $S_{\mathrm{T,LC}}$ are referred to as the global and local consumption speeds, respectively [2]. $\dot{m}_{\mathrm{r}}$ and $\rho_{\mathrm{r}}$ are the reactants mass flow rate and density. $\overline{A}_{\mathrm{f}}$ is the flame surface area estimated based on the mean progress variable ($\overline{c}$). In Equation (1b), $\Sigma$ is the flame surface density, $L_{\mathrm{\xi}}$ is the length of the mean progress variable contour, $I_{0}$ is the Bray-Cant [41] stretch factor, and $\eta$ and $\xi$ are the curvilinear axes, which are respectively normal and tangent to a given $\overline{c}$ contour (usually $\overline{c}=0.5$). As reviewed in Driscoll et al. [2, 3], there are two ongoing questions regarding the global and local consumption speeds of methane-air turbulent premixed flames in the literature. First, the DNS studies of [42, 43, 44] suggest that increasing the turbulence intensity increases the flame surface area such that the local consumption speed nearly equals the global consumption speed. However, several other studies, see for example  [3, 45, 46, 40, 47, 48, 21, 49], suggest that as the turbulence intensity increases, the global consumption speed increases, but the flame surface area and the local consumption speed plateau. For instance, Wabel et al. [40] show that, at $u^{\prime}/S_{\mathrm{L,0}}\approx 150$, the global consumption speed is about 5 times the local consumption speed. Wabel et al. [40] speculated that, as the preheat zone broadens with increasing the turbulence intensity, the turbulent diffusivity of the gas in this zone increases, enhancing the transport of reactants and as a result increase of the global consumption speed. Recently, Mohammadnejad et al. [11] showed that, in addition to the preheat zone, the reaction zone also broadens with increasing the turbulence intensity, which may also contribute to the increase of the gas diffusivity. Both Gülder [46] and Nivarti et al. [45] performed theoretical calculations and showed that the enhancement of gas diffusivity inside the flame region can allow for elaborating the reason for the increase of the global consumption speed with increase of the turbulence intensity and as a result the difference between the local and global consumption speed values reported in for example [3, 45, 46, 40, 47, 48, 21, 49]. However, in both Gülder [46] and Nivarti et al. [45], the preheat and reaction zones thicknesses are not measured; and as a result, it is not known whether the enhancement of gas diffusivity is due to broadening of only the preheat zone or both the preheat and reaction zones. This requires further investigations.

The second question is related to the behavior of the global consumption speed itself. As hypothesized by Damköhler [50], also known as the first Damköhler’s hypothesis, increasing $u^{\prime}/S_{\mathrm{L,0}}$ is expected to increase the global consumption speed linearly. Although, at relatively small turbulence intensities, the global consumption speed nearly follows a linear relation with the turbulence intensity [3, 51, 22, 50, 40, 21, 49, 52], at relatively intense turbulence conditions, the global consumption speed plotted against the turbulence intensity bends towards the axis of $u^{\prime}/S_{\mathrm{L,0}}$. This is referred to as the bending behavior in the literature, see for example [40, 21, 49, 52, 53, 26]. Several reasons such as nonlinear variation of gas diffusivity with turbulence intensity [3, 50, 51, 45] as well as local quenching [3] are suggested to contribute to the bending behavior. To our best of knowledge, a potential relation between the flame quenching (or local extinction) and the bending behavior has not been studied yet.

The internal structure of hydrogen-enriched turbulent methane-air premixed flames has been recently investigated in Mohammadnejad et al. [8] and Zhang et al. [54]. The study of Mohammadnejad et al. [8] shows that, for relatively moderate turbulence conditions ($u^{\prime}/S_{\mathrm{L,0}}$ up to 2.6), the preheat and reaction zone thicknesses of the hydrogen-enriched flames are similar to the laminar premixed flame counterparts. Mohammadnejad et al. [8] showed that hydrogen enrichment (up to 50%) does not significantly influence the preheat and reaction zone thicknesses normalized by the corresponding laminar flame counterparts. Compared to [8], whose tested turbulence intensity was limited to 2.6, Zhang et al. [54] investigated both the reaction and preheat zone thicknesses for larger turbulence intensities ($3\lesssim u^{\prime}/S_{\mathrm{L,0}}\lesssim 12.5$) with hydrogen-enrichment percentages examined up to 60%. They [54] showed that, similar to the methane-air turbulent premixed flames, increasing the turbulence intensity can lead to broadening of the hydrogen-enriched methane-air flames preheat zone. However, Zhang et al. [54] showed that the reaction zone thickness does not change by increasing the turbulence intensity. Although studies of both [8, 54] provide insight into the structure of hydrogen-enriched methane-air premixed flames at relatively small and moderate turbulence intensities, the preheat and reaction zone thicknesses of these flames pertaining to relatively intense turbulence conditions remain to be investigated.

Few experimental investigations [55, 56] have studied the effect of hydrogen-enrichment on local and global consumption speeds of turbulent premixed flames. Both Halter et al. [56] and Guo et al. [55] showed that increasing the hydrogen-enrichment percentage increases the global consumption speed normalized by the unstretched laminar flame speed. Guo et al. [55] showed this parameter also increases with increasing the turbulence intensity. The values of the local consumption speed are not reported in [56], however, Guo et al. [55] showed that increasing both the hydrogen-enrichment percentage and the turbulence intensity increase the local consumption speed normalized by the unstretched laminar flame speed. Although these studies [55, 56] elaborated the effect of hydrogen enrichment and turbulence intensity on the local/global consumption speeds, their investigations are performed for $u^{\prime}/S_{\mathrm{L,0}}\lesssim 5$. To our best of knowledge, values of the local and global consumption speeds for hydrogen-enriched premixed flame are not investigated at larger turbulence intensities.

This study aims to address three objectives. The first objective of the present study is to measure the preheat and reaction zone thicknesses of hydrogen-enriched methane-air premixed flames pertaining to relatively intense turbulence conditions. To our best of knowledge, such information is not available in the literature. The second objective of this study is to develop a consumption speed estimation framework that does not rely on the flamelet assumption. Then, using this framework, we aim to understand the reason for the difference between local and global consumption speeds reported in the literature. The third objective of the present study is to characterize the bending behavior of the global consumption speed and to develop a mathematical formulation that allows to estimate the amount of bending.

The setup used in this study is a Bunsen burner identical to that utilized in [11, 34]. This burner consists of an expansion section, a settling chamber (with 5 mesh screens), a contraction section, and a nozzle with exit diameter of 22.2 mm. The burner schematic and the burner nozzle technical drawing are shown in Figs. 1 and 2(a), respectively. A mixture of 75% hydrogen and 25% methane (volumetric percentages) is injected through the burner rim, generating a pilot flame, see Fig. 1. Similar to [11], the volumetric flowrate of the pilot flow is only 2% of the main flow, which is small, and as a result, the pilot flame is not expected to influence the main Bunsen flames. In order to generate a relatively wide range of turbulence intensities, three different turbulence generating mechanisms are used similar to [11]. Specifically, zero, one, or two perforated plates were utilized. For the conditions with no perforated plate (which is the first turbulence generating mechanism), turbulence is primarily produced by the mesh screens and the Kelvin-Helmholtz instability that is developed in the jet shear layer, similar to [11, 57, 58, 59]. For the second turbulence generating mechanism, a perforated plate, with the technical drawing shown in Fig. 2(c), is placed 44.4 mm upstream of the nozzle exit plane. This plate contains 1.8 mm diameter circular holes arranged in a hexagonal pattern, generating a blockage ratio of $60$%. For the third turbulence generating mechanism, one perforated plate (identical to that used for the second turbulence generating mechanism) is positioned 44.4 mm, and a second plate is located 30.4 mm upstream of the burner exit plane, see Fig. 2(b). The second perforated plate has holes positioned on a hexagonal arrangement with a diameter of 1.8 mm. The blockage ratio of this plate is $77$%.

Simultaneous Planar Laser-Induced Fluorescence (PLIF) of hydroxyl radical ($\mathrm{OH}$) and formaldehyde molecule ($\mathrm{CH_{2}O}$) as well as separate Stereoscopic Particle Image Velocimetry (SPIV) were performed in this study. The diagnostics used here is similar to that used in Mohammadnejad et al. [11]. However, compared to [11], the SPIV experiments in the present study are performed only for the non-reacting flow conditions and to characterize the background turbulent flow. Both PLIF and SPIV data acquisition rate were 1 Hz. A dual-head Nd:YAG pump laser (Quanta Ray PIV400) and a frequency-doubled dye laser (Sirah Precision Scan) are used for generating the PLIF signals. Two laser pulses with wavelengths of 532 nm and $355\pm 3$ nm are produced by the pump laser using harmonic generators. As shown in Fig. 1, the 532 nm beam is used to pump the dye laser to generate a 283 nm ($282.94\pm 0.005$ nm) beam. This beam is used to excite $\mathrm{Q_{1}}$(6) line of the OH A-X system (1,0) vibrational band [11, 8, 60, 61, 62, 63, 64, 65, 66, 67]. The 355 nm beam is used to excite the $\tilde{A}^{1}A_{2}-\tilde{X}^{1}A_{1}4_{0}^{1}$ vibronic manifold of $\mathrm{CH_{2}O}$ similar to [11, 8, 60, 62, 63, 64, 68, 69, 70, 67]. Using mirrors ($\mathrm{M_{1}-M_{3}}$) and cylindrical lenses ($\mathrm{L_{1}-L_{3}}$) shown in Fig. 1, the 282.94 nm and 355 nm beams are converted to coincident laser sheets with the thickness and height of about $250~{}\mathrm{\mu m}$ and 40 mm, respectively. Two PIMAX ICCD cameras (see, $\mathrm{C_{1}}$ and $\mathrm{C_{2}}$ in Fig. 1), each equipped with a Nikkor UV lens and a bandpass filter with center wavelength of $320\pm 20$ nm (for OH) and Coastal Optics UV lens and a Schott GG 395 longpass filter (for $\mathrm{CH_{2}O}$), are used to collect the PLIF signals. The projected spatial resolution of both cameras is $89~{}\mathrm{\mu m}$ per pixel. The knife edge technique [71, 72] was used to obtain the line spread functions of both cameras. It was calculated that the optical system effective resolution (the full width at half maximum of the line-spread function) is $263~{}\mathrm{\mu m}$. Two photodiodes, shown in Fig. 1, are utilized to collect pulse-to-pulse variations of the 282.94 nm and 355 nm beams energies, which are utilized for reduction of the PLIF data discussed in the next section. Further details regarding the utilized PLIF system can be found in [11, 8].

The SPIV system contains a dual-cavity, double-pulse, Nd:YAG laser (BSL Twins CFR PIV200) as well as $\mathrm{C_{3}}$ and $\mathrm{C_{4}}$ sCMOS cameras (LaVision Imager sCMOS), which are equipped with Scheimpflug adapters and Tokina lenses, see Fig. 1. The laser generates pairs of 532 nm beams with a 4–33 $\mu$s separation time (depending on the tested mean bulk flow velocity). $\mathrm{M_{1}-M_{3}}$ mirrors as well as $\mathrm{L_{1}-L_{3}}$ cylindrical lenses (see Fig. 1) are utilized to convert the PIV laser beams to laser sheets that are coincident with those of the PLIF lasers. The field of view and projected spatial resolution of the SPIV cameras are 84 mm$\times$70 mm and 29.3 $\mu$m/pixel, respectively. Olive oil is atomized using a TSI Six-Jet atomizer and is fed into the burner. PLIF and SPIV images are registered and mapped to one physical coordinate system using a three-dimensional LaVision Type 20 target plate. LaVision DaVis 8.4 software is used to calculate the velocity vectors. The details of the PIV vector calculation is identical to that described in [11].

A Cartesian coordinate system, presented in Fig. 3, is used in the present study. The origin of the coordinate system is at the intersection of the burner centerline and the nozzle exit plane. The $x-$ and $y-$axes are at the imaging plane and lie inside the laser sheets. The $z-$axis is normal to the $x-y$ plane. The PLIF and SPIV field of views, shown by a green dashed window in Fig. 3, are coincident. In the figure, the flow characterization region (FCR) highlights the location at which the non-reacting turbulent flow characteristics are averaged, which is similar to [11].

In total, 68 experimental conditions are tested, with the corresponding details tabulated in Table 1. Mixture of methane (grade 2 with 99% chemical purity) and hydrogen (grade 5 with 99.999% chemical purity) is used as the fuel in this study. The fuel-air equivalence ratio ($\phi$) for all tested conditions is 0.7. All tested conditions labels are provided in the first column of the table. In this column, U$\#$ and H$\#$ indicate the tested mean bulk velocity ($U$) and the hydrogen-enrichment percentage ($\mathrm{H_{2}}$%). Four mean bulk flow velocities of 5, 15, 25, and 35 m/s are examined, which are provided in the second column of the table. $\mathrm{H_{2}}\%$ varied from 0 to a maximum of $70$%, with the corresponding values provided in the third column of Table 1. In the label of each tested condition, T$\#$ indicates the number of utilized perforated plates. Conditions with zero (T0), one (T1), and two (T2) perforated plates are highlighted by the blue, green, and red colors in the table, respectively. Cantera simulations [73] with the GRI-Mech 3.0 mechanism are used to estimate the unstretched laminar flame speeds ($S_{\mathrm{L,0}}$). Laminar flame thickness is calculated using $\delta_{\mathrm{L}}=(\lambda/c_{\mathrm{p}})/(\rho_{0}S_{\mathrm{L,0}})$, with $\lambda$ and $c_{\mathrm{p}}$ being the gas thermal conductivity and specific heat (both estimated at 1500 K), and $\rho_{0}$ being the fuel and air mixture density (estimated at 300 K), similar to [11, 3, 9, 22].

Root Mean Square (RMS) of the velocity fluctuations along the $y$ ($u^{\prime}$), $x$ ($v^{\prime}$), and $z$ ($w^{\prime}$) directions as well as the integral ($\Lambda$), Taylor ($l_{\mathrm{T}}=\Lambda(u^{\prime}\Lambda/\nu)^{-1/2}$), and Kolmogorov ($\eta_{\mathrm{K}}=(\nu^{3}/\epsilon)^{1/4}=\Lambda(u^{\prime}\Lambda/\nu)^{-3/4}$) length scales are averaged within the red dotted-dashed window of Fig. 3, and the results are presented in the table. In the equations used for calculation of the Taylor and Kolmogorov length scales, $\nu$ is the fuel and air mixture kinematic viscosity (estimated at 300 K) and $\epsilon=u^{\prime 3}/\Lambda$ [3, 74]. The integral length scale, $\Lambda$, is estimated using the formulations provided in [74] and the procedures discussed in [11, 8]. Similar to [11], conditions pertaining to the first turbulence generating mechanism feature the Kelvin-Helmholtz instability. This instability is generated by the jet shear layer, can influence the values of the RMS velocity and integral length scale, and does not correspond to the background turbulent flow. Thus, using the Proper Orthogonal Decomposition technique [75] and similar to [11], the Kelvin-Helmholtz modes are identified and excluded from the non-reacting flow data pertaining to zero perforated plates, following the procedures discussed in [11]. The Reynolds ($Re_{\mathrm{T}}$) and Karlovitz ($Ka$) numbers are estimated using $Re_{\mathrm{T}}=u^{\prime}\Lambda/\nu$ and $Ka=(u^{\prime}/S_{\mathrm{L,0}})^{3/2}(\Lambda/\delta_{\mathrm{L}})^{-1/2}$, with the values provided in Table 1. The effective Lewis number ($L_{\mathrm{eff}}$) of the fuel and air mixture was estimated using the formulations provided in Mohammadnejad et al. [8], with the corresponding values presented in the last column of the table.

All test conditions are overlaid on the Borghi–Peters [22, 76] diagram, as shown in Fig. 4. The results in the figure suggest that flames with $U=5$ m/s and no perforated plate pertain to the wrinkled flames regime. Flames with $U=5$ m/s and one perforated plate as well as Flame U15H70T1 pertain to the corrugated flames regime. The rest of the tested conditions are positioned in the thin reaction zones regime. In Fig. 4, the newly-suggested border of $u^{\prime}\Lambda/(S_{\mathrm{L,0}}\delta_{\mathrm{L}})=180$ [9, 3], which separates the flames with thin and broadened preheat zones, is also presented. This border suggests that the conditions pertaining to $U=15,25$, and 35 m/s with 2 perforated plates as well as Flame U35H00T0 are expected to feature a broadened preheat zone. Flames U25H00T0, U25H40T0, and U35H40T0 are positioned near the border of $u^{\prime}\Lambda/(S_{\mathrm{L,0}}\delta_{\mathrm{L}})=180$. The rest of the tested conditions are expected to feature a relatively thin preheat zone. The predictions of the Borghi-Peters diagram related to thin and broadened preheat zones are discussed in detail later in the results section.

Representative planar laser-induced fluorescence images pertaining to the test conditions of U5H00T0 ($Ka=0.3$), U5H40T0 ($Ka=0.1$), U35H00T2 ($Ka=76.0$), and U35H40T2 ($Ka=40.4$) are presented in the first, second, third, and fourth rows of Fig. 6, respectively. The results in the first and second columns are representative hydroxyl and formaldehyde PLIF data. The third column is the heat release rate estimated utilizing the procedure discussed in Section 3. The corresponding preheat, reaction, and combustion products zones were obtained using the procedure discussed in Section 3, and these zones are shown by the cyan, red, and yellow colors in the last column of Fig. 6, respectively. For both pure and hydrogen-enriched methane-air flames with $Ka=0.3$ (U5H00T0) and 0.1 (U5H40T0), both preheat and reaction zones remain relatively thin, see the cyan and red regions in Figs. 6(d and h). However, for both pure and hydrogen-enriched flames, as the Karlovitz number increases to large values, the preheat and reaction zones broaden significantly. This observation is consistent with the results presented in [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21].

Following the procedure elaborated in [11] and discussed in Section 3, the thicknesses of both preheat and reaction zones were estimated for each frame of the tested conditions. The averaged preheat ($\delta_{\mathrm{P}}$) and reaction ($\delta_{\mathrm{F}}$) zone thicknesses are presented in Figs. 7(a) and (b), respectively. The results in Figs. 7(a) and (b) normalized by the corresponding laminar flame preheat ($\delta_{\mathrm{P,L}}$) and reaction ($\delta_{\mathrm{F,L}}$) zone thicknesses for matching fuel-air equivalence ratio and hydrogen-enrichment percentages are also shown in Figs. 7(c) and (d), respectively. Values of $\delta_{\mathrm{P,L}}$ and $\delta_{\mathrm{F,L}}$ are tabulated in Table 3 and correspond to test conditions with no turbulence generator and at the smallest tested mean bulk flow velocity (5 m/s). The condition associated with the largest tested turbulence intensity (that is U35H00T2) features the maximum uncertainty for estimation of the preheat and reaction zone thicknesses. In order to estimate the maximum uncertainty associated with $\delta_{\mathrm{P}}$ and $\delta_{\mathrm{F}}$, the test condition of U35H00T2 was repeated three times, the maximum deviation from the mean value was estimated, and the corresponding error bars were overlaid on Figs. 7(a) and (b), respectively. The probability density functions of the test conditions with the largest instantaneous preheat and reaction zone thicknesses, which correspond to U35H30T2 and U35H10T2 test conditions, are overlaid by the hot color bar contour on Figs. 7(a) and (b), respectively.

Overlaid on Fig. 7(c) are the normalized preheat zone thickness pertaining to the studies of Skiba et al. [9], Wang et al. [21], and Mohammadnejad et al. [11], which are shown by the dark red dotted-dashed, blue dotted, and black dashed curves respectively. The studies of [9, 21, 11] correspond to pure methane-air flames. Unlike the present study that the fuel-air equivalence ratio is fixed and equals 0.7, results of [9] pertain to the fuel-air equivalence ratios of 0.65, 0.85, and 1.05, and those of [21] pertain to $\phi=0.4$, 0.7, and 1.0. Comparison of our results with those of Skiba et al. [9] and Wang et al. [21] at a similar fuel-air equivalence ratio (0.65 and 0.7) suggests that the values of $\delta_{\mathrm{P}}/\delta_{\mathrm{P,L}}$ nearly collapse. The results show that, for both pure and hydrogen-enriched methane-air flames, the preheat zone features broadening, and this becomes more pronounced with increasing the turbulence intensity. The broadening of the preheat zone for hydrogen-enriched methane-air turbulent premixed flames has been recently reported in the study of Zhang et al. [54]. However, the tested conditions in [54] correspond to to a maximum turbulence intensity of about 12.5. The results of the present study suggest this broadening extends to $u^{\prime}/S_{\mathrm{L,0}}\approx 42$ (for 10% hydrogen-enrichment).

Driscoll et al. [3] and Skiba et al. [9] suggest the turbulent premixed methane-air flames with $u^{\prime}\Lambda/({S_{\mathrm{L,0}}\delta_{\mathrm{L}}})\gtrsim 180$ may feature preheat zone broadening. Variation of the normalized preheat zone thickness versus $u^{\prime}\Lambda/({S_{\mathrm{L,0}}\delta_{\mathrm{L}}})$ is presented in Fig. 8 for all tested conditions. As can be seen, test conditions with $u^{\prime}\Lambda/({S_{\mathrm{L,0}}\delta_{\mathrm{L}}})\gtrsim 180$ feature $\delta_{\mathrm{P}}/\delta_{\mathrm{P,L}}\gtrsim 1.5$, however, those with $u^{\prime}\Lambda/({S_{\mathrm{L,0}}\delta_{\mathrm{L}}})\lesssim 180$ feature $\delta_{\mathrm{P}}/\delta_{\mathrm{P,L}}\lesssim 1.5$. For pure methane-air turbulent premixed flames, our results agree with those of [9, 3]. The results presented in Fig. 8 suggest that the border of $u^{\prime}\Lambda/({S_{\mathrm{L,0}}\delta_{\mathrm{L}}})\approx 180$ can also be used for predicting the preheat zone broadening of hydrogen-enriched methane-air turbulent premixed flames.

The normalized reaction zone thickness reported in the studies of Skiba et al. [9], Zhou et al. [12], and Mohammadnejad et al. [11] are overlaid on Fig. 7(d), which are shown by the dark red dotted-dashed, blue dotted, and black dashed lines respectively. Results of the present study show that increasing the turbulence intensity increases the normalized reaction zone thickness as shown in Fig. 7(d). This increasing trend is similar for all tested hydrogen-enrichment percentages and agrees well with the results of both Zhou et al. [12] and Mohammadnejad et al. [11]. Compared to the results of the present study and those of [11, 12], Skiba et al. [9] suggest that the reaction zone thickness remains nearly constant and does not vary by increasing the turbulence intensity. For pure methane-air turbulent premixed flames, Mohammadnejad et al. [11] experimentally showed that the reason for the broadening of the preheat zone is due to penetration of turbulent eddies into this zone. Also, Mohammadnejad et al. [11] showed that the probability density function of the turbulent eddy size is significantly influenced by the utilized turbulence generating mechanism. As a result, they [11] speculated the reason for the controversial observation reported for the reaction zone thickness in their study and in [12] compared to that in [9] is due to the difference between the utilized turbulence generating mechanisms. Nevertheless, the broadening of the reaction zone for hydrogen-enriched methane-air flames shown in Fig. 7(b and d) at relatively large turbulence intensities is reported in the present study for the first time to the best knowledge of the authors. In essence, the results presented in Fig. 7 show that both the preheat and reaction zones of the pure and hydrogen-enriched methane-air turbulent premixed flames can feature broadening. For example, the preheat and reaction zone thicknesses can feature values that are respectively 6.3 and 4.9 times those of the laminar flame counterparts as shown by the probability density functions in Figs. 7(a and b).

The local consumption speed normalized by the unstretched laminar flame speed is estimated from Eq. (1b) and the results are presented in Fig. 9(a) for all tested conditions. Please note that in Eq. (1b), the stretch factor ($I_{0}$) is assumed to be unity, similar to [40], as the effective Lewis number is close to unity for the tested conditions (please see the last column of Table 1). The maximum uncertainty in estimation of the local consumption speed is related to test condition of U35H00T2, and is shown by the error bar in Fig. 9(a). The procedure for estimation of the error bar is similar to that discussed in Section 5.1. The results in Fig. 9(a) are also color coded based on the integral length scale of the test conditions and presented in Fig. 9(b). Overlaid on both Figs. 9(a) and (b) are the results extracted from Wabel et al. [40] and Wang et al. [21], which are shown by the red dotted-dashed and blue dotted lines, respectively. The integral length scale in the study of Wabel et al. [40] varies between 6.1 and 41 mm and that in Wang et al. [21] is 2.9 mm. As can be seen, variation of $S_{\mathrm{T,LC}}/S_{\mathrm{L,0}}$ is significantly dependent on the integral length scale and follows two trends. For $\Lambda\lesssim 4$ mm, our results suggest that the normalized local consumption speed plateaus at about 2 following [21]; and, for $\Lambda\gtrsim 4$ mm, $S_{\mathrm{T,LC}}/S_{\mathrm{L,0}}$ plateaus at about 5 following [9]. Plateau of the normalized local consumption speed with increase of the turbulence intensity (evident in Fig. 9) has been a matter of discussions in the literature. Several investigations [40, 21, 45, 3] suggested the reason for this plateau is saturation of the flame surface area with increase of $u^{\prime}/S_{\mathrm{L,0}}$ at large values of this parameter. However, the results presented in Fig. 7 as well as those in [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21] suggest that increasing the turbulence intensity broadens both the preheat and reaction zones, and as a result, the flamelet assumption may not hold for such conditions. This means that Eq. (1b), which is developed based on the flamelet assumption, may lead to inaccurate estimation of the local consumption speed. Instead, here, we propose a new formulation that does not rely on the flamelet assumption, which is discussed in the following.

In order to develop a formulation that does not take into account the flamelet assumption, first, similar to [40, 21, 93, 94], the OH PLIF images are binarized and the mean progress variable field ($\overline{c}_{\mathrm{OH}}$) is obtained by averaging the binarized images. Please note this binarization is only performed to calculate the mean progress variable fields, is conducted for flames that are reported to be both relatively thin [40] and thick [21], and does not necessarily imply a presumed structure for the reaction zone. Then, a curvilinear coordinate system ($\xi-\eta$), with $\xi$ and $\eta$ being respectively tangent and normal to $\overline{c}_{\mathrm{OH}}$ contours are constructed. The time-average heat release rate is integrated along $\eta-$axis, which is locally normal to the mean progress variable contours of $\overline{c}_{\mathrm{OH}}=0.5$. This integral depends on $\xi$ and is given by

where $\overline{HRR}_{\mathrm{T}}(\eta,\xi)$ is the mean heat release rate of turbulent flames estimated locally. $\eta_{\mathrm{min}}$ and $\eta_{\mathrm{max}}$ are the extents of the $\eta-$axis in the domain of investigation. For a relatively large domain of investigation, the left hand side of Eq. (2) does not change by extending the bounds of integration. The turbulent flame burning rate ($B_{\mathrm{T}}$) is defined as the spatially averaged (along $\xi$) value of $B_{\mathrm{T,\xi}}(\xi)$ and is given by

In Eq. (3), $\xi_{\mathrm{min}}$ and $\xi_{\mathrm{max}}$ correspond to the boundaries of $\xi-$axis in the domain of investigation. In the denominator of Eq. (3), $L_{\xi,\mathrm{T}}$ is the length of $\overline{c}_{\mathrm{OH}}=0.5$ contour and depends on the tested condition. Similarly, the laminar flame burning rate can be estimated from

where $\overline{HRR}_{\mathrm{L}}(\eta,\xi)$ and $L_{\xi,\mathrm{L}}$ are the mean laminar flame heat release rate estimated locally and the laminar flame length estimated using the length of $\overline{c}_{\mathrm{OH}}=0.5$. Dividing Eq. (3) by Eq. (4) leads to

It is assumed that the heat release rate is proportional to $I_{\mathrm{HRR}}$ estimated using the $\mathrm{OH}$ and $\mathrm{CH_{2}O}$ PLIF signals following the procedures provided in Section 3. Using this assumption, Eq. (5) can be simplified to

where $\overline{I}_{\mathrm{HRR,T}}$ and $\overline{I}_{\mathrm{HRR,L}}$ are the mean $I_{\mathrm{HRR}}$ corresponding to turbulent and laminar flames, respectively. Equation (6) allows for estimation of the turbulent flame burning rate normalized by that for the laminar counterpart and is independent of the internal flame structure. In fact, in the following, it is shown that for flames whose internal structure follows the flamelet assumption, Eq. (6) leads to values of burning rate that are similar to those obtained utilizing the formulation proposed by Driscoll [2], i.e. Eq. (1b).

Variation of $B_{\mathrm{T}}/B_{\mathrm{L}}$, estimated from Eq. (6), versus $u^{\prime}/S_{\mathrm{L,0}}$ is presented in Fig. 10(a). The maximum uncertainty associated with calculation of $B_{\mathrm{T}}/B_{\mathrm{L}}$ pertains to test condition of U35H00T2, which is shown by the error bar in the figure. Among the formulations proposed in the literature, see for example the review paper by Driscoll [2], derivation of the global consumption speed formulation does not necessarily depend on the flamelet assumption, and as a result, $S_{\mathrm{T,GC}}/S_{\mathrm{L,0}}$ can be compared to $B_{\mathrm{T}}/B_{\mathrm{L}}$. Since the field of view of the present investigation is relatively small, the normalized global consumption speed values cannot be estimated for the present investigation; however, $S_{\mathrm{T,GC}}/S_{\mathrm{L,0}}$ associated with the studies of [40, 21] are overlaid on Fig. 10 by the red dotted-dashed and blue dotted curves. Two trends can be observed for the variation of the normalized burning rate. In order to highlight this, the data presented in Fig. 10(a) is color coded based on the tested integral length scale and presented in Fig. 10(b). As evident in the figure, for $\Lambda\lesssim 4$ mm, the values of $B_{\mathrm{T}}/B_{\mathrm{L}}$ are close to the normalized global consumption speed pertaining to the study of Wang et al. [21]; however, for $\Lambda\gtrsim 4$ mm, the results are closer to those of Wabel et al. [40]. Comparison of the results presented in Fig. 10 with those in Fig. 9 suggests that the values of $B_{\mathrm{T}}/B_{\mathrm{L}}$ and $S_{\mathrm{T,GC}}/S_{\mathrm{L,0}}$ are similar, and they are both significantly larger than those of $S_{\mathrm{T,LC}}/S_{\mathrm{L,0}}$.

In order to study the reason for the disparity between the values of the normalized burning rate and the normalized local consumption speed, the following parameter is defined.

Variation of $J$ versus $u^{\prime}/S_{\mathrm{L,0}}$ for the present study is shown in Fig. 11 by the open black circular data symbol. As can be seen, the values of $J$ are larger than unity, and this parameter increases with increasing $u^{\prime}/S_{\mathrm{L,0}}$. Variation of the normalized global consumption speed divided by the normalized local consumption speed (i.e. $S_{\mathrm{T,GC}}/S_{\mathrm{T,LC}}$) is extracted from the studies of Wabel et al. [40] and Wang et al. [21], and their results are presented in Fig. 11 by the solid dark red and blue triangular data symbols, respectively. Agreeing, with our results, those of past investigations [40, 21] also suggest that the local consumption speed underpredicts the turbulent premixed flames burning velocity.

The potential underlying reason for the underprediction of the burning velocity by the local consumption speed has been a matter of discussion over the past decades, see for example [3]. Gülder [46] and Nivarti et al. [45] suggested that the reason for the observed difference between the local and global consumption speed values can be explained by taking into account the enhanced diffusivity of the reactants due to penetration of small scale eddies into the flame region. Using the mathematical formulations proposed by Nivarti et al [45] and Gülder [46] as well as the non-dimensional numbers related to test conditions of the present study ($Ka$, $u^{\prime}/S_{\mathrm{L,0}}$, and $Re_{\lambda}=u^{\prime}l_{\mathrm{T}}/\nu$), the ratio of the global to local consumption speed was estimated and the results are presented in the Fig. 11 by the green dot-and-circle as well as pink double circle data symbols, respectively. As can be seen, the ratio of the burning rate to local consumption speed in this study ($J$, the black circular data points) is close to that observed in the literature [40, 21] for the ratio of the global to local consumption speeds ($S_{\mathrm{T,GC}}/S_{\mathrm{T,LC}}$). Also, the results estimated based on the mathematical formulation proposed in Nivarti et al. [45] for $S_{\mathrm{T,GC}}/S_{\mathrm{T,LC}}$ agrees well with $J$ obtained from the present study. However, values of the global consumption speed divided by the local consumption speed obtained based on the Gülder’s formulation [46] deviates from those of the present study and other experimental results [40, 21] at large turbulence intensities. The reason for the deviation of the results of the present study and the prediction from the formulation provided by Gülder [46] is speculated to be possibly linked to the assumptions made for derivation of the formulation in [46]. Specifically, Gülder [46] assumed that the enhanced diffusivity occurs only at the Taylor length scale, while Nivarti et al. [45] considered enhancement of diffusivity at all length scales smaller than the laminar flame thickness.

The underlying reason for the increasing trend pertaining to variations of $J$ and $S_{\mathrm{T,GC}}/S_{\mathrm{T,LC}}$ with $u^{\prime}/S_{\mathrm{L,0}}$ is hypothesized to be linked to the internal structure of the tested flames. In order to investigate this hypothesis, variations of $J$ verses the preheat and reaction zone thicknesses are shown in Figs. 12(a) and (b), respectively. As can be seen, $\delta_{\mathrm{P}}/\delta_{\mathrm{P,L}}$ and $\delta_{\mathrm{F}}/\delta_{\mathrm{F,L}}$ are positively correlated with $J$. Such correlation is even more pronounced for the reaction zone thickness and follows a linear trend. Thus, a line was fit to the data in Fig. 12(b) using the least-square technique. The equation of the line is give by

Equation (8) suggests that, at the limit of $\delta_{\mathrm{F}}=\delta_{\mathrm{F,L}}$, where the flamelet assumption holds, the normalized burning rate estimated from Eq. (6) equals the normalized consumption speed estimated from Eq. (1b). However, increasing the reaction zone thickness leads to larger values of the normalized burning rate compared to the normalized local consumption speed. This means that thickening of the turbulent premixed flames is linked to the deviation of the normalized burning rate (and the normalized global consumption speed) from the corresponding normalized local consumption speed values. In fact, the reason for this deviation being linked to the internal flame structure was speculated by Wabel et al. [40] earlier. However, compared to Wabel et al. [40] who suggested the reason for the deviation between $S_{\mathrm{T,GC}}$ and $S_{\mathrm{T,LC}}$ may be due to the broadening of the preheat zone, our results show such deviation is linked to broadening of both the preheat and reaction zones. It is experimentally shown that [11, 95], with increasing the turbulence intensity, the eddies can penetrate into the preheat zone. This is expected to increase the gas turbulent diffusivity, which increases the burning rate [50, 3, 45, 46], and as a result, the larger than unity values of $J$ presented in Fig. 12(a). It is speculated that, similar to the preheat zone, the turbulent eddies may also penetrate into the reaction zone, leading to increase of $\delta_{\mathrm{F}}$ (as evident in Fig. 7(b) as well as the results presented in [11, 21, 12]), increasing the turbulent diffusivity, and as a result increasing the flame burning rate. However, the above speculation (that is penetration of turbulent eddies into the reaction zone) remains to be investigated experimentally. It is proposed that Eqs. (7) and (8) can be combined with Eq. (1b) in order to modify the formulation of the local consumption speed suggested by Driscoll et al. [2] accounting for the reaction zone thickening. This equation is given by

The values of the normalized burning rate estimated from Eq. (6) along with the prediction of Eq. (9) are presented by circular and square data symbols in Fig. 13 for all test conditions. The results in this figure suggest that the local consumption speed corrected for the effect of flame thickening (Eq. (9)) follows the values of the burning rate estimated from Eq. (6). This means that the reason for the deviation between the global (here burning rate) and local consumption speeds is the thickening of the reaction zone. This is not shown in the past investigations to our best of knowledge and is demonstrated here.

Although the above argument suggests the values of normalized burning rate can be reconciled using the reaction zone thickness and the integral of the flame surface density, the reason why the variation of $B_{\mathrm{T}}/B_{\mathrm{L}}$ versus $u^{\prime}/S_{\mathrm{L,0}}$ bends toward the horizontal axis with increasing the turbulence intensity, see Fig. 10, is unclear in the literature, has been a matter of discussions [3], and is investigated here. This behavior is referred to as the bending behavior and is reported in past investigations, see for example [40, 21, 49, 52, 53, 26]. Here, first, the bending behavior is characterized by estimating the difference between the values of the normalized burning rate and the prediction of the Damköhler’s first hypothesis [50] (shown by the dashed line in Fig. 10). This difference is referred to as $D$ and is estimated by

For example, at $u^{\prime}/S_{\mathrm{L,0}}=19.4$ and $\mathrm{H_{2}}\%=50$%, $D$ is about 4.6, which is shown by the double-sided arrow in Fig. 10(a). $D$ was estimated for all tested conditions in Table 1, and the results are presented in Fig. 14(a). Also overlaid on the figure is the difference between the values of the normalized global consumption speed pertaining to the studies of [40, 21] from the dashed line in Fig. 10(a). The estimated values of $D$ are also color coded based on the integral length scale of the tested conditions and shown in Fig. 14(b). As can be seen, $D$ is nearly zero for test conditions corresponding to small turbulence intensities, but it increases with increasing $u^{\prime}/S_{\mathrm{L,0}}$. Also, the results in Fig. 14(b) show that, at a given turbulence intensity, $D$ nearly decreases with increasing the integral length scale.

The underlying reason for the above observations is hypothesized to be linked to pronounced occurrence of local extinctions and incomplete combustion at relatively intense turbulence conditions. The hydroxyl radical is a post flame species and its concentration is expected to be positively related to completion of the combustion process [22, 96]. Thus, it is expected that generation of OH radical should be inversely related to incomplete combustion as well as the occurrence of local extinction events. In fact, occurrence of local extinctions has been investigated using the OH PLIF signal in the past, see for example the studies of [97, 98, 99, 100, 101]. In the present study, the number of pixels that feature relatively large values of OH PLIF signal (more than 15% of the maximum) were obtained and averaged for all acquired frames of a given test condition, which is referred to as $\overline{N}_{\mathrm{OH}}$. In order to facilitate comparison between different test conditions, $\overline{N}_{\mathrm{OH}}$ was normalized by the number of pixels that are expected to feature relatively large time-averaged OH PLIF signal, which is referred to as $N_{\mathrm{\overline{OH}}}$. Variations of $(\overline{N}_{\mathrm{OH}}/N_{\mathrm{\overline{OH}}})^{-1}$ versus $D$ as well as $u^{\prime}/S_{\mathrm{L,0}}$ are presented in Fig. 15(a), and (b), respectively. The maximum uncertainty for estimation of $(\overline{N}_{\mathrm{OH}}/N_{\mathrm{\overline{OH}}})^{-1}$ corresponds to the test condition of U35H00T2, is estimated using the procedure discussed in Section 5.1, and is shown by the error bar in Fig. 15(a). The results in the figure show that $(\overline{N}_{\mathrm{OH}}/N_{\mathrm{\overline{OH}}})^{-1}$ is highly correlated with $D$. This correlation suggests that the occurrence of local extinctions and incomplete combustion are related to the observed bending behavior.

The results presented in Fig. 15(b) suggest that increasing the turbulence intensity increases $(\overline{N}_{\mathrm{OH}}/N_{\mathrm{\overline{OH}}})^{-1}$. As this parameter increases, the possibility of local extinction occurrence increases, the fuel can potentially escape from the flame region, decreasing the generated heat release rate, and as a result the pronounced deviation of the burning rate from the prediction of the Damköhler’s first hypothesis shown in Fig. 14. It is of interest to develop a model that allows for predicting the effect of governing parameters on the occurrence of local extinctions, the bending behavior, and as a result the burning rate. To address these, first, the results presented in Fig. 15(b) are color-coded based on the corresponding tested integral length scale and presented in Fig. 15(c). As evident in the inset of the figure, $(\overline{N}_{\mathrm{OH}}/N_{\mathrm{\overline{OH}}})^{-1}$ does not change by changing the integral length scale at relatively small values of $u^{\prime}/S_{\mathrm{L,0}}$. However, at relatively large values of the turbulence intensity, decreasing the integral length scale increases the values of $(\overline{N}_{\mathrm{OH}}/N_{\mathrm{\overline{OH}}})^{-1}$ (compare the values of this parameters at $u^{\prime}/S_{\mathrm{L,0}}\approx 20$ and at different integral length scales). At relatively small turbulence intensities, the turbulent eddies cannot penetrate into the flame zone, local extinctions rarely occur, and the values of $(\overline{N}_{\mathrm{OH}}/N_{\mathrm{\overline{OH}}})^{-1}$ are relatively small and nearly independent of the tested integral length scales. However, at larger values of $u^{\prime}/S_{\mathrm{L,0}}$, eddies feature relatively large rotational kinetic energy and can penetrate into the flame zone [11]. At large values of $u^{\prime}/S_{\mathrm{L,0}}$, decreasing the integral length scale facilitates penetration of relatively small scale eddies into the flame zone, increasing the occurrence of local extinctions, and as a result the larger values of $(\overline{N}_{\mathrm{OH}}/N_{\mathrm{\overline{OH}}})^{-1}$ and $D$ shown in Figs. 15 and 14(b). In essence, the above argument suggests that both $(\overline{N}_{\mathrm{OH}}/N_{\mathrm{\overline{OH}}})^{-1}$ and $D$ should be positively (negatively) correlated with $u^{\prime}/S_{\mathrm{L,0}}$ ($\Lambda$). To investigate this, variations of $(\overline{N}_{\mathrm{OH}}/N_{\mathrm{\overline{OH}}})^{-1}$ and $D$ versus $({u^{\prime}}/{S_{\mathrm{L,0}}})({\Lambda}/{\delta_{\mathrm{L}}})^{-1}$ are presented in Figs. 16(a) and (b), respectively. The parameter $({u^{\prime}}/{S_{\mathrm{L,0}}})({\Lambda}/{\delta_{\mathrm{L}}})^{-1}$ is similar to that used for calculation of the Karlovitz number, however, different exponents for ${u^{\prime}}/{S_{\mathrm{L,0}}}$ and ${\Lambda}/{\delta_{\mathrm{L}}}$ are used here. This is because, when presented against $({u^{\prime}}/{S_{\mathrm{L,0}}})({\Lambda}/{\delta_{\mathrm{L}}})^{-1}$ (compared to the Karlovitz number), both $(\overline{N}_{\mathrm{OH}}/N_{\mathrm{\overline{OH}}})^{-1}$ and $D$ variations collapse. Using the least-square technique, a second order polynomial curve was fit to the variation of $D$ versus $({u^{\prime}}/{S_{\mathrm{L,0}}})({\Lambda}/{\delta_{\mathrm{L}}})^{-1}$, and it is obtained that

Combining Eqs. (10) and (11), it can be shown that the normalized burning rate is obtained from

In Eq. (12), Term I pertains to the Damköhler’s first hypothesis, and Term II highlights the negative contribution of bending to the normalized burning rate of turbulent premixed flames. To assess prediction of Eq. (12), first, values of the normalized burning rate color-coded based on the normalized integral length scale are presented in Fig. 17. Then, the values of $B_{\mathrm{T}}/B_{\mathrm{L}}$ obtained from Eq. (12) are overlaid on the figure for several values of $\Lambda/\delta_{\mathrm{L}}$ by the dotted-dashed curves. The prediction of Damköhler [50] is also shown by the black dashed line. The results in Fig. 17 suggest that the formulation proposed in Eq. (12) allows for prediction of the normalized burning rate. The results also show that, flames with larger $\Lambda/\delta_{\mathrm{L}}$ feature less pronounced bending behavior; and, at a fixed value of $u^{\prime}/S_{\mathrm{L,0}}$, increasing $\Lambda/\delta_{\mathrm{L}}$ increases the normalized burning rate. In fact, at the limit of $\Lambda/\delta_{\mathrm{L}}$ approaching infinity, Eq. (12) predicts that $B_{\mathrm{T}}/B_{\mathrm{L}}$ approaches the prediction of the Damköhler’s first hypothesis shown by the dashed line in Fig. 17. Also, for a fixed value of $\Lambda/\delta_{\mathrm{L}}$, the variation of the normalized burning rate versus $u^{\prime}/S_{\mathrm{L,0}}$ features a parabolic behavior. For extremely large values of $u^{\prime}/S_{\mathrm{L,0}}$, the effect of Term II in Eq. (12) becomes dominant, which is equivalent to pronounced occurrences of local extinctions, leading to decrease of $B_{\mathrm{T}}/B_{\mathrm{L}}$.

The bending behavior and decrease of the burning rate due to local extinctions are consistent with the results of past studies, see for example [26, 7, 102, 103, 104]. The implication of Eq. (12) prediction and the results presented in Fig. 17 is that the burning rate cannot increase indefinitely with increasing the turbulence intensity. However, at very large normalized integral length scales, the amount of bending is less pronounced. This conclusion is of relevance to DNS studies. Compared to experimental investigations, the DNS studies tend to utilize a relatively small domain of investigation and as a result a relatively small integral length scale, see for example studies of [32, 33, 44, 43]. Our results suggest, increasing the size of the domain of investigation and the integral length scale could potentially allow for the DNS investigations to capture relatively larger values of the burning rate.