4.1. Optimal parameters for source term model

For the systematic development of any model and, in particular, a model for the progress variable source term $\dot {\omega }_{i}$, it is important to first identify the relevant input parameters, with which the quantity of interest can be best parametrized. This is achieved in an optimal estimator analysis (Moreau, Teytaud & Bertoglio Reference Moreau, Teytaud and Bertoglio2006; Berger et al. Reference Berger, Kleinheinz, Mueller and Pitsch2018), in which the capability of a set of parameters $\boldsymbol {\psi }$, e.g. $\boldsymbol {\psi } = [C_{H_2O},Z]$, to parametrize a target quantity $Q$, e.g. $Q=\dot {\omega }_{C_{H_2O}}$, is quantified by an error norm referred to as irreducible error. The amount of scatter of $Q$ with respect to the conditional mean $\langle Q \mid \boldsymbol {\psi } \rangle$ is measured by the quadratic error norm and the parametrization is good if small irreducible errors are observed. Further details on the concept of the optimal estimator and technical details are provided by Moreau et al. (Reference Moreau, Teytaud and Bertoglio2006) and Berger et al. (Reference Berger, Kleinheinz, Mueller and Pitsch2018). In the following, the irreducible errors are conditionally averaged with respect to progress variable $C_{i}$, defined as

(4.1)\begin{equation} \epsilon_{{irr}^*}^2 = \langle \left(\dot{\omega}_{i} - \langle \dot{\omega}_{i} \mid \boldsymbol{\psi} \rangle\right)^2 \mid C_{i} \rangle, \end{equation}

where $\dot {\omega }_{i}$ represents the source term for either $Y_{H_2O}$ or $Y_{H_2}$. The irreducible error $\epsilon _{{irr}^*}^2$ is normalized by the maximum value of the conditionally averaged source term $\langle \dot {\omega }_{i} \mid C_{i} \rangle$, yielding

(4.2)\begin{equation} \epsilon_{{irr}}^2 = \frac{\epsilon_{{irr}^*}^2}{\max(\langle \dot{\omega}_{i} \mid C_{i} \rangle)^2}. \end{equation}

With this normalization, a value of unity indicates irreducible errors that are as large as the maximum value of the average source term within the flame front. In the following, only the analysis based on the hydrogen-based progress variable $C_{i}=C_{H_2}$ is presented since the analysis concerning the water-based progress variable yields similar findings.

Figure 2 shows the irreducible errors for different parametrizations of the progress variable source term with different sets of parameters evaluated by means of the DNS data. As the development of an LES combustion model is pursued, three different filter sizes are shown: the unfiltered fields, and two filter sizes $\varDelta$ that correspond to one and two flame thicknesses $l_{F}$. For the filtering procedure, a box-filter is used.

Figure 2.

Irreducible errors of progress variable source term using different sets of parameters. Analysis for (a) unfiltered fields and (b,c) filtered fields with two different filter sizes $\varDelta$. Note that $\varDelta = l_{F}$ and $\varDelta = 2 l_{F}$ correspond to filtering of 10 and 20 grid points in the DNS, respectively.

For the unfiltered fields, irreducible errors are shown if parametrizing the progress variable source term only by progress variable, e.g. $\boldsymbol {\psi } = [C_{H_2}]$, or by progress variable and mixture fraction, e.g. $\boldsymbol {\psi } = [C_{H_2},Z]$. Using only progress variable for parametrization, the average fluctuations of the source term are as large as 70 % of the mean value at a progress variable of $C_{H_2}=0.7$, indicating strong variations, which cannot be described solely by progress variable. Consistent with the analyses of Berger et al. (Reference Berger, Attili and Pitsch2022a), including mixture fraction significantly improves the parametrization of the source term as it intrinsically captures the effects of differential diffusion. For instance, the maximum of the conditionally averaged irreducible errors in figure 2 drops from $\max (\epsilon _{irr}^2) = 0.66$ for the parametrization only by progress variable to a value of $\max (\epsilon _{irr}^2) = 0.012$ if additionally including mixture fraction. Thus, the instantaneous progress variable source term $\dot {\omega }_{H_2}$ may be parametrized as $\dot {\omega }_{H_2}(C_{H_2},Z)$. This also indicates that the assumption of a one-step global reaction that is used to defined the mixture fraction in (2.2) does not negatively impact the parametrization of the detailed chemical reactions by progress variable and mixture fraction. The parametrization of the source term by the different parameter sets is visualized in figure 3, where two scatters of the DNS data are shown. If the source term is parametrized by progress variable only, significant scattering can be seen, while in figure 3(b), where the source term is parametrized by progress variable and mixture fraction, a good parametrization is achieved as the data collapse onto a manifold.

Figure 3.

Parametrization of source term $\dot {\omega }_{H_2}$ if using (a) progress variable and (b) progress variable and mixture fraction to condition the DNS data. Additionally, the conditional mean $\langle \dot {\omega }_{H_2} \mid C_{H_2} \rangle$ is shown in panel (a).

For the filtered fields, additional input parameters need to be considered to parametrize the filtered progress variable source term $\bar {\dot {\omega }}_{H_2}$. While the instantaneous progress variable source term $\dot {\omega }_{H_2}$ only depends on progress variable and mixture fraction, the filtered progress variable source term $\bar {\dot {\omega }}_{H_2}$ is also affected by the subfilter distribution, which can be represented by the higher-order moments of progress variable and mixture fraction, such as $\widetilde {C_{H_2}}$, $\tilde {Z}$, $\widetilde {{C_{H_2}^{{''}^2}}}$ etc. This can be formally written as

(4.3)\begin{equation} \bar{\dot{\omega}}_{H_2}= \bar{\dot{\omega}}_{H_2}\left(\widetilde{C_{H_2}},\tilde{Z},\widetilde{{C_{H_2}^{{''}^2}}},\widetilde{{Z^{{''}^2}}}, \widetilde{C_{H_2} Z},\ldots\right) . \end{equation}

Note that $\widetilde {{Q^{{''}^2}}}$ denotes the subfilter variance of a quantity $Q$. Thus, for the filtered fields, irreducible errors of several more different parameter sets are investigated in figure 2(b,c) as combinations of different subfilter moments are assessed.

While progress variable and mixture fraction describe the unfiltered source term very well, it is evident that their parametrization capability is significantly reduced for increasing filter sizes. In particular, the irreducible error with $\boldsymbol {\psi } = [\widetilde {C_{H_2}},\tilde {Z}]$ is close to $\boldsymbol {\psi } = [\widetilde {C_{H_2}}]$ in figure 2(c). Also, using only progress variable and any of the second order moments does not reduce irreducible errors significantly. Thus, three parameters are needed and the best parametrization is achieved with the set $\boldsymbol {\psi } = [\widetilde {C_{H_2}},\tilde {Z},\widetilde {{C_{H_2}^{{''}^2}}}]$. However, adding a fourth parameter, e.g. a second variance, does not yield any significant further reduction and, hence, is not needed. This is also quantitatively assessed by the maximum irreducible errors of each parametrization, which are provided in Appendix B. It is worth noting that using the variance of progress variable $\widetilde {{C_{H_2}^{{''}^2}}}$ in addition to $\widetilde {C_{H_2}}$ and $\tilde {Z}$ yields lower irreducible errors than using the variance of mixture fraction $\widetilde {{Z^{{''}^2}}}$, indicating a higher dependence of the subfilter source term fluctuations on the subfilter fluctuations of progress variable than of mixture fraction. Thus, an LES combustion model for thermodiffusively unstable flames should be based on progress variable, mixture fraction and the variance of progress variable. Similar conclusions are drawn from an optimal estimator analysis for the water-based progress variable, showing that $\bar {\dot {\omega }}_{H_2O}$ is well parametrized by the parameter set $\boldsymbol {\psi } = [\widetilde {C_{H_2O}},\tilde {Z},\widetilde {{C_{H_2O}^{{''}^2}}}]$.

4.2. Formulation of source term model

Based on the findings of the optimal estimator analysis, a model for the progress variable source term is formulated in this section. A presumed p.d.f. approach is chosen and the filtered progress variable source term $\bar {\dot {\omega }}_{H_2}$ is expressed as

(4.4)\begin{equation} \bar{\dot{\omega}}_{H_2}= \bar{\rho} \int \frac{1}{\rho} \dot{\omega}_{H_2}(C_{H_2},Z) \tilde{\mathcal{P}}(C_{H_2},Z) \,\text{d}C_{H_2} \,\text{d}Z. \end{equation}

