2.1.1. Simplified mixture fraction models

Initial efforts to capture differential diffusion effects were conducted by Vreman et al. (Reference Vreman, van Oijen, de Goey and Bastiaans2009), who computed an effective Lewis number for each controlling variable. Their approach accounted for differential diffusion between heat and fuel but did not capture preferential diffusion between individual species. Moreover, the results were found to be sensitive to the choice of controlling variables.

Another approach to capture differential and preferential diffusion was proposed by Regele et al. (Reference Regele, Knudsen, Pitsch and Blanquart2013), who extended the FPV model to incorporate constant, non-unity Lewis numbers. They employed a reduced mixture fraction that considers only fuel and oxidiser, and assumed one-step chemistry for the derivation of its transport equation. This offers the advantage of a simplified closure of the diffusion terms, however, only the diffusion effects of the major species are taken into account. Note that the mixture fraction transport equation includes a cross-diffusion term that accounts for diffusion induced by the gradient of the progress variable. However, no corresponding cross-diffusion term is included in the progress variable transport equation. The authors employed a flamelet database consisting of 1D FP flames with varying equivalence ratios. Schlup & Blanquart (Reference Schlup and Blanquart2019) extended the work by Regele et al. (Reference Regele, Knudsen, Pitsch and Blanquart2013) and relaxed the constant Lewis number assumption to include mixture-averaged transport and thermal (Soret) diffusion. The model was tested on a spherically expanding flame and a TD unstable FP flame, both showing good results. A three-dimensional (3D) simulation of a turbulent premixed lean hydrogen/air flame showed good agreement in regions without superadiabatic temperatures, but deviations were observed in superadiabatic regions. Yao & Blanquart (Reference Yao and Blanquart2024) applied this model in the framework of the large eddy simulations (LES) of a low-swirl burner with a lean premixed hydrogen/air mixture using a presumed probability density function (PDF) approach for turbulence-chemistry interaction (TCI), while Berger et al. (Reference Berger, Attili, Gauding and Pitsch2025) also applied the model in the LES of a lean hydrogen/air slot flame also using a presumed PDF approach for TCI.

Böttler et al. (Reference Böttler, Chen, Xie, Scholtissek, Chen and Hasse2022, Reference Böttler, Lulic, Steinhausen, Wen, Hasse and Scholtissek2023) adopted an approach conceptually similar to that of Regele et al. (Reference Regele, Knudsen, Pitsch and Blanquart2013) and Schlup & Blanquart (Reference Schlup and Blanquart2019), employing a reduced mixture fraction based on the major species ( $\textrm {H}_{2}$ , $\textrm {O}_{2}$ and ${\textrm {H}_{2}}\textrm {O}$ ). However, instead of directly transporting this reduced mixture fraction, they solved transport equations for the major species and reconstructed both the reduced mixture fraction and the progress variable from those species prior to the table lookup. In the transport equations of the major species, the diffusion coefficients of the transported species are well defined and can be tabulated directly. The correction velocity and thermal diffusion are neglected in the transport equations, as their closure is not straightforward within this framework. As the model was originally developed for forced ignition in hydrogen/air mixtures (Böttler et al. Reference Böttler, Chen, Xie, Scholtissek, Chen and Hasse2022), a transport equation for the enthalpy $h$ was included, resulting in a 3D manifold with $Y_{{c}}$ , $Z$ and $h$ as controlling variables. For the differential diffusion term in the enthalpy equation, the gradients of the non-transported species were approximated using the corresponding gradients in the 1D flamelets. Instead of 1D FP unstretched flamelets, the composition space model (CSM) (Scholtissek et al. Reference Scholtissek, Domingo, Vervisch and Hasse2019a , Reference Scholtissek, Domingo, Vervisch and Hasseb ) was used for generating the flamelet database, enabling the representation of flame structures from various canonical premixed flame configurations, including arbitrary combinations of strain and curvature (Böttler et al. Reference Böttler, Scholtissek, Chen, Chen and Hasse2021). The model was applied to a spherically expanding TD unstable lean hydrogen/air flame and compared with a second model, also based on the CSM, in which the mixture fraction and the curvature of the flamelets have been varied (Böttler et al. Reference Böttler, Lulic, Steinhausen, Wen, Hasse and Scholtissek2023). Instead of using curvature directly as a controlling variable, the H radical was employed as a second progress variable, resulting in a 3D manifold with a reduced mixture fraction and two progress variables. The newly proposed model showed improved accuracy for TD unstable flames, although discrepancies with DC simulations remain. Schepers & van Oijen (Reference Schepers and van Oijen2025) effectively combined the two aforementioned models and replaced the major species $\textrm {O}_{2}$ with the $\textrm {H}$ radical. This substitution is motivated by the significant contribution of the $\textrm {H}$ radical to the differential diffusion of enthalpy. To construct their model, they employed a flamelet database generated from 1D FP flames with varying enthalpy and mixture fraction levels. Their approach demonstrated improved accuracy when applied to a 2D TD unstable FP flame, while discrepancies remain, for example, in the linear regime of TD instabilities.

2.1.2. Full mixture fraction models

In contrast, the foundation for the framework presented in this study, originally proposed by Oijen & Goey (Reference Oijen van and Goey de2000) and de Swart et al. (Reference de Swart, Bastiaans, van Oijen, de Goey and Cant2010) within the FGM framework (van Oijen et al. Reference van Oijen, Donini, Bastiaans, ten Thije Boonkkamp and de Goey2016), is not based on a reduced mixture fraction definition. In their formulation, constant Lewis numbers were assumed (i.e. $\textit{Le}_i$ does not vary across the flame front). They demonstrated that, under the manifold assumption, i.e. that any quantity $\xi$ depends solely on the controlling variables, the diffusion terms in the transport equations of the controlling variable can be closed, as illustrated through the following considerations.

While diffusion terms could be computed directly during simulation using the manifold, this would require storing all species mass fractions and diffusion coefficients, causing significant memory and computational overhead. Instead, gradients are partially precomputed in the manifold, greatly reducing the number of variables retrieved from the table. For any manifold quantity $\xi (\mathcal{Y}_1, \ldots , \mathcal{Y}_n)$ (e.g. a species mass fraction $Y_k$ ), the total differential is defined as

(2.6) \begin{equation} \frac {\partial \xi }{\partial x_i} = \sum _{j=1}^{N_{{c}}} \left .\frac {\partial \xi }{\partial \mathcal{Y}_{\!j}}\right |_{\mathcal{Y}_1,\ldots ,\mathcal{Y}_{j-1},\mathcal{Y}_{j+1},\ldots ,\mathcal{Y}_{n}} \frac {\partial \mathcal{Y}_{\!j}}{\partial x_i}, \end{equation}

where gradients with respect to the controlling variables $\mathcal{Y}_{\!j}$ are precomputed within the manifold, while spatial gradients of the controlling variables are obtained during the simulation. (Note that in the following, $\left .{\partial \xi }/{\partial \mathcal{Y}_{\!j}}\right |_{\mathcal{Y}_1,\ldots ,\mathcal{Y}_{j-1},\mathcal{Y}_{j+1},\ldots ,\mathcal{Y}_{n}}$ is abbreviated as ${\partial \xi }/{\partial \mathcal{Y}_{\!j}}$ .) Various approaches exist for this purpose, and the method employed in this work for calculating the gradients with respect to the controlling variables is discussed in detail in the following section. Assuming constant Lewis numbers, the diffusion flux is defined as $Y_{k} V_{k,i} = - {\alpha }/{{\textit{Le}_k}}\,{\partial Y_k}/{\partial x_i}$ , which can be reformulated as $Y_k V_{k,i} = - \sum _{l=1}^{N_{{c}}} ({\alpha }/{\textit{Le}_k} \, {\partial Y_k}{/\partial \mathcal{Y}_l}) {\partial \mathcal{Y}_l} / \partial x_i = - \sum _{l=1}^{N_{{c}}} \varLambda _{Y_k,\mathcal{Y}_l} \partial \mathcal{Y}_{l} / \partial x_i$ using (2.6).

For an exemplary controlling variable $\mathcal{Y}_{\!j}$ , the transport equation is given by

(2.7) \begin{equation} \frac {\partial \rho \mathcal{Y}_{\!j}}{\partial t} + \frac {\partial \rho u_i \mathcal{Y}_{\!j}}{\partial x_i} = \frac {\partial }{\partial x_i}\! \left (\rho \sum _{k=1}^{N_{{s}}} \varsigma _{j,k} Y_k V_{k,i} \right ) + S_{\mathcal{Y}_{\!j}}, \end{equation}

