A Numerical Study on Premixed Turbulent Planar Ammonia/Air and Ammonia/Hydrogen/Air Flames: An Analysis on Flame Displacement Speed and Burning Velocity

The economic storage and transportation of ammonia (NH3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {NH}_3$$\end{document}), and its capability to be thermally decomposed to hydrogen (H2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {H}_2$$\end{document}) make it a potential carbon-free synthetic fuel for the future. To comprehend the fundamental characteristics of NH3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {NH}_3$$\end{document} as a primary fuel enriched with H2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {H}_2$$\end{document} under low turbulent premixed flame conditions, three quasi direct numerical simulations (quasi-DNS) with detailed chemistry and the mixture-averaged transport model are conducted under decaying turbulence herein. The Karlovitz number is fixed to 4.28 for all the test conditions. The blending ratio (α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document}), specifying the hydrogen concentration in the ammonia/hydrogen mixture, varies from 0.0 to 0.6. The results reveal that the mean value of the density-weighted flame displacement speed (Sd∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_{\textrm{d}}^{*}$$\end{document}) is similar to (higher than) the unstrained premixed laminar burning velocity (SL0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_{\textrm{L}}^{0}$$\end{document}) for NH3/\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {NH}_3/$$\end{document}air flame (NH3/H2/\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {NH}_3/\hbox {H}_2/$$\end{document}air flames). Furthermore, the performance of two extrapolation relations for estimating Sd∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_{\textrm{d}}^{*}$$\end{document} as linear and non-linear functions of flame front curvature is discussed thoroughly. The performances of both models are almost similar when evaluating the data near the leading edge of the flame. However, the non-linear one offers more accurate results near the trailing edge of the flame. The results show that the mean flame stretch factor increases with increasing the blending ratio, suggesting that the mean flamelet consumption velocity deviates from SL0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_{\textrm{L}}^{0}$$\end{document} by enriching the mixture with H2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {H}_2$$\end{document}. The mean value of the local equivalence ratio (ϕ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi$$\end{document}) for the turbulent NH3/\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {NH}_3/$$\end{document}air flame is almost equal to its laminar counterpart, while it deviates significantly for NH3/H2/\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {NH}_3/\hbox {H}_2/$$\end{document}air flames. In addition, the local equivalence ratio for the flame front with positive curvature values is higher than the negatively curved regions for NH3/H2/\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {NH}_3/\hbox {H}_2/$$\end{document}air flames due to H2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {H}_2$$\end{document} preferential diffusion. Furthermore, the results indicate that hydrogen is consumed faster in positively curved regions compared to the negatively curved zones due to enhanced reaction rates of specific chemical reactions.


