A Combustion Regime-Based Model for Large Eddy Simulation

Abstract

The aim of this work is to propose a unified (generalized) closure of the chemical source term in the context of Large Eddy Simulation able to cover all the regimes of turbulent premixed combustion. Turbulence/combustion scale interaction is firstly analyzed: a new perspective to look at commonly accepted combustion diagrams is provided based on the evidence that actual turbulent flames can experience locally several combustion regimes although global non-dimensional numbers would locate such flames in a single specific operating point of the standard combustion diagram. The deliverable is a LES subgrid scale model for turbulent premixed flames named Localized Turbulent Scales Model (LTSM). This is founded on the estimation of the local reacting volume fraction of a computational cell that is related to the local turbulent and laminar flame speeds and to the local flame thickness.

1. Introduction

Combustion regimes have been theoretically investigated for many years [1,2,3,4,5,6,7,8] by simply assuming the turbulence integral length scale, L, the associated turbulent velocity fluctuation, $u_{rms}^{\prime }$, the laminar flame speed, ${\mathcal{S}}_{\mathcal{L}}$, and the flame front thickness, ${\delta }_L$, as the main quantities characterizing turbulence–chemistry interaction. Such analysis, strictly valid for premixed flames but extendable to non-premixed flames too [9], leads to an order of magnitude definition of the combustion regimes in the so-called Klimov–Williams combustion diagram [1] that reports $log\left(u_{rms}^{\prime }/{\mathcal{S}}_{\mathcal{L}}\right)$ vs. $log\left(L/{\delta }_L\right)$. A brief history of combustion diagrams is provided in the Supplementary Materials.

The standard combustion diagram and its spectral evolution [10] do not include the effect of important flame physics such as heat losses, flame curvature, viscous dissipation, and transient dynamics, all affecting quenching. The effect of the Lewis number on quenching produced by stretching is also not considered [11]. In the literature, there are several numerical simulations [12,13] and experimental studies [14,15,16] on the effect of aerodynamic stretch on a laminar flame front (often a stagnation-point laminar flame [17]). These studies predict stretch factors that produce extinction (and therefore are called critical stretch factors), showing that quenching is favored by $Le>1$ and by flame non-adiabaticity.

In addition, it is observed that there can be several local sources of turbulence even for a simple single jet configuration, thus providing different and local turbulent macro-scales. For example, turbulence can be modified or produced by acoustic waves [18]. It can be damped by the flame front, but it can also be locally produced or enhanced by the flame front itself: the dilatational effects due to the heat release may even result into a non-negligible energy backscatter from small to large scales in combustion regimes in which the flame–transit time is long enough to allow for activation of the nonlinear convective mechanisms of the energy cascade, as recently confirmed through DNS data analysis [19] and revealed by some previous energy spectra [20]. In flames, there is also a spectrum of chemical times, depending on the local state of the flow in terms of pressure, temperature, and composition. This spectrum of chemical times is necessarily reduced just to one chemical time if ${\mathcal{S}}_{\mathcal{L}}$ and ${\delta }_L$ are assumed constant in the analysis of combustion diagrams [21]; this also implies to neglect the effect of conductivity (or diffusivity more generally) on local laminar flame speed [22]. Furthermore, standard theory assumes unity Prandtl number.

As a matter of fact, a real flame can experience several regimes depending on the local turbulent and chemical conditions. For example, Figure 1 reports on the standard combustion diagram the global states (the three big yellow marked points) of three ${\mathrm{CH}}_4/{\mathrm{H}}_2/\mathrm{Air}$ (equivalence ratio 0.7) premixed flames at 1, 10, and 40 bar, derived by DNS simulation data [23]. On the same diagram, combustion regimes associated with an instantaneous flowfield of the three flames are also reported; for these clouds of points, the quantities along the axis now have to be read as instantaneous quantities. Instantaneous local turbulent velocity and correlation length scales were calculated in the fresh mixture at the $c=0.1$ value of the progress variable; the instantaneous local laminar flame speed was estimated according to the expression of Bougrine [24]; the instantaneous local flame thickness was calculated as the path length normal to the iso-c surfaces from 0.1 to 0.99. It is observed how the actual combustion regimes are spread around the global state point, covering different regions of the combustion map. It is also observed that the instantaneous distributions of the Prandtl number (not reported) exhibit $Pr\ne 1$ in every flames, and the distribution becomes wider by increasing pressure [23].

A central point in turbulent combustion modeling is the closure of the filtered chemical source term. Some models were designed to work in specific combustion regimes, e.g., Zimont [25,26,27] and Peters [21,28] developed turbulent flame speed (${\mathcal{S}}_{\mathcal{T}}$) models for a highly turbulent premixed combustion regime characterized by “flamelets” thickened by small-scale turbulent eddies penetrating into their preheating zone; Bray proposed how to estimate the domain of countergradient transport through his own characteristic number [29,30,31] and specific models for this mechanism mainly associated with weak, medium-scale turbulence were suggested [32,33]; Clavin introduced the effect of flame front stretch rate and curvature on flame speed for weakly perturbed laminar flames [34,35], later extended by Matalon to the weak turbulence (Markstein) regime [36]; Minamoto suggested a closure for the filtered reaction rate in LES of MILD combustion [37]. An important feature of the models adopted for such turbulent combustion closure would be their self-adaptation to the local combustion regime.

Upon the above literature review, it is clear that there exist different models specific for each combustion regime. A unified model covering all the regimes is currently missing. This motivates the present work whose final objectives are to build such a general framework, here called Localised Turbulent Scales Model (LTSM), and to develop a subgrid Large Eddy Simulation model for the chemical source terms, self-adapting to the local combustion regimes.

To reach its objectives, this work will be based on theoretical analysis of turbulent combustion and some well known results in specific combustion regimes. The article consists of two parts: in the first one, a methodology to identify the different combustion regimes locally exhibited by turbulent premixed flames will be presented (LTSM); in the second one, a subgrid scale model for the filtered chemical source terms will be provided and validated.

2. LTSM Theoretical Framework

The methodology adopted for the development of the theoretical framework of the Localised Turbulent Scales Model (LTSM) is presented here. The goal of the first part of the present article is the identification of premixed combustion regimes based on comparison of orders of magnitude (comparing the physical length and velocity scales involved in the process and reckoning on their possible interaction) and some assumptions. In particular, it is aimed at identifying combustion regimes by specifying their boundaries in terms of the turbulent Reynolds and Damköhler numbers and removing the unity Prandtl number assumption. In fact, actual local Prandtl number distributions in flames show values lower than one in ideal gas conditions; values greater than one are also possible in real gas conditions [38].

Identification of combustion regimes means to understand the spatial structure and morphology of the reacting flow resulting from the interaction of turbulent and chemical scales. Hence, the methodology adopted will consist of identifying when the flame can be considered locally laminar or turbulent, when combustion can be assumed to be volumetric or distributed, how eddies smaller or larger than the flame front can affect it, and then which eddies may survive to a flame front without being dissipated before combustion can take place at their scale, and among these ”survived” eddies which may wrinkle, thicken, or quench the flame front.

2.1. The Combustion Regime Identification Strategy

The analysis will take into account characteristic times at the generic eddy scale ℓ; this eddy is assumed to belong to the Kolmogorov’s inertial range, $\left[\eta ,{\ell }_{\Delta }\right]$ (${\ell }_{\Delta }$ being the local macro-scale and $\eta$ the local dissipative scale) and thus it will have a characteristic velocity $u_{\ell }^{\prime }=u_{\Delta }^{\prime }{\left(\ell /{\ell }_{\Delta }\right)}^{1/3}$. Characteristic times of different processes will be considered, such as the turbulent time ${\tau }_{t_{\ell }}=\ell /u_{\ell }^{\prime }$, the viscous time ${\tau }_{{\nu }_{\ell }}={\ell }^2/\nu$, the chemical time ${\tau }_{ch}={\delta }_L/{\mathcal{S}}_{\mathcal{L}}=\alpha /{\mathcal{S}}_{\mathcal{L}}^2={\delta }_L^2/\alpha$ (the Mallard–Le Chatelier theory is adopted), $\alpha =\kappa /\left(\rho C_p\right)$ being the thermal diffusivity, $\kappa$ the thermal conductivity, $\rho$ the density, and $C_p$ the specific heat at constant pressure. As a reminder, the laminar flame speed in the Mallard–Le Chatelier theory is estimated as follows:

$${\mathcal{S}}_{\mathcal{L}}\approx {\left(\frac{\alpha }{{\tau }_{ch}}\right)}^{1/2},$$

(1)

the laminar flame thickness as follows:

$${\delta }_{\mathcal{L}}\approx {\left(\alpha {\tau }_{ch}\right)}^{1/2},$$

(2)

and that these two expressions imply the following equation:

$$\frac{{\delta }_{\mathcal{L}}{\mathcal{S}}_{\mathcal{L}}}{\alpha }=1.$$

(3)

The Prandtl number, the turbulent Reynolds, and Damköhler (first and second species) numbers will be defined from characteristic times and used to identify the boundaries of combustion regimes.

In the following sections, the scale ${\ell }^{*}$ of the smallest “surviving” eddy (not dissipated by the flame front) will be firstly identified, and then checked if it may wrinkle the flame front, thus becoming a “wrinkling” eddy of scale ${\ell }_w^{*}$. Then, it will be investigated if “survived” eddies of scale $\ell <{\delta }_L$ may thicken the flame front itself. Finally, comparing the newly defined length scales, combustion regimes will be naturally identified.