Here, $\tilde {\mathcal {P}}(C_{H_2},Z)$ is the density-weighted joint subfilter p.d.f. of progress variable and mixture fraction, and hence, to obtain the filtered source term $\bar {\dot {\omega }}_{H_2}$, the inverse density-weighting in (4.4) is performed. Favre-filtered quantities, such as the progress variable, are obtained as

(4.5) \begin{equation} \tilde{C}_{H_2} = \int C_{H_2} \tilde{\mathcal{P}}(C_{H_2},Z) \,\text{d}C_{H_2} \,\text{d}Z, \end{equation}

while the filtered density is given by

(4.6)\begin{equation} \frac{1}{\bar{\rho}} = \int \frac{1}{\rho} \tilde{\mathcal{P}}(C_{H_2},Z) \,\text{d}C_{H_2} \,\text{d}Z. \end{equation}

Similar to the previous section, the discussion is based on the modelling of $\bar {\dot {\omega }}_{H_2}$, while the model for the source term $\bar {\dot {\omega }}_{H_2O}$ of the water-based progress variable is discussed at the end of the section.

The following discussion of the model formulation is structured in three parts: first, a flamelet-based model to represent the unfiltered source term $\dot {\omega }_{H_2}(C_{H_2},Z)$ is discussed; second, a presumed p.d.f. model of $\tilde {\mathcal {P}}(C_{H_2},Z)$ is introduced and assessed by the DNS data; third, the model formulation based on a water-based progress variable is presented.

4.2.1. Modelling the unfiltered source term

First, a model for $\dot {\omega }_{H_2}(C_{H_2},Z)$ is required. For this, figure 4 shows the joint p.d.f. of progress variable and mixture fraction in the DNS, which displays the range of possible states in the turbulent flame, and compares different sets of representative flamelets with the turbulent flame states. In figure 4(a), strained flamelets with different strain rates $a$ and an equivalence ratio of $\phi =0.4$ are shown. These flamelets are generated in FlameMaster (Pitsch Reference Pitsch1998) by computing a counterflow flame, in which burned and unburned mixtures are injected towards each other. Both streams feature an equivalence ratio of $\phi =0.4$ and a pressure of $p=1$ bar, while the temperatures are set to $T_{u}=298\,{\rm K}$ in the unburned gas and to $T_{b}=1418\,{\rm K}$ in the burned gas, where the latter represents the adiabatic flame temperature. The outermost flamelets represent the unstretched flamelet (dotted black curve) and the conditions when extinction occurs (yellow curve), e.g. defined by a rapid reduction of reaction rates if strain rate is further increased. In figure 4(b), the joint p.d.f. is compared with a set of unstretched premixed flamelets with different equivalence ratios, ranging from $\phi =0.34$ to $\phi =1.0$, where the former is close to the flammability limit. It is evident that the strained flamelets can only parametrize a small portion of the entire flame state, while the set of unstretched flamelets covers almost the entire joint p.d.f. except for a small region of low mixture fraction and high progress variable values. However, extrapolation into this region does not pose any significant challenges as the source term of progress variable is close to zero due to the low mixture fraction values. Thus, the set of unstretched flamelets is used for modelling, as already proposed by Bastiaans et al. (Reference Bastiaans, Vreman and Pitsch2007) and Regele et al. (Reference Regele, Knudsen, Pitsch and Blanquart2013).

Figure 4.

Joint p.d.f. of progress variable and mixture fraction in the turbulent flame (a) overlaid by strained flamelets with different strain rates $a$ obtained from a premixed counterflow flame and (b) overlaid by unstretched flamelets with different equivalence ratios $\phi$.

Since the set of unstretched flamelets does not feature any overlaps, building a table to predict $\dot {\omega }_{H_2}(C_{H_2},Z)$ is straightforward. A comparison of the table prediction $\dot {\omega }_{H_2}^{FL}$ and the conditional average $\langle \dot {\omega }_{H_2} \mid C_{H_2},Z \rangle$ of the DNS data is shown in figure 5. While certain differences are seen in regions with high mixture fractions, which are characterized by high values of curvature, good qualitative agreement is observed overall and, in particular, at flame states of high probability, cf. figure 4. As will be shown quantitatively below in figure 12, low modelling errors are achieved with this model.

Figure 5.

(a) Conditionally averaged progress variable source term with respect to mixture fraction and progress variable and (b) model prediction of progress variable source by the set of unstretched premixed flamelets.

Analogously, a model based on the water-based progress variable can be built for the instantaneous source term $\dot {\omega }_{H_2O}$, yielding similarly good agreement with the DNS data and is discussed at the end of this section.

4.2.2. Modelling the subfilter p.d.f.

For LES, a model for the subfilter p.d.f. $\tilde {\mathcal {P}}(C_{H_2},Z)$ or $\tilde {\mathcal {P}}(C_{H_2O},Z)$ is required. The joint p.d.f. of the hydrogen-based progress variable and mixture fraction is shown in figure 4, while the p.d.f. with a water-based progress variable is shown in figure 6. Comparing the two p.d.f.s, a distinct difference is evident: in contrast to the hydrogen-based progress variable, the water-based progress variable features super-equilibrium values in the burned gas, e.g. $C_{H_2O} > 1$. In particular, $C_{H_2O}$ behaves similar to the temperature field, which features super-adiabatic values, cf. figure 1. Thus, $C_{H_2O}$ is not bounded like the hydrogen-based progress variable with $C_{H_2} \in [0,1]$. This makes a formulation of a p.d.f. for $C_{H_2O}$ challenging, so in the following, a presumed subfilter p.d.f. is only suggested for the hydrogen-based progress variable. However, a model for the water-based progress variable can be also formulated using the presumed subfilter p.d.f. for the hydrogen-based progress variable as discussed in the next section.

Figure 6.

Joint p.d.f. of progress variable $C_{H_2O}$ and mixture fraction $Z$ in the turbulent flame.

It is evident from figure 4 that progress variable and mixture fraction are not independent. This effect is also reflected by the strongly curved trajectories of the unstretched flamelets in figure 4(b), which feature a local minimum of mixture fraction for intermediate progress variable values. To transform mixture fraction and progress variable into a set of independent variables, which significantly simplifies the modelling of $\tilde {\mathcal {P}}(C_{H_2},Z)$, the joint p.d.f. of progress variable and a flamelet index is investigated. The latter is the value of the nominal equivalence ratio $\phi ^{FL}$ for each unstretched flamelet in figure 4(b) and, hence, is constant within a flamelet. Thus, in contrast to mixture fraction, $\phi ^{FL}$ is independent of progress variable. To determine $\phi ^{FL}$ locally from mixture fraction and progress variable, the flamelet index $\phi ^{FL}$ is tabulated as a function of these two parameters analogously to the source term in the previous section. Then, $\phi ^{FL}$ can be computed for each DNS data point and the joint p.d.f. of progress variable and the flame index is shown in figure 7. Conveniently, $\phi ^{FL}$ and $C_{H_2}$ are statistically almost independent. To highlight this aspect, figure 7 also shows the joint distribution of the local equivalence ratio $\phi$ of the DNS, which is determined as (Peters Reference Peters2000)

(4.7)\begin{equation} \phi = \frac{Z}{1-Z} \frac{1-Z_{st}}{Z_{st}} \end{equation}