Introduction
Decreasing greenhouse gas (GHG) emissions and increasing the efficiency of current combustion devices could be accompanied by utilizing carbon-free synthetic fuels.One of the promising carbon-free fuel components for reaching low emissions is hydrogen ( H 2 ) (Kob- ayashi et al. 2019).However, it should be emphasized that the economic storage and transportation of hydrogen are not resolved completely constituting a major obstacle in its sustainable utilization (Valera-Medina et al. 2018;Kobayashi et al. 2019).Considering this, other carbon-free fuel compounds should be examined as an alternative to pure H 2 .Due to its high hydrogen density, ammonia ( NH 3 ) has the potential of being utilized as a future carbon-free fuel since its transportation, storage, handling, distribution, and production are settled (Kobayashi et al. 2019).Since the boiling temperature and condensation pressure of pure ammonia are relatively similar to propane, ammonia transportation may not have significant challenges for transport ships (Kobayashi et al. 2019).However, it should be mentioned that the low burning velocity and heat of combustion of pure ammonia refrain its robust utilization as a main fuel.Therefore, improving its combustion properties with added hydrogen ( H 2 ) and/or enriched oxygen ( O 2 ) is a promising approach to the sustain- able utilization of NH 3 (Kobayashi et al. 2019).Despite the ongoing research in comprehending the fundamental mechanisms associated with premixed turbulent combustion, there still exist many unresolved problems regarding the properties of such flames.The flame displacement speed and the turbulent burning velocity are important fundamental properties of premixed turbulent flames, and the knowledge of these properties for novel fuels such as ammonia and its blends with either hydrogen or oxygen is essential for the design of modern combustion devices.
Examining the flame displacement speed ( S d ), as the local flame front speed relative to the local fluid flow velocity, helps in comprehending the local turbulence-chemistry interaction in-depth (Song et al. 2021).It should be noted that in premixed turbulent combustion, detailed knowledge about the flame displacement speed is required due to its wide application in the Flame Surface Density (FSD) modeling (Trouvé and Poinsot 1994).The flame displacement speed is considered as a local property, and its magnitude differs inside the flame brush (Giannakopoulos et al. 2015).Therefore, in order to compare the flame displacement speed values throughout the flame brush conclusively, the density-weighted flame displacement speed ( S * d ) which is defined by normalizing the flame displacement speed value with the ratio of the local density to the reactants density should be evaluated.For an unstrained premixed laminar planar flame, the density-weighted flame displacement speed across the flame brush is equal to the unstrained premixed laminar burning velocity ( S 0 L ) (Giannakopoulos et al. 2015).However, the aforementioned validity is disputed for the flames experiencing stretch which consists of the strain and curvature terms (Trouvé and Poinsot 1994).It should be mentioned that since premixed turbulent flames continuously experience different strain and curvature values, both the density-weighted and density-unweighted flame displacement speed values vary significantly across the flame brush (Trouvé and Poinsot 1994;Giannakopoulos et al. 2015).Gran et al. (1996) discussed the variation of S d and its components (reaction, normal diffusion, and tangential diffusion) with respect to the flame front curvature ( ) for a stoichiometric turbulent planar methane/ air flame based on 2-D direct numerical simulation (DNS) using detailed chemistry.They showed that the tangential diffusion component of the flame displacement speed contributes considerably to the flame displacement speed value for the positive and negative flame front curvature regions.Furthermore, their results indicated that S d is negatively correlated with the flame front curvature, and negative values for S d are observed in positively curved regions.Later, Echekki and Chen (1999) discussed the variation of S * d components with respect to the flame front curvature obtained from 2-D DNS results of a stoichiometric turbulent planar methane/air flame using detailed chemistry.They showed that the reaction (normal diffusion) component of S * d is always positive (negative) within the considered zone in their study, and the normal diffusion component of S * d does not vary with the flame front curvature.Furthermore, they showed that the flame displacement speed is negatively correlated with the flame front curvature.Chakraborty and Cant (2005) showed the negative correlation between S d and for premixed turbulent planar flames with different Lewis number values varying between 0.8 and 1.2 based on 3-D DNS using simple chemistry.They showed that the mean density-weighted flame displacement speed value divided by the unstrained premixed laminar burning velocity value is noticeably higher than unity throughout the flame brush when the Lewis number (Le) is equal to 0.8, and it does not deviate distinctly from the unity value when Le is equal to 1.0 and 1.2.Han and Huh (2008) showed a negative correlation between S d and for premixed turbulent planar flames with different turbulence intensity and Lewis number values based on 3-D DNS with simple chemistry.A similar negative correlation between S * d and was observed recently for a wide range of turbulence intensity values for premixed turbulent planar flames by Nivarti and Cant (2019) based on 3-D DNS using simple chemistry.They showed that the mean value of S * d is approximately equal to S 0 L for the test cases considered in their study by assuming the unity Lewis number for the species.Song et al. (2021) showed that the most probable value of S * d is similar to S 0 L for highly turbulent premixed hydrogen/air flames with an equivalence ratio of 0.7 based on 3-D DNS using detailed chemistry.More recently, Herbert et al. (2020) and Chakraborty et al. (2022) compared the performance of various extrapolation relations connecting S * d to the flame front curvature and stretch rate based on 3-D DNS with simple/detailed chemistry for thermo-diffusively neutral flame conditions under various turbulence intensity values.It should be noted that the extrapolation relations utilized in their studies were originally developed for weakly stretched premixed laminar flames (Chakraborty et al. 2022).They showed that the performance of the linear model based on the flame front curvature ( ) is exceptionally robust for their test cases compared to other relations such as LS model (linear model based on stretch) and NQ model (quasisteady non-linear model) (Herbert et al. 2020).
In addition to the flame displacement speed, the turbulent burning velocity ( S T ), as a global property, is considered as a fundamental property for premixed turbulent combustion (Peters 2000;Poinsot and Veynante 2012).Many attempts have been performed in the past to propose correlations for predicting the turbulent burning velocity, see e.g.Clavin and Williams (1979); Gülder (1991), and the detailed review paper by Lipatnikov and Chomiak (2002).The turbulent burning velocity results obtained from the proposed correlations in the literature seem to be highly sensitive to the flow configuration (Bilger et al. 2005).Therefore, the applications of these correlations seem to be restricted.It should be noted that the turbulent burning velocity and the flame surface area are connected through the well-known Damköhler's first hypothesis within the flamelet regime (Damköhler 1940;Driscoll 2008).According to this hypothesis, the ratio of the turbulent burning velocity to the unstrained premixed laminar burning velocity (S T ∕S 0 L ) is equal to the ratio of the wrinkled flame surface area to the unwrinkled flame surface area (A T ∕A L ) .Later, Bray and Cant (1991) modified Damköhler's first hypothesis by proposing a flame stretch factor ( I 0 ) defined as the ratio of the mean flamelet consumption velocity to the unstrained premixed laminar burning velocity (S � L ∕S 0 L ) .It should be noted that Damköhler's first hypothesis was validated, i.e. S � L = S 0 L , for statically premixed turbulent planar flames based on DNS, see e.g.Chakraborty and Cant (2005); Han and Huh (2008); Nivarti andCant (2017a, 2017b); Varma et al. (2021), andBrearley et al. (2020).According to the results presented in the aforementioned papers, this hypothesis seems to be valid when the Lewis number of the flame is set to the unity value suggesting that S � L = S 0 L .Very recently, Song et al. (2021) validated this hypothesis for highly turbulent lean H 2 ∕air premixed flames by showing that S T ∕S 0 L and A T ∕A L are highly correlated, and the mean flame stretch factor values are near the unity value.However, the validity of this hypothesis might be disputed when the fuel Lewis number deviates from the unity value.For example, Haworth and Poinsot (1992) showed that the ratio of the turbulent burning velocity to the unstrained premixed laminar burning velocity is higher (lower) than the ratio of the wrinkled flame surface area to the unwrinkled flame surface area when the reactants Lewis number is lower (higher) than the unity value.A similar observation was previously reported in the literature, see e.g.Trouvé and Poinsot (1994); Chakraborty and Cant (2005), and Han and Huh (2008).
It should be emphasized that the knowledge of the local flamelet consumption velocity as the flame speed where the reactants are consumed (Poinsot and Veynante 2012), and the fuel consumption rate help to improve our understanding of the behavior of premixed turbulent flames.Haworth and Poinsot (1992) and Rutland and Trouvé (1993) discussed the effect of the Lewis number on the local flamelet consumption velocity using DNS.They showed that S c is larger (smaller) than S 0 L on the flame surface with positively (negatively) curved regions when the Lewis number value is less than unity.The opposite trend is observed when the Lewis number value is higher than unity.Following the works presented in Haworth and Poinsot (1992) and Rutland and Trouvé (1993), Bell et al. (2007) showed that the local burning velocity is positively (negatively) correlated with the flame front curvature for hydrogen/air (propane/air and methane/air) flames based on 2-D DNS using detailed chemistry.They mentioned that higher local flamelet burning velocity values in negatively curved regions for the propane/air flames indicate that the fuel consumption rate is intensified in these zones.More recently, Lee et al. (2021) showed that the fuel consumption rate is enhanced in positively curved regions for premixed turbulent lean H 2 ∕air flames in the thin reaction zones regime propagating in forced turbulence based on 3-D DNS with detailed chemistry.In their recent 3-D DNS study with detailed chemistry, Rieth et al. (2021) showed the augmentation of hydrogen consumption rate in positively curved regions for preheated NH 3 ∕H 2 ∕N 2 ∕air flames developing in a shear layer within the broken reaction zone regime.They showed that chemical reactions such as H 2 + OH ⇔ H + H 2 O consuming H 2 (in the forward reaction) are enhanced in positively curved regions due to the local enhancement of the equivalence ratio.Furthermore, they noted that this augmentation is promoted under pressurized conditions.
To the best of the authors' knowledge, the flame displacement speed statistics and the validity of Damköhler's first hypothesis for premixed turbulent NH 3 ∕air and NH 3 ∕H 2 ∕ air flames are not explored thoroughly in the literature.In addition, the effect of H 2 enrichment on the local equivalence ratio, molecular diffusion flux, fuel consumption rate, and reactions consuming/producing molecular hydrogen across the flame brush needs to be discussed systematically.Understanding these details helps to explore the underlying physics behind the local/global characteristics of premixed turbulent NH 3 ∕ air and NH 3 ∕H 2 ∕air flames.Such an understanding gives an insight into the accurate combustion modeling of such flames.Hence, the present study focuses on examining these properties based on quasi direct numerical simulations (quasi-DNS) using detailed chemistry and the mixture-averaged transport model to consider the preferential diffusion.In this respect, the main objectives of the present work are (1) to study the density-weighted flame displacement speed statistics and the accuracy of extrapolation relations evaluating the density-weighted flame displacement speed, (2) to assess the validity of Damköhler's first hypothesis, (3) to examine the variation of the local equivalence ratio across the flame brush and investigate its change with the flame front curvature, and (4) to explore the influence of the positively/negatively curved regions on the molecular diffusion flux of H 2 ∕H , fuel consumption rate, and important reac- tions consuming/producing H 2 .The remainder of the paper is organized as follows.The numerical methods will be discussed thoroughly in Sect. 2. The results will be discussed in Sect.3. Finally, the summary of the main concluding remarks will be outlined in Sect. 4.