2.2. The Combustion Closure Strategy

The Favre filtered chemical source term in the energy and single species transport equations is modeled here as $\widetilde{{\omega }_i}\approx {\gamma }^{*}{\omega }_i^{*}$, ${\gamma }^{*}$ and ${\omega }_i^{*}$ being the local reacting volume fraction of the computational cell and the reaction rate of the $i-th$ chemical species.

Let us consider a local control volume ${\mathcal{V}}_{\Delta }$ with characteristic dimension ${\ell }_{\Delta }={\mathcal{V}}_{\Delta }^{1/3}$ and assume that it contains a flame front having an actual surface area ${\mathcal{A}}_{\mathcal{F}}$ and an actual thickness ${\delta }_{\mathcal{F}}$ and hence a reacting volume ${\mathcal{V}}_{\mathcal{F}}$. It is possible to define the fraction ${\gamma }^{*}$ of the local reacting volume ${\mathcal{V}}_{\mathcal{F}}$ as follows:

$${\gamma }^{*}=\frac{{\mathcal{V}}_{\mathcal{F}}}{{\mathcal{V}}_{\Delta }}\approx \frac{{\mathcal{A}}_{\mathcal{F}}{\delta }_{\mathcal{F}}}{{\mathcal{V}}_{\Delta }}\approx \frac{{\mathcal{S}}_{\mathcal{T}}}{{\mathcal{S}}_{\mathcal{L}}}{\mathcal{A}}_{\mathcal{L}}\frac{{\delta }_{\mathcal{F}}}{{\mathcal{V}}_{\Delta }}\approx \frac{{\mathcal{S}}_{\mathcal{T}}}{{\mathcal{S}}_{\mathcal{L}}}{\ell }_{\Delta }^2\frac{{\delta }_{\mathcal{F}}}{{\ell }_{\Delta }^3}=\frac{{\mathcal{S}}_{\mathcal{T}}}{{\mathcal{S}}_{\mathcal{L}}}\frac{{\delta }_{\mathcal{F}}}{{\ell }_{\Delta }}.$$

(4)

This expression has been obtained with two main assumptions. The first is that, within a wrinkled flame front, the iso-surfaces of the progress variable are parallel [39]. The second assumption is that the ratio between the turbulent and the laminar flame surface areas scales as the ratio between the associated flame speeds and is expressed as follows:

$$\frac{{\mathcal{A}}_{\mathcal{F}}}{{\mathcal{A}}_{\mathcal{L}}}\equiv \frac{{\mathcal{A}}_{\mathcal{T}}}{{\mathcal{A}}_{\mathcal{L}}}\approx \frac{{\mathcal{S}}_{\mathcal{T}}}{{\mathcal{S}}_{\mathcal{L}}}.$$

(5)

Subgrid flame front physics is synthesized in this ratio: the subgrid flame may be laminar or turbulent, wrinkled or not, thickened by turbulence or not, depending on the local conditions of the flow. The ${\gamma }^{*}$ estimation becomes the problem of estimating the characteristics of the local flame front (turbulent flame speed, laminar flame speed and thickness, turbulent or laminar) from the filtered conditions of the flow and depending on the related local premixed combustion regime.

An extinction or flame stretch factor ${\mathcal{G}}_{ext}\le 1$, introduced in Section 4, is also included in Equation (4) to take into account flame quenching due to subgrid scales. Hence, the local reacting volume fraction is finally modeled as follows:

$${\gamma }^{*}={\mathcal{G}}_{ext}\frac{{\mathcal{S}}_{\mathcal{T}}}{{\mathcal{S}}_{\mathcal{L}}}\frac{{\delta }_{\mathcal{F}}}{\Delta }.$$

(6)

3. Vortices/Flame Front Interaction

This section aims to answer two fundamental questions related to the interaction problem between vortices and a flame front. What is the smallest eddy with a life-time longer than the characteristic chemical time? In addition, then, among the surviving scales, what is the smallest eddy able to locally wrinkle a flame front? Answering the first question allows for defining the range of scales that can interact with a flame front, and eventually enter into it and produce local flame thickening. Answering the second question allows instead for defining the smallest turbulent scales that apply the highest strain rate and curvature onto the flame front.

3.1. The Smallest Surviving Eddy

The interaction between a premixed flame front and eddies has been widely analyzed in the literature. Results clearly show that the dissipative Kolmogorov scales $\eta$ cannot quench a flame front [22]. An estimate of the smallest turbulent scale that can affect a laminar flame front without being dissipated can be obtained by imposing that the turbulent $\ell -$scale Damköhler number of second species be greater than one, i.e.,

$$D{a}_{\ell }^{II}=\frac{{\tau }_{{\nu }_{\ell }}}{{\tau }_{ch}}=P{r}^{-1}{\left(\frac{\ell }{{\delta }_{\mathcal{L}}}\right)}^2\ge 1,$$

(7)

where ${\tau }_{{\nu }_{\ell }}={\ell }^2/\nu$ is the lifetime of the vortex of scale ℓ, ${\tau }_{ch}={\delta }_{\mathcal{L}}/{\mathcal{S}}_{\mathcal{L}}={\delta }_{\mathcal{L}}^2/\alpha$ is the chemical time, and $Pr=\nu /\alpha$ the Prandtl number.

It is stressed that the present analysis removes the $Pr=1$ assumption. Hence, the smallest surviving scale $l^{*}$ is estimated as follows:

$${\ell }^{*}=P{r}^{1/2}{\delta }_{\mathcal{L}}={\ell }_{\Delta }{\left(D{a}_{\Delta }^{I}R{e}_{\Delta }\right)}^{-1/2}={\ell }_{\Delta }D{a}_{\Delta }^{I{I}^{-1/2}},$$

(8)

where $R{e}_{\Delta }=u_{\Delta }^{\prime }{\ell }_{\Delta }/\nu$ is the turbulent Reynolds number defined in terms of the local length and velocity macro-scales, ${\ell }_{\Delta }$ and $u_{\Delta }^{\prime }$; $D{a}_{\Delta }^I={\ell }_{\Delta }/\left(u_{\Delta }^{\prime }{\tau }_{ch}\right)$ is the turbulent ${\ell }_{\Delta }-$scale Damköhler number of the first species.

The smallest surviving scale ${\ell }^{*}\in [\eta ,{\ell }_{\Delta }]$ for whatever Prandtl number if the following condition is met:

$${\ell }_{\Delta }\ge {\ell }^{*}\ge \eta \iff R{e}_{\Delta }^{-1}\le D{a}_{\Delta }^I\le R{e}_{\Delta }^{1/2}\Rightarrow D{a}_{\Delta }^I\ge R{e}_{\Delta }^{-1},$$

(9)

with ${\ell }^{*}={\ell }_{\Delta }$ for $D{a}_{\Delta }^I=R{e}_{\Delta }^{-1}$ and ${\ell }^{*}=\eta$ for $D{a}_{\Delta }^I=R{e}_{\Delta }^{1/2}$. It is observed that, for $Pr=1$, it is always ${\ell }^{*}\equiv {\delta }_{\mathcal{L}}$. Eddies $\ell \ge {\ell }^{*}$ will not be dissipated by the flame front: for whatever Prandtl number, no eddies will survive for $D{a}_{\Delta }^I<R{e}_{\Delta }^{-1}$, while the whole spectrum of eddies will for $D{a}_{\Delta }^I\ge R{e}_{\Delta }^{1/2}$.

An important observation, of course valid in the present order of magnitude framework, is as follows:

$${\ell }^{*}\le {\delta }_{\mathcal{L}}\iff Pr\le 1,$$

(10)

and this means that, only for $Pr\le 1$, typical of ideal gas combustion, eddy-scales smaller than the flame thickness may survive and affect the flame itself, e.g., entering into it (although this has not been proved yet) and thickening it.

Another observation is that the Taylor micro-scale, ${\lambda }_g={\ell }_{\Delta }R{e}_{\Delta }^{-1/2}$, could be used to estimate the local strain of the flow as $u_{\Delta }^{\prime }/{\lambda }_g$ in the eddy–flame interaction (flame strain modeling), according to the equation:

$${\lambda }_g\ge {\ell }^{*}\iff D{a}_{\Delta }^I\ge 1.$$

(11)

This relation is always verified in the range of existence of ${\ell }^{*}$ in (9) since $1>R{e}_{\Delta }^{-1}$ (assuming the flow locally turbulent, i.e., $R{e}_{\Delta }>1$).

3.2. The Smallest Wrinkling Eddy

Looking at the vortices/flame front interaction problem, some eddies will be destroyed by viscous dissipation before combustion can take place. The range of the surviving eddy-scales has been defined in Section 3.1. Among these scales, those smaller than the local flame front thickness may locally enter into it and thicken it, provided that $u_{\Delta }^{\prime }\ge {\mathcal{S}}_{\mathcal{L}}$, $u_{\Delta }^{\prime }$ being the local $rms$ velocity fluctuation. The “surviving” condition is not sufficient for an eddy to locally wrinkle a flame front. Extending this concept to the generic scale l, only eddies with characteristic (tangential) velocity $u_{\ell }^{\prime }\ge {\mathcal{S}}_{\mathcal{L}}$ will be able to locally wrinkle a flame front.

Given the smallest surviving eddy ${\ell }^{*}={\ell }_{\Delta }{\left(D{a}_{\Delta }^{I}R{e}_{\Delta }\right)}^{-1/2}$ within the inertial cascade, the smallest surviving and wrinkling eddy ${\ell }_w^{*}$ will be given by the following condition:

$$\frac{u_{{\ell }^{*}}^{\prime }}{{\mathcal{S}}_{\mathcal{L}}}\ge 1\Rightarrow \frac{u_{\Delta }^{\prime }}{{\mathcal{S}}_{\mathcal{L}}}{\left(\frac{{\ell }^{*}}{{\ell }_{\Delta }}\right)}^{1/3}=P{r}^{1/2}{\left(R{e}_{\Delta }D{a}_{\Delta }^{I^{-2}}\right)}^{1/3}\ge 1,$$

(12)

Hence, the smallest surviving scale ${\ell }^{*}$ becomes ${\ell }_w^{*}$ only if the following condition is met:

$$D{a}_{\Delta }^I\le P{r}^{3/4}R{e}_{\Delta }^{1/2}\Rightarrow \frac{u^{\prime }}{{\mathcal{S}}_{\mathcal{L}}}\ge P{r}^{-1/6}{\left(\frac{{\ell }_{\Delta }}{{\delta }_{\mathcal{L}}}\right)}^{1/3}.$$

(13)

For higher Damköhler numbers, the smallest surviving scale ${\ell }^{*}={\ell }_{\Delta }{\left(D{a}_{\Delta }^{I}R{e}_{\Delta }\right)}^{-1/2}$ decreases but is unable to wrinkle the flame front; in this higher $D{a}_{\Delta }^I$ range, the smallest surviving and wrinkling scale becomes constant and equal to the minimum value ${\left.{\ell }_w^{*}\right|}_{min}$ reached at $D{a}_{\Delta }^I=P{r}^{3/4}R{e}_{\Delta }^{1/2}$. Substituting this maximum Damköhler number for the smallest surviving and wrinkling scale into Equation (8), the minimum ${\ell }_w^{*}$ is obtained as follows:

$${\left.{\ell }_w^{*}\right|}_{min}=P{r}^{-3/8}\eta \ge \eta .$$

(14)

Since inner cut-off scales are generally expressed in terms of the Karlovitz number, defined in Equation (26), the relation (14) can be rewritten as follows:

$${\left.{\ell }_w^{*}\right|}_{min}/{\delta }_{\mathcal{L}}=P{r}^{-3/8}K{a}_t^{-1/2}.$$

(15)

This expression has the same $K{a}_t$ dependence as in [40], with $P{r}^{-3/8}$ as the proportionality coefficient.

In summary, ${\ell }_w^{*}$ decreases as ${\ell }_{\Delta }{\left(D{a}_{\Delta }^{I}R{e}_{\Delta }\right)}^{-1/2}$ by increasing $D{a}_{\Delta }^I$ up to $D{a}_{\Delta }^I=P{r}^{3/4}R{e}_{\Delta }^{1/2}$. Then, ${\ell }_w^{*}$ maintains its minimum value $P{r}^{-3/8}\eta$ for whatever $D{a}_{\Delta }^I\ge P{r}^{3/4}R{e}_{\Delta }^{1/2}$, implying that the highest local flame front curvature can be of the order ${\ell }_w^{{*}^{-1}}=P{r}^{3/8}/\eta$.

From what was found above and considering ${\left.u_{{\ell }^{*}}^{\prime }\right|}_{min}={\mathcal{S}}_{\mathcal{L}}$, the minimum scale Reynolds number for flame front wrinkling will be as follows:

$${\left.R{e}_{{\ell }_w^{*}}\right|}_{min}=\frac{{\mathcal{S}}_{\mathcal{L}}{\ell }^{*}}{\nu }\equiv \frac{{\mathcal{S}}_{\mathcal{L}}{\left.{\ell }_w^{*}\right|}_{min}}{\nu }=P{r}^{-1/2}\ge 1.$$

(16)

It is observed that simple chemistry DNS calculations [22,40] showed the following equation:

$$\frac{{\left.{\ell }_w^{*}\right|}_{min}}{{\delta }_{\mathcal{L}}}=0.2+5.5{\left(\frac{\epsilon {\delta }_{\mathcal{L}}}{{\mathcal{S}}_{\mathcal{L}}^3}\right)}^{-1/6}\Rightarrow {\left.R{e}_{{\ell }_w^{*}}\right|}_{min}\approx 5.5^{3/2}D{a}_{\ell }^{I^{-1/2}}P{r}^{-1},$$

(17)

with $\epsilon =u_{\ell }^{\prime 3}/\ell$. Experimental work [41,42] instead showed the following equation:

$$\frac{{\left.u_{{\ell }_w^{*}}^{\prime }\right|}_{min}}{{\mathcal{S}}_{\mathcal{L}}}=2.5{\left(\frac{{\left.{\ell }_w^{*}\right|}_{min}}{{\delta }_{\mathcal{L}}}\right)}^{-1}\Rightarrow {\left.R{e}_{{\ell }_w^{*}}\right|}_{min}\approx 2.5P{r}^{-1}.$$

(18)

Large differences are expected between experiments, simple chemistry DNS calculations, and theoretical results. However, the trends of these three approaches are qualitatively well reproduced, looking at the ${\left.R{e}_{{\ell }_w^{*}}\right|}_{min}=f\left(Pr\right)$ dependence.

3.3. The $\mathcal{T}$urbulence-$\mathcal{T}$hickened $\mathcal{R}$egime

In Section 3.1, it was shown that some scales smaller than the flame thickness may not be dissipated. These scales may affect the flame front itself thickening it. In particular, it was theorized that eddies ${\ell }^{*}\le {\delta }_{\mathcal{L}}$ exist only for $Pr\le 1$ and hence the $\mathcal{T}$urbulence-$\mathcal{T}$hickened $\mathcal{R}$egime cannot be experienced when $Pr>1$.

In the literature, there is a well known model for the turbulence-thickened regime: it is the Zimont’s model [25,26,27]. According to this model, the chemical time does not change through the whole thickening process: considering the flame front composed of a reaction layer, ${\delta }_r$, and a thermal (preheating) layer, ${\delta }_{tl}$, a constant ${\tau }_{ch}$ implies that ${\delta }_r$ is not affected by the flow, while ${\delta }_{tl}$ may be thickened by the survived eddies.

The thickening process can be modeled via a turbulent thermal diffusivity ${\alpha }^{*}\approx u^{*}{\ell }_{tk}^{*}$, $u^{*}$ being the characteristic velocity of the scale ${\ell }_{tk}^{*}\approx u^{*}{\tau }_{ch}$ that is the survived and thickening turbulent scale obtained via Equation (7) with a turbulent viscosity ${\nu }^{*}\approx u^{*}{\ell }_{tk}^{*}$. With these assumptions, the thickened flame front ${\delta }_{\mathcal{F}}^{*}$ can be estimated as follows:

$${\delta }_{\mathcal{F}}^{*}\approx {\left({\alpha }^{*}{\tau }_{ch}\right)}^{1/2}\approx {\left(u^{*}{\ell }_{tk}^{*}{\tau }_{ch}\right)}^{1/2}\approx u^{*}{\tau }_{ch}\approx {\ell }_{\Delta }D{a}_{\Delta }^{I^{-3/2}},$$

(19)

having obtained $u^{*}$ via the turbulence energy cascade, $u^{*}\approx u_{\Delta }^{\prime }{({\ell }_{tk}^{*}/{\ell }_{\Delta })}^{1/3}\approx u_{\Delta }^{\prime }D{a}_{\Delta }^{I^{-1/2}}$. It is observed that ${\ell }_{tk}^{*}$ has to be larger than or equal to ${\ell }^{*}$ given by Equation (8) and, in particular, it comes out that ${\ell }_{tk}^{*}\equiv {\delta }_{\mathcal{F}}^{*}$. The same result for ${\delta }_{\mathcal{F}}^{*}$ is obtained, as in Zimont’s model [25,27], by assuming ${\alpha }^{*}\approx u^{*}{\ell }_{tk}^{*}\approx u^{*}{\delta }_{\mathcal{F}}^{*}$, i.e., that the smallest surviving scale is of the order of the thickened flame front.

From Equation (19), the laminar thickened flame speed is as follows:

$${\mathcal{S}}_{\mathcal{L}}^{*}\approx {\left(\frac{{\alpha }^{*}}{{\tau }_{ch}}\right)}^{1/2}\approx {\left(\frac{u^{*}{\delta }_{\mathcal{F}}^{*}}{{\tau }_{ch}}\right)}^{1/2}\approx u_{\Delta }^{\prime }D{a}_{\Delta }^{I^{-1/2}}\equiv u^{*}.$$

(20)

Assuming ${\mathcal{S}}_{\mathcal{T}}/{\mathcal{S}}_{\mathcal{L}}\approx {\mathcal{A}}_{\mathcal{T}}/{\mathcal{A}}_{\mathcal{L}}$ and assuming the ratio between the turbulent and laminar flame surface areas, ${\mathcal{A}}_{\mathcal{T}}/{\mathcal{A}}_{\mathcal{L}}\approx D{a}_{\Delta }^{I^{3/4}}$ [25,27], the turbulent thickened flame speed can be scaled as

$${\left.{\mathcal{S}}_{\mathcal{T}}\right|}_{\mathcal{Z}}\approx {\mathcal{A}}_{\mathcal{Z}}{\mathcal{S}}_{\mathcal{L}}^{*}\frac{{\mathcal{A}}_{\mathcal{T}}}{{\mathcal{A}}_{\mathcal{L}}}\approx {\mathcal{A}}_{\mathcal{Z}}{\mathcal{S}}_{\mathcal{L}}^{*}D{a}_{\Delta }^{I^{3/4}}\approx {\mathcal{A}}_{\mathcal{Z}}{u}_{\Delta }^{\prime }D{a}_{\Delta }^{I^{1/4}},$$

