Inner Flame Front Structures and Burning Velocities of Premixed Turbulent Planar Ammonia/Air and Methane/Air Flames

Abstract

Ammonia ($\mathrm {NH_3}$) has attracted interest as a future carbon-free synthetic fuel due to its economic storage and transportation. In this study, quasi direct numerical simulations (quasi-DNS) with detailed-chemistry have been performed in 3-D to examine the flame thickness and assess the validity of Damköhler’s first hypothesis for premixed turbulent planar ammonia/air and methane/air flames under different turbulence levels. The Karlovitz number is systematically changed from 4.26 to 12.06 indicating that all the test conditions are located within the thin reaction zones combustion regime. Results indicate that the ensemble average values of the preheat zone thickness deviate slightly from the thin laminar flamelet assumption, while the reaction zone regions remain relatively intact. Following the balance equation of reaction progress variable gradient, normal strain rate and the tangential diffusion component of flame displacement speed variation in the normal direction to the flame surface are found to be responsible for thickening the flame. However, the sum of reaction and normal diffusion components of flame displacement speed variation in the normal direction to the flame surface is in charge of flame thinning for ammonia/air and methane/air flames. In addition, the validity of Damköhler’s first hypothesis is confirmed by indicating that the ratio of the turbulent burning velocity to the unstrained premixed laminar burning velocity is relatively equal to the ratio of the wrinkled to the unwrinkled flame surface area. Furthermore, the probability density functions of the density-weighted flame displacement speed show that the bulk of flame elements propagate identical to the unstrained premixed laminar flame.

1 Introduction

An improved comprehension of premixed turbulent combustion physics is required to enhance the efficiency of certain class of combustion devices such as gas-turbines for power generation and internal combustion engines for transportation (Peters 2000). Many attempts have been made in the past to decrease greenhouse gas (GHG) emissions through increasing the efficiency of such combustion devices (Kobayashi et al. 2019). However, meeting the low emission aims requires using carbon-free synthetic fuels in industrial sectors as well. Using hydrogen ($\mathrm {H_2}$) produced from green-energy origins seems to be one of the possible solutions for attaining this target. However, economic transportation and storage of pure hydrogen are still not resolved, which is an obstacle in its sustainable operation (Kobayashi et al. 2019). This emerges the need of investigating other fuel compounds as well. Among all hydrogen carrier compounds, ammonia ($\mathrm {NH_3}$) has the superiority of being utilized as a fuel in combustion devices since it has a very high hydrogen density and its transportation, storage, handling, distribution, and sustainable production are quite well established (Zamfirescu and Dincer 2008; Valera-Medina et al. 2018; Kobayashi et al. 2019). The boiling temperature and condensation pressure of ammonia are similar to propane, and it could be effortlessly utilized in marine engines running on propane (Valera-Medina et al. 2018; Kobayashi et al. 2019). It should be emphasized that using pure ammonia as a fuel has its own disadvantages due to its low unstrained premixed laminar burning velocity and heat of combustion (Kobayashi et al. 2019). While using ammonia results in eliminating carbon from combustion products, the level of nitrogen oxides ($\mathrm {NOx}$) emissions from each combustion device running on pure ammonia should be carefully addressed. Successful implementation of utilizing ammonia as a primary fuel for combustion devices has been reviewed and discussed thoroughly in Valera-Medina et al. (2018) and Kobayashi et al. (2019). Despite the active research in the field of premixed turbulent combustion, systematic numerical investigation on instantaneous inner flame front structures and burning velocities of premixed turbulent ammonia/air flames evolved under different turbulence intensities has not been carried out. The current research gap could be due to the low burning velocity of pure ammonia/air flames making it a less interesting fuel to study compared to other fuels. It has been argued that the fundamental properties of $\mathrm {NH_3}$ as a main fuel could be enhanced with adding hydrogen ($\mathrm {H_2}$) and/or enriching it with oxygen ($\mathrm {O_2}$) (Kobayashi et al. 2019). However, the level of $\mathrm {NOx}$ produced from these flames increases significantly (Li et al. 2015; Nozari and Karabeyoğlu 2015) requiring more research on its sustainable use. Considering this, since ammonia as a main fuel is utilized in some combustion devices, e.g. a swirl-stabilized burner reported in Somarathne et al. (2016), there is a need to address the fundamental properties of premixed turbulent ammonia/air flames in the literature.

One of the important characteristics of instantaneous inner flame front structures of premixed turbulent flames is the local flame thickness. It should be noted that the experimental and numerical results of local flame thickness evolution under different turbulence intensities are contradictory, and there is not a very precise answer for this fundamental issue (Tamadonfar and Gülder 2015; Yu et al. 2019). Several experimental studies show the thick preheat zone region for methane/air and propane/air flames stabilized on Bunsen-type and swirl-stabilized burners, see e.g. Mansour et al. (1998) and O’young and Bilger (1997). However, thin flame is observed for methane/air flames on Bunsen-type and V-shaped burners, see e.g. Soika et al. (1998) and Buschmann et al. (1996). de Goey et al. (2005) and Tamadonfar and Gülder (2015) showed that the ratio of the turbulent mean flame thickness to the unstrained premixed laminar flame thickness of methane/air flames is dependent on the equivalence ratio. Wabel et al. (2017) showed that the preheat zone region tends to be exceptionally widened under extremely high turbulence intensities, while the reaction zone region retains its thin characteristics. Since the two cases studied in Wabel et al. (2017), showing broadened preheat/thin reaction layers, were located in the predicted boundary for the broken reaction zones regime, they concluded that the boundary limit for this region should be re-examined in the Borghi-Peters regime diagram. More recently, Mohammadnejad et al. (2020) observed broadened preheat and reaction layers under high turbulent Karlovitz number Bunsen-type flames. They suggested that the penetration of energetic eddies into the preheat and reaction zones might be the active mechanism behind the observed trend. Furthermore, Driscoll (2008) argued that the formation of thick flames observed during experiments might be due to the misalignment of laser sheet to the flame front layers and/or the effect of strain rate on the flame front, and the penetration of small scale structures might not be the main reason for this observation. Aspden et al. (2016, 2019) observed a thick preheat zone region for turbulent planar methane/air and hydrogen/air flames under extremely high Karlovitz numbers using direct numerical simulation (DNS) and implicit large eddy simulation (ILES). They showed that the reaction zone region is still thin for both mixtures with increasing Karlovitz number. In their DNS studies, Sankaran et al. (2015) showed that the preheat layer tends to be thickened, but the reaction zone layer is not considerably altered for premixed turbulent methane/air Bunsen-type flames.

Since the majority of experimental studies available in the literature utilized planar imaging technique for evaluating the flame thickness, see e.g. Tamadonfar and Gülder (2015), due to extreme challenges in 3-D measurements, introducing a correction factor for 2-D projection data and validating it over DNS data may result in a better comprehension of limitations/challenges happening during measurements. Therefore, the experimental data could be used in a more robust manner for developing/assessing new models (Chakraborty and Hawkes 2011). In addition, Sankaran et al. (2007) explored the effects of strain rates and flame front curvature on local flame thickness by examining the conditional means of flame thickness and the balance equation of reaction progress variable gradient. They discussed that the normal strain rate and the flame displacement speed variation in the normal direction to the flame surface are responsible for controlling the scalar gradient. Their methodology helps to unlock the underlying active mechanisms behind flame thinning/broadening. A similar approach was recently implemented in Wang et al. (2017) showing the relative importance of normal strain rate and flame self-propagation on the flame thickness of a high Karlovitz number premixed jet flame. In their study, it has been shown that the variation of the flame displacement speed in the normal direction to the flame surface (normal strain rate) has a thickening (thinning) effect. A similar observation was reported by Nilsson et al. (2018) for premixed turbulent planar methane/air flames under high Karlovitz number test cases. Furthermore, Kasten et al. (2021) discussed the flame front orientation statistics and its impact on normal strain rate behavior. Since the latter has an impact on flame thinning/broadening, as mentioned in Dinkelacker et al. (1998); Chakraborty (2007); Sankaran et al. (2007); Wang et al. (2017); Nilsson et al. (2018), investigating the flame front orientation statistics can deepen our understanding on flame front topologies.