Flame configuration
Numerical simulations of premixed turbulent planar NH 3 ∕air and NH 3 ∕H 2 ∕air flames have been performed under the equivalence ratio of 1.0.It should be mentioned that for the combustion of ammonia, the NO x level is reduced for stoichiometry and slightly rich mixture (Kobayashi et al. 2019).This justifies the selection of the equivalence ratio of unity for the current study.The blending ratio ( ) indicating the hydrogen concentration in the ammonia/hydrogen mixture is evaluated as X H 2 ∕(X H 2 + X NH 3 ) , where X i is the mole fraction of species i.This ratio is systematically varied from 0.0 to 0.6 in this study to promote the combustion properties of NH 3 ∕air flames.The reactants gas temperature ( T r ) is set to 298 K under standard atmospheric pressure.This results in an unstrained premixed laminar burning velocity S 0 L (laminar thermal flame thickness th ) of approximately 7.01 (2.007), 30.37 (0.604), and 68.30 (0.376) cm/s (mm) for the conditions with = 0.0, 0.4, and 0.6, respectively.It should be noted that th is evalu- ated using the maximum temperature gradient.The adiabatic flame temperature ( T ad ) for the conditions with = 0.0, 0.4, and 0.6 are then equal to 2066, 2146, and 2179 K, respectively.
The Passot-Pouquet spectrum (Passot and Pouquet 1987) is utilized to initiate an isotropic homogeneous turbulent flow field in a rectangular cube with a domain size equal to 80 th × 10 th × 10 th .It should be noted that this spectrum is implemented in the latest version of a Python script which is originally developed by Saad et al. (2017).The domain has been discretized on a Cartesian grid of 1760 × 220 × 220 cells with uniform spacing.Flame-turbulence interactions occur under decaying turbulence.The prescribed turbulent flow is then superimposed on the solution of the 1-D unstrained premixed reference laminar flame.The grid spacing verifies that 22 computational cells reside across the flame thickness.In addition, the initial Kolmogorov length scale ≈ ΛRe

−3∕4 Λ
resides within at least one computation cell for all the test cases, where Λ and Re Λ are the integral length scale and turbulent Reynolds number, respectively.The initial non-dimensional turbulence intensity ( u � ∕S 0 L ), non-dimensional turbulent integral length scale ( Λ∕ th ), turbulent Reynolds number ( Re Λ = u � Λ∕ ), Karlovitz num- ber ( Ka = (u � ∕S 0 L ) 1.5 (Λ∕ th ) −0.5 ), and Damköhler number ( Da = ΛS 0 L ∕u � th ), are shown in Table 1, where u ′ and are the turbulence intensity and reactants kinematic viscosity, respectively.The reactants effective Lewis number ( Le eff ) is evaluated using the follow- ing expression (Ichikawa et al. 2015): where Le i corresponds to the Lewis number of species i.The effective Lewis number value, presented in Table 1, decreases from 0.97 to 0.69 with increasing the blending ratio from 0.0 to 0.6.It should be mentioned that all the test conditions lie within the thin reaction zones combustion regime since the Karlovitz number is higher than the unity value.According to Table 1, is increased from 0.0 to 0.6, while u � ∕S 0 L and Λ∕ th are kept con- stant.It should be borne in mind that AH000, AH040, and AH060 in Table 1 correspond to the test cases where is equal to 0.0, 0.4, and 0.6, respectively.