with the source term $S_{\mathcal{Y}_{\!j}}$ and the weighting factor $\varsigma _{j,k}$ of species $k$ associated with the controlling variable $\mathcal{Y}_{\!j}$ , which, for example, corresponds to the weighting factor $a_{k}$ in the case of a reaction progress variable. The diffusion term in the $\mathcal{Y}_{\!j}$ transport equation can be reformulated as

(2.8) \begin{align} - \rho \sum _{k=1}^{N_{{s}}} \varsigma _{j,k} Y_k V_{k,i} = \rho \sum _{k=1}^{N_{{s}}} \varsigma _{j,k}\! \left (-Y_k V_{k,i}\right ) - \frac {\kappa }{c_p} \frac {\partial \mathcal{Y}_{\!j}}{\partial x_i} + \frac {\kappa }{c_p} \frac {\partial \mathcal{Y}_{\!j}}{\partial x_i} \notag \\ = \underbrace {\rho \sum _{k=1}^{N_{{s}}} \varsigma _{j,k}\! \underbrace {\left ( -Y_k V_{k,i} - \frac {\kappa }{\rho c_{\!p}} \frac {\partial {Y}_{k}}{\partial x_i} \right )}_{j_{k,i}}}_{\textit{differential diffusion}} + \underbrace {\frac {\kappa }{c_{\!p}} \frac {\partial \mathcal{Y}_{\!j}}{\partial x_i}}_{\textit{unity Lewis diffusion}}, \end{align}

which is decomposed into a unity Lewis number contribution and a contribution accounting for differential diffusion, as proposed, for example, in the work of Nicolai et al. (Reference Nicolai, Dressler, Janicka and Hasse2022). Using (2.6) (and the closure of the diffusion flux proposed above), the differential diffusion flux is obtained as

(2.9) \begin{equation} \sum _{k=1}^{N_{{s}}}\varsigma _{j,k} j_{k,i} = \sum _{l=1}^{N_{{c}}} \varGamma _{\mathcal{Y}_{\!j},\mathcal{Y}_{l}}\frac {\partial \mathcal{Y}_l}{\partial x_i}, \end{equation}

where $\varGamma _{\mathcal{Y}_{\!j},\mathcal{Y}_l}$ denotes the contribution of the controlling variable $\mathcal{Y}_{l}$ on $\mathcal{Y}_{\!j}$ :

(2.10) \begin{equation} \varGamma _{\mathcal{Y}_{c},\mathcal{Y}_{l}} = \sum _{k=1}^{N_{{s}}} \varsigma _{j,k} \varLambda _{Y_{k},\mathcal{Y}_l}^{*} = \sum _{k=1}^{N_{{s}}} \varsigma _{j,k}\! \left ( \varLambda _{Y_{k},\mathcal{Y}_l} - \frac {\kappa }{\rho c_p} \frac {\partial {Y}_k}{\partial \mathcal{Y}_l} \right )\!. \end{equation}

Finally, the transport equation for the controlling variable $\mathcal{Y}_{\!j}$ reads

(2.11) \begin{equation} \frac {\partial \rho \mathcal{Y}_{\!j}}{\partial t} + \frac {\partial \rho u_i \mathcal{Y}_{\!j}}{\partial x_i} = \frac {\partial }{\partial x_i}\! \left (\rho \sum _{l=1}^{N_{{c}}} \varGamma _{\mathcal{Y}_{\!j},\mathcal{Y}_l}\frac {\partial \mathcal{Y}_l}{\partial x_i} + \frac {\kappa }{c_p} \frac {\partial \mathcal{Y}_{\!j}}{\partial x_i} \right ) + S_{\mathcal{Y}_{\!j}}. \end{equation}

Based on these considerations, it becomes evident that the diffusion terms consist of two components: normal diffusion aligned with the gradient of each controlling variable and cross-diffusion (or drift) terms directed along the gradients of the other controlling variables. These terms effectively account for multidimensional differential and preferential diffusion effects. For instance, in a 2D manifold parametrised by the progress variable and the mixture fraction, the transport equation for each variable includes not only its own diffusion term but also an additional cross-diffusion term involving the gradient of the other controlling variable.

Some notes regarding the diffusion terms are as follows:

1. (i) In the original approach of Oijen & Goey (Reference Oijen van and Goey de2000), de Swart et al. (Reference de Swart, Bastiaans, van Oijen, de Goey and Cant2010), however, it was additionally assumed that the controlling variables $\mathcal{Y}_i$ depend locally only on the progress variable $Y_{{c}}$ , i.e. $\mathcal{Y}_i = \mathcal{Y}^{\textrm{1D}}_{i}(Y_{{c}})$ . This assumption effectively decouples the equations (as the cross-diffusion term in the progress variable transport equation vanishes) and implies that differential diffusion occurs only in the direction of the gradient of the progress variable. Donini et al. (Reference Donini, Bastiaans, van Oijen and de Goey2015) introduced enthalpy as an additional controlling variable and adopted the same assumptions as de Swart et al. (Reference de Swart, Bastiaans, van Oijen, de Goey and Cant2010). Mukundakumar et al. (Reference Mukundakumar, Efimov, Beishuizen and van Oijen2021) eliminated the additional assumptions by employing a mathematical reformulation that avoids the need to compute gradients within the manifold. Instead, additional terms are stored in the manifold, whose gradients are evaluated during the simulation. However, this approach is only valid under the assumption of constant Lewis numbers, as it relies on the spatial gradients of the Lewis numbers being zero.

2. (ii) Abtahizadeh et al. (Reference Abtahizadeh, de Goey and van Oijen2015) and Nicolai et al. (Reference Nicolai, Dressler, Janicka and Hasse2022) also relaxed the assumptions and included cross-diffusion fluxes between all controlling variables in the context of the autoignition of $\textrm {CH}_{4}$ / $\textrm {H}_{2}$ flames (Abtahizadeh et al. Reference Abtahizadeh, de Goey and van Oijen2015) and turbulent stratified $\textrm {CH}_{4}$ / $\textrm {H}_{2}$ flames (Nicolai et al. Reference Nicolai, Dressler, Janicka and Hasse2022). However, as both studies also assumed constant, non-unity Lewis numbers, Pérez-Sánchez et al. (Reference Pérez-Sánchez, Fortes and Mira2025) recently extended the model to incorporate a mixture-averaged diffusion formulation that also accounts for contributions to the diffusion fluxes arising from the correction velocity. Their flamelet database consists of FP flames with varying equivalence ratios and the extended model was tested on a stratified hydrogen/air triple flame and showed very accurate results. It was subsequently applied to 2D TD unstable FP hydrogen/air flames (Fortes et al. Reference Fortes, Pérez-Sánchez, Both, Grenga and Mira2025), where it also yielded good agreement with DC simulations with minor deviations, for example, in numerically obtained dispersion relations.

The framework employed in this work should incorporate the least restrictive assumptions in the model derivation, while ensuring the robust computation of diffusion terms, as the FWI introduces an additional effect that increases the complexity and further challenges the model. The least restrictive assumptions are provided by the previously introduced framework (Nicolai et al. Reference Nicolai, Dressler, Janicka and Hasse2022; Fortes et al. Reference Fortes, Pérez-Sánchez, Both, Grenga and Mira2025; Pérez-Sánchez et al. Reference Pérez-Sánchez, Fortes and Mira2025), although these studies have focused exclusively on FP flames rather than on FWIs. Thus, the following section presents an extended modelling framework, which integrates a mixture-averaged diffusion model including thermal diffusion and accounts for heat losses via an additional controlling variable alongside the reaction progress variable and mixture fraction, usually employed for TC approaches focusing on TD unstable flames.

2.2.1. Closure of the diffusion terms

In this study, differential diffusion effects are taken into account by modelling the diffusive fluxes using the mixture-averaged approximation, including thermal (Soret) diffusion. Note that, in principle, any diffusion model can be used. In the supplementary material, the extension to the multicomponent diffusion model, including multicomponent thermal diffusion, is demonstrated and validated. Employing the mixture-averaged approximation (Hirschfelder, Bird & Curtiss Reference Hirschfelder, Bird and Curtiss1964), the diffusion flux in the transport equations is expressed as

(2.12) \begin{equation} \rho Y_k V_{k,i} = \rho Y_k V_{k,i}^{{D}} + \rho Y_k V_{k,i}^{{T}} + \rho Y_k V_{k,i}^{{C}}, \end{equation}