As mentioned earlier, the turbulent burning velocity is a key fundamental property of premixed turbulent flames (Peters 2000; Poinsot and Veynante 2012), and there exists numerous works in the literature determining a model for predicting the turbulent burning velocity, see e.g. Liu and Lenze (1989), Gülder (1991), Peters (1999), Fogla et al. (2017), and the comprehensive review paper by Lipatnikov and Chomiak (2002) reviewing the available empirical correlations in the literature. It should be noted that the turbulent burning velocity modeling is connected to the Damköhler’s first hypothesis relating the turbulent burning velocity to the flame surface area (Peters 1986, 1999). According to this well-known hypothesis (Damköhler 1940), the ratio of the turbulent burning velocity to the unstrained premixed laminar burning velocity is equal to the ratio of the wrinkled to the unwrinkled flame surface area. On one hand, this hypothesis was successfully confirmed using DNS studies of Chakraborty and Cant (2005), Hawkes and Chen (2006), Han and Huh (2008), Nivarti and Cant (2017), Ahmed et al. (2019), Klein et al. (2020), Brearley et al. (2020), and Varma et al. (2021) for statistically premixed turbulent planar flames with unity Lewis number. It should be emphasized that the primary fuels for these studies are mostly methane and hydrogen. More recently, Sabelnikov et al. (2019) assessed the validity of this hypothesis under high Karlovitz number test cases by analyzing the propagation of single-reaction waves in the constant density turbulent flow using DNS. On the other hand, Gülder (2007) declared that there exists a body of experimental data not supporting the Damköhler’s first hypothesis. Tamadonfar and Gülder (2015) showed that this hypothesis might not be quite credible for premixed turbulent Bunsen-type flames under certain experimental conditions using methane/air, ethane/air, and propane/air flames. More recently, Chakraborty et al. (2019) showed that the Damköhler’s first hypothesis might be invalid for flames having high non-zero mean flame front curvature such as Bunsen-type flames, in line with Tamadonfar and Gülder (2015), but it stays valid for statistically premixed turbulent planar unity Lewis number flames using DNS. The validity of this hypothesis might be disputed when the global Lewis number is deviated from a unity value. For example, Haworth and Poinsot (1992) and Trouvé and Poinsot (1994) reported that the ratio of the turbulent burning velocity to the unstrained premixed laminar burning velocity is higher (lower) than the ratio of the wrinkled to the unwrinkled flame surface area when the Lewis number is lower (higher) than unity using DNS with simple chemistry. A similar dependency between the turbulent burning and the Lewis number was previously discussed by Bell et al. (2007), Han and Huh (2008), and Chakraborty et al. (2014) based on DNS using simple/detailed chemistry, and it was reviewed comprehensively by Lipatnikov and Chomiak (2005).

To the best of the authors’ knowledge, certain aspects of instantaneous flame front inner structures as well as the validity of Damköhler’s first hypothesis and their dependence on turbulence characteristics for premixed turbulent ammonia/air flames are not explored in the literature. Therefore, the target of the current study is to address some of these issues by investigating the effect of turbulence intensity on the aforementioned characteristics of statistically premixed turbulent planar ammonia/air and methane/air flames using quasi direct numerical simulations (quasi-DNS) with detailed-chemistry. It should be emphasized that, to the best of our knowledge, performing such computationally expensive simulations were not feasible with the currently available open-source solvers until recently. In this work, the recently developed open-source dynamic load balancing model (DLBFoam) (Tekgül et al. 2021) is utilized for simulating these test conditions in the open-source CFD toolbox OpenFOAM. Since the majority of works in the literature utilize methane/air flames in their studies, it is of particular interest to use it as a benchmark fuel/air mixture to compare its results with ammonia/air flames. In this respect, the main objectives of the current work are (1) to study the underlying mechanisms behind the thinning/broadening of the flame thickness and present the budget terms associated with this phenomenon, (2) to explore the flame front orientation statistics, and (3) to assess the validity of Damköhler’s first hypothesis for ammonia/air and methane/air flames. The remainder of the paper is organized as follows. The numerical methods will be discussed thoroughly in Sect. 2. This will be followed by discussing the results in Sect. 3. Finally, the main concluding remarks will be summarized in Sect. 4.

2 Numerical Methods

2.1 Flame Configuration

In this work, the analysis is based on numerical simulations of statistically premixed turbulent planar methane/air and ammonia/air flames with an equivalence ratio of 0.7 and 1.0, respectively. Choosing the aforementioned equivalence ratio values is based on the fact that the fuel Lewis number for both mixtures is approximately near the unity value. This implies that the instability induced by thermo-diffusivity is of minor importance indicating that these flames are globally thermo-diffusively neutral. Another important reason for selecting the aforementioned stoichiometry for ammonia/air flames is that the total $\mathrm {NOx}$ emissions reduce significantly near stoichiometric and slightly rich mixtures with its minimum value around $\phi = 1.1$ depending on the reactants temperature (Kobayashi et al. 2019). The reactants gas temperature $T_{\mathrm {r}}^{}$ is set to 320 K under standard atmospheric pressure. This results in an unstrained premixed laminar burning velocity $S_{\mathrm {L}}^{0}$ of approximately 22.16 (7.99) cm/s for methane/air (ammonia/air) flames. The unstrained premixed laminar thermal flame thickness $\delta _{\mathrm {th}}^{}$ evaluated from the maximum temperature gradient is equal to 0.6375 (1.9156) mm for methane/air (ammonia/air) flames. The adiabatic flame temperature $T_{\mathrm {ad}}^{}$ for methane/air (ammonia/air) flame is approximately equal to 1854 (2079) K. For simplification, the unity Lewis number is considered for all species indicating that the mass diffusivity of all species is equal to the thermal diffusivity of the mixture, and the perfect gas law is assumed for all species. A similar simplification has been previously performed in the literature, see e.g. Boger et al. (1998); Chakraborty (2007); Han and Huh (2008); Dopazo et al. (2015); Brearley et al. (2020); Kasten et al. (2021); Varma et al. (2021).

The initial isotropic homogeneous turbulent flow field is initiated using the Passot-Pouquet spectrum (Saad et al. 2017). The initial turbulence intensity $u_{}^{\prime }$ and the initial integral length scale $\Lambda$ are specified following this spectrum on a cube with a side length equal to $26\delta _{\mathrm {th}}^{}$ ($20\delta _{\mathrm {th}}^{}$) for methane/air (ammonia/air) flames. Flame-turbulence interactions happen under decaying turbulence, similar to Boger et al. (1998); Hawkes and Chen (2004); Chakraborty (2007); Han and Huh (2008); Chakraborty and Hawkes (2011); Chakraborty et al. (2014); Nivarti and Cant (2017); Bhide and Sreedhara (2020); Brearley et al. (2020). More information on different turbulence generation methods for statistically premixed turbulent planar flame configurations have been extensively reviewed and discussed in Klein et al. (2017). The generated turbulent flow is then superimposed on the solution of 1-D unstrained premixed laminar flame corresponding to the equivalence ratio of 0.7 (1.0) for methane/air (ammonia/air) flame. The grid spacing guarantees that 13–17 computational cells reside across the flame thickness depending on the turbulence intensity of the studied case. This ensures that the spatial resolution is sufficient enough for resolving the flame front structure. In addition, the initial Kolmogorov length scale $\eta \approx \Lambda Re_{\Lambda }^{-3/4}$ is equal or larger than the grid resolution for all the test conditions, where $Re_{\mathrm {\Lambda }}^{}$ is the turbulent Reynolds number. It should be mentioned that the initial Kolmogorov length scale grows with time considering the decay of turbulence (Pope 2000). The initial non-dimensional turbulence intensity $u_{}^{\prime }/S_{\mathrm {L}}^{0}$, non-dimensional turbulent integral length scale $\Lambda /\delta _{\mathrm {th}}^{}$, turbulent Reynolds number $Re_{\mathrm {\Lambda }}^{} = u_{}^{\prime }\Lambda /\nu$, Karlovitz number $Ka = (u_{}^{\prime }/S_{\mathrm {L}}^{0})^{1.5}(\Lambda /\delta _{\mathrm {th}}^{})^{-0.5}$, and Damköhler number $Da = \Lambda S_{\mathrm {L}}^{0}/u_{}^{\prime }\delta _{\mathrm {th}}^{}$ are presented in Table 1 for the six test conditions studied in this work (cases A–H), where $\nu$ is the reactants kinematic viscosity. According to Table 1, $u_{}^{\prime }/S_{\mathrm {L}}^{0}$ is increased systematically from 4.0 to 8.0 for each fuel/air mixture, while $\Lambda /\delta _{\mathrm {th}}^{}$ is kept constant. Figure 1 shows all the numerical test conditions on a Borghi–Peters regime diagram for premixed turbulent combustion. These data indicate that all the test cases lie within the thin reaction zones regime, which is beyond the Klimov-Williams criterion, i.e. $Ka = 1$ (Peters 2000).

Table 1 Simulation parameters studied in the current work

Fig. 1

Numerical data on a Borghi–Peters regime diagram for premixed turbulent combustion (Borghi 1985; Peters 1999)

2.2 Reactive Flow Solver

Quasi direct numerical simulations (quasi-DNS) with detailed-chemistry for statistically premixed turbulent planar flames have been performed using the open-source C++ CFD toolbox OpenFOAM (Weller et al. 1998), where mass, momentum, species mass fractions, and energy conservation equations are solved with compressible PIMPLE algorithm. The terminology of quasi is borrowed from the recent work of Zirwes et al. (2020), where they reviewed and discussed the capability of OpenFOAM for conducting DNS. The detailed chemical kinetic mechanism of GRI-Mech 3.0 (Smith et al. 1999) with 53 species and 325 reactions is used for methane/air flames, and the chemical kinetic mechanism of Stagni et al. (2020) with 31 species and 203 reactions is utilized for ammonia/air flames. The performance of the aforementioned chemical kinetic mechanism against a wide range of experimental data is reported in Shrestha et al. (2021). The reactingFoam solver utilized in the current work is linked to the open-source library of pyJac (Niemeyer et al. 2017) into the OpenFOAM framework. The aforementioned library provides C subroutines for evaluating the analytical Jacobian of a prescribed chemical kinetic mechanism needed for solving a system of ordinary differential equations (ODE) for calculating the chemical source terms. These chemical source terms are the heat release rate for the energy equation, and the net production/consumption rate of each species for the species transport equations. In our study, every computational cell is considered as a stand-alone chemistry problem as where their state is defined through pressure, temperature, and species mass fractions, i.e. local thermochemical composition vector. It should be emphasized that since the chemical time scales are much smaller than the flow time scale, the chemical source terms introduce high stiffness to the system of ordinary differential equations. Therefore, an operator-splitting approach is utilized in this work to separate the calculation of the chemistry solution from the fluid flow solution. These source terms implying the rate of change of the local thermochemical composition are then evaluated by solving the system of ordinary differential equations using semi-implicit Euler method (Seulex) through integrating in each chemistry problem over the CFD time step (Hairer and Wanner 1996).

