Large Eddy Simulation of Combustion for High-Speed Airbreathing Engines

Abstract

Large Eddy Simulation (LES) has rapidly developed into a powerful computational methodology for fluid dynamic studies, between Reynolds-Averaged Navier–Stokes (RANS) and Direct Numerical Simulation (DNS) in both accuracy and cost. High-speed combustion applications, such as ramjets, scramjets, dual-mode ramjets, and rotating detonation engines, are promising propulsion systems, but also challenging to analyze and develop. In this paper, the building blocks needed to perform LES of high-speed combustion are reviewed. Modelling of the unresolved, subgrid terms in the filtered LES equations is highlighted. The main families of combustion models are presented, focusing on finite-rate chemistry models. The density-based finite volume method and the reaction mechanisms commonly employed in LES of high-speed H₂-air combustion are briefly reviewed. Three high-speed combustor applications are presented: an experiment of supersonic flame stabilization behind a bluff body, a direct connect facility experiment as a transition case from ramjet to scramjet operation mode, and the STRATOFLY MR3 Small-Scale Flight Experiment. Several combinations of turbulence and combustion models are compared. Comparisons with experiments are also provided when available. Overall, the results show good agreement with experimental data (e.g., shock train, mixing, wall heat flux, transition from ramjet to scramjet operation mode).

1. Introduction

During the past decade, Large Eddy Simulation (LES) [1,2,3], which is a compromise between the Direct Numerical Simulation (DNS) [4], and the Reynolds-Averaged Navier–Stokes (RANS) methods [5], has been extensively used to study high-speed combustion applications, e.g. [6,7,8,9,10,11,12,13,14,15,16,17]. In LES, turbulent scales larger than the grid spacing are resolved, while subgrid scales and the effect they have on the larger scales are modeled. This allows for capturing transient phenomena such as turbulence–chemistry, shock–shock, shock–boundary-layer, and shock–flame interactions. The increasing number of successful LES, with many demonstrating excellent qualitative and quantitative agreement with experimental data, makes LES a well-accepted method in both academia and industry. Due to the rapid advancement in high-performance computing, LES has become much more available, but physics modeling, reliability, and accuracy of LES remain to be firmly established for high-speed combustion applications as discussed in [18].

Regarding high-speed combustion applications, we refer to systems with flow velocities of the order of 1 km/s together with the additional complexities of the physics of turbulence, compressibility, shock- and rarefaction-waves, convective and radiative heat-transfer, and combustion reactions. With this complexity, it is not surprising that ramjet, scramjet, and dual-mode ramjet combustion and detonation combustion are less well understood than other propulsion systems such as jet-engines or solid propellant rockets. Many high-speed combustion model investigations depend on inviscid steady-state stream-tube formulations of compressible reacting flows, e.g., [19,20]. These are too simple, given the complex multi-physics nature of the engine types considered, for future air transportation systems but are likely to be used as engineering tools in the future. The complexity of high-speed combustion applications makes experiments very challenging and expensive to perform, usually providing only limited amounts of data. For the future, we envision an increasing degree of cooperation between experimental and computational work, with high-fidelity simulations and comprehensively instrumented experiments jointly providing a large amount of data. These data will then be processed simultaneously and in multiple ways to provide detailed insights in the complex turbulent flow physics and the chemistry that are fundamental to these novel engine cores.

Here, we describe and discuss the physical, chemical, and numerical aspects of LES to provide a comprehensive overview of the methodology and its advantages and disadvantages. In addition, we also provide a few examples of LES results and how they compare with experimental data to inform what can be achieved using LES in terms of accuracy, cost, and information based on selected results from different LES computations of high-speed combustion. This work ends with a brief description of the projected future of LES of high-speed combustion, as well as a summary of features of LES that need to be improved for the methodology to gain even more widespread use, in particular among industrial or engineering practitioners.

2. Large Eddy Simulation of High-Speed Airbreathing Propulsion

During the last decade, the interest in high-speed airbreathing propulsion has increased significantly, driven by its many potential scientific, economic, and strategic opportunities. Many aspects of this subject remain unsolved, and efforts from both the experimental and numerical sides have been devoted to achieving the goal of practical hypersonic flight. Here, we review a selection of relevant works up to the current state-of-the-art; for more details, the reader is referred to the review by Huang et al. [21].

Ingenito and Bruno [6] developed a finite-rate chemistry LES combustion model by incorporating a compressibility correction to the filtered chemical source terms. This model was validated against jet in cross-flow supersonic combustion data from experiments performed at NASA Dryden Research Center. Genin and Menon [7] used a combination of a dynamic subgrid flow model, a quasi-laminar combustion model, and a solver combining shock-capturing and low-dissipation schemes to study supersonic combustion behind a wedge. Koo et al. [8] combined a conventional Smagorinsky subgrid flow model with a Eulerian probability density function turbulent combustion model using a quadrature-based method. This approach was validated on the United Technologies Research Center (UTRC) scramjet model combustor. Fureby et al. [9] combined the mixed subgrid flow model with the partially stirred reactor turbulent combustion model to study supersonic combustion behind strut- and hypermixer injector configurations. In [9], an explicit low dissipation total variation diminishing numerical scheme was used.

Saghafian et al. [11] combined the dynamic Smagorinsky subgrid flow model with a Compressible Flamelet/Progress Variable (CFPV) model. To account for the strong variations in temperature and pressure, compressibility corrections were introduced in the modeled reaction rates. Comparisons with experiments for H₂-air combustion in a cross-flow configuration showed good agreement. Another notable effort is the HyShot II experiment, conducted in the High Enthalpy Shock Tunnel in Göttingen. Several research teams have successfully computed this case: Larsson et al. [12] used the Vreman subgrid flow model, a compressible CFPV combustion model, and a solution-adaptive numerical scheme involving a shock sensor and hybrid low dissipation/shock-capturing scheme. Nordin-Bates et al. [14] used a dynamic subgrid flow model and the partially stirred reactor turbulent combustion model to study the same case, showing good agreement with experimental data and providing detailed flow physics descriptions.

The HIFiRE 2 experiment in the Arc Heated Scramjet Test Facility at the NASA Langley Research Center has also been studied by several research teams: Lacaze et al. [15] used the mixed subgrid flow model, a quasi-laminar combustion model, and an all-Mach-number solver, whereas Saghafian et al. [22] used the previously developed model of Saghafian et al. [11]. Fulton et al. [13] used a hybrid RANS-LES method based on the k-ω turbulence model and the Lenormand subgrid model together with two different partially stirred reactor combustion models to examine supersonic H₂-air combustion in the University of Virginia ramp combustor. Another application is found in Hassan et al. [16] who focused on C₂H₄-air combustion in a cavity combustor tested at the US Air Force Research Laboratory. Here, a dynamic hybrid RANS-LES model and a quasi-laminar combustion model were used. Later, Ruan et al. [17] computed the same case utilizing a dynamic subgrid flow model and a quasi-laminar combustion model.

3. Large Eddy Simulation Modeling

LES is based on filtered mass, momentum, and energy equations. The low-pass filtering is used to separate the resolved scale flow (denoted by tildes and overbars, depending on whether density-weighted filtering is used) from the unresolved (or subgrid) scale flow. The filtering implies that the physics on scales larger than the filter width, ∆, can be resolved, whereas the physics on scales smaller than ∆ must be modeled. For a linear viscous mixture with Fourier heat conduction and Fickian diffusion, the filtered mass, momentum, and energy equations are