where $V_{k,i}^{{D}}$ is the mixture-averaged diffusion velocity, $V_{k,i}^{{T}}$ is the thermal diffusion velocity and $V_{k,i}^{{C}}$ is the correction velocity, which is required for the non-mass conservative mixture-averaged approximation. Thermal diffusion is considered and is modelled using a simplified formulation from Chapman & Cowling (Reference Chapman and Cowling1990). It is applied only to light species with molecular weights $W_k \lt 5$ (specifically H and $\textrm {H}_{2}$ ; see Reaction Design 2015; Schlup & Blanquart Reference Schlup and Blanquart2018; Howarth et al. Reference Howarth, Day, Pitsch and Aspden2024 for further details). Details on the formulation of the mixture-averaged diffusion velocity, thermal diffusion velocity and correction velocity are provided in the supplementary material. Overall, the diffusion flux $Y_kV_{k,i}$ of species $k$ in direction $i$ is given by

(2.13) \begin{align} Y_k V_{k,i} &= - D_{\textit{mix},k} \frac {\partial Y_k}{\partial x_i} - \frac {Y_k D_{\textit{mix},k}}{\overline {W}} \frac {\partial \overline {W}}{\partial x_i} - \frac {D_{\textit{therm},k}}{\rho T} \frac {\partial T}{\partial x_i} \notag \\ &\quad + Y_k \sum _{j=1}^{N_{{s}}} D_{\textit{mix},j} \frac {\partial Y_{\!j}}{\partial x_i} + \frac {Y_k}{\overline {W}} \frac {\partial \overline {W}}{\partial x_i} \sum _{j=1}^{N_{{s}}} D_{\textit{mix},j} Y_{\!j} \notag \\ &\quad + \frac {Y_k}{\rho T} \frac {\partial T}{\partial x_i} \sum _{j=1}^{N_{{s}}} D_{\textit{therm},j}. \end{align}

Based on (2.6), similar to de Swart et al. (Reference de Swart, Bastiaans, van Oijen, de Goey and Cant2010), Nicolai et al. (Reference Nicolai, Dressler, Janicka and Hasse2022), Pérez-Sánchez et al. (Reference Pérez-Sánchez, Fortes and Mira2025), the diffusion flux can be reformulated as

(2.14) \begin{align} Y_k V_{k,i} &= - \sum _{l=1}^{N_{{c}}} \Bigg ( D_{\textit{mix},k} \frac {\partial Y_k}{\partial \mathcal{Y}_l} + \frac {Y_k D_{\textit{mix},k}}{\overline {W}} \frac {\partial \overline {W}}{\partial \mathcal{Y}_l} + \frac {D_{\textit{therm},k}}{\rho T} \frac {\partial T}{\partial \mathcal{Y}_l} \notag \\ &\quad - Y_k \sum _{j=1}^{N_{{s}}} D_{\textit{mix},j} \frac {\partial Y_{\!j}}{\partial \mathcal{Y}_l} - \frac {Y_k}{\overline {W}} \frac {\partial \overline {W}}{\partial \mathcal{Y}_l} \sum _{j=1}^{N_{{s}}} D_{\textit{mix},j} Y_{\!j} \notag \\ &\quad - \frac {Y_k}{\rho T} \frac {\partial T}{\partial \mathcal{Y}_l} \sum _{j=1}^{N_{{s}}} {D_{\textit{therm},j}} \Bigg ) \frac {\partial \mathcal{Y}_l}{\partial x_i} \notag \\ &\quad = - \sum _{l=1}^{N_{{c}}} \varLambda _{Y_{k},\mathcal{Y}_l} \frac {\partial \mathcal{Y}_l}{\partial x_i}. \end{align}

This flux is included in the transport equations of the respective controlling variable (e.g. progress variable $Y_{{c}}$ ).

2.2.2. Extension to heat losses

To capture heat losses in flamelet manifolds, the flamelet database must be extended to account for heat losses within the flamelets. To distinguish these additional states, the enthalpy $h$ is typically employed as an additional controlling variable. However, in this study, temperature $T$ has proven to be more suitable for the construction of manifolds involving differential diffusion and heat losses, which is discussed in detail in § 4. The transport equation for temperature can be derived from the enthalpy transport equation and is given by

(2.15) \begin{align} c_{\!p} \frac {\partial \rho T}{\partial t} + c_{\!p} \frac {\partial \rho u_i T}{\partial x_i}\nonumber \\ &\quad = \frac {\partial }{\partial x_i}\! \left (\kappa \frac {\partial T}{\partial x_i}\right ) -\! \left ( \rho \sum _{k=1}^{N_{{s}}}c_{{{p}},k} Y_k V_{k,i} \right ) \frac {\partial T}{\partial x_i} + \dot {\omega }_T^{\prime }, \end{align}

where $c_{\!p,k}$ is the specific heat capacity of species $k$ at constant pressure, $\kappa$ is the thermal conductivity and the heat release rate is defined as $\dot {\omega }_T^{\prime } = - \sum _{k=1}^{N_{{s}}}h_k \dot {\omega }_k$ .

2.2.3. Controlling variable transport equations

Finally, the reformulated transport equations for the (up to) three controlling variables used in the manifolds employed in this study are provided.

The transport equation for the reaction progress variable $Y_{{c}}$ reads

(2.16) \begin{align} \frac {\partial (\rho Y_{c})}{\partial t} + \frac {\partial (\rho u_i Y_{c})}{\partial x_i} &= \frac {\partial }{\partial x_i}\! \underbrace {\left ( \rho \sum _{l=1}^{N_{c}} \varGamma _{Y_{c},\mathcal Y_l} \frac {\partial \mathcal Y_l}{\partial x_i} + \frac {\kappa }{c_p} \frac {\partial Y_{c}}{\partial x_i} \right )}_{= -\,\rho \sum _{k=1}^{N_{s}} a_k Y_k V_{k,i}} + \rho \dot {\omega }_{Y_{c}}, \end{align}

where $\dot {\omega }_{Y_{{c}}} = \sum _{k=1}^{N_{{s}}} a_k \dot {\omega }_k$ and $\varGamma _{Y_{{c}},\mathcal{Y}_l}$ denotes the contribution of the controlling variable $\mathcal{Y}_l$ on $Y_{{c}}$ :

(2.17) \begin{equation} \varGamma _{Y_{{c}},\mathcal{Y}_l} = \sum _{k=1}^{N_{{s}}} a_{k} \varLambda _{Y_{k},\mathcal{Y}_l}^{*} = \sum _{k=1}^{N_{{s}}} a_{k}\! \left ( \varLambda _{Y_{k},\mathcal{Y}_l} - \frac {\kappa }{\rho c_p} \frac {\partial Y_{{c}}}{\partial \mathcal{Y}_l} \right )\!. \end{equation}

For the transport equation of the elemental mass fraction $Z_e$ , it follows that

(2.18) \begin{align} \frac {\partial (\rho Z_e)}{\partial t} + \frac {\partial (\rho u_i Z_e)}{\partial x_i} &= \frac {\partial }{\partial x_i}\! \underbrace {\left ( \rho \sum _{l=1}^{N_{c}} \varGamma _{Z_e,\mathcal Y_l}\, \frac {\partial \mathcal Y_l}{\partial x_i} + \frac {\kappa }{c_p}\, \frac {\partial Z_e}{\partial x_i} \right )}_{=- \rho \sum _{k=1}^{N_{s}} \gamma _{e,k}\,\dfrac {W_e}{W_k}\,Y_k V_{k,i}} , \end{align}

where $\varGamma _{Z_{e},\mathcal{Y}_l}$ denotes the contribution of the controlling variable $\mathcal{Y}_l$ to $Z_{e}$ :

(2.19) \begin{equation} \varGamma _{Z_{e},\mathcal{Y}_l} = W_e\sum _{k=1}^{N_{{s}}} \frac {\gamma _{e,k}}{W_k} \varLambda _{Y_{k},\mathcal{Y}_l}^{*} = W_e \sum _{k=1}^{N_{{s}}} \frac {\gamma _{e,k}}{W_k}\! \left ( \varLambda _{Y_{k},\mathcal{Y}_l} - \frac {\kappa }{\rho c_p} \frac {\partial Y_k}{\partial \mathcal{Y}_l} \right )\!. \end{equation}

The transport equation of the temperature $T$ reads

(2.20) \begin{equation} c_{\!p} \frac {\partial \rho T}{\partial t} + c_{\!p} \frac {\partial \rho u_i T}{\partial x_i} = \frac {\partial }{\partial x_i}{\left (\kappa \frac {\partial T}{\partial x_i}\right )} + \rho\! \underbrace {\left ( \sum _{l=1}^{N_{{c}}} \varGamma _{T,\mathcal{Y}_l} \frac {\partial \mathcal{Y}_l}{\partial x_i} \right )}_{=\sum _{k=1}^{N_{{s}}} c_{\!p,k} Y_k V_{k,i} } \frac {\partial T}{\partial x_i} + \dot {\omega }_T^{\prime },\end{equation}