To overcome the computational load imbalance among processors, an open-source dynamic load balancing model (DLBFoam), developed recently by Tekgül et al. (2021), is utilized to distribute the computation of chemistry problems equally among all processors using the MPI communication protocol. The DLBFoam model may be the first open-source implementation of chemistry dynamic load balancing reported in the literature enabling such computationally expensive simulations. Inflow and outflow boundary conditions in the direction of mean flame propagation are utilized with zero mean inlet velocity, and periodic boundary conditions are employed for lateral directions illustrating an unbounded planar flame brush. For pressure, a non-reflecting boundary condition is implemented for both inflow and outflow. Time discretization has been performed using a second-order implicit backward Euler method as well as the cubic discretization scheme for spatial derivatives. Each simulation runs for at least two initial eddy turn over times $t_{\mathrm {e}}^{} = \Lambda / u_{}^{\prime }$ ensuring that the simulation time is equal/larger than the chemical time scale $t_{\mathrm {chem}}^{} = \delta _{\mathrm {th}}^{} / S_{\mathrm {L}}^{0}$. An example of the instantaneous flame front, temperature field, and velocity field after running the simulation for one initial eddy turn over time for case H is shown in Fig. 2. The interested reader is referred to Zirwes et al. (2020) for detailed information on the cubic discretization scheme implemented in OpenFOAM, and Tekgül et al. (2021) and Morev et al. (2022) for comprehensive explanation about the customized solver employed here.

Fig. 2

An instantaneous snapshot of flame front, temperature field, and velocity field after running the simulation for one initial eddy turn over time for case H. Fresh and burned gases are on the left and right sides of the domain, respectively

The accuracy of using different discretization schemes implemented in OpenFOAM on 3-D Taylor-Green vortex was tested/validated by Zirwes et al. (2020) on a cube with 512$^3$ cells. This problem is re-visited and discussed in Appendix A focusing on implementing the cubic discretization scheme for a cube with different cell sizes changing from 64$^3$, 128$^3$, 256$^3$, and 512$^3$ cells. Furthermore, the results of homogeneous isotropic turbulence under a decaying condition using the cubic discretization scheme are discussed in Appendix B.

3 Results and Discussion

3.1 Thermal Flame Thickness

In this section, the effect of turbulence intensity on local flame thickness of premixed turbulent methane/air and ammonia/air flames will be investigated. An overview of flame thinning/broadening could be observed in the temperature based progress variable slice presented in Fig. 3, where the progress variable is equal to $c = (T - T_{\mathrm {r}}^{})/(T_{\mathrm {ad}}^{} - T_{\mathrm {r}}^{})$, where T is the instantaneous temperature. It is observed that the flame surface is convoluted due to the turbulent structures, and the preheat zone region seems to be widened, while the reaction zone region retains its laminar-like structure. To quantify this effect, the thickening factor $\Theta (c)$ which has been previously introduced in Aspden et al. (2016) is utilized to evaluate the effect of turbulence intensity on local flame thinning/broadening. This parameter is defined for fuel/air mixture with the global Lewis number equal to unity using the following relation (Aspden et al. 2016):

$$\begin{aligned} \Theta (c) = \frac{\left\langle \nabla c (c)_{c = c^{*}} \right\rangle _{\mathrm {L}}}{\left\langle \nabla c (c)_{c = c^{*}} \right\rangle _{\mathrm {T}}}, \end{aligned}$$

(1)

where $c^{*}$ is the location where conditioning takes place, and $\mathrm {L}$ and $\mathrm {T}$ correspond to laminar and turbulent conditions, respectively. Figure 4 shows the variation of the heat release rate with respect to the progress variable for an unstrained premixed laminar methane/air flame with an equivalence ratio of 0.7. The heat release rate has been normalized by its maximum value. The preheat zone is located in the region where the progress variable is changing approximately from 0.03 to 0.55, and the reaction zone is where the progress variable is between 0.55 and 0.92. It should be mentioned that in the preheat zone region, the low-temperature radicals such as $\mathrm {CH_{3}O}$ for methane/air flames and $\mathrm {H_{2}O_{2}}$ for ammonia/air flames reach to their maximum values, while the heat release rate reaches to its maximum value in the reaction zone region. In order to be consistent with some of the previous studies in the literature, the preheat and reaction zone regions are defined as where $c^{*}$ is equal to 0.2 and 0.8, respectively. The preheat and reaction zone thicknesses have comparable values when conditioning the data at $c^{*}=0.2$ and 0.8, respectively. It should be added that the selected progress variable values for defining the preheat and reaction zone regions were previously utilized by Yu et al. (2019) in their recent work. However, different researchers selected different progress variable values when defining the preheat and reaction zone layers. For example, Dinkelacker et al. (1998) specified the preheat and reaction zones as where $c^{*}$ is equal to 0.25 and 0.50, respectively. Yuen and Gülder (2009) defined the preheat zone and reaction zone regions as where the progress variable is equal to 0.3 and 0.5, respectively. Nilsson et al. (2018) defined the preheat zone region as where $c^{*}=0.10$, and Attili et al. (2017) characterized the preheat and reaction zones at 900 K and 1800 K for their lean preheated methane/air flames under high pressure condition. The corresponding progress variable values were then 0.07 and 0.7 for the preheat and reaction zones, respectively. Sabelnikov et al. (2019) defined the local reaction zone thickness as where $c^{*}$ is connected with the maximum reaction rate value. The probability density functions (PDFs) of the preheat zone and reaction zone thickening factors for case H are shown in Fig. 5. It has been observed that the shape of each PDF is similar to a log-normal distribution. The highest probability of thickening factors for both regions are around one, while the preheat zone thickening factor seems to be skewed towards a higher value than unity. Experimental and numerical results reported in Dinkelacker et al. (1998); Yuen and Gülder (2009); Tamadonfar and Gülder (2015); Sabelnikov et al. (2019) showed a similar observation for PDFs of local flame thickness.

Fig. 3

Slice of a progress variable contour for case H captured at $t^{+} = 2$, where $t^{+}$ is the ratio of time to initial eddy turn over time

Fig. 4

Heat release rate with respect to the progress variable for a 1-D premixed laminar methane/air flame with an equivalence ratio of 0.7. The unburned, preheat, reaction, and hot product zones are shown herein. The heat release rate has been normalized by its maximum value

In order to understand the role of turbulence intensity on flame thickness, the ensemble average values of thickening factors for both methane/air and ammonia/air flames over different progress variable values are presented in Fig. 6. Results indicate that the flame fronts tend to be slightly widened in the preheat zone regions suggesting a moderate deviation from the thin laminar flamelet assumption for all the test cases. This deviation is more noticeable when evaluating the thickening factors at $c^{*}=0.05-0.10$. This might be due to the turbulent mixing resulting in strong temperature gradient in the leading edge of the flame. However, the turbulence structures might dissipate in the reaction zone region due to the existence of large heat release (Sankaran et al. 2007). Therefore, the reaction zone region is not disrupted for the test cases examined in this work. Previous experimental results of Tamadonfar and Gülder (2015) and Wabel et al. (2017) along with recent DNS studies of Aspden et al. (2016, 2019) showed a similar trend for methane/air flames. However, Aspden et al. (2019) showed that the thickening factor in the preheat zone region is less than unity for premixed turbulent hydrogen/air flames even at a low turbulent Karlovitz number test condition. Taking this observation into account, Aspden et al. (2016) concluded that the response of thickening factor to the turbulence intensity is dependent on the global Lewis number. Following their conclusion, the results presented in Fig. 6 show that the response of thickening factor to the turbulence intensity is the same for both fuel/air mixtures studied in this work since the global Lewis numbers are equal to unity value.

Fig. 5

Probability density functions of the preheat zone and reaction zone thickening factors for case H. The statistics are captured at $t^{+} = 2$

As mentioned earlier in Introduction, carrying out experimental measurements to evaluate 3-D scalar gradient of flame temperature for evaluating the thickening factor is an extremely challenging task. Only a few experimental measurements have succeeded in estimating the 3-D local flame thickness, see e.g. de Goey et al. (2005). The current quasi-DNS database is utilized to estimate the 2-D local thickening factor by excluding the temperature gradient in the third direction for evaluating $\Theta (c)$, Fig. 6. The 2-D local thickening factors overestimate the 3-D values by about 20–25%. This observation suggests that when detecting thick flames from 2-D imaging measurements, the experimental data should be interpreted attentively. Chakraborty and Hawkes (2011) proposed a correction factor of $4/\pi$ for correcting the 2-D flame surface density data into its actual 3-D value. Following their proposal, since the local temperature gradient is directly utilized to estimate the flame surface density and thickening factor, the proposed correction value could also be utilized to correct the 2-D local thickening factor estimation to its actual 3-D value. It is required to bear in mind that the proposed correction factor is based on isotropic condition of the angle between local flame surface and normal vectors on the measurement plane. Figure 6 shows that implementing the correction factor of $4/\pi$ to the gradient of 2-D temperature based reaction progress variable results in having an excellent quantitative agreement between the corrected 2-D local thickening factor and its actual 3-D value. The accuracy of this agreement increases significantly with increasing turbulent Reynolds number improving the isotropy assumption (Chakraborty and Hawkes 2011).

