Large Eddy Simulation of Supersonic Combustion Using the Eulerian Stochastic Fields Method

The development of supersonic combustion engines (scramjet-type) presents several challenges. Numerical simulations of scramjet engine combustion have become an attractive alternative to experimental investigations. However, supersonic combustion simulations still have several challenges: such as adequate modelling of shock/boundary layer and turbulence/chemistry interactions. The present work presents two large-eddy simulation probability density function (LES-PDF) novel formulations developed for high-speed applications. One is a conservative joint-scalar approach, where the joint-probability for reactive scalars and energy is solved, while the other is a joint velocity-scalar. The LES-PDF transport equations are solved using the Eulerian stochastic fields technique implemented in a density-based compressible solver. The performance of the models is verified through the simulations of two and three-dimensional supersonic reacting mixing layers and compared against DNS data from the literature. The results show that the joint-scalar formulation is accurate and robust, while the joint velocity-scalar closures require further development.


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 one-time one-point fine-grained Eulerian probability density function (PDF) is therefore defined as: 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 and applying a spatial filter to it [8]: where ( · ) is the filtering operator and 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 S α ≈ S α (p, ), 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: where Γ = μ/σ + μ sgs /σ sgs is the total diffusion coefficient. A sub-grid mixing timescale is defined as in [19]: where 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 α and C H 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 subgrid viscosity, μ sgs . The equations for the n th -set of Eulerian stochastic fields for the mass fractions and enthalpy can be obtained following [10]: where Y n α and H n are the stochastic variables for mass fractions of the chemical specie α and enthalpy, respectively, of the stochastic field n. The Wiener process dW n i is approximated by dt 1/2 γ , where γ = {−1, 1} is a dichotomic vector, ensuring dW i = 0. Filtered quantities are calculated from the average of the stochastic fields as Q ≈ Q .
The filtered pressure field is therefore calculated, when applying the ideal gas law, without neglecting its sub-grid fluctuations terms: where R u is the gas universal constant, MW α is the molecular weight of the chemical specie α and T n is the temperature of the n th 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].

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: where d, v i , ζ and Z α are the sample density, velocity, total energy and mass fraction, while ρ(x, t), u i (x, t), e t (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: The simplified Langevin model [2] is used to close the unknown terms for the velocity and density part. The tensor G ij is then defined as: where k sgs 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, P(x, t), instead of the filtered pressure, 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 P ≡ P( , U i , E t , Y α ) 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: The employed closure relation for the dissipation of the sub-grid kinetic energy, sgs , is the following: where the constant C is equal to 1.05. The micro-mixing constants C Y α 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, k sgs , can be directly obtained from the stochastic fields information: 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: 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).

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 (N 2 /O 2 /H 2 O) enters the domain at speed U 1 = 1151.6 m/s (Ma = 1.5) and temperature 1475 K (hereafter stream 1). The cold fuel stream (N 2 /H 2 ) inlet speed is U 2 = 669.1 m/s (Ma = 1.1) and temperature 545 K (stream 2). The vorticity thickness is defined as: where |·| 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 mediumweak compressible effects. To trigger the instability growth, a perturbation is introduced in the crosswise component of the velocity [16] proportional to the convective velocity In the two dimensional case, a maximum perturbation of 0.1U c is used, while in the three-dimensional case the perturbation is periodic in the span-wise direction and its maximum size is U c . For the two-dimensional case a domain length of L x = 350δ w,0 and L 2 = 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 L 3 = 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 H 2 /O 2 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  [15,19,20,27,28]. The results are compared with the two and three-dimensional DNS results of Ferrer [16].

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.
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 x 1 /δ w,0 = 250 the SPDF has a nonnegligible difference to the DNS results close to x 2 /δ 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. Figure 3 presents the normalised root mean square of the axial velocity fluctuation u 1 . 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 x 1 /δ w,0 = 250 and x 1 /δ w,0 = 300. 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.
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. 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 H 2 and O 2 profiles of the SPDF and VSPDF models show similar trends, although the hydrogen has been more consumed in the SPDF model. The H 2 O profile (indicative of reaction progress) shows that the combustion on the Smagorinsky model still occurs but remains constrained to 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  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 twodimensional, 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). 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 x 1 /δ 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 x 1 /δ 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. Figure 9 presents the averaged axial velocity at x 1 /δ 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 Fig. 9 Normalised averaged axial velocity at x 1 /δ w,0 = 300. Experimental data for non-reacting mixing layer of [29] ( ) and [30] ( ) Fig. 10 Normalised averaged velocity rms at x 1 /δ w,0 = 300. Experimental data for non-reacting mixing layer of [29] ( ) and [30] ( ) 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. Figure 10 shows the normalised averaged velocity rms at x 1 /δ 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 nonreactive 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].

Three-dimensional case
Chemical species mass fractions are shown in Figs. 11 (mean) and 12 (fluctuations). Major reactive species, H 2 , O 2 and H 2 O, are in good agreement with DNS results. The non-reactive species N 2 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 H O 2 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].
The mass fractions fluctuations profiles (Fig. 12) show similar behaviour. Predictions for H 2 , O 2 , H 2 O and N 2 are higher for all approaches than the DNS values. However, with

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 quasilaminar 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.