with the contribution $\varGamma _{T,\mathcal{Y}_l}$ of the controlling variable $\mathcal{Y}_l$ to the temperature $T$ defined as

(2.21) \begin{equation} \varGamma _{T,\mathcal{Y}_l} = \sum _{k=1}^{N_{{s}}} c_{\!p,k} \varLambda _{Y_k,\mathcal{Y}_l}. \end{equation}

To apply these transport equations in fully coupled simulations, the coefficients $\varGamma _{\mathcal{Y}_{\!j},\mathcal{Y}_l}$ must be computed from the gradients $\partial \xi / \partial \mathcal{Y}_k$ within the manifold. These gradients are efficiently evaluated during the construction of the tabulated flamelet manifold prior to the simulation, as detailed below.

2.2.4. Calculation of gradients in the manifold

The following section outlines the calculation of gradients of quantities $\xi$ in the manifold for a general manifold with controlling variables $\mathcal{Y}_1$ to $\mathcal{Y}_n$ , as required in (2.14). An evaluation on the final manifold is not straightforward, as it is typically constructed in terms of normalised controlling variables. Therefore, a robust, generalised method is required for the efficient calculation of gradients with respect to the controlling variables. While Pérez-Sánchez et al. (Reference Pérez-Sánchez, Fortes and Mira2025) proposed an approach for the controlling variables mixture fraction and progress variable, a general framework has not yet been established in the literature. Note that the proposed procedure also relies solely on the manifold assumption (i.e. that each quantity $\xi$ depends exclusively on the controlling variables) and introduces no further assumptions.

The controlling variables are normalised and stored on a regular grid to enable direct, search-free access and straightforward interpolation within the tabulated manifold. In this context, the normalised controlling variables ( $\mathcal{Y}_{\textit{norm},1}$ to $\mathcal{Y}_{\textit{norm},n}$ ) depend only on the preceding controlling variables through their respective normalisation bounds. They are independent of subsequent variables, as such dependencies would require iterative table access, something generally avoided in tabulated manifold approaches for performance reasons. Thus, a controlling variable $\mathcal{Y}_n$ is normalised as

(2.22) \begin{align} \mathcal{Y}_{\textit{norm},n} = \frac {\mathcal{Y}_n - \mathcal{Y}_{n,\textit{min}}(\mathcal{Y}_{1}, \ldots , \mathcal{Y}_{n-1})}{\mathcal{Y}_{n,max }(\mathcal{Y}_{1}, \ldots , \mathcal{Y}_{n-1}) - \mathcal{Y}_{n,\textit{min}}(\mathcal{Y}_{1}, \ldots , \mathcal{Y}_{n-1})}\, . \end{align}

Next, the manifold assumption $\xi = \xi (\mathcal{Y}_1,\mathcal{Y}_2,\ldots ,\mathcal{Y}_{n})$ is exploited. Applying the chain rule yields a coordinate transformation for any quantity $\xi$ (note that $\big |_{\mathcal{Y}_{\textit{norm},1},\ldots ,\mathcal{Y}_{\textit{norm},j-1},\mathcal{Y}_{\textit{norm},j+1},\ldots ,\mathcal{Y}_{\textit{norm},n}}$ is omitted here for brevity in the derivatives with respect to the normalised controlling variable $\mathcal{Y}_{\textit{norm},j}$ ):

(2.23) \begin{align} \frac {\partial \xi }{\partial \mathcal{Y}_{\textit{norm},j}} &= \frac {\partial \mathcal{Y}_1}{\partial \mathcal{Y}_{\textit{norm},j}} \frac {\partial \xi }{\partial \mathcal{Y}_{1}} + \frac {\partial \mathcal{Y}_2}{\partial \mathcal{Y}_{\textit{norm},j}} \frac {\partial \xi }{\partial \mathcal{Y}_{2}} \notag \\ &\quad + \ldots + \frac {\partial \mathcal{Y}_n}{\partial \mathcal{Y}_{\textit{norm},j}} \frac {\partial \xi }{\partial \mathcal{Y}_{n}} \, . \end{align}

This can be expressed in matrix–vector form for all controlling variables $\mathcal{Y}_i$ :

(2.24) \begin{equation} \underbrace { \begin{bmatrix} \dfrac {\partial \xi _i}{\partial \mathcal{Y}_{\textit{norm},1}} \\[12pt] \dfrac {\partial \xi _i}{\partial \mathcal{Y}_{\textit{norm},2}} \\[12pt] \vdots \\[12pt] \dfrac {\partial \xi _i}{\partial \mathcal{Y}_{\textit{norm},n}} \end{bmatrix} }_{\underline {g}_{\textit{norm}}} = \underbrace { \begin{bmatrix} \dfrac {\partial \mathcal{Y}_1}{\partial \mathcal{Y}_{\textit{norm},1}} & \dfrac {\partial \mathcal{Y}_2}{\partial \mathcal{Y}_{\textit{norm},1}} & \ldots & \dfrac {\partial \mathcal{Y}_n}{\partial \mathcal{Y}_{\textit{norm},1}} \\[12pt] \dfrac {\partial \mathcal{Y}_1}{\partial \mathcal{Y}_{\textit{norm},2}} & \dfrac {\partial \mathcal{Y}_2}{\partial \mathcal{Y}_{\textit{norm},2}} & \ldots & \dfrac {\partial \mathcal{Y}_n}{\partial \mathcal{Y}_{\textit{norm},2}} \\[12pt] \vdots & \vdots & \ddots & \vdots \\[12pt] \dfrac {\partial \mathcal{Y}_1}{\partial \mathcal{Y}_{\textit{norm},n}} & \dfrac {\partial \mathcal{Y}_2}{\partial \mathcal{Y}_{\textit{norm},n}} & \ldots & \dfrac {\partial \mathcal{Y}_n}{\partial \mathcal{Y}_{\textit{norm},n}} \end{bmatrix} }_{\underline {\underline {C}}} \underbrace { \begin{bmatrix} \dfrac {\partial \xi _i}{\partial \mathcal{Y}_1} \\[12pt] \dfrac {\partial \xi _i}{\partial \mathcal{Y}_2} \\[12pt] \vdots \\[12pt] \dfrac {\partial \xi _i}{\partial \mathcal{Y}_n} \end{bmatrix} }_{\underline {g}}\!. \end{equation}

This expression can be inverted to obtain $\underline {g}$ from $\underline {g}_{\textit{norm}}$ and $\underline {\underline {C}}$ , whose entries can be directly computed on the (normalised) manifold:

(2.25) \begin{equation} \underline {g} = \underline {\underline {C}}^{-1} \boldsymbol{\cdot }\underline {g}_{\textit{norm}}. \end{equation}

In this case, a linear system of equations must be solved. However, the relation in (2.22) can be exploited: if it holds for all controlling variables, the resulting matrix $\underline {\underline {C}}$ becomes upper triangular, since

(2.26) \begin{equation} \frac {\partial \mathcal{Y}_{\!j}}{\partial \mathcal{Y}_{\textit{norm},k}} = 0 \quad \text{for all }\quad k \gt j . \end{equation}

In this case, $\underline {g}$ can be directly determined from $\underline {g}_{\textit{norm}}$ and $\underline {\underline {C}}$ through back substitution.

Furthermore, the entries of the matrix $\underline {\underline {C}}$ can be simplified through analytical derivation. The diagonal entries can be determined analytically from (2.22) as

(2.27) \begin{equation} \frac {\partial \mathcal{Y}_{\!j}}{\partial \mathcal{Y}_{\textit{norm},j}} = \mathcal{Y}_{j,max } - \mathcal{Y}_{j,min }, \end{equation}

while the remaining non-diagonal entries can also be determined from a reformulation of (2.22):

(2.28) \begin{equation} \frac {\partial \mathcal{Y}_{\!j}}{\partial \mathcal{Y}_{\textit{norm},k}} = \mathcal{Y}_{\textit{norm},j}\! \left ( \frac {\partial \mathcal{Y}_{j,\textit{max}}}{\partial \mathcal{Y}_{\textit{norm},k}} - \frac {\partial \mathcal{Y}_{j,\textit{min}}}{\partial \mathcal{Y}_{\textit{norm},k}} \right ) + \frac {\partial \mathcal{Y}_{j,\textit{min}}}{\partial \mathcal{Y}_{\textit{norm},k}}. \end{equation}

