1 Introduction

Hypersonic aircraft technology can be achieved by the use of air-breathing scramjet engines. In these engines combustion occurs at supersonic speeds, which involves complex physical process such as shock-boundary layer interaction and turbulence-chemistry interaction, which occur under a wide range of scales. Computational Fluid Dynamics (CFD) simulations of such phenomena are still challenging and few numerical model exists [1] that can be used in high-speed and reacting flows.

Probability Density Function (PDF) methods can provide an elegant closure to the highly non-linear source term [2] occurring in such flows. PDF turbulence models do not require a very fine mesh, unlike quasi-linear approaches [3], and are therefore promising for turbulence-combustion modelling.

The first compressible PDF formulation was presented by Kollmann and co-workers [4, 5]. In their model the PDF description included variables, such as density, velocity, internal energy, mixture fraction and the relative rate of volume expansion. This representation allowed the fully closure of the reactive and convective terms. Nevertheless, probably due to its mathematical complexity, it was not been further explored.

Alternative compressible PDF models were developed by Delarue and Pope [6]. The authors presented a PDF formulation including velocity, internal energy, pressure and turbulence frequency. Nik et al. [7] expanded the model to account for chemical reactions. This model also provided fully closure for the reactive and convective terms in a Large Eddy Simulation (LES) framework. The presented PDF equation was then solved using a stochastic Lagrangian procedure, where Lagrangian (or notional) particles are tracked and their ensemble average represents the PDF. The process involved interpolation between the Eulerian mesh and the Lagrangian particles which lead to non-uniform spatial resolution of the PDF. To maintain uniform accuracy (and stability) additional approaches were introduced [8].

An alternative approach to the Lagrangian particles is the Eulerian formulation. In this method, first developed by Valiño [9] and Sabel’nikov and Soulard [10], the stochastic equations equivalent to the PDF are solved within a Eulerian framework. This system of equations can be easily implemented into Eulerian CFD solvers. The stochastic Eulerian approach has been used in compressible PDF simulations before [11], however, only mass fractions and enthalpy were considered and sub-grid pressure fluctuations were neglected. This approximation allowed the closure of the reactive term, as employed by [4,5,6]. This approach has also been commonly used by several authors [12,13,14] within the Lagrangian framework, however, its accuracy has not been quantitatively investigated yet.

This work investigates two new formulations based on a Eulerian compressible LES-PDF modelling, recently published in [15]. One formulation is a joint scalar PDF model [11], employing a filtered pressure closure as suggested by [13], which does not require sub-grid pressure fluctuations modelling. The other formulation is a joint velocity-scalar PDF model, resulting in a more complete approach as it can close both unresolved convective and source terms. Both LES-PDF approaches are solved using the Eulerian Monte-Carlo stochastic fields technique.

The performance of the models is investigated through two and three-dimensional supersonic reacting spatial mixing layers. Mixing layers show different degrees of mixing and turbulent strain, creating different combustion behaviour and therefore are prototype combustion benchmarks. The results are compared with predictions of the Smagorinsky model coupled with the quasi-laminar combustion (neglecting sub-grid fluctuations) and with literature DNS data [16].

2 The LES-PDF Modelling

2.1 Joint scalar PDF (SPDF)

The approach shown here has firstly been published in [15] and it is briefly revised. The joint scalar PDF model includes enthalpy and mass fractions variables into its sample space. The one-time one-point fine-grained Eulerian probability density function (PDF) is therefore defined as:

$$ f^{\prime}\left( \eta,\boldmath{Z};{x},t\right) = \delta \left( h\left( \boldmath{x},t\right) - \eta\right) \times \prod\limits_{\alpha=1}^{N_{s}} \delta (Y_{\alpha}\left( \boldmath{x},t\right) - Z_{\alpha}) $$
(1)

where η, Zα are the sample enthalpy and mass fractions, h(x, t) and Yα(x, t) are the real enthalpy and mass fractions, respectively. A non-closed filtered transport equation for LES-PDF can be obtained by deriving a transport equation for \(f^{\prime }\) and applying a spatial filter to it [8]:

$$ \begin{array}{@{}rcl@{}} \frac{\partial \overline{\rho}\widetilde{f}}{\partial t} + \frac{\partial \overline{\rho}\widetilde{u}_{i}\widetilde{f}}{\partial x_{i}} &=& \frac{\partial} {\partial x_{i}} \left( {\Gamma} \frac{\partial \widetilde{f}}{\partial x_{i}} \right) - \frac{\partial} {\partial Z_{\alpha}} \left( \overline{\rho} \widetilde{ S_{\alpha} \mid \boldmath{\Psi}} \widetilde{f} \right)\\ &&+ \frac{\partial} {\partial x_{i}} \left( \overline{\rho}\widetilde{u}_{i}\widetilde{f} - \overline{\rho} \widetilde{u_{i} \mid \boldmath{\Psi}} \widetilde{f} \right)- \frac{\partial} {\partial \eta}\left( \overline{\rho} {\widetilde{\frac{1}{\rho}\frac{D p}{D t}} \mid \boldmath{\Psi}} \widetilde{f} + \overline{\rho} \widetilde{ \frac{\tau_{ij}}{\rho} \frac{\partial u_{i}}{\partial x_{j}}} \mid \boldmath{\Psi}\widetilde{f} \right)\\ &&- \frac{\partial^{2}} {\partial \psi_{\alpha} \partial \psi_{\beta}} \left( \overline{\rho} \widetilde{{\Gamma} \frac{\partial \phi_{\alpha}}{\partial x_{i}} \frac{\partial \phi_{\beta}}{\partial x_{i}}} \mid \boldmath{\Psi} \widetilde{f} \right) \end{array} $$
(2)

where \(\left (\ \widetilde {\cdot } \ \right )\) is the filtering operator and \(\widetilde {f}\) is the filtered PDF. The proposed LES-PDF equation relies on the assumption that the reactive source term is a function of sample variables and the filtered pressure. In this way, the reactive source term is partially modelled as \(\widetilde {S}_{\alpha } \approx S_{\alpha }\left (\overline {p},\boldmath {\Psi }\right )\), where Ψ = [η, Zα], in a similar fashion as [11, 13], although it is still accurate closed to regions of small compressibility or low-Mach number flows. The Smagorinsky model [17] is used to close the convective terms. The several remaining unclosed terms are modelled with the IEM micro-mixing model [18]. The joint scalar LES-PDF equation obtained is:

$$ \begin{array}{@{}rcl@{}} \frac{\partial \overline{\rho}\widetilde{f}}{\partial t} + \frac{\partial \overline{\rho}\widetilde{u}_{i} \widetilde{f}} {\partial x_{i}} &= & \frac{\partial}{\partial x_{i}}\left( {{\varGamma}}^{\prime} \frac{\partial \widetilde{f}}{\partial x_{i}} \right) - \frac{\partial} {\partial Z_{\alpha}}\left( \overline{\rho} S_{\alpha}\left( \overline{p},\boldmath{\Psi}\right) \widetilde{f} - \frac{1}{2} \frac{C_{Y_{\alpha}}}{\tau_{sgs}} \overline{\rho} \left( Z_{\alpha} - \widetilde{Y}_{\alpha} \right) \widetilde{f} \right)\\ && - \frac{\partial}{\partial \eta}\left( \frac{D\overline{p}}{Dt}\widetilde{f} + \widetilde{\tau}_{ij} \frac{\partial \widetilde{u}_{i}}{\partial x_{j}}\widetilde{f} - \frac{1}{2} \frac{C_{H}}{\tau_{sgs}} \overline{\rho} \left( \eta - \widetilde{h} \right) \widetilde{f} \right) \end{array} $$
(3)

where \({{\varGamma }}^{\prime } = \mu /\sigma + \mu _{sgs}/\sigma _{sgs}\) is the total diffusion coefficient. A sub-grid mixing timescale is defined as in [19]:

$$ \tau_{sgs} = \left( \frac{\mu+\mu_{sgs}} {\overline{\rho}{\Delta}^{2}} \right){~}^{-1}\left( 1 - \exp(-\mathscr{R}^{2}) \right) $$
(4)

where \(\mathscr{R}\) is the ratio μsgs/μ (akin to a sub-grid Reynolds number). The time-scale reverts to a molecular one when the local flow is laminar. The model constants \(C_{Y_{\alpha }}\) and CH are equal to 2 following [20]. Equal diffusivity is assumed (unity Lewis number), where the Schmidt number σ and Prandtl number are equal to unity, while their equivalent sub-grid number, σsgs are equal to 0.7. The conventional Smagorinsky model is employed to obtain the sub-grid viscosity, μsgs. The equations for the nth-set of Eulerian stochastic fields for the mass fractions and enthalpy can be obtained following [10]:

$$ \begin{array}{@{}rcl@{}} \frac{ \partial \overline{\rho} \mathscr{Y}^{n}_{\alpha}}{\partial t} +\frac{\partial \overline{\rho}\widetilde{u}_{i} \mathscr{Y}^{n}_{\alpha}}{\partial x_{i}} &=& \frac{\partial} {\partial x_{i}} \left( {\varGamma}^{\prime} \frac{\partial \mathscr{Y}^{n}_{\alpha}}{\partial x_{i}} \right) + \overline{\rho} S^{n}_{\alpha}\left( \overline{p},\boldmath{\Psi}\right)- \frac{1}{2} \frac{C_{Y_{\alpha}}}{\tau_{sgs}} \overline{\rho} \left( \mathscr{Y}^{n}_{\alpha} - \widetilde{Y}_{\alpha} \right)\\ &&+ \left( 2 \overline{\rho} {\varGamma}^{\prime} \right){~}^{1/2} \frac{\partial \mathscr{Y}^{n}_{\alpha}}{\partial x_{i}}\frac{\text{d} {W^{n}_{i}}}{\text{d} t} \end{array} $$
(5)
$$ \begin{array}{@{}rcl@{}} \frac{\partial \overline{\rho}\mathscr{H}^{n}}{\partial t} +\frac{\partial \overline{\rho}\widetilde{u}_{i} \mathscr{H}^{n}}{\partial x_{i}} &=& \frac{\partial} {\partial x_{i}} \left( {\varGamma}^{\prime} \frac{\partial \mathscr{H}^{n}}{\partial x_{i}} \right) + \frac{D\overline{p}}{Dt} + \widetilde{\tau}_{ij} \frac{\partial \widetilde{u}_{i}}{\partial x_{j}}- \frac{1}{2} \frac{C_{H}}{\tau_{sgs}} \overline{\rho} \left( \mathscr{H}^{n} - \widetilde{h} \right)\\ &&+ \left( 2 \overline{\rho} {\varGamma}^{\prime} \right){~}^{1/2}\frac{\partial \mathscr{H}^{n}}{\partial x_{i}}\frac{\text{d} {W^{n}_{i}}}{\text{d} t} \end{array} $$
(6)

where \(\mathscr{Y}^{n}_{\alpha }\) and \(\mathscr{H}^{n}\) are the stochastic variables for mass fractions of the chemical specie α and enthalpy, respectively, of the stochastic field n. The Wiener process \(d{W^{n}_{i}}\) is approximated by dt1/2γ, where γ = {− 1,1} is a dichotomic vector, ensuring \(\left \langle dW_{i} \right \rangle = 0\). Filtered quantities are calculated from the average of the stochastic fields as \( \widetilde {Q} \approx \left \langle Q \right \rangle \).

The filtered pressure field is therefore calculated, when applying the ideal gas law, without neglecting its sub-grid fluctuations terms:

$$ \overline{p} = \overline{\rho} \widetilde{RT} \approx \overline{\rho} R_{u} \left[ \frac{1}{N_{F}} \sum\limits_{n=1}^{N_{F}} \left( \sum\limits_{\alpha=1}^{N_{s}} \frac{\mathscr{Y}^{n}_{\alpha}}{MW_{\alpha}} \right) T^{n} \right] $$
(7)

where Ru is the gas universal constant, MWα is the molecular weight of the chemical specie α and Tn is the temperature of the nth field. This approach has not been tried previously with the stochastic fields approach, only in Lagrangian formulations [13]. Accounting for sub-grid pressure fluctuations may improved accuracy in hydrocarbon fuels-pure oxygen systems [21].

2.2 Joint velocity-scalar PDF (VSPDF)

This new formulation, also published in [15], is a more complete approach. Here, the velocity components, along with density, mass fractions and total energy are included into the PDF sample space. This allows the exact closure of the convective and chemical source terms in the LES framework. The one-time one-point joint velocity-scalar fine-grained Eulerian PDF is defined as follows:

$$ \begin{array}{@{}rcl@{}} f^{\prime}(d,\boldmath{v},\zeta, \boldmath{Z};\boldmath{x},t)& = & \delta (\rho(\boldmath{x},t) - d) \prod\limits_{i=1}^{3} \delta (u_{i}(\boldmath{x},t) - v_{i})\\ &&\times\delta (e_{t}(\boldmath{x},t) - \zeta) \prod\limits_{\alpha=1}^{N_{s}} \delta (Y_{\alpha}(\boldmath{x},t) - Z_{\alpha}) \end{array} $$
(8)