Fig. 6

Ensemble average values of thickening factors for a case A, b case B, c case C, d case F, e case G, and f case H with respect to the progress variable. The statistics are captured at $t^{+} = 2$

3.2 Strain Rate and Curvature Dependence of Thermal Flame Thickness

One of the potential mechanisms behind flame thinning/broadening is the flow straining as discussed in Chakraborty (2007); Sankaran et al. (2007); Driscoll (2008). Therefore, examining the effect of strain rate on thickening factor is quite helpful in understanding this mechanism as suggested by Dinkelacker et al. (1998). The normal strain rate, $a_{\mathrm {N}}^{}$, and tangential strain rate, $a_{\mathrm {T}}^{}$, are evaluated using the following formulations (Candel and Poinsot 1990):

$$\begin{aligned} a_{\mathrm {N}}^{} = n_{j}^{}n_{k}^{}\left( \frac{\partial u_{j}^{}}{\partial x_{k}^{}} \right) , \end{aligned}$$

(2)

$$\begin{aligned} a_{\mathrm {T}}^{} = \nabla \cdot\vec {u} - a_{\mathrm {N}}^{}= \left( \frac{\partial u_{j}^{}}{\partial x_{j}^{}} \right) - n_{j}^{}n_{k}^{}\left( \frac{\partial u_{j}^{}}{\partial x_{k}^{}} \right) , \end{aligned}$$

(3)

where $\nabla .\vec {u}$ is the dilatation rate, $u_{j}^{}$ is the jth component of the velocity vector, and $n = -\nabla c / \mid \nabla c \mid$ is the normal vector on a specified progress variable pointing towards the reactants. The normal and tangential strain rates, and the dilatation rate across the flame brush for cases A, C, F, and H are shown in Fig. 7. Results show that the mean value of tangential strain rate is approximately constant across the flame brush. A similar observation was previously reported in Han and Huh (2008) and Dopazo et al. (2015) based on 3-D DNS using simple chemistry. It is clear that if the tangential strain rate (dilatation rate) prevails the dilatation rate (tangential strain rate), the compressive (extensive) normal strain rate is promoted. Furthermore, an increase (decrease) in the mean value of tangential (normal) strain rate is observed with increasing turbulence intensity, see e.g. Figure 7a, b, c and d. The mean value of the flow dilatation rate is positive across the flame brush, and it reaches a peak value inside the flame. This could be due to severe heat release across the flame brush (Dopazo et al. 2015). In their recent work, Kasten et al. (2021) reported a similar behavior using simple chemistry under forced turbulence in the reactants. A similar trend was previously reported in Dopazo et al. (2015) for a flame under decaying turbulence in the corrugated flamelets regime with an Arrhenius one-step chemistry. Results show that the maximum value of the dilatation rate for ammonia/air flames (Fig. 7c, d) is higher than the corresponding value for methane/air flames (Fig. 7a, b). This could be due to the existence of higher heat release factor $\tau$ for stoichiometric ammonia/air flames than lean methane/air flames, where $\tau$ is evaluated as $(T_{\mathrm {ad}} - T_{\mathrm {r}})/T_{\mathrm {r}}$. Since the mean value of tangential strain rate is approximately constant, and the mean dilatation rate value is higher for ammonia/air flames compared to methane/air flames under constant $u_{}^{\prime }/S_{\mathrm {L}}^{0}$, the mean normal strain rate value is higher for ammonia/air flames. This means that iso-scalar lines widening, i.e. flame thickening, due to normal strain rate is more prominent for ammonia/air flames than methane/air flames.

Fig. 7

Normalized mean values of the normal strain rate, tangential strain rate, and dilatation rate across the flame brush for a case A, b case C, c case F, and d case H. Each parameter has been normalized by the chemical time scale. The statistics are captured at $t^{+} = 2$

Figure 8 shows the joint PDFs between the normalized tangential (normal) strain rate and preheat/reaction zone thickening factors for cases F and H. The tangential and normal strain rates have been normalized by the chemical time scale. Results show that the preheat (reaction) zone thickening factor has a strong (weak) negative correlation with the tangential strain rate for ammonia/air flames at low to medium Karlovitz numbers. In addition, increasing the turbulence intensity results in having a broader thickening factor distribution being more pronounced in the preheat zone region than the reaction zone region. Another important observation is that a notable flame thinning happens under the action of positive local tangential strain rate across the flame brush with increasing turbulence intensity. This observation is consistent with earlier findings, see e.g. Chakraborty (2007) and Sankaran et al. (2007).

Results in Fig. 8e, g show that the local normal strain rate is mostly positive for low Karlovitz number case, i.e. case F. This indicates that the extensive normal strain rate is mainly responsible for widening the iso-scalar lines for this test condition. However, increasing $u_{}^{\prime }/S_{\mathrm {L}}^{0}$ results in inducing compressive normal strain rate across the flame which is more pronounced in the preheat zone region than the reaction zone region, Fig. 8f, h. It should be emphasized that the flame front then tends to become thinner under the action of compressive normal strain rate. The authors confirm that the above observation exists for methane/air flames examined in this study as well, and the joint PDFs are not presented here for conciseness.

Fig. 8

Joint PDFs between the normalized tangential strain rate and the preheat (reaction) zone thickening factor for a, c case F, and b, d case H. Joint PDFs between the normalized normal strain rate and the preheat (reaction) zone thickening factor for e, g case F, and f, h case H. The tangential (normal) strain rate has been normalized by the chemical time scale. The statistics are captured at $t^{+} = 2$

It is of particular importance to examine the dependency of flame front curvature on both strain rate components, i.e. tangential and normal terms, to deepen our comprehension on the role of these components on flame thickening. The flame front curvature, $\kappa$, is evaluated using the following formulation presented in Trouvé and Poinsot (1994) and Han and Huh (2008):

$$\begin{aligned} \kappa = 0.5 \times (\nabla \cdot n)_{c = c^{*}}. \end{aligned}$$

(4)

It should be mentioned that in this definition, the positive (negative) flame front curvature is convex (concave) towards the reactants (Trouvé and Poinsot 1994). The joint PDFs between the normalized tangential strain rate and the normalized flame front curvature, conditioned at $c^{*} = 0.2$ and $c^{*} = 0.8$, for cases F and H are shown in Fig. 9a–d. The flame front curvature has been normalized by the laminar flame thickness. Results show that tangential strain rate is negatively correlated with the flame front curvature. This finding was previously observed in Chakraborty and Cant (2004) and Chakraborty (2007) based on 3-D DNS with simple chemistry. This observation confirms that the incidence of negative tangential strain rate is more likely at positive flame front curvature. Chakraborty (2007) discussed that the defocusing of heat, flow divergence, and the acceleration of flow at positively curved regions results in lower $a_{\mathrm {T}}^{}$ value at these regions compared to the negatively curved regions. Furthermore, increasing $u_{}^{\prime }/S_{\mathrm {L}}^{0}$ results in having flame front elements with smaller radii experiencing higher values of tangential strain rate. This is more noticeable for the preheat zone region than the reaction zone region where the heat release might alleviate the flame front wrinkling.

The joint PDFs between the normalized normal strain rate and the normalized flame front curvature for cases F and H are presented in Fig. 9e–h. Negative correlation between normal strain rate, conditioned at $c^{*} = 0.2$, and flame front curvature could be observed for both cases. This observation indicates that at highly negative curved regions, the widening of iso-scalar lines due to the normal strain rate is more pronounced in comparison to the positively curved regions. However, there exists no clear correlation between the normal strain rate, conditioned at $c^{*} = 0.8$, and the flame front curvature for the studied cases indicating that the thickening of flame due to this mechanism is not sensitive to flame front curvature for the reaction zone region. As discussed earlier, increasing $u_{}^{\prime }/S_{\mathrm {L}}^{0}$ results in decreasing the normal strain rate significantly as well as widening the distribution implying that the compressive strain rate could be responsible for bringing the iso-scalar lines together, i.e. thinning the flame. For the sake of brevity, results for methane/air flames are not presented here.

Fig. 9

Joint PDFs between the normalized tangential strain rate and the flame front curvature for a case F conditioned at $c^{*}=0.2$, b case F conditioned at $c^{*}=0.8$, c case H conditioned at $c^{*}=0.2$, and d case H conditioned at $c^{*}=0.8$. Joint PDFs between the normalized normal strain rate and the flame front curvature for e case F conditioned at $c^{*}=0.2$, f case F conditioned at $c^{*}=0.8$, g case H conditioned at $c^{*}=0.2$, and h case H conditioned at $c^{*}=0.8$. The tangential and normal strain rates have been normalized by the chemical time scale, and the flame front curvature has been normalized by the laminar thermal flame thickness. All the statistics are captured at $t^{+} = 2$