These mathematical considerations enable the efficient calculation of the derivatives of any quantity $\xi$ with respect to the controlling variables $\mathcal{Y}_i$ on manifolds of arbitrary dimension. This aspect will be discussed in more detail in the context of the specific manifolds used in this study (§ 4).

5.5.1. Linear regime

The linear regime is analysed via numerical dispersion relations, focusing on the growth rate of the flame amplitude in response to sinusoidal perturbations of the form $A_0 \sin (2\pi y / \lambda )$ with varying wavelengths $\lambda = 2\pi / k$ , where the domain width exactly matches one wavelength of the imposed perturbation ( $L_y = \lambda$ ; see figure 3 c). For further details on numerically obtained dispersion relations, the interested reader is referred to Altantzis et al. (Reference Altantzis, Frouzakis, Tomboulides, Kerkemeier and Boulouchos2011, Reference Altantzis, Frouzakis, Tomboulides, Matalon and Boulouchos2012), Berger et al. (Reference Berger, Attili and Pitsch2022a ).

Figure 8. Dispersion relations from DC simulations and fully coupled TC simulations with the different manifolds. Symbols refer to the growth rates extracted from simulations, solid lines represent cubic spline fits to these growth rates. The growth rates $\omega$ and wavelengths $\lambda$ are normalised by the laminar flame time $\tau _{{l}} = \delta _{T,{l}} / s_{{l}}$ and laminar (thermal) flame thickness $\delta _{T,{{l}}}$ of the corresponding unstretched 1D flamelet. The grey dashed line indicates the growth rate of the DL (hydrodynamic) instability (Matalon Reference Matalon2007).

Figure 8 shows numerical dispersion relations obtained from DC simulations and TC simulations with the different manifolds. For the simulations at varying wavelengths $\lambda$ , the growth rate $\omega$ is determined by tracking the temporal evolution of the maximum displacement of the flame front, as detailed in Altantzis et al. (Reference Altantzis, Frouzakis, Tomboulides, Kerkemeier and Boulouchos2011), and shown as a function of the initial perturbation wavelength $\lambda$ . The manifolds whose flamelet databases include flamelets with varying equivalence ratios $\varphi$ of the unburnt mixture yield highly accurate results in the fully coupled TC simulations, closely matching those of the DC simulations. Only the maximum growth rate, $\omega _{\textit{max}} = \max (\omega (\lambda ))$ , and the cutoff wavelength, $\lambda _{{c}}$ , which marks the transition point where the dispersion relation changes sign at low wavelengths (i.e. high wavenumbers), are slightly underpredicted. In contrast, the manifolds M-EGR and M-HOQ significantly underestimate both the maximum growth rate $\omega _{\textit{max}}$ and the corresponding wavelength at which this maximum occurs, $\lambda _{\omega _{\textit{max}}} = {\arg}\,\max (\omega (\lambda ))$ . While manifold M-HOQ still provides an accurate prediction of the cutoff wavelength $\lambda _{{c}}$ , manifold M-EGR overestimates this value (i.e. underestimates the corresponding wavenumber). The manifolds that only capture heat losses (M-EGR and M-HOQ) are therefore not suitable for the accurate prediction of TD instabilities. This limitation is attributed to the absence of mixture variation in the flamelet databases of M-EGR and M-HOQ. Furthermore, this implies that enthalpy variation alone is not sufficient to capture TD instabilities.

5.5.2. Nonlinear regime

Next, the nonlinear regime is investigated. The setup for the nonlinear regime is similar to that used for the linear regime; however, the lateral domain size is set to $L_y = 100 \delta _{T,{l}} \gt 30 \lambda _{{c}}$ to ensure that no relevant length scales of the TD instabilities are suppressed (Berger et al. Reference Berger, Kleinheinz, Attili and Pitsch2019; Lapenna, Lamioni & Creta Reference Lapenna, Lamioni and Creta2021). The streamwise domain length is chosen to be sufficiently large ( $L_x = 150 \delta _{T,{l}}$ ) to guarantee that all species and temperature profiles (as well as the profiles of the controlling variables in the TC simulations) exhibit no gradients at the domain boundaries.

Figure 9. Temperature $T$ (fist row, left), elemental mass fraction $Z_{\textrm {H}}$ (first row, right), $\textrm {H}_{2}$ source term $\dot {\omega }_{\textrm {H}_{2}}$ (second row), species profiles of ${\textrm {H}_{2}}\textrm {O}$ (third row), $\textrm {H}$ (fourth row) and $\textrm {OH}$ (fifth row) for a snapshot of a 2D TD unstable flame in the nonlinear regime. The relative errors $\varepsilon (\xi )$ obtained from an a-priori lookup using the various manifolds are displayed on the right of each respective quantity $\xi$ .

A-priori analysis: first, an a-priori analysis is conducted for the different manifolds. Figure 9 shows various thermochemical quantities (e.g. temperature $T$ , ${\textrm {H}_{2}}\textrm {O}$ mass fraction $Y_{{\textrm {H}_{2}}\textrm {O}}$ ) for a snapshot of the DC simulation of the TD unstable flame in the nonlinear regime and the respective relative errors $\varepsilon (\xi )$ from an a-priori lookup with the manifolds, as defined in § 5.2. First, the manifolds that account for mixture variations (M-TD, M-EGR-TD and M-HOQ-TD) are discussed, followed by the manifolds without mixture variation in the flamelet database (M-HOQ and M-EGR).

At the top left of figure 9, in addition to the temperature from the DC simulation used as a reference, the relative temperature error $\varepsilon (T)$ for manifold M-TD is shown, exhibiting a maximum error of only $3\,\%$ . For the other manifolds, temperature serves as a controlling variable. The second row shows the $Y_{\textrm {H}_{2}}$ source term, $\dot {\omega }_{\textrm {H}_{2}}$ , together with the relative a-priori errors $\varepsilon (Y_{{\textrm {H}_{2}}\textrm {O}})$ for the manifolds including mixture variation (i.e. M-TD, M-EGR-TD and M-HOQ-TD). Local errors of up to $8\,\%$ occur in regions of strong positive curvature, demonstrating the manifolds’ consistency with the reference DC simulation and indicating good expected accuracy in the subsequent a-posteriori validation. The third row shows $Y_{{\textrm {H}_{2}}\textrm {O}}$ together with the corresponding relative a-priori errors. Errors below $0.2\,\%$ for the manifolds including mixture variation indicate very good agreement for this major species; similar accuracy is obtained for $Y_{\textrm {O}_{2}}$ , which is not shown for brevity. For the minor species $\textrm {H}$ and $\textrm {OH}$ (fourth and fifth rows, respectively), the manifolds M-TD, M-EGR-TD and M-HOQ-TD exhibit local mass fraction errors of up to $8\,\%$ , primarily in regions of high flame curvature, similar to the source term of $Y_{\textrm {H}_{2}}$ .

The relative error of $Z_{\textrm {H}}$ for manifold M-HOQ (top right of figure 9) clearly demonstrates its inability to predict local mixture variations resulting from differential and preferential diffusion. This limitation is also evident in all other quantities. For instance, the $Y_{\textrm {H}_{2}}$ source term $\dot {\omega }_{\textrm {H}_{2}}$ deviates by more than $50\,\%$ from the reference solution in some regions, and the mass fractions of minor species such as $\textrm {H}$ and $\textrm {OH}$ exhibit errors of similar magnitude. Even major species such as ${\textrm {H}_{2}}\textrm {O}$ exhibit relative errors exceeding $20\,\%$ at various locations. It is noted that the results from the a-priori lookup of the M-EGR manifold are similar to those of M-HOQ (as this manifold is likewise unable to capture the TD instability) and are therefore omitted for brevity.

Figure 10. Time-averaged relative error of the flame consumption speed $\varepsilon (s_{{c}})$ , obtained from an a-priori analysis of all manifolds with respect to the reference DC simulation.

Figure 10 shows the time-averaged relative error of the a-priori flame consumption speed $s_{{c}}$ , averaged over more than $100$ snapshots, as defined by Remiddi et al. (Reference Remiddi, Lapenna, Cavalieri, Schintu, Indelicato, Attili, Berger, Pitsch and Creta2024), in comparison with the reference DC simulation for each manifold.

The flame consumption speed $s_{{c}}$ is defined as

(5.2) \begin{equation} s_{{c}} = - \frac {1}{\rho _{{u}} Y_{\textrm {H}_{2},{u}} L_y} \int _{0}^{L_y}\int _{0}^{L_x}\dot {\omega }_{\textrm {H}_{2}} {\rm d}x\, {\rm d}y. \end{equation}