where d, vi, ζ and Zα are the sample density, velocity, total energy and mass fraction, while ρ(x, t), ui(x, t), et(x, t) and Yα(x, t) are the respective real fields. Following the same methods for the SPDF, it is possible to derive the following closed transport equation:

$$ \begin{array}{@{}rcl@{}} \frac{\partial \overline{\rho} \widetilde{f}}{\partial t} + \frac{\partial \overline{\rho} v_{i} \widetilde{f}}{\partial x_{i}} &=& - \frac{\partial}{\partial d}\left( - \overline{\rho}^{2} \frac{\partial \widetilde{u}_{i}}{\partial x_{i}} \widetilde{f} \right)\\ &&-\frac{\partial}{\partial Z_{\alpha}} \left( \frac{\partial \widetilde{J}_{\alpha,i}} {\partial x_{i}} \widetilde{f} + \overline{\rho} S_{\alpha}\left( \boldmath{\Phi}\right) \widetilde{f} - \frac{1}{2} C_{Y_{\alpha}} \frac{\epsilon}{k} \overline{\rho} \left( Z_{\alpha} - \widetilde{Y}_{\alpha} \right) \widetilde{f} \right)\\ &&- \frac{\partial} {\partial v_{i}} \left( -\frac{\partial \overline{p}}{\partial x_{i}} \widetilde{f} + \frac{\partial \widetilde{\tau}_{ij}}{\partial x_{i}} \widetilde{f} + \overline{\rho} G_{ij} \left( v_{j} - \widetilde{u}_{j} \right) \widetilde{f} \right)\\ &&+ \frac{\partial^{2}}{\partial v_{i} \partial v_{i}} \left( \frac{1}{2} C_{0} \epsilon_{sgs} \widetilde{f} \right)- \frac{\partial}{\partial \zeta}\left( \frac{\partial \widetilde{q}_{i}}{\partial x_{i}} \widetilde{f} - \frac{\partial \overline{p} \widetilde{u}_{i}}{\partial x_{i}}\widetilde{f} + \frac{\partial \widetilde{\tau}_{ij}\widetilde{u}_{j}}{\partial x_{i}} \widetilde{f}\right.\\ &&\left.- \frac{1}{2} C_{e_{t}} \frac{\epsilon_{sgs}}{k_{sgs}} \overline{\rho} \left( \zeta - \widetilde{e_{t}} \right) \widetilde{f}\right) \end{array} $$
(9)

The simplified Langevin model [2] is used to close the unknown terms for the velocity and density part. The tensor Gij is then defined as:

$$ G_{ij} = -\frac{\epsilon_{sgs}}{k_{sgs}}\left( \frac{1}{2} + \frac{3}{4}C_{0} \right)\delta_{ij} $$
(10)

where ksgs and 𝜖sgs are the sub-grid kinetic energy and its dissipation respectively. The IEM micro-mixing model [18] is also chosen to close the remaining unknown terms on the total energy and mass fractions part.

Two further approximations are necessary to derive the stochastic field equations. The first is to substitute the pressure terms by the stochastic pressure field, \(\mathscr{P}(\boldmath {x},t)\), instead of the filtered pressure, \(\overline {p}\), neglecting sub-grid pressure fluctuations. This results in a different momentum equation from the Burgers’ equation, preventing the occurrence of numerical shocks. The stochastic pressure can be directly obtained from stochastic variables \(\mathscr{P} \equiv \mathscr{P}(\varrho ,\mathscr{U}_{i},\mathscr{E}_{t},\mathscr{Y}_{\alpha })\) and present the same properties of the real pressure field.

The other approximation used is to neglect the stochastic difference on the right hand side of the continuity equation if written on conservative form, so it can ensure mass conservation for all set of stochastic fields. This approximation does not affect the first-moments and is exact if the density is constant. Although this approximation can be avoided in a Lagrangian framework, the SPDEs are aimed to be coupled with a Eulerian solver. This numerical approximation act as an additional force to prevent the fields from diverge striongly by adding or removing stochastic density. The proposed set of Eulerian SPDEs is therefore:

$$ \frac{\partial \varrho^{n}}{\partial t} + \frac{\partial \varrho^{n} \mathscr{U}^{n}_{i}} {\partial x_{i}} = 0 $$
(11)
$$ \frac{\partial \varrho^{n} \mathscr{U}^{n}_{i}}{\partial t} + \frac{\partial \varrho^{n} \mathscr{U}^{n}_{j} \mathscr{U}^{n}_{i}}{\partial x_{j}} = -\frac{\partial \mathscr{P}^{n}}{\partial x_{i}} + \frac{\varrho^{n}}{\overline{\rho}} \frac{\partial \widetilde{\tau}_{ij}}{\partial x_{i}} + \varrho^{n} G_{ij} \left( \mathscr{U}^{n}_{j} - \widetilde{u}_{j} \right) + \varrho^{n} \sqrt{C_{0} \frac{\epsilon_{sgs}}{\overline{\rho}}} \frac{\text{d} {W^{n}_{i}}}{\text{d} t} $$
(12)
$$ \frac{\partial \varrho^{n} \mathscr{Y}^{n}_{\alpha}}{\partial t} +\frac{\partial \varrho^{n} \mathscr{U}^{n}_{i} \mathscr{Y}^{n}_{\alpha}}{\partial x_{i}} = \frac{\varrho^{n}}{\overline{\rho}}\frac{\partial \widetilde{J}_{\alpha,i}} {\partial x_{i}} + \varrho^{n} S_{\alpha}\left( \boldmath{\Psi}\right) - \frac{1}{2} C_{Y_{\alpha}} \frac{\epsilon_{sgs}}{k_{sgs}} \varrho^{n} \left( \mathscr{Y}^{n}_{\alpha} - \widetilde{Y}_{\alpha} \right) $$
(13)
$$ \frac{\partial \varrho^{n} \mathscr{E}^{n}_{t}}{\partial t} + \frac{\partial \varrho^{n} \mathscr{U}^{n}_{i} \mathscr{E}^{n}_{t}}{\partial x_{i}} = \frac{\varrho^{n}}{\overline{\rho}}\frac{\partial \widetilde{q}_{i}}{\partial x_{i}} - \frac{\varrho^{n}}{\overline{\rho}}\frac{\partial \overline{p} \widetilde{u}_{i}}{\partial x_{i}} + \frac{\varrho^{n}}{\overline{\rho}}\frac{\partial \widetilde{\tau}_{ij}\widetilde{u}_{j}}{\partial x_{i}} - \frac{1}{2} C_{e_{t}} \frac{\epsilon_{sgs}}{k_{sgs}} \varrho^{n} \left( \mathscr{E}^{n}_{t} - \widetilde{e_{t}} \right) $$
(14)

The employed closure relation for the dissipation of the sub-grid kinetic energy, 𝜖sgs, is the following:

$$ \epsilon_{sgs} = C_{\epsilon} k_{sgs}^{3/2}/{\Delta} $$
(15)

where the constant C𝜖 is equal to 1.05. The micro-mixing constants \(C_{Y_{\alpha }}\) and \(C_{e_{t}}\) are set to 2 and the Langevin constant is set to 2.1, following LES-PDF convention [8]. The sub-grid kinetic energy, ksgs, can be directly obtained from the stochastic fields information:

$$ k_{sgs}= \frac{1}{2} \left( \frac{1}{N_{f}} \sum\limits_{n=1}^{N_{f}} (\mathscr{U}^{n}_{i}-\widetilde{u}_{i})^{2} \right) $$
(16)

At last, the filtered variables can be obtained from the average of the Eulerian stochastic fields. For a variable Q(x, t), it is possible to recover the filtered and the Favre-filtered values:

$$ \overline{Q} = \frac{1}{N_{f}}\sum\limits_{n=1}^{N_{f}} Q^{n} \ \ \widetilde{Q} = \frac{{\sum}_{n=1}^{N_{f}} \varrho^{n} Q^{n}}{{\sum}_{n=1}^{N_{f}} \varrho^{n}} $$
(17)

The developed Eulerian stochastic differential equations are equivalent to Eq. 9 with mild assumptions. The continuity and momentum equations resemble those of Azarnykh et al. [22]. This set of equations also does not generate the numerical shocks predicted in [23] if a conventional discretisation scheme is employed.

The novelty of the present VSPDF formulation stems from both the selection of independent variables and its solution method, which has advantages in high speed flows. The VSPDF approach includes a set of variables similar to Kollmann [5] but without a dilatation variable and with density instead of pressure, unlike the formulation of Nouri et al. [24]. This has benefits when combining the method with a density-based solver, which are typically used in high-speed flows. The use of Eulerian stochastic fields method to solve the VSPDF has not been used before for fully compressible and reactive flows. Moreover, the addition of an stochastic pressure, unlike other VPDF methods [6, 23] , allowed sub-grid pressure fluctuations as well as increase the stability of the SPDE system (11)-(14).