and contains the same information as mixture fraction. Here, $Z_{st}$ indicates the stoichiometric value of mixture fraction, which is $Z_{st}=0.028$. Note that $\phi$ is different from the flamelet index $\phi ^{FL}$, as $\phi$ represents the local composition while $\phi ^{FL}$ indicates the nominal equivalence ratio in the unburned gas that corresponds to the local composition. For instance, from figure 4(b), it is evident that mixture fraction and the local composition within a flamelet vary due to differential diffusion effects, but $\phi ^{FL}$ is constant for a flamelet by definition. Comparing the two distributions $\mathcal {P}(C_{H_2},\phi ^{FL})$ and $\mathcal {P}(C_{H_2},\phi )$ in figure 7, a significantly stronger correlation between $\phi$ and progress variable $C_{H_2}$ is seen. This is further emphasized in figure 8, where the two joint distributions $\mathcal {P}(C_{H_2}=C_0,\phi ^{FL})$ and $\mathcal {P}(C_{H_2}=C_0,\phi )$ are shown for constant values of progress variable $C_0$. It is evident that the distribution $\mathcal {P}(C_{H_2},\phi )$ strongly varies with progress variable, e.g. a shifting of the peak is seen. In contrast, the distribution $\mathcal {P}(C_{H_2}, \phi ^{FL})$ is less sensitive to the value of progress variable, e.g. the peak and width of the distribution are similar for all progress variable values. Further, it is interesting to note that the conditional average of the flamelet index $\langle \phi ^{FL} \mid C_{H_2} \rangle$ lies within the flame above the nominal equivalence ratio of $\phi =0.4$. This indicates that, on average, combustion occurs on a richer flamelet of $\phi ^{FL} \approx 0.5$, which is consistent with the analysis of Berger et al. (Reference Berger, Attili and Pitsch2022b), who showed that combustion occurs at richer mixtures due to the mean positive strain rate universally observed in turbulent flows (Rutland & Trouvé Reference Rutland and Trouvé1993; Luca et al. Reference Luca, Al-Khateeb, Attili and Bisetti2018; Chu et al. Reference Chu, Berger, Gauding, Attili and Pitsch2024).

Figure 7.

Comparison of (a) joint distribution of progress variable and local equivalence ratio and (b) the joint distribution of progress variable and the flamelet index.

Figure 8.

(a) Joint distribution of progress variable and local equivalence ratio and (b) joint distribution of progress variable and the flamelet index evaluated at constant values of progress variable.

Assuming statistical independence of $\phi ^{FL}$ and $C_{H_2}$, and similar to the model of Knudsen et al. (Reference Knudsen, Kolla, Hawkes and Pitsch2013), a model for the presumed p.d.f. is proposed,

(4.8)\begin{equation} \tilde{\mathcal{P}}\left(C_{H_2},\phi^{FL}\right) = \beta(C_{H_2}) \delta\left(\phi^{FL} - \tilde{\phi}^{FL}\right). \end{equation}

Similar to previous works (Galpin et al. Reference Galpin, Naudin, Vervisch, Angelberger, Colin and Domingo2008; Vreman et al. Reference Vreman, van Oijen, de Goey and Bastiaans2009; Pfitzner Reference Pfitzner2021), a $\beta$-p.d.f. is assumed for $C_{H_2}$ and a $\delta$-p.d.f. is assumed for the flamelet index $\phi ^{FL}$. The $\delta$-p.d.f. implies that the combustion process locally occurs on an unstretched flamelet with an equivalence ratio of $\tilde {\phi }^{FL}$, where $\tilde {\phi }^{FL}=\tilde {\phi }^{FL}(\tilde {C}_{H_2}, \widetilde {{C_{H_2}^{{''}^2}}}, \tilde {Z})$ varies with the local values of the filtered progress variable, the progress variable variance and mixture fraction. Further note that the assumption of an $\beta$-p.d.f. for $C_{H_2}$ is possible, as $C_{H_2}$ is strictly bound to $C_{H_2} \in [0,1]$. For the water-based progress variable, a $\beta$-p.d.f. cannot be assumed due to the super-unity values of $C_{H_2O}$. Note that if using the formulation of a hydrogen-based progress variable and mixture fraction, the information of the super-equilibrium regions in the post-flame region, where $C_{H_2}=1$, is provided by mixture fraction.

Using (4.8), the filtered source term can be expressed as

(4.9)\begin{align} \bar{\dot{\omega}}_{H_2} &= \bar{\rho} \int \frac{\dot{\omega}_{H_2}(C_{H_2},\phi^{FL})}{\rho(C_{H_2},\phi^{FL})} \beta(C_{H_2}) \delta(\phi^{FL} - \tilde{\phi}^{FL}) \,\text{d}C_{H_2} \,\text{d}\phi^{FL}, \end{align}

(4.10)\begin{align} &= \bar{\rho} \int \frac{\dot{\omega}_{H_2}(C_{H_2},\tilde{\phi}^{FL})}{\rho(C_{H_2},\tilde{\phi}^{FL})} \beta(C_{H_2}) \,\text{d}C_{H_2} . \end{align}

Note that this integral depends on $\widetilde {C_{H_2}}$, $\widetilde {{C_{H_2}^{{''}^2}}}$ and $\tilde {Z}$ since the $\beta$-p.d.f. is parametrized by $\widetilde {C_{H_2}}$ and $\widetilde {{C_{H_2}^{{''}^2}}}$, and $\tilde {\phi }^{FL}$ is implicitly determined by $\tilde {Z}$. The latter is given by

(4.11)\begin{equation} \tilde{Z} = \int Z(C_{H_2},\phi^{FL}) \beta(C_{H_2}) \delta(\phi^{FL} - \tilde{\phi}^{FL}) \,\text{d}C_{H_2} \,\text{d}\phi^{FL}. \end{equation}

Thus, to determine $\tilde {\phi }^{FL}$ for the computation of $\bar {\dot {\omega }}_{H_2}$ in (4.10), the integral for the filtered mixture fraction in (4.11) is computed for all different flamelets, and the value $\tilde {\phi }^{FL}$ is given by the flamelet that matches the value of $\tilde {Z}$ at a given value of $\widetilde {C_{H_2}}$ and $\widetilde {{C_{H_2}^{{''}^2}}}$.

In practice, a table can be built directly as a function of $\widetilde {C_{H_2}}$, $\tilde {Z}$ and $\widetilde {{C_{H_2}^{{''}^2}}}$. For this, (4.10) and (4.11) are evaluated for each flamelet for a grid of $\widetilde {C_{H_2}}$ and $\widetilde {{C_{H_2}^{{''}^2}}}$ values. Here, a linear spacing is used for $\widetilde {C_{H_2}}$ and a quadratic spacing for $\widetilde {{C_{H_2}^{{''}^2}}}$ to enhance the table resolution for small values. Since $\tilde {Z}$ monotonically increases with $\tilde {\phi }^{FL}$, the data can be interpolated from the ($\widetilde {C_{H_2}}$, $\widetilde {{C_{H_2}^{{''}^2}}}$, $\tilde {\phi }^{FL}$)-space to the ($\widetilde {C_{H_2}}$, $\widetilde {{C_{H_2}^{{''}^2}}}$, $\tilde {Z}$)-space.

Figure 9 shows three panels of the tabulated value of $\bar {\dot {\omega }}_{H_2}$, where either $\widetilde {{C_{H_2}^{{''}^2}}}$ or $\tilde {Z}$ are kept constant. In figure 9(a), the table is evaluated at zero progress variable variance and, hence, corresponds to the tabulated unfiltered source term. In figure 9(b), the same plot is shown at $\widetilde {{C_{H_2}^{{''}^2}}}=0.1$. As expected, the source term magnitude reduces with increasing progress variable variance if comparing with panel (a). The same trend can be seen in figure 9(c), where the variation of the source term with progress variable variance is explicitly shown for a fixed value of mixture fraction. Further note that in an LES, the relevant range of the table is given by the constraint $\widetilde {{C_{H_2}^{{''}^2}}}<\tilde {C}_{H_2}(1-\tilde {C}_{H_2})$. This represents the maximum possible progress variable variance. However, for numerical stability, the table is extrapolated to the entire parametric space, e.g. by setting the value of any lookup scalar in the extrapolated region to the value at $\widetilde {{C_{H_2}^{{''}^2}}}=\tilde {C}_{H_2}(1-\tilde {C}_{H_2})$. Similarly, the table is extrapolated in the direction of mixture fraction, the physical region of which is bound towards lower values by the flamelet solution with the lowest equivalence ratio.

Figure 9.

Lookup table prediction of $\bar {\dot {\omega }}_{H_2}$ as function of the filtered parameters ($\widetilde {C_{H_2}}$, $\widetilde {{C_{H_2}^{{''}^2}}}$, $\tilde {Z}$) for (a,b) two particular values of $\widetilde {{C_{H_2}^{{''}^2}}}$ and (c) one for a constant value of $\tilde {Z}$.

4.2.3. Lookup-table based on water-based progress variable

