Direct numerical simulation of a supersonic turbulent boundary layer with hydrogen combustion

Abstract

The complex behaviours of supersonic turbulent boundary-layer flows interacting with combustion is explored through the use of direct numerical simulations. The chosen flow model is a non-premixed hydrogen–air flame ignited within a three-dimensional supersonic turbulent flat-plate boundary layer operating at a Mach number of 2.33 and a friction Reynolds number of approximately 1000. The simulation involves a finite-rate model for the hydrogen–air reaction including 9 species and 19 steps with adiabatic non-catalytic wall conditions. The inlet flow consists of preheated air in the main stream and a hydrogen film injected in proximity to the cold wall, inducing mixing and ignition in the outer layer. The specific inlet configuration leads to two successive transition stages of distinct mechanisms, with the first stage related to the mixing-layer instability and the second one to the boundary-layer instability. The excess Reynolds stresses resulting from the transition exhibit a downstream decrease, with values being lower than canonical adiabatic profiles, particularly in the outer layer. This difference is attributed to the reduced local friction Reynolds numbers, associated with the non-classical wall-normal distributions of density and viscosity. The effect of combustion on the recovery process and skin friction is analysed in comparison with a supplemented chemistry-frozen setting. Velocity–temperature and velocity–species correlations are further examined. The wall-normal profiles of turbulent Prandtl number tend to classical non-reacting values and the turbulent Schmidt number is only slightly affected by the reaction. The strong Reynolds analogy linking velocity and temperature fluctuations is found to be invalid in the outer layer due to the presence of large-scale temperature fluctuations and the inhomogeneity of hydrogen gas distribution, as revealed by a spectrum analysis. A statistical analysis of elementary reactions at varying wall-normal distances is provided, highlighting the dominance of hydrogen atom depletion in the inner region and the prevalence of water vapour production in the outer region, together contributing to the chemical heat release across the boundary layer. Turbulence–chemistry interaction is assessed through a comparison of mean turbulent and laminar chemical heat release rates, showing significant differences of up to an order of magnitude. Using a simple diagnostic tool, the modelling of fluctuations associated with hydrogen gas and two intermediate products is found to be crucial for turbulent chemistry closure.

1. Introduction

Understanding and modelling turbulent flow dynamics in compressible boundary layers is crucial for advancing fundamental studies in the design of high-speed flight vehicles. In specific scenarios, finite-rate chemical reactions play a significant role, affecting both external and internal flows of vehicles. These reactions result in heat release or absorption, altering gas composition and transport properties, thereby influencing the dynamics of bounded flows. Three primary categories of chemical reactions are commonly explored in high-speed boundary-layer flows. The first involves thermochemical non-equilibrium effects related to high-enthalpy gas dissociation, occurring prominently at the blunt body's leading edge where the generated shock wave is highly compressed (Park 1993; Candler 2019). The second category involves leveraging carbon dioxide to delay turbulent boundary-layer transition by absorbing energy from acoustic disturbances (Johnson & Candler 1998; Leyva et al. 2009). The third category focuses on the supersonic combustion process near combustor wall surfaces within air-breathing propulsion systems, which is the focus of this paper.

1.1. Experimental and numerical studies of near-wall supersonic combustion

The ignition and sustainability of combustion play a critical role in designing air-breathing propulsion systems for high-speed vehicles (Urzay 2018). Given the millisecond-scale residence time of reactants within the combustion chamber, various fuel injection schemes are designed to improve the fuel–air mixing process and combustion stability (Seleznev, Surzhikov & Shang 2019). These schemes include steps (Burrows & Kurkov 1973; Suraweera, Mee & Stalker 2005), cavities (Voland et al. 1999; Micka & Driscoll 2009; Storch et al. 2011), pylons (Anderson & Gooderum 1974; Waidmann et al. 1994; Boyce et al. 2000) and throttles (Zabaikin 2003; Tretyakov 2012). Apart from pylons enabling fuel injection into the main gas flow outside the boundary layer, the other schemes mainly employ fuel injectors implemented through the wall surface of designed combustor structures. In these schemes, complex dynamics of turbulence–flame interaction occurs in supersonic turbulent boundary layers. For instance, in schemes with cavities, multiple fuel injectors are implemented inside or at the upstream and downstream of the cavity, allowing enhanced mixing and flame stabilisation through the recirculation bubble generated therein.

While experimental databases of supersonic combustors are invaluable for numerical simulation validation, obtaining velocity measurement data near the wall is challenging in supersonic flows (Seleznev et al. 2019). One such experimental set-up, the Burrows–Kurkov combustor (Burrows & Kurkov 1973), has been investigated extensively through joint numerical simulations. In the experiment, preheated air at Mach 2.44 was injected into the main stream of a backward-step combustor, where it mixed with a sonic hydrogen fuel stream horizontally injected from a vertical slot at the step. The air–hydrogen mixture downstream along the combustor wall surface, at a slight angle, led to ignition in the boundary-layer flow. The experimental database includes temperature and species concentration at the combustor exit, along with wall pressure measurements and ignition location records.

The Burrows–Kurkov combustor serves as a benchmark for validating numerical approaches, including Reynolds-averaged Navier–Stokes (RANS) equations (Drummond, Rogers & Hussaini 1987; Eklund & Stouffer 1994; Tretyakov 2012; Gao et al. 2015; Zuo et al. 2020) and large-eddy simulation (LES) (Edwards, Boles & Baurle 2012; Vyasaprasath et al. 2015; Liu et al. 2020; Wei et al. 2023). Edwards et al. (2012) employed a hybrid LES/RANS method with a mixing eddy-viscosity model. Simulations covered both two-dimensional (2-D) spanwise-periodic grids and three-dimensional (3-D) grids, including top and side wall surfaces. Laminar chemistry or Gaussian quadrature assumptions (Donde, Koo & Raman 2012) were employed for closing filtered chemical reaction rates. The simulations provided good agreement for flame anchoring positions and time-averaged profiles at the combustor exit, with minor deviations in the peak values of temperature and water vapour concentration profiles. The choice between seven and nine species for chemical reactions showed no significant differences, and the use of Gaussian quadrature closure slightly increased the peak values. Recent LES studies (Gao et al. 2015; Wei et al. 2023) with 3-D grids and spanwise periodic conditions employed seven-species reaction models. Both studies achieved good agreement with experimental data regarding time-averaged profiles at the combustor exit, although Gao et al. (2015) observed increased skin friction due to combustion, whereas Wei et al. (2023) noted reduced skin friction.

The effect of combustion on supersonic skin friction is a key focus of a similar supersonic combustor experiment conducted by Suraweera et al. (2005). The combustor, analogous to the Burrows–Kurkov combustor, features a step configuration without inclination of the wall surface. With a main-stream Mach number of 4.5, the results showed a 70–80 % drag reduction compared with no hydrogen stream injection. Prediction models for skin friction, considering chemical heat release (Stalker 2005; Barth, Wheatley & Smart 2013; Liu et al. 2017b), were proposed based on the von Kármán momentum integral of boundary layers. The drag reduction mechanism was attributed to a decrease in Reynolds stresses and thickening of the turbulent boundary layer due to combustion (Barth et al. 2013). Note that even with combustion, the heat flux along the wall surface remained at a similar order of magnitude as in cases without hydrogen injection. However, this observation may be limited by the experimental operation time of approximately 3 ms, during which the wall surface temperature was maintained at an ambient temperature of 300 K. In a real flight condition with an extended operation time, boundary-layer combustion might introduce an additional heat source, presenting potential challenges for wall surface heat protection. Therefore, although the concept of using boundary-layer combustion for skin-friction control has been explored (Clark & Bade Shrestha 2014; Xue et al. 2020, 2021; Lu, Zhang & Qin 2024; Qu et al. 2024), the focus of the present work is not on drag reduction as a control technique. Instead, in this paper we study the physical mechanism of how combustion affects the skin friction.

Validating the different skin-friction observations in the 3-D LES calculations of Gao et al. (2015) and Wei et al. (2023) is challenging due to the absence of skin-friction measurements in the Burrows–Kurkov experiment. Those differences in the calculation results are already expected, as the near-wall flow may not be adequately resolved in LES, and different sub-grid models could lead to distinct skin-friction estimations. Additional investigations into different modelling strategies in a reacting setting are necessary, for example, the closures for turbulent heat-flux and diffusion terms, which often rely on standard gradient-diffusion formulations, assuming homogeneous values for turbulent Prandtl number $Pr_t$ and turbulent Schmidt number $Sc_t$. Those assumptions have rarely been assessed in a high-speed near-wall direct numerical simulation (DNS) data with combustion.

1.2. DNS of high-speed boundary layers with chemical reactions

To the best of the authors’ knowledge, no DNS has been conducted on the Burrows–Kurkov or similar supersonic combustors that fully resolves near-wall turbulence. Our preliminary calculations estimate a friction Reynolds number $Re_\tau$ in the Burrows–Kurkov experiment exceeding 10 000, significantly surpassing existing DNS calculations of compressible boundary layers (Pirozzoli & Bernardini 2011; Wenzel et al. 2018; Zhang, Duan & Choudhari 2018; Cogo et al. 2022). The high $Re_\tau$ can be attributed, in part, to the isothermal cooled wall condition in both the experiments and LES calculations. The downstream wall of the backward step is treated as isothermal at 300 K, resulting in a recovery temperature ratio of approximately 0.09 relative to the preheated main stream. This value is notably lower than the recovery temperature ratio in the range from 0.15 to 0.25 used in simulations of compressible boundary layers with cold walls (Duan, Beekman & Martin 2010; Zhang, Duan & Choudhari 2017; Zhang et al. 2022; Xu, Wang & Chen 2023; Yu et al. 2024).

In closely related fields, Martin et al. (1998) and Martin & Candler (2001) conducted pioneering work simulating Mach 4 turbulent boundary layers coupled with one-step reversible reaction models involving both endothermic and exothermic reactions. Their findings showed that endothermic reactions reduce temperature fluctuation magnitudes, while exothermic reactions increase them. Duan & Martín (2011) simulated turbulent boundary layers at Mach 3.4 and 10.3 using a five-species and five-reaction mechanism for air dissociation (Park 1993). The study focused on assessing turbulence–chemistry interaction in high-speed boundary layers, quantified by the difference between turbulent and laminar reaction rates. A priori tests indicated that the turbulence–chemistry interaction might not strongly impact the turbulent flow field, as suggested by interaction Damköhler numbers at least one order smaller than unity. This conclusion was further supported by a posteriori tests, artificially replacing the instantaneous reaction rate with the laminar reaction rate in DNS calculations. The authors noted that, in contrast to air dissociation, turbulence–chemistry interaction in combustion flows could play a significant role in chemical reaction and flow statistics, which can be attributed to the influence of turbulent fluctuations of temperature and species on radical reactions in combustion. This interaction is often modelled through stochastic methods in mixing layers and free shear flows (Baurle, Hsu & Hassan 1995; Baurle & Girimaji 2003; Koo, Donde & Raman 2011; Donde et al. 2012).

There has been a notable increase in studies over the past 5 years focusing on the DNS of high-speed reacting boundary-layer flows, particularly addressing thermochemical non-equilibrium effects associated with gas dissociation (Candler 2019). Di Renzo & Urzay (2021) simulated a Mach 10 high-enthalpy boundary layer under isothermal cold wall conditions, revealing a non-monotonic wall-normal temperature profile with a peak value approximately four times induced by viscous heating. Air dissociation occurred at the temperature peak region, corresponding to an endothermic reaction, whereas recombination prevailed closer to the wall surface. Passiatore and colleagues conducted a series of studies on hypersonic boundary-layer simulations assessing thermochemical non-equilibrium effects (Passiatore et al. 2021, 2022, 2023). In simulations with an adiabatic hot wall (Passiatore et al. 2021), air dissociation occurred at the wall boundary where the temperature peaked. Comparison with artificially chemically frozen simulations revealed that the endothermic reaction slightly reduced mean temperature, kinetic energy and temperature fluctuation but had minimal effect on skin friction in regions of fully developed turbulence. In another study employing a two-temperature model for a cold-wall boundary layer at Mach 12.48 (Passiatore et al. 2022), turbulence transport was found to redistribute gas to a vibrationally under-excited state near the wall and an over-excited state farther away. Classical closure strategies for sub-grid heat flux, using turbulent Prandtl and the strong Reynolds analogy (SRA) and its variances, were extended to vibrational temperature fluctuations, demonstrating validity similar to translational temperature. Recent studies also involved simulations of reacting boundary layers with impinged shock waves (Volpiani 2021; Passiatore et al. 2023). In addition, simulations of Mach 4.5 turbulent boundary layers were conducted across a range of enthalpies to examine the air dissociation effect on wall heat transfer (Li et al. 2022).