3 Case Description

Two and three-dimensional supersonic reactive mixing layers are simulated in order to investigate the new models performance. The configuration is the same presented by Ferrer [16]. The hot oxidizer stream (N2/O2/H2O) enters the domain at speed U1 = 1151.6 m/s (Ma = 1.5) and temperature 1475 K (hereafter stream 1). The cold fuel stream (N2/H2) inlet speed is U2 = 669.1 m/s (Ma = 1.1) and temperature 545 K (stream 2). The vorticity thickness is defined as:

$$ \delta_{w} = \frac{U_{1}-U_{2}}{ \left\lvert \partial \widetilde{u}_{1} / \partial x_{2} \right\rvert_{max,2}} $$
(18)

where \(\left \lvert \cdot \right \rvert _{max,2}\) indicates the maximum of the function at the crosswise coordinate y. The initial vorticity thickness is δw,0 = 0.198 mm. The composition of the streams are described in Table 1. The convective Mach number for this test case is 0.35, resulting in medium-weak compressible effects. To trigger the instability growth, a perturbation is introduced in the crosswise component of the velocity [16] proportional to the convective velocity Uc = (c1U1 + c2U2)/(c1 + c2). In the two dimensional case, a maximum perturbation of 0.1Uc is used, while in the three-dimensional case the perturbation is periodic in the span-wise direction and its maximum size is Uc.

Table 1 Mass fractions at the inlet [16]
Full size table

For the two-dimensional case a domain length of Lx = 350δw,0 and L2 = 80δw,0 is used. A mesh size of 36000 cells (360 × 108) is employed, using a stretching grid that increases the mesh elements density towards the centre line. For the three-dimensional simulation, the planar domain is extended to L3 = 40δw,0. The mesh size for all simulations is of 2.5 million (360 × 108 × 64), also applying a stretching grid with more elements at the centreline. For the current LES, the filter width is approximately four times the Kolmogorov scale Δ ≈ 4ηk, where ηk is taken from the DNS [16] at 300δw,0 from the inlet. The mesh is relatively fine for a typical LES, with the filter width smaller that the Taylor micro-scale within the inertial sub-range.

Simulations are carried out using the in-house code CompReal. First order zero-gradient boundary conditions are used at the outlet, top and bottom of the domain. Periodic boundary conditions are applied in the span-wise direction for the three-dimensional case. The H2/O2 skeletal mechanism of Yetteret al. [25] is used to describe the hydrogen combustion using 9 species and 19 reactions. At the inlet the flow is supersonic and the boundary conditions are prescribed. In the original DNS [16] the chemical mechanism of O’Conaire et al. [26] was employed. This mechanism has the same number species as Yetter’s [25] but two more steps and has a wider pressure applicability (up to 87 bar). However, in the relatively low range of pressures in the present test-case, the difference between mechanisms is expected to be small and the slightly faster Yetter mechanism is therefore retained.

A hybrid DRP-HLLC numerical scheme is used, where the low-dispersive DRP scheme employs a 13 points stencil and fourth order accuracy in regions with small density/scalar gradients, and the HLLC-TVD second order is used in regions with sharp gradients. A fourth-order central differences discretization is used to obtain the remaining derivatives. All equations are integrated with an explicit third order Runge-Kutta scheme, whereas for the stochastic equations the order is reduced to strong 0.5-order and weak first-order accuracy. A time step of 5.0 × 10− 9 secs is employed to ensure that the acoustic Courant number does not increase above 0.2 for both cases. The simulations are performed using 8 stochastic fields in accord with other works in the literature [15, 19, 20, 27, 28]. The results are compared with the two and three-dimensional DNS results of Ferrer [16].

4 Results

4.1 Two-dimensional case

The two-dimensional mixing layer is not representative of real turbulence and therefore consequences for LES modelling should be handled with caution. Nevertheless, the proposed PDF models represent the unresolved contributions and therefore its relative behaviour can be compared. The instantaneous mass fraction of the radical OH is displayed in Fig. 1, where the mixing layer growth and the combustion process are observed for the different models. The Smagorinsky model predicts a slow growth of the mixing layer, probably due to the large associated turbulence damping. The SPDF model generates adequate mixing levels and increases the growth of vorticity thickness both by increasing combustion level and the pressure correction in Eq. 7. The VSPDF formulation follows similar trend. However, stochastic noise can be observed at the end of the domain. The vortices in the VSPDF model are more distorted than the SPDF and the OH net production is lower.