For modelling unclosed filtered terms, a presumed subfilter p.d.f. has been introduced in previous sections. The subfilter p.d.f. was formulated only for the hydrogen-based progress variable as an analogous formulation is not possible for a water-based progress variable due to the super-equilibrium values of $C_{H_2O}$, which cannot be represented by a $\beta$-p.d.f. for $C_{H_2O}$. Thus, to model unclosed filtered terms with a water-based progress variable, a different approach is followed. In particular, no presumed p.d.f. is explicitly formulated for $C_{H_2O}$, but the previously introduced, well-posed presumed p.d.f. for $C_{H_2}$ is used. First, the values of $\widetilde {C_{H_2O}}$ and $\widetilde {{C_{H_2O}^{{''}^2}}}$ are determined with the $C_{H_2}$-based model. These values are computed for each unstretched flamelet $\phi ^{FL}$ on a pre-specified grid of $\widetilde {C_{H_2}}$ and $\widetilde {{C_{H_2}^{{''}^2}}}$ as

(4.12)$$\begin{gather} \widetilde{C_{H_2O}} = \int C_{H_2O}(C_{H_2},\phi^{FL}) \beta(C_{H_2}) \delta(\phi^{FL} - \tilde{\phi}^{FL}) \,\text{d}C_{H_2} \,\text{d}\phi^{FL}, \end{gather}$$

(4.13)$$\begin{gather}\widetilde{C_{H_2O}^2} = \int C_{H_2O}(C_{H_2},\phi^{FL})^2 \beta(C_{H_2}) \delta(\phi^{FL} - \tilde{\phi}^{FL}) \,\text{d}C_{H_2} \,\text{d}\phi^{FL}, \end{gather}$$

where $\widetilde {{C_{H_2O}^{{''}^2}}}$ is determined from $\widetilde {{C_{H_2O}^{{''}^2}}}=\widetilde {{C_{H_2O}}^2}-\widetilde {C_{H_2O}}^2$. It is worth noting that in (4.12) and (4.13), the instantaneous values of $C_{H_2O}$ are determined from the unstretched flamelet manifold that is parametrized by $C_{H_2}$ and $\phi ^{FL}$. Thus, $C_{H_2O}$ is looked up analogously to other quantities, such as $\rho$ or $\dot {\omega }_{H_2}$, as a function of $C_{H_2}$. Similarly, the filtered progress variable source term $\bar {\dot {\omega }}_{H_2O}$ of $C_{H_2O}$ is computed for each flamelet as

(4.14) \begin{equation} \bar{\dot{\omega}}_{H_2O} = \bar{\rho} \int \frac{\dot{\omega}_{H_2O}(C_{H_2},\phi^{FL}) }{\rho(C_{H_2},\phi^{FL})} \beta(C_{H_2}) \delta(\phi^{FL} - \tilde{\phi}^{FL}) \,\text{d}C_{H_2} \,\text{d}\phi^{FL}. \end{equation}

Since the flame state can be uniquely mapped from the ($\widetilde {C_{H_2}}$, $\widetilde {{C_{H_2}^{{''}^2}}}$)-space to the ($\widetilde {C_{H_2O}}$, $\widetilde {{C_{H_2O}^{{''}^2}}}$)-space for each flamelet, the filtered source term can be unambiguously expressed as a function of $\widetilde {C_{H_2O}}$ and $\widetilde {{C_{H_2O}^{{''}^2}}}$. For each flamelet, this relationship can be determined by computing the conditional mean

(4.15)\begin{equation} \bar{\dot{\omega}}_{H_2O}\left(\widetilde{C_{H_2O}}, \widetilde{{C_{H_2O}^{{''}^2}}},\tilde{\phi}^{FL}\right) = \left\langle \bar{\dot{\omega}}_{H_2O} \mid \widetilde{C_{H_2O}}, \widetilde{{C_{H_2O}^{{''}^2}}} \right\rangle_{\tilde{\phi}^{FL}} . \end{equation}

The notation $\langle \cdots \rangle _{\tilde {\phi }^{FL}}$ denotes that the mapping procedure from $[\widetilde {C_{H_2}}, \widetilde {{C_{H_2}^{{''}^2}}}]$ to $[\widetilde {C_{H_2O}}, \widetilde {{C_{H_2O}^{{''}^2}}}]$ is performed for each flamelet $\tilde {\phi }^{FL}$ separately. Figure 10 shows the mapping for the particular flamelet with $\tilde {\phi }^{FL}=0.4$. The visible data points are linked to the pre-specified discretization of the table based on $\widetilde {C_{H_2}}$ and $\widetilde {{C_{H_2}^{{''}^2}}}$. Due to the unique mapping, the scatter becomes a well-defined manifold and an unambiguous determination of $\bar {\dot {\omega }}_{H_2O}$ as a function of $\widetilde {C_{H_2O}}$ and $\widetilde {{C_{H_2O}^{{''}^2}}}$ is possible.

Figure 10.

Mapping from the hydrogen-based to water-based combustion model formulation for the flamelet with $\tilde{\phi}^{FL}=0.4$. Scatter of filtered source term $\bar {\dot {\omega }}_{H_2O}$ from table lookup as a function of filtered progress variable $\widetilde {C_{H_2O}}$ and progress variable variance $\widetilde {{C_{H_2O}^{{''}^2}}}$.

Finally, the data are mapped from the ($\widetilde {C_{H_2O}}$, $\widetilde {{C_{H_2O}^{{''}^2}}}$, $\tilde {\phi }^{FL}$)-space to the ($\widetilde {C_{H_2O}}$, $\widetilde {{C_{H_2O}^{{''}^2}}}$, $\tilde {Z}$)-space, where $\tilde {Z}$ is computed according to (4.11). The resulting table is shown in figure 11 for two different values of $\widetilde {{C_{H_2O}^{{''}^2}}}$ and at constant mixture fraction. The characteristic super-equilibrium regions ($\widetilde {C_{H_2O}}>1$) are well recovered with this approach. Additionally, it is worth noting that for large mixture fraction values featuring local equivalence ratios with $\phi ^{FL}>0.4$, the hydrogen mass fraction in the unburned state, defined by $\widetilde {C_{H_2O}}=0$, is larger than the hydrogen mass fraction $Y_{H_2,u}$ in the unburned gas in the DNS. Such regions cannot be parametrized with a hydrogen-based progress variable $C_{H_2} \in [0,1]$ and can also not occur in the DNS, so table entries in these regions represent only an extrapolation for numerical stability.

Figure 11.

Look-up table prediction of $\bar {\dot {\omega }}_{H_2O}$ as function of the filtered input parameters ($\widetilde {C_{H_2O}}$, $\widetilde {{C_{H_2O}^{{''}^2}}}$, $\tilde {Z}$) for (a,b) two particular values of $\widetilde {{C_{H_2O}^{{''}^2}}}$ and (c) one for a constant value of $\tilde {Z}$.

4.3. Assessment of source term model

Figure 12 shows the assessment of four different combustion models for three different filter sizes. The following models are considered for the comparison:

1. (i) model $\phi =0.4$-Fl. uses a single flamelet at the nominal equivalence ratio, neglecting the effects of a subfilter p.d.f. and the fluctuations of mixture fraction;

2. (ii) model $\phi =0.4$-Fl.+p.d.f. uses a single flamelet at the nominal equivalence ratio and while it neglects the fluctuations of mixture fraction, it considers the effects of a subfilter distribution by a presumed $\beta$-p.d.f. for progress variable. This is a common model for turbulent premixed flames of conventional fuels when thermodiffusive instabilities are not present (Pfitzner Reference Pfitzner2021);

3. (iii) model $\phi$-Var. describes the local mixture fraction fluctuations by a series of unstretched premixed flamelets with varying equivalence ratio according to Bastiaans et al. (Reference Bastiaans, Vreman and Pitsch2007) and Regele et al. (Reference Regele, Knudsen, Pitsch and Blanquart2013). However, the effects of the subfilter distribution are neglected; and

4. (iv) model $\phi$-Var.+p.d.f. is the model described in the previous section and additionally considers the presumed p.d.f. of (4.8) to account for the complex subfilter fluctuations of mixture fraction and progress variable.

Figure 12.

Model assessment for $\bar {\dot {\omega }}_{H_2}$ using a hydrogen-based progress variable. Modelling errors of the progress variable source term for three different filter sizes (a) $\varDelta = 0$, (b) $\varDelta =l_{F}$ and (c) $\varDelta =2 l_{F}$ for the four different flamelet-based models.

Figure 12 shows the conditionally averaged modelling error as a function of progress variable. The modelling error is defined as