For the manifolds, the corresponding thermochemical quantities, such as the source term of $\textrm {H}_{2}$ , are taken from the a-priori lookups. The relative errors in flame consumption speed $\varepsilon (s_{{c}})$ also exhibit only minor deviations for the manifolds whose flamelet databases contain flamelets at varying equivalence ratios $\varphi$ (i.e. M-TD, M-EGR-TD and M-HOQ-TD); they show a slight underprediction of approximately $4\,\%$ relative to the DC reference solution. It is particularly noteworthy that the relative errors from an a-priori analysis using a manifold constructed directly from a DC simulation of a TD unstable flame fall within the same range, as shown in Remiddi et al. (Reference Remiddi, Lapenna, Cavalieri, Schintu, Indelicato, Attili, Berger, Pitsch and Creta2024). Additionally, it should be noted that the relative error in flame consumption speed in an a-priori analysis can be directly associated with the error in the reactivity factor (normalised mean local flame speed) $I_0 = ({s_{{c}} / s_{{l}}})/( {A / L_y})$ (Bray & Cant Reference Bray and Cant1991), since the flame front surface area (i.e. length in two dimensions) $A/L_y$ remains unchanged in the a-priori assessment. It should also be noted that local errors in the source term $\dot {\omega }_{\textrm {H}_{2}}$ can lead to deviations in the local flame speed, for example, in regions of strong positive and negative curvature, which may locally differ from the (mean) errors in the overall consumption speed and reactivity factor across the domain; see also the a-priori analysis of the source term $\dot {\omega }_{\textrm {H}_{2}}$ . The other two manifolds that only capture heat losses (M-EGR and M-HOQ) significantly underestimate the flame consumption speed $s_{{c}}$ , as their flamelet databases do not include flamelets with varying equivalence ratio as discussed previously.

In summary, the manifolds that capture mixture variation (M-TD, M-EGR-TD and M-HOQ-TD) exhibit very high accuracy in the a-priori analysis, whereas the manifolds that account only for heat losses (M-EGR and M-HOQ) show significant errors for the TD unstable flame.

A-posteriori analysis: in the following, the manifolds are evaluated in an a-posteriori assessment based on the fully coupled TC simulations, in comparison with the DC reference simulation. This serves as a validation step to demonstrate that, in particular, the 3D manifolds can accurately represent TD instabilities a-posteriori, which is an essential prerequisite for subsequently investigating the 2D HOQ configuration where FWIs and TD instabilities interact.

Physical space

Figure 11. Snapshots of the flame front in the nonlinear regime showing the normalised temperature $\varTheta$ (top row) and elemental mass fraction of hydrogen $Z_{\textrm {H}}$ (bottom row) in a section of the computational domain for the DC simulation and the fully coupled TC simulations.

Figure 11 presents snapshots of a zoom into the flame fronts in the nonlinear regime, comparing the results of the fully coupled TC simulations with the DC reference. The top row shows profiles of the normalised temperature $\varTheta$ , defined as

(5.3) \begin{equation} \varTheta = \frac {T - T_{{u}}}{T_{\textit{ad}} - T_{{u}}}, \end{equation}

where $T_{\textit{ad}}$ is the adiabatic flame temperature of the unburned mixture, and the bottom row shows profiles of the elemental hydrogen mass fraction $Z_{\textrm {H}}$ .

The analysis of the profiles demonstrates that all manifolds accounting for mixture variations (M-TD, M-EGR-TD and M-HOQ-TD) are capable of accurately representing the structure and dynamics of TD instabilities. This includes local mixture variation in positively and negatively curved regions of the flame front, as well as the presence of superadiabatic temperatures ( $\varTheta \gt 1$ ). Furthermore, the manifolds reproduce both the small-scale features of the instabilities (i.e. small-scale wrinkling) and the large-scale structures (i.e. flame fingers). For manifold M-HOQ, the findings from the linear regime and the a-priori analysis of the nonlinear regime are confirmed: M-HOQ is not capable of predicting TD instabilities as only DL instabilities develop. Results for manifold M-EGR are not shown, as they are analogous to those of M-HOQ and similarly fail to reproduce TD instabilities.

These findings are further supported and quantified by comparisons of flame front curvature and instability cell size between the DC and TC simulations, as detailed in the supplementary material.

Thermochemical state space

Figure 12. Joint PDFs of the $\textrm {H}_{2}$ reaction rate $\dot {\omega }_{\textrm {H}_{2}}$ (fuel consumption rate) and the normalised reaction progress variable $Y_{{c,\textit{norm}}}$ (top row), and of the elemental mass fraction $Z_{\textrm {H}}$ and $Y_{{c,\textit{norm}}}$ (bottom row) for the DC simulation and the fully coupled TC simulations employing manifolds M-TD (centre) and M-HOQ (right). The respective conditional means $\langle \dot {\omega }_{\textrm {H}_{2}} \mid Y_{{c,\textit{norm}}} \rangle$ and $\langle Z_{\textrm {H}} \mid Y_{{c,\textit{norm}}} \rangle$ (red lines) are included for all cases, and the corresponding conditional means from the DC simulation, $\langle \dot {\omega }_{\textrm {H}_{2}} \mid Y_{{c,\textit{norm}}} \rangle _{\textit{DC}}$ and $\langle Z_{\textrm {H}} \mid Y_{{c,\textit{norm}}} \rangle _{\textit{DC}}$ , are also shown in the TC results (dashed grey lines). The respective profiles from an unstretched flamelet at the unburnt mixture conditions are included (black lines).

Following the analysis in physical space, the fully coupled TC simulations are compared with the DC simulation in thermochemical state space. Figure 12 presents the joint distributions of the fuel consumption rate $\dot {\omega }_{\textrm {H}_{2}}$ and the progress variable $Y_{{c,\textit{norm}}}$ (top row), as well as the elemental mass fraction $Z_{\textrm {H}}$ (mixture fraction; bottom row), for the DC case and for the TC simulations employing manifolds M-TD and M-HOQ. Results for manifolds M-HOQ-TD and M-EGR-TD are not displayed because their state spaces are nearly identical to that of M-TD, while those for manifold M-EGR are omitted as they closely resemble that of M-HOQ.

The large scatter in the fuel consumption rate, which correlates with fluctuations in the mixture fraction (Schneider et al. Reference Schneider, Nicolai, Schuh, Steinhausen and Hasse2025a ), as evident in figure 12 (bottom row) and figure 11 (bottom row), is captured accurately by manifold M-TD (and likewise by manifolds M-EGR-TD and M-HOQ-TD). In contrast, the absence of the mixture fraction variation in the flamelet database of manifold M-HOQ (and M-EGR) prevents reproduction of these fluctuations, which also explains the large discrepancies observed in physical space.

The conditional means $\langle \dot {\omega }_{\textrm {H}_{2}} \mid Y_{{c,\textit{norm}}} \rangle$ and $\langle Z_{\textrm {H}} \mid Y_{{c,\textit{norm}}} \rangle$ obtained from the TC simulations employing the manifolds including mixture variation agree closely with those from the DC simulation, and deviate markedly from the profiles of a FP flamelet evaluated at the unburnt mixture conditions. For the manifolds that only capture heat losses, the opposite holds: good agreement with the profiles from the FP flamelet but large deviations from the conditional means of the DC simulation. The thermochemical states that deviate from the conditional mean are caused by discrepancies in the negatively curved regions of the DL instability (see figure 11). These discrepancies are presumably associated with manifold states where enthalpy variation, captured in the flamelet database, compensates for mixture variations.

This demonstrates that the manifolds incorporating mixture variations in the database are capable of accurately predicting the instability. Consequently, the global flame consumption speed $s_{{c}}$ and the reactivity factor $I_0$ are also expected to be accurately captured a-posteriori by these manifolds. In contrast, the manifolds that only capture heat losses are expected to reproduce only the 1D laminar flame speed and, therefore, fail to account for the enhancement associated with TD instabilities. This is discussed in detail in the following.

Global flame properties

Figure 13. (a): Temporal averages of the normalised flame consumption speed $s_{{c}} / s_{{l}}$ , the normalised flame surface area $A / L_y$ and the reactivity factor $I_0$ (mean values $\pm 1$ standard deviation $\sigma$ ) for the DC reference simulation and the fully coupled TC simulations. (b) Relative errors of the temporally averaged normalised flame consumption speed $\varepsilon (s_{{c}} / s_{{l}})$ of the fully coupled TC simulation with respect to the DC reference simulation. (c) Relative errors of the temporally averaged reactivity factor $\varepsilon (I_0)$ of the fully coupled TC simulation with respect to the DC reference simulation.