Reactive flow solver
The mass, momentum, species mass fractions, and energy conservation equations are solved with the compressible PIMPLE algorithm using the open-source C++ CFD toolbox OpenFOAM-v8 (Weller et al. 1998; The OpenFOAM Foundation 2020) to perform quasi direct numerical simulations (quasi-DNS) with detailed chemistry for premixed turbulent planar NH 3 ∕air and NH 3 ∕H 2 ∕air flames.Recently, Zirwes et al. (2020) thoroughly dis- cussed and reviewed the capability of OpenFOAM in performing DNS, and the terminology of quasi-DNS is adopted from their work.It should be noted that performing the 3-D quasi-DNS of chemically reacting flows with detailed chemistry and the mixture-averaged transport model is computationally very expensive.Hence, the test cases explored in the current study are limited to the three cases mentioned earlier in Table 1.
The chemical kinetic mechanism of Stagni et al. (2020) including 31 species and 203 reactions is utilized in this study.The performance of this mechanism over a wide range of experimental conditions against different chemical kinetic mechanisms is discussed by Shrestha et al. (2021).For example, in their measurements, Lhuillier et al. (2020) reported that the laminar burning velocity values of stoichiometric NH 3 ∕H 2 ∕air flames with blend- ing ratios of 0.0, 0.4, 0.6 were equal to 6.52, 28.17, and 66.38 (cm/s) when the reactants temperature and pressure are equal to 298 K and 1 atm, respectively.Very recently, Gotama et al. (2022) reported a laminar burning velocity value of 29.03 (cm/s) for the same mixture with a blending ratio of 0.4 under the atmospheric condition.Comparing (1) ,  2020) mechanism show that choosing this chemical kinetic mechanism is a reasonable choice for the conditions tested in this work.However, Gotama et al. (2022) showed that the performance of some of the chemical kinetic mechanisms developed in the literature deteriorates at certain conditions for NH 3 ∕H 2 ∕air flames.This suggests that the uncertainty of predicting the laminar burning velocity could be large using different chemical kinetic mechanisms developed in the literature.
The reactingDNS solver wherein the mixture-averaged transport model is utilized for describing the reaction-diffusion process precisely is used in this study.It should be emphasized that the successful implementation of this solver was previously shown in Zhong et al. 2018Zhong et al. , 2020a, b;, b;Zhang et al. 2020, andRocha et al. 2021.Utilizing this transport model leads to more accurate results compared to the unity/non-unity Lewis number simulations when simulating the mixtures containing high mass diffusivity species such as the molecular/atomic hydrogen.It should be mentioned that the viscosity and thermal conductivity for each species as well as the binary diffusion coefficient for each species pair are evaluated using the third-order logarithmic polynomial fitting equations.The Wilke formula (Wilke 1950) is used for estimating the viscosity of the mixture, and a combination averaging formula is utilized for estimating the mixture-averaged thermal conductivity (Mathur et al. 1967).The mixture-averaged diffusion coefficient for each species is evaluated from the methodology described by Bird et al. (1960).The thermodynamic properties for each species are evaluated from the polynomials fitted to NIST-JANAF thermochemical tables.The obtained data are then utilized to calculate the thermodynamic properties of the mixture.It should be mentioned that the thermal diffusion effect (Soret effect) is ignored in the current study.Aspden et al. (2011a) mentioned that the Soret effect is capable of shifting the quantitative characteristics of lean H 2 ∕air flames slightly.However, the qualita- tive trends remain intact.This effect is also neglected in some other studies as well investigating premixed turbulent CH 4 ∕air and H 2 ∕air flames, see e.g.Aspden et al. (2016) and Tanahashi et al. (2000).Very recently, Netzer et al. (2021) studied the thermal diffusion effect on wrinkled laminar NH 3 ∕H 2 ∕N 2 ∕air premixed flames.They showed that the correlations between the NO mass fraction and the local temperature, the OH mass fraction, and the local equivalence ratio are similar by either including or excluding the thermal diffusion effect from the simulations.They concluded that the diffusion of hydrogen has a greater impact on the NO mass fraction level than the thermal diffusion.
The open-source library of pyJac (Niemeyer et al. 2017) is then linked to the react-ingDNS solver.This library provides C subroutines for generating the analytical Jacobian of a prescribed chemical kinetic mechanism.The analytical Jacobian is required for solving a system of ordinary differential equations (ODE) which is used for calculating the chemical source terms, i.e. the heat release rate and the net consumption/production rate of each species in the energy and species transport equations, respectively.In this work, every computational cell is treated as a stand-alone chemistry problem where their state is defined through pressure, temperature, and species mass fractions, i.e. local thermochemical composition vector.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.Hence, an operator-splitting approach is utilized 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 estimated by solving the system of ordinary differential equations using the semi-implicit Euler method (Seulex) through integrating in each chemistry problem over the CFD time step (Hairer and Wanner 1996).
It should be noted that an open-source dynamic load balancing model (DLBFoam), developed recently by Tekgül et al. (2021), is used to overcome the computational load imbalance among the processors by distributing the computation of chemistry problems equally among all the processors using the MPI communication protocol.This model was previously used for investigating the spray-assisted dual-fuel ignition (Kannan et al. 2021;Karimkashi et al. 2022), spray-assisted tri-fuel ignition (Gadalla et al. 2021), and premixed turbulent methane/air and ammonia/air flames (Tamadonfar et al. 2022).
Inflow and outflow boundary conditions with a zero mean inlet velocity in the direction of the mean flame propagation are used.In the lateral direction, periodic boundary conditions are employed, showing an unbounded planar flame brush.For pressure, a nonreflecting boundary condition is applied for both inflow and outflow regions.A secondorder implicit backward Euler method and cubic scheme have been used for time and spatial discretizations, respectively.An interested reader is referred to Zirwes et al. (2020) for detailed information on the cubic discretization scheme implemented in OpenFOAM.Each simulation runs for two initial eddy turn-over times ( t e ) verifying that the simula- tion time ( t sim ) is equal/larger than the chemical time scale ( t chem ), where t e = Λ∕u � and t chem = th ∕S 0 L .It should be emphasized that the simulation time for premixed turbulent flames propagating under decaying turbulence is generally considered within 2-4 eddy turn-over times (Klein et al. 2017).Therefore, all the statistics are evaluated at two initial eddy turn-over times in this work.It should be emphasized that our investigations indicate that the observed qualitative trends herein are not sensitive to the choice of the simulation time within 2-3 eddy turn-over times.

Density-weighted flame displacement speed and its components
The instantaneous flame displacement speed ( S d ) is decomposed into three components using the following relation (Gran et al. 1996;Echekki and Chen 1999): where S d, r = ω∕|∇c| , S d, n = n.∇(Dn.∇c)∕ |∇c| , and S d, t = −2D are the reaction, nor- mal diffusion, and tangential diffusion components of S d , respectively.In Eq. 2, ω is the reaction rate of the reaction progress variable, is the density, c is the progress variable based on the temperature T evaluated as (T − T r )∕(T ad − T r ) , n = −∇c∕|∇c| is the normal vector on a specified progress variable pointing towards the reactants, D is the progress variable diffusivity, and is the flame front curvature estimated as 0.5 (∇.n) (Trouvé and Poinsot 1994;Han and Huh 2008).It should be noted that the positive (negative) flame front curvature is convex (concave) towards the reactants.It should be mentioned that when the reaction progress variable is defined based on the temperature value, an additional term due to the diffusion fluxes through different species appears in Eq. 2 which is one to two orders of magnitude smaller than the other terms (Giannakopoulos et al. 2015;Yuvraj et al. 2022).Therefore, we have not included this term for the sake of simplicity in this analysis.Furthermore, the velocity acceleration resulting from the thermal expansion across the flame is eliminated by using the density-weighted form of the flame displacement speed ( S * d = S d ∕ 0 ), where 0 is the reactants density (Song et al. 2021).Therefore, comparing (2) the density-weighted flame displacement speed to the unstrained premixed laminar burning velocity is meaningful (Giannakopoulos et al. 2015).
Figure 1 shows the mean values of density-weighted flame displacement speed and its components, i.e. S * d, r = S d, r ∕ 0 , S * d, n = S d, n ∕ 0 , and S * d, t = S d, t ∕ 0 , with respect to the progress variable for all the test conditions.Each component has been normalized by the unstrained premixed laminar burning velocity ( S 0 L ) of its reference condition, and the statistics are captured at two initial eddy turn-over times.The results show that the region where S * d, r is not equal to zero increases with increasing considering that the intense heat release happens over a broader region across the flame brush.Furthermore, S * d, n is positive (negative) in the region where c is approximately lower (higher) than 0.6.It is worth mentioning that with increasing , the peak values of S * d, r ∕S 0 L and S * d, n ∕S 0 L shift towards the leading edge of the flame since the location of maximum heat release rate moves towards the leading edge of the flame (not shown), and their values decrease significantly.Since the mean value of the flame front curvature is near zero for each progress variable, the ensemble average value of S * d, t is approximately equal to zero with respect to the progress variable for all the test conditions.This implies that the contribution of the mean value of S * d, t on the density-weighted flame displacement speed is insignificant.This observation was previously reported in the literature.For example, Chakraborty and Cant (2005) showed that the mean value of S d, t is near zero across the flame brush for premixed turbulent planar flames with different Lewis number values based on 3-D DNS using simple chemistry.
The results in Fig. 1 indicate that S * d deviates from S 0 L over a broad range of c with increasing the blending ratio.Furthermore, the density-weighted flame displacement speed ( S * d ) averaged over the entire flame brush is approximately equal to S 0 L when = 0.0 , and it is approximately 12% (22% ) higher than S 0 L when = 0.4 (0.6).This observation suggests that the local flame elements for pure NH 3 ∕air flame, on average, propagate identical to the unstrained premixed laminar flame, and they move faster than the reference laminar condition with increasing , i.e. adding H 2 to the NH 3 ∕air flame.In addition, Fig. 1 shows that the flamelets at the trailing edge of the flame brush propagate faster than its leading edge with enriching ammonia, as the main fuel, with hydrogen, i.e. decreasing the effective Lewis number value.It should be mentioned that the dependency of the flame displacement speed variation across the flame brush with the Lewis number has already been discussed thoroughly by Trouvé and Poinsot (1994).They showed that when the Lewis number is extremely low, i.e. equal to 0.3, the flame elements propagate much faster than the laminar flame counterpart.In contrast to the observation reported in Fig. 1, Chakraborty and Cant (2005) showed that the mean value of S * d at the leading edge of the flame brush is larger than its mean value at its trailing edge for premixed turbulent flame simulations under different Lewis number values.Furthermore, it is required to study the effect of the flame front curvature on the density-weighted flame displacement speed in the following section since its value varies in positively/negatively curved regions.