(21)

${\mathcal{A}}_{\mathcal{Z}}=0.5$ being an empirical constant [26,43]. Assuming that the chemical time is unchanged, i.e., ${\delta }_{\mathcal{F}}^{\mathcal{T}}/{\left.{\mathcal{S}}_{\mathcal{T}}\right|}_{\mathcal{Z}}={\tau }_{ch}$, the thickness of the turbulence-thickened flame front is as follows:

$${\left.{\delta }_{\mathcal{F}}^{\mathcal{T}}\right|}_{\mathcal{Z}}\approx {\left.{\mathcal{S}}_{\mathcal{T}}\right|}_{\mathcal{Z}}{\tau }_{ch}\approx {\mathcal{A}}_{\mathcal{Z}}{\ell }_{\Delta }D{a}_{\Delta }^{I^{-3/4}}.$$

(22)

In terms of Prandtl and turbulent Reynolds numbers, ${\left.{\mathcal{S}}_{\mathcal{T}}/{\mathcal{S}}_{\mathcal{L}}\right|}_{\mathcal{Z}}$ becomes the following equation:

$${\left.\frac{{\mathcal{S}}_{\mathcal{T}}}{{\mathcal{S}}_{\mathcal{L}}}\right|}_{\mathcal{Z}}\approx {\mathcal{A}}_{\mathcal{Z}}{\left(PrR{e}_{\Delta }\right)}^{1/2}D{a}_{\Delta }^{I^{-1/4}}.$$

(23)

Imposing the conditions reported in

Table 1. Constraints to determine the validity range of the $\mathcal{T}urbulence$–$\mathcal{T}hickened$ flame $\mathcal{R}egime$.

List	Condition	Type of Condition		Range
1	${\ell }_{\Delta }\ge {\ell }^{*}\ge \eta$	surviving scale existence	⇒	$D{a}_{\Delta }^I\ge R{e}_{\Delta }^{-1}$
2	${\ell }^{*}\le {\delta }_{\mathcal{L}}$	potential thickening scales	⇒	$Pr\le 1$
3	${\delta }_{\mathcal{L}}\le {\ell }_{\Delta }$	no volumetric combustion 1	⇒	$D{a}_{\Delta }^I\ge P{r}^{-1}R{e}_{\Delta }^{-1}$
4	${\delta }_{\mathcal{L}}\ge \eta$	minimum thickening scale 1	⇒	$D{a}_{\Delta }^I\le P{r}^{-1}R{e}_{\Delta }^{1/2}$
5	$u_{\Delta }^{\prime }\ge {\mathcal{S}}_{\mathcal{L}}$	velocity fluctuation	⇒	$D{a}_{\Delta }^I\le P{r}^{-1}R{e}_{\Delta }$
6	${\ell }_{tk}^{*}\ge {\ell }^{*}$	minimum thickening scale 2	⇒	$D{a}_{\Delta }^I\le R{e}_{\Delta }^{1/2}$
7	${\delta }_{\mathcal{F}}^{*}\ge {\delta }_{\mathcal{L}}$	thickening effect 1	⇒	$D{a}_{\Delta }^I\le P{r}^{1/2}R{e}_{\Delta }^{1/2}$
8	${\mathcal{S}}_{\mathcal{L}}^{*}\ge {\mathcal{S}}_{\mathcal{L}}$	thickening effect 1	⇒	$D{a}_{\Delta }^I\le P{r}^{1/2}R{e}_{\Delta }^{1/2}$
9	${\delta }_{\mathcal{F}}^{\mathcal{T}}\ge {\delta }_{\mathcal{F}}^{*}$	thickening effect 2	⇒	$D{a}_{\Delta }^I\ge 1$
10	${\left.{\mathcal{S}}_{\mathcal{T}}\right|}_{\mathcal{Z}}\ge {\mathcal{S}}_{\mathcal{L}}^{*}$	thickening effect 2	⇒	$D{a}_{\Delta }^I\ge 1$
11	${\delta }_{\mathcal{F}}^{\mathcal{T}}\le {\ell }_{\Delta }$	no volumetric combustion 2	⇒	$D{a}_{\Delta }^I\ge 1$
12	${\left.\frac{{\mathcal{S}}_{\mathcal{T}}}{{\mathcal{S}}_{\mathcal{L}}}\right|}_{\mathcal{Z}}\frac{{\delta }_{\mathcal{F}}^{*}}{{\ell }_{\Delta }}\le 1$	reacting volume fraction	⇒	$D{a}_{\Delta }^I\ge P{r}^{2/7}R{e}_{\Delta }^{2/7}$

, the $\mathcal{T}urbulence$ - $\mathcal{T}hickened$ flame $\mathcal{R}egime$ validity ranges are obtained in terms of the local Damköhler number, $D{a}_{\Delta }^I$. It is observed that the Zimont constant ${\mathcal{A}}_{\mathcal{Z}}$ was not considered in this analysis, specifically in conditions 9–12.

In particular, the constraint number 12 in Table 1 requires that the volume fraction ${\gamma }^{*}$ in Equation (4) be less than one. Having adopted a model for the reacting volume fraction, subgrid flame front wrinkling and curvature effects are synthesized in the quantity ${\gamma }^{*}$. Considering the Zimont’s model [25,26,27], the expression for the local reacting volume fraction ${\gamma }^{*}$ in the turbulent thickened flame regime becomes as follows:

$${\gamma }^{*}=\frac{{\mathcal{S}}_{\mathcal{T}}}{{\mathcal{S}}_{\mathcal{L}}}\frac{{\delta }_{\mathcal{F}}}{{\ell }_{\Delta }}={\left.\frac{{\mathcal{S}}_{\mathcal{T}}}{{\mathcal{S}}_{\mathcal{L}}}\right|}_{\mathcal{Z}}\frac{{\delta }_{\mathcal{F}}^{*}}{{\ell }_{\Delta }},$$

(24)

with ${\delta }_{\mathcal{F}}^{*}$ and ${\left.{\mathcal{S}}_{\mathcal{T}}/{\mathcal{S}}_{\mathcal{L}}\right|}_{\mathcal{Z}}$ given by Equations (19) and (21), respectively.

Imposing the constraints in Table 1, the $\mathcal{T}urbulence$ - $\mathcal{T}hickened$ flame $\mathcal{R}egime$ validity ranges are obtained in terms of the local Damköhler number, $D{a}_{\Delta }^I$, according to the following equation:

$$P{r}^{1/2}R{e}_{\Delta }^{1/2}\ge D{a}_{\Delta }^I\ge P{r}^{2/7}R{e}_{\Delta }^{2/7}for1\ge Pr\ge R{e}_{\Delta }^{-1}.$$

(25)

4. The Extinction Factor

The wrinkling scales, and in particular the smallest one ${\ell }_w^{*}$ that imposes the highest strain onto the flame front, can stretch a flame front up to local quenching, thus becoming quenching scales.

To have local quenching, the thickness of the flame front has to be smaller than the local integral macroscale. Hence, if ${\delta }_{\mathcal{L}}\ge {\ell }_{\Delta }$ ($\mathcal{V}olumetric\mathcal{R}egime$) no quenching due to eddy flame stretching is possible: no ${\ell }_q^{*}$ exist for $D{a}_{\Delta }^I\le P{r}^{-1}R{e}_{\Delta }^{-1}$ (see Section 5). This is a first boundary of the quenching region in the combustion diagram: this boundary seems to be confirmed by the quenching model described below that compares well with experimental and DNS numerical data.

The Karlovitz number is related to the effect of aerodynamic stretching on the flame front and is defined as the ratio between the chemical time and the stretching time and is presented as follows:

$$K{a}_t=\frac{{\tau }_{ch}}{{\tau }_{st}}\equiv K{a}_{{\eta }_t}=\frac{{\delta }_L/{\mathcal{S}}_{\mathcal{L}}}{\eta /u_{\eta }}=Pr{\left(\frac{{\delta }_L}{\eta }\right)}^2\approx {\left(\frac{{\delta }_L}{\eta }\right)}^2.$$

(26)

By defining the Karlovitz number as in Equation (26), it is assumed that the $\eta$ scales are the most important or effective in stretching a flame. The smallest wrinkling eddy applies the highest strain on the flame front, hence it may be the most effective to produce quenching. Considering that the Karlovitz number defined in Equation (26) can also be expressed as $K{a}_{\eta }=R{e}_{\Delta }^{1/2}/D{a}_{\Delta }^I$, the existence condition for the smallest wrinkling scale in Equation (13) can be rewritten as $K{a}_{\eta }\ge P{r}^{-3/4}$ that is only a bit higher than the Klimov–Williams criterion.

This existence condition of the smallest wrinkling eddy is a necessary but not sufficient condition for quenching. In fact, Peters [6] observed that actual flames resist strain more than predicted by the Klimov–Williams (Karlovitz) criterion.

In the literature, there are some models that provide an extinction or stretch factor ${\mathcal{G}}_{ext}\le 1$ that can be effectively used to “damp” locally the reaction rate, thus mimicking the effect of subgrid flame stretching.

An extinction or stretch factor ${\mathcal{G}}_{ext}$ was firstly introduced by Bray [44] and then adopted in models [26,27,45]. However, quenching predictions based on the Bray’s stretch factor (BSF curves in Figure 2) or on the Klimov–Williams criterion discussed in Section 1 (KWC line) both do not agree with those found from both experiments and direct numerical simulations that provide, in fact, similar trends [22,42,46].