(4.16)\begin{equation} \epsilon^2_{{mod.}} = \frac{ (\dot{\omega}_{H_2} - \dot{\omega}_{H_2}^{mod})^2}{\max \left[\langle \dot{\omega}_{H_2} \mid C \rangle^2\right]}, \end{equation}

where the same normalization as for the irreducible errors is chosen. This definition compares the absolute modelling errors with the maximum average progress variable source term within the flame. Note that the value of the normalization is constant for one particular filter size but decreases for increasing filter size, as filtering reduces the maximum value of the progress variable source term.

For the unfiltered DNS data, a significant reduction of modelling errors is observed by the flamelet manifold model $\phi$-Var. with respect to the single-flamelet model $\phi =0.4$-Fl., which neglects the effects of differential diffusion and yields an average deviation larger than 100 % of the filtered source term value within the DNS. In contrast, the modelling errors for the $\phi$-Var. model are almost zero, indicating an accurate representation of the flame state by the series of unstretched premixed flamelets. In particular, the maximum value of the conditional mean modelling errors in figure 12 drops from 1.53 for the $\phi =0.4$-Fl. model to 0.032 for the $\phi$-Var. model. For the filtered fields, the single-flamelet model $\phi =0.4$-Fl. again shows high modelling errors, which does not significantly change for the $\phi =0.4$-Fl.+p.d.f. model that additionally considers a subfilter distribution. These two models entirely neglect the effects of local mixture fraction fluctuations, so their failure is expected. It is worth noting that a value of unity indicates a significant error, as the modelling errors are as large as the maximum value of the conditionally averaged progress variable source term. Further, the modelling errors of model $\phi$-Var. are seen to increase with increasing filter size. This results from a strong overprediction of the progress variable source term by the model, as it does not consider the effect of filtering. In particular, filtering leads to a reduction of the peak reaction rates. In contrast, the newly proposed model $\phi$-Var.+p.d.f., which accounts for the effects of a subfilter distribution by means of the presumed p.d.f. in (4.8), yields a significant reduction of modelling errors, leading to maximum modelling errors below 0.11 for both filter sizes, as shown in Appendix B.

For the source term model based on the water-based progress variable, as discussed in § 4.2.3, figure 13 shows the a priori modelling errors for $\dot {\omega }_{H_2O}$ for three different filter sizes. Large errors are observed for the models $\phi =0.4$-Fl. and $\phi =0.4$-Fl.+p.d.f. as differential diffusion effects are not considered. A significant reduction of modelling errors is only observed if using the flamelet-manifold model in conjunction with the presumed subfilter p.d.f. However, comparing with figure 12, the hydrogen-based formulation of the flamelet-manifold model with subfilter p.d.f. appears to yield lower a priori modelling errors than the water-based model variant.

Figure 13.

Model assessment for $\bar {\dot {\omega }}_{H_2O}$ using a water-based progress variable. Modelling errors of the progress variable source term for three different filter sizes (a) $\varDelta = 0$, (b) $\varDelta =l_{F}$ and (c) $\varDelta =2 l_{F}$ for the four different flamelet-based models.

To highlight the good model prediction of the $\phi =0.4$-Fl.+p.d.f. model for flames, where differential diffusion effects are not present, figure 14 shows the modelling errors of the $\phi =0.4$-Fl. and $\phi =0.4$-Fl.+p.d.f. models if evaluated for the DNS case of Berger et al. (Reference Berger, Attili and Pitsch2022b), where differential diffusion effects have been artificially suppressed. For the unfiltered field, the single-flamelet model yields a good representation of the DNS as there is only moderate impact of turbulence on the flame, e.g. due to turbulent strain rate and curvature, cf. also the discussion by Berger et al. (Reference Berger, Attili and Pitsch2022b). It is also evident that for an increasing filter size, the consideration of a $\beta$-distribution for the subfilter p.d.f. is a good assumption that leads to negligibly small modelling errors. This analysis highlights that the failure of the $\phi =0.4$-Fl.+p.d.f. model in figure 12 for the DNS case, which features strong differential diffusion effects, is specifically linked to the presence of thermodiffusive effects and differential diffusion.

Figure 14.

Model assessment for $\bar {\dot {\omega }}_{H_2O}$ using a water-based progress variable for the TurbStable DNS that does not feature differential diffusion effects. Only the $\phi =0.4$-Fl. and $\phi =0.4$-Fl.+p.d.f. models are considered.

4.4. Model for progress variable variance

The lookup tables require the variance of progress variable as an input. Hence, an algebraic correlation or transport equation for $\widetilde {{C_{H_2}^{{''}^2}}}$ or $\widetilde {{C_{H_2O}^{{''}^2}}}$ is needed. As will be shown in detail, $\widetilde {{C_{i}^{{''}^2}}}$ is found to be well parametrized by progress variable $\widetilde {C_{i}}$, mixture fraction $\tilde {Z}$ and the magnitude of the gradient of progress variable $| \boldsymbol {\nabla } \widetilde {C_{i}} |$. Thus, in this work, a data-driven model based on the DNS data is employed for $\widetilde {{C_{i}^{{''}^2}}}$ to allow for a rigorous analysis of the LES as errors linked to the modelling of $\widetilde {{C_{i}^{{''}^2}}}$ can be ruled out. However, note that such a model may not be applicable to other conditions and flame configurations, as it is specific to the present DNS case. In the following, the modelling of $\widetilde {{C_{i}^{{''}^2}}}$ based on the DNS data is discussed.

For an a priori analysis of filtered quantities, first, an appropriate filter size needs to be chosen. However, this is not known a priori and may even vary within an LES due to different local conditions. For instance, filtering the DNS data with large filter sizes leads to low gradients of the filtered fields, so the local gradient information in the LES may be regarded as a measure of the local filter size. Since a model for $\widetilde {{C_{i}^{{''}^2}}}$ should be valid at different filter sizes and the filter size is not known a priori, the DNS data are filtered at several different filter sizes, ranging from $\varDelta =3\delta _{DNS}$ to $\varDelta =40\delta _{DNS}$, where $\delta _{DNS}=70~{\rm \mu} {\rm m}$ is the spatial discretization of the DNS data. To avoid limitations arising from a too narrow range of considered filter sizes, the range of filter sizes is chosen sufficiently large. The data are then used to perform an optimal estimator analysis to identify suitable parameters for modelling $\widetilde {{C_{i}^{{''}^2}}}$ across different filter sizes. For the optimal estimator analysis, the following variables are considered: the filtered progress variable and mixture fraction, and the gradient of the filtered progress variable. The former two parameters indicate the local flame state, while the gradient of the filtered progress variable is expected to parametrize the effect of different filter sizes.

Figure 15 shows the irreducible errors for the parametrization of $\widetilde {{C_{H_2}^{{''}^2}}}$ with different sets of parameters. The irreducible errors are defined as

(4.17)\begin{equation} \epsilon^2_{{irr}} = \frac{ \left(\widetilde{{C_{H_2}^{{''}^2}}} - \left\langle \widetilde{{C_{H_2}^{{''}^2}}} \mid \boldsymbol{\psi} \right\rangle\right)^2}{\max\left[\left\langle \widetilde{{C_{H_2}^{{''}^2}}} \mid C_{H_2} \right\rangle\right]}, \end{equation}

where $\boldsymbol {\psi }$ represents the set of parameters considered for the optimal estimator analysis. The irreducible errors are normalized with the maximum value of the conditional mean variance $\langle \widetilde {{C_{H_2}^{{''}^2}}} \mid C_{H_2} \rangle$ in the DNS. Note that in this analysis, conditional averaging, denoted by $\langle \rangle$, is performed for the ensemble of data that includes all considered filter sizes. High irreducible errors are observed if only using progress variable, while a significant reduction is achieved if using all three parameters. It is worth noting that if performing the optimal estimator analysis only for one particular filter size, a very good parametrization is already achieved using $\widetilde {C_{H_2}}$ and $\tilde {Z}$. However, the filter size is not known a priori and figure 15 indicates that if using data with different filter sizes, $| \boldsymbol {\nabla } \tilde {C}_{H_2}|$ is also needed for an adequate parametrization. To visualize the good parametrization of $\widetilde {{C_{H_2}^{{''}^2}}}$ with $\widetilde {C_{H_2}}$, $\tilde {Z}$ and $| \boldsymbol {\nabla } \widetilde {C_{H_2}}|$, figure 16 shows the scatter of $\widetilde {{C_{H_2}^{{''}^2}}}$ with respect to $\widetilde {C_{H_2}}$ and $| \boldsymbol {\nabla } \widetilde {C_{H_2}}|$ at a given value of $\tilde {Z}$ for filtered data at various filter sizes. It is seen that the scatter aligns on a low-dimensional manifold. Thus, in the LES, the progress variable variance $\widetilde {{C_{H_2}^{{''}^2}}}$ is modelled by the following conditional average based on the DNS data:

(4.18)\begin{equation} \widetilde{{C_{H_2}^{{''}^2}}}_{model} = \left\langle \widetilde{{C_{H_2}^{{''}^2}}} \mid \widetilde{C_{H_2}}, \tilde{Z}, | \boldsymbol{\nabla} \widetilde{C_{H_2}}| \right\rangle_{DNS}. \end{equation}

Figure 15.

Parametrization of the progress variable variance considering filtered DNS data with different filter sizes: irreducible errors for four different sets of parameters.

Figure 16.

Parametrization of the progress variable variance: scattered filtered DNS data at different filter sizes for a given value of mixture fraction.

For the water-based progress variable, an analogous model is formulated and a similar good parametrization of $\widetilde {{C_{H_2O}^{{''}^2}}}$ by $\widetilde {C_{H_2O}}$, $\tilde {Z}$ and $| \boldsymbol {\nabla } \widetilde {C_{H_2O}}|$ is found.

5.1. LES configuration

The computational domain of the DNS of Berger et al. (Reference Berger, Attili and Pitsch2022b) is considered for the LES and a detailed description can be found from Berger et al. (Reference Berger, Attili and Pitsch2022b). The domain is periodic in the spanwise direction ($z$), open boundary conditions are prescribed at the outlet in streamwise direction ($x$) and slip conditions are imposed at the boundaries in cross-wise direction ($y$). The inlet velocities of the central jet are taken from an auxiliary fully developed turbulent channel flow simulation and the same boundary conditions as in the DNS are applied (e.g. no filtered turbulent channel is computed). Laminar coflows are applied outside of the central jet, featuring a uniform velocity, which is set to 15 % of the bulk velocity $U$ of the central jet. An overview of the simulation parameters of the DNS of Berger et al. (Reference Berger, Attili and Pitsch2022b) and the LES of this work can be found in tables 1 and 2, respectively. Note that the LES domains are larger in axial direction compared with the DNS as longer flames are obtained in the LES.

Table 2.

Simulation parameters of the DNS, which corresponds to the TurbUnstable case of Berger et al. (Reference Berger, Attili and Pitsch2022b), and the different LES of this work: slot width $H$ and bulk velocity $U$ of central jet inflow; velocity $U_{Coflow}$ of laminar coflow; grid resolution $\varDelta$; domain size $L_{i}$; number of grid points $N_{i}$ in each direction; and jet Reynolds number $Re$.

5.2. Governing equations and numerical methods

The flow is modelled by the filtered Navier–Stokes equations in the low-Mach-number limit (Tomboulides, Lee & Orszag Reference Tomboulides, Lee and Orszag1997) and the fluid is assumed to be an ideal gas. The subfilter-stress term in the momentum equation is closed by a dynamic Smagorinsky model (Pierce & Moin Reference Pierce and Moin1998) that uses Lagrangian averaging (Meneveau, Lund & Cabot Reference Meneveau, Lund and Cabot1996). Due to the low-Mach-number formulation, a Poisson equation is solved for the pressure that enforces mass conservation. Further details can be found from Knudsen et al. (Reference Knudsen, Kolla, Hawkes and Pitsch2013). For a rigorous assessment of the combustion model, four different formulations are investigated a posteriori in this work, which are referred to as C-model, C-Z-model(H $_2$O), C-Z-model(H $_2$) and Finitechem LES; an overview of the models is given in table 3.

Table 3.

Comparison of LES combustion models: parametrization of local flame state and subfilter closure. Note that the presumed p.d.f. for the two models C-Z-model(H $_2$) and C-Z-model(H $_2$O) is identical, cf. (4.10) and (4.14). In the Finitechem LES model, transport equations for all species of the detailed chemical reaction mechanism of Burke et al. (Reference Burke, Chaos, Ju, Dryer and Klippenstein2012) are solved.

The C-model only solves the filtered transport equations of the water-based progress variable. Rather than using a manifold of flamelets with different equivalence ratio, the lookup table is created from one single unstretched premixed flamelet with a nominal equivalence ratio of $\phi =0.4$, which corresponds to the value of unburned mixture in the DNS. The C-Z-model(H $_2$O) and C-Z-model(H $_2$) solve the filtered transport equations of progress variable and mixture fraction in (2.13) and (2.14), where the former uses the formulation with a water-based progress variable and the latter employs a hydrogen-based progress variable. These two models represent the previously introduced new models that represent the combustion process with a manifold of unstretched flamelets of different equivalence ratio and model the subfilter fluctuations via a presumed p.d.f. In all of these three models, the fluid properties, such as density, viscosity, diffusion coefficients for molecular diffusion and the Soret effect, and the progress variable source term are computed from the lookup tables. The temperature field, needed to compute the terms related to the Soret effect, is also extracted from the same lookup tables. Finally, the Finitechem LES model uses a detailed chemical reaction mechanism (Burke et al. Reference Burke, Chaos, Ju, Dryer and Klippenstein2012), the same used in the DNS, and solves transport equations for all species. Details on this model are described by Berger et al. (Reference Berger, Attili and Pitsch2022b). This model does not consider any subfilter closure on purpose and simply represents an under resolved DNS. Thus, it allows to assess the relevance of subfilter modelling for the present flames. Analogously to Knudsen et al. (Reference Knudsen, Kolla, Hawkes and Pitsch2013), the subfilter scalar flux in all models is closed using a turbulent diffusivity computed from a dynamic Smagorinsky-type model (Pierce & Moin Reference Pierce and Moin1998).

A semi-implicit finite difference code, based on the Crank–Nicolson time advancement scheme and an iterative predictor corrector scheme, is employed (Desjardins et al. Reference Desjardins, Blanquart, Balarac and Pitsch2008). Spatial and temporal staggering is used to increase accuracy and stability. The Poisson equation for the pressure is solved by the preconditioned conjugate gradient HYPRE solver (Falgout & Yang Reference Falgout and Yang2002). Momentum equations are discretized with a fourth-order scheme. In the species and temperature equations, the convective term is discretized with a fifth-order WENO scheme (Jiang & Shu Reference Jiang and Shu1996) and the diffusion operator is discretized with fourth-order central differences. The time integration of the chemical source terms employs a table lookup. Note that the numerical code for the LES is the same as for the DNS.

The mesh of the LES is equidistant in all directions and is eight times coarser compared with the mesh of the DNS. Thus, the grid spacing is roughly as large as the laminar unstretched premixed flame thickness, which is $l_{F}= 0.714\,\text {mm}$ if determined based on the maximum temperature gradient. The grid resolution is chosen as coarse as possible to challenge the subfilter models, while it still needs to be fine enough to resolve the mean flow in the LES. The time step is controlled by the Courant–Friedrichs–Lewy (CFL) condition based on the local velocities and a maximum CFL number of 0.5 is allowed. Further, the maximum time step is limited to $\Delta t = 1\,\mathrm {\mu }{{\rm s}}$ in the LES to ensure an accurate time integration of the chemical reactions. The LES are run for at least four flow through times $t_{Flow}$, defined by the ratio of domain length $L_{x}$ and bulk velocity $U$ as $t_{Flow} = L_{x}/U$, and statistics are collected after a statistically steady state is reached.

5.3. LES results and model assessment