Effect of the flame front curvature on the density-weighted flame displacement speed
Demonstrating the density-weighted flame displacement speed ( S * d ) in terms of the local flame front curvature ( ) is required in the context of Flame Surface Density (FSD) modeling (Herbert et al. 2020;Chakraborty et al. 2022).Figure 2 2022) for premixed turbulent planar methane/air and hydrogen/air flames, respectively.Furthermore, increasing from 0.0 to 0.6 results in a higher data scatter which is more pronounced near the trailing edge of the flame, i.e. d values is insignificant.This observation suggests that the turbulent structures are not very strong due to the absence of negative S * d values for these flame conditions.This scenario is in line with low turbulence intensity results of lean premixed H 2 ∕air flames presented by Song et al. (2021).However, they showed that the probability of observing negative S * d values increases significantly under extremely high turbulence intensity test conditions near the leading edge of the flame.
Following the methodology presented in Herbert et al. (2020) and Chakraborty et al. (2022), the performances of two extrapolation relations for estimating the density-weighted flame displacement speed ( S * d ), namely LC and N3P models, are discussed for premixed turbulent NH 3 ∕air and NH 3 ∕H 2 ∕air flames.The former, developed by Markstein (1951), is a linear model based on the flame front curvature.In this model, the ratio of the densityweighted flame displacement speed to the laminar burning velocity (S * d ∕S 0 L ) is evaluated using 1 − L M , where L M is the Markstein length.The latter, recently developed by Wu et al. (2015), is a non-linear extrapolation model with three terms, and th , where C is a free-parameter.Using linear and non-linear regression analysis, extrapolation results for two different progress variable values, i.e. c = 0.2 and 0.8, are presented in Fig. 2. Table 2 summarizes the values of L M ∕ th and C estimated using S * d and from our quasi-DNS data, correlation coefficient r between S * d acquired from our quasi-DNS and predicted data from extrapolation relations, and the goodness of fit g.o.f. for all the test conditions.An interested reader is referred to Herbert et al. (2020) and Chakraborty et al. (2022) for detailed explanation on r and g.o.f evaluations.The results indicate that the correlation coefficient obtained from the N3P model is slightly higher than the corresponding value evaluated from the LC model, and its value decreases with increasing c for each model.A further observation reveals high g.o.f.when the data are conditioned at c = 0.2 , and the accuracy of these models decreases significantly at higher progress variable values, i.e. c = 0.8 .It should be noted that the LC model is incapable of capturing the non-linearity relationship between and S * d leading to its weakened performance, especially at a higher progress variable value, i.e. c = 0.8 .Considering the data presented in Table 2, the performance of the N3P model is better than the LC model in general.However, the enhanced performance of this model is due to the additional freeparameter C making this model more complex for modeling purposes.A similar observation was thoroughly discussed very recently in Herbert et al. (2020) and Chakraborty et al. (2022) for thermo-diffusively neutral premixed turbulent planar flames based on 3-D DNS with simple/detailed chemistry.

Assessing the validity of Damköhler's first hypothesis
Examining the validity of Damköhler's first hypothesis for novel fuel blends is required due to its importance in premixed turbulent combustion modeling.Following Chakraborty et al. 2019, the ratio of the turbulent burning velocity ( S T ) to the so-called modified burn- ing velocity ( S ′ L ) is defined as follows: where ∑ gen = �∇c� is the flame surface density, A T = ∫ V |∇c|dV is the wrinkled flame sur- face area, and A L is the cross-sectional area normal to the direction of mean flame propa- gation, and ( S d ) s = ( S d )|∇c|∕|∇c| is the surface-averaged value of S d .Evaluating the wrinkled flame surface area ( A T ) based on the integration of the flame surface density over the flame brush has a number of advantages over calculating the flame surface area using a specific iso-surface as discussed in Klein et al. (2020).Among these advantages, the independency of the flame surface area value using the former method from any threshold values is desirable, and this method is utilized for evaluating the wrinkled flame surface area ( A T ) in the current study.It has been shown by Boger et al. (1998) that ( S d ) s is approxi- mately equal to 0 S 0 L by assuming unity Lewis number for all species using simple chemistry.This results in reducing Eq. 3 to Damköhler's first hypothesis relation as follows: Table 3 summarizes S T ∕S 0 L , A T ∕A L , and I 0 for all the test cases in this work, where the mean flame stretch factor I 0 is defined as the ratio of S T ∕S 0 L to A T ∕A L .The results show that both S T ∕S 0 L and the ratio between S T ∕S 0 L and A T ∕A L increase with increasing the blend- ing ratio.Our results suggest that Damköhler's first hypothesis remains relatively valid for case AH000, where the ratio of S T ∕S 0 L to A T ∕A L is near unity.It should be mentioned that the validity of this hypothesis was confirmed for the unity Lewis number flames based on DNS with simple/detailed chemistry (Chakraborty and Cant 2005;Han and Huh 2008;Nivarti and Cant 2017a, b;Brearley et al. 2020;Varma et al. 2021, andTamadonfar et al. 2022).Furthermore, this hypothesis might become less valid with increasing to 0.4 and 0.6, i.e. cases AH040 and AH060.It should be noted that with increasing , the effective Lewis number ( Le eff ) decreases from 0.97 to 0.69.As a result, this observation suggests that increasing (decreasing Le eff ) results in a faster diffusion of reactants into the reaction zone, while the rate at which the heat is conducted out from this zone is slower.Therefore, the enhancement of turbulent burning velocity is higher than the creation of the wrinkled flame surface area for these test conditions.This observation is in line with previous DNS results reported in the literature, see e.g.Haworth and Poinsot (1992); Rutland and Trouvé ( 1993); Trouvé andPoinsot (1994), andChakraborty et al. (2014).The results presented in Table 3 suggest that the mean flamelet consumption velocity deviates from the laminar burning velocity with increasing the blending ratio, implying the deviation from Damköhler's first hypothesis with adding H 2 to pure NH 3 ∕air flames.This means that the flame elements, on average, burn faster than their laminar counterpart by enriching the mixture with H 2 .It should be highlighted that the turbulent burning velocity values extracted from simulations of turbulence in a box are lower than the corresponding values obtained from experiments under the same turbulence intensity, as mentioned by Steinberg et al. (2021).
In addition, comparing the turbulent burning velocity results obtained from the current quasi-DNS data and the experimental measurements might not be straightforward due to the existence of several constraints.One of the reasons, described in the review paper by Driscoll et al. (2020), is the presence of periodic boundary conditions in the numerical simulations where these boundaries do not exist in experiments.Since the boundary conditions and the wrinkling process for the current numerical setups are not identical to any experimental measurement systems, direct comparison between their results might be inconclusive.Another reason why it is challenging to compare the current quasi-DNS results with experimental data is the use of decaying turbulence, resulting in low simulation time compared to the eddy turn-over time.