Fig. 1
figure 1

Snapshot of the \(\widetilde {Y}_{OH}\) radical at t = 250μs

Full size image

The averaged axial velocity at different downstream positions is shown in Fig. 2. The self-similarity of the velocity profiles is noticeable. Overall, SPDF and VSPDF results agree better with DNS data than Smagorinsky. Up to x1/δw,0 = 250 the SPDF has a non-negligible difference to the DNS results close to x2/δw = 0.5 (close to the flame front), with the SPDF close second. In the last station, there is a difference between SPDF and VSPDF, which can be attributed to the stronger combustion of the SPDF approach. The Smagorinsky model presents very small momentum changes and remains mostly undisturbed throughout the domain. The traditional Smagorinsky model represents the convective part by adding sub-grid viscosity, which is not enough to simulate the mixing unless more points are used.

Fig. 2
figure 2

Normalised time-averaged axial velocity

Full size image

Figure 3 presents the normalised root mean square of the axial velocity fluctuation \(u_{1}^{\prime \prime }\). The Smagorinsky model does not generate accurate results, with very low turbulent fluctuations. The SPDF and the VSPDF present good qualitative agreement with the DNS data. The VSPDF shows a close match with the DNS data at at x1/δw,0 = 250 and x1/δw,0 = 300. Both PDF approaches show more than four times the amount of fluctuations than the simple Smagorinsky.

Fig. 3
figure 3

Normalised time-averaged axial velocity fluctuation rms

Full size image

Similar behaviour is observed with the crosswise velocity fluctuations in Fig. 4, where PDF models agree well with the data downstream and the Smagorinsky model presents very small fluctuations.

Fig. 4
figure 4

Normalised time averaged crosswise velocity fluctuation rms

Full size image

Cross-correlations of the velocity components are shown in Fig. 5. Both models present good agreement with DNS data with the VSPDF showing slightly better results. Overall the velocity fluctuations are better represented by the VSPDF model, which provides an exact closure for the convective term. The Smagorinsky model suffers from the absence of an accurate combustion model and the required fine resolution to accurately reproduce the DNS velocity data.

Fig. 5
figure 5

Normalised time-averaged cross-correlation velocity fluctuation rms

Full size image

Figure 6 presents the chemical species time-averaged profiles for the Smagorinsky, SPDF and VSPDF models (two-dimensional DNS species data not available [16] ). Both models enhance mixing level with larger spread of species than the Smagorinsky model. The H2 and O2 profiles of the SPDF and VSPDF models show similar trends, although the hydrogen has been more consumed in the SPDF model. The H2O profile (indicative of reaction progress) shows that the combustion on the Smagorinsky model still occurs but remains constrained to the centreline. The Smagorinsky predictions of the HO2 radical (indicative of pre-ignition) shows an unphysical double peak, characteristic of a flame front.

Fig. 6
figure 6

Normalised time-averaged chemical species concentration at x1/δw,0 = 300

Full size image

The two-dimensional reactive mixing layer showcases the ability of the stochastic models to simulate supersonic flows and use a complex chemical mechanism for the hydrogen combustion. Both PDF approaches show a good agreement with the DNS, while the Smagorinsky model does not generate the correct level of fluctuations.

4.2 Three-dimensional case

Figure 7 shows a snapshot of OH mass-fraction distribution. Similar to the two-dimensional case, the comparison to the Smagorinsky model is performed as well. Unlike two-dimensional, the Smagorinsky models shows a “turbulent” structure. The SPDF approach presents a larger defined vortical growth than Smagorinsky, while the VSPDF model has the earliest transition to turbulence and the highest OH concentration (unlike the two-dimensional case).

Fig. 7
figure 7

Instantaneous \(\widetilde {Y}_{OH}\) distribution (only values above YOH,1 shown) . Spatial coordinates normalised by δw,0

Full size image

The vorticity thickness growth is presented in Fig. 8, comparing the DNS [16] and present LES approaches. The SPDF and the Smagorinsky models grow at a similar rate until approximately x1/δw,0 = 150, when combustion occurs the results diverge. The SPDF and Smagorinsky follow closely to the DNS data. The VSPDF’s growth rate is well above the DNS value from x1/δw,0 = 125 onwards. This suggests that the VSPDF sub-grid modelling introduces excessive momentum mixing, which is not damped by the Langevin or micro-mixing.