A model that can be more effectively used to predict quenching is the so-called quenching cascade model [46] that compares quite well with experimental and direct numerical simulation data on quenching ([42], pp. 212–214). The main assumptions at the base of the model are:

1

However strong the vortex strain may be, the time required to quench the flame is the same and it is given by ${\tau }_{ch}={\delta }_{\mathcal{F}}/{\mathcal{S}}_{\mathcal{L}}$;

2

Only eddies with large characteristic velocity $u_{\ell }^{\prime }$ (at least $\ge {\mathcal{S}}_{\mathcal{L}}$) are able to quench the flame;

3

If the mean eddy life-time is smaller than the chemical time ${\tau }_{ch}$, such eddies cannot quench the flame;

4

A temporal sequence of approximately ${\tau }_{ch}/{\tau }_{\ell }$ consecutive eddies of size ℓ with large $u_{\ell }^{\prime }$ must exist to induce actual quenching;

5

The mean eddy time-scale ${\tau }_{\ell }$ is estimated considering a multifractal description of turbulence with intermittency;

6

A discrete set of scales is recursively generated.

With these assumptions, the quenching cascade model estimates the total fraction of flame surface undergoing quenching during a time ${\tau }_{ch}$.

Three ${\mathcal{G}}_{ext}$ isolines of the quenching cascade model are shown in Figure 2, where they are compared with the associated isolines of Bray’s stretch factor, the Klimov–Williams criterion $K{a}_{\eta }=1$ and the $K{a}_{\eta }=100$ line suggested by Peters. In particular, the four red lines should be compared. Above them, extinction is predicted at different levels by these models. Since the extinction region predicted by the quenching cascade model has been validated with experimental data [22,42,46], it is assumed as a reference model: hence, high velocity fluctuations are required to produce localized flame quenching. While the Bray’s stretch factor and the Klimov–Williams criterion underestimate extinction, the Peters’ criterion seems to overpredict the flames’ capability of resisting strain.

5. Premixed Combustion Regimes

Thinking about the interaction between a flame front and turbulent eddies, it is possible to identify different combustion regimes through the comparison of the previously defined length scales (${\ell }_{\Delta }$, $\eta$, ${\ell }^{*}$, ${\ell }_w^{*}$, ${\ell }_{tk}^{*}$, ${\delta }_{\mathcal{L}}$) involved in the process.

At first, the local laminar flame front ${\delta }_{\mathcal{L}}$, the turbulent macro-scale ${\ell }_{\Delta }$, and the turbulent dissipative scale $\eta$ are compared. The regime associated with ${\delta }_{\mathcal{L}}\ge {\ell }_{\Delta }$ is named $\mathcal{V}olumetric\mathcal{R}egime$, ${\mathcal{V}}_{\mathcal{R}}$. The regime associated with ${\delta }_{\mathcal{L}}<\eta$ is named $\mathcal{W}rinkled$$\mathcal{R}egime$, ${\mathcal{W}}_{\mathcal{R}}$. These two combustion regimes are experienced, respectively, for low and high Damköhler numbers of first species, as shown in

Table 2. The three main premixed combustion regimes in the LTSM framework based on the comparison between turbulent length scales and laminar flame front thickness.

Regime	Scale Condition		${\mathbf{Da}}_{\Delta }^{\mathit{I}}$ Condition
${\mathcal{V}}_{\mathcal{R}}$	${\delta }_{\mathcal{L}}\ge {\ell }_{\Delta }$	⇒	$D{a}_{\Delta }^I\le P{r}^{-1}R{e}_{\Delta }^{-1}$
${\mathcal{TTC}}_{\mathcal{R}}$	${\ell }_{\Delta }>{\delta }_{\mathcal{L}}\ge \eta$	⇒	$P{r}^{-1}R{e}_{\Delta }^{1/2}\ge D{a}_{\Delta }^I>P{r}^{-1}R{e}_{\Delta }^{-1}$
${\mathcal{W}}_{\mathcal{R}}$	${\delta }_{\mathcal{L}}<\eta$	⇒	$D{a}_{\Delta }^I>P{r}^{-1}R{e}_{\Delta }^{1/2}$

.

Between ${\mathcal{V}}_{\mathcal{R}}$ and ${\mathcal{W}}_{\mathcal{R}}$, there is a more complex regime consisting of other subregimes. In fact, when ${\ell }_{\Delta }>{\delta }_{\mathcal{L}}\ge \eta$, some eddies may not be dissipated by the flame front, some may corrugate it, and in particular conditions of the Prandtl number, some may even thicken it: the role played by these scales is shown in Figure 3. The subregime characterized by ${\ell }_w^{*}>{\ell }^{*}$ is here called $\mathcal{C}orrugated\mathcal{R}egime$, ${\mathcal{C}}_{\mathcal{R}}$. The subregime characterized by ${\ell }_w^{*}={\ell }^{*}$ is named $\mathcal{T}hickened\mathcal{R}egime$, ${\mathcal{T}}_{\mathcal{R}}$. As discussed in Section 3.3, in certain ranges of the Prandtl number, the flame front may be even thickened by some eddies: this subregime is here named $\mathcal{T}urbulence-\mathcal{T}hickened\mathcal{R}egime$, ${\mathcal{TT}}_{\mathcal{R}}$. The intermediate combustion regime in Table 2 is named ${\mathcal{TTC}}_{\mathcal{R}}$ (since it includes both the ${\mathcal{TT}}_{\mathcal{R}}$ and ${\mathcal{C}}_{\mathcal{R}}$ regimes), and it nearly corresponds to the pocket flames region defined by Barrère [2]: the region with $Da\le 1$ corresponds to the $\mathcal{W}$ell-$\mathcal{S}$tirred $\mathcal{R}$eactor; the region including the $\mathcal{T}hickened$ and the $\mathcal{T}urbulence$-$\mathcal{T}hickened$ regimes corresponds to the Distributed Reaction Zone.

The premixed combustion regimes described in Table 2 are described with more details in Figure 4. Depending on the local Prandtl number, there can be four possible conditions: three of them relating to $Pr\le 1$ and only one to $Pr>1$.

The Most Probable Condition for $Pr\le 1$

Considering that common gases have $Pr\in [0.2,1]$ and that a sufficiently high Reynolds number is required (i.e., ${\ell }_{\Delta }/\eta \ge 10$ at least) to experience turbulence/combustion interaction, among the three conditions characterized by $Pr\le 1$, the third one is the most physical. In fact, in the first case, the Prandtl number would be very low even for high Reynolds; in the second case, the Prandtl number would be too low for gaseous combustion at sufficiently high Reynolds numbers. For this reason, and since the condition with $Pr>1$ appears “simpler”, the condition $1\ge Pr\ge R{e}_{\Delta }^{-6/13}$ will be analyzed in more detail in the following. The condition $1\ge Pr\ge R{e}_{\Delta }^{-6/13}$ implies that turbulence/combustion interaction will be assumed to happen for $R{e}_{\Delta }\ge P{r}^{-13/6}$, being $Pr\le 1$. For $R{e}_{\Delta }\le P{r}^{-13/6}$, the flame will be considered locally laminar.

The boundary lines of each regimes are reported in

Table 3. Main boundaries of the LTSM combustion regimes shown in Figure 5.

${\mathcal{V}}_{\mathcal{R}}$−${\mathcal{T}}_{\mathcal{R}}$	$D{a}_{\Delta }^I={\left(PrR{e}_{\Delta }\right)}^{-1}$	${\ell }_{\Delta }/{\delta }_{\mathcal{L}}=1$
${\mathcal{T}}_{\mathcal{R}}$−${\mathcal{TT}}_{\mathcal{R}}$	$D{a}_{\Delta }^I={\left(PrR{e}_{\Delta }\right)}^{2/7}$	$u_{\Delta }^{\prime }/{\mathcal{S}}_{\mathcal{L}}={\left({\ell }_{\Delta }/{\delta }_{\mathcal{L}}\right)}^{5/9}$
${\mathcal{TT}}_{\mathcal{R}}$−${\mathcal{C}}_{\mathcal{R}}$	$D{a}_{\Delta }^I={\left(PrR{e}_{\Delta }\right)}^{1/2}$	$u_{\Delta }^{\prime }/{\mathcal{S}}_{\mathcal{L}}={\left({\ell }_{\Delta }/{\delta }_{\mathcal{L}}\right)}^{1/3}$
${\mathcal{C}}_{\mathcal{R}}$−${\mathcal{W}}_{\mathcal{R}}$	$D{a}_{\Delta }^I=P{r}^{-1}R{e}_{\Delta }^{1/2}$	$u_{\Delta }^{\prime }/{\mathcal{S}}_{\mathcal{L}}=Pr{\left({\ell }_{\Delta }/{\delta }_{\mathcal{L}}\right)}^{1/3}$
${\mathcal{LV}}_{\mathcal{R}}$/${\mathcal{LF}}_{\mathcal{R}}$	$R{e}_{\Delta }\ge P{r}^{-13/6}$	$u_{\Delta }^{\prime }/{\mathcal{S}}_{\mathcal{L}}=P{r}^{-7/6}{\left({\ell }_{\Delta }/{\delta }_{\mathcal{L}}\right)}^{-1}$
$D{a}_{\Delta }^I$ line	$D{a}_{\Delta }^I=\left({\ell }_{\Delta }/u_{\Delta }^{\prime }\right)/\left({\delta }_{\mathcal{L}}/{\mathcal{S}}_{\mathcal{L}}\right)$	$u_{\Delta }^{\prime }/{\mathcal{S}}_{\mathcal{L}}=\left({\ell }_{\Delta }/{\delta }_{\mathcal{L}}\right)D{a}_{\Delta }^{I^{-1}}$
$K{a}_{\eta }$ line	$K{a}_{\eta }=\left({\delta }_{\mathcal{L}}/{\mathcal{S}}_{\mathcal{L}}\right)/\left(\eta /u_{\eta }\right)$	$u_{\Delta }^{\prime }/{\mathcal{S}}_{\mathcal{L}}=P{r}^{1/3}{\left(L/{\delta }_{\mathcal{L}}\right)}^{1/3}K{a}_{{\eta }_t}^{2/3}$