Local equivalence ratio across the flame brush
It should be mentioned that the local flamelet velocity is dependent on the local equivalence ratio ( ) (Netzer et al. 2021), and investigating the distribution of local across the flame brush can deepen our understanding of flame propagation for the test conditions examined in this study.Following the methodology presented in Chen and Bilger (2004); Day et al. (2009); Netzer et al. (2021), and Lee et al. (2021), the local equivalence ratio ( ) is evaluated using the following expression: where Z F and Z O are evaluated using 3 shows examples of instantaneous 2-D slices of local equivalence ratio ( ) for all the cases studied in this work.It is observed that the local changes significantly (insignificantly) across the flame brush for cases AH040 and AH060 (case AH000), and its value is strongly dependent on the flame front curvature.The underlying mechanism behind the aforementioned trend will be addressed thoroughly in the next sub-section of this manuscript.
In order to have an in-depth comprehension of the local equivalence ratio distribution across the flame brush, the scatter plot of the local equivalence ratio ( ) for each turbulent case as well as the mean values for both the laminar (dashed black line) and turbulent (solid black line) cases are presented in Fig. 4, resulting in the following observations.First, the local equivalence ratio decreases from the global equivalence ratio of unity within the flame brush for the laminar flames.This deviation from the unity value is less noticeable for case AH000 due to the fact that the effective Lewis number for this particular test case is near the unity value.However, there still exists the effect of hydrogen diffusion in this system which is capable of changing Z F ∕Z O .A similar observation was previously reported by Aspden et al. (2011a) for lean CH 4 ∕air flame where the fuel Lewis number is near the unity value.Second, this deviation enhances by adding H 2 to NH 3 ∕ air flame due to the diffusion of the fuel component, namely H 2 , into the flame.Third, scatter variation of the local equivalence ratio ( ) from a reference mean value increases with increasing .Fourth, increasing the blending ratio ( ) results in growing the difference between the mean values of local for turbulent and laminar flames across the flame brush.For example, the mean value of ∕ L conditioned at c = 0.75 increases from 1.016 to 1.064 when increasing the blending ratio from 0.0 to 0.6, where L corresponds to the  et al. (2021).Finally, the local returns to its global value, i.e. = 1 in this study, near the product regions where the effect of chemical reactions is insignificant.