To fully capture the underlying mechanism behind flame thinning/broadening over the flame front curvature, one should assess the influence of flame self-propagation on flame thickness as well since the flow straining is not the only active mechanism for this purpose. As mentioned earlier, the balance equation of $\mid \nabla c \mid$ was examined by Sankaran et al. (2007) to recognize the role of flame front curvature on flame thickness of premixed turbulent methane/air Bunsen-type flames by investigating the turbulent flow straining and flame self-propagation terms. Following their methodology helps to comprehend the mechanism behind flame thinning/broadening in more depth for premixed turbulent planar ammonia/air and methane/air flames as well. According to Sankaran et al. (2007), the balance equation of $\mid \nabla c \mid$ is written as follows:

$$\begin{aligned} \frac{D}{Dt}\mid \nabla c \mid = -\left( a_{\mathrm {N}}^{} + \frac{\partial S_{\mathrm {d}}^{} }{\partial n} \right) \mid \nabla c \mid , \end{aligned}$$

(5)

where $S_{\mathrm {d}}^{}$ is the flame displacement speed. As described by Gran et al. (1996) and Echekki and Chen (1999), the flame displacement speed could then be decomposed into its reaction $S_{\mathrm {d, r}}^{} = {\dot{\omega }} / \rho \mid \nabla c \mid$, normal diffusion $S_{\mathrm {d, n}}^{} = n.\nabla (\rho D n. \nabla c) / \rho \mid \nabla c \mid$, and tangential diffusion $S_{\mathrm {d, t}}^{} = -2D\kappa$ components where ${\dot{\omega }}$ is the reaction rate of the reaction progress variable, $\rho$ is the density, and D is the progress variable diffusivity. Sankaran et al. (2007) discussed that a positive (negative) valued term on the R.H.S. of Eq. 5 operates as a source (sink) for $\mid \nabla c \mid$ resulting in flame thinning (broadening). The budget terms acting on preheat zone and reaction zone regions with respect to the normalized flame front curvature are shown in Fig. 10. All the source terms have been normalized by the chemical time scale, and the flame front curvature has been normalized with the laminar thermal flame thickness. Results show that positive values of normal strain rate act as a sink of $\mid \nabla c \mid$ for the preheat zone region, and it is negatively correlated with the flame front curvature, see Fig. 10a, b, e and f. In the reaction zone region, the normal strain rate does not show any obvious correlation with the flame front curvature for both methane/air and ammonia/air flames indicating that the thickening effect induced by the normal strain rate is independent of the flame front curvature at this region. Furthermore, since the local normal strain rate value for ammonia/air flame is higher than the corresponding value for methane/air flame under constant $u_{}^{\prime }/S_{\mathrm {L}}^{0}$ and $\kappa \delta _{\mathrm {th}}$, the thickening effect induced by normal strain rate is higher for ammonia/air flame than methane/air flame. In addition, increasing $u_{}^{\prime }/S_{\mathrm {L}}^{0}$ for each fuel/air mixture results in decreasing the mean normal strain rate value indicating that the thickening effect induced by normal strain rate decreases.

Furthermore, the kinematic restoration $\partial S_\mathrm {{d, r+n}}^{}/ \partial n$ acts as a source of $\mid \nabla c \mid$ for the preheat zone region resulting in flame thinning for all the test cases studied here. Since the absolute value of $\partial S_\mathrm {{d, r+n}}^{}/ \partial n$ under a constant $u_{}^{\prime }/S_{\mathrm {L}}^{0}$ is higher for ammonia/air flame than the corresponding value for methane/air flame, the preheat zone thinning induced by kinematic restoration is more pronounced for ammonia/air flame than methane/air flame. Furthermore, increasing $u_{}^{\prime }/S_{\mathrm {L}}^{0}$ weakens the thinning effect induced by kinematic restoration slightly for each fuel/air mixture. It has been observed that in the reaction zone region, the effect of kinematic restoration on flame thinning is significant (insignificant) for ammonia/air (methane/air) flames with negative flame front curvature, while its effect is more prominent for positive flame front curvature regions. A further observation reveals that $\partial S_\mathrm {{d, t}}^{}/ \partial n$ primarily acts as a sink for highly positive curved regions, and it is positively correlated with the flame front curvature. A similar observation was previously reported for highly positive curved regions for Bunsen-type burner in Sankaran et al. (2007). With increasing $u_{}^{\prime }/S_{\mathrm {L}}^{0}$, the local value of $\partial S_\mathrm {{d, t}}^{}/ \partial n$ under a constant $\kappa$ value increases indicating the thickening effect activated by $\partial S_\mathrm {{d, t}}^{}/ \partial n$ improves. In addition, the thickening effect caused by this term is higher for ammonia/air flame than methane/air flame under a constant $u_{}^{\prime }/S_{\mathrm {L}}^{0}$ value. In general, the net balance of all source terms reported in Eq. 5 results in flame thickening for the preheat zone region for both methane/air and ammonia/air flames. For the reaction zone region, the net balance of all terms is positively correlated with the flame front curvature for both fuel/air mixtures. It should be noted that in their recent work, Wang et al. (2017) showed negative (positive) values of $a_{\mathrm {N}}^{}$ ($\partial S_{\mathrm {d}}^{} / \partial n$) for all progress variable magnitudes for a premixed methane/air jet flame under extremely high turbulent Karlovitz number. In addition, Nilsson et al. (2018) reported that the kinematic restoration has the thickening effect, while the normal strain rate has the thinning effect for lean premixed turbulent planar methane/air flames under high Karlovitz number regime.

Fig. 10

Budget of source terms in Eq. 5 and scatter of $\left( a_{\mathrm {N}}^{} + \partial S_{\mathrm {d}}^{}/\partial n \right)$ as a function of normalized flame front curvature. The data are conditioned at $c^{*}=0.2$ for a case A, b case C, e case F, and f case H, and conditioned at $c^{*}=0.8$ for c case A, d case C, g case F, and h case H. The source terms have been normalized by the chemical time scale, and the flame front curvature has been normalized by the laminar thermal flame thickness

3.3 Flame Front Orientation Characteristics

The flame front orientation statistics have been shown to influence the behavior of normal strain rate (Kasten et al. 2021). Since the latter term has a significant effect on flame thinning/broadening, according to Eq. 5, it is of great importance to study the flame front orientation statistics as well to deepen our perception about the cases examined in this study. According to Kasten et al. (2021), $\theta _{i} = \mathrm {cos}^{-1}({\hat{e}}_{i}.n)$ is the angle between normal to the flame and the corresponding eigenvectors of the strain rate tensor for $i=1$, 2, and 3, where i corresponds to the principle strain rate directions. The mean values of $\mid \mathrm {cos}(\theta _{1}) \mid$, $\mid \mathrm {cos}(\theta _{2}) \mid$, and $\mid \mathrm {cos}(\theta _{3}) \mid$ with respect to the progress variable is presented in Fig. 11a, b for methane/air and ammonia/air flames, respectively. Results show that n and ${\hat{e}}_{1}$, i.e. the most extensive principle strain rate eigenvector, are primarily oriented with each other for low turbulence intensity flame conditions, that is, cases A and F. This alignment is more pronounced near $c=0.6$, where the thermal expansion due to heat release is more dominant across the flame brush. Furthermore, increasing the turbulence intensity from $u_{}^{\prime }/S_{\mathrm {L}}^{0} = 4$ to 8 reduces this trend, while the alignment between n and ${\hat{e}}_{2}$ (${\hat{e}}_{3}$) improves across the flame brush. The current observation is consistent with previous findings reported in Chakraborty and Swaminathan (2007) and Kasten et al. (2021) discussing the dominant alignment of n and ${\hat{e}}_{1}$ for flames with $\tau Da \gg 1.0$.

To comprehend the orientation of flame fronts with respect to the flame front curvature, Fig. 11c, d shows the mean values of $\mid \mathrm {cos}(\theta _{i}) \mid$ with respect to the normalized flame front curvature for ammonia/air flames conditioned at $c^{*} = 0.2$ and 0.8, respectively. Results show that, for the preheat zone region, the alignment of n and ${\hat{e}}_{1}$ (${\hat{e}}_{2}$ and ${\hat{e}}_{3}$) decreases (increases) with increasing the flame front curvature, and the differences between $\mid \mathrm {cos}(\theta _{i}) \mid$ values decreases near highly positive curved regions. This observation suggests that the collinear alignment between n and ${\hat{e}}_{1}$ is more pronounced near the negative curved regions where the focusing of heat is more noticeable than the regions with positive flame front curvature. Furthermore, n aligns collinearly with ${\hat{e}}_{1}$ in the reaction zone region where the influence of heat release is quite dominant, see Fig. 11d. A further observation indicates that the mean values of $\mid \mathrm {cos}(\theta _{i}) \mid$ become roughly constant with increasing turbulence intensity in the reaction zone region indicating that $\mid \mathrm {cos}(\theta _{i}) \mid$ becomes insensitive to the flame front curvature.

Fig. 11

Mean values of $\mid \mathrm {cos}(\theta _{i}) \mid$ with respect to the progress variable for a cases A–C, and b cases F–H. Mean values of $\mid \mathrm {cos}(\theta _{i}) \mid$ with respect to the normalized flame front curvature for ammonia/air flames conditioned at the c preheat zone region, and the d reaction zone region. The flame front curvature has been normalized by the laminar flame thickness. All the statistics are captured at $t^{+} = 2$