1.3. Contributions for the present study

The literature review reveals that supersonic combustor calculations involving boundary-layer combustion have primarily employed RANS or LES methods. Although these calculations are tested against experimental databases, such as the well-known Burrows–Kurkov combustor, validation against near-wall measurements remain scarce. Conducting DNS of supersonic boundary-layer combustion can be valuable for the community, providing insights into the turbulence dynamics in this complex scenario and offering a database to assess closure models for turbulence and chemistry. In a closely related field, an increasing number of DNS studies over the past 5 years have focused on air dissociation in high-enthalpy, high-speed turbulent boundary layers. These simulations contribute valuable insights into physical processes such as turbulence–chemistry interactions and offer numerical support for validating closure models.

In this study, we perform DNS of a non-premixed hydrogen–air flame ignited in a 3-D turbulent flat-plate boundary layer with a main-stream Mach number of 2.33. A model wall jet profile is specified at the laminar inlet, featuring a Poiseuille flow profile close to the wall for the sonic hydrogen stream and a similar boundary-layer profile for the preheated air main stream farther away from the wall surface. These profiles are connected with a no-slip velocity at the interface, analogous to an infinitely thin nozzle outlet separating the jet stream and the main stream. The prescribed inlet profile induces flow transition (Liu, Wang & Piao 2017a), leading to turbulence transition downstream, without resorting to other laminar-to-turbulent transition strategies such as blow and suction (Pirozzoli, Grasso & Gatski 2004; Di Renzo & Urzay 2021; Passiatore et al. 2022). The mixing of the hydrogen and air stream induces auto-ignition and sustained turbulent combustion in the boundary-layer flow. A finite-rate model with 9 species and 19 reactions is employed for hydrogen–air combustion.

The wall–jet inlet configuration resembles that used in Burrows–Kurkov combustors, and the flow parameters have been adjusted from the original experimental values to ensure that the simulated turbulence aligns with DNS requirements (Poggie, Bisek & Gosse 2015). Nevertheless, we note that the objective of the present calculation is not to reproduce the results of the Burrows–Kurkov combustor. Instead, the described model captures some common characteristics of boundary-layer combustion in high-speed flows. These include flame propagation at the downstream edge of cavities (Storch et al. 2011; Lin et al. 2023) and combustion induced by injectors inclined with the wall surface (Belanger & Hornung 1992; Lin et al. 2010; Bao et al. 2014). In all these experimental configurations, there is the presence of a flame interacting with high-speed near-wall turbulence, a phenomenon that has been scarcely investigated using DNS.

This paper is the first step of our exploration into this DNS database, with an emphasis on the flow statistics in the transition and recovery processes, the correlations of velocity, temperature and species fluctuations, and the chemical statistics at various wall-normal distances. The structure of the paper is as follows. In § 2, the computational set-up is described. In § 3, 2-D and 3-D visualisations of turbulence flow fields are provided. In § 4, the transition behaviour associated with the specific wall–jet set-up is analysed, followed by an examination of the recovery process in comparison with the canonical turbulent boundary-layer flows. In § 5, we evaluate the turbulent Prandtl number, the turbulent Schmidt numbers and the SRA. Section 6 provides statistics on chemical reactions and the turbulent-chemistry modelling for RANS calculations. Conclusions and perspectives are given in § 7.

2. Numerical settings

2.1. Governing equations

The flow equations and chemical reaction models employed in this study are precisely the same as those employed by Yan et al. (2022), using the OpenCFD-Comb code. Therefore, we present the governing equations in a concise form, and interested readers can refer to Yan et al. (2022) for more comprehensive details.

The compressible Navier–Stokes equations for multi-species reacting flows are formulated in 3-D coordinates $(x, y, z)$, where the components represent the streamwise, wall-normal and spanwise directions, respectively. The conservative variables are chosen as ($\rho, \rho \boldsymbol {u}, \rho E, \rho _i$), where $\rho$ is the density, $\boldsymbol {u}=(u, v, w)$ denotes the streamwise, wall-normal and spanwise velocity components, $E$ represents the total energy per unit volume and $\rho _i=\rho Y_{i}$ is the density associated with the $i$th species. Here, $Y_{i}$ denotes the mass fraction of the $i$th species. The conservative equations are expressed as

(2.1)$$\begin{gather} \frac{\partial \rho}{\partial t} +\boldsymbol{\nabla}\boldsymbol{\cdot}\left (\rho \boldsymbol{u} \right ) =0, \end{gather}$$

(2.2)$$\begin{gather}\frac{\partial \rho\boldsymbol{u}}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} \left (\rho\boldsymbol{u} \boldsymbol{u} \right )={-} \boldsymbol{\nabla} p + \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\tau}, \end{gather}$$

(2.3)$$\begin{gather}\frac{\partial \rho E}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} \left ( \rho \left (E+p \right ) \boldsymbol{u} \right ) = \boldsymbol{\nabla} \left (\boldsymbol{\tau} \boldsymbol{\cdot} \boldsymbol{u} + \boldsymbol{q} \right ) + \dot{\omega}_T, \end{gather}$$

(2.4)$$\begin{gather}\frac{\partial \rho_i}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} \left ( \rho Y_i \boldsymbol{u} \right )= \boldsymbol{\nabla} \boldsymbol{\cdot}\left (\rho D_i \boldsymbol{\nabla} Y_i \right ) + \dot{\omega}_i. \end{gather}$$

The equations are closed by the Dalton's law, expressed as

(2.5)\begin{equation} p = \sum_{i=1}^{N} \rho_i \frac{R_0}{M_i} T, \end{equation}

where $p$ denotes the pressure, $T$ is the temperature, $R_0$ represents the universal gas constant and $M_i$ is the molar mass of the $i$th species. The molecular stress tensor $\boldsymbol {\tau }$ and the conductive heat flux $\boldsymbol {q}$ are given by

(2.6)$$\begin{gather} \boldsymbol{\tau} =\mu \left ( \boldsymbol{\nabla}\boldsymbol{u} + \boldsymbol{\nabla}\boldsymbol{u}^{\mathrm{T}} - \tfrac{2}{3} \left (\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} \right ) \boldsymbol{I} \right ), \end{gather}$$

(2.7)$$\begin{gather}\boldsymbol{q}=\kappa \boldsymbol{\nabla} T + \rho \sum_{i=1}^{N} D_i h_i \boldsymbol{\nabla} Y_i. \end{gather}$$

The thermodynamic variables, including specific heat capacity $C_p$, enthalpy $h$ and entropy $S$, and the transport coefficients, such as molecular viscosity $\mu$, thermal conductivity $\kappa$ and binary diffusion coefficients $D_i$ of the $i$th species, are evaluated based on average models that depend on the temperature and mass fraction of each species (Yan et al. 2022).

A 9-species, 19-step chemistry model from Li et al. (2004) is employed for hydrogen–air reaction. We denote the number of species and reactions as $N=9$ and $M=19$, respectively. The species involved are $\mathrm {O}$, $\mathrm {O}_{2}$, $\mathrm {H}_{2}$, $\mathrm {H}_2 {\mathrm {O}}$, $\mathrm {OH}$, $\mathrm {H}$, $\mathrm {HO}_{2}$, $\mathrm {H}_2 \mathrm {O}_{2}$ and $\mathrm {N}_{2}$. In practice, the continuity equation (2.1) and the species transport equations (2.4) associated with the first eight species are solved, whereas the mass fraction of $\mathrm {N}_{2}$ is obtained as $Y_{\mathrm {N}_{2}}=1-\sum _{i=1}^{i=8} Y_{i}$. The contribution of $j$th reaction step to $i$th species is computed as