. Such information is used to build the LTSM combustion diagram in Figure 5 for $1\ge Pr\ge R{e}_{\Delta }^{-6/13}$ both in the $u_{\Delta }^{\prime }/{\mathcal{S}}_{\mathcal{L}}vs.{\ell }_{\Delta }/{\delta }_{\mathcal{L}}$ and $D{a}_{\Delta }^{I}vs.PrR{e}_{\Delta }$ logarithmic planes. As a reminder, in the LTSM framework, these diagrams refer to local quantities and not to global ones as in the standard combustion diagram (in fact, the turbulent scales considered are $u_{\Delta }^{\prime }$ and ${\ell }_{\Delta }$ instead of $u_{rms}^{\prime }$ and L).

In the $u_{\Delta }^{\prime }/{\mathcal{S}}_{\mathcal{L}}vs.{\ell }_{\Delta }/{\delta }_{\mathcal{L}}$ plane, there are two boundaries that can change depending on the local Prandtl number: the first is the laminar region boundary, $R{e}_{\Delta }=P{r}^{-13/6}$, that corresponds to the $R{e}_{\Delta }=1$ line of the standard combustion diagram but shifted upwards depending on the $Pr<1$ value; the second is the lower boundary of the $\mathcal{C}orrugated\mathcal{R}egime$ (that is also the top boundary of the $\mathcal{W}rinkled\mathcal{R}egime$), that shifts downwards depending on the $Pr<1$ value (the ${\mathcal{C}}_{\mathcal{R}}$ regime disappears for $Pr=1$). The lower boundary of the $\mathcal{T}urbulence-\mathcal{T}hickened\mathcal{R}egime$ (that is also the upper boundary of the $\mathcal{C}orrugated\mathcal{R}egime$) is the $K{a}_{{\eta }_t}=1$ line with $Pr=1$ (that is the quenching line in the Klimov–Williams criterion).

Figure 6 shows the main scales of turbulence and combustion and the regimes for $1\ge Pr\ge R{e}_{\Delta }^{-6/13}$, i.e., for $R{e}_{\Delta }\ge P{r}^{-13/6}$, with $Pr\le 1$. It is observed that, in this range of Prandtl number, the conditions ${\ell }^{*}={\lambda }_g$ and ${\ell }_w^{*}={\ell }_{w,min}^{*}$ reported in the top frame of Figure 6 cannot be reached since $R{e}_{\Delta }^{-1}<R{e}_{\Delta }^{-6/13}$. Focusing on the eddy length scales and considering that the flame front destroys eddies smaller than ${\ell }^{*}$, it is observed that the turbulent energy cascade is not fully active everywhere: where ${\ell }^{*}>\eta$ the energy cascade stops at scale ${\ell }^{*}$.

6. Modeling the Reacting Volume Fraction

In the Large Eddy Simulation framework, the local turbulent macro-scale ${\ell }_{\Delta }$ will be assumed equal to the characteristic size of the computational cell (filter size), estimated through its volume, i.e., $\Delta ={\mathcal{V}}_{\Delta }^{1/3}$; its associated velocity, $u_{\Delta }^{\prime }$, will be assumed equal to the subgrid velocity derived through the specific eddy viscosity model adopted in the simulation. In Section 2.2, the strategy adopted for turbulent combustion closure was introduced: the Favre filtered chemical source term in the energy and single species transport equations is modeled as $\widetilde{{\omega }_i}\approx {\gamma }^{*}{\omega }_i^{*}$, where the local reacting volume fraction of the computational cell is modeled as in Equation (6).

In Equation (6), the laminar flame speed is estimated as ${\mathcal{S}}_{\mathcal{L}}\approx {\left(\alpha /{\tau }_{ch}\right)}^{1/2}$. The turbulent flame speed and the flame front thickness will be modeled depending on the combustion regimes. In particular, the laminar flame thickness will be estimated as ${\delta }_{\mathcal{L}}\approx {\left(\alpha {\tau }_{ch}\right)}^{1/2}$ but not in the $\mathcal{T}urbulence$-$\mathcal{T}hickened$$\mathcal{R}egime$ where its thickened counterpart defined in Equation (19) will be assumed. The local filtered chemical time required to estimate laminar quantities can be calculated as ${\tau }_{ch}=\rho C_{p}T/\left|\Delta {\mathcal{H}}_R\right|$, where $\Delta {\mathcal{H}}_R={\sum }_{i=1}^{N_s}{\mathcal{H}}_i\dot{{\omega }_i}$ is the heat of reaction, $N_s$ being the number of chemical species, ${\mathcal{H}}_i=h_{f_i}^0\left(T_r\right)+\Delta h_{s_i}\left(T\right)$ the enthalpy of the i-th chemical species, $h_{f_i}^0\left(T_r\right)$ its formation enthalpy at the reference temperature $T_r=298.15\mathrm{K}$, $\Delta h_{s_i}={\int }_{T_r}^T{C}_{p_i}\left(T\right)dT$ its sensible enthalpy, and $\dot{{\omega }_i}$ its reaction rate.

In conclusion, the local $Pr$, $R{e}_{\Delta }$, $D{a}_{\Delta }^I$ are first calculated. Based on these characteristic numbers, subgrid turbulent combustion is taken into account by modeling ${\gamma }^{*}$ by means of the expressions summarized in

Table 4. Summary of the LTSM combustion regimes shown in Figure 5, their existence conditions, and the related expressions for the reacting fine structures volume ${\gamma }^{*}$.

${\mathcal{V}}_{\mathcal{R}}$	$\begin{array}{c}D{a}_{\Delta }^I\le {\left(PrR{e}_{\Delta }\right)}^{-1}\\ \forall R{e}_{\Delta }>1andPr>0\end{array}$	${\mathit{\gamma }}^{*}=1$
${\mathcal{T}}_{\mathcal{R}}$	$\begin{array}{c}{\left(PrR{e}_{\Delta }\right)}^{2/7}\ge D{a}_{\Delta }^I\ge {\left(PrR{e}_{\Delta }\right)}^{-1}\\ ifR{e}_{\Delta }\ge P{r}^{-13/6}andPr\le 1\\ or\\ P{r}^{-1}R{e}_{\Delta }^{1/2}\ge D{a}_{\Delta }^I\ge {\left(PrR{e}_{\Delta }\right)}^{-1}\\ ifR{e}_{\Delta }\ge 1andPr>1\end{array}$	$\begin{array}{c}{\gamma }^{*}=\frac{{\delta }_{\mathcal{L}}}{\Delta }={\left(PrR{e}_{\Delta }D{a}_{\Delta }^I\right)}^{-1/2}\le 1\hfill \\ {\left.{\gamma }^{*}\right|}_{max}=1\hfill \\ {\left.{\gamma }^{*}\right|}_{min}={\left(PrR{e}_{\Delta }\right)}^{-9/14}\le 1\hfill \end{array}$
${\mathcal{TT}}_{\mathcal{R}}$	$\begin{array}{c}{\left(PrR{e}_{\Delta }\right)}^{1/2}\ge D{a}_{\Delta }^I\ge {\left(PrR{e}_{\Delta }\right)}^{2/7}\\ ifR{e}_{\Delta }\ge P{r}^{-13/6}andPr\le 1\end{array}$	$\begin{array}{c}{\gamma }^{*}={\left.\frac{{\mathcal{S}}_{\mathcal{T}}}{{\mathcal{S}}_{\mathcal{L}}}\right|}_{\mathcal{Z}}\frac{{\delta }_{\mathcal{F}}^{*}}{\Delta }={\mathcal{A}}_{\mathcal{Z}}{\left(PrR{e}_{\Delta }D{a}_{\Delta }^{I^{-7/2}}\right)}^{1/2}<1\hfill \\ {\left.{\gamma }^{*}\right|}_{max}={\mathcal{A}}_{\mathcal{Z}}<1\hfill \\ {\left.{\gamma }^{*}\right|}_{min}={\mathcal{A}}_{\mathcal{Z}}{\left(PrR{e}_{\Delta }\right)}^{-3/8}<1\hfill \end{array}$
${\mathcal{C}}_{\mathcal{R}}$	$\begin{array}{c}P{r}^{-1}R{e}_{\Delta }^{1/2}\ge D{a}_{\Delta }^I\ge {\left(PrR{e}_{\Delta }\right)}^{1/2}\\ ifR{e}_{\Delta }\ge P{r}^{-13/6}andPr\le 1\end{array}$	$\begin{array}{c}{\gamma }^{*}={\left.\frac{{\mathcal{S}}_{\mathcal{T}}}{{\mathcal{S}}_{\mathcal{L}}}\right|}_{\mathcal{C}}\frac{{\delta }_{\mathcal{L}}}{\Delta }\approx {\left.\frac{{\mathcal{S}}_{\mathcal{T}}}{{\mathcal{S}}_{\mathcal{L}}}\right|}_{\mathcal{Z}}\frac{{\delta }_{\mathcal{L}}}{\Delta }={\mathcal{A}}_{\mathcal{Z}}D{a}_{\Delta }^{I^{-3/4}}<1\hfill \\ {\left.{\gamma }^{*}\right|}_{max}={\mathcal{A}}_{\mathcal{Z}}{\left(PrR{e}_{\Delta }\right)}^{-3/8}<1\hfill \\ {\left.{\gamma }^{*}\right|}_{min}={\mathcal{A}}_{\mathcal{Z}}{\left(P{r}^{-2}R{e}_{\Delta }\right)}^{-3/8}<1\hfill \end{array}$
${\mathcal{LV}}_{\mathcal{R}}$	$\begin{array}{c}{\delta }_{\mathcal{L}}\ge \Delta \\ ifR{e}_{\Delta }<P{r}^{-13/6}andPr\le 1\\ orifR{e}_{\Delta }\le 1andPr>1\end{array}$	$\begin{array}{c}{\gamma }^{*}=1\hfill \end{array}$
${\mathcal{LF}}_{\mathcal{R}}$	$\begin{array}{c}{\delta }_{\mathcal{L}}<\Delta \\ ifR{e}_{\Delta }<P{r}^{-13/6}andPr\le 1\\ orifR{e}_{\Delta }\le 1andPr>1\end{array}$	$\begin{array}{c}{\gamma }^{*}=\frac{{\delta }_{\mathcal{L}}}{\Delta }<1\hfill \end{array}$

