A Numerical Investigation of Mixing Models in LES-FMDF for Compressible Reactive Flows

Abstract

The filtered mass density function (FMDF) model has been employed for large-eddy simulations (LES) of compressible high-speed turbulent mixing and reacting flows. However, the mixing model remains a pressing challenge for FMDF methods, especially for compressible reactive flows. In this work, a temporal development mixing layer with two different convective Mach numbers, ${\mathrm{M}}_{\mathrm{c}}=0.4$ and ${\mathrm{M}}_{\mathrm{c}}=0.8$, is used to investigate the mixing models. A simplified one-step reaction and a real hydrogen/air reaction are employed to study the mixing and turbulence-chemistry interaction. Two widely used mixing models, interaction by exchange with the mean (IEM) and Euclidean minimum spanning tree (EMST), are studied. Numerical results indicate that no difference is observed between the IEM and EMST models in simple reaction flows. However, for hydrogen/air reactions, the EMST model can predict the reaction more accurately in high-speed flow. For mixing models in compressible reactive flows, the requirement of localness preservation tends to be more essential as the convective Mach number increases. With the increase of compressibility, the sensitivity of the mixing model coefficient is reduced significantly. Therefore, the appropriate mixing model coefficient has a wider range. Results also indicate that a large error may result when using a fixed mixing model coefficient in compressible flows.

1. Introduction

The complex high-speed turbulent combustion process needs to be better understood to further develop the scramjet engine [1,2,3,4]. However, the accurate prediction of combustion in high-speed compressible flows remains a challenge due to the turbulence-chemistry interaction (TCI). On one hand, turbulence can promote the mixing of oxygen and fuel, increase the flame surface area, and strengthen the chemical reaction. On the other hand, the fluctuation disturbance caused by turbulence may also extinguish the local combustion flame. In addition, in the case of high Mach numbers, the stronger compression effect and even shock further complicate the TCI. Therefore, an accurate estimation of the TCI is important for the improvement of the simulation of turbulence combustion [5].

The transportable probability density function (TPDF) model can accurately close the source term of the chemical reaction with no assumptions of a low-dimensional flow pattern. Its ability to calculate the complex TCI has been demonstrated in many previous studies [6,7,8]. There are limited studies, however, on the application of such methods in high-speed compressible flows. At present, the high-speed source term, as well as mixing model, are major challenges in compressible high-speed reactive flows. The former attempts to model the effects of pressure and viscous dissipation on energy equation, which cannot be ignored in high-speed flows [9,10,11]. The latter is to close the conditional average/filtering of the molecular diffusion term. It is also the most important and the weakest link of the TPDF method in simulating both low- and high-speed turbulent reaction flow [12]. The model responsible for simulating the molecular diffusion process is called the mixing model. At present, the most widely used mixing models in the field of low-speed flow include the interaction by exchange with the mean (IEM) [13], modified curl (MC) [14], Euclidean minimum spanning tree (EMST) [15], multi-mapping closed (MMC) [16], parameterized scalar profile (PSP) [17], and shadow-position mixing model (SPMM) [18].

In previous studies, the mixing model has been widely investigated for incompressible reactive flows [19]. Mitarai et al. [20] investigated the performance of three mixing models, namely IEM, MC, and EMST, which are implemented in both Reynolds-averaged simulation and large-eddy simulation (LES). The results obtained from the simulation of a fully developed isotropic decaying turbulence showed that the EMST mixing model achieved significantly better results than the IEM and MC mixing models. Renfeng et al. [21] also studied these three mixing models in a series of non-premixed piloted methane/air jet flames using the PDF/ISAT with the GRI3.0 mechanism and determined that the burning index (BI) depended sensitively on the chemical mechanism used, as well as on the mixing model coefficient ${\mathrm{C}}_{\mathsf{\varphi }}$. Here, for the ${\mathrm{C}}_{\mathsf{\varphi }}$ values, the mean and RMS mixture fraction calculated by EMST were in good agreement with the experimental data, which was not observed in the case of the IEM and MC models. It also suggested that the relatively good predictions obtained using the IEM or MC models may be associated with the adjustment of mixing model coefficients. In addition, some numerical results [20] show that compared with that of RANS/PDF, the effect of mixing models in LES/FDF was not obvious.

In recent years, Zhou et al. [22] used a RANS-PDF method to study piloted premixed jet burner. Two flames with different combustion regimes and two mixing models, EMST and IEM, were discussed in detail. It was observed that the localness requirement of mixing in the composition space was essential for combustion in the flamelet regime, whereas it was dispensable for combustion in the broken reaction zone regime. Krisman et al. [23] used direct numerical simulation (DNS) datasets to provide both the initial conditions and inputs needed for the three mixing models (IEM, MC, and EMST) of the runs, allowing the study to focus on the mixing model. From these studies, it can be observed that the EMST model generally had an advantage in predicting extinction and reignition in a one-dimensional non-premixed jet flame, and the results suggested that the optimal value ${\mathrm{C}}_{\mathsf{\varphi }}$ was case-dependent.

Although there have been several studies to help us understand the performance of mixing models [24] and their applicable conditions, it has not yet been determined whether the characteristics and conclusions of these models can be applied to high-speed compressible flows. Firstly, it is an urgent need to estimate the impact of the mixing model coefficient ${\mathrm{C}}_{\mathsf{\varphi }}$ when using the LES/FDF method to calculate the high-speed reaction flow. The assumption ${\mathrm{C}}_{\mathsf{\varphi }}$ = 2.0 is valid only when both the scalar and velocity spectra are at equilibrium, which is widely used for RANS/PDF. However, in LES, the turbulence timescale is based on the sub-grid quantities, and a wide range of values has been used in previous LES/FDF [25,26] calculations. However, in high-speed flow, there is still a lack of further understanding of the value range and sensitivity of the mixing model coefficient. Secondly, there is also an immediate need to estimate the impact of the different mixing models on high-speed flow. Previous studies related to the LES/FDF method widely used the IEM model [10,11,21] and obtained acceptable simulation results. As for EMST, studies [22,27] of high-speed premixed flame only observed that the localness requirement of mixing in the composition space was essential for combustion in the flamelet regime. Therefore, it is necessary to explore the performance and applicability of different mixing models in high-speed flows.

In this work, a large eddy simulation-filtered mass density function (LES-FMDF) method [28,29] is used, which is a density-based method that considers a high-speed source term. The importance of the mixing formulation is illustrated by comparing the predictions from the IEM and EMST models in the compressible temporally developing reactive mixing layers.

The objective of this study is to investigate the performance of two mixing models in a compressible reactive flow. This work has been organized into the following contents: first, the transport equation of the FMDF and its solution using the Lagrangian particle/Eulerian mesh (LPEM) method are briefly discussed in Section 2. A brief introduction to two mixing models, IEM and EMST, is also included. To determine the interaction between turbulence and chemistry, a simple unreleased heat reaction was introduced to study the scalar mixing performance of different mixing models. Next, a real hydrogen/air reaction is used to study the impact of different mixing models on a strong TCI condition. Moreover, the effect of the compressible and model parameters is discussed in the following section.

2. Materials and Methods

2.1. Lagrangian PDF Approach

In the hybrid LES/FMDF problem, a Lagrangian particle/Eulerian mesh (LPEM) methodology is used to solve the velocity, pressure, and scalar fields. The filtered equations in the LES are as follows:

$$\frac{\mathsf{\partial }{\langle \mathsf{\rho }\rangle }_{\mathrm{l}}}{\mathsf{\partial }\mathrm{t}}+\frac{\mathsf{\partial }{\langle \mathsf{\rho }\rangle }_{\mathrm{l}}{\langle {\mathrm{u}}_{\mathrm{j}}\rangle }_{\mathrm{L}}}{\mathsf{\partial }{\mathrm{x}}_{\mathrm{j}}}=0,$$

(1)