Figure 13(a) shows temporal averages of the normalised flame consumption speed $s_{{c}} / s_{{l}}$ , the normalised flame surface area $A / L_y$ and the reactivity factor (mean values with $\pm 1$ standard deviation $\sigma$ ), for the DC simulation and the fully coupled TC simulations. The flame surface area is defined as the length of the isoline corresponding to $Y_{{c,\textit{norm}}} = 1 - Y_{\textrm {H}_{2}} / Y_{\textrm {H}_{2},{u}} = 0.9$ , following Howarth & Aspden (Reference Howarth and Aspden2022). Manifolds M-TD, M-EGR-TD and M-HOQ-TD accurately predict all three quantities, with relative errors in the global consumption speed $\varepsilon (s_{{c}}/s_{{l}})$ below $4\,\%$ (see figure 13 b). The relative errors of the reactivity factor $\varepsilon (I_0)$ are below $5\,\%$ (see figure 13 c), which is particularly relevant for the subsequent HOQ analysis (Schneider et al. Reference Schneider, Nicolai, Schuh, Steinhausen and Hasse2025a , Reference Schneider, Nicolai, Schuh, Steinhausen and Hasseb ). The manifolds that include mixture variations in the flamelet database are therefore capable of predicting the TD instability with high accuracy. As discussed in detail in the previous sections, the manifolds M-EGR and M-HOQ are, by contrast, unable to capture the TD instability and, consequently, fail to predict the enhanced global flame consumption speed $s_{{c}}/s_{{l}}$ , the increased flame surface area $A / L_y$ and the elevated local reactivity $I_0$ .

To summarise the findings from the 2D FP flame validation setup, the manifolds M-TD, M-EGR-TD and M-HOQ-TD accurately predict TD unstable flames in both the linear and nonlinear regimes, whereas the manifolds M-EGR and M-HOQ fail to capture the TD instability.

5.6.1. A-priori analysis

Figure 14 presents representative wall-normal profiles (at $y/\delta_{T,{l}} = 20$ ) of selected thermochemical quantities, analogous to the 1D HOQ analysis, for the DC simulation and a-priori lookups using manifolds M-EGR-TD and M-HOQ-TD at different time steps during the quenching process.

Figure 14. Profiles of temperature $T$ (first row), ${\textrm {H}_{2}}\textrm {O}$ mass fraction $Y_{{\textrm {H}_{2}}\textrm {O}}$ (second row), $\textrm {H}$ mass fraction $Y_{\textrm {H}}$ (third row), ${\textrm {H}_{2}}\textrm {O}_{2}$ mass fraction $Y_{\textrm {H}_{2}\textrm {O}_{2}}$ (fourth row) and heat release rate $\dot {\omega }^{\prime }_{T}$ (fifth row) along a wall-normal line ( $y/\delta _{T,{l}} = 20$ ) in the near-wall region, obtained from the 2D HOQ DC simulation and from a-priori lookups using manifolds M-EGR-TD and M-HOQ-TD, at different time steps $t$ .

Very similar results to those of the 1D HOQ case are observed. The mass fractions of major species such as ${\textrm {H}_{2}}\textrm {O}$ , as well as most other species including radicals like $\textrm {H}$ , can be predicted with high accuracy by both manifolds throughout the entire quenching process. For manifold M-EGR-TD, however, noticeable deviations appear in the near-wall region at later stages of the quenching process in the mass fraction profiles of species that are particularly sensitive to quenching dynamics, such as $Y_{\textrm {H}_{2}}{\textrm {O}_{2}}$ . In contrast, manifold M-HOQ-TD predicts these profiles with high accuracy as well. The heat release rate $\dot {\omega }^{\prime }_T$ profiles also show minor deviations when predicted using manifold M-EGR-TD. Since this quantity enters the transport equation for temperature $T$ in the fully coupled simulations, larger deviations can be expected in the a-posteriori analysis for the fully coupled TC simulations with manifold M-EGR-TD, compared with M-HOQ-TD, similar to what was observed in the 1D HOQ case. Nevertheless, both manifolds are capable of accurately predicting the individual effects of TD instabilities, differential diffusion, quenching (i.e. heat losses), as well as the interactions between these phenomena, in the a-priori assessment.

Note that these results also apply to other wall positions (i.e. different $y$ coordinates). An additional example, in which the flame front impinges on the wall at a different angle, is provided in the supplementary material.

5.6.2. A-posteriori analysis

Figure 15 shows the profiles of normalised temperature $\varTheta$ and elemental mass fraction $Z_{\textrm {H}}$ in the near-wall region at various time steps during the quenching process for the DC simulation and the corresponding TC simulations. The profiles show very good agreement for all time steps, especially also in the near-wall region where the flame quenches and a temperature boundary layer develops. Both manifolds are also able to reproduce the near-wall mixture shift, which was analysed and discussed in detail for the 1D HOQ case (see third row in figure 6). Furthermore, only minor differences between the two different TC simulations can be observed based on these profiles. This demonstrates that both manifolds also perform well a-posteriori for the 2D HOQ of a TD unstable flame. Thus, both manifolds capture the interaction between differential diffusion, the associated TD instabilities and FWIs. To further assess the performance of the two manifolds, the wall heat flux and quenching distance are analysed and compared across the simulation results in the following.

Figure 15. Snapshots of the normalised temperature $\varTheta$ (a) and the elemental mass fraction $Z_{\textrm {H}}$ (b) at different time steps of the 2D HOQ simulations are shown for the DC case (first row) and the fully coupled TC simulations using manifolds M-EGR-TD (second row) and M-HOQ-TD (third row).

Figure 16. Wall heat flux $\varPhi$ , normalised by the quenching wall heat flux of the 1D HOQ case (DC) $\varPhi _{{q,\textrm{1D}}}$ , shown for the same time steps as in figure 15 for the DC simulation (a) and the fully coupled TC simulations employing manifolds M-EGR-TD (b) and M-HOQ-TD (c).

Figure 16 shows the normalised wall heat flux $\varPhi / \varPhi _{{q,\textrm{1D}}}$ for the DC simulation and the fully coupled TC simulations employing manifolds M-EGR-TD and M-HOQ-TD. Note that, at the first time step, the flame has not yet reached the wall in any of the simulations, and the wall heat flux $\varPhi$ is therefore zero. The discussion first focuses on the analysis of the physical phenomena using the DC simulation, followed by an assessment of the reproduction of these phenomena in the TC simulations.

As seen in figure 15, a single ‘flame finger’ first approaches and quenches head-on at the lower wall, which is accompanied by an increase in the wall heat flux. At the edges of this region, the wall heat flux locally exceeds the quenching wall heat flux of the 1D HOQ $\varPhi _{{q,\textrm{1D}}}$ , which correlates with the locally increased mixture fraction (equivalence ratio). Subsequently, the flame propagates in both positive and negative $y$ directions, locally quenching in a SWQ manner. The wall heat flux is highest at the propagating edges, where the flame burns close to and quenches at the wall, whereas in regions where the flame has already extinguished, the hot products are cooled by the wall, resulting in a continued decrease in wall heat flux over time. At later times, additional parts of the flame, particularly in the upper region, also reach the wall and locally quench; in some areas still resembling SWQ, in others resembling HOQ, and in some regions, the flame quenches under varying local angles. Combined with the locally varying mixture fraction, this results in strong spatial and temporal fluctuations of the wall heat flux. Finally, once the flame is fully extinguished, the wall heat flux continues to decrease as the burned gases are further cooled by the wall. For a more detailed analysis of the 2D HOQ case based on DC simulations, the interested reader is referred to Schneider et al. (Reference Schneider, Nicolai, Schuh, Steinhausen and Hasse2025a , Reference Schneider, Nicolai, Schuh, Steinhausen and Hasseb , Reference Schneider, Rong, Steinhausen, Hasse and Nicolaic ).

Both manifolds are able to qualitatively predict the wall heat flux very well. The characteristic spatial and temporal evolution, namely the initial increase in the lower region due to the quenching of a ‘flame finger’, followed by propagation along the wall in a SWQ manner associated with elevated wall heat fluxes at the propagating flame edges, the subsequent cooling of the burned gases by the wall, and the quenching of local flame structures (such as additional ‘flame fingers’) in the upper region, is well captured by the representation of TD instability features in the fully coupled TC simulations. As a result, the spatio-temporal variability of the wall heat flux is accurately reproduced, including the elevated wall heat fluxes compared with the 1D HOQ case, which are partly associated with local mixture variations. Accordingly, both manifolds, M-EGR-TD and M-HOQ-TD, predict the wall heat fluxes within the correct range of values; however, given that the quenching process is inherently unsteady, a precise quantitative comparison based on individual time steps alone is neither meaningful nor appropriate. Therefore, in the following, a statistical comparison of the wall heat flux between the different simulations is presented.