Effect of the flame front curvature on the local equivalence ratio
Since the local equivalence ratio tends to vary under curved flame components, examining the effect of the flame front curvature ( ) on the local equivalence ratio ( ) is required.Considering this, the joint probability density functions (PDFs) between the normalized local equivalence ratio ( ∕ L ) and the normalized flame front curvature ( th ) are shown in Fig. 5 by sampling the data over c = 0.75 iso-surface.The results show a strong (mild) positive correlation between ∕ L and th for cases AH040 and AH060 (case AH000).It should be noted that selecting different progress variable values for sampling the data do not change the trend observed in Fig. 5.This observation indicates that the local equivalence ratio ( ) for flame elements with positive curvature values is higher than the negatively curved regions (see Fig. 3).This could be due to the enhancement of preferential diffusion of H 2 from the reactants into the positively curved reaction zone regions with increasing the blending ratio.A similar correlation between and was reported very recently based on 3-D DNS with detailed chemistry in the thin reaction zones regime for premixed turbulent H 2 ∕air flames propagating in forced turbulence (Lee et al. 2021), and in the broken reaction zone regime for preheated NH 3 ∕H 2 ∕N 2 ∕air flames developing in a shear layer (Rieth et al. 2021).
To reveal the underlying mechanism behind the observation reported in Fig. 5, it is required to address the diffusive fluxes in detail.Following Gran et al. (1996); Echekki and Chen (1999); Poinsot and Veynante (2012), and Lee et al. (2021), the total molecular diffusion flux MD T is evaluated as follows: where D m is the diffusivity of species m evaluated by where Y m is the mass fraction of species m, D mj is the binary diffusion coefficient between two species m and j, and MW is the molecular weight of the mixture.For simplicity, D m corresponds to the value obtained from the reference laminar flame.It should be emphasized that using D m from either the simulation results or the reference laminar flame does not change the results discussed herein.In Eq. 6, n m is the unit vector normal to an iso-sur- face evaluated as It should be noted that the first (second) term on the R.H.S. of Eq. 6 is the normal diffusion (curvature) term, and the last term corresponds to the diffusion flux due to the gradient of the molecular weight of the mixture.
It should be emphasized that the molecular diffusion fluxes of H 2 and H are examined due to their high diffusivity compared to the oxidizer.Figure 6 shows the mean values of the normalized molecular diffusion flux terms of H 2 and H versus the normalized flame front curvature ( th ) sampled at c = 0.75 iso-surface for all the test conditions.The flux terms have been normalized with the total diffusion flux value (L.H.S. of Eq. 6) evaluated for the corresponding laminar flame.The joint PDFs between the normalized total diffusion flux and the normalized flame front curvature are shown along the mean values of each diffusion flux term introduced earlier in Fig. 6 to present an improved perception of the diffusion flux distribution.The normal diffusion flux term, curvature diffusion flux term, diffusion flux due to the gradient of the molecular weight of the mixture, and total diffusion terms, shown in Fig. 6, are represented by MD 1 , MD 2 , MD 3 , and MD T , respectively.
The results indicate that the normal diffusion terms for both H 2 and H are not cor- related with the variation of the flame front curvature, and the diffusion flux values due to the gradient of the molecular weight of the mixture are quite negligible.However, the curvature and total molecular diffusion flux terms are positively (negatively) correlated with the flame front curvature for H 2 (H), except that the molecular hydrogen diffusion flux terms do not show any clear correlation with the flame front curvature for pure ammonia/air flame, i.e. case AH000.This observation indicates that the molecular hydrogen diffuses actively towards the flame front with positive curvature values for NH 3 ∕air flames enriched with H 2 , i.e. cases AH040 and AH060.Furthermore, the defocusing of molecular hydrogen from the negatively curved zones is clearly observed for case AH040.It should be noted that the defocusing of molecular hydrogen from the negatively curved regions for case AH060 is noticeable if the data are sampled at lower progress variable values (not shown for the sake of brevity).However, it should be mentioned that sampling the data at different progress variable values does not change the positive correlation dependency between the molecular diffusion flux terms and the flame front curvature for ammonia/hydrogen/air flames.The local enrichment of molecular hydrogen due to the preferential diffusion, resulting in higher H 2 concentra- tion in positively curved regions, describes the positive correlation between the local equivalence ratio and the flame front curvature.Furthermore, the diffusion of H towards the unburned zone is amplified in negatively curved zones, and the defocusing of H could be observed when the flame front curvature is positive.A similar observation was reported very recently by Lee et al. (2021) for lean premixed turbulent H 2 ∕air flames (located within the thin reaction zones regime), and by Rieth et al. (2021) for preheated H 2 ∕air flames and NH 3 ∕H 2 ∕N 2 ∕air flames (located within the thin and broken reaction zones regimes).
It should be mentioned that the variation of the total diffusion flux for the molecular hydrogen under curved regions has an effect on the fuel consumption rate.Figure 7 shows the joint PDFs between the hydrogen/ammonia production rate and the flame front curvature sampled at c = 0.75 iso-surface for case AH040.The fuel pro- duction rate has been normalized with the corresponding absolute value of the reference laminar flame.It should be borne in mind that the negative (positive) production rate value corresponds to the consumption (production) of the species.The results show that the molecular hydrogen H 2 is consumed at a faster rate in positively curved Fig. 7 a Joint PDF between the normalized hydrogen production rate and the normalized flame front curvature ( th ) for case AH040.b Joint PDF between the normalized ammonia production rate and the normalized flame front curvature ( th ) for case AH040.The data are sampled at c = 0.75 iso-surface.The fuel production rate for each fuel component has been normalized by the corresponding absolute value of the reference laminar flame obtained at c = 0.75 iso-surface.The negative value for the production rate implies that the species is being consumed regions compared to the zones with the negative flame front curvature values, Fig. 7a, because the total diffusion flux value is higher under the former regions, see Fig. 6c.This enhanced fuel component consumption rate in positively curved regions could be considered an indication of a higher local burning velocity value at the aforementioned zones.A similar observation of the enhanced fuel consumption rate for the high mass diffusivity species in positively curved regions was previously reported in the literature, see e.g. Lee et al. (2021) and Rieth et al. (2021).Moreover, Bell et al. (2007) reported that the local burning velocity for a premixed lean H 2 ∕air flame is higher in positively curved regions than in negatively curved zones, indicating the fuel consumption rate is enhanced in positively curved zones.It should be noted that the negative correlation between the molecular hydrogen production rate and the flame front curvature is not observed anymore for the current test case when sampling the data at a higher progress variable value, i.e. c > 0.85 .The current observation regarding the correlation between the H 2 production rate and the flame front curvature is also valid when the blending ratio value is equal to 0.6, i.e. case AH060.Furthermore, the results indicate that the consumption rate of NH 3 is not dependent on the flame front curvature for the current test case, i.e. Figure 7b.This observation suggests that the fuel consumption rate is not dependent on the flame front curvature when the species Lewis number is relatively near the unity value, i.e.NH 3 species.It should be emphasized that the observed trend between the NH 3 production rate and the flame front curvature is the same for other test cases studied in this work (not shown for the sake of brevity).This observation was previously shown by Aspden et al. (2016) for premixed turbulent lean CH 4 ∕air flames with the global Lewis number of around unity, where they reported the fuel consumption rate is uniform along the flame front for low Karlovitz number flames.
Furthermore, detecting the reactions responsible for the trend observed in Fig. 7a is of particular importance.Among 203 reactions existing in the utilized chemical kinetic mechanism, 23 reactions include the molecular hydrogen species.By investigating all these reactions, two reactions, i.e. reaction 2 ( H 2 + O ⇔ H + OH) and reaction 3 ( H 2 + OH ⇔ H + H 2 O ), are found to be responsible for the trend shown earlier in Fig. 7a, and other reactions such as reaction 78 ( NH + H ⇔ N + H 2 ) do not actually contribute to the formation of this trend.Figure 8 shows the joint PDFs between the normalized H 2 pro- duction rates of these reactions and the normalized flame front curvature for case AH040.The data are sampled at c = 0.75 iso-surface, and the hydrogen production rate for each chemical reaction has been normalized by the mean absolute value of the total hydrogen production rate evaluated for the turbulent condition at c = 0.75 iso-surface.The results show that the H 2 production rates of reactions 2 and 3 have a mild negative correlation with the flame front curvature.One of the potential reasons for the trend detected in Fig. 8a-b could be the augmentation of the total diffusion flux of H 2 in positively curved regions, see Fig. 6c.This enhancement results in increasing the hydrogen concentration in positively curved regions compared to the negatively curved zones, leading to a higher local equivalence ratio.Hence, the reactions containing molecular hydrogen react intensely with O and OH radicals in positively curved zones depleting these radical pools from these regions compared to the negatively curved zones, as shown in Fig. 9.A similar observation was recently reported by Rieth et al. (2021) for preheated NH 3 ∕H 2 ∕N 2 ∕air flames developing in a shear layer, where the conversion rate of H 2 + OH ⇔ H + H 2 O was shown to be intensi- fied in positively curved regions for the test cases at 1 and 10 atm pressure.In addition, reaction 78 (as an example for all other reactions not contributing to this trend) is responsible for producing H 2 (not consuming it), and it does not show a clear correlation with the flame front curvature for this case.It should be noted that the contribution of some of the reactions such as reaction 155 ( H + HONO ⇔ H 2 + NO 2 ) in producing H 2 are very insig- nificant for the current test case (not shown), and their examination is not discussed herein.It is worth noting that since the molecular hydrogen does not exist as a fuel component for NH 3 ∕air flame, i.e. case AH000, and the total diffusion flux of H 2 is independent of the flame front curvature, shown earlier in Fig. 6a, the trend observed in Fig. 8a-b does not hold for this case (not shown for the sake of brevity).