. These expressions guarantee ${\gamma }^{*}\le 1$; the identity $\Delta \equiv {\ell }_{\Delta }$ was assumed in deriving their non-dimensional number dependence. It is observed that the Zimont expression for ${\mathcal{S}}_{\mathcal{T}}$ has been assumed not only in the $\mathcal{T}urbulence$-$\mathcal{T}hickened$$\mathcal{R}egime$ but also in the $\mathcal{C}orrugated$$\mathcal{R}egime$, due to the lack of reliable experimental data [43] and the superiority of the Zimont model with respect to other models [43]. In the $\mathcal{W}rinkled$$\mathcal{R}egime$, subgrid hydrodynamic effects and Lewis number effects on ${\mathcal{S}}_{\mathcal{T}}$ have been neglected for the time being. However, it is observed that this regime appears to be of minor importance in industrial applications [43]. Work is currently going on by the present authors to define, in this low strain regime (Markstein regime), a $\mathcal{X}$ factor taking into account thermo-diffusive effects (induced by reactants’ Lewis number) that increase or decrease the turbulent flame speed and flame wrinkling: ${\mathcal{S}}_{\mathcal{T}}/{\mathcal{S}}_{\mathcal{L}}=1-\mathcal{X}$, with $\mathcal{X}>0$ or $\mathcal{X}<0$ depending on local curvature, strain, and Markstein number sign.

7. Model Validation

The Localised Turbulent Scale Model is validated here against data coming from a Direct Numerical Simulation performed by these authors [47]. In particular, the time average and rms fluctuation of some quantities resulted from the LES simulation were compared with their DNS counterpart.

7.1. The Test Case

The test case consists of an unconfined and atmospheric Bunsen flame developing along the streamwise direction (z). This premixed flame is produced by three adjacent rectangular slot burners whose size is undefined in the spanwise direction (x) (due to the periodic boundary conditions imposed in the simulation) and that are separated along the transversal direction (y) by means of two $0.17\mathrm{mm}$ thick walls. The central slot burner injects a fresh mixture of methane, hydrogen, and air, while the two side burners inject hot combustion products of the same central mixture.

The reactant mixture, with an equivalence ratio $\mathsf{\Phi }=0.7$ and with 0.2 mole fraction of hydrogen, is injected from the central slot with a bulk velocity of $100\mathrm{m}{\mathrm{s}}^{-1}$ and at $600\mathrm{K}$. The velocity of the coflow stream is $25\mathrm{m}{\mathrm{s}}^{-1}$. The central jet Reynolds number is 2253, based on the width of the jet, $1.2\mathrm{mm}$, its bulk velocity, and the kinematic viscosity $5.327·{10}^{-5}{\mathrm{m}}^2{\mathrm{s}}^{-1}$. Homogeneous isotropic turbulence is forced at the inlet. Such turbulence is artificially produced by means of a synthetic turbulence generator implemented from [48]. In particular, the spatial correlation length scales and velocity fluctuations provided as input to this generator are: $L_{zz}=0.8\mathrm{mm}$, $L_{xx}=L_{yy}=L_{zz}/2=0.4\mathrm{mm}$, $u_z^{\prime }=u_x^{\prime }=u_y^{\prime }=12\mathrm{m}{\mathrm{s}}^{-1}$ with no shear stresses (the Reynolds stress tensor is diagonal) [49].

7.2. The Numerical Set-Up

The computational domain is three-dimensional with 61 nodes in the spanwise direction (x), extending from −1.5 to $1.5\mathrm{mm}$, and along which periodic boundary conditions are forced. The computational domain has four zones: a central injection zone extending from $-4\mathrm{mm}$ to 0 in the streamwise direction (z) with a width (along y) of $1.2\mathrm{mm}$ and with $60\times 35$ nodes ($zy$); two surrounding zones extending from $-0.4\mathrm{mm}$ to 0 in the streamwise direction and from the central slot external wall (wall thickness $0.17\mathrm{mm}$) up to $18\mathrm{mm}$ outward in the transversal direction (y) with $8\times 89$ nodes ($zy$); a main mixing and reacting zone downstream of the injection extending from 0 to $24\mathrm{mm}$ in the streamwise direction and from $-18\mathrm{mm}$ to $18\mathrm{mm}$ in the transversal direction with $241\times 221$ nodes ($zy$). The whole computational domain has 56,785 × 61 = 3,463,885 nodes.

Since the LES results are going to be compared against DNS numerical data, it is useful to quantify how much the LES computational grid is coarser than the DNS grid. Here, we simply report the main quantities of the two grids: the $\Delta y^{LES}/\Delta y^{DNS}$ spacing ratio is in the range 1.42–6.25, and the $\Delta z^{LES}/\Delta z^{DNS}$ is in the range 3.14–8.21, with the higher values reached in the flame front region; the ratio between the aspect ratios, $(\Delta y^{LES}/\Delta z^{LES})/(\Delta y^{DNS}/\Delta z^{DNS})$, is in the range of 0.17–2.00.

The simulation was performed using the same in-house HeaRT code previously used for the DNS, which solves the fully compressible reactive Navier–Stokes equations in their conservative formulation. The fully explicit third-order accurate TVD low-storage Runge–Kutta scheme of Shu and Osher [50] is used for temporal integration. A staggered formulation of the transported quantities is adopted for the spatial integration: a second-order centered scheme is used for all diffusive fluxes and momentum convective terms, while scalar convective terms are modeled through the AUSM method coupled to a third-order QUICK interpolation. No artificial damping or filtering was necessary for ensuring a stable solution.

The ideal gas equation of state is assumed. The same models adopted to estimate molecular transport properties in the DNS calculations are here adopted. The Hirschfelder and Curtiss [51] approximate formula is adopted to model molecular transport in a multicomponent mixture, taking into account preferential diffusion effects. The Soret thermo-diffusive effect is also considered. Molecular transport properties for individual species, except their binary mass diffusivities, are accurately calculated through kinetic theory by using a priori the software library provided by A. Ern (EGlib) [52,53]. The calculated values are stored in a look-up table from 200 to 5000 K every 100 K. Values at intermediate temperatures are calculated at run-time by linear interpolation and then used to derive mixture-average properties through the Wilke’s formula with Bird’s correction for viscosity ([54], pp. 23–29) and Mathur’s expression for thermal conductivity ([54], pp. 274–278). Binary mass diffusion coefficients are calculated by means of kinetic theory expressions ([54], pp. 525–528) at run-time. Thermo-diffusion coefficients required to model the Soret effect are calculated at run-time by using the EGlib routines. The chemical mechanism adopted is a skeletal mechanism having 17 transported species and 58 reactions; in particular, the transported species are: ${\mathrm{CH}}_4$, ${\mathrm{H}}_2$, ${\mathrm{O}}_2$, ${\mathrm{N}}_2$, $\mathrm{OH}$, $\mathrm{O}$, $\mathrm{H}$, ${\mathrm{HO}}_2$, ${\mathrm{H}}_2{\mathrm{O}}_2$, ${\mathrm{CH}}_3$, ${\mathrm{CH}}_2$, $\mathrm{CH}$, ${\mathrm{CH}}_2\mathrm{O}$, $\mathrm{HCO}$, $\mathrm{CO}$, ${\mathrm{CO}}_2$, ${\mathrm{H}}_2\mathrm{O}$ [55].

The walls are assumed viscous and adiabatic. Velocity, temperature, and mass fractions are fixed at the inlet, with turbulent velocity fluctuations artificially generated [48]. Periodic boundary conditions are applied in the spanwise direction, while improved partially non-reflecting outflow boundary conditions are imposed at the other open boundaries of the computational domain [56,57]. The subgrid scale model adopted for the turbulence closure is the dynamic Smagorinsky model [58]. This is used not only to calculate the eddy viscosity ${\mu }_t$, but also the instantaneous and local subgrid velocity fluctuation required in the LTSM model, i.e., $u_{\Delta }^{\prime }={\mu }_t/\left(C_S\rho \Delta \right)$ where $C_S$ is the Smagorinsky constant calculated dynamically.