$$\{\begin{array}{l}{\partial }_t(\overline{\rho })+\nabla \cdot (\overline{\rho }\tilde{\mathit{v}})=0\hfill \\ {\partial }_t(\overline{\rho }{\tilde{Y}}_i)+\nabla \cdot (\overline{\rho }\tilde{\mathit{v}}{\tilde{Y}}_i)=\nabla \cdot ({\overline{\mathit{J}}}_i-{\mathit{b}}_i)+{\overline{\dot{w}}}_i\hfill \\ {\partial }_t(\overline{\rho }\tilde{\mathit{v}})+\nabla \cdot (\overline{\rho }\tilde{\mathit{v}}⨂\tilde{\mathit{v}})=-\nabla \overline{p}+\nabla \cdot (\overline{\mathit{S}}-\mathit{B})\hfill \\ {\partial }_t(\overline{\rho }\tilde{E})+\nabla \cdot (\overline{\rho }\tilde{\mathit{v}}\tilde{E})=\nabla \cdot (-\overline{p}\tilde{\mathit{v}}+\overline{\mathit{S}}\tilde{\mathit{v}}+\overline{\mathit{h}}+{\Sigma }_{i=1}^N(h_{i,f}^{\theta }{\overline{\mathit{J}}}_i)-{\mathit{b}}_E)+{\Sigma }_{i=1}^N({\overline{\dot{w}}}_i{h}_{i,f}^{\theta })\hfill \end{array}$$

(1)

in which $\rho$ is the density, $\mathit{v}$ the velocity, $Y_i$ the species mass -fractions, ${\mathit{j}}_i \approx D_i\nabla Y_i$the species mass flux vectors, $D_i$ the species diffusivities, ${\dot{w}}_i=P_{ij}{\dot{w}}_j$ the species formation rates, $P_{ij}$ the stoichiometric matrix, and ${\dot{w}}_j$ the reaction rates. In the momentum equation (1₃), $p=\rho RT$ is the pressure, $R$ the (composition-dependent) gas constant, $T$ the temperature, and $\mathit{S} \approx 2\mu {\mathit{D}}_D$ the viscous stress tensor with $\mu$ the viscosity and ${\mathit{D}}_D=\frac{1}{2}{(\nabla \mathit{v}+\nabla {\mathit{v}}^T)}_D$ the deviatoric part of the rate-of-strain tensor. In the energy equation (1₄), $E={\sum }_{i=1}^N(Y_i{\int }_{T_0}^T{C}_{p,i}dT)-p/\rho +v^2/2$ is the total energy, being the sum of the sensible energy and the kinetic energy, $-p\mathit{v}$ the pressure work, $\mathit{S}\mathit{v}$ the viscous work, $\mathit{h} \approx \kappa \nabla T$ the heat flux, $\kappa$ the thermal diffusivity, and $\sum _{i=1}^N({\dot{w}}_i{h}_{i,f}^{\theta })$ the chemical heat release. The influence of the unresolved (or subgrid) flow is collected in the subgrid stress tensor $\mathit{B}=\overline{\rho }(\widetilde{\mathit{v}⨂\mathit{v}}-\tilde{\mathit{v}}⨂\tilde{\mathit{v}})$ and the subgrid flux vectors ${\mathit{b}}_i=\overline{\rho }(\widetilde{\mathit{v}{Y}_i}-\tilde{\mathit{v}}{\tilde{Y}}_i)$ and ${\mathit{b}}_E=\overline{\rho }(\widetilde{\mathit{v}E}-\tilde{\mathit{v}}\tilde{E})+(\overline{p\mathit{v}}-\overline{p}\tilde{\mathit{v}})-(\overline{\mathit{S}\mathit{v}}-\overline{\mathit{S}}\tilde{\mathit{v}})$. Here, $D_i$ and $\kappa$ are obtained from the viscosity, $\mu$, using Schmidt, ${\mathrm{Sc}}_i$, and Prandtl, $\mathrm{Pr}$, numbers [23], respectively. Equation (1) forms a set of $N$ + 5 equations, where $N$ is the number of species, but due to the constraint $\sum _i{Y}_i=1$, one species (here $Y_{N2}$) can be discarded.

3.1. LES Subgrid Flow Modeling

The subgrid stress tensor and flux vectors, $\mathit{B}$, ${\mathit{b}}_E$, and ${\mathit{b}}_i$, in the filtered transport Equation (1) can be closed using a wide variety of models. Sagaut [1] provides a comprehensive review of subgrid models for non-reactive and incompressible flows. Most often these models are extended to compressible and reactive flows using dimensional arguments. This is questionable, but in the absence of experimental or direct numerical simulation data, this constitutes today’s standards. Based on literature surveys, the most common subgrid models for high-speed reactive flows are functional subgrid viscosity models [1], of the form $\mathit{B}=-2{\mu }_k{\tilde{\mathit{D}}}_D+\frac{2}{3}\overline{\rho }k\mathit{I}$, where ${\mu }_k$ is the subgrid viscosity and $k$ the subgrid kinetic energy, ${\mathit{b}}_i=\left({\mu }_k/S{c}_t\right)\nabla Y_i$, and ${\mathit{b}}_E=\left({\mu }_k/P{r}_t\right)\nabla E$, where $S{c}_t \approx 0.7$ and $P{r}_t \approx 0.8$ are the turbulent Schmidt and Prandtl numbers, respectively. These models are by construction purely dissipative, and different techniques can be used to provide the necessary closure expressions for $k$ and ${\mu }_k$. The most used model is the Smagorinsky (SMG) model [24], by which ${\mu }_k=c_D\overline{\rho }{\mathsf{\Delta }}^2\left|\tilde{\mathit{D}}\right|$ and $k=c_I{\mathsf{\Delta }}^2{\left|\tilde{\mathit{D}}\right|}^2$, in which $c_D$ and $c_I$ are model coefficients that can be estimated from the inertial subrange behavior or dynamically. A model with higher flexibility and accuracy is the Localized Dynamic k-equation Model (LDKM) [25], in which a modeled transport equation, ${\partial }_t\left(\overline{\rho }k\right)+\nabla \cdot \left(\overline{\rho }k\tilde{\mathit{v}}\right)=-\mathit{B}\cdot \tilde{\mathit{D}}+\nabla \cdot \left({\mu }_k\nabla k\right)-c_{\epsilon }\overline{\rho }{k}^{3/2}/\mathsf{\Delta }$, is solved for $k$, and the subgrid viscosity is modeled as ${\mu }_k=c_k\overline{\rho }\mathsf{\Delta }\sqrt{k}$. Here, the model coefficients, $c_k$ and $c_{\epsilon }$, are evaluated dynamically using scale similarity. Another model, designed for high-speed flows with shocks and rarefaction waves, is the Hyperviscosity (HV) model, $\mathit{B}=-2{\mu }_k{\tilde{\mathit{D}}}_D+\overline{\rho }\left(\beta \left(\nabla ·\tilde{v}\right)+\frac{2}{3}k\right)\mathit{I}$ [26]. This model aims to improve the modeling of shocks using combinations of shear and bulk viscosities so that ${\mu }_k=c_{\nu }\overline{\rho }{\mathsf{\Delta }}^4\overline{\left|({\nabla }^2\left|\tilde{\mathit{D}}\right|)\right|}$ and $\beta =c_{\beta }{\mathsf{\Delta }}^4\overline{\left|{\nabla }^2\tilde{\mathit{D}}\right|}$, in which $c_{\nu }=0.025$ and $c_{\beta }=5.0$ are model coefficients assuming that $k$ can be absorbed into the spherical part.

3.2. LES Combustion Modeling

The filtered species formation rates in the species transport equations, ${\overline{\dot{w}}}_i$, require separate treatment to describe the effect of the combustion chemistry on the flow and how the turbulence affects the combustion chemistry. The main reason for this is the inherent non-linearities in species mass-fractions and temperature in the reaction rates, the wide range of scales involved, and the large number of terms that arise due to the reaction mechanism. Several methods to model the filtered species formation rates exist, and they can generally be divided into the following.

Flamelet models assume that the region of combustion is thin compared to the dimensions of the flow and behaves like an interface between fuel and oxidizer (in non-premixed combustion) or between reactants and products (in premixed combustion) [27]. If scale separation applies, the flow and chemistry can be decoupled such that the species transport Equation (1₂) may be replaced by a small number of transport equations, usually for a mixture fraction, $z$, and a progress variable, $c$. The chemistry can then be computed and tabulated using canonical flames, resulting in a flamelet library, parameterized by $z$ and $c$, and by additional variables such as $p$ or the scalar dissipation rate, $\chi$, so that $Y_i^{lib}=Y_i^{lib}\left(z,c,\dots \right)$. The filtered mass fractions, ${\tilde{Y}}_i$, can then be computed from $\tilde{z}$, $\tilde{c}$, … and $Y_i^{lib}$, using a statistical approach by which ${\tilde{Y}}_i={\int }_0^1{\int }_0^1\mathcal{P}(z,c)Y_i^{lib}(z,c)$, where $\mathcal{P}=\mathcal{P}(z,c)$ is a Probability Density Function (PDF). In practice a presumed PDF is often used, in which Beta and Gaussian functions are used for $z$ and $c$, respectively. This requires additional models for the variances of the mixture fraction and progress variable ${z^{″}}^2$ and ${c^{″}}^2$. The source term of the progress variable, ${\overline{\dot{w}}}_c$, needs also to be modeled, for example from analytical expressions for the flame surface density, $\mathsf{\Sigma }$, wrinkling, $\Xi$, and the laminar flame speed, $s_u$ [28], or using tabulated values based on the flamelet library [29]. Flamelet models are efficient, but ${\overline{\dot{w}}}_c$ may require additional modeling to account for compressibility effects depending on the formulation, e.g., [30].

Finite-rate chemistry models are based on the species transport equation (1₂) using closure models for the filtered species formation rates, ${\overline{\dot{w}}}_i$, explicitly considering the underlying chemical reaction mechanism and associated thermal and transport properties. An assortment of LES finite-rate chemistry models of varying complexity developed based on different assumptions are available including the Thickened Flame Model (TFM) [31], the Eddy Dissipation Concept (EDC) model [32], the Fractal Model (FM) [33], the Partially Stirred Reactor (PaSR) model [34], the Filtered Density Function (FDF) model [35], the Eulerian Stochastic Fields (ESF) model [36], and the more advanced Linear Eddy Model (LEM) [37]. Based on the experiences from previous investigations, e.g., [38], we hereafter use the EDC, FM, PaSR, and ESF models as typical representative models of this class of LES combustion models.

The EDC model [32], developed by Magnussen [39], following the ideas of Chomiak [40], which were in turn based on observations by Batchelor and Townsend [41], suggests that at high Re numbers, turbulent fine structures are not uniformly distributed but concentrated in smaller localized regions. By defining the filtered mass fractions, temperature, and species formation rates as ${\tilde{Y}}_i={\gamma }^{*}{Y}_i^{*}+\left(1-{\gamma }^{*}\right)Y_i^0$, $\tilde{T}={\gamma }^{*}{T}^{*}+\left(1-{\gamma }^{*}\right)T^0$, and ${\overline{\dot{w}}}_i={\gamma }^{*}{\dot{w}}_i^{*}+\left(1-{\gamma }_i^{*}\right){\dot{w}}_i^0$, in which * and ⁰ denote values in the fine structures and surroundings, respectively, and ${\gamma }^{*}$ is the fine structure volume fraction, subgrid equations of mass and energy can be formulated such that $\overline{\rho }(Y_i^{*}-{\tilde{Y}}_i)=(1-{\gamma }^{*}){\tau }^{*}{\dot{w}}_i^{*}$ and $\overline{\rho }\sum _{i=1}^N(Y_i^{*}{h}_i^{*}-{\tilde{Y}}_i{\tilde{h}}_i)=(1-{\gamma }^{*}){\tau }^{*}\sum _{i=1}^N{h}_{i,f}^{\theta }{\dot{w}}_i^{*}$, in which ${\tau }^{*}$ is the subgrid mixing time scale. The fine structure reacting volume fraction and the subgrid time scale can be estimated from a model of the turbulent cascade process, first proposed in [39] and later modified for LES in [32], such that ${\gamma }^{*} \approx 1.02{(\nu /\mathsf{\Delta }{k}^{1/2})}^{3/4}$ and ${\tau }^{*}=1.24{(\mathsf{\Delta }\nu /k^{3/2})}^{1/2}$.

The FM model [33] is based on the same basic assumptions as the EDC model but uses a somewhat simpler formulation, ${\overline{\dot{w}}}_i \approx {\gamma }^{*}\widetilde{\dot{w_i}}$, in which $\widetilde{\dot{w_i}}$ denote the species formation rates evaluated at the local LES state, i.e., $\widetilde{\dot{w_i}}={\dot{w}}_i(\overline{\rho },{\tilde{Y}}_i,\tilde{T})$. The fine structure volume fraction is in this model estimated using fractal theory so that ${\gamma }^{*}={\gamma }_N{(\mathsf{\Delta }/{\ell }_K)}^{D_3-3}$, in which ${\gamma }_N=N_{\ell ,K}/N_T$ is the ratio of the number of Kolmogorov scales, $N_{\ell ,K}$, to the overall number of scales, $N_T$, generated locally, whereas $D_3$ is the fractal dimension of the turbulence. Given theoretical estimates for Kolmogorov turbulence for $D_3$, $N_{\ell ,K}$, and $N_T$ [33], the FM model is closed. For high Re numbers the fine structure volume fraction asymptotically approaches the limit value ${\gamma }^{*} \approx 0.314$.

The PaSR model [34] is also based on the same assumptions as the EDC and FM models but considers the filtered species formation rates to be obtained using a bi-modal PDF, distinguishing between the fine structure and surroundings, so that $\mathcal{P}={\gamma }^{*}\delta ({\psi }^{*})+(1-{\gamma }^{*})\delta ({\psi }^0)$. Here, $\psi \in \left\{\rho ,Y_i,T\right\}$ denotes the thermochemical state. The resulting rates are ${\overline{\dot{w}}}_i=\int \mathcal{P}\left(\psi \right){\dot{w}}_i\left(\psi \right)d\psi ={\gamma }^{*}{\dot{w}}_i^{*}+(1-{\gamma }^{*}){\dot{w}}_i^0$, which can be simplified as ${\overline{\dot{w}}}_i \approx {\gamma }^{*}\widetilde{\dot{w_i}}$ since ${\dot{w}}_i$ are most often negligible outside of the fine structure regions. The modeling of ${\gamma }^{*}$ is based on a geometrical analysis [32], resulting in ${\gamma }^{*} \approx {\tau }_c/\left({\tau }^{*}+{\tau }_c\right)$, where ${\tau }_c \approx {\delta }_u/s_u$ is a global representative chemical time scale, whereas ${\tau }^{*}$ is based on the fact that the fine structure area-to-volume ratio is determined by the dissipative length scale, ${\ell }_D={\left(\nu /\left(v^{\prime }/\mathsf{\Delta }\right)\right)}^{1/2}$, and that the velocity influencing this is the Kolmogorov velocity scale, $v_K$, such that ${\tau }^{*}={\ell }_D/v_K=\sqrt{{\tau }_K{\tau }_{\mathsf{\Delta }}}$. Here, ${\tau }_{\mathsf{\Delta }}=\mathsf{\Delta }/v^{\prime }$ represents turbulent dissipation and small-scale mixing and ${\tau }_K$ is the Kolmogorov time scale.

The Eulerian Stochastic Fields (ESF) model [36] is a stochastic rather than deterministic method of representing the probability density function class of combustion models. In the ESF approach, the PDF, $𝒫$, is represented by an ensemble of $N_s$ stochastic fields with each field encompassing the $N$ scalars, ${\xi }_{\alpha }^n$ for $1\le n\le N_s$ and $1\le \alpha \le N$. In the present work, the Itô formulation of the stochastic integral is adopted, and the stochastic fields evolve according to

$$\begin{gathered} \overline{\rho }d{\xi }_{\alpha }^n+\overline{\rho }\tilde{\mathit{v}}\cdot \nabla {\xi }_{\alpha }^{n}dt=\nabla \cdot \left(\sigma \nabla {\xi }_{\alpha }^n\right)dt+\overline{\rho }\sqrt{\frac{2\sigma }{\overline{\rho }}}\nabla {\xi }_{\alpha }^n\cdot d{\mathit{W}}_{\alpha }^n \\ -\frac{1}{2}\overline{\rho }\frac{C_d}{{\tau }_k}({\xi }_{\alpha }^n-{\tilde{\varphi }}_{\alpha })dt+{\dot{w}}_{\alpha }^n\left({\xi }_{\alpha }^n\right)dt, \end{gathered}$$

(2)

where $\sigma$ is the total diffusion, and $d{\mathit{W}}_{\alpha }^n$ the increments of a Wiener process, different for each field but independent of the spatial location. The stochastic fields given by (2) form an equivalent stochastic system to the PDF, smooth over the scale of the filter width. Following [36], eight stochastic fields are solved for each scalar variable, making this approach very expensive.

4. Numerical Methods for High-Speed Flows

Numerical methods for LES of high-speed combustion with shocks, contact discontinuities, and rarefaction waves are both important and challenging. Finite Volume Methods (FVM) prevail in fluid dynamics, in particular for realistic geometries. A specific requirement is the need for low-dissipation schemes to capture most of the resolved turbulence. These schemes are, however, poorly suited to capture shocks and other flow discontinuities, and to circumvent this, artificial dissipation and diffusion techniques [42] and hybrid schemes [43], which combine two schemes and blend different flux reconstruction algorithms, have been developed. Next, we briefly describe the density-based FVM employed in the example studies presented in Section 5.

The reactive LES Equations (1) are discretized using Gauss theorem in conjunction with a multi-step time-integration scheme. The discretized continuity Equation (1₁) reads

$$\frac{1}{\Delta t}\left({\rho }_P^{n+1}-{\rho }_P^n\right)+\frac{1}{d{V}_P}{{\displaystyle \sum }}_f{\left[\rho \mathit{v}\right]}_f^n·d{\mathit{A}}_f=0$$

(3)

in which the subscripts denote the cell center location P and cell faces f, the superscripts denote the time indices, $d{V}_p$ the control volume P, $d{\mathit{A}}_{\mathrm{f}}$ the area of face f with direction along the surface normal, ∆t the time-step, and the summation runs over all cell faces of P. Equation (1₃) is discretized similarly but using an operator splitting approach with an inviscid and a viscous step:

$$\begin{gathered} \frac{1}{\Delta t}\left({\left(\rho \mathit{v}\right)}_P^{*}-{\left(\rho \mathit{v}\right)}_P^n\right)+\frac{1}{d{V}_P}{{\displaystyle \sum }}_f{\left[\left(\rho \mathit{v}\otimes \mathit{v}\right)+p\mathit{I}\right]}_f^n·d{\mathit{A}}_f=0, \\ \frac{1}{\Delta t}({\left(\rho \mathit{v}\right)}_P^{n+1}-{\left(\rho \mathit{v}\right)}_P^{*})=\frac{1}{d{V}_P}{{\displaystyle \sum }}_f {[{\mu }_{eff}\nabla \mathit{v}]}_f^{*}d{\mathit{A}}_f+\frac{1}{d{V}_P}{{\displaystyle \sum }}_f {[{\mu }_{eff}(\nabla {\mathit{v}}^T-\frac{2}{3}tr\left(\mathit{D}\mathit{I}\right))]}_f^{*}d{\mathit{A}}_f, \end{gathered}$$

(4)

in which ${\mu }_{eff}$ is the sum of molecular and subgrid viscosities. The energy Equation (1₄) is discretized by a similar operator splitting approach such that

$$\frac{1}{\Delta t}\left({\left(\rho E\right)}_P^{*}-{\left(\rho E\right)}_P^n\right)+\frac{1}{d{V}_P}{{\displaystyle \sum }}_f{[\rho \mathit{v}\left(E+\frac{p}{\rho }\right)]}_f^n·d{\mathit{A}}_f=\frac{1}{d{V}_P}{{\displaystyle \sum }}_f{\left[\mu \left(\nabla \mathit{v}+\nabla {\mathit{v}}^T\right)\mathit{v}\right]}_f^n·d{\mathit{A}}_f,$$

$$\frac{1}{\Delta t}({\left(\rho e\right)}_P^{n+1}-{\left(\rho e\right)}_P^{*})=\frac{1}{d{V}_P}{{\displaystyle \sum }}_f{\left[-\frac{{\kappa }_{eff}}{c_V}\nabla e\right]}_f^{*}·d{\mathit{A}}_f,$$

(5)

in which $e=E-{\mathit{v}}^2/2$ is the internal energy and κ_eff the sum of molecular and subgrid thermal diffusivities. The species transport Equations (1₂) are discretized similarly so that

$$\frac{1}{\Delta t}({(\rho Y_i)}_P^{*}-{(\rho Y_i)}_P^n)+\frac{1}{d{V}_P}{{\displaystyle \sum }}_f{\left[\rho \mathit{v}{Y}_i\right]}_f^n·d{\mathit{A}}_f=0,$$

$$\frac{1}{\Delta t}({(\rho Y_i)}_P^{n+1}-{(\rho Y_i)}_P^{*})=\frac{1}{d{V}_P}{{\displaystyle \sum }}_f{\left[\rho D_{eff}\nabla Y_i\right]}_f^n·d{\mathit{A}}_f=\frac{1}{\Delta t}({(\rho {\widehat{Y}}_i)}_P^{n+1}-{(\rho {\widehat{Y}}_i)}_P^{*}),$$

(6)

$${\displaystyle \sum }_{i=0}^m{\left[{\alpha }_i\rho {\widehat{Y}}_i\right]}_P^{n+i}=\Delta t{\displaystyle \sum }_{i=1}^m{\left[{\beta }_i{\dot{w}}_i\right]}_P^{n+i}$$

in which D_eff is the sum of molecular and subgrid mass diffusivities. In (6₃), the combustion chemistry is treated using Strang splitting [44] and integrated using a Rosenbrock solver [45]. The convective fluxes are reconstructed using the Kurganov central scheme [46,47], whereas the diffusive fluxes are reconstructed using local linear interpolation. The inviscid equations are solved explicitly whereas the viscous equations are solved using Gauss–Seidel [48].

5. Combustion Chemistry

For high-speed combustion, hydrogen (H₂) and small hydrocarbon fuels, such as ethylene (C₂H₄), are favorable due to faster ignition, whereas larger hydrocarbons, such as kerosene, are favorable from an energy contents point of view. Most research studies on high-speed combustion are performed with H₂, but also C₂H₄ is sometimes considered. Following [49,50], the choice of chemical reaction mechanism is generally found to be important for accurate and reliable LES as confirmed in [38] for LES of high-speed combustion. Based on these observations, it is then natural to study this aspect in greater detail also for other high-speed combustion applications.

For H₂-air combustion, Figure 1 compares experimental data and mechanism predictions for the laminar flame speed, s_u, ignition time, τ_ign, and extinction strain rate, σ_ext, for the global mechanism of Marinov (M1) [51]; the skeletal mechanisms of Eklund and Stouffer (ES7) [52], Baurle and Girimaji (BG7) [53], and Davidenko et al. (D7) [54]; the detailed mechanisms of Jachimovski (J20) [55], Alekseev et al. (K30) [56], and Wang et al. (USCII) [57]; as well as the pathway-centric mechanism of Zettervall and Fureby, (Z22) [50]. For τ_ign, Figure 1a, all mechanisms examined perform similarly above 1000 K, but with deviations in the crossover regime, 850 K < T < 1000 K, and below. The reaction pathways of the global and skeletal mechanisms, all neglecting the low temperature of H₂O₂ and HO₂ chemistry, are unable to reproduce the bending of τ_ign in the crossover region. For the detailed mechanisms, this bending appears too strong, resulting in too long ignition delay times for T < 900 K. The pathway-centric Z22 mechanism is, however, able to better reproduce the τ_ign data across the whole range of temperatures. Concerning s_u, in Figure 1b, M1 overpredicts s_u whereas ES7, BG7, and D7 all underpredict s_u. J20, K30, USCII, and Z22 all show good agreement with the experimental data. For the laminar flame temperature, T_fl (not shown), all mechanisms show good agreement with each other and with the experimental data. For σ_ext, Figure 1c, M1 overpredicts σ_ext, whereas ES7, BG7, and D7 underpredict σ_ext. J20, K30, USCII, and Z22 all show acceptable agreement with each other and with the experimental data, suggesting that these mechanisms can represent the effects associated with fluid-dynamic strain. The main difference between the detailed (J20, K30, USCII) and pathway-centric (Z22) mechanisms is the more detailed breakdown of H₂ for Z22 and the resulting larger radical pool, rendering notably improved predictions of τ_ign over a wide range of combustion conditions.

6. Examples of Supersonic Combustion LES

In this section we show results for a few selected high-speed combustor applications to demonstrate the relative importance of the different modeling aspects and to provide information about the cost and accuracy of LES for simple and advanced combustor configurations. This section also serves the additional purpose of demonstrating which type of information can be extricated from LES of high-speed combustion applications of varying complexity.

6.1. The Waidmann Supersonic Combustion Experiments

A case that has been used many times for validation of numerical simulation models for supersonic combustion is the experiment by Waidmann et al. [64,65,66,67], which also provides insight into supersonic flame stabilization behind a bluff body. In this case, preheated air is expanded through a Laval nozzle into a combustor section at Ma ≈ 2.0 and T = 340 K. The combustor section has a width of 40 mm, a height of 50 mm, and an upper-wall divergence angle of 3° to compensate for the boundary layer growth. A 32 × 6 mm² wedge-shaped strut is placed 100 mm downstream of the nozzle, and parts of combustor side-walls are made of quartz to allow optical access. H₂ at T = 250 K is injected sonically through fifteen 1.0 mm diameter injectors at the strut base. The equivalence ratio is φ ≈ 0.034, and the experimental studies included Schlieren imaging, OH Laser Induced Fluorescence (OH-LIF), Laser Doppler Velocimetry (LDV), and Coherent Anti Stokes Raman Scattering (CARS) to characterize velocity and temperature.

Many numerical simulations of this case have been presented using RANS [67], and hybrid RANS/LES and LES [7,68,69,70,71,72,73,74,75]. By comparing, it appears as if the agreement between predictions and experiments has improved over time, and particularly after switching from RANS to hybrid RANS/LES and LES. Simulations are now able to capture the dynamics of the shear layers and the wake. The studies presented here originate primarily from [75] and use hexahedral grids with ∼7.4 and ∼59.2 million cells, respectively, with refinement around the strut and towards the walls. According to [75], the finer grids do not significantly influence the results besides providing more details of the smallest explicitly resolved flow scales.

Figure 2a presents volumetric renderings of the density gradient, $\nabla \rho$, for non-reacting LES using the HV, SMG, and LDKM subgrid models. Included is also an experimental Schlieren image to facilitate qualitative comparison of the flow, mixing, and shock structures. Oblique shocks are formed at the tip of the wedge that reflect off the upper and lower walls before impinging on the H₂-filled wake. Together with the bent expansion fans, originating from the base of the wedge, this results in a shock-train. The upper and lower wall boundary layers are influenced by the reflected shocks, resulting in increased wall temperature and thickening of the boundary layer. The wedge boundary layer separates at the base, and shear-layers are formed that roll-up into Kelvin–Helmholtz vortices. Comparing the LES predictions with the Schlieren image shows good agreement for the shock-train and mixing pattern. An effect of the subgrid modeling is seen in the $\nabla \rho$ images, suggesting that the HV and LDKM models are more accurate than the SMG model. This is seen in the last crossing of the reflected shocks. Figure 2b compares cross-sectional profiles of the time-averaged axial velocity, $\langle v_x\rangle$, at five cross sections, x/h = 13.0, 20.83, 22.0, 26.16, and 34.5. Results from LES with the SMG, HV, and LDKM models are included together with RANS results based on the k-ω model. Only small differences are found in $\langle v_x\rangle$, but with the HV and LDKM models showing better agreement with the experimental data. These differences are, however, more pronounced for the rms axial velocity fluctuations, $v_x^{rms}$, as detailed in [75].

Figure 3 shows volumetric renderings of $\nabla \rho$ from the side together with the heat-release rate, $Q$, for six combinations of reaction mechanisms and combustion models. In the left column we compare the D7, J20, and Z22 reaction mechanisms using the LDKM turbulence model and PaSR combustion model, whereas in the right column we compare the PaSR, TFM, and ESF combustion models using the Z22 reaction mechanism and the LDKM turbulence model. Experimental Schlieren and OH-LIF images included are also included in Figure 3 to facilitate qualitative comparisons of the flow, mixing, and shock-train. When combustion is present, the expansion fans at the wedge corners almost disappear, and the recompression shocks become weaker compared to the non-reacting case in Figure 2a. Also, the recirculation region becomes longer and wider and acts as a virtual flame holder, and both the peaks reverse velocity and base pressure increase compared to the non-reacting case. The shear-layers from wedge corners are more noticeable in the reacting case, and ignition occurs in these as a result of shock impingement. Concerning the reaction mechanism, we find that the D7 mechanism ignition occurs in the shear-layers because of the too short ignition delay time, see Figure 1, resulting in too early and very intense heat-release. This also inhibits the formation of the Kelvin–Helmholtz vortices observed in the Schlieren and OH-LIF images, affecting the downstream flow. The J20 and Z22 mechanisms, on the other hand, result in delayed ignition and flame stabilization by recirculation of hot gases in agreement with the Schlieren and OH-LIF images. Finally, concerning the influence of the LES combustion model, we find very little difference between the predictions from the LES-ESF and LES-PaSR models, whereas the LES-TFM model results in a somewhat less intense combustion region.

Figure 4a,b shows profiles of $\langle v_x\rangle$ and the time-averaged temperature, $\langle T\rangle$, from the reacting cases at five cross sections, x/h = 13.0, 20.83, 22.0, 26.16, and 34.5. At the first location, all LES show reverse $\langle v_x\rangle$ on each side of the jets and overall good agreement with the experimental profile, although with the predicted jets being stronger than the measured ones. At the second location, the predicted $\langle v_x\rangle$ defect is narrower than the measured one and slightly asymmetric [63,67,70]. At the fourth and fifth locations, a virtually flat $\langle v_x\rangle$ profile is observed, and good agreement is obtained for all LES. At the first location, the temperature $\langle T\rangle$ in the shear-layers is overpredicted by D7, well predicted by J20 and slightly underpredicted by Z22. This is consistent with Figure 3 showing that D7 ignites too early, whereas J20 and Z22 ignite further downstream in the shear-layers allowing the flame to anchor in the wake of the strut. At the second location, the J20 and Z22 LES results show good agreement with the experimental profile, whereas the D7 results in a weak bi-modal $\langle T\rangle$ profile as a result of the intense burning in the shear-layers. At the third and fourth sections, no experimental data are available, and here we find that the LES results with Z22 result in wider $\langle T\rangle$ profiles compared to D7 and J20. At the fifth location, the LES results with the Z22 reaction mechanism show the best agreement with the experimental data. Regarding the sensitivity to the LES combustion model, we find that the LES-PaSR, EDC, TFM, and ESF models behave similarly, but with the LES-TFM model resulting in slightly delayed and narrower $\langle T\rangle$ profiles as a result of the lower heat release in the TFM model as seen by comparing Figure 3f with Figure 3d,e.

To understand more the flame structure and the flame stabilization mechanisms, Figure 5 shows images of the axial velocity, $v_x$, second invariant of the velocity gradient tensor, ${\lambda }_2$, colored by the rate-of-strain $\left|\tilde{\mathit{D}}\right|$, mixture fraction, $z$, reaction progress variable, $c$, extended Takeno Flame Index, TFI [76], and the explosive eigenvalue, ${\lambda }_e$, of the Chemical Explosion Mode Analysis, CEMA [77], overlaid on semi-transparent images of the density gradient, $|\nabla \rho |$. The Mach 2 flow through the combustor is responsible for the shock- and rarefaction waves, the shear-layers, and the recirculation region behind the strut. The vortical structure of the wake consists of interacting spanwise and longitudinal vortex structures as illustrated by λ₂, and high values of $\left|\tilde{\mathit{D}}\right|$ persist in the strut shear-layers. Comparing with Figure 2a, the strut-wake is affected by the volumetric expansion caused by the heat-release, $Q$, in the shear layers, Figure 3, due to the combustion of injected H₂. The reflected shock-waves interact with the wake shear-layers and the heat release regions downstream of the strut, promoting shock-induced combustion assisting the flame-holding process which otherwise is dominated by recirculation of hot combustion products from downstream as evident from the mixture fraction, $z$, and reaction progress variable, $c$, distributions. The highest value of $z$ is found in the H₂ jets, and in the recirculation region downstream of the strut, whereas the stoichiometric value, $z_{st}$, mainly follows the topological structure of $|\nabla \rho |$ and ${\lambda }_2$, being entrained in the vortex field on the edge of the wake. On the contrary, the highest value of $c$ is found in the wake, dominated by OH and H₂O. Based on $z$ and $c$, the flame lift-off distance, defined as the distance from the strut at which $c$ abruptly rises representing complete combustion, is about four strut-heights. Based on the TFI, the flame is a diffusion flame but with rich and lean premixed zones dominating in the recirculation region and along the flame edge, respectively. The CEMA explosive eigenvalue, ${\lambda }_e$, reveals positive values of ${\lambda }_e$ in the near-wake region, composed of partially reacted mixtures. Large positive values of ${\lambda }_e$ are found in the reacting parts of the shear-layers; large negative values of ${\lambda }_e$ are found in the upstream parts of the shear-layers, and intermediate values of ${\lambda }_e$ are found in the downstream parts of the shear layers.

6.2. The Fotia and Driscoll Direct Connect Facility Experiments

The Fotia and Driscoll direct connect facility experiments [78,79] offer an opportunity to study transition from ramjet to scramjet operation and to further validate the LES models. A constant-area isolator is fed by a two-dimensional Ma 2.2 nozzle, and H₂ is injected at the end of the isolator through a 2.49 mm diameter injector on the centerline, 44.5 mm upstream of the leading edge of a cavity with a slanted downstream wall. A 4° diverging area combustor section starts just after the cavity and discharges into a large-diameter exhaust. An electric heater combined with a hydrogen fueled vitiator was used to provide air stagnation temperatures between $T_0$ = 1040 K and 1500 K at stagnation pressures between $p_0$ = 420 kPa and 590 kPa. For the transition studies, $T_0$ and $p_0$ were fixed at 1400 K and 448.2 kPa, respectively, with the mode transition being facilitated by varying the H₂ inflow, resulting in global equivalence ratios ranging between $\varphi$ = 0.34 (ramjet mode) and 0.19 (scramjet mode), with the three intermediate equivalence ratios of $\varphi$ = 0.29, 0.26, and 0.22, representing different stages of the transition from ramjet to scramjet mode.

The computational set-up starts at the Laval nozzle and ends at the combustor exit. Hexahedral grids with ~11 million cells (baseline grid) and ~85 million cells (fine grid), both being refined at the walls and in the cavity, are used. Around 87% and 94% of the kinetic energy was resolved rendering both grids appropriate for LES.

Figure 6 presents perspective views of the LES predictions of (a) the ramjet Case F1 with $\varphi$ = 0.34, (b) the transition Case F3 with $\varphi$ = 0.26, and (c) the scramjet Case F5 with $\varphi$ = 0.19 in terms of contours of the refractive index gradient, $\nabla n$, iso-surfaces of the second invariant of the velocity gradient, λ₂ (gray), and volumetric renderings of the H₂ mass-fraction (green), and the temperature, $T$ (warm colors). Based on $\nabla n$, we find that the shock-train begins directly after the throat but then develops differently along the isolator depending on $\varphi$: for high values of $\varphi$, Case F1, the shock-train is short whereas for low values of $\varphi$, Case F5, the shock-train is longer. After this region, the shock-train gradually disintegrates due to the combined effects of shock-boundary layer interaction, shock–shock interactions, and downstream volumetric expansion due to exothermicity. As observed by comparing the $\nabla n$ distributions, the volumetric expansion increases with increasing $\varphi$ and is higher in the ramjet mode than in the scramjet mode.

The vortical flow, characterized by coherent λ₂ structures, appears to be initiated simultaneously in the upper and lower boundary layers of the isolator just after the end of the shock-train. This behavior suggests that the isolator boundary layer separates in this region, further enforcing the breakdown of the shock-train. The vortical flow then rapidly spreads across the full isolator cross-section, increasing the transverse mixing and simultaneously widening the boundary layer. Just upstream of the cavity, H₂ is injected perpendicular to the main flow direction. This jet-in-cross flow configuration enhances both the transverse and spanwise mixing. Downstream of the H₂ injection and along the cavity, the coherent λ₂ flow structures are modified by the vortices from the jet-in-crossflow, and the vortices shed from the cavity leading edge and the heat release from the combustion. An instantaneous image of λ₂ appears chaotic, but averaging λ₂ reveals a horseshoe vortex pair sweeping around the H₂ jet plume and being slightly diverted down into the cavity before breaking up near the end of the cavity. On average, the H₂ plume itself consists of a counter-rotating vortex pair and ring-like vortices surrounding the pair. In addition, boundary layer flow structures (longitudinal and hairpin vortices as well as streaks) also evolve after separation in the isolator, interacting with the plume vortex structures. The refractive index gradient, $\nabla n$, shows how the density changes as a consequence of the longer axial velocity, $v_x$, penetration, and reduced heat release with decreasing $\varphi$ from Case F1 to F5. Note that $\nabla n$ is influenced by the H₂ injection, the shock-train, as well as the volumetric expansion from heat release. Moreover, the effective plume angle decreases as the isolator flow speed increases due to the change in volumetric expansion over the cavity. Burning occurs in the cavity and along the borders of the H₂ plume. The differences between the λ₂ distributions are associated with the differences in the formation and break-up of the isolator shock-train, being pushed downstream with diminishing $\varphi$, and the heat-release which also decreases with diminishing $\varphi$. The temperature, $T$, increases along the isolator as a result of the increase in pressure, $p$, due to the blockage produced by the volumetric expansion from the heat release. After H₂ injection, a large increase in $T$ is observed around the H₂ plume and in the cavity. The cavity acts as a flame-holder since hot gases from the trailing edge of the cavity are re-circulated towards the leading edge where it heats up the shear layer above the cavity. The shear layer then interacts both with the H₂ rich plume and the hot air from the isolator surrounding the plume. The heat release, $Q$, occurs primarily on the outer edges of the H₂ rich plume but to some extent also in the trailing edge of the cavity where H₂ is injected.

Figure 7 compares time-averaged and rms wall pressures, ${\langle p\rangle }_{wall}$ and $p_{wall}^{rms}$, from experiment and LES. Angle brackets denote time averaging and the superscript ‘rms’ denotes rms fluctuations. Pressure acts as a marker of the combustion dynamics in a thermally choked flow since any change in $Q$ will result in a change in the shock-train and pressure rise. Reasonable to good agreement between LES and experimental data is found for ${\langle p\rangle }_{wall}$, in particular for the trends associated with $\varphi$, and for the transition from ramjet to scramjet mode. This implies that the mixing and combustion are better predicted for case 5 compared to case 1. The peak of $p_{wall}^{rms}$ occurs in the isolator due to the breakdown of the shock-train. Based on $p_{wall}^{rms}$, the ramjet case F1 is the steadiest since the reaction-zone is mainly located in the low-speed upstream part of the cavity shear-layer. The scramjet case F5 is also steady with the reaction zone moved further aft in the shear layer. The intermediate case F3, representing transition from ramjet to scramjet mode of operation, is the most unsteady case particularly in the isolator where the shock train dissolves.

6.3. The Small-Scale Flight Experiment (SSFE)

In the EU FP7 project LAPCAT [80] and the EU H2020 project STRATOFLY [81], numerical and experimental studies of a conceptual design for antipodal hypersonic passenger travel have been performed. In LAPCAT, a flight vehicle model was developed based on a waverider design, LAPCAT-II MR2, having a dorsal-mounted, dual-mode ramjet engine [82]. In STRATOFLY this model was further refined, resulting in the STRATOFLY MR3 model. For propulsion studies, a downscaled version of the LAPCAT-II MR2 model was developed by Langener et al. [83] to be used in the Small-Scale Flight Experiment (SSFE) at the DLR High Enthalpy piston-driven shock tunnel at Göttingen (HEG) [84]. Figure 8a shows a CAD image of the small scale STRATOFLY MR3 model inserted in the test section of the HEG. The wings are partially cut to fit the HEG test section, and the combustor is modified to compensate for some of the adverse scaling effects in subscale tests [83], by means of optimizing the combustor layout and increasing the equivalence ratio to $\varphi$ ≈ 1.0. With this set-up, free-flight tests were carried out at conditions corresponding to Ma 7.34 flight at an altitude of 28 km [85,86]. At these conditions $T$ = 264 K, $\rho$ = 0.0272 kg/m³, and $v$ = 2398 m/s. In the SSFE, the downscaled (1:60) STRATOFLY MR3 model is suspended in the test-section by two kevlar/aramide cords. Each cord is arranged in a V-configuration passing the SSFE model and is then attached on either side to the clamping mechanism. To allow free flight tests, the cords are cut exactly when the shock enters the test-section, allowing for ~4 ms of test-time before the model is caught by a cradle. The STRATOFLY MR3 model includes a H₂ injection system, pressure transducers, and a transparent upper wall enabling ultra-fast measurements with temperature sensitive paint or heat-flux measurements [87] and line-of-sight integrated species measurements using laser absorption spectroscopy [86]. Several combined, theoretical, numerical, and experimental studies have recently been published [88,89], providing a reasonable understanding of the complex flow and combustion physics.

Figure 8b shows the SSFE model used for simulations, consisting of an inlet, an isolator, a combustor, and an exhaust nozzle. Fuel injection is performed on the top of the semi-struts in the front of the combustor at the top and bottom of the central injection strut. The simulation encompasses the HEG test-section with the SSFE model. A hexahedral grid of ~34 million cells was used for the baseline case, and finer grid with ~272 million cells was used to study grid sensitivity. Both grids contain the same geometry, with low resolution in the freestream, but are refined along the engine flow path and towards the combustor walls and struts. Approximately 88% and 93% of the kinetic energy was resolved, rendering both grids acceptable.

Figure 9 shows exploded views from the side and top of the volumetric rendering of the (a) axial velocity, $v_x$, (b) temperature, $T$, (c), pressure, $p$, (d) heat-release rate, $Q$, and (e) species mass fractions, $Y_i$, from LES of the SSFE experiments. High axial velocities are found across the intake with the exception of a thin boundary layer adjacent to the intake wall. A counter-rotating vortex pair develops some distance downstream of the intake lip due to the laminar to turbulent transition and continues to develop until it reaches the combustor intake. Here, $v_x$ is relatively uniform with the exception of the boundary layer and the imprint of the counter-rotating vortex pair. The flow separates at the lower wall, between and just downstream of the semi-struts. This separation is shallow and starts at the junction between the downstream part of the intake and the upstream part of the combustor, following the break-up of the counter-rotating vortex pair from the intake. The H₂ injected from the semi-struts creates two large counter-rotating vortices that dominate the flow in the remainder of the combustor. The flow between the semi-struts and the full-strut appears very complex, dominated by the large counter-rotating vortices from the struts, the three-dimensional shock-train through the combustor, the boundary layers, and the small-scale turbulence developing from the interacting large-scale processes. Here, only limited volumetric expansion is observed since the H₂ is mixing and just about to ignite, primarily due to shock impingement. Just after the full strut, a significant decrease in $v_x$ is observed due to large volumetric expansion and increased pressure as a result of the H₂ being fully ignited, resulting in a substantial heat release. Further downstream and in the nozzle, $v_x$ recovers and increases to levels sufficiently high to provide thrust. The temperature, $T$, is low all along the intake, except at the walls where it increases slightly due to friction and in the combustor inlet where $T$ increases slightly due to compression. After the semi-struts, where H₂ is injected, $T$ increases once the injected H₂ has mixed with the compressed air from the intake. The ignition is due to shock compression and causes the edges of the large-scale vortices created by the semi-struts to burn intensely at high $T$. When additional H₂ is injected sideways from the full-strut, ignition occurs throughout the cross-section, and $T$ increases to around 2500 K across the whole cross-section. Further downstream and in the nozzle, $T$ gradually decreases due to the expansion in the combustor and nozzle where the flow is dominated by the large-scale vortices which generate mixing and small-scale turbulence.

The pressure, $p$, shows clear imprints of the three-dimensional shock train developing in the combustor, and particularly downstream of the full strut, where the effects of the large-scale vortical structures start to mutually interact with the pressure through baroclinic torque and vortex stretching, influencing mixing at all resolved scales. The high $p$ region ends where the constant area combustor changes into the diverging area combustor section, where the flow starts to accelerate. The shock-train is still present along the diverging section of the combustor and into the nozzle but it rapidly decreases in strength. The $p$-waves influence the chemical reaction rates which are enhanced in the constant area combustor, promoting heat release, $Q$, which occurs on the edges of the H₂ plumes once the H₂ has mixed properly with the ingested air. The majority of $Q$ is found just after H₂ has been injected sideways from the full strut, resulting in complete combustion across the combustor. The species distributions, H₂ in green, OH in magenta, H₂O in blue, and remaining radicals in orange, essentially follow the flow through the constant area and diverging area combustor sections and leave a feeble trace of unburnt H₂, OH, and radicals towards the upper wall. Most species are deeply embedded in the primary longitudinal vortex structures, interfered only by the shock trains and the chemical reactions. Due to the high turbulence levels, mixing is extremely efficient and promotes combustion.

Figure 10 shows images from (a) experiments and (b) LES of the wall heat flux on the lower combustor wall (left panel) and on the intake surface (right panel), as well as (c) the wall pressure along the upper combustor surfaces as indicated. Regarding the wall heat flux, the same colormaps are employed for the experimental and LES images on each section, but different colormaps are used for the upper combustor wall and intake surface. The measured wall heat flux image for the intake surface is obtained from [89] whereas those of the lower combustor wall are provided by DLR as part of STRATOFLY. For the intake, the wall surface heat flux increases downstream of the flat portion between the inlet lip and the downstream curved intake due to the compression of the flow. Imprints of the vortex pair close to the center plane and observed in Figure 10a,b can be seen in both the experimental and computational images. For the lower combustor wall, the differences between the experimental and computational results are larger, but the trends are similar with a region of high wall heat flux between the two semi struts and the full strut. Regarding the wall pressure, $p_w$, normalized by the free stream total pressure, $p_t$, in Figure 10c, results from the non-reacting and reacting cases are included together with the Rayleigh line. For the reacting case, two peaks at x ≈ 0.47 m and 0.65 m can be seen; these correspond to the pressure increase due to injection at the semi-struts and full struts, respectively. The increase in $p_w$ between x ≈ 0.55 m and 0.90 m is due to the volumetric expansion caused by combustion and is reasonably reproduced by the LES. Both the experimental and computed peak pressures exceed the maximum theoretical pressure rise due to H₂-air combustion, as indicated by the Rayleigh line. This suggests the presence of a strong shock-train throughout the combustor as seen in Figure 9c.

Figure 11 compares cross sectional profiles of the time-averaged axial velocity, $\langle v_x\rangle$ temperature, $\langle T\rangle$, and mass fractions $\langle Y_{H2}\rangle$, $\langle Y_{OH}\rangle$, $\langle Y_{H2O}\rangle$, and $\langle Y_{NO}\rangle$ at x = 0.395, 0.550, 0.700, 0.800, 0.900, and 1.100 m along the combustor. The highest $\langle v_x\rangle$ occurs over the intake and the nozzle whereas lower $\langle v_x\rangle$ prevail in the combustor as a consequence of the complex flow pattern and the volumetric expansion due to the exothermicity. The wakes of the semi struts appear clearly in $\langle v_x\rangle$ at x = 0.550 m and of the full strut at x = 0.700 m, after which the $\langle v_x\rangle$ deficit has disappeared due to the volumetric expansion that slowly increases $\langle v_x\rangle$ further downstream. $\langle T\rangle$ is close to the freestream value over the intake, x = 0.395 m, with some increase observed in the boundary layers due to friction heating. Just downstream of the semi struts, at x = 0.550 m, $\langle T\rangle$ increases along the edges of the injected $\langle Y_{H2}\rangle$ plume as it ignites. This is particularly evident on the outer edges of the $\langle Y_{H2}\rangle$ plumes. Further downstream, at x = 0.550 m to 0.900 m, the full combustor cross-section is ignited, resulting in high $\langle T\rangle$ across the full engine duct. Further downstream, at x = 1.100 m, $\langle T\rangle$ decreases somewhat due to the widening of the nozzle and the flow acceleration. The imprints of the $\langle Y_{H2}\rangle$ injection from the two semi struts and the combustion occurring on the edges of the $\langle Y_{H2}\rangle$ plumes from two semi-struts can easily be seen. The rapid spread of combustion across the cross-section is also evident just downstream of the full strut. $\langle Y_{OH}\rangle$ is formed where combustion, heat release, and volumetric expansion occur, and high values of $\langle Y_{OH}\rangle$ are found at the boundaries of the $\langle Y_{H2}\rangle$ plumes, at x = 0.550 m, and just after the full strut, x = 0.700 m, where the sideways injected H₂ promotes complete combustion across the full width of the engine duct. $\langle Y_{OH}\rangle$ gradually diminishes further downstream, at x = 0.80 m and 0.90 m as combustion completes at x = 1.10 m, where $\langle Y_{OH}\rangle$ is zero. $\langle Y_{H2O}\rangle$ is non-existing over the intake but forms as soon as combustion has started on the edges of the $\langle Y_{H2}\rangle$ plumes at x = 0.550 m and continues to form with downstream distance from H₂ injection. Based on $\langle Y_{H2O}\rangle$, complete combustion has occurred just after x = 0.800 m, where the computed $\langle Y_{H2O}\rangle$ profiles are virtually flat across the engine duct. Comparing with experimentally measured line-of-sight integrated $\langle Y_{H2O}\rangle$ values, excellent agreement is observed between x = 0.550 m and 1.100 m. Regarding NO, very good agreement between the experimentally measured line-of-sight integrated $\langle Y_{NO}\rangle$ values and the computed $\langle Y_{NO}\rangle$ values is observed. This demonstrates reciprocal validation of the LES computations and the experiments, increasing the value of the analysis.

7. Conclusions

Predictive modeling of high-speed combustion is becoming increasingly important for the development of novel airbreathing propulsion systems such as ramjet, scramjet, dual-mode ramjet, and rotating detonation engines. However, the interaction of turbulence with physical processes, such as jets, sprays, vaporization, chemical kinetics, and convective and radiative heat-transfer, is a great challenge. The increase in computational capacity in the past decades has made many of these flows accessible to numerical simulation, and major advances in modeling turbulent flow have now been made as we are entering the era of sufficient computer power for Large Eddy Simulation (LES) and similar methods. The definition of LES is to numerically solve the spatially filtered flow equations, whereby the large (energetic) scales of the flow, which are directly affected by boundary conditions, are fully resolved while the small (less energetic) scales are subject to sub-filter modelling. Restricting the difficult sub-filter modelling problem to only the smallest scales does provide a clear improvement over traditional RANS methods, but it also comes at a significant computational cost because the solutions to the transient spatially filtered equations need to be sampled for a long time to obtain converged statistics.

Today, LES modelling of high-speed combustion has already come a long way, and there are many examples of successful simulations in the literature. These range from simple canonical set-ups to large complex geometries as exemplified in this chapter. Some of the advances include the ability to predict, for example, shock-trains, boundary layers, shock-induced ignition, jet in cross-flow with separation shocks and complex vortical structures, thermal choking, acoustic feedback, extinction events, and detonation waves. However, for historical and practical reasons, the subgrid models used in LES are still based on the ideas inherited from RANS models, and the improved predictions that are observed in LES are mainly due to resolving the large turbulence scales. The subgrid models themselves still provide some limitations, and improvements will be necessary to meet the future demands on accuracy, speed, and reliability of predictions in practical applications. In addition, some of the additional physics involved in more complex engineering simulations may not be described with enough detail by commonly used models.

In the authors’ opinion, the most critical areas that need to be improved upon for high-speed combustion LES are the following: (i) Models for turbulent stresses and transport terms, specifically adapted for compressible flows. The accuracy of these models needs to be investigated when there are shock–shock, shock–boundary-layer, and shock–interface interactions present in the flow. Here there is also a need for accurate DNS simulations for validation of any proposed closure model. While the exact influence of turbulence on high-speed combustion in general can be debated, it is expected to be of high importance, at least in the boundary layers and for the mixing in jet shear layers and in cavities. (ii) Models for the filtered reaction source terms in compressible chemically reactive flows. The unresolved reaction source terms are a major problem of reacting RANS, and it is unclear if the additional resolution of LES provides any improvement here since the reactions take place on the smallest scales which are still not resolved. While many models exist, there is still no consensus of a generally applicable accurate model. Furthermore, the development of reaction model closures will be tightly coupled with the exploration of how high-speed turbulence affects combustion. (iii) Models of liquid sprays including break-up and vaporization. This is a notoriously difficult problem already in incompressible flow, and for airbreathing high-speed engines, there will likely be interactions between droplets and shocks, expansion fans, and large drag forces. (iv) Thermodynamics and transport phenomena in high-speed, high-temperature flows, including deflagration-to-detonation transition, thermal radiation, non-equilibrium thermodynamics behind shocks, ionization, dissociation, and excitation of chemical species.

As indicated, there is a large range of important topics that may be targeted in future work. In terms of LES of high-speed combustion, it is suggested to perform wide studies where a combustor is studied in a wide range of operating conditions (fuels, Ma numbers, equivalence ratios, temperatures, injection velocities, etc.) to further probe the design space. Secondly, there are several phenomena of high-speed combustion that deserve further focused studies, including the effect of turbulence on high-speed flames and detonations, thermoacoustic instabilities in supersonic combustion chambers, and shock interactions with turbulent boundary layers. Third, the different communities using and developing LES should start to engage in more cross-disciplinary activities for development of concept designs as well as systematic validation of models. The future of LES for high-speed combustion applications is promising, in particular if the LES results are analyzed together with experimental data. This, in turn, puts pressure on the development of both experimental facilities and diagnostics, redirecting the experimental work towards fewer but more carefully instrumented experiments. Only by combining experimental and simulation data can the next generation of airbreathing engines be developed to the standard required.