Concluding remarks
The density-weighted flame displacement speed statistics and the validity of extrapolation relations for estimating this quantity, the Damköhler's first hypothesis, and the influence of curved flame front regions on the local equivalence ratio and molecular diffusion fluxes were examined based on quasi-DNS results with detailed chemistry and the mixture-averaged transport model of low turbulent premixed planar NH 3 ∕air and NH 3 ∕H 2 ∕air flames in decaying turbulence under atmospheric condition.The blending ratio determining the hydrogen concentration in the ammonia/hydrogen mixture was varied from 0.0 to 0.6 in this study.The initial Karlovitz number was fixed to 4.28 for all the test cases indicating that the numerical data are located within the thin reaction zones combustion regime.It 1 3 should be emphasized that the test cases were limited to the three cases since performing the 3-D quasi-DNS of chemically reacting flows with detailed chemistry and the mixtureaveraged transport model is computationally very expensive.The main findings are summarized as follows: 1.The mean value of the density-weighted flame displacement speed ( S * d ) was similar to the unstrained premixed laminar burning velocity ( S 0 L ) for NH 3 ∕air flame, while it deviated from S 0 L for NH 3 ∕H 2 ∕air flames.This observation suggested that the flame elements, on average, moved identically to the unstrained premixed laminar burning velocity for NH 3 ∕air flame, and they propagated faster than S 0 L for NH 3 ∕H 2 ∕air flames.Furthermore, two extrapolation relations for evaluating S * d , i.e.LC and N3P models, were examined.The performances of both models were almost similar when sampling the data at c = 0.2 , while the N3P model performed better at c = 0.8 .This improved performance was due to the additional free-parameter that existed in the N3P model.2. The results showed that Damköhler's first hypothesis remained valid for the NH 3 ∕air flame indicating that the ratio of the turbulent burning velocity to the unstrained premixed laminar burning velocity was relatively equal to the ratio of the wrinkled flame surface area to the unwrinkled flame surface area.However, this hypothesis did not maintain when increasing the blending ratio ( ) to 0.4 and 0.6.Furthermore, the ratio of the turbulent burning velocity to the unstrained premixed laminar burning velocity was promoted by increasing the blending ratio.3. The mean value of the local equivalence ratio ( ) across the flame brush for the turbulent NH 3 ∕air flame was almost equal to its laminar counterpart, while it was higher for NH 3 ∕H 2 ∕air flames.The local equivalence ratio for positively curved regions was higher than the negatively curved regions for NH 3 ∕H 2 ∕air flames.This observation could be due to the local enrichment of H 2 in positively curved regions due to the preferential diffusion.4. The results indicated that the hydrogen production rate was negatively correlated with the flame front curvature, and the ammonia production rate was not dependent on the flame front curvature for NH 3 ∕air flames enriched with H 2 .This observation sug- gested that the fuel consumption rate was intensified in positively curved regions.This enhanced hydrogen consumption rate in positively curved zones could be interpreted as a higher local burning velocity value at these zones.5.The results showed that hydrogen production rates through reactions H 2 + O ⇔ H + OH and H 2 + OH ⇔ H + H 2 O have a mild negative correlation with the flame front curva- ture while other reactions such as NH + H ⇔ N + H 2 did not show any clear correlation with the flame front curvature for the test cases with enriched hydrogen in ammonia/air flames.

Fig. 1
Fig. 1 Mean values of normalized density-weighted flame displacement speed S * d ∕S 0 L (right y-axis) and its components, i.e. S * d, r ∕S 0 L , S * d, n ∕S 0 L , and S * d, t ∕S 0 L (left y-axis), with respect to the progress variable (c) for a case AH000, b case AH040, and c case AH060 shows the scatter of S * d ∕S 0 L with respect to the normalized flame front curvature ( th ).Since the statistical behavior of this property might vary depending on the choice of the progress variable value, the data are conditioned/analyzed at two different progress variable values of c = 0.2 and 0.8 for the sake of completeness.The results show a negative correlation between S * d and for all the test conditions/progress variable values studied in this work.This implies that the flame elements in negatively curved regions propagate faster than the regions with the positive flame front curvature.A similar observation was previously reported in the literature based on 2-D/3-D DNS with simple chemistry, see e.g.Nivarti and Cant (2019);Herbert et al. (2020), andTrivedi andCant (2022), and detailed chemistry, see e.g.Echekki and Chen (1999) andChakraborty et al. (

Fig. 2
Fig. 2 Scatter of S * d ∕S 0 L as well as the values of S * d ∕S 0 L evaluated using LC and N3P models with respect to the normalized flame front curvature ( th ) for case AH000 conditioned at a c = 0.2 and d c = 0.8 , case AH040 conditioned at b c = 0.2 and e c = 0.8 , and case AH060 conditioned at c c = 0.2 and f c = 0.8

Fig. 3
Fig. 3 Instantaneous 2-D slices of the local equivalence ratio ( ) for a case AH000, b case AH040, and c case AH060 at t sim = 2t e .The white, gray, and black lines correspond to the isosurfaces of c = 0.2 , 0.5, and 0.8, respectively

Fig. 5
Fig. 5 Joint PDFs between the normalized local equivalence ratio ( ∕ L ) and the normalized flame front curvature ( th ) sampled at c = 0.75 iso-surface for a case AH000, b case AH040, and c case AH060

Fig. 6
Fig. 6 Mean values of the diffusion flux terms for the molecular hydrogen for a case AH000, c case AH040, and e case AH060 with respect to th .Mean values of the diffusion flux terms for the atomic hydrogen for b case AH000, d case AH040, and f case AH060 with respect to th .The data are sampled at c = 0.75 iso-surface for the test cases.The normal diffusion, curvature diffusion, diffusion due to the gradient of the molecular weight of the mixture, and total diffusion terms are represented as MD 1 , MD 2 , MD 3 , and MD T , respectively.The diffusion flux terms have been normalized with the total diffusion flux value evaluated for the corresponding laminar flame.The joint PDFs between the normalized total diffusion flux and the normalized flame front curvature are shown for each test case

Fig. 8
Fig. 8 Joint PDFs between the normalized H 2 production rate of a reaction 2, b reaction 3, and c reaction 78 and the normalized flame front curvature ( th ) for case AH040.The data are sampled at c = 0.75 isosurface.The H 2 production rate for each reaction has been normalized by the mean absolute value of the hydrogen production rate evaluated from the turbulent test case obtained at c = 0.75 iso-surface

Table 1
experimental data with the laminar burning velocity values of 7.01, 30.37, and 68.30 (cm/s) obtained numerically fromStagni et al. ( these

Table 2
Values of L M ∕ th , C, correlation coefficients r between S * d acquired from quasi-DNS and predicted ones from LC and N3P extrapolation relations, and the goodness of fit g.o.f. for all the test conditions at t sim = 2t e

Table 3
Values of S T ∕S 0 L , A T ∕A L , and I 0 for all the test conditions at t sim = 2t e