(2.8)\begin{equation} \dot{\omega}_{i,j}= \left ( \chi'_{i,j}- \chi_{i,j} \right ) \phi_{j}, \end{equation}

where the elementary reaction rate of $j$th step is expressed as $\phi _{j}=\varphi _{j}- \varphi '_{j}$. The forward and backward reaction rates of $j$th step are given by $\varphi _{j}$ and $\varphi '_{j}$, and their associated stoichiometric coefficients are denoted as $\chi _{i,j}$ and $\chi '_{i,j}$, respectively. The production rate of the $i$th species $\dot {\omega }_i$ in (2.4) is given by

(2.9)\begin{equation} \dot{\omega}_i= W_i \sum_{j=1}^{M} \dot{\omega}_{i,j}. \end{equation}

The heat release rate in (2.3) is given by

(2.10)\begin{equation} \dot{\omega}_T={-}\sum_{i=1}^{N} \Delta h_{f,i}^o \dot{\omega}_i, \end{equation}

where $\Delta h_{f}^o$ is the standard formation enthalpy corresponding to the $i$th species.

2.2. Numerical methods

Calculations are conducted using the OpenCFD-Comb code, an enhanced version of the compressible DNS code OpenCFD that integrates equations for multi-species reacting flows. The original OpenCFD code has been exclusively employed in the investigation of compressible turbulent boundary layers over the past two decades (Li, Fu & Ma 2006; Li et al. 2010; Zhang et al. 2012; Zhu et al. 2017; Xu et al. 2021, 2023; Zhao et al. 2024). The reacting version has been applied in recent studies to simulate a supersonic jet flame (Fu et al. 2019) and a turbulent mixing process involving the Richtmyer–Meshkov instability with the presence of hydrogen combustion (Yan et al. 2022). Both studies used the 9-species, 19-step model from Li et al. (2004) for hydrogen–air reaction, the same chemical model employed in the present study. The convection terms in the governing equations are discretised using an optimised six-order monotonicity-preserving scheme (OMP6) (Li, Leng & He 2013), consistently applied in the aforementioned studies. An eighth-order centred difference scheme is employed to discretise the viscous terms. Time-stepping is carried out using a third-order Runge–Kutta method, with the time stepping of chemical reaction terms performed through the Strang splitting scheme (Strang 1968; Ren & Pope 2008; Fu et al. 2019; Yan et al. 2022).

2.3. Flow configuration

The simulated flow involves a non-premixed wall–jet system consisting of a hydrogen stream and a vitiated air stream with water vapour. The cold hydrogen stream is injected close to the wall at Mach 1, and the preheated air stream is injected as the main stream at Mach 2.33. The flow configuration is illustrated in figure 1, and the flow parameters are detailed in table 1. Inlet profiles of streamwise velocity, temperature and species concentrations are also shown in figure 2. In both streams, species concentrations are assumed to be homogeneous. The hydrogen and air stream profiles are connected at a wall-normal distance $h_0=4$ mm, where $h_0$ represents the height of the hydrogen stream. The interface at $h_0=4$ mm corresponds to an infinitely thin injection nozzle, where the velocity is imposed with a no-slip condition and the temperature is set to 300 K. A Poiseuille profile is prescribed for the hydrogen stream inlet, with the temperature reaching a maximum of 300 K at the wall $y=0$ and at the stream interface $y=h_0$. The minimum temperature in the profile, $T_c=254$ K, is prescribed at the jet centre $y=h_0/2$ (Peter & Kloker 2022). The main stream is set as a self-similar solution of the compressible Blasius boundary layer using the ideal gas law without chemical reaction. The prescribed boundary-layer thickness is $\delta _{m,0}=3.2$ mm. Characteristic non-reflecting conditions are applied at the inflow, outflow, and free-stream boundaries (Poinsot & Lelef 1992). The wall surface is treated as no-slip, adiabatic and non-catalytic. The last two conditions are interpreted as $\partial T / \partial y =0$ and $\partial Y_i /\partial y =0$ at the wall surfaces, allowing variations of temperature and species concentrations in the streamwise and spanwise directions.

Figure 1. Flow configuration.

Table 1. Parameters of inlet profile.

Figure 2. Inlet profiles of velocity and temperature (a) and species concentrations (b). The entire calculation domain extends to $y^*=y/h_0=22$ and the constant values at the outer flow regions are not displayed.

Two major differences exist in the simulated flow configuration compared with the Burrows–Kurkov experiments and corresponding LES calculations. The first is the use of a laminar model wall–jet inlet, in contrast to the turbulent air main stream in the experiments. The exact inlet condition for the turbulent main stream in the Burrows–Kurkov set-up is unknown, and key results, such as the flame anchoring position, are reported to be highly sensitive to the thickness of the inlet turbulent boundary layer (Edwards et al. 2012; Liu et al. 2020; Wei et al. 2023). In these studies, various turbulent inflow generation methods have been used, which can lead to potential differences in results. Using a laminar inlet can prevent this problem and enhance the reproducibility by future studies. At the same time, the use of a laminar inlet can reduce the computational costs by eliminating the need for a precursor domain to generate turbulence. Although the laminar inlet may result in slower hydrogen–air mixing compared with the turbulent air main stream, adjustments to flow parameters, such as reducing main-stream boundary-layer thickness and increasing main-stream temperature, have been made to ensure turbulent mixing and ignition. The effective calculation length is also extended to one half downstream with respect to the Burrows–Kurkov combustor exit.

The second distinction involves the use of an adiabatic wall condition instead of the isothermal cold wall in the experiment and the LES calculations. This choice is motivated by the necessity to reduce the friction Reynolds number to around 1000, making DNS feasible. As mentioned in a previous section, the cold wall condition with a recovery temperature ratio of 0.09 leads to strong shear near the wall surface associated with high friction Reynolds number. The adiabatic wall condition is implemented by imposing a Neumann condition for the temperature at the wall surface, allowing its evolution in the streamwise and spanwise directions. This set-up enables temperature fluctuations at the wall and facilitates the quantification of the chemical heat introduced by combustion, in comparison with a chemistry-frozen case. Note that this approach has been applied in simulating film cooling experiments to assess the efficiency of cooling along streamwise directions in a worst-case heat transfer scenario (Peter & Kloker 2022).

The size of the effective computation domain is $L_x \times L_y \times L_z = 150h_0 \times 22h_0 \times 10h_0$. The streamwise direction is further extended downstream to $300h_0$, incorporating a sponge layer with coarse grids to dampen reflections at the outlet boundary. The complete calculation domain is discretised with $N_x \times N_y \times N_z = 2206 \times 585 \times 356$, where $N_x$, $N_y$ and $N_z$ denote the grid numbers in each direction. In the streamwise direction, the grid spacing varies linearly from $x=0$ to $x=75h_0$, and the spacing is set to be uniform from $x=75h_0$ to the end of the effective domain at $x=150h_0$. The stretching function from Passiatore et al. (2021) is employed in the wall-normal direction to refine the near-wall region, expressed as

(2.11)\begin{equation} \frac{y(j)}{L_y}=(1-\alpha) {\left ( \frac{j-1}{N_y-1} \right )}^3+\alpha \frac{j-1}{N_y-1}, \end{equation}

where $\alpha =0.25$ and $j \in [1,N_y]$. The grid is uniform in the spanwise direction.

The wall jet experiences an instability transition due to the prescribed inlet profile before entering a turbulent regime. The statistics are calculated through temporal and spanwise averages, except for the spectral analysis carried out in the spanwise direction in § 5, where only the temporal average is taken. The temporal average considers 6000 snapshots, nearly equivalent to the time duration of the flow passing through $750h_0$ in the streamwise direction. The temporal and spanwise average of a generic flow variable $\alpha$ is denoted as $\bar {\alpha }$, and the associated fluctuation is expressed as $\alpha '$. The Favre average is denoted as $\tilde {\alpha }$, and the associated fluctuation is expressed as $\alpha ''$.

Calculations are conducted with identical flow parameters and grids in both reacting and supplemented non-reacting settings. In the latter, the chemical reaction is artificially frozen, and only the transport equations of $\mathrm {H}_{2}$, $\mathrm {O}_{2}$, $\mathrm {H}_2{\mathrm {O}}$ and $\mathrm {N}_{2}$ are implemented. Table 2 presents the resulting boundary-layer properties and grid sizes for the reacting case at three downstream locations, whereas table 3 provides the corresponding results for the non-reacting case. The superscript ${\cdot }^*$ denotes the normalisation with the jet height $h_0$, whereas the superscript ${\cdot }^+$ denotes the normalisation with the viscous length scale $\delta _\nu =\bar {\mu }_w/(\bar {\rho }_w u_\tau )$, in which the subscript ${\cdot }_w$ denotes the variables at the wall. The friction velocity is denoted as $u_\tau =\sqrt {\bar {\tau }_w / \bar {\rho }_w}$ where $\bar {\tau }_w$ represents the averaged shear stress at the wall. The superscript ${\cdot }^\star$ denotes the semi-local scale, normalised with $\bar {\nu } / \sqrt {\bar {\tau }_w /\bar {\rho }}$. The Reynolds number based on the streamwise position $x$ is given by $Re_x=\rho _{\infty } U_{\infty } x/\mu _{\infty }$. The friction Reynolds number is defined as $Re_\tau =\delta / \delta _\nu$, where $\delta$ denotes the local boundary-layer thickness. The Reynolds number based on the boundary-layer thickness $\delta$ is denoted as $Re_\delta =\rho _{\infty } U_{\infty } \delta /\mu _{\infty }$. The Reynolds number based on the momentum thickness is defined as $Re_\theta =\rho _{\infty } U_{\infty } \theta /\mu _{\infty }$ with $\theta =\int _{0}^{\delta } ( \bar {\rho }\tilde {u}/\rho _0 u_0 )(1-\tilde {u}/u_0)\,{{\rm d}y}$. The friction Mach number is defined by $Ma_\tau =\bar {u}_\tau / c_w$ where $c_w$ is the gas sound velocity at the wall surface. The boundary-layer shape factor is defined as $H=\delta _d/\theta$ where $\delta _d$ denotes the local displacement thickness.

Table 2. Grid spacing and boundary-layer properties at three streamwise positions for the reacting case.

Table 3. Grid spacing and boundary-layer properties at three streamwise positions for the non-reacting case.

For the reacting case, the designed grid meets the recommended DNS requirements proposed by Poggie et al. (2015). In the supplemented non-reacting case, the first grid space at the wall $y^+_w$ is slightly higher than unity due to the higher values of $Re_\tau$. As the non-reacting case serves as a reference, we assume that this does not lead to qualitative differences in the comparison between the two cases. Additional discussions on the choice of flow parameter values and on the mesh adequacy are provided in Appendix B.

3. Flow visualisation

Spanwise sections of instantaneous flow fields are presented in figure 3. The prescribed inflow profile leads to large-scale vortices reminiscent of Kelvin–Helmholtz billows, a consequence of strong shear associated with velocity and density differences at the interface $y^*=1$. Further downstream, the mixing of preheated air and the hydrogen stream induces ignition in the reacting case, evidenced by the maximum temperature field reaching around 2200 K in figure 3(c). The velocity fields show that ignition results in a larger boundary-layer thickness, approximately 1.5 times that of the non-reacting case at $x^*=150$ (cf. tables 2 and 3). This observation aligns with theoretical analyses and numerical simulations of relevant supersonic combustors (Barth et al. 2013; Xue et al. 2021), where the increased boundary-layer thickness is attributed to the presence of flame surfaces in the boundary layer.

Figure 3. Instantaneous spanwise sections for the reacting (a,c,e,g,i,k) and non-reacting (b,d,f,h,j) cases.

In the reacting case, the mass fraction fields of hydrogen (figure 3e) and oxygen (figure 3g) reveal an excess of hydrogen in the boundary layer, with all oxygen consumed downstream. The variation of $Y_{\mathrm {O}_{2}}$ in the boundary layer is non-monotonic along the streamwise direction in the reacting case (figure 3g): the mixing process enhances the mass fraction of oxygen close to the wall, whereas ignition consumes oxygen from the streamwise position around $x^*=40$. Accordingly, the chemical reaction leads to the production of water vapour (figure 3i) and intermediate products such as $\mathrm {OH}$ (figure 3k). The distribution of the radical $\mathrm {OH}$ serves as an indicator of the flame surface position. It can be observed from figure 3(k) that the flame surface lies in the outer layer beneath the edge of the boundary layer, relatively far from the inner region where the boundary-layer shear is maximum. The flame surface becomes thicker as it moves downstream, consistent with higher temperature in the flame surface in figure 3(c). Conversely, in the non-reacting case, a pure mixing process of the cold hydrogen film with the preheated air occurs, with temperature and species concentrations displaying monotonic evolution in the boundary layer along the streamwise direction. The mass concentration field of $\mathrm {OH}$ for the non-reacting case is not shown in figure 3, as it is not produced.

3-D vortex structures are shown in figures 4 and 5 for the reacting and non-reacting cases, respectively. The identified vortex structures (green) are identified using the $\lambda _{c,i}$ criterion, an approach shown to be directly applicable to compressible wall-bounded turbulence (Zhou et al. 1999; Kolár 2009). In both reacting and non-reacting scenarios, 2-D instability modes manifest near the inlet at the stream interface $y^*=1$, associated with shear instability. These 2-D modes experience a transition into large-scale 3-D vortex structures around $x^*=10$, subsequently breaking into significantly smaller scales downstream at approximately $x^*=60$. In the reacting case, the vortex structures are lifted further away from the wall, possibly due to chemical effects such as thermal expansion and baroclinic torque generation. An abrupt breaking of lifted vortices is observed around $x^*=105$ (enlarged view in figure 4), after which the vortex structures are found in close proximity to the wall, with an apparent separation from the background high-temperature region.

Figure 4. 3-D visualisation for the reacting case. Green: vortex structures identified through $\lambda _{c,i}$-criterion. Background: temperature field. The entire domain from $x^*=0$ to $x^*=150$, is displayed alongside an enlarged view of vortex structures from $x^*=95$ to $x^*=130$.

Figure 5. 3-D visualisation for the non-reacting case. Green: vortex structures identified through $\lambda _{c,i}$-criterion. Background: temperature field. The entire domain from $x^*=0$ to $x^*=150$, is displayed alongside an enlarged view of vortex structures from $x^*=95$ to $x^*=130$.

4. Flow statistics

4.1. Transition process: interaction of mixing-layer and boundary-layer instability

We begin by characterising the streamwise evolution of streamwise kinetic energy $K_x=\overline {\rho u'' u''}$, and the heat release rate $\bar {\dot {\omega }}_T$. The streamwise kinetic energy integrated in the wall-normal direction in both reacting and non-reacting cases is presented in figure 6. Here $K_x$ experiences an immediate increase from the inlet, reaching a peak at $x^*=7.6$ before a rapid decay from the peak value. The decay slows down around $x^*=15$, followed by a slight increase in kinetic energy in both the reacting and non-reacting cases. The chemical heat release rate $\bar {\dot {\omega }}_T$ integrated in the wall-normal direction is represented by the red line in figure 6. We identify a streamwise position $x^*=21.7$, which corresponds to 1L' of the maximum heat release rate in the streamwise direction. It can be inferred that upstream of $x^*=21.7$, there is no combustion effect on the flow dynamics and the streamwise statistics of the reacting and non-reacting cases should closely match except for numerical and averaging errors. Further downstream, ignition occurs and modulates the flow dynamics, resulting in distinct evolutions of $K_x$ between the reacting and non-reacting cases.

Figure 6. Streamwise evolution of wall hydrogen mass fraction $\tilde {Y}_{\mathrm {H}_{2},w}$, skin friction $C_f$, kinetic energy $K_x$ and heat release rate $\bar {\dot {\omega }}_T$. The kinetic energy and the heat release rate are integrated along the wall-normal direction. Here $\tilde {Y}_{\mathrm {H}_{2},w}$, $C_f$ and $K_x$ are shown both for the reacting and non-reacting cases. The maximum value of each curve is normalised to unity.

Figure 7 presents the profiles of wall-normal velocity $\tilde {u}$ and streamwise kinetic energy $K_x$ in the range $0< x^*<21.7$ before the ignition tunes in. A budget analysis of turbulent production associated with $K_x$ is provided (Fang et al. 2020; Yu et al. 2022), which can be expressed as

(4.1)\begin{equation} P_{K_x}=\underbrace{-\overline{\rho u'' u''} \frac{\partial \tilde{u}}{\partial x}}_{{P_{11}^x}}\underbrace{-\overline{\rho u'' v''}\frac{\partial \tilde{u}}{\partial y}}_{P_{11}^y }, \end{equation}

where $P_{K_x}$ represents the total production of $K_x$, $P_{11}^x$ represents the production by streamwise evolution of mean flow and $P_{11}^y$ is the production by shear stress.

Figure 7. (a) Mean velocity profiles. (b) Mean kinetic energy profiles. (c–f) Turbulent production budgets $P_{K_x}=P_{11}^x+P_{11}^y$ at four streamwise positions. The same legend is used in (a,b) and the same line style is employed in (c–f).

The mean velocity profile experiences a rapid distortion between $x^*=1.9$ and $x^*=7.6$, where the double-shear mixing layer between the hydrogen stream and main stream is significantly attenuated. At $x^*=1.9$ very close to the inlet, two peaks appear in the kinetic energy profile, corresponding to the inner and outer shear layers on either side of the mixing layer, respectively. By examining the sign of $P_{11}^x$ in figure 7(c), we identify that the inner shear layer accelerates, whereas the outer shear layer decelerates. At $x^*=7.6$, where the kinetic energy reaches its maximum, the kinetic energy profile exhibits a larger peak on the outer side than on the inner side, with the latter being shifted towards the wall surface.

Since $x^*=7.6$, $K_x$ rapidly decays and attains a local maximum in the streamwise direction. At $x^*=15.0$, when compared with the kinetic profile at $x^*=7.6$, both inner and outer peaks in the energy profile are no longer discernible, because the distorted mean velocity profile has recovered to a monotonic one without shear regions far from the wall. This change is responsible for the rapid decay of $K_x$. The peak value of $P_{K_x}$ is located in the near-wall region dominated by $P_{11}^y$, indicating the onset of a boundary-layer instability. This is further supported by examining the streamwise evolution of skin friction $C_f$ and wall hydrogen mass fraction $\tilde {Y}_{\mathrm {H}_{2},w}$ in figure 6. In the range $x^*<7.6$, despite the strong increase in $K_x$, there is almost no change of $C_f$ and $\tilde {Y}_{\mathrm {H}_{2},w}$, suggesting that the near-wall region remains in a laminar regime. From $x^*>7.6$, $\tilde {Y}_{\mathrm {H}_{2},w}$ suddenly decreases, which indicates a significant portion of $\mathrm {H}_{2}$ is lifted by the strong flow motion near the wall. Concurrently, the skin friction increases, demonstrating the induced transition in the boundary layer. An outer layer peak is reestablished at $x^*=21.7$, as the near-wall fluctuations develops to the outer region.

In summary, the transition process associated with the prescribed wall–jet profile is characterised by two successive stages, each with distinct instability mechanisms. The first range $0< x^*<7.6$ corresponds to the mixing-layer-induced instability, marked by a sudden increase of kinetic energy without changes in the wall surface properties, and by a rapid distortion of mean velocity profiles. A second transition, occurring at $x^*>7.6$, is primarily governed by the boundary-layer transition, which is characterised by the increase in $C_f$. By the streamwise position $x^*=21.7$, $K_x$ has already evolved into a dual-peak profile, before the combustion tunes in.

4.2. Recovery process: outer-layer modulation by local density and viscosity

We examine the flow evolution from $x^*=21.7$ to the downstream end $x^*=150$, where combustion is active. The mean velocity profiles under the viscous scale, using the van Driest transformation, are presented in figure 8(a). A canonical adiabatic profile from Pirozzoli & Bernardini (2011) at $Ma=2$ and $Re_\tau =1000$, flow condition close to the presented profiles at $x^*=125$ and $x^*=150$ (cf. table 2), is provided for comparison. The profiles at $x^*=21.7$ and $x^*=50$ deviate from the canonical one, and the logarithmic law ($u^+=({1}/{\kappa } )\ln (y^+)+B$ where $\kappa =0.41$ and $B=5.2$). As the flow moves downstream, both profiles at $x^*=125$ and $x^*=150$ show good agreement with the canonical one until a wall-normal distance of $y/\delta \approx 0.3$. In the outer layer, the obtained profiles exhibit lower values than the canonical one. The effect of compressible transformation is further examined using that proposed in Griffin et al. (2021), where the transformed velocity, denoted as $u^+_{GFM}$, is represented under the semi-local scaling $y^\star$ in figure 8(b). A good match with the linear and logarithmic laws is also achieved. Note that a relatively good match in the inner layer can still be achieved without any compressible transformation due to the slow variation of near-wall density, as shown in figure 29(a) in Appendix A. This may be due to both the adiabatic wall condition (Modesti & Pirozzoli 2016) and also the prominent presence of hydrogen near the wall.

Figure 8. Mean velocity profiles using (a) the van Driest transformation and (b) the compressible transformation from Griffin, Fu & Moin (2021), in the reacting case at four streamwise positions. The linear and logarithmic laws are superposed, together with the canonical adiabatic profile at $Ma=2$ and $Re_\tau =1000$ from Pirozzoli & Bernardini (2011) in (a).

The evolution of Reynolds stress distributions is presented in figure 9 under the viscous scale $u^+$, with values normalised by the associated wall stresses $\bar {\rho }_w u_\tau ^2$. At the upstream position $x^*=21.7$, the profiles display significantly higher magnitudes than the canonical profiles, particularly in the outer region. This finding is likely to be attributed to the two-stage transition associated with the prescribed wall–jet inlet: a substantial amount of turbulent kinetic energy has been generated in the outer layer during the first transition stage linked to the mixing-layer instability, with additional energy supplemented from the near-wall turbulent fluctuations in the second stage. As the flow progresses downstream, the magnitudes of Reynolds stresses decrease. At streamwise positions $x^*=125$ ($Re_\tau =1067$) and $x^*=150$ ($Re_\tau =1051$), the Reynolds stresses recover to values that are visually lower than those of the canonical profile ($Re_\tau =1000$), with significant differences in the outer region.

Figure 9. Reynolds stresses normalised by the associated wall stresses at four streamwise positions in the reacting case. The canonical adiabatic profiles at $Ma=2$ and $Re_\tau =1000$ from Pirozzoli & Bernardini (2011) are superposed.

To investigate whether this difference is specifically linked to the combustion effect, the peak of normalised Reynolds shear stress $\bar {\tau }^+_{t,xy}=-\overline {\rho u'' v''}/(\bar {\rho }_w u_\tau ^2)$ in the wall-normal direction, denoted as $\bar {\tau }_{t,xy}^{pk}=\max _{y}(\bar {\tau }^+_{t,xy})$, is tracked along the streamwise for both the reacting and non-reacting cases (Ceci et al. 2022). Figure 10(a) shows that $\bar {\tau }_{t,xy}^{pk}$ decays in both cases along slightly different trajectories before reaching plateau regimes around $x^*\approx 70$. A correlation of $\bar {\tau }_{t,xy}^{pk}$ with $Re_\tau$ is provided in Ceci et al. (2022) for supersonic adiabatic boundary layers, taking the form $\tau _{t,xy}^{pk}=1-B_1 Re_{\tau }^{-6/7}$ with $B_1=13.62$. The nearly constant values achieved for both the reacting and non-reacting cases are found to be below the reference curve, suggesting that the primary driving factor for the negative difference in peak values is not limited to the combustion effect.

Figure 10. (a) Streamwise evolution of the peak value of Reynolds shear stress, $\bar {\tau }_{t,xy}^{pk}$, for both the reacting and non-reacting cases. A reference correlation from Ceci et al. (2022) is superposed. (b) Wall-normal profiles of local friction Reynolds number prefactor, $\sqrt {\bar {\rho }/\bar {\rho }_w} (\bar {\mu }_w/\bar {\mu })$, at different streamwise positions, complemented by a canonical adiabatic model profile based on Pirozzoli & Bernardini (2011).

We attribute this difference to the modulation effect of mean density and viscosity, associated with the inlet set-up of cold hydrogen jet stream and preheated air main stream. Typically, canonical wall-normal profiles are labelled by a prescribed friction Reynolds number $Re_\tau$ at the wall surface, along with a specific thermal wall condition, such as an adiabatic wall or a cold wall. The underlying assumption is that the wall-normal distribution of the associated density and viscosity profiles, which affect the local Reynolds number at a certain wall distance, should closely resemble those of the canonical one. If this assumption does not hold, simply using $Re_\tau$ is insufficient to accurately match a flow profile at the same $Re_\tau$ but with an arbitrary density and viscosity distribution. To characterise the local modulation of density and viscosity, we introduce the local friction Reynolds number (Patel, Boersma & Pecnik 2016), which is defined as

(4.2)\begin{equation} Re_{\tau}^{{\star}}=\sqrt{\frac{\bar{\rho}}{\bar{\rho}_w}}\frac{\bar{\mu}_w}{\bar{\mu}}Re_{\tau}. \end{equation}

Compared with $Re_{\tau }$, $Re_{\tau }^{\star }$ can represent the effect of local density and viscosity at a specific wall distance, while keeping the same characteristic velocity and length scale. Figure 10(b) presents the profile of the modulation prefactor $\sqrt {\bar {\rho }/\bar {\rho }_w}\bar {\mu }_w/\bar {\mu }$ in both the reacting and non-reacting cases. For comparison, a canonical adiabatic model is constructed using the density profile from Pirozzoli & Bernardini (2011) at $Re_\tau =1000$ and $Ma=2$, assuming that temperature is proportional to the inverse of density and that dynamic viscosity follows a $\frac {2}{3}$ power law against temperature. As shown by the red curve, the modulation prefactor $\sqrt {\bar {\rho }/\bar {\rho }_w}\bar {\mu }_w/\bar {\mu }$ increases with wall distance due to the decreasing temperature from the wall. In both the reacting and non-reacting cases, the profiles of this prefactor are found to be below the canonical curve, particularly in the outer layer. Since the reduced $Re_{\tau }^{\star }$ leads to locally reduced turbulent fluctuations, the magnitudes of Reynolds stresses in the present setting should be lower than those of the canonical adiabatic one at the same prescribed $Re_{\tau }$, which aligns with the observation at $x^*=125$ and $x^*=150$ in figure 9. A simple explanation can be given based on the temperature profiles in figure 10(c,d): with the presence of a cold hydrogen stream, the temperature does not evolve monotonically in the wall-normal direction, significantly deviating from a canonical adiabatic temperature profile where the temperature monotonically decreases from the wall surface.

The local density and viscosity depend on the local temperature and species concentrations, which can be influenced significantly by the combustion process. We examine the effect of combustion on the recovery process, and interpret the distinct recovery trajectories concerning $\bar {\tau }^{pk}_{t,xy}$ in figure 10(a). At upstream positions, $\bar {\tau }^{pk}_{t,xy}$ decays more slowly in the reacting case and exhibits larger peak magnitudes, whereas downstream, it reaches a plateau of smaller magnitudes, as shown by the zoomed plot in figure 10(a). The wall-normal profiles of density, viscosity and $Re_{\tau }^{\star }$ at $x^*=50$ and $x^*=150$ are presented in figure 11(a,b), along with those of temperature and species in figure 11(c,d). At $x^*=50$, we observe larger $Re_{\tau }^{\star }$ in the reacting case than in the non-reacting case, consistent with the larger $\bar {\tau }^{pk}_{t,xy}$ at the same position in figure 10(a). The values of $Re_\tau$ are close at the wall, and $\sqrt {\bar {\rho }/\bar {\rho }_w}$ exhibits a relatively more pronounced difference between the reacting and non-reacting cases, relative to that of $\bar {\mu }_w/\bar {\mu }$. Moreover, $\sqrt {\bar {\rho }/\bar {\rho }_w}$ scales with $1-\bar {Y}_{\mathrm {H}_{2}}$ at different wall-normal distances. Therefore, the slower decay of Reynolds shear stresses at $x^*=50$ is primarily attributed to the modulation of hydrogen distribution by combustion. Effectively, considering that hydrogen is very light, the depletion of hydrogen by combustion in the flow results in an increase in density, leading to larger values of $Re_{\tau }^{\star }$.

Figure 11. (a,b) Mean profiles of local friction Reynolds number $Re_{\tau }^{\star }$, $\sqrt {\bar {\rho }/\bar {\rho }_w}$ and $\bar {\mu }_w/\bar {\mu }$. (c,d) Mean profiles of temperature $\bar {T}$ and mass fraction of all the species except hydrogen $1-\bar {Y}_{\mathrm {H}_{2}}$. Both reacting and non-reacting cases are presented, at two streamwise positions $x^*=50$ (a,c) and $x^*=150$ (b,d).

At the downstream position $x^*=150$, the chemical heat release significantly raises the wall temperature, leading to lower $Re_\tau$ at the wall and subsequently lower $Re_{\tau }^{\star }$ in the flow. An attenuated peak value of Reynolds shear stress is consistently found in comparison to the non-reacting case.

In summary, the influence of combustion on the recovery process is attributed to its modulation of mean density and viscosity profiles through the chemical heat release and the depletion of hydrogen. The present finding aligns with those in reacting mixing layers, where the reduction of turbulent kinetic energy and Reynolds stresses around the flame surface was attributed to the mean density effect modulated by chemical heat release, compared with the inert case (Mahle et al. 2007).

4.3. Effect of combustion on the skin friction

The effect of combustion on the wall surface is assessed by comparing the skin friction $C_f$, the wall temperature $T_w$, the wall dynamic viscosity $\mu _w$ and the wall shear strain $({\rm d}U/{{\rm d}y})_w$ in the reacting and non-reacting cases. The results are presented in figures 12(a)–12(c). The skin friction undergoes nearly the same initial increasing stage, and from $x^*\approx 50$, the reacting case shows a more pronounced decreasing trend than the non-reacting one. At the downstream end $x^*=150$, the skin friction in the reacting case is approximately 80 % of that in the non-reacting case. The adiabatic condition $\partial T / \partial y = 0$ allows for variations in wall surface temperature. The chemical heat release results in a significantly higher wall temperature with respect to the non-reacting case, subsequently leading to higher viscosity at the wall $\mu _w$. As $C_f$ is proportional to the product of $\mu _w$ and $({\rm d}U/{{\rm d}y})_w$, the lower value of $C_f$ with combustion is because of a reduced wall shear stress. The relationship between $C_f$ and $Re_\theta$ is shown in figure 12(d) using samples from the downstream positions $125< x^*<150$. A correlation $C_f=A_3 Re_\theta ^{-B_3}\times 10^{-3}$ for canonical adiabatic boundary layer at $Ma=2$ is added for comparison, where $A_3=17.4$ and $B_3=0.245$ (Ceci et al. 2022). It is not surprising to find that $C_f$ in both the reacting and non-reacting cases is lower than the canonical one, as the wall-normal flow profile significantly differs from a canonical adiabatic profile with smaller local Reynolds numbers.

Figure 12. (a) Skin friction $C_f$. (b) Wall temperature $T_w$. (c) Wall shear strain $({\rm d}U/{{\rm d}y})_w$ and viscosity $\mu _w$. (d) Relation between $C_f$ and $Re_\theta$ in the streamwise range $125< x^*<150$. A correlation of the canonical adiabatic supersonic boundary layer from Ceci et al. (2022) is superposed.

The combustion-related modulation on the skin friction is quantified by the Renard–Deck decomposition extended to a compressible boundary-layer setting (Renard & Deck 2016; Fan, Li & Pirozzoli 2019), which is expressed as

(4.3)\begin{align} C_f&= \underbrace{\frac{2}{\rho_{\infty}u_{\infty}^3}\int_0^\infty\bar{\tau}_{v,xy}\frac{\partial \tilde{u}}{\partial y} {{\rm d}y}}_{C_{f,V}} +\underbrace{\frac{2}{\rho_{\infty}u_{\infty}^3}\int_0^\infty\bar{\tau}_{t,xy}\frac{\partial \tilde{u}}{\partial y} {{\rm d}y}}_{C_{f,T}}\nonumber\\ &\quad + \underbrace{\frac{2}{\rho_{\infty}u_{\infty}^3}\int_0^\infty (\tilde{u}-u_\infty)\left ( \bar{\rho}\left ( \tilde{u}\frac{\partial \tilde{u}}{\partial x}+\tilde{v}\frac{\partial \tilde{u}}{\partial y} -\frac{\partial}{\partial x}(\bar{\tau}_{v,xx}+\bar{\tau}_{t,xx}) \right )\right )}_{C_{f,G}}. \end{align}

The terms $C_{f,V}$, $C_{f,T}$ and $C_{f,G}$ represent the wall-normal contributions from viscous shear stress, Reynolds shear stress and streamwise variation, respectively. The results of this decomposition in both the reacting and non-reacting flows are presented in figure 13(a). The contribution from streamwise variation $C_{f,G}$ exhibits similar trends in both cases, accounting for 18 % of total $C_f$ at $x^*=150$ in the reacting case. The viscous and Reynolds shear stresses contribute equally over a large range of streamwise positions, each of which accounts for 41 % of total skin friction at $x^*=150$. Both $C_{f,V}$ and $C_{f,T}$ show smaller values than the non-reacting case, with the Reynolds shear stress contribution $C_{f,T}$ exhibiting a larger difference. These results align with the findings that shear strains and Reynolds shear stresses in the flow are reduced due to combustion.

Figure 13. (a) Wall skin friction using Renard–Deck decomposition. Here $C_{f,V}$, $C_{f,T}$ and $C_{f,G}$ denote the budgets from viscous shear stress, Reynolds shear stress and streamwise variation, respectively. (b) Splitting Reynolds shear stress’ budget $C_{f,T}$ into a canonical portion $C_{f,T_c}$ and a non-canonical portion $C_{f,T_n}$. Black, reacting; blue, non-reacting.

To further explain the considerable difference of $C_f$ between the present configuration and the canonical adiabatic correlation in figure 12(d), following Yu et al. (2023a), we decompose the wall-stress-normalised Reynolds shear stress $\bar {\tau }^+_{t,xy}=\bar {\tau }_{t,xy}/(\bar {\rho }u^2_\tau )$ into a canonical portion $\bar {\tau }^+_{c,xy}=\bar {\tau }_{c,xy}/(\bar {\rho }u^2_\tau )$ and a non-canonical portion $\bar {\tau }^+_{n,xy}$, as expressed by

(4.4)\begin{equation} \bar{\tau}^+_{t,xy}=\bar{\tau}^+_{c,xy}+\bar{\tau}^+_{n,xy}. \end{equation}

Given the relatively small variation of $Re_\tau$ in the streamwise direction (cf. tables 2, 3), we assume a constant canonical Reynolds shear stress $\bar {\tau }^+_{c,xy}$ at each streamwise position (Yu et al. 2023a). The canonical adiabatic profile shown in figure 9(d) is employed, transformed to $\bar {\tau }_{c,xy}$ and $y/\delta$ using local wall quantities. The resulting $\bar {\tau }_{t,xy}$, $\bar {\tau }_{c,xy}$ and $\bar {\tau }_{n,xy}$ are presented in figure 14. At $x^*=21.7$, the non-canonical portion is characterised by large positive magnitudes due to the specific transition process, which was addressed as the mixing-layer-induced portion in a shock-impingement scenario (Yu et al. 2023a). As the mixing-layer-induced portion gradually diminishes downstream, the non-canonical portion becomes negative due to the modulation of mean velocity and density profiles in the outer layer. The skin-friction budget from Reynolds shear stress can be further decomposed into a canonical portion $C_{f,Tc}$ and a non-canonical portion $C_{f,Tn}$, expressed as

(4.5)\begin{equation} C_{f,T}=\underbrace{\frac{2}{\rho_{\infty}u_{\infty}^3}\int_0^\infty\bar{\tau}_{c,xy}\frac{\partial \tilde{u}}{\partial y} {{\rm d}y}}_{C_{f,Tc}}+\underbrace{\frac{2}{\rho_{\infty}u_{\infty}^3}\int_0^\infty\bar{\tau}_{n,xy}\frac{\partial \tilde{u}}{\partial y} {{\rm d}y}}_{C_{f,Tn}}. \end{equation}

The associated decomposition results are presented in figure 13(b). At $x^*=30$, the mixing-layer-induced portion accounts for approximately half of the $C_{f,Tc}$. This effect nearly disappears by $x^*=70$. Further downstream, the modulation effect results in a negative contribution to $C_{f,Tc}$, which is more important in the reacting case: $C_{f,Tn}$ accounts for $-24\,\%$ and $-5\,\%$ of $C_{f,T}$, respectively, at the downstream end $x^*=150$.

Figure 14. Splitting the Reynolds shear stress (a) into a canonical portion (b) and a non-canonical portion (c). Wall-normal distributions under $y/\delta$ at four streamwise positions in the reacting case are presented. All the shear stresses are normalised by the free-stream variables $\rho _{\infty } U^2_{\infty }$.

5. Temperature and species statistics

5.1. Mean profiles, fluctuations and spanwise spectra

The profiles of averaged temperature and their corresponding fluctuation counterparts are provided in figure 15(a,b). Non-monotonic mean temperature profiles are observed in both the reacting and non-reacting cases. The elevated temperature near the wall aligns with the adiabatic condition, where the viscous heat generated within the intense shear region of the turbulent boundary layer is constrained by the solid wall. As one moves away from the wall, the averaged temperature decreases due to the presence of the cold hydrogen stream. Further moving outward, in the reacting case, a rise in temperature occurs in the outer layer corresponding to the chemical heat release, followed by a decrease towards the main-stream temperature beyond the boundary layer. Conversely, in the non-reacting case, an increase in temperature in the outer layer is attributed to turbulent mixing with the preheated main stream. In both scenarios, the averaged temperature rises downstream, consistent with the observed increase in wall temperature in figure 12(b). At $x^*=150$, the temperature peak for the reacting case is identified at $y^+=795$, whereas the peak associated with maximal temperature fluctuation in figure 15(b) is found at $y^+=576$, approximately corresponding to the maximal gradient of the averaged profile. For the non-reacting case, the temperature fluctuation peaks are also correlated with the maximal growth of the average profile. Notably, the non-reacting case exhibits higher peak values in temperature fluctuations with respect to the reacting case due to the larger temperature gradient in the former. The mean profiles of density and pressure fluctuations, together with transport property and Mach-number statistics are provided in Appendix A.

Figure 15. Averaged profiles of static temperature (a) along with the root-mean-squares of their fluctuation counterparts in (b). The profiles are evaluated at three streamwise positions in both reacting (black) and non-reacting (blue) cases.

Figure 16 shows the mean mass fraction of each species. In the non-reacting scenario, only four species with non-zero values prescribed at the inlet are displayed. In the reacting cases, all nine species are presented, including intermediate products $\mathrm {O}$, $\mathrm {OH}$, $\mathrm {H}$ and $\mathrm {HO}_{2}$. The non-reacting case exhibits nearly uniform profiles of $\mathrm {O}_{2}$ and $\mathrm {H}_2{\mathrm {O}}$ due to pure turbulent mixing at the three downstream positions. A higher concentration of $\mathrm {H}_{2}$ is observed in proximity to the wall, with a diminished presence of $\mathrm {N}_2$. Chemical reactions induce a redistribution of species by consuming $\mathrm {O}_{2}$ at the wall, leading to the region being filled with produced $\mathrm {H}_2\mathrm {O}$. The peaks associated with the averaged intermediate product profiles primarily fall within the range $780< y^+<840$, corresponding to the maximum of the averaged temperature profile. An exception involves the profile of $\mathrm {H}$ with prominent distributions observed at the wall. The presence of $\mathrm {H}$ and $\mathrm {OH}$ at the wall is further investigated in § 6. The mean fluctuations of species mass fraction are provided in figure 31 in Appendix A.

Figure 16. Averaged profiles of species mass fraction $\tilde {Y}_i$: (a) O; (b) ${\rm O}_2$; (c) ${\rm H}_2$; (d) ${\rm H}_2{\rm O}$; (e) OH; ( f) H; (g) ${\rm HO}_2$; (h) ${\rm H}_2{\rm O}_2$; (i) ${\rm N}_2$. The profiles are evaluated at three streamwise positions in the reacting (black) and non-reacting cases (blue).

A spectral analysis is performed in the periodic spanwise direction. The spanwise Fourier transform associates a generic fluctuation $\alpha '(x,y,z,t)$ with its spectral counterpart $\hat {\alpha }(x,y,k_z,t)$, where $k_z$ denotes spanwise wavenumbers. The spanwise spectra $E_{\alpha \alpha }(x,y,k_z)=\overline {\frac {1}{2}\hat {\alpha }\hat {\alpha }^\star }$ at different wall-normal distance are computed for streamwise velocity, temperature and hydrogen mass fraction fluctuations at the streamwise location $x^*=150$. The obtained spectra together with the spanwise snapshots under linear scale $y^*$ are presented in figure 17. The premultiplied velocity spectra $k_z E_{uu}$ reveal dual peaks in figure 17(a,g). The first peak is very close to the wall, aligning with the peak of Reynolds normal stresses in figure 9(a). The second peak is located in the outer layer, corresponding to larger spanwise wavelengths. This bimodal structure is often observed in confined flows with high Reynolds numbers (Yu & Xu 2022). The related velocity snapshots (figure 17b,h) consistently show that velocity fluctuations are prominent in the inner region. Large-scale velocity structures visually exist in the outer region, but their associated fluctuation energy is considerably smaller.

Figure 17. Premultiplied spectra $k_z E_{uu}$, $k_z E_{TT}$ and $k_z E_{Y_{\mathrm {H}_{2}}Y_{\mathrm {H}_{2}}}$ for reacting (a–f) and non-reacting (g–l) cases at $x^*=150$. The linear scale $y^*$ is employed, and $\lambda ^*_z$ represents the wavelength corresponding to $k_z$ under linear scale. The normalisation of $E_{uu}$ and $E_{TT}$ ensures their maximum value is unity. The white dashed lines indicate the boundary-layer thickness. Instantaneous streamwise sections of $u$, $T$ and $Y_{\mathrm {H}_{2}}$ are provided using $y^*$ and $z^*$ as the wall-normal scale accordingly.

The temperature spectra $k_z E_{TT}$ and mass fraction of hydrogen spectra $k_z E_{Y_{\mathrm {H}_{2}}Y_{\mathrm {H}_{2}}}$, exhibit particular features due to the wall–jet configuration. In the reacting case (figure 17c), the temperature peak is situated at the wall-normal position approximately halfway through the boundary-layer thickness. This peak also corresponds to a wavelength approximately half that of the boundary-layer thickness. In the associated snapshot (figure 17d), chemical reaction induces local heat release, causing spatial inhomogeneity throughout the wall-normal direction within the boundary layer. The spectra of hydrogen mass fraction $k_z E_{Y_{\mathrm {H}_{2}}Y_{\mathrm {H}_{2}}}$ (figure 17e) share close resemblance with the temperature spectra. The cold flow region beneath $y^*=5$ and in the spanwise region $z^*<2$ and $z^*>8$ (periodic) in (figure 17d) correlates with excess hydrogen (figure 17f). In the outer region above $y^*=5$, the chemically heated flow is mixed with the preheated main stream, resulting in temperature fluctuations at the boundary layer's edge. In both reacting and non-reacting cases, the temperature spectra resemble the hydrogen spectra, indicative of the cold flow region strongly linked to hydrogen. A closer examination shows that there also exists fluctuations in the temperature spectra (figure 17c,i) near the wall, aligning with the peak in the velocity spectra. However, the near-wall temperature fluctuations are considerably smaller than those observed in the outer layer and visually disappear upon normalisation.

5.2. Turbulent Prandtl numbers

We define the correlation coefficients between two arbitrary variables $\alpha$ and $\beta$ as

(5.1)\begin{equation} R_{\alpha,\beta}=\frac{\overline{\alpha \beta}}{\sqrt{\overline{\alpha^2}}\sqrt{\overline{\beta^2}}}. \end{equation}

The correlation coefficients between velocity and temperature fluctuations are examined in figures 18(a)–18(d). Previous studies have indicated that the velocity–temperature correlation values depend on the gradient of averaged temperature profiles (Duan et al. 2010; Di Renzo & Urzay 2021). Positive $R_{u'',T''}$ is mostly associated with positive gradient of averaged temperature, which is attributed to ejection and sweep events in the turbulent boundary layers, and vice versa. This observation is validated in the current calculation, as shown in (c,d).

Figure 18. (a,b) Correlation coefficients between the fluctuations of wall-normal velocity and temperature $R_{v'',T''}$. (c,d) Correlation coefficients between the fluctuations of streamwise velocity and temperature $R_{u'',T''}$. (e,f) Turbulent Prandtl numbers $Pr_{t}$. (a,c,e) Reacting case (black). (b,d,e) Non-reacting case (blue). The axes $R_{v'',T''}=0$ and $R_{u'',T''}=0$ are superposed in (a–d) and $Pr_{t}=0.9$ is superposed in (e,f).

In the reacting case, $R_{u'',T''}$ exhibits two zero points. The first corresponds to the local temperature minimum linked to the cold hydrogen stream, whereas the second corresponds to the local temperature maximum related to the flame surface. Conversely, in the non-reacting case, $R_{u'',T''}$ has only one zero point inside the boundary layer, corresponding to the local temperature minimum associated with the cold hydrogen stream. The wall-normal velocity–temperature correlation $R_{v'',T''}$ shows an opposite sign compared with $R_{u'',T''}$.

The turbulent Prandtl number $Pr_t$, associating the Reynolds shear stress $\overline {\rho u'' v''}$ (figure 9d) with the turbulent heat transport $\overline {\rho v'' T''}$, is defined as

(5.2)\begin{equation} Pr_t=\frac{\overline{\rho u'' v''} \partial \tilde{T} / \partial y}{\overline{\rho v'' T''} \partial \tilde{u} / \partial y}. \end{equation}

As $R_{v'',T''}$ involves in the denominator of $Pr_t$ with a closely similar form, consequently, apart from the grid points very close to the wall surface, the wall-normal profiles of turbulent Prandtl number exhibit two discontinuities corresponding to the two zero points of $R_{v'',T''}$ in the reacting case, and one discontinuity in the non-reacting case. Such discontinuous turbulent Prandtl number distributions associated with non-monotonic temperature profiles were recently observed in hypersonic high-enthalpy boundary-layer flows (Di Renzo & Urzay 2021) and supersonic boundary layers with film cooling (Peter & Kloker 2022). We argue that such a discontinuous profile may not raise specific concerns in the closure of turbulent heat transport terms using $Pr_t$: as both $\partial \tilde {T} / \partial \tilde {y}$ and $\overline {\rho v'' T''}$ are close to zero around the discontinuous location, using a finite value of $Pr_t$ is sufficient to obtain a turbulent heat transport term $\overline {\rho v'' T''}$ close to zero.

The obtained turbulent Prandtl numbers are presented under inner scale $y^+$ in figure 18(e,f). For both the reacting and non-reacting cases, $Pr_t$ is found to be around 1.1 in the logarithmic region and varies from 0.5 to 0.9 in the outer layer, as shown in figure 18(e,f). The observed trends are consistent with recent DNS results of film cooling (figure 6 in Peter & Kloker 2022) where $Pr_t\approx 1$ was found close to the wall followed by a gradual decrease to $Pr_t\approx 0.7$ in the outer layer. Although the chemical reaction introduces an extra discontinuity close to the boundary edge, it does not lead to dramatic change of turbulent Prandtl numbers in the inner region and in a large portion of the outer layer region, with respect to the non-reacting case.

5.3. Turbulent Schmidt numbers

Turbulent Schmidt numbers, commonly employed to model species transport fluxes $\overline {\rho v'' Y''_i}$, are evaluated in both reacting and non-reacting scenarios at three streamwise positions. For each species $i$, it is defined as

(5.3)\begin{equation} Sc_{t,i}=\frac{\overline{\rho u'' v''} \partial \tilde{Y}_i / \partial y}{\overline{\rho v'' Y''_i} \partial \tilde{u} / \partial y}. \end{equation}

The wall-normal profiles of $Sc_{t,i}$ are presented under inner scale $y^+$ in figure 19. In the non-reacting case, turbulent Schmidt numbers for all four species exhibit an increasing trend from the wall to the edge of the boundary layer, with values around 0.8 in a large portion of the outer layer. This increasing trend and the obtained values align with the DNS film cooling results (panels $x=30$ and $x=60$ in figure 26 of Peter & Kloker 2022).

Figure 19. Turbulent Schmidt numbers $Sc_{t,i}$ with inner scale $y^+$: (a) O; (b) ${\rm O}_2$; (c) ${\rm H}_2$; (d) ${\rm H}_2{\rm O}$; (e) OH; ( f) H; (g) ${\rm HO}_2$; (h) ${\rm H}_2{\rm O}_2$; (i) ${\rm N}_2$. The profiles are evaluated at three streamwise positions in the reacting (black) and non-reacting (blue) cases.

Chemical reactions slightly alter the values of turbulent Schmidt numbers. The values associated with $\mathrm {O_2}$, $\mathrm {H}_{2}$, $\mathrm {H_2O}$ and $\mathrm {N_2}$ are found to be larger than unity in the outer layer close to the edge. Here $Sc_{t,\mathrm {H_2O}}\approx 1$ is identified in the inner layer where the mass fraction of $\mathrm {H}_2\mathrm {O}$ is prominent. The intermediate products $\mathrm {O}$, $\mathrm {OH}$, $\mathrm {HO}_{2}$ and $\mathrm {H}_2\mathrm {O}_{2}$ exhibit values around 0.8 in the outer regions linked to the peaks in their mass fraction profiles. The zero points in $\overline {\rho v'' Y''_i}$ introduce discontinuities in the turbulent Schmidt numbers, leading to distorted profiles for those intermediate products at around $y/\delta =0.8$. $Sc_{t,\mathrm {H}}$ is identified in the range from 0.8 to 1.5 in the outer region.

5.4. Strong Reynolds analogy

In this section, we explore the validity of the SRA and its subsequent extensions, which relate the temperature and velocity fluctuations. The original version of SRA proposed by Morkovin (1962) is written as

(5.4)\begin{equation} \frac{T''_{rms}/\tilde{T}}{(\gamma-1) Ma^2 (u''_{rms}/\tilde{u})} =1. \end{equation}

The formula is derived for an ideal gas with a constant $C_p$ based on the assumption that the fluctuation of total temperature is zero. In the present multi-species setting, $C_p$ depends on the temperature and species concentrations. If we maintain $C_p$ in the formula and assume $\tilde {C}_p T''_{rms} \gg C''_p \tilde {T}$, the expression can be written as

(5.5)\begin{equation} \frac{T''_{rms} \tilde{C}_p}{u''_{rms} \tilde{u}}=1. \end{equation}

The evaluation of the left-hand side of (5.4) and (5.5) is presented in figure 20. Both the original SRA using constant $C_p$, and the adapted one using a variable $C_p$, exhibits a nearly constant value within the range $10< y^+<200$. The adapted one shows a closer value to unity. In the outer region, both of these expressions significantly deviate from the unity.

Figure 20. SRA and its extensions at different streamwise positions for the reacting (a,c,e) and non-reacting (b,d,f) cases. SRA and SRA-$C_p$ represent the left-hand side of (5.4) and (5.5), respectively.

Other versions, including ESRA (Cebeci 2012), GSRA (Gaviglio 1987), RSRA (Rubesin 1990) and HSRA (Huang, Coleman & Bradshaw 1995), are also assessed. Readers may refer to Duan et al. (2010) and Yu et al. (2023b) for the detailed expressions. Figure 20 depicts the ratio between the left-hand-side and right-hand-side values for these various formulations. The GSRA, RSRA and HSRA exhibit severe deviations from unity, whereas the values for ESRA fall outside the displayed range. The validity of underlying assumptions associated with those extended versions is further examined. In the ESRA, a similarity relation between the total temperature $T_t$ and the streamwise velocity is assumed, following $(\tilde {T}_t-\tilde {T}_w)/(\tilde {T}_{t,\infty }-\tilde {T}_w)=\tilde {u}/u_\infty$. However, the temperature profile is not monotonic due to the wall–jet inlet setting and this similarity is lost. The associated ratio $(\tilde {T}_t-\tilde {T}_w)/(\tilde {T}_{t,\infty }-\tilde {T}_w)/(\tilde {u}/u_\infty )$ is examined in figure 21(a), which is significantly different from unity. In the GSRA, RSRA and HSRA, the mixing length assumption is employed for both velocity and temperature fluctuations, leading to the proportional relation $u''(\partial \tilde {T} / \partial y) \propto T'' (\partial \tilde {u} / \partial y)$. The associated ratio $(u'' / T'') (\partial \tilde {T} / \partial \tilde {u} )$ is examined in figure 21(b), which also deviates from a constant. Since the underlying assumptions of those extended versions cannot hold, their results show no improvement with respect to the original SRA.

Figure 21. Examination of the assumptions in SRA extensions: (a) $(\tilde {T}_t-\tilde {T}_w)/(\tilde {T}_{t,\infty }-\tilde {T}_w)/(\tilde {u}/u_\infty )$ for ESRA; (b) $(u'' / T'') (\partial \tilde {T} / \partial \tilde {u} )$ for GSRA, RSRA and HSRA. The results of the reacting case are shown at three streamwise positions.

To comprehend the relative validity of the original SRA in the inner region and its invalidity in the outer region, we additionally calculate a linear coherence spectrum $\xi _{u,T}^2(x,y,k_z)$, defined as

(5.6)\begin{equation} \xi_{u,T}^2=\frac{\overline{\hat{T} \hat{u}^\blacktriangle}}{\overline{\hat{T} \hat{T}^\blacktriangle}\overline{\hat{u} \hat{u}^\blacktriangle}}, \end{equation}

to examine the spectral correlation between temperature and velocity fluctuations (Baars, Hutchins & Marusic 2016; Yu & Xu 2021). The complex conjugate is represented by ${\cdot }^\blacktriangle$ and $\bar {{\cdot }}$ exceptionally denotes only the average over time in the definition of spectra. The linear coherence spectra of the reacting case at $x^*=150$, as shown in figure 22, allows us to measure the temporal dependence of velocity and temperature fluctuations at various spanwise wavelengths $\lambda _z$. The coherence coefficient is found at around 0.5 for $y^+ < 200$ across a broad spectrum of wavelengths, with the maximum coherence value of 0.561, indicating a partial correlation between velocity and temperature fluctuations. Consistently, the left-hand side of the SRA identity attains nearly constant values at wall-normal positions $y^+ < 200$ in figure 20.

Figure 22. Linear coherence spectra $\xi _{u,T}^2(x,y)$ in the reacting case at $x^*=150$. The logarithmic scale $y^+$ is used, and $\lambda ^+_z$ represents the wavelength corresponding to $k_z$. The white dashed line indicates the boundary-layer thickness.

In the region $y^+ > 200$ (equivalent to $y^*>1.62$), the coherence between temperature and velocity is lost due to large-scale temperature fluctuations associated with turbulent mixing and local chemical heat release. The inhomogeneity in temperature distribution is inherently linked to the distribution of hydrogen gas. However, the chemically active large-scale scalar fluctuations associated with temperature and hydrogen do not induce, in turn, important large-scale velocity fluctuations. The spectral coherence is consequently lost in the outer layer.

6. Reaction statistics

6.1. Chemical reaction rates and analysis of elementary reactions

Figure 23 presents the mean reaction rates of each species $\bar {\dot {\omega }}_i$, where positive and negative values represent chemical production and consumption, respectively. These reaction rates are displayed without normalisation in figure 23, which reveals a gradual decrease in chemical reaction rates downstream. Since the reaction rate only depends on the local temperature and mixture composition while the temperature rises downstream, the reducing rate values are likely attributed to the depletion of reactant concentrations in the mixture. The peak rate values of species are predominantly situated in the outer layer, aligning with the temperature peak. Apart from the outer region, significant reaction rates are observed at the wall for radicals $\mathrm {OH}$ and $\mathrm {H}$, along with non-negligible near-wall rates in $\mathrm {H}_{2}$ and $\mathrm {H}_2\mathrm {O}$.

Figure 23. Mean reaction rates of each species $\bar {\dot {\omega }}_i$: (a) O; (b) ${\rm O}_2$; (c) ${\rm H}_2$; (d) ${\rm H}_2{\rm O}$; (e) OH; ( f) H; (g) ${\rm HO}_2$; (h) ${\rm H}_2{\rm O}_2$; (i) ${\rm N}_2$. The axes corresponding to zero reaction rate are superposed (grey dashed).

We explore the elementary reaction rates that govern the chemical production and consumption of species in the near-wall region. By substituting (2.8) into (2.9) and taking the Reynolds average, the mean reaction rate of the $i$th species is expressed as follows:

(6.1)\begin{equation} \bar{\dot{\omega}}_{i}= W_i \sum_{j=1}^{M} \left ( \chi'_{i,j}- \chi_{i,j} \right ) \bar{\phi}_{j}. \end{equation}

From this expression, the elementary reaction rates of the $j$th reaction step $\bar {\phi }_{j}$ contribute to the reaction rates of the $i$th species by multiplying the associated stoichiometric coefficients $\chi _{i,j}'- \chi _{i,j}$ and the molecular mass $W_i$. All 19 elementary reaction rates are examined, most of which show peak values in the outer layer corresponding to the maximum heat release rate, and near-zero values close to the wall. It is found that the elementary reactions 3, 5 and 8 show non-negligible values close to the wall, as presented in figure 24. The corresponding chemical reaction steps are detailed in table 4. The reaction rate of elementary reaction 3 is negative at the wall, indicating the consumption of $\mathrm {H}$ and $\mathrm {H_2O}$, and the production of $\mathrm {OH}$ and $\mathrm {H}_{2}$. Elementary reaction 5 results in the consumption of $\mathrm {H}$ and the production of $\mathrm {H}_{2}$. Meanwhile, elementary reaction 8 favours the forward reaction, leading to the consumption of $\mathrm {OH}$ and $\mathrm {H}$ with the production of $\mathrm {H_2O}$. A summary of those results provides clear insights into the near-wall chemistry. All three elementary reactions contribute to the consumption of $\mathrm {H}$. Reactions 3 and 5 significantly contribute to the production of $\mathrm {H}_{2}$. Although elementary reactions 3 and 8 promote the evolution of $\mathrm {H_2O}$ and $\mathrm {OH}$ in opposite directions, reaction 3 exhibits a rate magnitude twice that of reaction 5. With identical stoichiometric coefficients in both reactions, their combined effect results in near-wall consumption of $\mathrm {H_2O}$ and production of $\mathrm {OH}$, consistent with the observed negative and positive rates, respectively.

Figure 24. Mean elementary reaction rates $\bar {\phi }_{j}$ of reaction steps 3, 5 and 8: (a) $\bar {\phi }_{3}$; (b) $\bar {\phi }_{5}$; (c) $\bar {\phi }_{8}$.

Table 4. Elementary reactions 3, 5 and 8 from 9 species 19 steps hydrogen–air reaction model (Li et al. 2004). In both reactions 5 and 8, $\mathrm {M}$ represents a third body being an arbitrary species.

6.2. Chemical heat release rate budget

By taking the Reynolds average of (2.10), the mean chemical heat release rate can be expressed as

(6.2)\begin{equation} \bar{\dot{\omega}}_T={-}\sum_{i=1}^{N} \bar{\dot{\omega}}_{T,i}, \end{equation}

where $\bar {\dot {\omega }}_{T,i} = \Delta h_{f,i}^o \bar {\dot {\omega }}_i$ represents the heat release budget resulting from the mass variation of ith species due to chemical reaction. Figure 25 presents $\bar {\dot {\omega }}_{T,i}$ for each species. The magnitudes of the heat release budget associated with $\mathrm {H_2O}$ and $\mathrm {H}$ are at least an order of magnitude higher than those of the other species. In figure 26, the budget profiles of $\bar {\dot {\omega }}_{T,\mathrm {H_2O}}$, $\bar {\dot {\omega }}_{T,\mathrm {H}}$ and their sum are presented, alongside the total mean heat release rate $\overline {\dot {\omega }_T}$. At all three downstream positions investigated, the sum of $\bar {\dot {\omega }}_{T,\mathrm {H_2O}}$ and $\bar {\dot {\omega }}_{T,\mathrm {H}}$ closely matches $\bar {\dot {\omega }}_T$. The production of $\mathrm {H}_2\mathrm {O}$ dominates the peak of the heat release rate in the outer layer. The depletion of $\mathrm {H}_2\mathrm {O}$ in the near-wall region due to elementary reaction 3 results in slight heat absorption. However, the depletion of $\mathrm {H}$ prevails at the wall, leading to heat release due to the positive standard formation enthalpy of $\mathrm {H}$. We conclude that the chemical heat release is dominated by $\mathrm {H}_2\mathrm {O}$ in the outer layer and by $\mathrm {H}$ at the wall, resulting in prominent heat release rates across the entire boundary layer. Note that these results are obtained under adiabatic non-catalytic conditions for the wall, and the validity of the results concerning near-wall reactions may not hold with other wall conditions.

Figure 25. Budgets of heat release rate from each species $\bar {\dot {\omega }}_{T,i}$: (a) O; (b) ${\rm O}_2$; (c) ${\rm H}_2$; (d) ${\rm H}_2{\rm O}$; (e) OH; ( f) H; (g) ${\rm HO}_2$; (h) ${\rm H}_2{\rm O}_2$; (i) ${\rm N}_2$.

Figure 26. Heat release rate budgets associated with $\bar {\dot {\omega }}_{T,\mathrm {H}}$, $\bar {\dot {\omega }}_{T,\mathrm {H_2O}}$ and $\bar {\dot {\omega }}_{T,\mathrm {H}}+\bar {\dot {\omega }}_{T,\mathrm {H_2O}}$, along with the total heat release rate $\bar {\dot {\omega }}_T$ at different streamwise positions: (a) $x^*=100$; (b) $x^*=125$; (c) $x^*=150$.

6.3. Turbulence–chemistry interaction

We investigate the turbulence–chemistry interaction by evaluating closure strategies for the mean turbulent heat release rate $\bar {\dot {\omega }}_T(T,Y_i)$. The simplest closure strategy is the laminar chemistry, represented by $\dot {\omega }_T(\bar {T},\bar {Y}_i)$, where the averaged fields of temperature and species concentration are input into the nonlinear expression of the chemical source terms. In RANS calculations, where the average is taken both in time and spanwise directions, using the laminar chemistry can lead to important differences against $\bar {\dot {\omega }}_T(T,Y_i)$ (Poinsot & Veynante 2005).

The mean laminar heat release rate is evaluated at $x^*=68$ and $x^*=150$, where the first position corresponds to the streamwise maximum of $\bar {\dot {\omega }}_T(T,Y_i)$ in figure 6. Considerable differences to more than an order of magnitude are revealed between the turbulent and laminar heat release rates at both positions, as shown in figures 27 and 28. These results indicate that the laminar closure would introduce substantial errors in relevant RANS calculations.

Figure 27. Reference turbulence heat release rate $\bar {\dot {\omega }}_T(T,Y_i)$ (turbulent), heat release rate closed by laminar chemistry $\dot {\omega }_T(\bar {T},\bar {Y}_i)$ (laminar), and mean heat release relate obtained by assumed methods $\breve {\dot {\omega }}_T|_{\bar {Y}_p}$ and $\breve {\dot {\omega }}_T|_{\bar {T}}$ (assumed) at $x^*=68$: (a–g) $\breve {\dot {\omega }}_T|_{\bar {Y}_p}$ for each species, (a) O, (b) ${\rm O}_2$, (c) ${\rm H}_2$, (d) ${\rm H}_2{\rm O}$, (e) OH, ( f) H, (g) ${\rm HO}_2$ and (h) ${\rm H}_2{\rm O}_2$; (i) $\breve {\dot {\omega }}_T|_{\bar {T}}$. S. I. units are used for the rate values.

Figure 28. Reference turbulence heat release rate $\bar {\dot {\omega }}_T(T,Y_i)$ (turbulent), heat release rate closed by laminar chemistry $\dot {\omega }_T(\bar {T},\bar {Y}_i)$ (laminar), and mean heat release relate obtained by assumed methods $\breve {\dot {\omega }}_T|_{\bar {Y}_p}$ and $\breve {\dot {\omega }}_T|_{\bar {T}}$ (assumed) at $x^*=150$: (a–g) $\breve {\dot {\omega }}_T|_{\bar {Y}_p}$ for each species, (a) O, (b) ${\rm O}_2$, (c) ${\rm H}_2$, (d) ${\rm H}_2{\rm O}$, (e) OH, ( f) H, (g) ${\rm HO}_2$ and (h) ${\rm H}_2{\rm O}_2$; (i) $\breve {\dot {\omega }}_T|_{\bar {T}}$. S. I. units are used for the rate values.

The reasons behind this considerable difference requires further exploration: among all the variables in the expression of chemical source terms, which variable's fluctuation contributes most strongly to this difference? Conversely, which species’ fluctuation must be modelled for an accurate closure of the chemical source terms, while the modelling of which species is less necessary? Drawing inspiration from another closure strategy in RANS, the assumed probability density function (p.d.f.) approach (Baurle et al. 1995; Baurle & Girimaji 2003), we propose a simple diagnostic method to examine the importance of modelling a specific fluctuation.

The assumed p.d.f. approach employs an a priori assumed form of the p.d.f. for the temperature and species concentrations to close the chemical source terms. Taking $\bar {\dot {\omega }}_T(T,Y_i)$ as an example, the chemical source term can be represented at each grid point ($x, y$) by integrating the heat release rate of each realisation with the associated joint p.d.f. (j.p.d.f.) of temperature and species concentrations. In this context, the mean heat release rate can be written as

(6.3)\begin{equation} \bar{\dot{\omega}}_T (T,Y_i) = \int \dot{\omega}_T(\hat{T},\hat{Y}_i ) \boldsymbol{P}(\hat{T},\hat{Y}_i ) \,{\rm d} \hat{T} \,{\rm d} \hat{Y}_i, \end{equation}

where $\hat {T}$ and $\hat {Y}_i$ denote the sample space of $T$ and $Y_i$, and $\boldsymbol {P}(\hat {T},\hat {Y}_i )$ represents their associated j.p.d.f. When statistical independence between the variables is assumed, $\boldsymbol {P}(\hat {T},\hat {Y}_i )$ can be represented by the product of individual p.d.f.s for each variable, or the product of several individual p.d.f.s with a j.p.d.f. For example, it is commonly assumed that the temperature follows a Gaussian distribution, which is independent of the concentrations of species assumed to follow a multivariate $\beta$ function (Baurle & Girimaji 2003).

To assess the importance of modelling the fluctuation of species $Y_p$, we consider a scenario where we possess only the knowledge of its mean value $\bar {Y}_p$, while we have complete knowledge of the j.p.d.f. for all other variables, $\boldsymbol {P}(\hat {T},\hat {Y}_1,\ldots,\hat {Y}_{p-1}, \hat {Y}_{p+1},\ldots,\hat {Y}_{N})$, from our DNS data. We define and evaluate an ‘assumed’ mean heat release rate $\breve {\dot {\omega }}_T|_{\bar {Y}_p}$, expressed as follows:

(6.4)\begin{equation} \breve{\dot{\omega}}_T|_{\bar{Y}_p}=\int \dot{\omega}_T(\hat{T},\hat{Y}_i ) \boldsymbol{P}(\hat{T},\hat{Y}_1,\ldots,\hat{Y}_{p-1},\hat{Y}_{p+1},\ldots,\hat{Y}_{N})\delta(\hat{Y}_p-\bar{Y}_p) \,{\rm d} \hat{T} \,{\rm d} \hat{Y}_i, \end{equation}

where $\delta (\hat {Y}_p-\bar {Y}_p)$ is a delta function assuring that the prescribed $Y_p$ can only take its mean value. In a more concise manner, $\breve {\dot {\omega }}_T|_{\bar {Y}_p}$ can be regarded as

(6.5)\begin{equation} \breve{\dot{\omega}}_T|_{\bar{Y}_p}=\bar{\dot{\omega}}_T (T,Y_1,\ldots,Y_{p-1},\bar{Y}_{p},Y_{p+1},\ldots,Y_N), \end{equation}

where $\bar {Y}_{p}$ is inserted into the expression of the heat release rate, replacing each realisation of $Y_{p}$. Similarly, we can assess the importance of modelling temperature fluctuations by introducing $\breve {\dot {\omega }}_T|_{\bar {T}}$, defined as

(6.6)\begin{equation} \breve{\dot{\omega}}_T|_{\bar{T}}=\bar{\dot{\omega}}_T (\bar{T},Y_i). \end{equation}

The results of $\breve {\dot {\omega }}_T|_{\bar {Y}_p}$ for each species excluding $\text {N}_2$, along with $\breve {\dot {\omega }}_T|_{\bar {T}}$, are presented in figures 27 and 28 for $x^*=68$ and $x^*=150$, respectively. The difference between the assumed mean heat release (red dotted) and the turbulent mean heat release (black) indicates the importance of modelling a specific fluctuation. At both positions, the absence of modelling the fluctuations from $\text {H}_2$ leads to the most substantial difference, with the largest discrepancy in the outer layer. This could be attributed to the strong inhomogeneity of $\text {H}_2$ in the outer layer, as evidenced by the instantaneous streamwise section of $\text {H}_2$ in figure 17( f). If the fluctuation of $\text {H}_2$ is not modelled, it can be assumed that $\text {H}_2$ is homogeneously redistributed in the spanwise direction. At local positions where no ignition occurs, the temperature and relevant species concentrations are likely to be sufficient, leading to an excess amount of heat release due to the redistributed $\text {H}_2$. Conversely, the absence of modelling the fluctuations of $\text {O}_2$ results in an underestimation of turbulent heat release rate. At local positions where no ignition occurs, the temperature and species concentrations may be insufficient to generate heat release with the redistributed $\text {O}_2$; at local positions where there is substantial ignition, the amount of $\text {O}_2$ is reduced when redistributed, leading to a reduced heat release.

The no-modelling of the two intermediate products, $\text {OH}$ and $\text {H}$, also results in significant differences with respect to the turbulent mean heat release. In contrast, the no-modelling of $\mathrm {H}_2 \mathrm {O}$, $\mathrm {HO}_{2}$ and $\mathrm {H}_2 \mathrm {O}_{2}$ shows smaller differences, suggesting that modelling the fluctuations of these three species is relatively less critical. The effect of modelling temperature fluctuations is shown by $\breve {\dot {\omega }}_T|_{\bar {T}}$, which reveals a relatively smaller influence at $x^*=150$ compared with $x^*=68$. This can be attributed to the reduced temperature fluctuations at the downstream location compared with the upstream, as evaluated in figure 15(b).

7. Conclusions

We have conducted an analysis of a supersonic turbulent boundary layer with the presence of hydrogen combustion using DNS. The study has primarily focused on the flow dynamics and evaluate the modelling strategies for temperature, species fluctuations and the employed chemical reactions, which are influenced significantly by the wall–jet inlet set-up and the combustion process.

The wall–jet inlet, which includes a cold hydrogen stream near the wall and a preheated air stream far from the wall, triggers a particular transition process. This process is initiated by the mixing-layer-induced instability in the outer region, followed by the development of near-wall fluctuations. The wall–jet transition results in excess Reynolds stresses and turbulent fluctuations, particularly in the outer layer, leading to a streamwise recovery process. Downstream, the Reynolds stresses reduce to values below the canonical adiabatic profiles, attributed to the lower local friction Reynolds number $Re_{\tau }^{\star }$, which is modulated by the mean density and viscosity profiles. The influence of combustion on the recovery process is tracked by the streamwise peak values of Reynolds shear stress in both the reacting and non-reacting cases. Upstream, the combustion effect is primarily manifested by the depletion of hydrogen, leading to an increase of density and the corresponding $Re_{\tau }^{\star }$. Consequently, the Reynolds shear stress peak in the reacting case recovers more slowly than in the non-reacting case. Downstream, the peak attains a plateau region where the reacting case exhibits a lower peak value, which is attributed to the higher temperature due to combustion, resulting in a smaller Reynolds number. We note that although the peak values have reached a plateau, it does not imply that the flow has reached an equilibrium state (Ceci et al. 2022): at an infinitely downstream position, the hydrogen should be entirely depleted, leading to no chemical reaction, and the multi-species gas should be perfectly mixed. In such a case, the temperature of wall surface would reach a constant adiabatic temperature, leading to a closer match with the canonical adiabatic profile.

The impact of combustion on skin friction is examined through a skin friction decomposition in both the reacting and non-reacting cases. With combustion present, the budgets from both the viscous and Reynolds shear stresses are reduced, resulting in skin friction that is approximately 80 % of the non-reacting case.

The turbulent Prandtl number and turbulent Schmidt numbers have been employed extensively in the LES and RANS modelling. We have evaluated their validity using the present reacting DNS database. The turbulent heat transport term $\overline {\rho v'' T''}$ exhibits zero points in wall-normal profiles associated with local maxima and minima in temperature profiles, introducing discontinuities in $Pr_t$. Nevertheless, for both reacting and non-reacting cases, $Pr_t\approx 1.1$ is found in the inner layer and $Pr_t\approx 0.5 \sim 0.9$ in the outer layer, which is consistent with the values and trends in a non-reacting scenario (Peter & Kloker 2022). Chemical reactions lead to slight increase of $Sc_{t,i}$ for $\mathrm {O}_{2}$, $\mathrm {H}_{2}$, $\mathrm {H}_2\mathrm {O}$ and $\mathrm {N}_{2}$ in the outer layer, whereas relatively minor changes are observed for the other species.

The SRA and its extensions have been assessed in both reacting and non-reacting cases. The original correlation identity, with a constant $C_p$, is relatively satisfied in the inner layer, whereas it does not work in the outer layer. Using a varying $C_p$ slightly improves the results in the range of inner layer, but leads to no improvement in the outer layer. A spanwise spectral analysis has been conducted, where the spanwise spectra of temperature and $\mathrm {H}_{2}$ are characterised by peak values in the outer layer, associated with inhomogeneity from the scalar mixing and local heat release. The large-scale scalar fluctuations do not generate significant velocity fluctuations at the same scales, leading to the invalidity of SRA in the outer layer. The spectral linear coherence spectrum between the velocity and temperature fluctuations further supports the relative success of SRA in the inner region, with a coherence value around 0.5 comparable to canonical wall-bounded turbulent flows (Yu & Xu 2021). However, this spectral coherence is lost in the outer region. The extended versions of SRA, including ESRA, GSRA, RSRA and HSRA, exhibit no improvement and even more severe deviations from unity than the original version. This is because these extensions are primarily based on similarity relations between the velocity and temperature profiles, which are not satisfied within the present mixing setting.

Reaction statistics at various wall-normal distances, a focus of this investigation, have been conducted by averaging the reaction source terms, elementary reaction rates and the chemical heat release rates. The near-wall reaction mechanism, particularly involving $\mathrm {H}_{2}$, $\mathrm {H}_2\mathrm {O}$, $\mathrm {H}$ and $\mathrm {OH}$, has been explored through elementary reaction rates. Elementary reactions 3, 5, and 8 among the 19 reactions in Li et al. (2004) are identified as important to the near-wall reaction. The budget analysis of chemical heat release rate highlights the dominant contributions of $\mathrm {H}_2\mathrm {O}$ in the outer layer and $\mathrm {H}$ at the wall to the overall heat release. These findings also suggest that employing a reduced reaction scheme, such as the one-step global reaction $\mathrm {H}_{2}+2\mathrm {O}_{2} \Leftrightarrow 2\mathrm {H}_2\mathrm {O}$, would be insufficient to capture the near-wall production and consumption of species when using an adiabatic and non-catalytic wall condition.

Turbulent chemistry closure in a RANS scenario has been examined by comparing the laminar and turbulent chemistry. The differences between the two are found to be more than an order of magnitude. We use a simple diagnostic method to assess the necessity of modelling the fluctuations from each variable in the expression of chemical source terms, including species concentrations and temperature. The method involves inserting the mean value of a specific variable into the expression of local heat release at each point, while maintaining the local information of other variables, before averaging to obtain an ‘assumed’ mean heat release rate. It has been shown that the absence of modelling the fluctuations of $\mathrm {H}_2$ results in the most severe deviation from the turbulent heat release rate, whereas the absence of modelling $\mathrm {H}$ and $\mathrm {OH}$ also results in relatively considerable deviations.

In conclusion, this study has offered insight into the statistical characteristics of wall-bounded turbulent reacting flows, aligning with recent studies on thermochemical non-equilibrium effects in high-speed turbulent boundary layers. The results from DNS have given practical recommendations for closure models, including the wall-normal profiles of turbulent Prandtl and Schmidt numbers evaluated in a multi-species reacting scenario, the invalidity of SRA in the outer layer due to scalar mixing and suggestions for turbulent chemistry closure.

Further work with this DNS database will focus on the following three orientations. First, since the present work shows the prescribed wall–jet inlet is associated with particular transition behaviours, complementary calculations with inlets of a laminar jet stream and a turbulent main stream will be conducted, to compare the transition and mixing behaviour with the present laminar setting. Second, as the current paper concentrates on flow dynamics and reaction statistics, discussions on the ignition mechanism, which is of significant practical interest, have been omitted. The effects of inflow conditions on the ignition process will be further investigated. Third, as scalar mixing and chemical heat release primarily occur in the outer layer, a more detailed analysis of vorticity and scalar transport near the edge of turbulent boundary layer is expected, where compressible effects associated with thermal expansion or baroclinic torque may play an important role.