Figure 17 shows a snapshot of all LES that have been performed for the TurbUnstable DNS case of Berger et al. (Reference Berger, Attili and Pitsch2022b) that features thermodiffusive instabilities. In addition to the four different combustion models, the filtered DNS is also shown for comparison. Note that the filter size $\varDelta$, used to filter the DNS, is not equal to the LES grid, but is chosen as $\varDelta = 17 \delta _{DNS}$ as solving the filtered transport equations in LES leads to an implicit filtering with an unknown filter size a priori. The filter size is chosen such that the LES can be compared with a filtered field with similar local flame thickness and scalar gradients, which is further discussed in the following. From figure 17, significant differences among the different LES combustion models are evident. The C-model yields a significantly longer flame compared with the other LES and the DNS, which is expected as it does not capture the effects of differential diffusion, which have been shown to have a leading order impact on the flame dynamics and reaction rates (Berger et al. Reference Berger, Attili and Pitsch2022b). For the same reason, this combustion model is also not capable of reproducing the characteristic super-adiabatic temperatures that result from differential diffusion effects (Giannakopoulos et al. Reference Giannakopoulos, Gatzoulis, Frouzakis, Matalon and Tomboulides2015). In contrast, all other combustion models feature super-adiabatic temperatures as they intrinsically capture differential diffusion effects, e.g. by transporting two scalars such as progress variable and mixture fraction for the C-Z-models and in the Finitechem LES-model, each species possesses a different diffusion coefficient. However, differences among these three models are evident, e.g. the Finitechem LES model yields a significantly longer flame and the super-adiabatic temperature regions in the shear layer that start at the nozzle and evolve in the burned gas are only very weakly visible compared with the two flames with the C-Z-model and the DNS.

Figure 17.

Instantaneous snapshots of the temperature distribution for the different LES and the filtered DNS data.

To quantitatively highlight the difference in flame height, figure 18 shows the axial fuel flux $\mathcal {F}$ along the axial direction:

(5.1)\begin{equation} \mathcal{F} = \frac{1}{\rho_{u} U H Y_{H_2,u} }\int \langle \rho u Y_{H_2} \rangle_{z,t} \,\text{d}y, \end{equation}

where $\rho$ represents the gas density, $u$ the axial velocity, $Y_{H_2}$ is the mass fraction of hydrogen, and the operator $\langle \cdots \rangle _{z,t}$ represents averaging in the homogeneous $z$-direction and in time. For convenience, the fuel flux is normalized by its value at the inlet.

Figure 18.

Fuel mass flux along axial direction ($x$) for the different LES and the DNS.

From figure 18, a significantly different fuel consumption and, hence, flame length is evident among the different flames. Consistent with the flame images in figure 17, the LES with the C-model significantly overpredicts the flame length of the DNS. In contrast, the two LES with the C-Z-model approach the DNS in terms of fuel consumption and flame length sufficiently well even though slightly longer flames are obtained. Noteworthy, the Finitechem LES case yields a longer flame than the two cases with the C-Z-model, indicating that the missing subfilter modelling has a significant effect on the LES prediction.

Figure 19 shows profiles of the averaged values of the filtered temperature field $\langle \tilde {T} \rangle _{z,t}$ and its fluctuations $\sigma _{T}$ at different heights above the burner for the different cases. The fluctuations are defined as

(5.2)\begin{equation} \sigma_{T} = \sqrt{\langle( \tilde{T} - \langle \tilde{T} \rangle_{z,t})^2 \rangle_{z,t}}. \end{equation}

Note that $\sigma _{T}$ of the DNS is computed based on the filtered fields to allow for a direct comparison. Consistent with the fuel consumption and flame lengths in figure 18, different temperature profiles are observed at different heights above the burner with the DNS profiles approaching the burned state the fastest, e.g. at $x/H\approx 6$, the mean temperature is already close to the adiabatic flame temperature and only small fluctuations $\sigma _{T}$ are visible. As already pointed out, temperature overshoots are not visible for the LES with the C-model and only weakly present for the Finitechem LES case. For instance, at $x/H=1$, the Finitechem LES case yields $T_{max} = 1498\,{\rm K}$, which is 74 K lower compared with the overshoots in the DNS. However, differences are also observed among the two C-Z-models: the C-Z-model(H $_2$) features higher temperature overshoots with $T_{max} = 1650\,{\rm K}$ at $x/H=1$ compared with the DNS with $T_{max} = 1572\,{\rm K}$, while the temperature peak values of $T_{max} = 1593\,{\rm K}$ of the C-Z-model( $H_2O$) compare well with the DNS at $x/H=1$. Further, the C-Z-model( $H_2$) features an earlier widening of the jet compared with the C-Z-model( $H_2O$) model as visible from the wider temperature profile at $x/H = 3$ at $\langle \tilde {T} \rangle _{z,t} \approx 1200\,{\rm K}$. Similar trends are visible from the temperature fluctuations, e.g. higher temperature fluctuations are observed for the C-Z-model(H $_2$) than for the C-Z-model(H $_2$O), whereas the latter is close to the DNS values for $x/H\le 3$. The other two cases reveal significantly different fluctuation profiles due to the longer flames, e.g. the fluctuations monotonically increase for the C-model up to $x/H=10$ as the flame penetrates sufficiently far downstream in this case.

Figure 19.

(a–d) Mean temperature profiles and ( f–h) temperature fluctuations along lateral direction ($y$) at different heights above the burner for the different LES and DNS.

Figure 18 indicates a significantly different fuel consumption of the different flames, which is directly linked to a different consumption speed $s_{c}$. The latter is defined as (Poinsot & Veynante Reference Poinsot and Veynante2005)

(5.3)\begin{equation} s_{c} ={-} \frac{1}{\rho_{u} A_0 Y_{H_2,u}} \int_{\mathcal{V}} \dot{\omega}_{H_2} \,\text{d} V . \end{equation}

Here, $Y_{H_2,u}$ and $\rho _{u}$ are the hydrogen mass fraction and density in the unburned gas, $\dot {\omega }_{H_2}$ represents the fuel consumption rate, and $A_0$ refers to a reference flame surface area, which is the surface area of the Favre-averaged hydrogen mass fraction field. Due to the spatial inhomogeneity of the flow in the axial direction $x$, it is appropriate to define a local turbulent flame speed (Attili et al. Reference Attili, Stefano, Denker, Bisetti and Pitsch2020). Thus, the streamwise direction is divided into a number of equally sized volumes $\mathcal {V}$ to evaluate the volumetric integral of the source term in (5.3). Each of these extends along the entire spanwise $z$ and cross-wise direction $y$. The size of these volumes in the streamwise direction is as small as possible but finite such that converged statistics are obtained. It is worth noting that for the LES with the C-Z-model( $H_2O$) and C-model, which use a water-based progress variable, a consumption speed could be similarly defined based on the water production rate $\dot {\omega }_{H_2O}$, but no significant differences are observed if using such a definition, so for simplicity, all simulations are compared according to (5.3). For the cases based on flamelet models, the source term $\dot {\omega }_{H_2}$ is then determined from a table-lookup based on either $(\tilde {C}_{H_2},\widetilde {C^{''2}_{H_2}} , \tilde {Z} )$ or $(\tilde {C}_{H_2O},\widetilde {C^{''2}_{H_2O}} , \tilde {Z} )$. For the Finitechem LES case, $\dot {\omega }_{H_2}$ is evaluated with the detailed chemical mechanism without considering any subfilter model.

Figure 20 shows the consumption speed of the different cases along the axial direction. Consistent with the fuel consumption, the DNS case features the highest fuel consumption, similar values are observed for the two C-Z-models, and the low fuel consumption rates of the Finitechem LES and C-model cases lead to significantly longer flames compared with the other cases. To further analyse this behaviour, the consumption speed is decomposed as

(5.4)\begin{equation} s_{c} = s_{L} \frac{A}{A_0} I_0. \end{equation}

Here, $s_{L}$ is the laminar unstretched burning velocity, the ratio ${A}/{A_0}$ represents the effect of flame wrinkling on the turbulent flame speed, where $A$ denotes the averaged turbulent flame surface area, and $I_0$ is referred to as stretch factor and accounts for variations of the local reactivity from the unstretched laminar burning velocity and variations of the local flame thickness as introduced, e.g. by Bray (Reference Bray1990) and Bray & Cant (Reference Bray and Cant1991). To compute the averaged turbulent flame surface area, the instantaneous flame surface area $A_{inst.}$ in each volume $\mathcal {V}$ is determined from the hydrogen mass fraction according to Vervisch et al. (Reference Vervisch, Bidaux, Bray and Kollmann1995); technical details can be found from Berger et al. (Reference Berger, Attili and Pitsch2022b). Then, $A_{inst.}$ is averaged in time to yield $A$. It is worth noting that also for the C-Z-model( $H_2O$) and C-model, the flame surface area has been computed based on the hydrogen mass fraction field, which is determined from the look-up tables. This is important as the flame surface area determined from the water mass fraction field possesses spurious contributions due to long tails into the burned gas, cf. Berger et al. (Reference Berger, Attili and Pitsch2022a) and Howarth et al. (Reference Howarth, Hunt and Aspden2023). The stretch factor $I_0$ is then computed from (5.4).

