1. Introduction
High-speed flows, encompassing objects entering a planetary atmosphere or engaged in atmospheric hypersonic flight, pose a challenging topic of interest for manifold configurations. In the case of such vehicles, the bow shock wave ahead of the body and the boundary layers developing around it transform kinetic energy into thermal and chemical energy, potentially generating thermally imperfect gas effects, and dissociating nitrogen and oxygen molecules (Anderson Reference Anderson2006; Bertin & Cummings Reference Bertin and Cummings2006; Candler Reference Candler2019; Urzay & Di Renzo Reference Urzay and Di Renzo2021). Heat transfer mechanisms in the boundary layer might differ substantially from conventional fluid mechanics, requiring examination of limiting cases to understand various physical processes and compare with experiments. Investigating the intricate dynamics of boundary layers is paramount for the design and performance optimisation of high-speed aircraft and spacecraft. In particular, during the later stages of re-entry trajectories or under suborbital flight conditions with higher free-stream densities, the flow behaviour becomes more intricate due to the potential transition from a laminar to a turbulent flow regime. This transition can be induced by factors such as free-stream disturbances, erosion and ablation of thermal protection systems, or surface defects. As a result, there is a significant increase in skin friction and heat flux at the walls. For these reasons, laminar-to-turbulent transition (Marxen, Iaccarino & Magin Reference Marxen, Iaccarino and Magin2014; Mortensen Reference Mortensen2015; Kline, Chang & Li Reference Kline, Chang and Li2019), together with fully turbulence dynamics more recently (Di Renzo & Urzay Reference Di Renzo and Urzay2021; Passiatore et al. Reference Passiatore, Sciacovelli, Cinnella and Pascazio2021, Reference Passiatore, Sciacovelli, Cinnella and Pascazio2022; Li et al. Reference Li, Yu, Sun, Liu and Yuan2022; Di Renzo, Williams and Pirozzoli Reference Di Renzo, Williams and Pirozzoli2024; Sciacovelli et al. Reference Sciacovelli, Cannici, Passiatore and Cinnella2024), has been identified as a subject of major concern for the accurate prediction of the aerothermodynamic loads on vehicles flying at hypersonic speeds. The complexity of experimentally studying these phenomena makes numerical simulations of high-enthalpy high-speed flows a crucial, albeit complicated, part of this endeavour. From a modelling standpoint, our existing knowledge on supersonic/hypersonic cryogenic turbulent boundary layers must be extended to chemically reacting flows.
Even without considering high-enthalpy effects, high-speed turbulent flows and, specifically, wall-bounded ones, are an active topic of research. A defining characteristic of these flows is the high kinetic energy of the free-stream gas, which has a profound impact on the near-wall turbulence. When such highly energetic flows are decelerated by a wall, they generate significant thermal fluxes and intense pressure waves, presenting substantial technical challenges in selecting appropriate surface materials. Moreover, the pronounced temperature variations lead to significant changes in both density and viscosity. Density fluctuations in response to changes in pressure (intrinsic compressibility effects) are instead usually neglected if Morkovin’s hypothesis holds (Morkovin Reference Morkovin1962), then turbulence statistics in compressible flows can be approximately collapsed onto their incompressible counterparts by taking into account mean property variations (Van Driest Reference Van Driest1951; Zhang et al. Reference Zhang, Bi, Hussain and She2014; Trettel & Larsson Reference Trettel and Larsson2016; Volpiani et al. Reference Volpiani, Iyer, Pirozzoli and Larsson2020; Griffin, Fu & Moin Reference Griffin, Fu and Moin2021; Hasan et al. Reference Hasan, Larsson, Pirozzoli and Pecnik2024). The validity of Morkovin’s hypothesis has been stressed lately by trying to increase wall-cooling and Mach numbers with the aim of ‘breaking’ compressible transformations. Indeed, if the Mach number is significant, intrinsic compressibility becomes important. This may cause the failure of classical incompressible transformations and compressibility corrections should be developed based on intrinsic compressibility rather than variable-property effects (Hasan et al. Reference Hasan, Larsson, Pirozzoli and Pecnik2023).
It is evident that to draw general conclusions and develop the necessary corrections, a comprehensive database of high-speed compressible turbulent boundary layers must be established. This database should encompass a broad range of Mach numbers and wall-to-recovery temperature ratios. While significant progress is being made in this regard for cryogenic configurations (Huang et al. Reference Huang, Nicholson, Duan, Choudhari and Bowersox2020; Xu et al. Reference Xu, Wang, Wan, Yu, Li and Chen2021; Cogo et al. Reference Cogo, Salvadore, Picano and Bernardini2022; Huang, Duan & Choudhari Reference Huang, Duan and Choudhari2022; Wenzel, Gibis & Kloker Reference Wenzel, Gibis and Kloker2022; Cogo et al. Reference Cogo, Baù, Chinappi, Bernardini and Picano2023; Xu, Wang & Chen Reference Xu, Wang and Chen2023; Gibis et al. Reference Gibis, Sciacovelli, Kloker and Wenzel2024), some practical challenges limit the extension of this knowledge to high-enthalpy flows. Specifically, these limitations are connected with the similarity properties of high-enthalpy laminar hypersonic boundary layers. Currently, direct numerical simulations of high-enthalpy turbulent boundary layers are obtained by computing the entire transition process from a laminar state. For this reason, most of the computational cost derives from the use of very large domains, due to the difficulty in achieving a fully developed state. Depending on the desired path to turbulence for the flow to follow, free-stream disturbances, localised wall defects, or suction and blowing strips can be used as tripping. The latter method enables guided excitation of instability waves and is the most commonly used approach to achieve modal transition. However, significant distances are required for these waves to become nonlinearly unstable and trigger the breakdown to turbulence, particularly in hypersonic, wall-cooled boundary layers. Moreover, a common challenge in these configurations is the reliability and repeatability of results, as the primary source of discrepancy across different studies lies in the attainment of a turbulent equilibrium state downstream of the transitional region.
Simulating the entire flow evolution, including the laminar and transitional sections, becomes computationally unfeasible when a sufficiently large, well-developed turbulent region is desired. The most straightforward solution to overcome this issue is to initiate the simulations directly in the turbulent regime. However, this approach comes with its own challenges, such as selecting an appropriate turbulence seeding method and determining a minimal distance from the inflow at which the flow properties become decorrelated (Ceci et al. Reference Ceci, Palumbo, Larsson and Pirozzoli2022). Nonetheless, it enables the execution of fully turbulent simulations, covering a range of moderate to high Reynolds numbers. This approach is commonly used for all the cryogenic studies mentioned previously. When it comes to high-enthalpy boundary layers, where chemical dissociation is triggered, the main drawback lies in the peculiar self-similar properties of these flows, which have, until now, discouraged the use of a turbulent inflow.
Self-similar solutions, in which flow characteristics remain invariant under specific scaling transformations, provide crucial insights into the fundamental nature of the flow and form the basis for developing inflow turbulence generation methods. The introduction of additional characteristic time scales due to chemistry forces the solution to be only locally self-similar (Liñán & Da Riva Reference Liñán and Da Riva1962), or at least this is our understanding to date, which is limited to the laminar regime. This change in the nature of the similarity solutions holds profound importance, and its understanding could pave the way towards a thorough description of the dynamics of hypersonic reacting boundary layers. A preliminary investigation has been performed by Passiatore & Di Renzo (Reference Passiatore and Di Renzo2023), where the impact of different chemical rates and their interplay with diffusion phenomena has been highlighted in laminar boundary layers. The same understanding has not yet been fully established for turbulent flow fields, meaning that no comprehensive investigation into similarity behaviour has been conducted.
The present work aims at bridging this gap with the exploration of the self-similar behaviour of turbulent boundary layers with chemical reactions. To investigate different finite-rate chemistry regimes, a parametric study should be conducted by varying free-stream stagnation enthalpy, wall temperature and Reynolds number. Alternatively, these parameters can be held constant while adjusting the Damköler number by modifying the chemical time scale. To facilitate this, a simplified dissociation–recombination mechanism is employed, with reaction rates deliberately controlled to induce strong chemical activity even under moderate free-stream stagnation enthalpy. Consequently, three direct numerical simulation (DNS) studies are performed with varying reaction rates. Given the motivation described previously and in the view of forthcoming studies, mastering the fundamental aspects of similarity and its implications for chemically reacting turbulent boundary layers is essential. Therefore, these DNS cases are used to examine the behaviour of species transport equations under varying reaction rates and to explore scenarios where the chemical leading edge is displaced relative to the flow leading edge, assessing the potential recovery of a self-similar state.
The remainder of the work is organised as follows. Section 2 reports the governing equations and physical models, with particular attention to the dissociation–recombination mechanism and numerical method. In § 3, the description of the test cases studied in this work is reported, together with some statistics of the reference cases. The results of this work are discussed in § 4 and the main findings are summarised in § 5.
2. Formulation
2.1. Governing equations and physical models
In this work, the Navier–Stokes conservation equation for compressible multi-component flows
\begin{align}&\qquad\qquad\qquad\quad \frac {\partial \rho u_i}{\partial t} + \frac {\partial \left ( \rho u_i u_j + p \delta _{ij} \right )}{\partial x_j} = \frac {\partial \tau _{ij}}{\partial x_j}, \end{align}
\begin{align}& \frac {\partial \rho E}{\partial t} + \frac {\partial \left [\left (\rho E + p \right )u_j\right ]}{\partial x_j} = \frac {\partial (u_i \tau _{ij} - q_j)}{\partial x_j} -\frac {\partial }{\partial x_j}\left (\sum _{n=1}^{\textit{NS}} \rho _n u_{n,j}^D h_n \right ) ,\end{align}
\begin{align}&\qquad\quad \frac {\partial \rho _n}{\partial t} + \frac {\partial (\rho _n u_j)}{\partial x_j} = -\frac {\partial \rho _n u_{n,j}^D }{\partial x_j} + \dot {\omega }_n \quad (n=1,\ldots ,\textit{NS}) ,\end{align}
are solved. In the preceding equations,
$\rho _n$
is the partial density of the
$n$
th species, with
$\rho =\sum _{n=1}^{{NS}} \rho _n$
being the mixture density, NS is the total number of species,
$u_i$
(
$i=1,2,3$
) are the components of the velocity vector in the direction
$x_i$
and
$p$
is the pressure. The system of conservation equations is closed with the equation of state for a multi-component mixture of ideal gases,
\begin{align} \nonumber\\[-24pt]p = R_0 T \sum _{n=1}^{{NS}} \frac {\rho _n}{W_n} = R_0 \rho T \sum _{n=1}^{{NS}} \frac {Y_n}{W_n}, \end{align}
where
$Y_n$
and
$W_n$
are the mass fraction and molar mass of the
$n$
th species, respectively. Here,
$R_0 = 8.314$
J mol−1 K−1 is the universal constant of gases.
In the momentum conservation equation,
$\tau _{ij}$
denotes the viscous stress tensor
\begin{equation} \tau _{ij} = \mu \left (\frac {\partial u_i}{\partial x_j} + \frac {\partial u_j}{\partial x_i} \right ) -\frac {2}{3}\mu \frac {\partial u_k}{\partial x_k}\delta _{ij}, \end{equation}
with
$\mu$
the mixture dynamic viscosity and
$\delta _{ij}$
the Kronecker symbol. Additionally,
$q_j = -\lambda ({\partial T}/{\partial x_j})$
represents the heat flux in the energy conservation equation, where
$\lambda$
is the mixture thermal conductivity and
$T$
the temperature. Sutherland’s law is used to compute viscosity, while the thermal conductivity of the mixture is evaluated with the constant Prandtl number assumption
$Pr = c_p \mu / \lambda = 0.71$
. In the previous expression,
$c_p = \sum _{n=1}^{{NS}}Y_n c_{p,n}$
is the mixture averaged specific heat capacity at constant pressure and
$c_{p,n}$
is the constant-pressure heat capacity of the
$n$
th species. In the energy equation, the stagnation energy is transported,
$E=e + {1}/{2}u_i u_i$
, with
$e$
being the specific internal energy of the mixture defined as
\begin{equation} e = h - p/\rho = \sum _{n=1}^{{NS}}Y_n h_n - p/ \rho , \end{equation}
where
$h$
is the specific enthalpy of the mixture and
$h_n$
is the partial specific enthalpy of species
$n$
given by
\begin{equation} h_n = h_{n,\textit{ref}} + \int _{T_{{ref}}}^T c_{p,n}(T')\, {\rm d} T'. \end{equation}
In this formulation,
$h_{n,{ref}}$
, which is a reference value of the specific enthalpy taken at the reference temperature
$T_{{ref}}$
, and
$c_{p,n}$
are evaluated using the nine-coefficient NASA polynomials (McBride, Zehe & Gordon Reference McBride, Zehe and Gordon2002). In the species conservation equation, the diffusion velocity is defined as
\begin{equation} u_{n,j}^D = {- D \frac {\partial Y_n}{\partial x_j} .} \end{equation}
In the notation,
$D$
is the mixture-averaged mass diffusivity, computed using a constant Schmidt number assumption
$Sc = \mu / ( \rho D ) = 1.0$
.
In the previous formulation, several hypotheses were adopted to model high-temperature chemically reacting gases. Transport properties were determined using low-temperature laws, such as Sutherland’s law and the assumption of a constant Prandtl number. Additionally, mass diffusivity was calculated under the constant Schmidt number hypothesis. These simplifications, along with the neglect of vibrational non-equilibrium effects in the current configuration, represent significant approximations. However, the primary objective of this study is to evaluate the behaviour of dissociation–recombination mechanisms exclusively. Therefore, these modelling choices were made to isolate the effects of non-equilibrium chemistry.
2.2. Chemical mechanism
To simplify the problem and isolate the effects of the chemical mechanism, we are considering a two-species mixture composed of O
$_2$
and O. It can be envisioned as any two-species mixture, subject to the simple dissociation–recombination mechanism, such as
\begin{equation} \textrm{O}_2 \Longleftrightarrow 2 \textrm{O}. \end{equation}
The chemical source term of O
$_2$
is computed following the law of mass action
\begin{equation} \dot {\omega }_{{O}_2}= k_f \left ( \frac {\rho Y_{{O}_2}}{W_{{O}_2}} \right ) - k_b \left (\frac {\rho Y_{{O}}}{W_{{O}}} \right )^2, \end{equation}
whereas
$\dot {\omega }_{{O}} = -\dot {\omega }_{{O}_2}$
. In (2.10),
$k_{f}$
is the forward reaction rate coefficient, computed as
\begin{equation} k_f=A T^n \operatorname {exp}\left (\frac {-E_a}{R_0 T}\right ), \end{equation}
and
$k_b = \displaystyle k_f/K_{{eq}}$
is the backward reaction rate coefficient, with
$K_{{eq}}$
being the equilibrium constant computed as a function of Gibbs free energy. With
$K_{{eq}}$
unequivocally defined by the thermodynamics of the mixture, one can only manipulate the coefficients
$A$
,
$n$
and
$E_a$
in (2.11), corresponding to the pre-exponential factor, temperature exponent and activation energy of a reaction, respectively, to trigger different scenarios. In § 3, we will characterise the values of
$E_a$
and
$A$
for the cases under investigation, while the temperature exponent
$n$
is always equal to −1.5. To have an idea of the order of magnitude of these parameters in a five-species air mixture (Park Reference Park1989), the reaction corresponding to the dissociation of molecular oxygen is characterised by
$A= 2 \times 10^{15}$
$\rm m^3\, (mol \, s)^{-1}$
and
$E_a/R_0=59\,500$
K, in contrast to the dissociation of nitrogen, which has
$A=7 \times 10^{15}$
$\rm m^3\, (mol \, s)^{-1}$
and
$E_a/R_0=1\,13\,200$
K. As a result, molecular oxygen dissociation has a significant reaction rate coefficient at
$T={2550} \, \textrm {K}$
, while molecular nitrogen becomes active at
$T \sim 4000\, \textrm{K}$
. With these magnitudes in mind, the remainder of the paper will present results in which these parameters are varied to trigger specific dissociation rates at moderate stagnation enthalpy values.
2.3. Numerical methods
The transport equations are approximated using finite-difference schemes on a collocated Cartesian computational grid. The convective fluxes are discretised with a hybrid method consisting of a sixth-order skew-symmetric scheme (Pirozzoli Reference Pirozzoli2010) combined with sixth-order targeted essentially non-oscillatory (TENO) shock capturing (Fu, Hu & Adams Reference Fu, Hu and Adams2016). The skew-symmetric formulation is deployed in the smooth regions of the flow, while the TENO reconstruction is used for stencils that cross shock waves or steep gradients. A directional shock sensor controls the activation of the TENO scheme as described by Williams, Di Renzo & Moin (Reference Williams, Di Renzo and Moin2022). The conservation equations are advanced in time using the strong-stability-preserving third-order Runge–Kutta scheme described by Gottlieb, Shu & Tadmor (Reference Gottlieb, Shu and Tadmor2001).
The aforementioned algorithm is implemented in the hypersonics task-based research (HTR) solver (Di Renzo, Fu & Urzay Reference Di Renzo, Fu and Urzay2020; Di Renzo & Pirozzoli Reference Di Renzo and Pirozzoli2021; Di Renzo Reference Di Renzo2022). The HTR solver is based on the task-based environment Legion, which enables efficient memory management along with improved mapping of the workload on computer nodes at runtime (Bauer et al. Reference Bauer, Treichler, Slaughter and Aiken2012). This leads to high scalability in GPU-based supercomputers and allows for the exploitation of heterogeneous architectures (CPU + GPU) without making machine-tailored modifications of the solver. The HTR solver is written in a combination of C++ and the Regent programming language, which is specifically designed to interact with Legion (Slaughter et al. Reference Slaughter, Lee, Treichler, Bauer and Aiken2015)
3. Problem statement
3.1. Computational set-up and flow parameters
The analysed configuration is a spatially evolving, zero-pressure-gradient boundary layer developing over a flat plate. The prescribed edge conditions correspond to
$M_e = U_e/a_e = 10$
based on the edge values of the velocity
$U_e$
and frozen speed of sound
$a_e$
. The edge values of pressure and temperature are
$p_e = 4600 \, \textrm {Pa}$
and
$T_e = 500\, \textrm {K}$
, and pure molecular oxygen is considered at the free stream. The plate’s wall is isothermal at a temperature
$T_w=5 T_e$
. The value of
$T_w/T_e$
has been selected to achieve different regimes of chemistry/transport interaction by just varying the pre-exponential factor
$A$
. While this parameter is not expected to play a major role on the findings of this work because of the simplified chemistry model, its influence should be addressed in future works whereby more complicated chemical processes are considered. In addition, the wall is assumed to be non-catalytic, in such a way that the wall-normal diffusion velocity is infinitely fast at the surface for both species. The computational domain is a rectangle adjacent to the flat plate. The extent of the domain is
$L_x \times L_y \times L_z = 2000\delta ^*_{{in}} \times 100\delta ^*_{{in}} \times 20 \pi \delta ^*_{{in}}$
, with
$x$
,
$y$
and
$z$
denoting the streamwise, wall-normal and spanwise directions, respectively, and
$\delta ^*_{{in}}$
is the displacement thickness of the boundary layer at the inlet section. This displacement thickness corresponds to
$Re_{\delta ^*_{{in}}} = \rho _e U_e \delta ^*_{{in}} / \mu _e = 15\,000$
, where
$\mu _e$
is the dynamic viscosity computed in the edge conditions. The computational grid is discretised with
$N_x \times N_y \times N_z = 8000 \times 480 \times 512$
, for a total of approximately 2 billion grid points. The grid points in the
$x$
- and
$z$
-directions are uniformly distributed, whereas the stretching function
$y_j = 100\delta ^*_{{in}}{\tanh}(\alpha y_j)/{\tanh}(\alpha )$
with
$j=1,2,\ldots, 480$
and
$\alpha = 1.88$
warrants higher wall-normal resolution near the wall. One-dimensional non-reflecting boundary conditions (Poinsot & Lele Reference Poinsot and Lele1992; Okong’o & Bellan Reference Okong’o and Bellan2002) are imposed at the outflow boundary and along the upper surface of the domain. At the inflow, the profiles of temperature, velocity and molar fractions deriving from the locally self-similar form of the laminar boundary-layer equations are imposed (Di Renzo & Urzay Reference Di Renzo and Urzay2021). These equations have been specialised with the same transport properties and chemical mechanism of the present work, similar to Passiatore & Di Renzo (Reference Passiatore and Di Renzo2023).
3.2. Forcing technique
To achieve a fully turbulent state, downstream of the inflow, a suction and blowing strip is used to excite the instability waves in the boundary layer, as in most studies on controlled boundary layer transition in the hypersonic regime. The forcing term consists of a distribution of wall-normal velocity applied on a strip of the wall boundary condition as
\begin{equation} v_w =f(x)g(z)\sum _{m=1}^{N_{{mode}}} A_m\sin (\omega _m t + \beta _m z) ,\end{equation}
where the shape functions are the same as those used by Franko & Lele (Reference Franko and Lele2013), namely
\begin{align}&\qquad\qquad\qquad\qquad f(x) = \exp \left [\displaystyle -\frac {1}{2} \left (\frac {x-x_c}{\sigma }\right )^2\right ], \end{align}
\begin{align}& g(z) = 1 + 0.1\left \{ \exp \left [-\left (\frac {z-z_c-z_w}{z_w}\right )^2\right ]- \exp \left [-\left (\frac {z-z_c+z_w}{z_w}\right )^2\right ]\right \}, \end{align}
where
$x_c=97 \delta ^*_{{in}}$
is the centre of the Gaussian distribution with variance
$\sigma =0.75\delta ^*_{{in}}$
,
$z_c=L_z/2$
and
$z_w = 2 \pi \delta ^*_{{in}}$
. The Reynolds number at the forcing location is
$Re_x \approx 1.5 \times 10^6$
, with
$Re_x= \rho _e u_e x / \mu _e$
being the Reynolds numbers based on the distance from the leading edge
$x$
. A two-dimensional fundamental mode is forced with amplitude corresponding to
$8\, \%$
of the edge velocity. A collection of two-dimensional, three-dimensional and stationary waves are added on top of the fundamental one, with an amplitude corresponding to
$0.4\, \%$
of the edge velocity. The characteristics of the modes are given in table 1. It is now well known that the receptivity of the boundary layer with respect to the imposed waves will be affected by chemistry (Malik & Anderson Reference Malik and Anderson1991; Hudson, Chokani & Candler Reference Hudson, Chokani and Candler1997; Perraud et al. Reference Perraud, Arnal, Dussillols and Thivet1999; Marxen et al. Reference Marxen, Iaccarino and Magin2014; Passiatore et al. Reference Passiatore, Gloerfelt, Sciacovelli, Pascazio and Cinnella2024). In fact, the transitional region extends differently in the cases studied in this work, due to the different chemical activity. Nevertheless, the investigation of the transition behaviour is beyond the scope of the present work and the following analysis is performed in the region where the transition to turbulence has occurred for all the cases under investigation.
Parameters of the modes excited by the forcing function: angular frequency, amplitude and spanwise wave number.
3.3. Presentation of the test cases
We consider three DNS cases in which the activation energy is kept fixed, whereas the pre-exponential factor
$A$
is varied. In this way, dissociation effects become active at similar temperatures in all the presented cases and the comparison is based on the rate of the dissociation–recombination mechanism. Accordingly, these cases differ from the inflow profiles (computed according to the chemical rate) and reaction rates uniquely. The activation energy corresponds to
$E_a/R_0 =50 \times 10^3$
K, whereas the reaction rate
$A$
(with dimension
$[\rm {m^3\, (mol \def\lumina\hspace {0.5mm} s)^{-1}]}$
, omitted hereafter for conciseness) is equal to
$10^{15}$
for the slowest chemistry,
$10^{16}$
for the medium chemistry, and
$10^{17}$
for the fast chemistry.
During our analysis, we compare the reference simulations with companion semi-frozen simulations, where chemical reactions are artificially activated downstream in the fully turbulent region. This is made to shift the origin of the hydrodynamic field with the origin of the thermochemical field. For these cases, the velocity and temperature inflow profiles are kept identical to those of the corresponding reference case, while the inlet molar fractions are forced to be
$X_{{O}_2}=1$
and
$X_{{O}}=0$
. The chemical source terms are artificially set to zero up to
$Re_{x_0} \approx 18 \times 10^6$
, where a fully turbulent state is achieved, and activated downstream, as it will be shown later. These simulations will be referred to as ‘FR’ cases. Additionally, for one mechanism, we have also computed a case for which the atomic oxygen molar fraction is set to zero at the inflow uniquely, but chemistry is immediately switched on from the first computational point. This case, which is similar to the chemical non-equilibrium case presented by Passiatore et al. (Reference Passiatore, Sciacovelli, Cinnella and Pascazio2021), will be named ‘FR Inflow’. For one analysis, we will use two laminar companion simulations at
$A=10^{15}$
, obtained by deactivating the suction and blowing forcing. The first corresponds to the reference case, while the second represents the FR Inflow case. For the sake of clarity, we summarise all the considered cases in table 2.
Summary of the the considered test cases.
Before diving into the analysis of the presented cases, in figure 1, we report the inflow temperature and molar fraction profiles (used for the reference cases), plotted against the similarity variable. Multiple considerations are in order. The profiles of temperature are influenced by the intensity of the chemical source terms. Dissociation reactions are endothermic in the forward direction, draining energy from the boundary layer and resulting in higher temperatures when chemical activity is lower. All the profiles exhibit a non-monotonic behaviour due to wall cooling, but as the rate increases, the profiles are flatter and have lower temperature peaks. As for the evolution of the species molar fractions, we have spanned a concentration range of dissociated species that goes from 8 % up to 16 %, approximately aligned with boundary-layer results present in the literature for the five-species air mixture. The most important finding when going to such high chemical intensity is the non-monotonicity of the profiles according to the rate of the reaction. Since the wall temperature is sufficiently low to induce an equilibrium constant of the chemical reaction that is biased towards the non-dissociated side, the near-wall recombination rate becomes increasingly more intense as the pre-exponential factor increases. Similar insights will be observed in the DNS results.
Inflow similarity profiles, plotted versus the similarity variable: (a) atomic oxygen molar fractions and (b) normalised temperature.
Based on the inflow profiles, we can univocally compute the Damköler number of the three reference simulations, as
\begin{equation} Da = \delta ^*_{{in}}/(U_e t_{c,{in}}), \qquad t_{c,{in}} = \frac { \rho _{e}}{ \dot {\omega }_{O\textit{max},{in}}}. \end{equation}
In the previous equation, the ratio
$\delta ^*_{{in}}/U_e$
is the same for the three cases; therefore,
$Da$
differs in the three cases because of the chemical characteristic time scales uniquely. Accordingly, the computed Damköler numbers are
$4.2 \times 10^{-5}$
,
$1.2 \times 10^{-4}$
and
$1.8 \times 10^{-4}$
for
$A=10^{15}$
,
$A=10^{16}$
and
$A=10^{17}$
, respectively. It is worth noting that the Damköler numbers achieved in these cases are similar to those measured in previous studies involving the five species air mixture (Di Renzo & Urzay Reference Di Renzo and Urzay2021; Passiatore et al. Reference Passiatore, Sciacovelli, Cinnella and Pascazio2021). Moreover, the Damköler number is directly proportional to the rate coefficient of the dissociation–recombination reaction (therefore,
$Da$
and
$A$
can be alternatively used to label each chemical mechanism) and that higher maximum concentrations of atomic oxygen are encountered for higher values of
$Da$
, as observed in figure 1. In general, the correspondence between reaction rate and maximum level of dissociation is not always true. As already observed by Passiatore & Di Renzo (Reference Passiatore and Di Renzo2023), as the reaction rate increases (or the activation energy decreases), chemical activity becomes increasingly more intense and the profiles of atomic molar fraction tend to closely follow the evolution of the temperature profile. In these conditions, the maximum atomic oxygen percentage exhibits a nonlinear dependence on the Damköhler number, driven by the interaction between transport phenomena and chemical production. This nonlinearity becomes even more pronounced with wall cooling, where recombination effects play a significant role in the near-wall region.
3.4. Data collection and boundary layer statistics
The results presented here are obtained by time- and spanwise-averaging DNS solution snapshots of the flow field in the statistically steady state. In particular, the averages are performed by sampling the flow field every ten time steps for more than one computational flow-through time, reaching an adequate statistical convergence of the results. The time- and spanwise-averaged value of any variable
$\phi$
is represented by the overline symbol
$\overline {\phi }$
, whereas
$\phi ' = \phi - \overline {\phi }$
denotes the corresponding fluctuations. Similarly, density-weighted or Favre averages are denoted by the tilde symbol
$\widetilde {\phi }$
, with
$\phi '' = \phi - \widetilde {\phi }$
being the fluctuations around the average. The analysis will mainly focus on the comparison of thermochemical quantities among the different chemical mechanisms and against the corresponding FR simulations in the fully turbulent region. For the sake of completeness, we report in table 3 some boundary layer properties at different stations (used hereafter in the analysis) for the case with
$A=10^{15}$
, to give an idea of the Reynolds number ranges and resolution. These properties slightly change when changing the reaction rate; therefore, only one case is provided. Figure 2 shows the transitional and turbulent structures, coloured by the atomic molar fraction, for the same case. In the following, superscript
$\bullet ^+$
denotes normalisation with respect to the viscous length scale
$l_v = \overline {\mu }_w/(\overline {\rho }_w u_\tau )$
. Unless otherwise specified, the wall-normal evolution of statistics is displayed in inner semi-local units
$y^{*}=\overline {\rho } u_\tau ^{*}y/\overline {\mu }$
, with
$u_\tau ^{*}=\sqrt {\overline {\tau }_w/\overline {\rho }}$
.
Boundary layer properties at five selected streamwise stations, for the case with
$Re_\tau = \rho _w u_\tau \delta / \overline {\mu }_w$
. In the table,
$Re_\theta =\rho _e u_e \theta / \mu _e$
is the Reynolds number based on the local momentum thickness
$\theta =\int _0^\delta \rho u/\rho _e u_e(1-\rho u/\rho _e u_e)\text{ d}y$
and
$Re_{\delta _2}=\mu _e/\overline {\mu }_w Re_\theta$
the incompressible-scaled counterpart.
$Re_\tau = \rho _w u_\tau \delta / \overline {\mu }_w$
is the friction Reynolds number and
$Ma_\tau = u_\tau /\overline {c}_w$
the friction Mach number.
$H=\delta ^{*}/\theta$
is the shape factor and
$B_q=q_w/(\overline {\rho }_w u_\tau \overline {h}_w)$
is the dimensionless heat flux. Lastly,
$\Delta x^+$
,
$\Delta y^+_{{w}}$
,
$\Delta y^+_\delta$
and
$\Delta z^+$
denote the grid sizes in inner variables in the
$x$
-direction,
$y$
-direction at the wall and the boundary layer edge, and in the
$z$
-direction, respectively. The bold Reynolds numbers represent the streamwise stations used in the analysis.
Isosurfaces of Q-criterion corresponding to 0.3, coloured with the local values of atomic oxygen molar fraction, for the case with
$A=10^{15}$
. The side planes are coloured by the magnitude of the density gradient. The three vertical lines mark the streamwise stations used in the analysis.
(a) Skin friction coefficient, (b) wall heat flux, (c) wall-normal profiles of Griffin et al. (Reference Griffin, Fu and Moin2021) velocity transformation and (d) Reynolds stresses, for the three reference cases.
Profiles of Favre-averaged atomic oxygen molar fraction: (a) streamwise evolution at the wall and (b) wall-normal profiles in the fully turbulent region, for the three reference cases. The dash-dotted grey line represents the equilibrium composition computed using the average thermodynamic state of case
$10^{16}$
, at
$Re_{\delta _2}=3000$
.
Before diving into the main findings of this work, an overview of the statistics of the reference cases is presented in this section. The trend of the skin friction coefficient, defined as
$C_f= {2\overline {\tau }_w}/{\rho _e u_e^2}$
, and normalised wall heat flux, with
$q_w = \overline {\lambda \partial T/\partial y}_{|_w}$
, is shown in figure 3(a,b), whereas the GFM transformed velocity (Griffin et al. Reference Griffin, Fu and Moin2021) and Reynolds stresses at
$Re_{\delta _2}=3000$
are shown in figure 3(c,d) for the three considered Damköler numbers. Of note, the wall shear stress in the fully turbulent region, velocity transformation and Reynolds stresses present negligible differences across the three cases, and also when comparing the reference case and the corresponding FR ones (not shown). This result, which is consistent with multiple works in the literature, highlights the weak coupling between the hydrodynamic and chemical fields in hypersonic boundary layers. However, the impact of the chemical mechanism on the wall heat flux is more evident. A nearly frozen chemistry increases the heat flux, as clearly shown by the
$A=10^{15}$
case. The other two mechanisms show a very similar value of heat flux in the fully turbulent region. Figure 4 shows the evolution of the atomic oxygen molar fraction at the wall and the wall-normal profile of the same quantity in a station in the fully turbulent region. The evolution of the case at
$A=10^{15}$
is very similar to the trend of atomic oxygen or nitric oxide in the five-species air mixture computations found in the literature (Di Renzo & Urzay Reference Di Renzo and Urzay2021; Passiatore et al. Reference Passiatore, Gloerfelt, Sciacovelli, Pascazio and Cinnella2024). Specifically, dissociation is enhanced as the flow evolves in the laminar region because of the increased residence time of the gas in the high-temperature region of the boundary layer and more products are transported to the wall from the temperature peak. As the flow transitions to turbulence, an inversion of the trend with respect to the laminar evolution is observed, due to the enhanced mixing with the colder, external layers. Although molecular oxygen dissociation continues downstream, the overall mixture composition remains far from the equilibrium one. The other two mechanisms have an evolution more similar to the wall heat flux. Indeed, dissociation slightly decreases moving downstream in the laminar region, whereas after the transition, the asymptotic value is higher. Also, despite the last case with
$A=10^{17}$
having the highest reaction rate, the amount of chemical products at the wall by the end of the domain is smaller than for the
$A=10^{16}$
case because faster recombination occurs at the wall. As a matter of fact, when extracting the wall-normal distributions of
$X_{{O}}$
in figure 4(b), the profile of
$A=10^{17}$
is non-monotonic (similar to the inflow profiles). Another way to envision this aspect is to compare these profiles against the corresponding equilibrium composition profile evaluated using the mean mixture temperature of case
$A=10^{16}$
. With the data statistically averaged, this equilibrium composition is an approximate, yet meaningful, result. For a simple two-species mechanism, the equilibrium composition is computed as
$K_{{eq}} = (\rho _{{O}}/W_{{O}})^2/(\rho _{{O}_2}/W_{{O}_2})$
,
$K_{{eq}}$
being unequivocally defined by the thermodynamic state of the mixture. The equilibrium profile shows that the atomic oxygen at the wall should not be much higher than 10 %, whereas both
$A=10^{16}$
and
$A=10^{17}$
show a higher concentration at the wall, justifying the negative production rate that will be observed later in figure 11. The intersection of these profiles with the equilibrium one corresponds approximately to the location where the production rate switches sign. All cases are out of equilibrium, and their differing behaviours provide valuable insights into various flow regimes.
4. Discussion
4.1. Budgets of species conservation equation
To investigate the similarity properties of chemically reacting boundary layers, we observed the budgets of the species conservation equation. The interplay between turbulent transport and chemical non-equilibrium is investigated by considering the averaged
$n$
th species conservation equation, which reads
\begin{equation} \frac {\partial \overline {\rho } \widetilde {Y}_n}{\partial t} + \frac {\partial }{\partial x_j} \left (\overline {\rho } \widetilde {u}_j \widetilde {Y}_n \right ) = \frac {\partial }{\partial x_j} \big(\overline { \rho _n u^D_{n,j}}\big) + \overline {\dot {\omega }}_n + \frac {\partial }{ \partial x_j} \overline {\rho u''_j Y''_n}. \end{equation}
On the left-hand side, the time derivative is equal to zero in statistically stationary flows and the second term represents the mean convection. On the right-hand side, the first term represents the mean molecular diffusion, the second is the averaged production rate, and the last is the turbulent transport.
Budgets of species conservation equation for atomic oxygen at
$Re_x = 7 \times 10^{6}$
(laminar region). Mechanism at (a)
$A=10^{15}$
, (b)
$A=10^{16}$
and (c)
$A=10^{17}$
. The terms are normalised with
$\overline {\mu }/(\overline {\rho } \tau _w)$
. The legend is as follows:
$\overline {\rho } \widetilde {u} \partial \widetilde {Y_O}/\partial x$
(
),
$\overline {\rho } \widetilde {v} \partial \widetilde {Y_O}/\partial y$
(
),
$\partial /\partial x (\overline {\rho u'' Y_O''})$
(
),
$\partial /\partial x (\overline {\rho _O u^D_{O,x}})$
(
),
$\partial /\partial y (\overline {\rho v'' Y_O''})$
(
),
$\partial /\partial y (\overline {\rho _O u^D_{O,y}})$
(
),
$\overline {\dot {\omega }_O}$
(
).
Figure 5 shows the profiles of the terms in (4.1) for atomic oxygen extracted at a station in the laminar region of the flow at
$Re_x = 7 \times 10^{6}$
. Here, the turbulent transport and streamwise diffusion terms are negligible and all the other terms in the equations are of the same order of magnitude. This result is entirely consistent with the locally self-similar theory for laminar hypersonic boundary layers, whereby streamwise convection plays a significant role in the species transport. Therefore, the production rates will approximately scale with the residence time of the gas on the flat plate and history effects will be relevant for this flow.
Budgets of species conservation equation for atomic oxygen at
$Re_{\delta _2} = 3000$
, in (a,c,e) inner scaling and (b,d,f) outer scaling. Mechanism at (a,b)
$A=10^{15}$
; (c,d)
$A=10^{16}$
; and (e,f)
$A=10^{17}$
. The terms are normalised with
$\overline {\mu }/(\overline {\rho } \tau _w)$
. The legend is as follows:
$\overline {\rho } \widetilde {u} \partial \widetilde {Y_O}/\partial x$
(
),
$\overline {\rho } \widetilde {v} \partial \widetilde {Y_O}/\partial y$
(
),
$\partial /\partial x (\overline {\rho u'' Y_O''})$
(
),
$\partial /\partial x (\overline {\rho _O u^D_{O,x}})$
(
),
$\partial /\partial y (\overline {\rho v'' Y_O''})$
(
),
$\partial /\partial y (\overline {\rho _O u^D_{O,y}})$
(
),
$\overline {\dot {\omega }_O}$
(
).
Figure 6 shows the profiles of the budget terms for atomic oxygen for the last streamwise location reported in table 3 in inner scaling (panels a,c,e) and outer scaling (panels b,d,f) for
$A=10^{15}$
,
$A=10^{16}$
and
$A=10^{17}$
in the first, second and third row, respectively. Close to the wall, regardless of the reaction rates, chemical source terms are balanced by the wall-normal diffusion. Turbulent transport grows as well as the molecular diffusion moving away from the wall and they both contribute to balance the chemical source terms. By the edge of the boundary layer, the molecular diffusion decreases and the budgets show a balance between wall-normal turbulent transport and chemical source terms. From the visualisation in outer scaling, one can better appreciate that diffusion plays an important role only in the region very close to the wall, while most of the species transport is dominated by wall-normal turbulent transport. The peak chemical production is around
$y^* \approx 10$
for all the mechanisms considered. From this location, dissociated species are transported away from the peak, primarily towards the wall by molecular diffusion and away from the wall by turbulent mixing. This behaviour is consistent with previous studies involving more complex chemistry, and reflects the relative dominance of different transport mechanisms depending on the wall-normal position and local gradients. The location of the peak production defines the centre of chemical activity and thus serves as a meaningful characteristic time scale. We can identify three regions. Between the wall and the peak of production, chemical rates are balanced by the wall-normal diffusion (chemical viscous layer); for
$ 15 \lesssim y^*$
, chemical source terms are balanced by the wall-normal turbulent transport (chemical outer layer). Between these two regions, chemical source terms are balanced by the combination of molecular diffusion and turbulent transport, defining a chemical buffer layer. Another important aspect refers to the variation in the streamwise direction. No streamwise flux contribution to the budget is comparable with the source terms, meaning that the state of the boundary layer depends only on its local Reynolds number and history effects are negligible. For this reason, one could infer that, as opposed to the well-known locally self-similar approximation for the laminar regime, the solution of the boundary layer is approximately globally self-similar in the well-developed fully turbulent region. In these conditions, the chemical source terms should not scale with the residence time of the gas on the flat plate, as suggested by the laminar theory, but with a local time scale proportional to turbulent mixing characteristic times. As an example to corroborate this statement, we compare two normalisation methods for chemical production rates in the slow chemistry case, at two Reynolds numbers extracted in the fully turbulent region. Figure 7(a) depicts the chemical source profiles normalised with the distance from the leading edge and the edge velocity. This scaling fails to collapse the profiles. However, the normalisation based on the local inner units shown in figure 7(b), which is better for considering local turbulent mixing, improves the collapse. This confirms that production rates depend on local turbulent mixing when a fully turbulent state is reached and that an approximate global self-similar solution of the species turbulent boundary layer can be sought.
Production rates of the
$A=10^{15}$
case normalised by (a) outer variables and (b) inner ones at two different Reynolds numbers.
Streamwise evolution of Favre-averaged atomic oxygen molar fraction of
$A=10^{15}$
in the reference and FR Inflow cases. Profiles extracted for (a,c) laminar and (b,d) turbulent calculations at (a,b) the wall, (c)
$y/\delta = 0.5$
and (d)
$y/\delta = 0.05$
.
4.2. History effects
Another perspective on the dual nature of chemistry in laminar versus turbulent regimes is through the quantification of history effects. At this point, we will use the two laminar simulations of
$A=10^{15}$
, described in table 2. Figure 8 shows the streamwise evolution of the atomic oxygen molar fraction for both the reference and frozen-chemistry inflow profiles together with their companion laminar simulations. The first portion of the streamwise domain is not shown due to the oscillations in composition generated by the suction and blowing forcing in the transitional calculations. Notably, the dissociated layer in the laminar regime is much thicker than in the turbulent region. In the latter, it extends only up to
${\approx}0.07 \delta$
, whereas the laminar species layer thickness nearly coincides with the momentum boundary layer one. For this reason, comparisons are made at different heights from the wall. The influence of the frozen-chemistry inflow profile is significant in the laminar calculations (figures 8
a and 8
c). At the wall, the frozen-chemistry inflow laminar simulation is able to approach the reference molar fraction value only by the end of the computational domain. At
$0.5 \delta$
, where the influence of the streamwise convective fluxes is large, as shown in figure 5, the gap between the two laminar calculations persists throughout the entire computational domain. In contrast, the profiles extracted from the turbulent calculations, shown in figures 8(b) and 8(d), demonstrate that the influence of the frozen-inflow profile is confined to the beginning of the computational domain before the transition to turbulence takes place. As soon as the turbulent transition occurs, the data for the FR Inflow case immediately collapses onto the original computation. This is encouraging, as it indicates that any deviation from the baseline in the fully turbulent region is quickly absorbed, allowing the reference to be perfectly recovered and history effects to be completely lost.
Snapshot of
$X_{{O}}$
in an instantaneous
$xy$
-plane for the three reference cases and for the corresponding FR ones. On the same planes, the density gradient is shown with a cutoff below 0.1 and above 0.25. The
$y$
-axis has been stretched for visualisation purposes.
4.3. Impact of delayed chemical activation in the turbulent flow
The results from the previous section indicate that the gap in oxygen molar fraction between the reference case and the FR Inflow case rapidly decays shortly after the turbulent transition. In this section, to further characterise the coupling between the hydrodynamic and chemical fields, the activation of chemical reactions is directly shifted to the fully turbulent region. In this way, the gap between the reference and FR solution is maximised in the turbulent region and it is easier to characterise the decay of the error between the two solutions. Figure 9 presents an instantaneous
$xy$
-plane, extracted after the transitional zone, for reference and FR cases, at
$A=10^{15}$
,
$A=10^{16}$
and
$A=10^{17}$
, arranged sequentially. The contours depict atomic oxygen and the normalised density gradient, with thresholds set below 0.1 and above 0.25. The superimposed density gradient contours indicate that the FR cases exhibit behaviour similar to the reference cases regarding boundary layer thickness and turbulent scales regardless of chemical activity. Once chemistry is initiated at
$Re_{x_0} \approx 18 \times 10^6$
, a chemical boundary layer gradually forms in all three cases. Beyond the initial transient phase, the evolution of
$X_{{O}}$
demonstrates strong qualitative agreement between the reference and FR cases. From a quantitative standpoint, we report the streamwise evolution of the atomic oxygen molar fraction in figure 10, at four different heights: at the wall, at
$y/\delta = 0.003$
, at
$y/\delta = 0.03$
and at
$y/\delta = 0.05$
. Apart from the laminar and transitional region in which no meaningful comparison can be done between the reference and FR cases, after
$Re_{x_0}$
, several considerations can be made. For the three cases, as soon as chemistry is switched on, a very steep jump in the atomic concentration is observed. After the gap between the atomic oxygen molar fraction decreases, the behaviour is nearly asymptotic and a faster reaction rates results in a smaller gap. However, no chemical mechanism is able to perfectly collapse with the reference case within this streamwise domain. Figure 11 shows the evolution of the averaged mass fraction of the atomic oxygen and its fluctuations. The first streamwise station is extracted right after the activation of the production rates for the FR cases (at the same
$Re_x$
to match the same computational station among the three cases). Indeed, at
$Re_x = 18.6 \times 10^6$
, the gap between the mean composition profiles extracted in the reference cases and the FR ones is large, as expected. As the overall structure of the species dissociation layer is rapidly recovered after the chemistry activation and the statistics of the velocity fluctuations are similar in the reference and FR cases, the turbulent mixing of molecular and atomic oxygen of FR cases rapidly tends to approach the reference cases, as shown by the Favre-averaged fluctuation profiles of
$\displaystyle \sqrt {\widetilde {Y^{\prime \prime 2}_{{O}}}}$
shown in figures 11(b), 11(d) and 11(f). The low concentration of atomic oxygen has the effect of increasing the mean temperature, as shown in figure 12(a), and shifting the dissociation reaction towards its forward direction. Regarding the source terms, the values of FR cases are higher than the reference ones in figure 12(b), since the flow in FR cases is trying to compensate for the lack of dissociation. At
$Re_{\delta _2} = 2500$
, the profiles of the mass fraction have relaxed towards the reference ones, without perfectly collapsing them, whereas the source terms are almost perfectly superposed. The collapse of the source terms is favoured by the weak coupling between the hydrodynamic fields of the gas composition. In fact, the Sutherland law and constant Prandtl number assumptions used in this formulation contribute to decreasing the influence of the gas composition on the momentum and energy transport. For this reason, the mean temperature field is rapidly recovered by the FR cases after the dissociation chemistry is activated, as shown by figures 12(c) and 12(e). Similar profiles of species fluctuations in conjunction with similar profiles of velocity fluctuations determine the rapid recovery of the wall-normal species turbulent fluxes, which in turn set the chemical production profiles (figures 12
d and 12
f) through the species transport budget discussed in the previous subsection.
Streamwise evolution of Favre-averaged atomic oxygen molar fraction at (a) the wall, (b)
$y/\delta\,{=}\,0.003$
, (c)
$y/\delta = 0.03$
and (d)
$y/\delta = 0.05$
, for the three reference cases and the corresponding FR ones.
Wall-normal profiles of (a,c,e) Favre-averaged mass fraction and (b,d,f) mass fraction fluctuations of atomic oxygen, at three different Reynolds numbers in the reference and FR cases.
The rapid collapse of the chemical source terms of the FR cases on the reference ones is probably favoured in the presented set-up by the weaker coupling of the hydrodynamic and thermal field with the gas composition field for the sake of simplicity considered in this study. The inclusion of models that evaluate the transport properties of the mixture using the mixture composition will add additional couplings in the system and potentially delay the recovery of the dissociation rate field in the FR cases. Future studies will require the assessment of these aspects in configurations involving more complex thermochemistry.
Profiles of (a,c,e) mean temperature and (b,d,f) normalised chemical production rates at three different Reynolds numbers in the reference and FR cases.
Profiles of temperature fluctuations at two Reynolds numbers for the reference and FR cases.
Figure 13 shows the wall-normal profiles of temperature fluctuations in the analysed test cases. As discussed by Duan & Martín (Reference Duan and Martín2011) and Passiatore et al. (Reference Passiatore, Sciacovelli, Cinnella and Pascazio2021), the intensity of the chemical reaction modulates the amplitude of the temperature fluctuations by converting fluctuations of the internal energy into changes of the gas composition. In fact, the case with slower chemistry shows the most intense temperature fluctuations and
$\widetilde {T''^2}$
decreases as the pre-exponential increases for
$10 \lesssim y^* \lesssim 100$
. The peak of
$\widetilde {T''^2}$
is located for all cases at
$ y^* \sim 10^3$
, whereby the fluid coming from the free stream is mixed by turbulence with the hot dissociated gas generated within the boundary layer. The temperature fluctuations of the FR cases are very similar to the reference profiles for
$y^* \lesssim 100$
, where chemistry is active, and show lower values than the reference close to the outer peak as consequence of the different composition profiles achieved by the gas. As the profiles are sampled more downstream the activation of the chemistry, the
$\widetilde {T''^2}$
profiles become more similar between the FR and corresponding reference case.
Relative error of atomic oxygen molar fraction
$\varDelta ^r_{\widetilde {X_{{O}}}} = ({\widetilde {X_{{O}}}} - {\widetilde {X_{{O}}}}^{{FR}})/{\widetilde {X_{{O}}}}$
at (a) the wall and (b)
$y/\delta = 0.05$
, plotted versus the
$\hat {x}$
coordinate, normalised by the boundary layer thickness at
$x=x_0$
. In panel (b), the dash-dotted line represents the same
$\varDelta ^r_{\widetilde {X_{{O}}}}$
computed with the cases
$A=10^{15}$
and
$A=10^{15}$
, FR inflow.
4.4. Relaxation length
In the attempt to quantify the time (or length) required for the FR cases to collapse with the evolution of the reference cases – or, alternatively, for the difference in atomic oxygen molar fraction between the two cases to reach zero – we consider the portion of computational domain starting at
$x_0$
(distance from the leading edge corresponding to
$Re_{x_0} \approx 18 \times 10^6$
). We can define the relative error in the molar fraction of atomic oxygen as
\begin{equation} \varDelta ^r_{\widetilde {X_{{O}}}} = \frac {\left({\widetilde {X_{{O}}}} - {\widetilde {X_{{O}}}}^{\textit{FR}}\right)}{{\widetilde {X_{{O}}}}_{|_{x_0}}} \end{equation}
and we can plot it versus the
$x$
-coordinate in the computational domain, normalised by a reference boundary-layer thickness. Specifically, if our computational domain virtually starts at
$x_0$
, we can define the length of the streamwise domain from
$x_0$
to the outflow as
$\hat {L}_x$
. Accordingly, we will denote the streamwise coordinate as
$\hat {x}$
. By normalising
$\hat {x}$
with the boundary layer thickness at
$x_0$
,
$\delta _{x_0}$
, we can observe how many ‘inflow’
$\delta _{x_0}$
are needed for the relative error to approach zero. In other words, if this represents the domain of a new simulation starting at
$Re_{x_0}$
, we can estimate how many
$\delta _{x_0}$
are required for the solution to recover the baseline, depending on the chemical mechanism, particularly under the frozen-chemistry inflow profile assumption. This trend is depicted in figure 14. For the fastest chemistry,
$100\delta _{x_0}$
is nearly sufficient for the flow field to reach a thermochemically physical state at the wall, although larger discrepancies remain in the outer region. This trend can be explained by considering that the error between the baseline and the FR simulation comprises the errors introduced in the convection, diffusion and chemical terms. At the wall, with the temperature fixed, infinitely fast diffusion and convection equal to zero, the relative error in the molar fraction is entirely due to the diffusion-dominated near-wall region; therefore, the reduced number of effective degrees of freedom favours a faster recovery of the composition. For the slowest chemistry, the error is larger throughout the boundary layer and an intermediate solution is found for
$A=10^{16}$
. Figure 14(b) also shows the relative error computed with the FR Inflow case for the
$A=10^{15}$
scenario (dash-dotted line). This profile has small errors already at
$x_0$
, as previously shown in § 4.2, and relative differences with the reference case become nearly negligible quickly thereafter.
Going back to the definition of the virtual origin at
$x_0$
, one can envision that this
$x_0$
is the inflow boundary of a fully turbulent simulation. In this set-up, the temperature and velocity profiles are prescribed to align with exact DNS results. However, a key simplification is made for the atomic oxygen profile, which is initially assumed to be zero at the inflow, effectively neglecting dissociation effects. While this approximation ensures good accuracy in the description of the hydrodynamic fields, it introduces a significant thermochemical limitation, as it represents the worst-case scenario where inflow dissociation is entirely disregarded. Nevertheless, the earlier discussion, particularly the comparison against the FR inflow simulation for the
$A=10^{15}$
case, suggests that avoiding the worst-case scenario (
$X_{{O}} = 0$
) and, instead, prescribing a partially dissociated mixture could lead to a faster relaxation of the solution to a physical state. In fact, a comparison between the solid and dash-dotted black lines of figure 14(b) indicates that the relative error at the end of the domain decreases by approximately 67 %. In other words, even if the exact molar fraction evolution at
$x_0$
cannot be provided, introducing a tentative profile with some level of dissociation may expedite the solution’s relaxation to the reference results.
Another way to approach this problem is by plotting the relative error against a chemically relevant length scale. In the previous section, the analysis of the budget terms in the species conservation equation provided valuable insight into how production rates are governed by chemical time scales near the peak of production. Based on this, a chemical length scale can be defined by considering the maximum value of the chemical source term in the wall-normal direction at the location
$x_0$
, where chemistry becomes active in the FR case, multiplied by a characteristic velocity at the same streamwise position:
\begin{equation} l_c = (\rho _{\mid \dot {\omega }_{O\, max }}/\dot {\omega }_{O\, max }) u_\tau \quad \textrm{at} \; x =x_0. \end{equation}
The friction velocity is considered here as a characteristic velocity consistently with the analysis performed in § 4.1. It is worth noting that the maximum chemical source term is also justified by the non-catalytic-wall assumption. In non-catalytic walls, local chemical reactions are neglected and changes in the composition are only due to the fast diffusion transport from the peak of dissociation happening in the viscous sub-layer (Di Renzo & Urzay Reference Di Renzo and Urzay2021). For this reason, the relevant chemical production rate is obtained at the near-wall dissociation peak. The streamwise coordinate
$\hat {x}$
can be normalised by the length scale obtaining
\begin{equation} \Delta x_c = \hat {x}/l_c. \end{equation}
The relative error, as defined in (4.2), is plotted against the streamwise coordinate normalised by the chemical length scale in figure 15 for six wall-normal stations. Figure 6(a,c,e) provides context regarding the distance from the wall, and the balance between production rates and transport at each location. This scaling successfully collapses results across all chemical mechanisms, indicating that the relaxation process is governed by the same dynamics regardless of the reaction rate. Furthermore, all selected wall-normal locations show a consistent collapse, demonstrating that whether in the viscous sublayer, buffer layer or outer region, the dominant mechanisms driving the relaxation process remain the same. Only very close to the edge of the chemical outer layer, the three curves start to deviate slightly. For the fastest chemistry, the relative error decays over spatial scales smaller than the overall domain length. As the pre-exponential factor decreases, the characteristic time scale increases, eventually becoming comparable to the flow crossing time.
To conclude, knowing a characteristic velocity (e.g.
$u_\tau$
) and the peak rate of chemical production (i.e. Damköler number), one can estimate a priori what is the required extension of the computational domain necessary to have a relative error that reaches zero if, at the inflow, the hydrodynamic field is the exact one and a frozen-chemistry assumption is adopted. Of course, the length of this domain can be reduced either if the reaction rate is faster or if a better inflow profile is used.
Relative error of atomic oxygen molar fraction
$\varDelta ^r_{\widetilde {X_{{O}}}}\,{=}\,({\widetilde {X_{{O}}}}\,{-}\,{\widetilde {X_{{O}}}}^{{FR}})/{\widetilde {X_{{O}}}}$
at (a) the wall, (b)
$y^*\,{=}\,2$
, (c)
$y^* = 4$
, (d)
$y^* = 8$
, (e)
$y^* = 15$
and (f)
$y^* = 30$
, plotted against the streamwise coordinate defined in (4.4).
4.5. Constant shear-layer equations
To conclude the discussion in the previous section, we use the constant shear-layer equations to estimate a possible candidate profile to be imposed at the inflow of a potentially fully turbulent calculation starting at
$x_0$
. This is a reasonable choice given the approximate global similarity in the fully turbulent region. We have proved that
$X_{{O}}=0$
at
$x_0$
is a poor approximation, which nevertheless provides a good collapse to the physical profile for the fast chemistry in the prescribed computational domain. This means that a better approximation of
$X_{{O}}$
at
$x_0$
will necessarily provide a better solution in approximately
$100 \delta _{x_0}$
, even for the slower mechanisms. The formulation of these equations for a chemically reacting flow can be obtained by imposing the constant shear layer flow conditions and using an eddy-viscosity model to represent the Reynolds stress tensor (Di Renzo & Urzay Reference Di Renzo and Urzay2019; Urzay & Di Renzo Reference Urzay and Di Renzo2021). In particular, for the species partial density equations, referring to the atomic oxygen uniquely, the formulation reads as
\begin{equation} \frac {{d}}{\textrm{d}y}\left (-\overline { \rho _{{O}} u^D_{{O},y}} + \frac {\mu _t}{Sc_t} \frac {{\rm d} \widetilde {Y_{{O}}}}{\text{dy}} \right ) + \overline {\dot {\omega }_{{O}}} = 0. \end{equation}
Turbulent Schmidt number for the three cases at
$Re_x = 18.6 \times 10 ^6$
.
Some adjustments to the diffusive term and production rate should be made. When solving this equation, we do not usually have access to the diffusive flux or the averaged chemical source terms. Therefore, they can be estimated as follows:
\begin{equation} \overline { \rho _{{O}} u^D_{{O},y}} \approx -\overline {\rho _O} \overline {D} \frac {\text{d} \widetilde {Y_{{O}}}}{\text{dy}} \end{equation}
and
\begin{align} \nonumber\\[-24pt]\overline {\dot {\omega }_{{O}}} \approx \dot {\omega }_{{O}} (\widetilde {T},\overline {\rho _O},\overline {\rho _{O_2}} ). \end{align}
While the impact of the first approximation is expected to be negligible, using a mean-field closure of
$\overline {\dot {\omega }_{{O}}}$
may be a significant approximation (Williams, Di Renzo & Moin Reference Williams, Di Renzo and Moin2024, Reference Williams, Di Renzo and Moin2025). In the discussed formulation, the turbulent viscosity is computed with a mixing length model based on a semi-local scaling
\begin{equation} \mu _t = \kappa \overline {\rho } y \sqrt {\frac {\tau _w}{\overline {\rho }}} \left [ 1 - \text{exp}\left (\frac {- y^*}{A^+} \right ) \right ]^2, \end{equation}
with
$\kappa =0.41$
and the van Direst damping factor
$A^+=17$
. As for the turbulent Schmidt number
$Sc_t$
, a constant value close to one is usually adopted. To verify this approximation, the
$Sc_t$
has been computed from the DNS as
\begin{equation} Sc_t=\frac {\overline {\rho u''v''} \partial \widetilde {Y}_{{O}}/\partial y}{\overline {\rho v''Y_{{O}}''} \partial \widetilde {u}/\partial y} \end{equation}
and presented in figure 16 for the three different chemistries. The turbulent Schmidt number in the outer region of the boundary layer is close to unity; therefore, a turbulent Schmidt number equal to 0.9 is considered.
(a,c,e) Comparison between exact terms from the budgets of species equations and approximated terms computed with the constant shear-layer equations and (b,d,f) comparison between DNS profiles of atomic oxygen and solution of constant shear-layer (CSL) (4.5). (a,b) Slow chemistry, (c,d) medium chemistry and (e,f) fast chemistry, at
$Re_{x_0}$
.
An a priori analysis of the constant shear layer equations has been carried out using the DNS data. Figures 17(a), 17(c) and 17(e) show the results of this a priori analysis performed at
$Re_{x_0}$
(virtual turbulent inflow) for the slow, medium and fast chemistry, respectively. As expected, the molecular diffusion term perfectly collapses with the exact one, whereas the error induced by the eddy-viscosity model and by neglecting the fluctuations due to the nonlinearity of
$\dot {\omega }$
are important. Of note, for the fast chemistry, the production terms computed with the averaged field show a recombination in the outer part of the boundary layer which is not visible in the DNS results. That said, (4.5) can be solved to obtain the profiles of atomic oxygen mass fraction imposing the mean temperature profile and
$\tau _w$
from the DNS. In practical application, the shear stress at the inflow will be an input of the problem and the thermal inflow profile could be estimated using methods proposed over the years for low-temperature compressible boundary layers (Hasan et al. Reference Hasan, Larsson, Pirozzoli and Pecnik2024). The solution of the partial density equation is provided in figure 17(b,d,f) and compared with the corresponding DNS profiles. The estimation of the shear-layer equations shows an important discrepancy, which is attenuated, at least in the outer part of the boundary layer, when increasing the rate of the mechanism. Nevertheless, for the goal of the present work, this profile may be a good approximation of inflow chemical composition profiles for a direct numerical simulation. The discrepancy at the wall will be rapidly attenuated, as previously shown by the evolution of the relative error, and at increasing
$\delta$
, this solution is for sure better than imposing a non-dissociated profile at the turbulent inflow of the simulation. It can be expected that a physical solution could be recovered shortly downstream in inflow if an appropriate turbulent inflow is used for the hydrodynamics. It is worth specifying that while the proposed profiles demonstrate potential as mean inflow profiles for turbulent simulations, they are not yet sufficient on their own to set-up a fully turbulent simulation. At this stage, it remains unclear which method should be used or how to superimpose fluctuations onto the primitive field for chemically reacting flows, but this is under investigation.
5. Conclusion
This work presents a comprehensive investigation into the similarity properties of hypersonic turbulent boundary layers, using a two-species dissociation–recombination mechanism in direct numerical simulation (DNS) studies across various hypersonic scenarios. We focused on three DNS cases, where the reaction mechanism is manipulated maintaining a constant activation energy while varying the pre-exponential factor
$A$
of the Arrhenius reaction to amplify or reduce the dissociation/recombination reaction rate. This approach allows us to initiate dissociation effects at the same temperature across all cases, enabling a comparative analysis based on the velocity of the dissociation–recombination mechanism uniquely. Additionally, we contrast these reference simulations with companion simulations, wherein chemical reactions are activated downstream within the fully turbulent region at
$Re_{x_0}$
. These simulations represent numerical set-ups in which the computational domain virtually starts at
$Re_{x_0}$
in a fully turbulent state, with the dynamic quantities accurately represented at the inflow and the chemical inflow profile deriving from a frozen-chemistry assumption. The hydrodynamic field in the fully turbulent region remains mostly unchanged by different chemical reaction rates and is consistent with the field obtained under the frozen-chemistry assumption. The temperature profiles are slightly altered, with faster reaction rates presenting the smallest temperature peaks. The wall heat flux is notably influenced by chemical mechanisms, with higher heat fluxes achieved by the nearly frozen chemistry, while the exothermic recombination in the
$A=10^{17}$
case aligns it closer to the intermediate reaction rate. The evolution of atomic oxygen at the wall shows enhanced dissociation in laminar flow at increasing Reynolds number, which reverses in turbulence due to mixing with cooler layers for
$A=10^{15}$
. In contrast, the opposite behaviour is observed for the other two cases, where the asymptotic value of
$X_{{O}}$
at the wall in the fully turbulent region is greater than in the laminar flow. The higher wall concentrations in the
$A=10^{17}$
case lead to negative production rates and enhanced recombination at the wall, with a similar, albeit to a lesser extent, effect observed for
$A=10^{16}$
. When comparing the reference simulations against the FR ones, we found that the chemical source terms rapidly converge to the reference values within a small transitional region beginning at
$Re_{x_0}$
. Meanwhile, the profiles of molar fraction significantly reduce their gap in the same region, subsequently decreasing asymptotically at a slope dependent on the Damköhler number. The behaviour of the chemical mechanism in the fully turbulent region is approximately globally self-similar, in contrast to what is observed in the laminar region. Depending on the distance from the wall, the source terms are balanced by laminar wall-normal transport, turbulent wall-normal transport or a combination of both. The relative error of the molar fraction due to delayed activation of the chemistry can be plotted versus a chemical length scale that collapses all chemical mechanisms. In this way, it is possible to assess how many thicknesses of a hypothetical inflow boundary layer with
$X_{{O}}=0$
are necessary to recover the physical behaviour of a simulation that considers an inflow with consistent chemical dissociation. At the same time, a tentative profile with some level of dissociation will expedite the solution’s alignment with the reference results. Since the source terms are uniquely dependent on the wall-normal variable, as shown in the budget evolution, a good approximation of a local tentative profile for the species can be obtained by solving the constant shear-layer equations.
This work represents a first step towards a broader investigation of relaxation phenomena in high-enthalpy boundary layers. By employing a simplified two-species model and basic transport properties, we created a controlled framework to isolate and characterise the influence of chemical reactions on the flow recovery process. These results provide a foundation for ongoing studies, where a full five-species air mixture, along with detailed multi-component transport models, will be incorporated to better capture the complexity of real high-enthalpy environments. Moving to a realistic air chemistry will introduce new challenges: each reaction will have different activation energies and the coexistence of multiple Damköhler numbers will lead to overlapping chemical regimes within the same test case. Similarly, more accurate transport models, accounting for species-dependent viscosity, thermal conductivity and differential diffusion, will introduce stronger couplings between the chemical, thermal and hydrodynamic fields. These complexities will make it more difficult to systematically control the flow regimes, but will also enable a more faithful description of recovery processes relevant to practical applications. While the simplified assumptions adopted here reduce physical realism, they provide a necessary baseline that will guide the interpretation of the more intricate behaviour expected in the fully coupled case. In particular, we anticipate that the increased coupling in the real gas model may modify relaxation lengths and flow recovery rates, although previous studies suggest that the overall influence of chemistry on the main hydrodynamic structures will remain relatively moderate. The extension of the present work to more complex boundary conditions involving catalysis, ablation or radiation will also be investigated in future work.
Acknowledgements.
We acknowledge the CINECA award under the ISCRA initiative, for the availability of high-performance computing resources and support.
Funding.
This material is based upon work supported by the Air Force Office of Scientific Research under award number FA8655-24-1-7001.
Declaration of interests.
The authors report no conflict of interest.
Appendix A. Verifications on equilibrium turbulence
To prove that that an equilibrium turbulent state of the boundary layer is achieved, we have adopted some metrics already used in the literature on cryogenic boundary. In figure 18, we report three metrics. In panel (a), we check that the residual of the compressible equivalent of the von Kármán equation is smaller than a prescribed tolerance, the relative error being defined as
\begin{align} \nonumber\\[-24pt]E_{VK} = \frac {C_f/2-\text{d}\theta /\text{d}x}{C_f/2}. \end{align}
In panel (b), we verify that the friction coefficient collapses onto a power-law relationship versus
$Re_\theta$
,
\begin{align} \nonumber\\[-24pt]C_f = A Re_\theta ^{-B}. \end{align}
In panel (c), we exploit the trend for the peak of turbulent shear stresses that, as suggested by Chen, Hussain & She (Reference Chen, Hussain and She2019), would collapse onto the following correlation:
\begin{align} \nonumber\\[-18pt]\tau _{xy}^{pk} = 1 - C Re_\tau ^{-6/7}. \end{align}
The fitting coefficients
$A$
,
$B$
and
$C$
are provided in the figure caption. As shown, all diagnostics exhibit reasonable agreement with the scaling laws, which gives us confidence that no history effects from the forcing are present.
Evolution of (a)
$E_{VK}$
, (b) skin friction coefficient and (c) peak of shear stresses. In panel (a), dashed grey lines indicate the
$\pm 4\, \%$
error band; in panel (b), dashed grey line indicates (A2) with A = 0.0065 and B = 0.268 as fitting coefficients; in panel(c), dashed grey line indicates (A3) with C = 18.8 as fitting coefficient.