$$\frac{\mathsf{\partial }{\langle \mathsf{\rho }\rangle }_{\mathrm{l}}{\langle {\mathrm{u}}_{\mathrm{i}}\rangle }_{\mathrm{L}}}{\mathsf{\partial }\mathrm{t}}+\frac{\mathsf{\partial }}{\mathsf{\partial }{\mathrm{x}}_{\mathrm{j}}}\left[{\langle \mathsf{\rho }\rangle }_{\mathrm{l}}{\langle {\mathrm{u}}_{\mathrm{j}}\rangle }_{\mathrm{L}}{\langle {\mathrm{u}}_{\mathrm{i}}\rangle }_{\mathrm{L}}+{\langle \mathrm{p}\rangle }_{\mathrm{l}}{\mathsf{\delta }}_{\mathrm{i}\mathrm{j}}-{\langle {\mathsf{\tau }}_{\mathrm{i}\mathrm{j}}\rangle }_{\mathrm{L}}+{\mathsf{\tau }}_{\mathrm{i}\mathrm{j}}^{\mathrm{s}\mathrm{g}\mathrm{s}}\right]=0,$$

(2)

$$\frac{\mathsf{\partial }{\langle \mathsf{\rho }\rangle }_{\mathrm{l}}{\langle {\mathrm{e}}_{\mathrm{t}}\rangle }_{\mathrm{L}}}{\mathsf{\partial }\mathrm{t}}+\frac{\mathsf{\partial }}{\mathsf{\partial }{\mathrm{x}}_{\mathrm{j}}}\left[{\langle {\mathrm{u}}_{\mathrm{j}}\rangle }_{\mathrm{L}}\left({\langle \mathsf{\rho }\rangle }_{\mathrm{l}}{\langle {\mathrm{e}}_{\mathrm{t}}\rangle }_{\mathrm{L}}+{\langle \mathrm{p}\rangle }_{\mathrm{l}}\right)+{\langle {\mathrm{q}}_{\mathrm{i}}\rangle }_{\mathrm{l}}-{\langle {\mathrm{u}}_{\mathrm{i}}\rangle }_{\mathrm{L}}{\langle {\mathsf{\tau }}_{\mathrm{i}\mathrm{j}}\rangle }_{\mathrm{L}}+{\mathrm{H}}_{\mathrm{j}}^{\mathrm{s}\mathrm{g}\mathrm{s}}+{\mathsf{\sigma }}_{\mathrm{j}}^{\mathrm{s}\mathrm{g}\mathrm{s}}\right]=0,$$

(3)

where the density ${\langle \mathsf{\rho }\rangle }_{\mathrm{l}}$, Favre-filtered velocity vector ${\langle {\mathrm{u}}_{\mathrm{i}}\rangle }_{\mathrm{L}}$, pressure ${\langle \mathrm{p}\rangle }_{\mathrm{l}}$, and total energy ${\langle {\mathrm{e}}_{\mathrm{t}}\rangle }_{\mathrm{L}}$ comprise the four primary transport variables. In addition, ${\mathrm{q}}_{\mathrm{i}}$ is the heat flux, ${\mathsf{\delta }}_{\mathrm{i}\mathrm{j}}$ is the Kronecker delta function, ${\mathsf{\tau }}_{\mathrm{i}\mathrm{j}}$ the viscous stress tensor, and ${\mathsf{\tau }}_{\mathrm{i}\mathrm{j}}{}^{\mathrm{s}\mathrm{g}\mathrm{s}}$, ${\mathrm{H}}_{\mathrm{i}}{}^{\mathrm{s}\mathrm{g}\mathrm{s}}$, and ${\mathsf{\delta }}_{\mathrm{i}}{}^{\mathrm{s}\mathrm{g}\mathrm{s}}$ denote the sub-grid scale (SGS) stress, SGS enthalpy, and SGS viscosity, respectively. Typically, the SGS stress is closed by the turbulent viscosity assumption:

$${\mathsf{\tau }}_{\mathrm{i}\mathrm{j}}^{\mathrm{s}\mathrm{g}\mathrm{s}}=-2{\langle \mathsf{\rho }\rangle }_{\mathrm{l}}{\mathrm{v}}_{\mathrm{t}}\left({\langle {\mathrm{S}}_{\mathrm{i}\mathrm{j}}\rangle }_{\mathrm{L}}-\frac{1}{3}{\langle {\mathrm{S}}_{\mathrm{k}\mathrm{k}}\rangle }_{\mathrm{L}}{\mathsf{\delta }}_{\mathrm{i}\mathrm{j}}\right)+\frac{2}{3}{\langle \mathsf{\rho }\rangle }_{\mathrm{l}}{\mathrm{k}}^{\mathrm{s}\mathrm{g}\mathrm{s}}{\mathsf{\delta }}_{\mathrm{i}\mathrm{j}},\mathsf{\cdots }\mathsf{\cdots }{\mathrm{v}}_{\mathrm{t}}={\mathrm{C}}_{\mathsf{\mu }}\Delta \sqrt{{\mathrm{k}}^{\mathrm{s}\mathrm{g}\mathrm{s}}},$$

(4)

where ${\mathsf{\upsilon }}_{\mathrm{t}}$ is the turbulent viscosity. The SGS turbulent kinetic energy ${\mathrm{k}}^{\mathrm{s}\mathrm{g}\mathrm{s}}$ can be obtained from the SGS turbulent kinetic energy transport equation proposed by Yoshizawa [30]; $\Delta =\sqrt[3]{{\Delta }_{\mathrm{x}}\mathsf{\cdot }{\Delta }_{\mathrm{y}}\mathsf{\cdot }{\Delta }_{\mathrm{z}}}$ represents the filter scale.

The compressible FMDF transport equation is given as [28]:

$$\begin{array}{ll}\frac{\mathsf{\partial }}{\mathsf{\partial }\mathrm{t}}{\mathrm{F}}_{\mathrm{L}\mathsf{\varphi }}+\frac{\mathsf{\partial }}{\mathsf{\partial }{\mathrm{x}}_{\mathrm{i}}}\left[{\langle {\mathrm{u}}_{\mathrm{i}}\rangle }_{\mathrm{L}}{\mathrm{F}}_{\mathrm{L}\mathsf{\varphi }}\right]& =\frac{\mathsf{\partial }}{\mathsf{\partial }{\mathsf{\psi }}_{\mathsf{\alpha }}}\left[-{\langle \frac{1}{\mathsf{\rho }}\frac{\mathsf{\partial }}{\mathsf{\partial }{\mathrm{x}}_{\mathrm{i}}}\left(\mathsf{\rho }\mathrm{D}\frac{\mathsf{\partial }{\mathsf{\varphi }}_{\mathsf{\alpha }}}{\mathsf{\partial }{\mathrm{x}}_{\mathrm{i}}}\right)|\mathsf{\psi }\rangle }_{\mathrm{l}}{\mathrm{F}}_{\mathrm{L}\mathsf{\varphi }}\right]\\ & +\frac{\mathsf{\partial }}{\mathsf{\partial }{\mathsf{\psi }}_{\mathsf{\alpha }}}\left[-{\langle {\mathrm{S}}_{\mathsf{\alpha }}|\mathsf{\psi }\rangle }_{\mathrm{l}}{\mathrm{F}}_{\mathrm{L}\mathsf{\varphi }}\right]-\frac{\mathsf{\partial }}{\mathsf{\partial }{\mathrm{x}}_{\mathrm{i}}}\left[\left({\langle {\mathrm{u}}_{\mathrm{i}}|\mathsf{\psi }\rangle }_{\mathrm{l}}-{\langle {\mathrm{u}}_{\mathrm{i}}\rangle }_{\mathrm{L}}\right){\mathrm{F}}_{\mathrm{L}\mathsf{\varphi }}\right]\end{array}$$

(5)