Figure 20.

Consumption speed of the different LES and the filtered DNS along the axial direction $x$.

Figure 21 shows the increase of flame surface area along the axial direction for the different simulations. Consistent with the visual impression from figure 17, low flame wrinkling is seen for the C-model. Among the LES, the highest flame wrinkling is obtained in the C-Z-model( $H_2$), which visually also shows the most small-scale structures in figure 17, yielding similar flame wrinkling rates as in the filtered DNS. Figure 22 shows the stretch factor for all cases. For this, it is worth noting that the value of the filtered DNS ($I_0 \approx 6-5$), which is depicted in figure 22, is higher than the values of the unfiltered DNS ($I_0 \approx 5-4$, cf. Berger et al. Reference Berger, Attili and Pitsch2022b). This is a direct consequence of the filtering procedure: (i) filtering does not affect the consumption speed $s_{c}$ as it only redistributes the local source term; (ii) filtering leads to a reduction of flame surface area. Thus, the effects of subfilter wrinkling are accounted for in the stretch factor. While none of the LES can reproduce the high stretch factors of the DNS, the C-Z-model( $H_2O$) approaches closest the DNS and yields the highest values among the LES cases. From figures 21 and 22, it is evident that while the two C-Z-models yield similar overall fuel consumption, the local flame dynamics differ among the two cases, e.g. the C-Z-model( $H_2$) features higher flame wrinkling and the C-Z-model( $H_2O$) features higher overall reaction rates. However, both cases show a significant enhancement of flame wrinkling and reaction rates compared with the C-model, which does not account for differential diffusion effects, indicating a significant impact of the proposed model on both contributions.

Figure 21.

Flame surface area of the different LES and the filtered DNS along the axial direction $x$.

Figure 22.

Stretch factor of the different LES and the filtered DNS along the axial direction $x$.

To understand the variations of the stretch factor, the local flame state is further investigated in figure 23, which shows the joint distribution among the hydrogen-based progress variable and mixture fraction. The C-model is not shown as it does not feature any effects of differential diffusion and, hence, mixture fraction fluctuations. Among the different LES, the conditional mean $\langle \tilde {Z} \mid \tilde {C}_{H_2} \rangle$ of the C-Z-model( $H_2O$) qualitatively resembles the filtered DNS best as the C-Z-model( $H_2$) and Finitechem LES cases yield a strong overprediction compared with the filtered DNS at high progress variable values. Further, the narrow distribution of mixture fraction around the conditional mean in the Finitechem LES case indicates a reduced impact of differential diffusion effects in this case, which is consistent with the low super-adiabatic temperatures observed in figure 17. In contrast, significant mixture fraction fluctuations are visible for the C-Z-model( $H_2$), which are even larger towards the burned gas than the fluctuations in the DNS and explain the higher super-adiabatic temperatures compared with the DNS in figure 17. The C-Z-model( $H_2O$) represents an intermediate state, which features a sufficient widening of the joint distribution while extremal values of mixture fraction are similar to the DNS.

Figure 23.

Local flame state in the LES and filtered DNS characterized by joint distribution of $\tilde {C}_{H_2}$ and $Z$. (a) Filtered DNS. (b) Finitechem LES. (c) C-Z-model(${\rm H}_2{\rm O}$). (d) C-Z-model(${\rm H}_2$).

As reaction rates sensitively depend on mixture fraction and progress variable, a good representation of the local flame state is critical. In figure 24, the conditional mean fuel consumption rate $\langle \bar {\dot {\omega }}_{H_2} \mid \tilde {C}_{H_2} \rangle$ is depicted for all cases. It is evident that the C-model features the lowest source term as it does not account for the effects of differential diffusion. Note that for the C-model and C-Z-model( $H_2O$), the values of $\bar {\dot {\omega }}_{H_2}$ and $\tilde {C}_{H_2}$ are looked up as a function of $\tilde {C}_{H_2O}$, $\widetilde {C^{''2}_{H_2O}}$ and $\tilde {Z}$. All other LES tend to under predict the fuel consumption rate at low progress variable values, while they feature higher reaction rates towards the burned gas compared with the DNS. Both trends are linked to the joint distributions depicted in figure 23 as the DNS features the largest conditional mean mixture fraction values for small progress variable values, while the LES feature higher mixture fraction values towards the burned gas. However, to explain the differences in the stretch factor in figure 22, it is important to note that the source term is volumetrically integrated in (5.3). Thus, variations of the local flame thickness also affect the stretch factor. As flame thickness is directly linked to the local gradients within the flame, figure 25 shows the gradient of temperature conditionally averaged with respect to temperature for all LES and the filtered DNS. The filter size of the DNS has been chosen such that the magnitude of the maximum temperature gradient of the filtered field approximately matches the maximum values of the LES, or in other words, to match the implicit filter size characteristic of the LES field. Thus, note that a different filter size would lead to a different mean temperature gradient in figure 25 and conditionally averaged fuel consumption rate in figure 24. However, the particular filter size of $\varDelta = 17 \delta _{DNS}$ is chosen as the resulting temperature gradient and, hence, flame structure resembles the LES fields the closest for the comparison. The smallest gradients are observed for the C-model as a result of the small source term in figure 24. Further, the highest gradients among the LES are visible for the C-Z-model( $H_2$), which explains the lower value of the volumetric integral and, hence, stretch factor compared with the C-Z-model( $H_2O$) in figure 22. Most noteworthy, the DNS features a significant thickening of the mean flame front at low temperatures, which is not present in any LES. In conjunction with the relatively high reaction rates of the DNS for low progress variable values in figure 24, this leads to the high stretch factor values in figure 22.

Figure 24.

Conditionally averaged source term $\bar {\dot {\omega }}_{H_2}$ with respect to $\tilde {C}_{H_2}$ for all LES and the filtered DNS.

Figure 25.

Conditionally averaged temperature gradients with respect to temperature for all cases.

5.4. LES results of jet flame without differential diffusion effects

Finally, an additional LES of the DNS of Berger et al. (Reference Berger, Attili and Pitsch2022b), where differential diffusion effects are artificially suppressed, has been performed. The computational set-up of this DNS features also a slot burner configuration and the operational conditions are the same in terms of jet Reynolds number, Karlovitz number, unburned temperature, equivalence ratio and pressure as those discussed in § 5.1. However, the Lewis numbers of all species have been set to unity and the Soret effect has been neglected to suppress differential diffusion effects; further technical details of the simulations are discussed by Berger et al. (Reference Berger, Attili and Pitsch2022b). Performing LES of this flame with the C-model allows an assessment of the capability of the employed modelling framework to predict turbulent jet flames. As no differential diffusion is present, good predictability of the C-model is expected. It is worth noting that the unstretched laminar flamelet used to generate the table for the C-model in this case has been also computed with unity Lewis numbers and no Soret effect to be consistent with the DNS.

Figure 26 shows two snapshots of the temperature field of the filtered DNS and LES. Analogously to the previous discussion, the filter size is chosen as $\varDelta = 17 \delta _{DNS}$ such that gradients in LES and filtered DNS are similar in magnitude. Additionally, figure 27 shows the fuel mass flux along the axial distance. Generally, good agreement of the filtered DNS and C-model is observed even though the flame length is overpredicted in the LES. In particular, the C-model overpredicts the flame length of the DNS without differential diffusion effects by approximately 25 %, while for the DNS with differential diffusion effects, the C-model yields a flame, which is approximately 270 % longer than the DNS. Here, the flame length is defined as the axial distance, where 99 % of the fuel is consumed, e.g. at $\mathcal {F}=0.01$ in figures 18 and 27. This indicates that the employed modelling framework is suitable to predict turbulent premixed jet flames and the aforementioned discrepancies of the C-model to predict the DNS with strong thermodiffusive instabilities can be clearly linked to the fact that differential diffusion effects are not considered.

Figure 26.

Instantaneous snapshots of the temperature fields for the TurbStable DNS without differential diffusion effects and the corresponding LES.

Figure 27.

Fuel mass flux along axial direction ($x$) for the TurbStable DNS without differential diffusion effects and the corresponding LES.