Figure 17. Mean values $\pm$ one standard deviation $\sigma$ of the quenching wall heat flux $\varPhi _{{q}}$ (a) and quenching distance $x_{{q}}$ (b) obtained from 2D HOQ simulations using DC and TC, with manifolds M-EGR-TD and M-HOQ-TD.

Figure 17 shows the mean quenching wall heat flux $\overline {\varPhi }_{{q}}$ (panel a) and mean quenching distance $\overline {x}_{{q}}$ (panel b), each presented with $\pm$ one standard deviation ( $\sigma$ ), for the 2D HOQ simulations using DC and fully coupled TC simulations employing the M-EGR-TD and M-HOQ-TD manifolds. Compared with the 1D HOQ (dashed grey lines in figure 17), the mean wall heat flux increases due to TD instabilities in the DC simulations. Moreover, it exhibits noticeable scatter around the mean value as a result of locally varying flame structures (see Schneider et al. (Reference Schneider, Nicolai, Schuh, Steinhausen and Hasse2025a , Reference Schneider, Nicolai, Schuh, Steinhausen and Hasseb ) for further details). This behaviour is qualitatively captured in the TC simulations employing either manifold.

Quantitatively, both manifold M-HOQ-TD and M-EGR-TD predict the mean quenching wall heat flux with good agreement. The mean quenching distance is likewise well predicted by both manifolds, although it is slightly overestimated in the TC simulation with M-HOQ-TD and to a greater extent with M-EGR-TD. With manifold M-EGR-TD, the standard deviation of the quenching wall heat flux is slightly underestimated, while that of the quenching distance is overestimated. For M-HOQ-TD, the latter is also slightly overestimated, but to a lesser extent. Overall, both manifolds perform very well in terms of predicting the global quenching properties. To further assess the performance of the manifolds with respect to the temporal evolution of the quenching process, the wall heat flux profiles as a function of time are analysed in the following.

Figure 18. (a) Distribution of the wall heat flux $\varPhi$ over the relative time $t - t_{{q}}$ for the 2D HOQ configuration for the DC simulation and the TC simulations employing manifolds M-EGR-TD and M-HOQ-TD. (b) Conditional means of the wall heat flux $\langle \varPhi \,|(t-t_{{q}}) \rangle$ for the same simulations.

Figure 18(a) shows distributions of the wall heat flux $\varPhi$ over the relative time $(t - t_{{q}})$ , as well as the conditional means of the wall heat flux $\langle \varPhi \,|(t-t_{{q}}) \rangle$ for the DC and TC simulations (for manifolds M-EGR-TD and M-HOQ-TD). To facilitate comparison, the conditional means of the wall heat flux, $\langle \varPhi \mid (t - t_{{q}}) \rangle$ , are also shown in a combined plot in figure 18(b).

As illustrated in the preceding analysis of figure 17, it is evident that both TC simulations (employing manifolds M-EGR-TD and M-HOQ-TD) accurately capture the increase in wall heat flux, as well as the variations observed in the wall heat flux profiles, which are attributed to TD instabilities. However, the extreme values, i.e. the maximum quenching wall heat flux, are captured more accurately by manifold M-HOQ-TD than by manifold M-EGR-TD. In addition to the previously discussed deviations in the quenching wall heat fluxes at $(t - t_{{q}} = 0)$ , similar to the 1D HOQ case, the shape of the wall heat flux profile predicted by M-EGR-TD does not perfectly match that of the DC simulation. The wall heat flux profile is somewhat broader, and the post-quenching wall heat flux remains elevated relative to the 1D HOQ. This behaviour is not observed in the DC simulation and appears only to a smaller extent in the TC simulation employing manifold M-HOQ-TD, and is presumably linked to varying rates of heat loss at the wall (Efimov et al. Reference Efimov, de Goey and van Oijen2019). These observations are fully consistent with the results of the 1D HOQ analysis.

In summary, both manifolds, M-EGR-TD and M-HOQ-TD, demonstrate high accuracy in capturing the HOQ of a TD unstable hydrogen/air flame, with M-HOQ-TD exhibiting slightly superior performance compared with M-EGR-TD. From the analysis of the 1D HOQ and the a-priori analysis of the 2D HOQ, it becomes apparent that the local thermochemical state is better predicted in the near-wall region during quenching with regard to some minor species, such as ${\textrm {H}_{2}}{\textrm {O}_{2}}$ and $\textrm {HO}_{2}$ . However, the overall performance of M-EGR-TD with respect to global quantities such as the wall heat flux and also most species distributions, including major species like ${\textrm {H}_{2}}\textrm {O}$ and radical species like $\textrm {H}$ , is very good. In this context, it is also apparent that an accurate prediction of the reactivity factor $I_0$ leads to a reliable estimation of the mean wall heat flux. Consequently, the quenching of lean premixed hydrogen/air flames depends strongly on the local mixture composition, irrespective of the underlying physical mechanism (EGR or HOQ) by which heat losses are introduced in the flamelets used to construct the underlying manifold database. This behaviour contrasts with that observed for hydrocarbon fuels, where manifolds based on HOQ flamelets, compared with those derived from flamelets with varying enthalpy in the initial mixture (e.g. through EGR), perform significantly better in predicting pollutant formation (primarily CO), and also better in capturing global quantities of the quenching process, such as the wall heat flux. This can be explained as follows. When comparing the thermochemical state spaces of $\textrm {H}_{2}$ and $\textrm {CH}_{4}$ , the one corresponding to $\textrm {H}_{2}$ is less complex in the context of the quenching process. This is, for instance, reflected in the a-priori analysis of the 1D and 2D HOQ cases, where only minor deviations are observed in the species profiles of ${\textrm {H}_{2}}\textrm {O}_{2}$ and $\textrm {HO}_{2}$ when using manifold M-EGR-TD, compared with the DC simulations and manifold M-HOQ-TD. However, since these minor species have mass fractions several orders of magnitude smaller than those of the major species, their deviations have a negligible impact on the overall thermochemical state and, consequently, on thermal conductivity.

For $\textrm {CH}_{4}$ , on the other hand, the prediction of especially CO, but also other species, near the wall is known to be inaccurate (Ganter et al. Reference Ganter, Heinrich, Meier, Kuenne, Jainski, Rißmann, Dreizler and Janicka2017, Reference Ganter, Straßacker, Kuenne, Meier, Heinrich, Maas and Janicka2018). Because CO constitutes a significantly larger fraction of the total mass in $\textrm {CH}_{4}$ flames (Efimov et al. Reference Efimov, de Goey and van Oijen2019; Steinhausen et al. Reference Steinhausen2020, Reference Steinhausen, Zirwes, Ferraro, Scholtissek, Bockhorn and Hasse2023; Schneider et al. Reference Schneider, Steinhausen, Nicolai and Hasse2024), compared with species such as ${\textrm {H}_{2}}\textrm {O}_{2}$ and $\textrm {HO}_{2}$ , also in the near-wall region, deviations in its prediction exert a much stronger influence on the local thermochemical state and, consequently, on the thermal conductivity and, thus, on the wall heat flux. This effect is particularly relevant in fully coupled simulations, where the thermal conductivity feeds back into the transport equations of the controlling variables via the diffusion term. Additionally, it should also be noted that the HOQ manifolds were primarily developed to improve the prediction of CO. As a secondary effect, the wall heat flux is captured slightly more accurately with the HOQ-based manifold, although it is already predicted with good accuracy (showing only small deviations in the peak values) when using the EGR-based manifolds (Efimov et al. Reference Efimov, de Goey and van Oijen2019; Schneider et al. Reference Schneider, Steinhausen, Nicolai and Hasse2024). For further details, especially on the importance of accurately predicting enthalpy gradients at the wall to improve overall near-wall prediction accuracy, see Ganter et al. (Reference Ganter, Heinrich, Meier, Kuenne, Jainski, Rißmann, Dreizler and Janicka2017, Reference Ganter, Straßacker, Kuenne, Meier, Heinrich, Maas and Janicka2018), Efimov et al. (Reference Efimov, de Goey and van Oijen2019). Based on these observations, in the context of lean $\textrm {H}_{2}$ combustion, manifold M-EGR-TD is generally recommended for accurately predicting the quenching process in TD unstable flames, as it is easier to construct, being based on FP rather than HOQ flames. If an accurate prediction of extreme values or local thermochemical states of certain minor species is required, manifold M-HOQ-TD should be preferred. In this context, highly unsteady processes, such as turbulent flames or flashback (Gruber et al. Reference Gruber, Chen, Valiev and Law2012), require further investigation using both manifolds.