3.4 Assessing the Validity of Damköhler’s First Hypothesis

As previously discussed, the assessment of Damköhler’s first hypothesis validity is the subject of ongoing research since the turbulent burning velocity and flame surface area are basically connected through this hypothesis (Damköhler 1940; Peters 1986, 1999), and it is of particular importance to examine the validity of this hypothesis for novel fuels. Following the formulation presented in Chakraborty et al. (2019), the ratio of the turbulent burning velocity $S_{\mathrm {T}}^{}$ to the so-called modified burning velocity (mean local flamelet consumption velocity according to Tamadonfar and Gülder 2015) $S_{\mathrm {L}}^{'}$ is defined as follows:

$$\begin{aligned} \frac{S_{\mathrm {T}}^{}}{S_{\mathrm {L}}^{'}} = \frac{\int _{V}{\dot{\omega }}dV/\rho _{0}A_{\mathrm {L}}}{\int _{V}\overline{(\rho S_{\mathrm {d}})}_{s}^{}\sum _{\mathrm {gen}}dV/\rho _{0}\int _{V}\sum _{\mathrm {gen}}dV} = \frac{A_{\mathrm {T}}^{}}{A_{\mathrm {L}}^{}}, \end{aligned}$$

(6)

where $\rho _{0}$ is the reactants density, $\sum _{\mathrm {gen}} = \overline{\mid \nabla c \mid }$ is the flame surface density, $A_{\mathrm {T}}^{} = \int _{V}\mid \nabla c \mid dV$ is the wrinkled flame surface area, and $A_{\mathrm {L}}^{}$ is the cross-sectional area normal to the direction of mean flame propagation, and $\overline{(\rho S_{\mathrm {d}})}_{s}^{} = \overline{(\rho S_{\mathrm {d}})\mid \nabla c \mid } / \overline{\mid \nabla c \mid }$ is the surface-averaged value of $\rho S_{\mathrm {d}}$. An alternative method for evaluating the wrinkled flame surface area $A_{\mathrm {T}}^{}$ is to estimate the reaction zone iso-surface area by excluding the contribution from the preheat zone region. It has been shown by Boger et al. (1998) that $\overline{(\rho S_{\mathrm {d}})}_{s}^{}$ is approximately equal to $\rho _{0} S_{\mathrm {L}}^{0}$ by assuming unity Lewis number for all species using simple chemistry. This results in reducing Eq. 6 to Damköhler’s first hypothesis relation as follows:

$$\begin{aligned} \frac{S_{\mathrm {T}}^{}}{S_{\mathrm {L}}^{0}} = \frac{A_{\mathrm {T}}^{}}{A_{\mathrm {L}}^{}}. \end{aligned}$$

(7)

It should be emphasized that Chakraborty et al. (2019) showed that Eq. 7 relation retains its validity for statistically premixed turbulent planar flames where the mean flame front curvature disappears and the correlation between local flame front curvature and $\mid \nabla c \mid$ is insignificant. More recently, Klein et al. (2020) showed that the validity of Eq. 7 might also be dependent on the choice of species utilized for defining the reaction progress variable as well. In this work, the temperature field is utilized to evaluate the parameters presented in Eqs. 6 and 7. In addition, since the 3-D experimental measurements of temperature gradient for evaluating the wrinkled flame surface area is a very challenging task, it is very helpful to examine the evaluation of flame surface area based on 2-D projection of actual quasi-DNS temperature data as well. This methodology helps in interpreting the planar measurements obtained from experiments more accurately. Figure 12 shows the time evolution of non-dimensional turbulent burning velocity and flame surface area for all the conditions studied in this work. Results in Fig. 12 leads to the following observations. First, increasing the turbulence intensity results in enhancing the turbulent burning velocity. For example, $S_{\mathrm {T}}^{}/S_{\mathrm {L}}^{0}$ increases from 1.39 to 1.84 with increasing $u_{}^{\prime }/S_{\mathrm {L}}^{0}$ from 4 to 8 for ammonia/air flames at $t_{}^{+} = 2$, where $t^{+}$ is the ratio of time to initial eddy turn over time. A similar trend was previously reported in the literature, see e.g. Han and Huh (2008); Tamadonfar and Gülder (2014, 2015); Aspden et al. (2016); Nivarti and Cant (2017); Brearley et al. (2020). It should be mentioned that the ratio of the turbulent burning velocity to the unstrained premixed laminar burning velocity values are low when they are compared with the available experimental data in the literature under the same turbulence intensity. However, it should be emphasized that direct comparison between the turbulent burning velocity values obtained from the current quasi-DNS results and the experimental data might not be straightforward due to various reasons. One of the reasons for this difference could be the existence of periodic boundary conditions in the numerical test cases in which these conditions do not happen during the experiments (Driscoll et al. 2020). Furthermore, it is recommended to compare the results obtained from one flame category (envelop, oblique, flat, and spherical flames) to that particular category since each flame category has different boundary conditions and flame wrinkling process (Bilger et al. 2005; Driscoll 2008). This suggests that comparing the burning velocity values obtained from the current setup with the corresponding values measured from the experimental setups might be impractical. Second, the results show that the magnitudes of $S_{\mathrm {T}}^{}/S_{\mathrm {L}}^{0}$ for methane/air and ammonia/air flames are approximately equal over time under a constant $u_{}^{\prime }/S_{\mathrm {L}}^{0}$ value. Third, the magnitudes of $A_{\mathrm {T}}^{}/A_{\mathrm {L}}^{}$ enhances with increasing $u_{}^{\prime }/S_{\mathrm {L}}^{0}$. This is due to an enhancement in flame front wrinkling originated from turbulence structures (Brearley et al. 2020). DNS results of planar flames using simple/detailed chemistry reported in Han and Huh (2008); Nivarti and Cant (2017); Brearley et al. (2020); Song et al. (2021), and the experimental work performed on Bunsen-type burner by Tamadonfar and Gülder (2015) showed a similar trend. Furthermore, the magnitude of flame surface area derived from 2-D projection of temperature data reduces compared to its actual 3-D value since the temperature gradient data in one direction is absent in the evaluation. To improve the aforementioned observation, the correction factor of $4/\pi$ proposed in Chakraborty and Hawkes (2011) is utilized to estimate the actual flame surface area from 2-D projection of temperature data. For all the test conditions, the corrected flame surface area overpredicts the actual 3-D flame surface area. This overprediction, after one initial eddy turn over time, improves with increasing turbulent Reynolds number. Therefore, using the aforementioned correction factor for estimating the flame surface area gives a reasonable result under high turbulent Reynolds number, while it does not provide suitable results for low turbulent Reynolds number cases, i.e. cases A and F. It should be borne in mind that this correction factor was developed under isotropy assumption (Chakraborty and Hawkes 2011). The successful implementation of using this correction factor for evaluating the flame surface area from 2-D projection data was previously reported for freely-propagating premixed turbulent planar hydrogen/air flames (Klein et al. 2020) and for premixed turbulent Bunsen-type flames (Chakraborty et al. 2019). Finally, it should be mentioned that evaluating the wrinkled flame surface area from the volume integration of the flame surface density or direct measurement of the reaction zone iso-surface area lead to almost similar results after about one initial eddy turn over time.

Fig. 12

Non-dimensional turbulent burning velocity and flame surface area with respect to normalized time for a case A, b case B, c case C, d case F, e case G, and f case H. The flame surface area was estimated from 3-D (volume integration of the flame surface density and direct measurement of the reaction zone iso-surface area), 2-D, and 2-D with a correction factor temperature data. The time has been normalized by the initial eddy turn over time

The ratio of $S_{\mathrm {T}}^{}/S_{\mathrm {L}}^{0}$ to $A_{\mathrm {T}}^{}/A_{\mathrm {L}}^{}$ with respect to time is presented in Fig. 13 for all the test conditions. It should be noted that in this figure, the wrinkled flame surface area is evaluated by integrating the flame surface density over the volume. Results show that $S_{\mathrm {T}}^{}/S_{\mathrm {L}}^{0}$ and $A_{\mathrm {T}}^{}/A_{\mathrm {L}}^{}$ are almost equal to each other after about one initial eddy turn over time. This observation suggests the validity of Damköhler’s first hypothesis, Eq. 7, for both methane/air and ammonia/air flames examined in this work. This finding is in agreement with previous findings reported in Han and Huh (2008); Chakraborty et al. (2014); Nivarti and Cant (2017); Ahmed et al. (2019); Sabelnikov et al. (2019); Brearley et al. (2020); Klein et al. (2020); Varma et al. (2021); Song et al. (2021). However, the validity of this hypothesis might not be maintained for premixed turbulent planar flames with a global Lewis number smaller or higher than unity (Haworth and Poinsot 1992; Trouvé and Poinsot 1994; Bell et al. 2007; Han and Huh 2008; Chakraborty et al. 2014), and for flames with significantly non-zero mean flame front curvature, i.e. Bunsen-type flames, see, DNS results of Chakraborty et al. (2019), and the experimental works of Daniele et al. (2013) and Tamadonfar and Gülder (2015). Haworth and Poinsot (1992), Trouvé and Poinsot (1994), Bell et al. (2007), Han and Huh (2008), and Chakraborty et al. (2014) showed that the ratio of $S_{\mathrm {T}}^{}/S_{\mathrm {L}}^{0}$ to $A_{\mathrm {T}}^{}/A_{\mathrm {L}}^{}$ is higher (smaller) than unity for premixed turbulent planar flames with a global Lewis number smaller (larger) than unity. For example, Chakraborty et al. (2014) reported that the ratio of $S_{\mathrm {T}}^{}/S_{\mathrm {L}}^{0}$ to $A_{\mathrm {T}}^{}/A_{\mathrm {L}}^{}$ is around 3.5 (0.85) when $Le = 0.34$ ($Le = 1.20$) using DNS with simple chemistry. It could be concluded from Fig. 13 that the validity of Damköhler’s first hypothesis might not be dependent on the fuel type, and it might be more related to the global Lewis number of the mixture. A further analysis shows that the ratio of $S_{\mathrm {T}}^{}/S_{\mathrm {L}}^{0}$ to $A_{\mathrm {T}}^{}/A_{\mathrm {L}}^{}$ in which $A_{\mathrm {T}}^{}$ is evaluated based on the 2-D projection of temperature data with a correction factor underpredicts the actual ratio for low turbulent Reynolds number test conditions, i.e. cases A and F. However, this underprediction improves drastically with increasing turbulent Reynolds number. This observation indicates that when using the aforementioned correction factor for getting an estimation of actual flame surface area from two-dimensional measurement, the results need to be interpreted cautiously especially for the test conditions with low turbulent Reynolds number.

Fig. 13

The ratio of the non-dimensional turbulent burning velocity to the non-dimensional flame surface area with respect to time for a cases A and F, b cases B and G, and c cases C and H. The flame surface area was evaluated using 3-D temperature data and 2-D temperature data with a correction factor. The time has been normalized by the initial eddy turn over time

Table 2 $S_{\mathrm {L}}^{\prime }/ S_{\mathrm {L}}^{0}$ values for all the studied conditions

According to Chakraborty et al. (2019), the ratio of the modified burning velocity to the unstrained laminar burning velocity $S_{\mathrm {L}}^{\prime }/ S_{\mathrm {L}}^{0}$ is close to unity for statistically premixed turbulent planar flames with unity Lewis number ensuring the validity of Damköhler’s first hypothesis. Therefore, evaluating the aforementioned parameter could be beneficial in confirming the validity of this hypothesis for the test conditions studied in this work as well. Table 2 shows $S_{\mathrm {L}}^{\prime }/ S_{\mathrm {L}}^{0}$ values for all the test conditions. Results show that this ratio is approximately equal to one for all the cases. In other words, the modified burning velocity $S_{\mathrm {L}}^{\prime }$ is almost equal to the unstrained premixed laminar burning velocity $S_{\mathrm {L}}^{0}$ implying that the mean local flamelet consumption velocity is approximately equal to the unstrained premixed laminar burning velocity for all the test cases. A similar trend was previously reported under forced turbulence by Ahmed et al. (2019) confirming the validity of Damköhler’s first hypothesis using DNS with simple chemistry. However, the current observation does not suggest that all local flame components burn similar to the unperturbed laminar flame. More recently, Song et al. (2021) investigated the local behavior of highly turbulent premixed hydrogen/air flames by examining the statistics of density-weighted flame displacement speed $S_{\mathrm {d}}^{*} = \rho S_{\mathrm {d}}^{} / \rho _{0}^{}$ using DNS. Following their methodology helps to comprehend the behavior of local flame front speed relative to the flow velocity for ammonia/air and methane/air flames in more depth. Figure 14a, b, c and d shows the PDFs of $S_{\mathrm {d}}^{*} / S_{\mathrm {L}}^{0}$ conditioned at $c^{*} = 0.2$ and 0.8 for methane/air and ammonia/air flames. Results show that the PDF distribution becomes wider with increasing $u_{}^{\prime }/S_{\mathrm {L}}^{0}$ from 4 to 8 for both regions and fuel/air mixtures. A similar observation was previously reported for hydrogen/air flames in Song et al. (2021). Observing negative $S_{\mathrm {d}}^{*} / S_{\mathrm {L}}^{0}$ values at the preheat zone region might be due to the strong turbulent mixing in this region. Furthermore, the probability of observing negative $S_{\mathrm {d}}^{*} / S_{\mathrm {L}}^{0}$ values are zero at the reaction zone region for both methane/air and ammonia/air flames due to the dissipation of turbulence structures at this region. A similar finding was recently reported for premixed turbulent hydrogen/air flames under extremely intense turbulence conditions (Song et al. 2021). To fully understand the behavior of $S_{\mathrm {d}}^{*}$ across the flame brush, mean values of $S_{\mathrm {d}}^{*} / S_{\mathrm {L}}^{0}$ with respect to the progress variable for methane/air and ammonia/air flames are presented in Fig. 14e, f, respectively. Results show that the mean values of $S_{\mathrm {d}}^{*} / S_{\mathrm {L}}^{0}$ are almost equal to unity across the flame brush except near the trailing edge of the flame indicating that most of ammonia/air and methane/air flame elements propagate identical to the laminar flame. Results clearly show that the probability of finding negative values for $S_{\mathrm {d}}^{*}$ decreases with increasing the progress variable.

Fig. 14

PDFs of the normalized density-weighted flame displacement speed for methane/air (ammonia/air) flames conditioned at the a, c preheat zone region, and b, d reaction zone region. Mean values of normalized density-weighted flame displacement speed with respect to the progress variable for e methane/air and f ammonia/air flames. The error bars indicate two standard deviations of the mean

4 Concluding Remarks

The instantaneous inner flame front structures and burning velocities of premixed turbulent ammonia/air and methane/air flames were investigated based on detailed-chemistry quasi-DNS results of statistically stationary premixed turbulent planar thermo-diffusively neutral flames under different turbulence levels. The Karlovitz number was systematically changed from 4.26 to 12.06. The simulations were performed using the open-source CFD toolbox OpenFOAM, and the recently developed open-source dynamic load balancing model (DLBFoam) was used to overcome the computational load imbalance between processors. It is worth noting that the numerical simulations with up to 85 million computational cells should be considered as very large-scale simulations in the OpenFOAM framework. Such reacting flow simulations have been scarcely documented in the literature, and they are completely enabled by the recent DLBFoam package. The key findings in this study are summarized as follows:

• 1. The probability density functions of the preheat zone and reaction zone thickening factors had the log-normal distribution. The ensemble average values of preheat zone thickening factor were slightly more than unity suggesting a moderate deviation from the thin laminar flamelet assumption for both ammonia/air and methane/air flames. In addition, mean values of reaction zone thickening factor were approximately equal to unity indicating that this region was not altered with turbulence. The correction factor of $4/\pi$ proposed in Chakraborty and Hawkes (2011) was utilized to correct the 2-D local flame thickness to its 3-D value giving satisfying results under high turbulent Reynolds number test conditions.

• 2. The tangential strain rate was approximately constant across the flame brush, and it increased with increasing turbulence intensity. In addition, the normal strain rate decreased with increasing turbulence intensity. The mean value of dilatation rate was positive across the flame brush, and it attained a peak value inside the flame brush. The mean normal strain rate and dilatation rate values for ammonia/air flame were higher than the corresponding values for methane/air flame under a constant $u_{}^{\prime }/S_{\mathrm {L}}^{0}$ value due to possibly higher heat release factor.

• 3. In general, the tangential strain rate was negatively correlated with the thickening factor. The normal strain rate was mostly positive for low Karlovitz number cases resulting in flame front widening. Increasing $u_{}^{\prime }/S_{\mathrm {L}}^{0}$ induced compressive normal strain rate across the flame resulting in bringing the iso-scalar lines together, i.e. flame thinning.

• 4. The balance equation of $\mid \nabla c \mid$ was examined to explore the active mechanisms for flame thinning/broadening. In general, $a_{\mathrm {N}}^{}$ acted as a sink for $\mid \nabla c \mid$, and $\partial S_\mathrm {{d, r+n}}^{}/ \partial n$ operated as a source of $\mid \nabla c \mid$ for the conditions tested in this study. When the flame front curvature is positive, $\partial S_\mathrm {{d, t}}^{}/ \partial n$ acted as a sink for $\mid \nabla c \mid$. The differences between values of these terms were discussed and compared for ammonia/air and methane/air flames.

• 5. The flame front orientation statistics were analyzed for all the test cases. Results showed that n and ${\hat{e}}_{1}$ were mainly oriented with each other for low Karlovitz number test cases, and this tendency decreased with increasing $u_{}^{\prime }/S_{\mathrm {L}}^{0}$. These statistics were found to be independent on the fuel type. The dependency of flame front orientation statistics with the flame front curvature was discussed as well.

• 6. The turbulent burning velocity and wrinkled flame surface area of ammonia/air and methane/air flames were examined. Results showed that both of these parameters were increased with enhancing $u_{}^{\prime }/S_{\mathrm {L}}^{0}$. The correction factor proposed in Chakraborty and Hawkes (2011) for evaluating the actual flame surface area from 2-D projection of actual data gave satisfying results under high turbulent Reynolds number flame conditions.

• 7. The validity of Damköhler’s first hypothesis was examined for ammonia/air and methane/air flames. Results showed that the ratio of the turbulent burning velocity to the unstrained premixed laminar burning velocity $S_{\mathrm {T}}^{}/S_{\mathrm {L}}^{0}$ was approximately equal to the ratio of the wrinkled to the unwrinkled flame surface area $A_{\mathrm {T}}^{}/A_{\mathrm {L}}^{}$ ensuring the validity of this hypothesis. Further analysis revealed that $S_{\mathrm {L}}^{\prime }$ was almost equal to $S_{\mathrm {L}}^{0}$ approving the validity of Damköhler’s first hypothesis.

• 8. Mean values of density-weighted flame displacement speed $S_{\mathrm {d}}^{*}$ were almost equal to the unstrained premixed laminar burning velocity $S_{\mathrm {L}}^{0}$ for ammonia/air and methane/air flames under different Karlovitz numbers indicating that the bulk of flame elements propagated similar to the unstrained premixed laminar flame for the test conditions examined in this study.

Appendices

Appendix A Taylor-Green Vortex

One of the fundamental test cases for the assessment of turbulent flow simulation methodologies is the Taylor-Green vortex (T.G.V.), as mentioned by DeBonis (2013), which is a canonical problem in fluid mechanics. Examining the T.G.V. enables us to study the vortex dynamics, turbulent transition and its decay, and the process involved in energy dissipation (Taylor and Green 1937). Zirwes et al. (2020) examined the accuracy of different discretization schemes implemented in OpenFOAM using the T.G.V. test case. They argued that the cubic scheme could estimate the time progress of dissipation rate accurately, and the differences between peak value of the dissipation rate, evaluated directly from the kinetic energy, was less than 1% when comparing with the spectral DNS solution. However, they showed that using the second order central difference scheme underestimates the peak dissipation rate by about 5%. This problem is revisited here by simulating the T.G.V. case for different grid resolutions using the cubic discretization scheme and comparing dissipation rate results obtained from two methodologies with the 3-D pseudo-spectral reference solution (van Rees et al. 2011).

The flow conditions and the post-processing approach utilized in this section was initiated by the AIAA First International Workshop on High-Order Methods in Computational Fluid Dynamics (DeBonis 2013). This problem is defined on a tri-periodic cube with a side length equal to 2$\pi$. The initial distribution for the velocity and pressure fields are thoroughly described in DeBonis (2013) and Zirwes et al. (2020). According to the AIAA workshop, the Reynolds number was set to 1600. The grids are equally spaced in the Cartesian environment of 64$^3$, 128$^3$, 256$^3$, and 512$^3$ cells. The simulations are solved by using the pimpleFoam solver available in OpenFOAM-v8 by solving the incompressible Navier–Stokes equations.

The integrated kinetic energy, $E_{k}$, is derived using the following relation (DeBonis 2013):

$$\begin{aligned} E_{k} = \frac{1}{\rho _{0} V}\int _{V}\rho \frac{u_{i} u_{i} }{2}dV. \end{aligned}$$

(A1)

The kinetic energy could then be differentiated in time to evaluate the kinetic energy dissipation rate, $\epsilon (E_{k})$, using the following relation (DeBonis 2013):

$$\begin{aligned} \epsilon (E_{k}) = - \frac{\partial E_{k}}{\partial t}. \end{aligned}$$

(A2)

Furthermore, the kinetic energy dissipation rate could also be evaluated by integrating the enstrophy, $\epsilon (\zeta )$, through a constant using the following formulation (DeBonis 2013):

$$\begin{aligned} \epsilon (\zeta ) = 2 \nu \frac{1}{\rho _{0}V}\int _{V}\rho \frac{\omega _{i} \omega _{i} }{2}dV. \end{aligned}$$

(A3)

The evolution of kinetic energy over time for all the conditions and the pseudo-spectral solution of van Rees et al. (2011) are shown in Fig. 15a. Results indicate that the differences between kinetic energy is not significant for grids with 128$^3$, 256$^3$, and 512$^3$ cells compared to the reference solution. However, the domain with 64$^3$ cells contains less kinetic energy. This observation is consistent with previous findings reported in DeBonis (2013). The directly computed and enstrophy based dissipation rates over time are shown in Fig. 15b, c. It is worth mentioning that the values of $\epsilon (E_{k})$ and $\epsilon (\zeta )$ should be identical (DeBonis 2013). However, there is a large discrepancy between these two values for the solution of coarse grids, i.e. 64$^3$ and 128$^3$ cells, while this difference enhances significantly with refining the mesh size. The energy residual which is a practical quantity to represent this difference, as previously discussed in Vuorinen and Keskinen (2016), is shown in Fig. 15d.

Fig. 15

Normalized integrated kinetic energy a, directly computed kinetic energy dissipation rate (b), enstrophy based kinetic energy dissipation rate (c), and energy residual d with respect to $t / t_{c}^{}$ for grids of 64$^3$, 128$^3$, 256$^3$, and 512$^3$ cells, where $t_{c}^{}$ is the characteristic time scale. The reference data refers to the pseudo-spectral solution reported in van Rees et al. (2011)

It has been argued in DeBonis (2013) that the differences between $\epsilon (E_{k})$ and $\epsilon (\zeta )$ existed in coarse meshes (64$^3$ and 128$^3$ cells) might be due to both the numerical dissipation and dissipation due to the vorticity, while the former one reduces with improving the resolution. This observation is confirmed by looking into the contours of vorticity norm in this study as well, Fig. 16. It is clear that the vorticity structures are very well defined for domains with 256$^3$ and 512$^3$ cells, Fig. 16c, d, respectively. However, the vorticity structures are not well resolved for the domains with 64$^3$ and 128$^3$ cells. The current observation is in a very good agreement with vorticity contours presented in DeBonis (2013). This observation confirms the accuracy of using the cubic spatial discretization scheme implemented in OpenFOAM.

Fig. 16

Contours of vorticity norm for grids with 64$^3$, 128$^3$, 256$^3$, and 512$^3$ cells at $t / t_{c}^{} = 8$ (still transitioning), where $t_{c}^{}$ is the characteristic time scale. The contours are on $x = \pi$ plane, $y = \pi$ to $2\pi$, and $z = \pi$ to $2\pi$

Appendix B Homogeneous Isotropic Turbulence under Decaying Turbulence

The intention of this part is to validate isotropy and homogeneity of the flow using OpenFOAM. As previously mentioned in Sect. 2.1, the homogeneous isotropic velocity field is initialized using the Passot-Pouquet spectrum (Saad et al. 2017) on a tri-periodic cube with a side length equal to case B. The values of initial turbulence intensity and integral length scale are equal to the magnitudes reported for case B. The simulation is conducted under decaying turbulence, i.e. similar to the conditions examined in this study. Figure 17 shows the instantaneous velocity and vorticity fields for the aforementioned condition after running the simulation for two initial eddy turn over times.

Fig. 17

Instantaneous a velocity, and b vorticity fields resembling case B under the non-reacting test condition after running the simulation for two initial eddy turn over times

Figure 18 shows the cloud data of (u, v), (u, w), and (v, w) velocity component pair after running the simulation for two initial eddy turn over times. Results show that the data form a symmetric pattern indicating that there is not a directional preference in the flow field. Therefore, the average values of uv, uw, and vw are equal to zero. Satisfying this criteria ensures that the turbulence field is isotropic (Kundu et al. 2016).

Fig. 18

Scatter plot of (u, v) (left), (u, w) (middle), and (v, w) (right) velocity component pair after running the simulation for two initial eddy turn over times

Further analysis is required to validate the homogeneity of the flow. The velocity-derivative skewness $S_{u}$ and flatness $F_{u}$ factors are evaluated using the following expressions (Tao et al. 2020):

$$\begin{aligned} S_{u} = \frac{\left\langle \left[ \left( \partial u / \partial x \right) ^{3} + \left( \partial v / \partial y \right) ^{3} + \left( \partial w / \partial z \right) ^{3} \right] / 3 \right\rangle }{\left\langle \left[ \left( \partial u / \partial x \right) ^{2} + \left( \partial v / \partial y \right) ^{2} + \left( \partial w / \partial z \right) ^{2} \right] / 3 \right\rangle ^{3/2}}. \end{aligned}$$

(B4)

$$\begin{aligned} F_{u} = \frac{\left\langle \left[ \left( \partial u / \partial x \right) ^{4} + \left( \partial v / \partial y \right) ^{4} + \left( \partial w / \partial z \right) ^{4} \right] / 3 \right\rangle }{\left\langle \left[ \left( \partial u / \partial x \right) ^{2} + \left( \partial v / \partial y \right) ^{2} + \left( \partial w / \partial z \right) ^{2} \right] / 3 \right\rangle ^{2}}. \end{aligned}$$

(B5)

The time evolution of velocity derivative skewness and flatness factors are shown in Fig. 19. The skewness (flatness) factor reaches a value of approximately $-0.5$ (3.5) which is consistent with previous findings reported in Tao et al. (2020). This observation affirms that using the current discretization scheme offered in OpenFOAM and extensively utilized by Zirwes et al. (2020) in their recent work could be a promising approach in simulating homogeneous isotropic turbulence under decaying condition.

Fig. 19

Time evolution of velocity-derivative skewness and flatness factors under decaying turbulence