7.3. The Results

The actual velocity fluctuation at the end of the injection channel is $u^{\prime }\approx 12\mathrm{m}{\mathrm{s}}^{-1}$, and the turbulent length scale is $L_t\approx 1.05\mathrm{mm}$. Using these quantities and the kinematic viscosity of the fresh mixture, the central jet turbulent Reynolds number is 236, and the Kolmogorov length scale is $\eta \approx 17.42\mathsf{\mu }$m. The adiabatic flame temperature is $2072\mathrm{K}$. The laminar flame speed and flame front thickness (estimated from kinetics simulations including the Soret effect) at these conditions are ${\mathcal{S}}_{\mathcal{L}}=1.01\mathrm{m}{\mathrm{s}}^{-1}$ and ${\delta }_{\mathcal{F}}=0.378\mathrm{mm}$, respectively. Hence, $u_{rms}^{\prime }/{\mathcal{S}}_{\mathcal{L}}=11.88$ and $L_t/{\delta }_{\mathcal{F}}=2.78$; $K{a}_{\eta }=471$, $D{a}_{L_t}^I=0.24$, thus ideally locating this flame into the broken reaction regime of the standard combustion diagram, well above the Peter’s boundary $K{a}_{\eta }=100$ of the thin reaction zone and where turbulence is expected to strongly influence premixed flame structure.

However, such a conclusion based on this commonly accepted procedure is wrong. In fact, the distribution of reacting cells reported for an instantaneous DNS flowfield in Figure 7. To build this diagram, a progress variable was defined as the normalized sum of ${\mathrm{CO}}_2$, $\mathrm{CO}$, ${\mathrm{H}}_2\mathrm{O}$, and ${\mathrm{H}}_2$ mass fractions. Then, the iso-surface $c=0.1$ (in the fresh mixture) was selected; on this surface, the auto-correlation length scales ${\ell }_{\Delta }$ and rms velocity fluctuations $u_{\Delta }^{\prime }$ were calculated: these quantities characterize the local eddies interacting with the local flame front. On the same surface, the local laminar flame speed ${\mathcal{S}}_{\mathcal{L}}$ was also estimated through the method developed in [24]; the local flame front thickness ${\delta }_{\mathcal{L}}$ was estimated as the curvilinear length of the curvilinear segment normal to the iso-c surfaces starting from $c=0.1$ and ending at $c=0.99$. shows that the flame experiences all the combustion regimes depicted within the LTSM framework. Most of the flame belongs to the thin reaction zone and only the reacting regions close to the flame tip ($z>1.5\mathrm{cm}$) enter the quenching region with an increased probability of localized extinctions for the farthest regions from the inlet.

The comparison between LES and DNS data are very good in terms of averaged data, and good in terms of rms fluctuations. In fact, the average streamwise velocity and temperature flowfields look very similar. For brevity reasons, only a brief analysis comparing transversal and axial distributions of main quantities is provided in the following; more results are provided in the Supplementary Materials.

In the LES case, the temporal statistics were calculated by means of 22,033 transversal planes, sampled in $4.4\mathrm{ms}$. Since the DNS is more computationally expensive, it was run for an admissible CPU time collecting 152 time-samples of the complete three-dimensional flowfields in $1.73\mathrm{ms}$; to enhance the convergence of the otherwise poor temporal statistics, every single plane (140) in the spanwise direction of the 152 3D-samples was assumed as an independent event, thus increasing the number of samples to 21,280. Inspecting contour maps of rms velocity fluctuations (not show here) built according to this strategy, it was observed that the profiles of the DNS data exhibit less symmetry with respect to the central jet axis than LES data, and this was more visible moving downstream. This was attributed to the lack of a fully converged statistics of the DNS rms data: since the tip of the flame is the most unstable region of the flow, exhibiting strong transversal oscillations and localized extinctions, more samples would be required. Hence, since the flow is topologically symmetric with respect to the central axis of the jet, the two transversal half-planes were considered as independent events in both LES and DNS, thus doubling the number of samples. However, these statistics are not still fully converged for the whole transversal plane, especially for rms fluctuations far downstream from injection. This is clear comparing the LES and DNS average and rms fluctuation of the streamwise velocity and temperature distributions along the central jet axis (see Supplementary Materials). This suggests to limit the comparison of transversal profiles up to 6 h, i.e., $7.50\mathrm{mm}$, above the injection: the heights 2 h and 5 h are reported in the paper, while the others are available in the Supplementary Materials.

Looking at streamwise velocity profiles in Figure 8 (top), it is observed that the average turbulent flat profiles quickly evolves into the one characteristic of a free jet, with peak velocity fluctuation located in the two shear-layers. The agreement between LES and DNS average data are very good, while that of rms fluctuations is good, showing lower LES peaks in the shear-layers up to 4 h and an LES profile translated towards a bit higher values at 6 h. Similar comments can be drawn for the transversal velocity profiles in Figure 8 (bottom) that, in general, exhibit a worse agreement.

Looking at temperature profiles in Figure 9 (top), it is observed that the agreement between averaged LES and DNS data are also very good, with a maximum error of ∼50 K over temperatures of the order of $1000\mathrm{K}$ along the central axis in the most downstream sections. The location of rms temperature peaks is well predicted, but the amplitude of the LES peaks exhibits an error of ∼70 K over temperatures of the order of $600\mathrm{K}$ in the middle sections, 3 h and 4 h. Along the central axis, the rms agreement deteriorates starting from the 4 h section and increases moving more downstream.

Mass fraction of main species are then compared in terms of their time averages at the same locations in Figure 9 (bottom), showing a good agreement too. It is observed that the ${\mathrm{H}}_2\mathrm{O}$ profiles are not reported since they are very similar to the ${\mathrm{CO}}_2$ profiles, both in shape and value. The worst agreement is achieved for the OH radical species.

It is observed that both temperature and velocity (especially the transversal one) profiles at the first two heights reveal that the dynamics of the flow downstream of the two walls separating the central jet from the side hot jets is not well captured by the LES simulation. This disagreement can be attributed to the poor spatial resolution of the LES grid to describe the thin wake regions of the two separating walls. Despite this, it is concluded that the transversal spreading of the flame is well reproduced at different heights from the injection, both in terms of velocity and temperature distributions. This implies that the coupling between the adopted turbulence and combustion models, i.e., the dynamic Smagorinsky and LTSM models, is effective in capturing the correct behavior of the present jet flame.

8. Discussion

This work provided a brief overview on how combustion regimes are introduced in standard and spectral combustion diagrams. Evidence that different regimes are experienced locally by the same turbulent flame has been given, reporting actual combustion diagrams for different DNS flame simulations. The actual combustion diagrams derived by DNS data are built on local quantities: auto-correlation length scales, rms velocity fluctuations, laminar flame speed, and flame front thickness. In particular, for one of these DNS cases, it is shown that the actual turbulent flame not only experiences several combustion regimes, but also that the global nominal combustion regime predicted according to the common methodology based on the global Karlovitz and Damköhler numbers is completely wrong.

Hence, a new perspective to look at combustion diagrams has been suggested. The analysis adopted consists of the identification of the main scales involved in turbulent combustion physics, in defining their role and their possible interaction. The result of this reckoning is a model describing combustion regimes in terms of turbulence and combustion scale interaction: this theoretical framework has been named LTSM, Localized Turbulence Scales Model. Adopting such a model as a subgrid scale model for Large Eddy Simulation implies assuming the local LES filter size as the local macro-scale ${\ell }_{\Delta }$, and the subgrid velocity fluctuation derived by the specific eddy viscosity model adopted as its associated local velocity scale.

Based on this framework, a closure for the filtered chemical source term in the context of Large Eddy Simulation has been suggested. The turbulent combustion subgrid scale model consists of modeling the local reacting volume fraction ${\gamma }^{*}$ as in the Eddy Dissipation Concept, but the model depends on the combustion regime locally experienced. More specifically, ${\gamma }^{*}$ depends on the turbulent and laminar flame speeds ratio, and on the flame thickness; in terms of non-dimensional numbers, it depends on the Prandtl, Reynolds, and first species Damköhler numbers.

Finally, the suggested subgrid model has been validated by comparing LES and DNS data related to an unconfined and atmospheric slot burner flame with two adjacent hot coflows. Average, rms velocity and temperature fluctuations, and average mass fractions (${\mathrm{CH}}_4$, ${\mathrm{H}}_2$, ${\mathrm{CO}}_2$, $\mathrm{CO}$, $\mathrm{OH}$) compare very well against DNS data. This, coupled with the interesting LES/DNS grid ratios adopted, confirms the goodness of the suggested subgrid model for the filtered chemical source term.

9. Conclusions

This work is an attempt to respond to the lack of a unified subgrid turbulent combustion model for Large Eddy Simulation, as shown from the literature review reported in Section 1, and this represents the novelty and value of the contribution of the paper. The Localised Turbulent Scales Model (LTSM) with its closure of filtered chemical source terms is the suggested model, able to self-adapt to the local combustion regimes taking into account the different conditions of the flow (turbulence, thermodynamics, composition, diffusive transport, acoustics). Although its reported validation against DNS data reveals it as a promising subgrid model, further comparison against experimental data and other combustion models would contribute to its robustness, applicability, and reliability; this will be the objective of future works.