$${\mathsf{\varphi }}_{\mathsf{\alpha }}=\{\begin{array}{l}{\mathrm{Y}}_{\mathsf{\alpha }},\hspace{1em}\mathsf{\alpha }=1,\mathsf{\dots },{\mathrm{N}}_{\mathrm{s}}\hfill \\ {\mathrm{h}}_{\mathrm{s}},\hspace{1em}\mathsf{\alpha }={\mathrm{N}}_{\mathrm{s}}+1\hfill \end{array},$$

(6)

$${\mathrm{S}}_{\mathsf{\alpha }}=\{\begin{array}{ll}\frac{{\dot{\mathsf{\omega }}}_{\mathsf{\alpha }}}{\mathsf{\rho }},& \mathsf{\alpha }=1,\mathsf{\dots },{\mathrm{N}}_{\mathrm{s}}\\ \frac{{\dot{\mathsf{\omega }}}_{\mathrm{T}}}{\mathsf{\rho }}+\frac{1}{\mathsf{\rho }}\left[\frac{\mathsf{\partial }\mathrm{p}}{\mathsf{\partial }\mathrm{t}}+{\mathrm{u}}_{\mathrm{i}}\frac{\mathsf{\partial }\mathrm{p}}{\mathsf{\partial }{\mathrm{x}}_{\mathrm{i}}}+{\mathsf{\tau }}_{\mathrm{i}\mathrm{j}}\frac{\mathsf{\partial }{\mathrm{u}}_{\mathrm{i}}}{\mathsf{\partial }{\mathrm{x}}_{\mathrm{j}}}\right],& \mathsf{\alpha }={\mathrm{N}}_{\mathrm{s}}+1\end{array},$$

(7)

where, ${\mathrm{Y}}_{\mathsf{\alpha }}$ is the mass fraction, ${\mathrm{h}}_{\mathrm{s}}$ is the sensible enthalpy, ${\dot{\mathsf{\omega }}}_{\mathsf{\alpha }}$ is the mass rate of the production of species α.

Next, a gradient-diffusion hypothesis is used to model the SGS velocity fluctuation conditioned on the scalar term:

$$\left({\langle {\mathrm{u}}_{\mathrm{i}}|\mathsf{\psi }\rangle }_{\mathrm{l}}-{\langle {\mathrm{u}}_{\mathrm{i}}\rangle }_{\mathrm{L}}\right){\mathrm{F}}_{\mathrm{L}\mathsf{\varphi }}=-{\langle \mathsf{\rho }\rangle }_{\mathrm{l}}{\mathrm{D}}_{\mathrm{t}}\frac{\mathsf{\partial }\left({\mathrm{F}}_{\mathrm{L}\mathsf{\varphi }}/{\langle \mathsf{\rho }\rangle }_{\mathrm{l}}\right)}{\mathsf{\partial }{\mathrm{x}}_{\mathrm{i}}}.$$

(8)

This work attempts to close the conditional micromixing term on the RHS of (5) using different mixing models; for instance, an IEM model result is given as follows:

$$\begin{array}{ll}\frac{\mathsf{\partial }}{\mathsf{\partial }{\mathsf{\psi }}_{\mathsf{\alpha }}}\left[-{\langle \frac{1}{\mathsf{\rho }}\frac{\mathsf{\partial }}{\mathsf{\partial }{\mathrm{x}}_{\mathrm{i}}}\left(\mathsf{\rho }\mathrm{D}\frac{\mathsf{\partial }{\mathsf{\varphi }}_{\mathsf{\alpha }}}{\mathsf{\partial }{\mathrm{x}}_{\mathrm{i}}}\right)|\mathsf{\psi }\rangle }_{\mathrm{l}}{\mathrm{F}}_{\mathrm{L}\mathsf{\varphi }}\right]& =\frac{\mathsf{\partial }}{\mathsf{\partial }{\mathrm{x}}_{\mathrm{i}}}\left({\langle \mathsf{\rho }\rangle }_{\mathrm{l}}{\langle \mathrm{D}\rangle }_{\mathrm{L}}\frac{\mathsf{\partial }\left({\mathrm{F}}_{\mathrm{L}\mathsf{\varphi }}/{\langle \mathsf{\rho }\rangle }_{\mathrm{l}}\right)}{\mathsf{\partial }{\mathrm{x}}_{\mathrm{i}}}\right)\\ & +\frac{\mathsf{\partial }}{\mathsf{\partial }{\mathsf{\psi }}_{\mathsf{\alpha }}}\left({\mathsf{\Omega }}_{\mathrm{m}}\left({\mathsf{\psi }}_{\mathsf{\alpha }}-{\langle {\mathsf{\varphi }}_{\mathsf{\alpha }}\rangle }_{\mathrm{L}}\right){\mathrm{F}}_{\mathrm{L}\mathsf{\varphi }}\right)\end{array}$$

(9)

where ${\mathsf{\Omega }}_{\mathrm{m}}$, the SGS mixing frequency, is given by ${\mathsf{\Omega }}_{\mathrm{m}}={\mathrm{C}}_{\mathsf{\phi }}\left({\langle \mathrm{D}\rangle }_{\mathrm{L}}+{\mathrm{D}}_{\mathrm{t}}\right)/\langle \mathsf{\rho }\rangle {\Delta }^2$, where $\mathrm{D}$ is the molecular diffusivity, and ${\mathrm{D}}_{\mathrm{t}}$ is the turbulent diffusivity coefficient. Furthermore, applying all the above model results to Equation (5), the model FMDF transport equation can be obtained as follows:

$$\begin{array}{ll}\frac{\mathsf{\partial }}{\mathsf{\partial }\mathrm{t}}{\mathrm{F}}_{\mathrm{L}\mathsf{\varphi }}+\frac{\mathsf{\partial }}{\mathsf{\partial }{\mathrm{x}}_{\mathrm{i}}}\left[{\langle {\mathrm{u}}_{\mathrm{i}}\rangle }_{\mathrm{L}}{\mathrm{F}}_{\mathrm{L}\mathsf{\varphi }}\right]& =\frac{\mathsf{\partial }}{\mathsf{\partial }{\mathrm{x}}_{\mathrm{i}}}\left({\langle \mathsf{\rho }\rangle }_{\mathrm{l}}\left({\langle \mathrm{D}\rangle }_{\mathrm{L}}+{\mathrm{D}}_{\mathrm{t}}\right)\frac{\mathsf{\partial }\left({\mathrm{F}}_{\mathrm{L}\mathsf{\varphi }}/{\langle \mathsf{\rho }\rangle }_{\mathrm{l}}\right)}{\mathsf{\partial }{\mathrm{x}}_{\mathrm{i}}}\right)\\ & +\frac{\mathsf{\partial }}{\mathsf{\partial }{\mathsf{\psi }}_{\mathsf{\alpha }}}\left({\mathsf{\Omega }}_{\mathrm{m}}\left({\mathsf{\psi }}_{\mathsf{\alpha }}-{\langle {\mathsf{\varphi }}_{\mathsf{\alpha }}\rangle }_{\mathrm{L}}\right){\mathrm{F}}_{\mathrm{L}\mathsf{\varphi }}\right)\\ & -\frac{\mathsf{\partial }}{\mathsf{\partial }{\mathsf{\psi }}_{\mathsf{\alpha }}}\left[{\widehat{\mathrm{S}}}_{\mathsf{\alpha }}(\mathsf{\psi }){\mathrm{F}}_{\mathrm{L}\mathsf{\varphi }}\right]\end{array}$$

(10)

In the modeled FMDF equation (Equation (10)), the velocity and pressure fields are not known and need to be obtained by solving the governing equations using finite-difference (FD) methods. Meanwhile, Equation (10) is typically solved through a Lagrangian Monte Carlo approach. Here, the multidimensional FMDF transport equation can be converted into stochastic diffusion processes:

$$\mathrm{d}{\mathrm{x}}_{\mathrm{i}}^{*}={\left({\langle {\mathrm{u}}_{\mathrm{i}}\rangle }_{\mathrm{L}}+\frac{1}{{\langle \mathsf{\rho }\rangle }_{\mathrm{l}}}\frac{\mathsf{\partial }{\langle \mathsf{\rho }\rangle }_{\mathrm{l}}\left({\langle \mathrm{D}\rangle }_{\mathrm{L}}+{\mathrm{D}}_{\mathrm{t}}\right)}{\mathsf{\partial }{\mathrm{x}}_{\mathrm{i}}}\right)}^{*}\mathrm{d}\mathrm{t}+\sqrt{2{\left({\langle \mathrm{D}\rangle }_{\mathrm{L}}+{\mathrm{D}}_{\mathrm{t}}\right)}^{*}}\mathrm{d}{\mathrm{W}}_{\mathrm{t}},$$

(11)

$$\mathrm{d}{\mathsf{\varphi }}_{\mathsf{\alpha }}^{*}={\mathrm{M}}^{*}\mathrm{d}\mathrm{t}+{\mathrm{S}}_{\mathsf{\alpha }}\left({\mathsf{\varphi }}_{\mathsf{\alpha }}^{*}\right)\mathrm{d}\mathrm{t},\mathsf{\alpha }=1,\mathsf{\dots },{\mathrm{N}}_{\mathrm{s}}+1,$$

(12)

where ${\mathrm{x}}_{\mathrm{i}}^{*}$ and ${\mathsf{\varphi }}_{\mathsf{\alpha }}^{*}$ denote the position and scalar value of a general MC particle, respectively. It should be noted that ${\mathrm{M}}^{*}$ denotes the molecular diffusion process and is determined by mixing models.

2.2. Mixing Models

2.2.1. IEM Model

The IEM model, also known as the linear mean-square estimation (LMSE) model [31], is one of the oldest and simplest deterministic mixing models. In this model, the component change of particles is only related to the mean value without considering other particles, and the scalar value of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ particle evolves according to the following formula:

$$\frac{\mathrm{d}{\mathsf{\varphi }}^{(\mathrm{i})}}{\mathrm{d}\mathrm{t}}=-{\mathsf{\Omega }}_{\mathrm{m}}\left({\mathsf{\varphi }}^{(\mathrm{i})}-\langle \mathsf{\varphi }|{\mathrm{x}}^{(\mathrm{i})}\rangle \right),$$

(13)

where ${\mathsf{\Omega }}_{\mathrm{m}}$ denotes the scalar mixing frequency, and $\langle \mathsf{\varphi }|{\mathrm{x}}^{(\mathrm{i})}\rangle$ is the mean of ${\mathsf{\varphi }}^{(\mathrm{i})}$ conditional on x. Generally, the scalar mixing frequency is linearly related to the time scale of the turbulence τ:

$${\mathsf{\Omega }}_{\mathrm{m}}=\frac{1}{{\mathsf{\tau }}_{\mathsf{\varphi }}}=\frac{{\mathrm{C}}_{\mathsf{\varphi }}}{\mathsf{\tau }}.$$

(14)

The model parameter ${\mathrm{C}}_{\mathsf{\varphi }}$ represents the correct decay rate of variance. Its role in this mixing model is to determine the degree of mixing between each particle and the cell-averaged based on the particle’s local mixing frequency value. Owing to its simple numerical process, the IEM model is widely applied to hybrid Euler/Lagrange simulations; However, because the IEM is a deterministic model, it does not contain PDF shape features. Therefore, its initial PDF shape will always be maintained and cannot relax to a Gaussian distribution. Since the change of the composition is directly related to the mean without being affected by the other particles, this model does not satisfy the requirement of localness in the composition space. Another significant IEM limitation is that the mixing rate is the same for all the components, implying that the difference in diffusion is not considered.

2.2.2. EMST Model

Using the EMST model [15] ensures that mixing occurs only between pairs of adjacent particles in the component space. Only adjacent particles are allowed to mix in the composition space, which not only maintains the three most important characteristics of the mixing model but also overcomes the problem of insufficient local characteristics in the component space of the IEM and MC models.

In the EMST approach, it is observed that a particle ensemble containing ${\mathrm{N}}_{\Delta }$ and each particle carries ${\mathrm{N}}_{\mathrm{s}}+1$ scalars. The matrix form of the EMST model evolution can be expressed as:

$$\frac{\mathrm{d}{\mathsf{\varphi }}_{\mathsf{\alpha }}^{(\mathrm{i})}}{\mathrm{d}\mathrm{t}}=-\frac{1}{{\mathsf{\omega }}^{(\mathrm{i})}}{\mathrm{M}}_{\mathrm{i}\mathrm{j}}{\mathsf{\varphi }}_{\mathsf{\alpha }}^{(\mathrm{j})},\mathsf{\cdots }\mathsf{\cdots }\mathrm{i}=1,\mathsf{\dots },{\mathrm{N}}_{\mathrm{mix}}\left({\mathrm{N}}_{\mathrm{mix}}<{\mathrm{N}}_{\Delta }\right),$$

(15)

where the elements of the interaction matrix ${\mathrm{M}}_{\mathrm{i}\mathrm{j}}$ are given by

$${\mathrm{M}}_{\mathrm{i}\mathrm{j}}=\{\begin{array}{l}-{\mathsf{\alpha }}_{\Delta }{\displaystyle \mathsf{\sum }_{\mathrm{v}=1}^{{\mathrm{N}}_{\mathrm{m}\mathrm{i}\mathrm{x}}-1}{\mathrm{B}}_{\mathrm{v}}}\left({\mathsf{\delta }}_{\mathrm{i}{\mathrm{m}}_{\mathrm{v}}}{\mathsf{\delta }}_{\mathrm{j}{\mathrm{n}}_{\mathrm{v}}}+{\mathsf{\delta }}_{\mathrm{j}{\mathrm{m}}_{\mathrm{v}}}{\mathsf{\delta }}_{\mathrm{i}{\mathrm{n}}_{\mathrm{v}}}\right),\mathrm{i}\mathsf{\ne }\mathrm{j}\hfill \\ {\mathsf{\alpha }}_{\Delta }{\displaystyle \mathsf{\sum }_{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{m}\mathrm{i}\mathrm{x}}}{\displaystyle \mathsf{\sum }_{\mathrm{v}=1}^{{\mathrm{N}}_{\mathrm{m}\mathrm{i}\mathrm{x}}-1}{\mathrm{B}}_{\mathrm{v}}}}\left({\mathsf{\delta }}_{\mathrm{i}{\mathrm{m}}_{\mathrm{v}}}{\mathsf{\delta }}_{\mathrm{j}{\mathrm{n}}_{\mathrm{v}}}+{\mathsf{\delta }}_{\mathrm{j}{\mathrm{m}}_{\mathrm{v}}}{\mathsf{\delta }}_{\mathrm{i}{\mathrm{n}}_{\mathrm{v}}}\right),\mathrm{i}=\mathrm{j}\hfill \end{array},$$

(16)

where $\mathsf{\delta }$ is the Kronecker function, the incidence matrix ${\mathrm{B}}_{\mathrm{v}}$ represents the coefficient of the ${\mathrm{v}}^{\mathrm{th}}$ extension tree that connects ${\mathrm{m}}_{\mathrm{v}}$ and ${\mathrm{n}}_{\mathrm{v}}$ particles, and ${\mathsf{\alpha }}_{\Delta }$ denotes the model parameters controlling the scalar variance dissipation. These initial model parameters of EMST need to be reconsidered under the condition of a higher mixing model coefficient ${\mathrm{C}}_{\mathsf{\varphi }}$.

It can be observed that at any given moment, there are only some particles ${\mathrm{N}}_{\mathrm{mix}}\left({\mathrm{N}}_{\mathrm{mix}}<{\mathrm{N}}_{\Delta }\right)$ that are mixing in the whole ensemble. The number of particles in the mixing state is determined by the aging process in the EMST model, which divides the fraction of particles that participate in the mixing process by half, usually during the initialization.

3. Results and Discussion

3.1. Numerical Scheme

The temporal development mixing layer is a classic viscous flow structure [32,33,34] as shown in Figure 1. It usually develops from two opposing parallel flows. This case includes basic phenomena and interactions in supersonic combustion, such as compressible effects, combustion heat release, turbulent mixing, shock waves, and flow transition. In addition, the computational grid can simulate a high Reynolds number at a small computational cost using periodic boundary continuation. In summary, it is a type of verification example of high-speed turbulent combustion that can explain the problem effectively. To explore the influence of the compressibility effect, two mixing layers with convective Mach numbers ${\mathrm{M}}_{\mathrm{c}}$ of 0.4 and 0.8 are selected for a comparative study.

The initial velocity of the mixing layer satisfies the condition of a hyperbolic tangent. To accelerate the calculation process, it is usually necessary to superimpose a small disturbance [35] on the basic flow field to stimulate the inherent instability of the flow. The initial pressure is given in

Table 1. Initialization parameters of the reaction (i).

Parameter	Fuel	Oxidant
Pressure (Pa)	29,778.54	29,778.54
Temperature (K)	288.15	288.15
Reynolds number	2800
Damköhler number	0.5
${\mathrm{Y}}_{\mathrm{F}}$	1.0	0.0
${\mathrm{Y}}_{\mathrm{O}}$	0.0	0.1
${\mathrm{Y}}_{\mathrm{P}}$	0.0	0.0

and

Table 2. Initialization parameters of the reaction (ii).

Parameter	Fuel (Top)	Oxidant (Bottom)
Pressure (Pa)	94,232.25	94,232.25
Temperature (K)	500	1200
Density (kg/m3)	0.386	0.274
${\mathrm{Y}}_{{\mathrm{H}}_2}$	0.05	0.0
${\mathrm{Y}}_{{\mathrm{O}}_2}$	0.0	0.278
${\mathrm{Y}}_{{\mathrm{H}}_2\mathrm{O}}$	0.0	0.0
${\mathrm{Y}}_{\mathrm{H}}$	0.0	0.0
${\mathrm{Y}}_{\mathrm{O}}$	0.0	0.0
${\mathrm{Y}}_{\mathrm{OH}}$	0.0	0.0
${\mathrm{Y}}_{{\mathrm{HO}}_2}$	0.0	0.0
${\mathrm{Y}}_{{\mathrm{H}}_2{\mathrm{O}}_2}$	0.0	0.0
${\mathrm{Y}}_{{\mathrm{N}}_2}$	0.95	0.722

, and the initial density is determined by the ideal gas state equation. Half of the initial vortex thickness ${\mathsf{\delta }}_{\mathsf{\omega }0}$ is selected as the characteristic length.

$${\mathsf{\delta }}_{\mathrm{w}0}={\frac{1}{2}\frac{\Delta \mathrm{u}}{{(\mathsf{\partial }\mathrm{u}/\mathsf{\partial }\mathrm{y})}_{\max }}|}_{\mathrm{t}={\mathrm{t}}_0}.$$

(17)

Periodic boundary conditions are taken in the x-direction, and simple no-reflection boundary conditions are applied in the y-direction. To match the periodic boundary, the flow direction length of the computational domain is comprised of two disturbance waves.

$$0\mathsf{\le }\mathrm{L}\mathrm{x}\mathsf{\le }2\mathsf{\times }\left(2\mathsf{\pi }/\mathsf{\alpha }\right).$$

(18)

To reduce the influence of the boundary treatment on the internal flow field, the calculation length of the y-direction takes the same scale as the flow direction.

$$-{\mathrm{L}}_{\mathrm{x}}/2\mathsf{\le }\mathrm{L}\mathrm{y}\mathsf{\le }{\mathrm{L}}_{\mathrm{x}}/2.$$

(19)

The computational domain is shown in Figure 1. The x and y directions of the flow direction are uniformly gridded. The grid sizes of LES and LES-FMDF are $121\mathsf{\times }181$ for ${\mathrm{M}}_{\mathrm{c}}=0.4$ and $201\mathsf{\times }361$ for ${\mathrm{M}}_{\mathrm{c}}=0.8$. To obtain accurate comparison results, the DNS grid size is set as $601\mathsf{\times }601$ (${\mathrm{L}}_{\mathrm{x}}={\mathrm{L}}_{\mathrm{y}}=2.8\mathsf{\times }{10}^{-2}$ m) for ${\mathrm{M}}_{\mathrm{c}}=0.4$ and $1201\mathsf{\times }1201$ (${\mathrm{L}}_{\mathrm{x}}={\mathrm{L}}_{\mathrm{y}}=5.6\mathsf{\times }{10}^{-2}$ m) for ${\mathrm{M}}_{\mathrm{c}}=0.8$.

In addition, two chemical reactions are considered in this case:

(i)

A simplicity, single-step chemical reaction without heat release;

(ii)

A hydrogen and oxygen reaction with 9 species and 21 elementary reaction steps, i.e., the O’ Conaire-2004 [36] reaction mechanism.

The former one (i) ignores the reaction heat release. Therefore, the effect of combustion on turbulence is not considered, which is convenient for analyzing and studying the mixing characteristics under high-speed conditions. In particular, the single-step irreversible reaction of fuel $\mathrm{F}$ and oxidant $\mathrm{O}$ directly to the product $\mathrm{P}$ is given as follows:

$$\mathrm{F}+\mathrm{r}\mathrm{O}\stackrel{{\mathrm{k}}_{\mathrm{f}}}{\mathsf{\to }}\left(1+\mathrm{r}\right)\mathrm{P}.$$

(20)

It is assumed that the activation energy of this reaction is zero, the stoichiometric coefficient $\mathrm{r}$ is 1, and the reaction rate is described by the law of mass action

$$\dot{\mathsf{\omega }}={\mathrm{k}}_{\mathrm{f}}{\mathrm{Y}}_{\mathrm{F}}{\mathrm{Y}}_{\mathrm{O}},$$

(21)

where ${\mathrm{k}}_{\mathrm{f}}=\mathrm{D}\mathrm{a}\mathsf{\times }{\mathrm{U}}_0/\mathrm{L}$ denotes the constant of the reaction rate, and the rate of reaction (i) is adjusted by the Damköhler number Da. The specific calculation settings are listed in Table 1.

In fact, the actual compressible reaction mixing layer is often accompanied by a large amount of energy release and density changes. Therefore, a true hydrogen/air reaction (ii) is essential for the further study of strong TCI between different mixing models. The detailed calculation parameters and initialization settings are listed in Table 2.

It should be noted that the case of the temporally developing mixing layer in this study adopts two convective Mach numbers, and two kinds of reaction estimated two mixing models: IEM and EMST. Additionally, under each focus, many kinds of mixing model coefficients and a series of instantaneous and statistical results are involved. The computational overhead caused by the 3D grid includes 7200 (base CPU hours) $\mathsf{\times }$ 2 (Mach numbers) $\mathsf{\times }$ 2 (mixing models) $\mathsf{\times }$ 2 (reactions) $\mathsf{\times }$ 6 (mixing model coefficient) and a total of 345,600 CPU hours, which is unacceptable. Therefore, a two-dimensional grid is adopted here, focusing on the development of the reaction mixing layer under two-dimensional perturbation. A value of 0.8 is considered as an optional maximum initial Mach number in this condition. At higher initial convective Mach numbers, the mixing layer development will be dominated by three-dimensional features [21]. Meanwhile, although the initial convective Mach number of the mixing layer is 0.8, local supersonic regions will gradually form after a period of development.

In the present work, a set of compressible LES-FMDF code [28,29] is used. A finite-difference scheme is used in LES code. Here, a WENO-CU6 scheme is adopted for the inviscid fluxes and the second-order centered scheme for viscous fluxes. A third-order Runge–Kutta (RK3) scheme is employed for time discretization. The simple reaction case simulations in this work use CFL = 0.5. The hydrogen/air reaction case simulations use CFL = 0.2. The same discretization schemes are used in DNS and all DNS simulation use CFL = 0.2. Moreover, in MC procedure, particle statistical information is obtained at the center of each finite-difference grid by ensemble averaging. The transfer of information from the LES to the MC particles is accomplished via a fourth-order Lagrangian interpolation.

3.2. Verification of LES-FMDF Methodology

To test the accuracy of the compressible LES-FMDF method used in this study, first, a comparison of the conventional LES, LES-FMDF, and direct numerical simulation (DNS) are introduced here. LES uses the finite rate reaction model to directly handle chemical reactions whereas the FMDF approach handles the SGS molecular mixing using the IEM model. Figure 2 shows the instantaneous results of the product mass fraction ${\mathrm{Y}}_{\mathrm{P}}$ calculated using the three methods in the case where ${\mathrm{M}}_{\mathrm{c}}=0.8$. The FMDF method has significant advantages over the mixing problem. The particles not only provide finer transport but also handle the molecular diffusion in the SGS scale, which is ignored in the LES method. Thus, the scalar interface calculated by FMDF has similar clarity as DNS. In the LES, the mixing occurs by default at a large scale, and the product mass fraction is seriously overestimated, as shown in Figure 2c.

The problem of LES overestimation of the mixture is more evident in the mean of the total product mass fraction as shown in Figure 3 and Figure 4. Before ${\mathrm{t}}^{*}=20$, when the vortex has not been rolled up, the results of LES shows a strong mixing occurring near the mixing layer. In contrast, the changes in the total product mass fraction predicted by the FMDF are in good agreement with those predicted by DNS over a long period.

Overall, the LES-FMDF method is proven to simulate scalar fields better than LES in compressible flows.

3.3. Comparison of Different Mixing Models in the Simplicity Reaction (i)

To study the effects of compressibility, we simulate two temporally developing mixing layers with two convective Mach numbers ${\mathrm{M}}_{\mathrm{c}}$ of 0.4 and 0.8. The former is a common subsonic flow problem with a slight compressibility, whereas the latter is a transonic flow with a strong compressibility. With the development of the flow, the latter case will produce local shock waves because of severe compression and expansion. To ensure better comparability between the two flows, the moment before the occurrence of the first vortex pair is selected for display. Figure 5 shows the instantaneous fields of the generated products ${\mathrm{Y}}_{\mathrm{P}}$ at a dimensionless time of 114 under the case of ${\mathrm{M}}_{\mathrm{c}}=0.8$ and 53 in the case of ${\mathrm{M}}_{\mathrm{c}}=0.4$.

It is worth noting that the cases of ${\mathrm{M}}_{\mathrm{c}}=0.8$ and ${\mathrm{M}}_{\mathrm{c}}=0.4$ have also evolved similar scalar field structures. Therefore, under the premise of no special explanation, the results at dimensionless moments ${\mathrm{t}}^{*}=114$ for ${\mathrm{M}}_{\mathrm{c}}=0.8$ and ${\mathrm{t}}^{*}=53$ for ${\mathrm{M}}_{\mathrm{c}}=0.4$ are adopted in the following paper. Figure 5 also shows a large difference as the compressibility effect increases. Not only does the gradient of parameters such as velocity, pressure, and scalar in the flow field become steeper, but there are local shock waves near the vortex in the case of ${\mathrm{M}}_{\mathrm{c}}=0.8$.

In Figure 6, the mean of the total product mass fraction ${\mathrm{P}}_{\mathrm{t}}$ is plotted versus the dimensionless time for different values of the convective Mach number. It is clearly shown that the increased compressibility infers reduced mixing [37]. To further evaluate the capability of different mixing models, the “Reynolds averaged” statistical moment of product mass fractions is calculated by the IEM and EMST methods. Hence, the “Reynolds averaged” [38,39] values ${\overline{\mathrm{Y}}}_{\mathrm{p}}$ are based on the x-direction for averaging. The resolved partial and subgrid partial Reynolds averaged values of the second moment are defined as $\overline{\mathrm{R}\left(\mathrm{a},\mathrm{b}\right)}$ and $\overline{\mathsf{\tau }\left(\mathrm{a},\mathrm{b}\right)}$, respectively.

$$\mathrm{R}(\mathrm{a},\mathrm{b})=(\langle \mathrm{a}\rangle -{\overline{\langle \mathrm{a}\rangle }}_{\mathrm{x}})(\langle \mathrm{b}\rangle -{\overline{\langle \mathrm{b}\rangle }}_{\mathrm{x}}),$$

(22)

$$\mathsf{\tau }(\mathrm{a},\mathrm{b})=\langle \mathrm{a},\mathrm{b}\rangle -\langle \mathrm{a}\rangle \langle \mathrm{b}\rangle .$$

(23)

Filtered DNS data are averaged over the LES grid, and the result is called filtered DNS (FDNS). Here, several simulations with a larger number of values of the mixing model coefficient ${\mathrm{C}}_{\mathsf{\varphi }}$ are performed. Both the IEM and EMST models in Figure 7 use the optimal mixing model coefficient ${\mathrm{C}}_{\mathsf{\varphi }}^{\mathrm{opt}}$ under the current calculation results. As shown in Figure 7, there is no significant difference between the IEM and EMST models with ${\mathrm{C}}_{\mathsf{\varphi }}^{\mathrm{opt}}$. In the present simulation, with no reactive exothermic coupling, the EMST model exhibits no obvious advantage.

The Reynolds-averaged profiles for the SGS second moments of the filtered fuel mass fraction ${\overline{\mathsf{\tau }}}_{\mathrm{L}}({\mathrm{Y}}_{\mathrm{F}},{\mathrm{Y}}_{\mathrm{F}})=\overline{{\langle {\mathrm{Y}}_{\mathrm{F}}{\mathrm{Y}}_{\mathrm{F}}\rangle }_{\mathrm{L}}}-\overline{{\langle {\mathrm{Y}}_{\mathrm{F}}\rangle }_{\mathrm{L}}}\overline{{\langle {\mathrm{Y}}_{\mathrm{F}}\rangle }_{\mathrm{L}}}$ are shown in Figure 8, which reflects the degree of unmixing in the sub-grid scale. These results also indicate that there is no significant difference between the IEM and EMST models under appropriate ${\mathrm{C}}_{\mathsf{\varphi }}$. Both the IEM and EMST models can always find a “better result” by adjusting the ${\mathrm{C}}_{\mathsf{\varphi }}$ in-mixing problem. In the calculations shown in Figure 8, ${\mathrm{C}}_{\mathsf{\varphi }}^{\mathrm{opt}}$ is most likely overestimating the molecular mixing (SGS scalar pulsation is lower than the results of the FDNS) to fit the mean of the product mass fraction ${\overline{\mathrm{Y}}}_{\mathrm{P}}$ of FDNS.

However, it is important to note that the optimal ${\mathrm{C}}_{\mathsf{\varphi }}$ required by EMST is higher than that required by the IEM model. This may be because the EMST model only works on some particles in which the age property is positive, as opposed to the IEM model, which mixes all particles in each cell. To be precise, the role of the mixing model coefficient in the EMST is not only a constant to represent the decay of variances, but also determines the rate at which particles entered and left the mixing state. Under the condition that the initial control parameters of the EMST model are not changed, it can be found in this study that the scalar variance decay rate of the EMST is weaker than that of IEM under the same mixing model coefficient, and the EMST often needs greater compensation here.

In the simplified reaction (i) considered, the mixture fraction is denoted as ${\mathrm{Z}}_{\mathrm{mix}}$:

$${\mathrm{Z}}_{\mathrm{m}\mathrm{i}\mathrm{x}}=\frac{{\mathrm{Y}}_{\mathrm{F}}-{\mathrm{Y}}_{\mathrm{O}}+1}{2},$$

(24)

where ${\mathrm{Y}}_{\mathrm{F}}$ and ${\mathrm{Y}}_{\mathrm{O}}$ are the mass fractions of the fuel and oxidant, respectively. Figure 9 shows the probability distribution function (PDF) of the mixture fraction obtained from the IEM and EMST models within the core mixing zone. Here, the core zone is set to $\mathrm{y}/{\mathsf{\delta }}_{\mathsf{\omega }0}\mathsf{\in }\left[-2.5,2.5\right]$ when the convective Mach number is 0.4, and $\mathrm{y}/{\mathsf{\delta }}_{\mathsf{\omega }0}\mathsf{\in }\left[-5,5\right]$ when the convective Mach number is 0.8. As shown in Figure 9, with appropriate ${\mathrm{C}}_{\mathsf{\varphi }}$, even the PDF distributions are nearly identical for the IEM and EMST models.

To further study the sensitivity to ${\mathrm{C}}_{\mathsf{\varphi }}$, $\mathsf{\partial }\overline{{\mathrm{Y}}_{\mathrm{P}}}/\mathsf{\partial }{\mathrm{C}}_{\mathsf{\varphi }}{|}_{{\mathrm{C}}_{\mathsf{\varphi }}^{\mathrm{opt}}}$ is considered, which is defined as the derivative of the ${\mathrm{C}}_{\mathsf{\varphi }}^{\mathrm{opt}}$ Reynolds averaged profiles of the product mass fraction. In Figure 10, it can be seen that the EMST model is less dependent on ${\mathrm{C}}_{\mathsf{\varphi }}$ because the $\mathsf{\partial }\overline{{\mathrm{Y}}_{\mathrm{P}}}/\mathsf{\partial }{\mathrm{C}}_{\mathsf{\varphi }}{|}_{{\mathrm{C}}_{\mathsf{\varphi }}^{\mathrm{opt}}}$ of the EMST model has a low value. Meanwhile, it can be seen that this sensitivity is local, and there are three obvious peaks of IEM at ${\mathrm{M}}_{\mathrm{c}}=0.4$. On the contrary, the EMST does not change significantly in the mixing range.

3.4. Comparison of Different Mixing Models in the Hydrogen/Air Reaction (ii)

For a real turbulent combustion process, the strong coupling relationship between molecular mixing, turbulent mixing, and chemical reactions cannot be ignored. The reaction considered in this section is a hydrogen/air reaction (ii) with nine components, and the detailed initialization settings are listed in Table 2 of Section 3.1. Figure 11 shows the instantaneous fields of the OH mass fraction ${\mathrm{Y}}_{\mathrm{OH}}$ at different convective Mach numbers, which are good indicators of the intensity of the reaction. As shown in Figure 11, the scalar field at ${\mathrm{M}}_{\mathrm{c}}=0.8$ shows more turbulent structures than that at ${\mathrm{M}}_{\mathrm{c}}=0.4$. Meanwhile, at higher convective Mach numbers, the reaction zone is more discrete and smaller.

To understand the structure of the flame of the compressible mixing layer, the flame index [40] is used here as follows:

$${\mathrm{F}}_{\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{x}}=\frac{\mathsf{\nabla }{\mathrm{Y}}_{{\mathrm{H}}_2}}{\left|\mathsf{\nabla }{\mathrm{Y}}_{{\mathrm{H}}_2}\right|}\mathsf{\cdot }\frac{\mathsf{\nabla }{\mathrm{Y}}_{{\mathrm{O}}_2}}{\left|\mathsf{\nabla }{\mathrm{Y}}_{{\mathrm{O}}_2}\right|}.$$

(25)

Figure 12 shows the flame index ${\mathrm{F}}_{\mathrm{index}}$ identification DNS results at different convection Mach numbers. It is also observed that the high-speed mixing layer presents an obvious partially premixed combustion mode. With the increase of the convective Mach number, the spatial and temporal scale of turbulence is further reduced, and the premixed region (${\mathrm{F}}_{\mathrm{index}}>0$), which is distributed in continuous layers at ${\mathrm{M}}_{\mathrm{c}}=0.4$, is cut apart by turbulence. When the convective Mach number is 0.8, only a small amount of premixed region exists in the center of the vortex core.

To illustrate the capability of different mixing models on the interaction of turbulence chemistry, Figure 13 shows a series of Reynolds average statistical moments of OH obtained from the IEM and EMST mixing models (with ${\mathrm{C}}_{\mathsf{\varphi }}$ = 30 in ${\mathrm{M}}_{\mathrm{c}}=0.4$; with ${\mathrm{C}}_{\mathsf{\varphi }}$ = 150 in ${\mathrm{M}}_{\mathrm{c}}=0.4$). The selection of the mixing model coefficient is shown in Appendix A. Although the results of the two mixing models deviate from those of FDNS, it is clearly shown that with the use of the detailed 9 species and 21 elementary reaction steps of the hydrogen/air reaction mechanism, the EMST model is superior to the IEM model in each case. This advantage is reflected in both the solvable scale and the SGS. On the contrary, the IEM model is found to always underestimate the reaction compared with the EMST model. At a high speed, the interface of the diffusion flame gradually tends to wrinkle to a greater degree, and the EMST model can highlight its advantages before the flame mode remains a flamelet. Because it is necessary to consider the localness of mixing on both sides of oxidizer and fuel. If the localness is ignored, the effect of a strong TCI magnifies the bias of a local weak heat release.

Furthermore, the scatter plots of the OH components in the mixture fraction space are plotted in Figure 14. The mixing fraction is often defined in terms of the mass fraction of the hydrogen atoms ${\mathrm{Z}}_{\mathrm{H}}$ in a reaction.

$${\mathrm{Z}}_{\mathrm{H}}={\mathrm{m}}_{{\mathrm{H}}_2}{\mathrm{Y}}_{{\mathrm{H}}_2}+{\mathrm{m}}_{{\mathrm{H}}_2\mathrm{O}}{\mathrm{Y}}_{{\mathrm{H}}_2\mathrm{O}}+{\mathrm{m}}_{\mathrm{H}}{\mathrm{Y}}_{\mathrm{H}}+{\mathrm{m}}_{\mathrm{O}\mathrm{H}}{\mathrm{Y}}_{\mathrm{O}\mathrm{H}}+{\mathrm{m}}_{\mathrm{H}{\mathrm{O}}_2}{\mathrm{Y}}_{\mathrm{H}{\mathrm{O}}_2}+{\mathrm{m}}_{{\mathrm{H}}_2{\mathrm{O}}_2}{\mathrm{Y}}_{{\mathrm{H}}_2{\mathrm{O}}_2},$$

(26)

$${\mathrm{Z}}_{\mathrm{m}\mathrm{i}\mathrm{x}}=\frac{{\mathrm{Z}}_{\mathrm{H}}-{\mathrm{Z}}_{\mathrm{H}}{}^{\mathrm{o}\mathrm{x}\mathrm{i}\mathrm{d}\mathrm{a}\mathrm{n}\mathrm{t}}}{{\mathrm{Z}}_{\mathrm{H}}{}^{\mathrm{f}\mathrm{u}\mathrm{e}\mathrm{l}}-{\mathrm{Z}}_{\mathrm{H}}{}^{\mathrm{o}\mathrm{x}\mathrm{i}\mathrm{d}\mathrm{a}\mathrm{n}\mathrm{t}}},$$

(27)

where ${\mathrm{m}}_{\mathrm{XX}}$ denotes the hydrogen molecular weight fraction of species ${\mathrm{Z}}_{\mathrm{H}}{}^{\mathrm{oxidant}}$, and ${\mathrm{Z}}_{\mathrm{H}}{}^{\mathrm{fuel}}$ represent the mass fraction of hydrogen atoms in the lower layer oxidant and the upper layer fuel in the initial flow field (Table 2), and ${\mathrm{Z}}_{\mathrm{st}}$ represents the stoichiometric mixture fraction (${\mathrm{Z}}_{\mathrm{st}}$ = 0.41).

As shown in Figure 14, the conditional mean of ${\mathrm{Y}}_{\mathrm{OH}}$ using the IEM model is always lower than that of the FDNS near the stoichiometric mixture fraction. Although the EMST model also underestimates the intensity of combustion near ${\mathrm{Z}}_{\mathrm{st}}$, it is in better agreement with the overall FDNS results. This also suggests that the performance of the IEM model was not as good as that of EMST in this case.

As a supplement, the information of the reaction progress variables, C, at the particle level, is given to further investigate the differences between the mixing models. The progress variable C is defined as the mass of the hydrogen element in the species of ${\mathrm{H}}_2\mathrm{O}$ over the total mass of the hydrogen element in the mixture:

$$\mathrm{C}=\frac{{\mathrm{m}}_{{\mathrm{H}}_2\mathrm{O}}{\mathrm{Y}}_{{\mathrm{H}}_2\mathrm{O}}}{{\mathrm{Z}}_{\mathrm{H}}},$$

(28)

where ${\mathrm{m}}_{{\mathrm{H}}_2\mathrm{O}}$ represents the molar mass fraction of H in the species ${\mathrm{H}}_2\mathrm{O}$. Figure 15 shows the PDFs of the progress variable in the core mixing zone. This means that the progress variables at four moments calculated by DNS in case of ${\mathrm{M}}_{\mathrm{c}}=0.4$ are 0.51 for $\mathrm{t}=60$ and 0.59 for $\mathrm{t}=70$, and those in case of ${\mathrm{M}}_{\mathrm{c}}=0.8$ are 0.49 for $\mathrm{t}=120$ and 0.56 for $\mathrm{t}=130$ in

Table 3. Mean progress variables at four moments in DNS, IEM, and EMST simulations in different cases.

The Mean Progress Variables $\overline{\mathbf{C}}$
Case 1: ${\mathbf{M}}_{\mathbf{c}}=0.4$				Case 2: ${\mathbf{M}}_{\mathbf{c}}=0.8$
	DNS	IEM	EMST		DNS	IEM	EMST
t = 60	0.51	0.48	0.52	t = 120	0.49	0.45	0.51
t = 70	0.59	0.54	0.59	t = 130	0.56	0.51	0.56

. Therefore, the PDFs compared with the mean properties are very similar. More details of the comparison are provided in the following table.

For the EMST model, the probabilities are large at both edges of the progress variable (C < 0.1, C > 0.9), but low at the middle. In contrast, for the IEM model, the PDFs of the progress variable for 0.1 < C < 0.5 are slightly larger owing to the non-local mixing. It is also observed that, with an increased convective Mach number, the requirement of localness preservation becomes even more important. For low values of the progress variable (0.1–0.3), IEM predicts that the PDFs will be slightly higher than EMST at a convective Mach number of 0.4, as shown in Figure 15b. However, at a convective Mach number of 0.8, a large number of particles driven by the IEM model are located in a poor reaction environment, as shown in Figure 15d.

Figure 16 shows the location of lower values of progress variable C, i.e., 0.1 < C < 0.3 in instantaneous results of temperature, and lower C values are often associated with lower temperatures. The EMST model performs slightly better at predicting the temperature field as ${\mathrm{M}}_{\mathrm{c}}=0.4$. The distribution of the lower progress variable marked by the black circle shows few differences, as shown in Figure 16b,c. However, there are a large number of particles of poor progress variable values in the IEM model when ${\mathrm{M}}_{\mathrm{c}}=0.8$. Particularly, in the vortex core region, the IEM shows an obvious low temperature. In addition, Figure 12 shows that when the convective Mach number is 0.8, the flame presents an obvious diffusion combustion mode. It is speculated that the EMST model highlight its advantages and produce few particles of poor progress variable values in Figure 16f.

The sensitivity of the different mixing models to the ${\mathrm{C}}_{\mathsf{\varphi }}$ in the real reaction is shown in Figure 17. It should be noted that ${\overline{\mathrm{E}}}_{\mathrm{OH}}$ is adopted as a value of the average deviation between the results ${\overline{\mathrm{Y}}}_{\mathrm{OH}}$ obtained by the different mixing models and the FDNS results under different mixing model coefficients. The value is calculated as follows:

$${\overline{\mathrm{E}}}_{\mathrm{O}\mathrm{H}}=\frac{1}{{\mathrm{N}}_{\mathrm{y}}}{\displaystyle \mathsf{\sum }_{\mathrm{i}=1}^{{\mathrm{N}}_{\mathrm{y}}}\left[\left|{\overline{\mathrm{Y}}}_{\mathrm{O}\mathrm{H}}({\mathrm{y}}_{\mathrm{i}})-{\overline{\mathrm{Y}}}_{\mathrm{O}\mathrm{H}}^{\mathrm{F}\mathrm{D}\mathrm{N}\mathrm{S}}\left({\mathrm{y}}_{\mathrm{i}}\right)\right|\right]}.$$

(29)

Figure 17 clearly shows that with the use of the EMST model, the deviation is small over the entire range of ${\mathrm{C}}_{\mathsf{\varphi }}$. This again proves that the EMST model is more capable of achieving better results than the IEM model in the present simulation. Similar to the simple reaction (i), Figure 18 shows that the sensitivity to ${\mathrm{C}}_{\mathsf{\varphi }}$ in the case of ${\mathrm{M}}_{\mathrm{c}}=0.8$ is almost an order of magnitude smaller than that in ${\mathrm{M}}_{\mathrm{c}}=0.4$. Under the action of strong TCI, the poor characteristics of IEM model in high-speed reaction cannot be compensated by adjusting the mixing model coefficient, which is different from the conclusion in Section 3.3 of the decoupling simple reaction. Moreover, when the convective Mach number is 0.4, the optimal ${\mathrm{C}}_{\mathsf{\varphi }}$ value of the two mixing models is 30, whereas the optimal ${\mathrm{C}}_{\mathsf{\varphi }}$ value is as high as 150 when the convective Mach number is 0.8. With the enhancement of the compressibility effect, the appropriate mixing model coefficient has a wider selection range, which makes it difficult to determine the optimal value. Even as the convective Mach number increased, the assumption that ${\mathrm{C}}_{\mathsf{\varphi }}$ is constant may have led to a greater error. Hence, it is more necessary to establish a dynamic mixing timescale model for LES-FMDF under high-speed reaction flow.

4. Conclusions

In this study, a compressible LES-FMDF method is employed for the temporally developing reactive mixing layer. This investigation is on the effects of two different mixing models in a compressible reactive flow. The following conclusions are drawn based on the results of the compressible flow:

• Under the appropriate mixing model coefficient, there is no significant difference between the IEM model and the EMST in the simple reaction. In addition, the probability density distributions of the core mixing regions of the two models are approximately the same. However, the effect of EMST intermittency under a high mixing model coefficient should be further discussed.
• In a real hydrogen/air reaction with a significant effect of TCI, the EMST model can predict the reaction more accurately. In addition, the requirement of localness preservation tends to be more essential as the convective Mach number increases.
• For both the IEM and EMST models for the hydrogen/air reaction, a larger value of ${\mathrm{C}}_{\mathsf{\varphi }}$ (60–180) is required to obtain the appropriate results when the convective Mach number increases. With the enhancement of compressibility, the sensitivity of the mixing model to the mixing model coefficient is greatly reduced, and the appropriate mixing model coefficient has a wider selection range, which makes it difficult to determine the optimal value.
• In the study of the simple reaction, the EMST model has a lower sensitivity than the IEM model, regardless of whether ${\mathrm{M}}_{\mathrm{c}}=0.4$ or ${\mathrm{M}}_{\mathrm{c}}=0.8$. The same conclusion is reached for the hydrogen/air reaction.

Overall, the results show that for the high-speed flows considered here, with an increase in the convective Mach number, the mixing process is gradually dominated by large-scale turbulent structures. In the simple reaction without heat release, the high-speed mixing layer is still the typical diffusion flame mode. However, after considering the real reaction, the combustion mode tends to the partially premixed diffusion flame, and the degree of non-premixed increases with the increase of the convective Mach number. In the case of strong turbulence and diffusion flame state, it is necessary to maintain the localness of the component space during the mixing process. In addition, for the mixing model coefficient, it can be seen from the results that higher convective Mach numbers require higher ${\mathrm{C}}_{\mathsf{\varphi }}$ values. This also indicates that it may not be correct for the LES-FMDF to continue using the sub-grid turbulence frequency and fixed mixing model coefficient to model the mixing frequency for a compressible high-speed flow. It is more essential to further establish a dynamic mixing timescale model to adapt to a high-speed flow.