Fig. 8
figure 8

Vorticity thickness growth

Full size image

Figure 9 presents the averaged axial velocity at x1/δw,0 = 300 together with experimental data for non-reacting cases. As expected, the reacting cases reduce the axial-momentum due to transverse dilatation. The Smagorinsky model has the higher velocity at the end of the domain, mostly because of the delay in the overall combustion rate, and both the SPDF and VSPDF results follow closely the DNS data. It is worth noting that the SPDF model uses the same sub-grid closure for the convective term as quasi-laminar Smagorinsky, but captures better the momentum statistics through improvement of the filtered pressure coupling and the representation of sub-grid combustion processes.

Fig. 9
figure 9

Normalised averaged axial velocity at x1/δw,0 = 300. Experimental data for non-reacting mixing layer of [29] (\(\vartriangle \)) and [30] (\(\triangledown \))

Full size image

Figure 10 shows the normalised averaged velocity rms at x1/δw,0 = 300. The Smagorinsky model generates considerably better results compared to its two-dimensional counterpart (see Fig. 3). All models show reasonable qualitative agreement with the DNS data. The VSPDF generates higher fluctuation levels, closer to experimental data for non-reactive mixing layers, while the SPDF shows broader distribution of fluctuations, reaching beyond the vorticity thickness. DNS and model fluctuations levels are (in average) 10 % below experimental non-reacting data, which indicates that the fluctuation levels are not strongly correlated to the combustion in the test case [16].

Fig. 10
figure 10

Normalised averaged velocity rms at x1/δw,0 = 300. Experimental data for non-reacting mixing layer of [29] (\(\vartriangle \)) and [30] (\(\triangledown \))

Full size image

Chemical species mass fractions are shown in Figs. 11 (mean) and 12 (fluctuations). Major reactive species, H2, O2 and H2O, are in good agreement with DNS results. The non-reactive species N2 is a good indicator of the mixing levels, and although all models agree qualitatively with DNS, the SPDF approach diffuses more (also seen in major species profiles). Observing the OH distribution, SPDF and VSPDF models generate a thicker flame than the DNS. The HO2 radical DNS profile could not be reproduced by any model. This may be attributed to the different chemical mechanism used in the simulations and DNS [26].

Fig. 11
figure 11

Normalised chemical species mass fractions at x1/δw,0 = 300

Full size image
Fig. 12
figure 12

Averaged chemical species mass fractions fluctuations at x1/δw,0 = 300

Full size image

The mass fractions fluctuations profiles (Fig. 12) show similar behaviour. Predictions for H2, O2, H2O and N2 are higher for all approaches than the DNS values. However, with good qualitative agreement. The SPDF model shows more diffusive behaviour overall, with the highest fluctuations. The SPDF N2 fluctuation profile, presents two peaks, not seen in the DNS or in the other models. The OH and HO2 fluctuations are somewhat different from the DNS results, where the models predict lower fluctuation levels. The results, however, are of the same order of magnitude (even in the small HO2) and comparable to the DNS data.

5 Conclusions

Two and three-dimensional supersonic reactive mixing layer simulations have been performed using two new LES-PDF formulations using the Eulerian stochastic fields. The models have been compared with a conventional Smagorinsky approach with a quasi-laminar treatment of the chemistry (neglecting sub-grid scalar fluctuations). All models use standard values for the sub-grid micro-mixing and Smagorinsky constants and no attempt has been done to calibrate or adjust the results to the DNS data.

The present SPDF model provides closure to the filtered pressure in the case of reactive ideal gases ensuring a strong coupling between the thermochemical states. The model improves the mixing predictions compared with conventional Smagorinsky model at the same resolution; even if both use the same sub-grid transport model. The SPDF approach show reasonable agreement with DNS data while the Smagorinsky model greatly underestimates the velocity fluctuations. This is probably due to the coupling of the pressure and the chemical state.

The VSPDF model increases the turbulence levels and at the same time reduce the amount of numerical diffusion in the scalars and a -priori maintains the closest coupling between thermodynamic and velocity. The results in the planar mixing layer are the closest to experimental data, however, it overestimates the growth in the vorticity layer in the three-dimensional simulation, compared to the SPDF or Smagorinsky approaches. These results suggest that the simplified Langevin model may need to account for sub-grid compressibility effects (for example following the modified model of [6] ). The present work complements the experimental validation in [15], using the benchmark configuration of a supersonic lifted burner [31], and shows that PDF approaches are a promising tool for supersonic combustion.