Dynamic large-eddy simulation of hypersonic transition delay over broadband wall impedance

Abstract

Three-dimensional numerical simulations of hypersonic boundary layer transition delay due to porosity representative of carbon-fibre-reinforced carbon-matrix ceramics (C/C) were carried out on a 7$^{\circ {}}$ half-angle cone for unit Reynolds numbers $Re_m=2.43 \times 10^6$–$6.40\times 10^6\ \text {m}^{-1}$, at the free-stream Mach number $M_\infty =7.4$, for both sharp and 2.5 mm nose tip radii. A broadband time-domain impedance boundary condition was used to model the acoustic effects of the porous surface on the flow field. A quasi-spectral sub-filter-scale dynamic closure was adopted to stabilize the computations upon turbulent breakdown under extreme cooling conditions, with wall-to-adiabatic temperature ratio of $T_{w}/ T_{ad} \simeq 0.08 $, while accurately recovering the growth rates of the unstable modes present in the early transition stages. Good agreement is observed with the reference experimental data, both in terms of the predicted extent of the transition delay and the measured second-mode frequency spectrum. The latter is strongly modulated by the formation of near-wall low-temperature three-dimensional streaks. Pressure disturbances concentrate in corridors of locally thickened boundary layer, with frequencies lower than what predicted by linear theory. Here, trapped wavetrains are formed, which can persist long into the turbulent region. Finally, it is shown that the presence of a porous wall simply shifts the onset of turbulence downstream, without affecting its structure.

1. Introduction

1.1. Background and relevance

Laminar-to-turbulent transition has critical implications for the design, performance and safe operation of hypersonic vehicles. At such high speeds, aerodynamic and thermal loads can increase dramatically upon the onset of turbulence. Reed et al. (1997) studied a low-Earth orbit hypersonic flight vehicle in fully laminar and fully turbulent flow conditions, and concluded that the requirements for the thermal protection systems (TPS) double in the latter case due to a factor of five increase in wall-heat flux caused by laminar-to-turbulent transition. A different study, developed during the National Aerospace Plane project and led by the Defense Science Board committee (NT Force 1992), reached similar conclusions regarding the operational performance of such a vehicle, reporting a threefold decrease in payload-to-weight ratio when comparing fully laminar with fully turbulent conditions. These studies emphasize the importance of laminar-to-turbulent transition delay for sustainability of high-speed flight.

In the low-disturbance environment typical of flight conditions, transition to turbulence over smooth surfaces is initiated primarily by free-stream fluctuations or atmospheric particles (Fedorov 2003; Browne, Al Hasnine & Brehm 2021). This receptivity process is followed by the modal growth of unstable waves, whose behaviour is typically accurately described by linear stability theory (LST) (Morkovin, Reshotko & Herbert 1994; Fedorov 2011). These then trigger nonlinear interactions and, ultimately, breakdown to turbulence (Morkovin 1969). In supersonic and hypersonic flows, the dynamics of the first and second modal instabilities is of great importance for transition. Its contribution is even higher in the case of slender and two-dimensional geometries where effects of other flow instabilities such as Görtler vortices and cross-flow disturbances are minimal.

In canonical incompressible boundary layers, longitudinal disturbances are always more unstable than oblique waves. That is no longer true when approaching boundary-layer-edge Mach numbers of 1.0. Mack's first mode – the compressible equivalent of the Tollmien–Schlichting unstable waves that govern transition in low-speed flow – is most unstable when oriented between 50 and 60$^{\circ }$ with respect to the flow direction. Only one unstable solution is present in the incompressible regime, but when a compressible flow is considered, a numerable set of additional unstable modes are possible. As the Mach number increases, the relative importance of the first of these additional transitional modes – i.e. Mack's second mode (Mack 1984) – increases and becomes dominant over the first. For adiabatic wall conditions, this happens around Mach numbers of $M = 4$ (Mack 1969). The presence of cooled walls destabilizes the second mode while stabilizing the first mode. At such conditions, two-dimensional disturbances govern the boundary layer stability dynamics once again. Mack's studies on the evolution of these and higher modes (Mack 1990) in two-dimensional boundary layers underscores the importance of second-mode attenuation in any attempt to delay transition in high-speed flows.

1.2. Transition in high-speed boundary layers

When the perturbation amplitudes are sufficiently large to initiate mode-to-mode interactions, nonlinear breakdown to turbulence can occur. The interaction between two oblique instability waves with equal but opposite wave angles is referred to as oblique breakdown, whereas resonance occurs when a two-dimensional or axisymmetric wave reaches an amplitude level sufficient to excite a secondary oblique instability at the same frequency (fundamental resonance) or half of the main mode frequency (subharmonic resonance).

The oblique breakdown path was first seen as a viable path to reach transition by Bestek, Thumm & Fasel (1993) in direct numerical simulations of a flat-plate boundary layer transition at Mach 1.6. The aforementioned path was also highlighted by the subsequent study of Chang & Malik (1994), who performed a nonlinear parabolic stability equation assessment of flows at low supersonic Mach numbers. In a temporal evolution simulation of a Mach 4.5 boundary layer, Adams & Kleiser (1996) reported that subharmonic resonance mechanisms had the strongest signature before breakdown. Since the simulation was performed with a temporally evolving parallel flow set-up, the behaviour of a spatially developing boundary layer was not accurately represented. Franko & Lele (2013) performed a direct numerical simulation (DNS) of a spatially developing flat-plate boundary layer and concluded that fundamental resonance was much stronger than subharmonic resonance. Additionally, when compared against the first-mode oblique breakdown transitional path to turbulent flow for the set-up studied in Franko & Lele (2013), both of the previously mentioned mechanisms develop slower and need an increased streamwise extent to induce a fully developed turbulent flow.

A formidable amount of high-speed flow transition research has focused on slender cones (Schneider 2004) due to their geometric simplicity and low cross-sectional blockage in ground testing. Early transition experiments (Demetriades 1974; Kendall 1975) reported the second mode as the most important unstable wave at $M_{\infty }$ numbers of 8 and 8.5 over cones with half-angles equal to $6^\circ$ and $4^\circ$, respectively. Stetson & Kimmel (1993) later performed experiments at Mach 8 over a 7$^{\circ }$ half-angle cone and reported the existence of higher harmonics with no evidence of subharmonic instability. The same flow conditions were used in a numerical study performed by Pruett & Chang (1995), who investigated different possible transition routes, concluding that oblique breakdown was more likely to occur than subharmonic resonance.

More recently, Sivasubramanian & Fasel (2014) used a broadband pulse to study a natural transition scenario in a hypersonic conical boundary layer with $T_w \approx T_{{ad}}$ and concluded that both fundamental resonance and second-mode oblique breakdown are equally viable transition paths at such near adiabatic conditions. A controlled-disturbance simulation of fundamental resonance under the same flow conditions was performed by Sivasubramanian & Fasel (2015), who carried out a detailed analysis of the nonlinear interactions leading up to turbulent breakdown. Some discrepancies between controlled-disturbance transition simulations and experiments were observed. Later, Hader & Fasel (2018) proposed a pseudo-random forcing that would remove possible unwanted biases towards a particular breakdown mechanism. The proposed strategy confirmed the relevance of fundamental resonance as a nonlinear breakdown path. This conclusion was supported by Chynoweth et al.'s (2014) transition study on a flared cone. The pseudo-random perturbation technique also led to an improved matching against the mean wall-heat flux measured experimentally in comparison with what was predicted by introducing controlled disturbances. A pseudo-random perturbation technique inspired by such previous work is also used in the current study, albeit in a grid-independent fashion, to mimic a natural transition scenario.

1.3. Transition delay over porous walls

Malmuth et al. (1998) proved through LST the capability of ultrasonically absorptive coatings (UACs) to attenuate the second mode. Experiments performed in the CalTech T-5 hypersonic wind tunnel (Fedorov et al. 2001; Rasheed et al. 2002) provided the first validation of this principle by using a cone with a half-angle of $5^\circ$ with two different surfaces: one smooth and impermeable, and the other perforated with regularly spaced micro-holes. They reported that the porous surface was capable of doubling the transitional Reynolds number.

The development of UACs continued with the goal of combining the acoustic absorption and TPS characteristics into one solution. Knowing that an irregular porous structure is typical of TPS materials used for hypersonic vehicles, Fedorov et al. (2003, 2006) conducted experiments in the T-326 Mach 6 wind tunnel at the Institute of Theoretical and Applied Mechanics (ITAM) in Novosibirsk first over a porous surface with regular microstructure, and then over a felt-metal surface. These experiments were the first to produce data capable of quantitatively demonstrating the growth rate attenuation of the second-mode waves via UACs, also confirmed by stability analysis. More recently, an error in the stability analysis performed by Fedorov et al. (2003) was corrected in Tritarelli, Lele & Fedorov (2015), demonstrating the presence of a shift of second-mode instabilities to lower frequencies over felt-metal porous structures.

Chokani et al. (2005) performed a bispectral analysis on this same data with the goal of identifying nonlinear mechanisms that could be triggered by the use of porous coatings; this was motivated by the observation of a weak enhancement of first-mode waves. The experiments revealed the occurrence of nonlinear phase locking involving first-mode waves exclusively present over porous walls and not over an impermeable surface. This effect was also shown to be small, not overcoming the overall stabilizing role of the porous coating. More recently, Lukashevich, Morozov & Shiplyuk (2016) performed new experiments in ITAM's tunnel, where it was shown that placing the porous insert in the region where the second mode is stable leads to an increase in the disturbance signal amplitude; the opposite is observed when the insert is placed in regions of second-mode growth. Angles of attack below $1^{\circ }$ were studied experimentally in the same tunnel by Morozov et al. (2018), who reported that the porous coatings were able to attenuate second-mode disturbances on either the windward or leeward side of the cone.

Motivated by the need for any realizable boundary layer control technology to be symbiotic with TPS, Wagner et al. (2013) pioneered the use of porous carbon–carbon (C/C) to control second-mode waves. This material represents the intermediate state in the fabrication process of the C/C silicon carbide already employed on hypersonic vehicles (Turner et al. 2006; Weihs, Longo & Turner 2008). The C/C material exists after the pyrolysis of the green body composed of carbon fibre reinforced plastic and before the liquid silicon infiltration that leads to the final product (Dittert & Kütemeyer 2017; Wagner et al. 2019). Experiments at approximately Mach 7.5 have been conducted in the DLR High Enthalpy Shock Tunnel Göttingen (HEG) demonstrating an increase of the laminar portion of the boundary layer as a result of second-mode attenuation over the porous surface. For these flow conditions, Wartemann et al. (2015) has performed an analysis based on the parabolized stability equations, confirming numerically the effectiveness of porous surfaces based on C/C in achieving transition delay. More recently, Running et al. (2023) performed transition experiments at Mach 6 over a flat plate at a negative angle of attack using a silicon-carbide foam and reported that the porous material had a strong ability to attenuate the instability mechanisms inside the hypersonic boundary layer.

Previous pore-resolving numerical studies of attenuating, cancelling or reinforcing second-mode instabilities in high-speed flow have focused on geometrically regular porosity over temporally developing boundary layers (Brès, Colonius & Fedorov 2008; De Tullio & Sandham 2010; Brès et al. 2013). Two-dimensional simulations investigating the influence of regular and irregular porous surfaces on the disturbances of spatially developing boundary layers over a flat plate (Egorov, Fedorov & Soudakov 2008; Wang & Zhong 2010, 2011) and over an axisymmetric cone (Lukashevich et al. 2012) were also performed based on the analytical relation developed by Fedorov et al. (2001), which couples the wall-normal velocity disturbance ($u'_n$) with the pressure disturbance ($p'$) at the boundary

(1.1)\begin{equation} u'_n(x,t) = p'(x,t) \mathrm{Real}(A_n) - \frac{1}{\omega} \frac{\partial p' (x,t)}{\partial t} \mathrm{Imag}(A_n), \end{equation}

where $A_n$ is the wall-normal admittance, being the wall admittance defined as the inverse of the wall-normal impedance

(1.2)\begin{equation} Z_{n}(\omega) = \frac{1}{\rho_0 a_0} \frac{\hat p(\omega)}{\hat u_n(\omega)} = A_n^{{-}1}, \end{equation}

where $\hat p(\omega )$ and $\hat u_n(\omega )$ are the Fourier transform of the pressure and wall-normal velocity fluctuations at the surface, $\omega$ is the angular frequency and $\rho _0 a_0$ is the base impedance. Equation (1.1) is effectively an impedance boundary condition capable of imposing a complex wall impedance only at a fixed frequency.

While previous numerical simulations on the effects of porous walls on boundary layers were limited to the investigation of the stability and receptivity frequency by frequency, the current work models the porous walls effects over the whole frequency spectrum, enabling simulations of the three-dimensional boundary layer transition process, from second-mode instability growth to three-dimensional turbulence breakdown. The present work is a significant extension of Sousa et al. (2023), who limited their investigation to axisymmetric configurations.

A time-domain impedance boundary condition (TDIBC) capable of imposing broadband acoustic effects from a porous surface in a high-fidelity calculation was first demonstrated by Scalo, Bodart & Lele (2015) in low-speed ($M_b<0.5$) compressible turbulent channel flow, relying on the causality principles outlined by Fung & Ju (2001). Initially, the technique was used with a single pole, concentrating the frequency response around a resonant frequency akin to Helmholtz resonance. It was shown that, for sufficiently high porosities, large-scale rollers with axis of rotation oriented in the spanwise direction were induced in the near-wall region, significantly altering the Reynolds stress profiles. More recently, a similar flow set-up was used by Chen & Scalo (2021) to study the porosity effects on supersonic and hypersonic turbulent channel flows. It was revealed that sufficiently high porosity can trigger the presence streamwise-travelling waves in the near-wall region.

Further development of the TDIBC technique was done by Lin, Scalo & Hesselink (2017), who implemented multipole reconstruction of the frequency response of piezoelectric energy extraction of a thermoacoustic engine, and by Douasbin et al. (2018), who devised an iterative fitting technique for the multipoles necessary to model the reflection coefficient of a laminar combustion chamber. In this context, Sousa et al. (2019) used a similar iterative fitting technique for the impedance frequency response of a porous wall and showed the attenuation effect on a broadband second-mode disturbance advecting over an axisymmetric conical boundary layer. Ultimately, in the current work, also done by Sousa et al. (2023), additional improvement on the fitting of the complex impedance response was done by using the vector fitting procedure described by Gustavsen & Semlyen (1999).

1.4. Manuscript outline

The current work aims at investigating the influence of microscale distributed wall porosity on the full transition path to turbulence by replicating the experiments performed in the HEG wind tunnel by Wagner et al. (2013) and Wagner (2014) in a cost-effective manner through the utilization of a novel dynamic large-eddy simulation (LES) technique outlined below. The present manuscript is based on the most recently available experimental dataset, which can be found in Wagner et al. (2019): a technical report AFRL-AFOSR-UK-TR-2020-0025 published by the Defense Technical Information Center (DTIC).

To the authors’ knowledge this is the first set of three-dimensional hypersonic boundary layer transition simulations performed at these flow conditions. The short duration of the experimental runs does not allow for the surface to heat up, leading to wall temperature adiabatic temperature ratios of approximately 0.1. The sharp gradients within the laminar boundary layer become more pronounced when turbulence sets in. These drive strong shock-like density and pressure fluctuations that need to be explicitly numerically dissipated, especially when using spectral or quasi-spectral numerics, to ensure numerical stability of the calculations at any level of resolution. This should be done without compromising solution accuracy in areas where the solution is well resolved, e.g. second-mode growth in the presence of a laminar base state. In this manuscript, we employed the quasi-spectral viscosity (QSV) method as presented in Sousa & Scalo (2022a) to address this requirement and reduce the resolution demands in the fully turbulent region. The QSV-LES method offers the benefits of the rapid convergence of spectral numerics in the transitional region, while removing the need for DNS-like resolution in the turbulent region. More broadly, the QSV method offers a dynamic sub-filter-scale (SFS) closure addressing both hydrodynamic SFS turbulence modelling and shock capturing needs under a unified mathematical framework for a high-order finite-difference discretization of the filtered compressible Navier–Stokes equations.

For the present dynamic LES transitional simulations to faithfully capture the flow dynamics observed in the experiments, it is imperative that the added numerical dissipation is capable of deactivating itself in laminar flow regions. This capability was demonstrated in prior publications by the same authors (Sousa et al. 2023), where axisymmetric simulations carried out with active QSV-LES closure at the same HEG wind tunnel conditions (Wagner et al. 2013; Wagner 2014) predicted second-mode wave frequencies and growth rates matching the linear stability predictions and companion DNS calculations. The dynamic nature of the QSV closure is further demonstrated here (see discussion regarding figures 9 and 10).

The computation methodology is first presented in § 2. This includes the computational-domain discretization, the governing equations solved, the dynamic LES scheme and the boundary conditions used in the simulations. Special attention is given to the TDIBC used to introduce the acoustic effects of the C/C surfaces without the need to resolve its intricate geometric structure.

In § 3, three-dimensional simulations of the full transition path to turbulence for different nose tip bluntness and Reynolds numbers are performed, with results compared against experimental data. Ultimately, the physical aspects of the flow are analysed and the effects of wall porosity on the transitional and turbulent regions are assessed.

The present effort is a significant extension of the previous work published by the same authors (Sousa et al. 2019, 2023) which focused only on axisymmetric simulations of linearly perturbed hypersonic boundary layers over sharp and blunt cones. Improvements in the dynamic LES modelling allowed three-dimensional high-fidelity simulations of the complete transition from laminar to turbulent flow and supported comparisons against experimental measurements.

2. Problem formulation

2.1. Computational set-up

The numerical study performed in this work is based on experiments conducted by Wagner et al. (2013) in the HEG. A complete description of such facility can be found in Hannemann et al. (2018). Three-dimensional (§ 3) dynamic LESs are carried out for the hypersonic boundary layer flow on a sharp cone at Reynolds numbers $Re_m = 2.43 \cdot 10^6\ \text {m}^{-1}$ and $4.06 \cdot 10^6\ \text{m}^{-1}$; and on a 2.5 mm blunt tip cone at $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$ and $6.40 \cdot 10^6\ \text {m}^{-1}$. Due to the short duration of the test run, the surface of the model stays at room temperature ($T_w=300$ K), effectively yielding an isothermal boundary condition. It is important to note that, although peak temperatures in the hypersonic boundary layer at the aforementioned experimental conditions can reach approximately 800 K, the current numerical studies have assumed a caloric perfect gas.

A list of all the available experimental flow conditions is reported in table 1. The reader is invited take note of the colours associated with each free-stream Reynolds number, which are used (at times with different degrees of shading or intensity) consistently throughout the manuscript to plot quantities at the corresponding flow conditions. The only exception are the flooded contour plots, where perceptually uniform sequential colour maps were used to improve readability.

Table 1. Free-stream flow parameters from Wagner et al. (2013) with corresponding colour scheme and symbols used throughout the manuscript.

In the case of the sharp cone, the initial flow conditions and inflow parameters are informed by combining the Taylor & Maccoll (1933) inviscid solution with a viscous solution for the boundary layer. The latter is derived by applying the Mangler transformation for bodies of rotation (Mangler 1948) to the compressible boundary layer similarity solution for a perfect gas over a flat plate (Illingworth 1949), as done by Lees (1956) for the flow over a cone. On the other hand, the blunt-nose simulations required an external precursor calculation performed by the DLR FLOWer code – a three-dimensional parallel hybrid multi-grid code, validated for hypersonic flows (Kroll & Fassbender 2006).

To avoid the computational costs associated with the very thin boundary layer near the leading edge of the cone, the inlet of the computational domain in this study is located at $x_i=0.045$ m downstream from the sharp cone tip, while the outlet is located at $x=L_x$ m (see table 2) further downstream, as shown in figure 1. The streamwise extent of the computational domain $L_x$ is chosen based on previous experimental results as the smallest distance needed to be able to capture the full transition path to turbulence for both impermeable and porous surface cases. The grid arrangement In the wall-normal direction follows the strategy used by Sousa et al. (2023) where the computational domain is placed entirely below the shock and wall-normal grid stretching is performed as a function of the streamwise distance from the nose tip with the intent of accurately resolving a growing boundary layer.

Table 2. Streamwise extent $(L_x)$ of the computational domain (figure 1) for the different flow conditions and cone geometries considered.

Figure 1. Computational set-up comprises a $7^\circ$ half-angle cone with $R_{tip} = 0$ (sharp) and $2.5$ mm (blunt) tip radii. The $A^{-}$ represents the incident wave, which drives the convolution process at the impedance boundary (§ 2.3), resulting in the output wave $A^{+}$. The impedance boundary starts at $X = 0.182$ m matching the experiments conducted by Wagner et al. (2013).

At the inlet and top boundaries, a Dirichlet boundary condition is used to input the flow velocity calculated by the analytical procedure or by a precursor simulation. Outlet conditions are homogeneous Neumann for all flow quantities. Periodic boundary conditions are imposed at the azimuthal extremes of the domain. Lastly, at the wall, Dirichlet conditions are used to impose no slip, no penetration and, when active, the TDIBCs used to model the porous surface. This approach allows us to model the acoustic response of a C/C-composed surface without the need to resolve its intricate geometrical structure and is a way to retain high resolution on the flow side by removing the grid resolution requirements that would be needed to capture the acoustic wave propagation in the pores.

Sponge layers are used at all non-solid boundaries. The inflow and outflow sponges extend into the domain by 0.03$L_x$, while the top boundary sponge is 5 % of the local domain height. The inlet sponge is used to weakly impose inflow boundary conditions, while the outlet sponge prevents the formation of upstream-travelling disturbances propagating in the subsonic portion of the boundary, while suppressing flow fluctuations convecting downstream before they interact with the outflow computational boundary.

2.2. Filtered governing equations in curvilinear coordinates

High-order structured compact-finite-difference simulations are carried out with the objective of capturing the perturbation evolution with minimum numerical dissipation and high resolving power for a given amount of points per wavelength. Additionally, a dynamic SFS model via the QSV closure (Sousa & Scalo 2022a), capable of performing both shock capturing and turbulence modelling, is leveraged to guarantee numerical stability upon three-dimensional turbulent breakdown under high cooling ratios ($T_w/T_r \simeq 0.1$) while preserving the second-mode wave dynamics, as demonstrated in Sousa et al. (2023). This results in the ability of the model to achieve convergence of the transition location, defined as the point where the mean heat flux at the surface departs from the laminar correlation (discussed in § 3.2), while also exhibiting the expected grid-dependent behaviour of a SFS closure in the turbulent region.

The filtered compressible Navier–Stokes equations in generalized curvilinear coordinates are presented hereafter following previous work by Jordan (1999) and Nagarajan, Lele & Ferziger (2007), who developed incompressible and compressible SFS closures, respectively. We assume the existence of a known, invertible mapping between $\boldsymbol {y}$, the physical Cartesian reference frame, and $\boldsymbol {x}$, the contravariant curvilinear coordinate system

(2.1)\begin{gather} x^i = x^i(y^1,y^2, y^3), \end{gather}

(2.2)\begin{gather}y^i = y^i(x^1,x^2, x^3), \end{gather}

where $x^i$ and $y^i$ are the $i$th coordinate of each respective system of reference. We then define a Favre filtering operator

(2.3)\begin{equation} \check f = \frac{\overline{J \rho f}}{\overline{J \rho}}, \end{equation}

which accounts for the Jacobian of the transformation, $J$, defined as the determinant of the Jacobi matrix ($J_{ij} = \partial y^i/\partial x^j$). The filtered governing equations can then be defined as

(2.4)\begin{gather} \frac{\partial \overline{J \rho}}{\partial t} + \frac{\partial }{\partial x^j}( \overline{J \rho} \check v^j)= 0, \end{gather}

(2.5)\begin{gather}\frac{\partial \overline{J \rho}\check v^i}{\partial t} + \frac{\partial}{\partial x^j}(\overline{J \rho} \check v^i \check v^j + \overline{J p} g^{ij} - J \check \sigma^{ij} + \overline{J \rho} \tau^{ij})={-}\varGamma^i_{qj}(\overline{J \rho} \check v^q \check v^j + \overline{J p} g^{qj} - J \check \sigma^{qj} + \overline{J \rho} \tau^{qj}), \end{gather}

(2.6)\begin{gather}\frac{\partial \overline{J E}}{\partial t} + \frac{\partial}{\partial x^j}(\overline{J(E + p)}\check v^j + J \check Q^{j})= \frac{\partial}{\partial x^k}( J \check \sigma^{ij}g_{ik}\check v^k )- \frac{\partial \overline{J \rho} C_p q^{j} }{\partial x_j}, \end{gather}

and the relevant SFS terms are

(2.7a,b)

Here, $C_p$ denotes the specific heat capacity at constant pressure and $v^i$ are the velocity components in contravariant curvilinear coordinate system. Sub-filter contributions resulting from the nonlinearities involving the dependency of molecular viscosity or conductivity on temperature have been neglected following Vreman, Geurts & Kuerten (1995), who showed those are negligible compared with other SFS terms. Following Nagarajan, Lele & Ferziger (2003), the SFS kinetic energy advection, , and the SFS turbulent heat dissipation, $\epsilon = \partial (\overline { \sigma ^{ij} v^i}) / \partial x^j - \partial \check (\sigma ^{ij} \check v^i) / \partial x^j$, terms are also neglected.

The tensors responsible for mapping a curvilinear physical space into a Cartesian reference space are the covariant and contravariant metric tensors, respectively,

(2.8a,b)\begin{equation} g_{ij} = \frac{\partial y^i \partial y^j}{\partial x^k \partial x^k} \quad \text{and} \quad g^{ij} = \frac{\partial x^k \partial x^k}{\partial y^i \partial y^j}, \end{equation}

as well as the Christoffel symbol of the second kind

(2.9)\begin{equation} \varGamma^i_{qj} = \frac{\partial x^i}{\partial y^l} \frac{\partial^2 y^l}{\partial x^q \partial x^j}. \end{equation}

In the derivation of these equations, the metric tensors and Christoffel symbols are assumed to be varying slowly over the spatial support of the filter kernel, therefore leading to additional negligible sub-filter flux terms.

In the curvilinear frame of reference the total energy, the viscous stress tensor and the heat flux vector are described by slightly modified relations described below

(2.10)\begin{gather} \frac{\overline{J p}}{\gamma - 1} = \overline{J E} - \frac{1}{2} \overline{J \rho} g_{ij} \check v^i \check v^j - \frac{1}{2} \overline{J \rho} g_{ij} \tau^{ij}, \end{gather}

(2.11)\begin{gather}\check \sigma^{ij} = \mu \left( g^{jk} \frac{\partial \check v^i}{\partial x^k} + g^{ik} \frac{\partial \check v^j}{\partial x^k} - \frac{2}{3} g^{ij}\frac{\partial \check v^k}{\partial x^k} \right), \end{gather}

(2.12)\begin{gather}\check Q^j ={-}k g^{ij} \frac{\partial \check T}{\partial x^i}. \end{gather}

Ultimately, the QSV closures for $\tau ^{mn}$ and $q^{n}$ when applied to generalized curvilinear coordinates are

(2.13) $$\begin{gather} \tau^{mn} ={-} C_{\tau^{mn}} \sum_{k=1}^3 \frac{1}{2} \left( g^{nk} \mathcal{D}^{mk}\frac{\partial \ddot v^m}{\partial x^k} + g^{mk} \mathcal{D}^{nk}\frac{\partial \ddot v^n}{\partial x^k} \right) , \end{gather}$$

(2.14)$$\begin{gather}q^{n} ={-} C_q \sum_{k=1}^3 g^{nk} \mathcal{D}^{kk}\frac{\partial \ddot T}{\partial x^k}, \end{gather}$$

where there is no summation implied over repeated $m$ and $n$ indices. Operations over the $k$ index are explicitly regulated by the summation symbol. The above relation is an improved notation over Sousa & Scalo (2022a). The double-dot superscript indicates spectrally modulated quantities, defined in curvilinear coordinates as

(2.15)\begin{gather} \frac{\partial \ddot v^m}{\partial x^k} = \frac{\partial \check v^m}{\partial x^k}* (1 - \tilde{G}_{{qsv}}), \end{gather}

(2.16)\begin{gather}\frac{\partial \ddot T}{\partial x^k} = \frac{\partial \check T}{\partial x^k}* (1 - \tilde{G}_{{qsv}}), \end{gather}

where $*$ is the convolution operator and,

(2.17)\begin{equation} \hat{\tilde{G}}_{{qsv}} = w[\beta \hat{\tilde{G}}_{{Pade}}(\alpha) + (1-\beta) \hat{\tilde{G}}_{{Fejer}}], \end{equation}

is the QSV's filter modulation transfer function, comprising a weighted average between a Padé and a Fejér filter. Additionally, $\mathcal {D}^{mn} = \upsilon ^m(\check v^n) \ell ^n$ is the dissipation magnitude tensor, comprising a length scale related to the physical grid spacing ($\ell ^n$) and the sub-filter velocity scale

(2.18)\begin{equation} \upsilon^m(\check v^n) = \sqrt{ \frac{2 E^m_{k_c}(\check v^n )}{\bar{\varDelta}^m}}, \end{equation}

where $\bar {\varDelta }^m$ is related to the computational grid spacing. Again, there is no summation implied over the $m$ and $n$ indices. The estimation of the energy content near the grid cutoff is performed via the operation

(2.19)\begin{equation} E^m_{k_c} (\check v^n) = \left[ (1 - \widetilde{G^m}_{{Pade}}(\alpha))* \frac{(\check v^n)^2}{2} \right]\frac{A^m(\alpha_{k_c})}{A^m(\alpha)}. \end{equation}

Equations (2.17) and (2.19) reflect typographical corrections pointed out by the corrigendum (Sousa & Scalo 2022b) of the original manuscript (Sousa & Scalo 2022a) outlining the QSV methodology. The full description of the spectral modulation operation done via $\tilde {G}_{{qsv}}$ can be found in Sousa & Scalo (2022a). The default values used in the aforementioned publication are also used in the current work, i.e. the coefficient $C_{\tau ^{mn}}$ is equal to 1 when $m=n$ and set to 0.6 when $m \neq n$, $C_{q} = 0.8$, $w = 0.85$, $\alpha =0.45$ and $\beta =0.8$.

The proposed closures for the SFS terms are not proper tensors representing stresses in the physical space, as the amplitudes of the different components that collectively define $\tau ^{mn}$ and $q^n$ are dependent on semi-local operations along a certain direction and therefore are not invariant with respect to a rotation of coordinates. However, implementing the SFS closure directly in the computational space makes the QSV methodology independent to the specific physical geometry, making it robust to geometrical complexities in the physical space. As a final remark, the results gathered in the current manuscript as well as in Sousa & Scalo (2022a) demonstrate QSV's capabilities on both Cartesian and curvilinear grids. Nonetheless, the consequences of the proposed closure not transforming like a proper tensor in physical space will be studied in future research.

The filtered compressible Navier–Stokes equations are discretized in space using a sixth-order accurate compact-finite-difference method with a staggered grid arrangement (Lele 1992) and in time using a third-order accurate explicit Runge–Kutta method. The numerical strategy used yields spectral-like resolution properties, being able to solve the perturbation evolution inside the boundary layer with minimum numerical dissipation and with high accuracy for a given amount of points per wavelength (Nagarajan et al. 2003). Additionally, the QSV closure has a dynamic property, being only active near areas where nonlinear effects are relevant and preserving the structure of ultrasonic transitional waves.

Moreover, in order for the present transitional simulations to faithfully capture the flow dynamics observed in the experiments they seek to emulate, it is imperative that the added numerical dissipation model is capable of adjusting its intensity. On the same note, the QSV-LES technique is able to dynamically modulate its magnitude and that facilitates the numerical computation of the transitional dynamics of the second-mode growth with minimal interference. This property was demonstrated by Sousa et al. (2023). In that manuscript, axisymmetric simulations of a hypersonic laminar boundary layer over a conical model were carried out. These simulations involved applying a single pulse perturbation under the experimental flow conditions reported in table 1. Even though they employed an active QSV model, the obtained second-mode growth rates closely aligned with the expected values predicted by linear stability theory.

The generic curvilinear coordinate system transformation applied in the scope of this work maps the physical-domain Cartesian coordinates $(y^1 = X, y^2 = Y, y^3 = Z)$ into (and from) the cylindrical coordinate system $(x^1 = X, x^2 = r, x^3 = \theta )$, where the equations are ultimately solved (see figure 1). Additionally, the boundary layer coordinates $x$ and $y$, representing streamwise and wall-normal directions, respectively, are also used, this being a natural way of displaying results of the flow over a conical geometry.

2.3. Time-domain impedance boundary conditions

In this work, the acoustic effect of the porous walls is modelled via the imposition of a wall-normal broadband complex impedance boundary condition, $Z_{n}(\omega )$, defined as the ratio of complex amplitudes of pressure and wall-normal velocities

(2.20)\begin{equation} Z_{n}(\omega) = \frac{1}{\rho_0 a_0} \frac{\hat p(\omega)}{\hat u_n(\omega)}, \end{equation}

normalized by the base impedance, $\rho _0 a_0$, also evaluated at the wall. Its inverse $Y_n(\omega ) = Z_{n}(\omega )^{-1}$ is referred to as the admittance. Impedance in the tangential directions can also be defined, relating instantaneous shear with the slip velocity at the interface.

Equation (2.20) assumes linearity of the acoustics inside the pores, which is justified by the following reasons: first, the thermoviscous dissipation experienced by the waves travelling inside the pores is sufficiently high to ensure rapid wave amplitude attenuation as supported by the conclusions drawn in Sousa et al. (2023), where a relatively thin surface layer of the C/C material is found to contribute to most of the absorption; second, the low density and pressure environment considered inhibits wave steepening (Thirani, Gupta & Scalo 2020), making acoustic nonlinearities in the pores unlikely to be important.

Acoustic reflections off an impedance boundary can be predicted with only wall-normal impedance information, relying on the assumption of locally acting surface (Kinsler et al. 1999). Such a hypothesis is valid in the current case because the C/C material in question is highly anisotropic (figure 2), with pores primarily aligned in the wall-normal direction and with small orifice sizes compared with the boundary layer thickness. Moreover, a parametric study allowing various degrees of wall-slip velocities demonstrated that the inviscid nature of the second-mode instability waves makes tangential impedance effects not as relevant as they would be in otherwise viscous dominated instabilities (Fedorov et al. 2003; Fedorov & Malmuth 2008). Therefore, we chose to model the acoustic effects of a porous C/C surface on the flow as a linear normal impedance boundary condition.

Figure 2. Reflected-light microscopy generated images (a) and a highlight of the intrinsic pores (b) of the C/C material analysed in this study.

Imposing directly the normal impedance at the boundary is impractical since its magnitude approaches infinity as the frequency goes to zero. To overcome this, different formulations using either the reflection coefficient ($\hat R$) or the wall softness ($\hat S$), defined as

(2.21a,b)\begin{equation} \hat{R}(\omega) = \frac{1-Z_{n}(\omega)}{1+Z_{n}(\omega)} ; \quad \hat{S}(\omega) = \hat{R}(\omega) + 1, \end{equation}

were developed, because these quantities are bounded in amplitude across the frequency spectrum when a passive boundary is considered. A more intuitive physical understanding of these complex-valued quantities can be achieved by formulating them as a function of the incident, $\hat {A}^-(\omega )$, and reflected, $\hat {A}^+(\omega )$, waves from a given surface (see figure 1) in the frequency domain via

(2.22a,b)\begin{equation} \hat{A}^-(\omega)\hat{R}(\omega) = \hat{A}^+(\omega) \quad \mathrm{and} \quad \hat{A}^-(\omega)\hat{S}(\omega) = \hat{A}^+(\omega) + \hat{A}^-(\omega) , \end{equation}

with all quantities evaluated at the impedance surface. Equation (2.22a,b) shows that the reflection coefficient can be interpreted as a transfer function ($\hat {R}$) seeing the incoming wave ($\hat {A}^-$) as the input and the reflected wave ($\hat {A}^+$) as the output. Similarly, the wall softness coefficient can be understood as a transfer function that returns the combined response of both incoming and outgoing waves at the surface given a certain perturbation input. The magnitude of either the complex softness or reflection coefficient is directly related to the amount of energy absorbed by an impedance surface and its argument, or angle, to the phase shift imparted by it.

The absorption coefficient, related to the acoustic power loss caused by thermoviscous dissipation inside the pores, is defined as

(2.23)\begin{equation} \beta(\omega) = 1 - | \hat R(\omega) |^2.\end{equation}

To apply a complex impedance in a numerical simulation, one needs to transform it to the time domain. The TDIBC is implemented in the current numerical scheme as a function of the wall softness ($\hat S$) based on the formalism of Fung & Ju (2004) and the numerical integration strategy in a compressible Navier–Stokes solver by Scalo et al. (2015), and subsequent work by Lin, Scalo & Hesselink (2016) and Douasbin et al. (2018). The complex wall softness $\hat S(\omega )$ is fitted with a set of complex pairs of poles $(p_k,p_k^*)$ and residues $(\mu _k,\mu _k^*)$

(2.24)\begin{equation} \hat{S}(\omega) = \frac{2}{1 + Z_{n}(\omega)} \approx \sum_{k=1}^{n_o} \left[ \frac{\mu_k}{s - p_k} + \frac{ \mu_k^{*}}{s - p_k^{*}} \right], \end{equation}

in such a way that their frequency-domain response is constrained so that it satisfies the causality, passivity and reality constraints (Rienstra 2006). In (2.24), $n_o$ is the number of oscillators, each described by a conjugate pair of residues and poles, with $s=j\omega$ being the complex Laplace coordinate. The fit in the frequency domain is performed following the vector fitting procedure described by Gustavsen & Semlyen (1999), where a least square approximation is performed iteratively to relocate the poles and residues and accurately represent the complex response of the function to be fitted in a given interval. Figure 3 shows a representative result of this fitting procedure using the C/C material acoustic properties as an example. A more thorough discussion on the impedance eduction method adopted in the current study can be found in Wagner et al. (2018a) and Sousa et al. (2023).

Figure 3. Real (a) and imaginary (b) components of the complex wall softness, $\hat S (\omega )$, representative of the C/C surface at wall pressure conditions corresponding to $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$ ($\circ$) and the multipole reconstruction of the complex signal following the vector fitting procedure (——–). The details of the model parameters and impedance eduction technique driven by ultrasonic bench test data can be found in Sousa et al. (2023).

The implementation of TDIBC in a Navier–Stokes solver uses the relation between the complex wall softness coefficient, $\hat {S}(\omega )$, with the incident, $\hat {A}^-(\omega )$, and reflected, $\hat {A}^+(\omega )$, waves as described in (2.21a,b). A (trivial) inverse Fourier transform of (2.24) is then taken

(2.25)\begin{equation} S(t) \approx \mathcal{H}(t) \sum_{k=1}^{n_o} [ \mu_k \,{\rm e}^{p_k t} + \mu_k^{*} \,{\rm e}^{p_k^{*} t}], \end{equation}

where $\mathcal {H}(t)$ is the Heaviside function, allowing us to construct the following convolution integral:

(2.26)\begin{equation} {{A}^+}(t) ={-}{{A}^-}(t) + \int_0^{\infty} S(\tau){{A}^-}(t - \tau) \, {\rm d} \tau, \end{equation}

carried out in the positive semi-infinite temporal axis $\tau = [0,\infty ]$ (Fung & Ju 2004), hence respecting causality since it only involves values of ${A}^-$ prior to the current time $t$. The time-domain definition of the incident and reflected waves is

(2.27)\begin{equation} \left.\begin{array}{@{}l@{}} A^-(t) = v'(t) + p'(t) / \rho_0 a_0 \\[3pt] A^+(t) = v'(t) - p'(t) / \rho_0 a_0 \end{array}\right\}, \end{equation}

where $v'$ and $p'$ are the fluctuating values of wall-normal velocity evaluated at the impedance boundary.

To numerically carry out the convolution integral shown in (2.26) and couple it with a Navier–Stokes solver, it needs to be evaluated at $t + \Delta t$, where $\Delta t$ is a finite time interval used to advance the numerical solution in time, as

(2.28)\begin{equation} {{A}^+}(t+\Delta t) ={-}{{A}^-}(t+\Delta t) + \sum_{k=1}^{n_o} [ s_k (t + \Delta t) + s_k^* (t + \Delta t)], \end{equation}

where

(2.29)\begin{equation} s_k (t + \Delta t) = \int_0^{\infty} \mu_k \,{\rm e}^{p_k t} {{A}^-}(t + \Delta t- \tau) \, {\rm d} \tau . \end{equation}

Using the internal addition property of the integral in (2.29), Scalo et al. (2015) showed that

(2.30)\begin{equation} s_k (t + \Delta t) = z_k s_k (t) + \int_0^{\Delta t} \mu_k \,{\rm e}^{p_k t} {{A}^-}(t + \Delta t- \tau) \, {\rm d} \tau, \end{equation}

where $z_k = {\rm e}^{p_k t}$. Following Fung & Ju (2004), the integral in (2.30) is evaluated using a trapezoidal quadrature rule, which is second-order accurate in time. The resulting relation for the contribution of each pole and residue pair can be written as

(2.31)\begin{equation} s_k (t + \Delta t) = z_k s_k (t) + \mu_k \Delta t [ w_{k0} {{A}^-}(t + \Delta t) + w_{k1} {{A}^-}(t) ], \end{equation}

where

(2.32a,b)\begin{equation} w_{k0} = \frac{z_k - 1}{p_k^2\Delta t^2} - \frac{1}{p_k \Delta t} \quad \mathrm{and} \quad w_{k1} ={-}\frac{z_k - 1}{p_k^2\Delta t^2} + \frac{z_k}{p_k \Delta t}. \end{equation}

After the numerical evaluation of the convolution integral and after summing the contribution of each pole and residue pair, the outgoing wave $A^+$ at time $t + \Delta t$ can be determined and, ultimately, the wall-normal velocity can then be recovered and imposed at each time step as

(2.33)\begin{equation} v'(t+\Delta t) = \tfrac{1}{2}[ A^-(t+\Delta t) + A^+(t+\Delta t)]. \end{equation}

As a final remark, a higher-order TDIBC implementation based on auxiliary differential equations was developed by Chen & Scalo (2021) and applied to the investigation of wall porosity effects on supersonic and hypersonic turbulent channels. In practice, no observable difference was found between the lower- and higher-order formulations for the current application of hypersonic boundary layer simulations, possibly due to small time step constraints. Ultimately, the second-order accurate approach was preferred due to its simplicity.

3. Transitional boundary layer simulations

In this section, a computational domain with an azimuthal extent of $8^\circ$ is used to simulate transitional hypersonic boundary layers over both sharp- and blunt-nosed cones, with impermeable and permeable walls. Such a width was initially chosen to capture at least two streamwise streaks (see figures 15 and 16), the largest transitional structures present in the current set-up. An auxiliary simulation with a 18-degree-sized domain in the azimuthal direction was performed for a sharp cone at the free-stream conditions corresponding to $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$ and the resulting transitional structures qualitatively matched with the more restrictive domain (Camillo et al. 2023).

3.1. Three-dimensional forcing

In this subsection, the three-dimensional forcing strategy used in the transitional simulations conducted in the current manuscript is discussed. One of the primary objectives of this research is to emulate the experiments conducted in the HEG wind tunnel (Wagner et al. 2013) to gain insights into the dynamics of hypersonic boundary layer flow. To accomplish this, we make use of measurements of the disturbance signal within the HEG wind tunnel as documented by Wagner et al. (2018b) (see figure 4). A simple curve fit yields

(3.1)\begin{equation} A(f) = 0.00448\,{\rm e}^{{-}0.02009f} + 0.00131, \end{equation}

where $f$ here is in kHz and $A(f)={ASD}(p'_w/\bar {p}_w)$ is the amplitude spectral density of the instantaneous wall pressure normalized by its average. It is assumed that such a relation is approximately Reynolds-number independent, as implied by the experimental data.

Figure 4. Amplitude spectral density (ASD) of wall-pressure signal measured by Wagner et al. (2018b) at the surface of a wedge probe normalized by its mean value for different operating flow conditions shown in solid lines (——–). The fit given by (3.1) is shown in dashed lines (– – –). Flow conditions are similar to the ones reported in table 1, used in the current work. Different colours indicate different Reynolds numbers as per the legend.

Figure 4 reveals that the pressure disturbance levels in the HEG tunnel decay monotonically with frequency. In the simulations we impose a disturbance field composed of 29 discrete modes in the 50 kHz–890 kHz range, with a spacing of 30 kHz, that is

(3.2)\begin{equation} f_m = 50 + 30m [\mathrm{kHz}] \quad \textrm{with} \ m \in [0,1,\ldots,28]. \end{equation}

This decision was taken with the intention of incorporating the high-energy/low-frequency portion of noise spectrum observed in the wind tunnel data, while also prudently constraining the lowest frequency to limit the period of the lowest harmonic, which sets the lower bound for the time window required for collecting meaningful statistics. The adopted averaging time window is $\bar {T} = 140\unicode{x2013}180\ \mathrm {\mu }$s, equivalent to 7–9 times the longest time scale present in the disturbance field. The sensitivity of the statistically averaged surface heat flux to $\bar {T}$ is shown in figure 5. A short averaging window fails to remove the signature of the second-mode waves; as the averaging time window increases, low wavenumber oscillations associated with the presence of turbulent spots in the simulations can be observed. Finally, for $\bar {T} > 140\ \mathrm {\mu }$s statistical convergence is achieved.

Figure 5. Statistical convergence analysis showing spanwise and time-averaged wall-heat flux profiles obtained for progressively increasing time averaging windows ($\bar {T}$) shown with different colours (see legend). The longest averaging period shown in each panel is highlighted with a thicker line. Impermeable-wall data (solid line) are compared against porous-wall cases (dashed lines). Data from sharp cone simulations at $Re_m = 2.43 \cdot 10^6\ \text {m}^{-1}$ (a,c) and $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$ (b,d) are compared against blunt-nose cases at $Re_m = 6.40 \cdot 10^6\ \text {m}^{-1}$ (e,g) and $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$ (f,h).

The spatial distribution of the disturbance field, on the other hand, is obtained via spatial Fourier filtering of a series of one-time generated pseudo-random three-dimensional scalar fields. This results in an analytically defined three-dimensional spatial profile with a filtered wavenumber content such that it is resolvable in every direction by the coarsest grid considered in this study, as explained below. This ensures that the imposed pseudo-random disturbance field remains independent of the grid resolution, hence not vitiating the grid-convergence analysis. Furthermore, the intent of removing high wavenumbers in the forcing region is not related to an attempt at ensuring that such content will not be presented within the simulation. High wavenumber content will naturally be generated by nonlinear interactions as the disturbance field advects downstream and the flow transitions from laminar to turbulent. When relatively high wavenumbers (near the grid cutoff) start to occur in the computational domain, however, the dynamic QSV-LES model detects them and increases the magnitude of the added eddy viscosity to ensure an accurate representation of the energy cascade to small scales within the grid resolution in question.

In the current study, the smallest spatial length scale of the disturbance field was set to $\delta _x = \delta _y = 0.5$ mm in the wall-normal and streamwise directions and to $\delta _z = 0.5$ degrees in the azimuthal direction. The precise control of the wavenumber spectrum of the disturbance field is achieved via sharp-spectral filtering of a broadband signal with the previously mentioned minimum allowed length scale in each direction. These scales remain resolvable even at the coarsest grid level considered in the current study, which is $2048\times 96\times 48$.

More details are provided in the following. An initial field is generated by seeding a designated volume with a uniformly distributed pseudo-random scalar field valued between $-1$ and $1$. In the current set of simulations, the disturbance volume has streamwise, wall-normal and azimuthal extents of $\ell _x = 0.075 L_x$, $\ell _y = 0.002$ m and $\ell _{\theta } = 8.0$ degrees, respectively, where $L_x$ is the overall computational domain along the cone surface (see table 2). The disturbance volume is bounded by the inflow plane upstream, by the cone's surface on the bottom wall and by the computational-domain periodic boundaries in the azimuthal direction, as shown in figure 6.

Figure 6. Visualization of forcing volume with dimensions, $\ell _x = 0.075 L_x$, $\ell _y = 0.002$ m and $\ell _{\theta } = 8.0$ degrees, with $L_x$ being the computational-domain extent in the streamwise direction (see table 2). Hyperbolic tangent functions are used to smoothly transition to the unperturbed regions over a quarter of the forcing volume's extent in both the streamwise and wall-normal directions.

A fast Fourier transform is then applied in three dimensions to the previously generated pseudo-random scalar field. Then, the Fourier coefficients with higher wavenumbers than the defined thresholds of $\kappa _x = 1/\delta _x$, $\kappa _y = 1/\delta _y$ and $\kappa _z = 1/\delta _z$, corresponding to spatial modes not resolvable on the coarsest grid considered, are set to zero. The resulting sharp-spectral filtered three-dimensional field of complex Fourier coefficients, comprising a finite number of spatial modes in each direction, can be inverse transformed and evaluated on any grid. As such, it constitutes an analytically defined field denoted as $\varPhi ({\boldsymbol x})$. This process is repeated for each of the temporal modes in (3.2), yielding $m$ statistically independent pseudo-random fields, $\varPhi _m({\boldsymbol x})$.

Finally, the constructed disturbance field ($\phi '$) can be defined as the superposition of $m$ individual scalar fields ($\phi '_m$)

(3.3)\begin{equation} \phi' = \sum_{m = 0}^{28} \phi'_m({\boldsymbol x}, t), \end{equation}

where

(3.4)\begin{equation} \phi'_m({\boldsymbol x}, t) = A(f_m) \Re ( \varPhi_m({\boldsymbol x}) \,{\rm e}^{{-}j2{\rm \pi} f_m t}). \end{equation}

An unsteady three-dimensional pressure disturbance field is then generated via

(3.5)\begin{equation} p'({\boldsymbol x},t) = p_{rms}\frac{\phi'}{\sigma(\phi')}, \end{equation}

where $\sigma (\phi ')$ is the standard deviation of $\phi '$ obtained by averaging over $100\ \mathrm {\mu }$s, and $p_{rms}$ is a user-defined value of the root mean squared pressure amplitude. A corresponding density fluctuation field is obtained via the isentropic relation

(3.6)\begin{equation} \rho'({\boldsymbol x},t) = p'/a_0^2, \end{equation}

where $a_0$ is the local speed of sound, which depends on the base temperature. This achieves a coupled pressure and density field with acoustic-like scaling.

In the forcing region defined in this manner, the total pressure and density are hardly imposed via

(3.7)\begin{gather} \rho = \rho'({\boldsymbol x},t)+\rho_0({\boldsymbol x}), \end{gather}

(3.8)\begin{gather}p = p'({\boldsymbol x},t)+p_0({\boldsymbol x}), \end{gather}

where $p_0({\boldsymbol x})$ and $\rho _0({\boldsymbol x})$ are the background steady laminar pressure and velocity fields deriving from companion precursor or analytical (Blasius and Taylor Maccoll) calculations, sustained by the inlet sponge. In the forcing region, the total instantaneous pressure and density are treated as known quantities in the momentum (2.5) and energy (2.6) equations, which are solved to obtain the instantaneous momentum and energy field. A blending coefficient ($\beta$) is used to smoothly shift from a forcing region into the computational domain downstream, where the added perturbations freely evolve according the governing equations obtained via combination of two tangent hyperbolic functions in the streamwise, $x$, and wall-normal, $y$, directions as

(3.9)\begin{equation} \beta(x,y) =\beta_x \beta_y, \quad \mathrm{where} \ \beta_{\gamma} = \left[ 0.5 \left( 1 - \tanh\left(\alpha \frac{\gamma - r_{\gamma}}{\ell_{\gamma}} \right) \right) \right], \end{equation}

where $\gamma$ is a placeholder variable that should be substituted by either the streamwise direction, $x$ or the wall-normal direction $y$. In the current study, the centre of the transition zones was set to $r_x=0.5 \ell _x + x_i$ and $r_y=0.5 \ell _y$, $x_i = 0.045$ m being the inlet plane location. The parameter $\alpha$ is chosen so that $\beta _\gamma = 0.001$ when $\gamma = \ell _{\gamma }$, i.e. $\alpha = 2 \tanh ^{-1}(1-0.002)$. By choosing this value we guarantee a smooth transition between 1 and 0 inside the forcing volume. This smooth transition allows the thus-created perturbation field to convect out of the forcing region seamlessly and it is then left to evolve naturally along the cone's surface.

The parameter $p_{rms}$ is a scalar parameter that governs the overall magnitude of the imposed disturbance and it was chosen by conducting multiple transitional simulations over a sharp cone with impermeable walls at $Re_m = 4.06 \cdot 10^6\ {\rm m}^{-1}$ with different $p_{rms}$ amplitudes. The transition location, defined as the streamwise distance where the mean heat flux departs from the laminar correlation values, in each of these simulations was compared against experimental observations. The choice of $p_{rms}/p_\infty = 3\,\%$ led to a predicted transition location of $X \approx 0.29$ m, close to the experimentally observed value of $X \approx 0.31$ m at the same flow conditions and nose tip geometry. The same dimensionless ratio of $p_{rms}/p_\infty$ was applied to all other geometries and Reynolds numbers simulated.

The adopted forcing technique is not reflective of the exact physical mechanisms that are present in the HEG wind tunnel experiments since it ignores the interaction between the tunnel noise and the shock, as well as the boundary layer receptivity process. However, it does initiate the transition process without any spatial or temporal bias, mimicking natural transition in the best way possible for boundary-layer-resolved simulations carried out under the shock. Despite these important caveats, there is very good agreement between the simulated three-dimensional transitional flow dynamics (capturing the boundary layer and hence, second-mode frequency modulation, due to the transitional streaks) and the surface pressure spectrum measurements (see figures 13 and 14). A more thorough study of the receptivity process, also accounting for the interaction of the free-stream noise with the attached shock wave, is left for future work.

3.2. Comparison against measured wall-heat flux

To identify the transition from laminar to turbulent conditions, we plot the streamwise profile of the Stanton number ($St$)

(3.10)\begin{equation} St = \frac{q_w}{\rho_e c_p U_e (T_{{ad}} - T_w)}, \end{equation}

which is directly proportional to the heat flux at the wall ($q_w$), with fully laminar and turbulent values predictable via semi-empirical correlations. The subscript $e$ indicates values taken at the boundary layer’s edge, the subscript $w$ refers to values at the wall and, lastly, ${ad}$ indicates adiabatic wall conditions. The adiabatic wall temperature is estimated by the use of the recovery factor $r = Pr^{\alpha }$, where $\alpha = 1/2$ if the flow is laminar, and $\alpha = 1/3$ at turbulent conditions (White & Corfield 2006)

(3.11)\begin{equation} T_{{ad}} = \left(1+r\frac{\gamma-1}{2}M_e^2\right)T_e. \end{equation}

The laminar correlations for the friction coefficient ($C_f$) and for the Stanton number ($St$) are obtained from the solution of the compressible Blasius similarity equations over a sharp cone (Lees 1956)

(3.12)\begin{gather} C_f = \frac{\mu_w \left.\dfrac{\partial u}{\partial y}\right|_w}{\dfrac{1}{2}\rho_e U_e^2} = \frac{\mu_w f''(0) \left.\dfrac{\partial \eta}{\partial y}\right|_w }{\dfrac{1}{2}\rho_e U_e} =\sqrt{6} \mathcal{C}_w Re_x^{{-}1/2} f''(0), \end{gather}

(3.13)\begin{gather}St = \frac{q_w}{\rho_e U_e C_p (T_{{ad}} - T_w)} = \sqrt{\frac{3}{2}} \frac{1}{Re_x^{1/2}Pr} \frac{T_e}{T_{{ad}} - T_w} \mathcal{C}_w g'(0), \end{gather}

where $\mathcal {C}$ is the Chapman–Rubesin parameter defined as $\mathcal {C} = \rho \mu / \rho _e \mu _e$, $\eta$ is the wall-normal similarity variable for boundary layers over axisymmetric bodies and

(3.14a,b)\begin{equation} f' = U/U_e \quad \mathrm{and} \quad g = T/T_e = \rho_e/\rho, \end{equation}

are the similarity functions that describe the profiles of the non-dimensional velocity and temperature, respectively. Since both of these quantities scale with $Re_x^{-1/2}$, one can compute their ratio

(3.15)\begin{equation} \frac{St}{C_f}\frac{T_{{ad}} - T_w}{T_e} = \frac{g'(0)}{2 f''(0) Pr}. \end{equation}

This operation is known as the Reynolds analogy, which is a good approximation for the ratio between heat flux and skin friction for both laminar and turbulent cases. Using this notion, the turbulent $St$ number profile is estimated via the van Driest correlation for the skin friction coefficient of a turbulent flat plate adapted to hypersonic flows (Franko & Lele 2013)

(3.16a–c)\begin{gather} a=\sqrt{r\frac{\gamma-1}{2}M_e^2\frac{T_e}{T_w}}; \quad b = \left(\frac{T_{{ad}}}{T_w} - 1 \right); \quad A = \frac{2a^2 - b}{\sqrt{b^2 + 4a^2}}; \end{gather}

(3.17a,b)\begin{gather}B = \frac{b}{\sqrt{b^2 + 4a^2}}; \quad \frac{\sin^{{-}1}A+\sin^{{-}1}B} {\sqrt{\left(\dfrac{T_{{ad}}}{T_e} - 1 \right) C_f}} = 4.15\log_{10}\left(\frac{Re_x C_f \mu_e}{\mu_w}\right) + 1.7, \end{gather}

which one can solve iteratively. Additionally, White & Corfield (2006) report that the turbulent skin friction of a cone can be related to that of a flat plate through

(3.18)\begin{equation} \frac{C_{f,{cone}}}{C_{f,{plate}}} = (2+d)^{({d}/({1+d}))} , \end{equation}

with values of $d$ between $1/8$ and $1/4$. This leads to a 10 %–15 % increase in turbulent heat transfer from a laminar flat plate to a sharp cone. In this work, the value of $m =1/8$ is used.

Analysing the laminar correlation for the Stanton number, (3.13), one can conclude that, for a given free-stream Mach number ($M_e$) and wall temperature ($T_w$), the Stanton number only depends on $Re_x$ and the similarity variable $g$. This means that the Stanton number can be collapsed if plotted vs $Re_x$, which carries all spatial dependencies. Similar conclusions can be reached for the turbulent correlation and the Reynolds analogy.

Although $Re_x$ is a common way of displaying the spatial variation of the heat flux over the surface of the cone, when used in transitional simulations it can falsely convey the idea that the transition onset is anticipated for a higher free-stream Reynolds number ($Re_m$). For this reason heat-flux results in figure 7 (and onwards) are displayed vs the absolute axial cone coordinate $X$ and wall-heat-flux values are normalized by the corresponding laminar value at $X = 0.275$ m. The use of such a coordinate system allows us to display the anticipation of the turbulence onset that occurs for both sharp- and blunt-nose geometries as $Re_m$ increases, and the transition delay effect due to the small (according to the classification by Softley 1969) nose tip bluntness of 2.5 mm. However, this results in a dependency of the fully turbulent heat-flux correlation value with $Re_m$.

Figure 7. Average heat flux from experimental measurements (symbols) and simulation results (lines) plotted vs the axial streamwise distance from the cone tips. Results from coarser to finer grid refinement levels ($2048\times 96\times 48$, $2560\times 120\times 64$ and $3072\times 144\times 80$) are shown with progressively darker shades of the colour associated with each flow condition (see table 1). For each cone geometry the heat flux is plotted normalized by the laminar value corresponding to the first measurement location at $X = 0.275$ m. Laminar and turbulent correlations (dotted lines) are normalized accordingly. Impermeable walls: solid lines and filled symbols; porous walls: dashed lines and open symbols. Note that experimental measurements over porous walls are not available for a sharp cone.

The experimental data shown in figure 7 and used in the rest of manuscript study can be found in the technical report by Wagner et al. (2019), where we note the following editorial correction: the experimental heat-flux curves measured for the sharp-tipped model are shown in figures 2–7 (section 2's seventh figure) of the report, where subplot a) was mislabelled as $R_{{tip}} = 0.1$ mm. The measurements obtained over a 2.5 mm nosed model are shown in subplot 2-7 b) and figures 5–11 (section 5's twelfth figure).

Various grid refinement levels are reported in figure 7, where it can be observed that the results for the normalized averaged heat-flux distribution along the streamwise direction exhibits the expected grid-dependent behaviour for a dynamic LES. The first step up in grid resolution from a coarse ($2048\times 96\times 48$) to an intermediate ($2560\times 120\times 64$) level yields a larger variation in the predicted transition location and turbulent heat-flux level than the following step from an intermediate ($2560\times 120\times 64$) to a fine ($3072\times 144\times 80$) grid. That is observed for the sharp-nosed cone simulations at $Re_m = [2.43 , 4.06] \cdot 10^6\ \text {m}^{-1}$, as well as the blunt-nosed cone runs at $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$, with either an impermeable or porous wall. The blunt-nosed simulation with impermeable walls at $Re_m = 6.40 \cdot 10^6 \text {m}^{-1}$ exhibits the largest degree of grid sensitivity. However, upon comparing the results from simulations utilizing the finest grid resolution with the experimentally measured values, it becomes evident that comparable heat-flux ramp-up rates and transitional region spatial extents are achieved even at the highest $Re_m$ case. The sharp cone case at $Re_m = 2.43 \cdot 10^6\ \text {m}^{-1}$ transitions from laminar to turbulent over a region of approximately 25 cm in length for both simulated and experimentally measured data; the same can be observed for the $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$ case, transitioning over 12 cm for sharp and 20 cm for blunt cone geometries. For the highest $Re_m$ case of $6.40 \cdot 10^6\ \text {m}^{-1}$ for the blunt cone only, the extent of the transitional region is 15 cm in both simulations and experiments, with predicted turbulent heat-flux values consistently below the experimental values.

The predicted transition onset location for the blunt impermeable runs is located upstream of the experimentally observed value, while the extent of the transition delay due to the porous walls is predicted within the same order of magnitude. Ultimately, the structure of the transitional region is directly affected by the amplitude and spectrum of the fluctuations introduced in the near-tip region through the three-dimensional forcing procedure. Better matching could be achieved if several more (rather costly) simulations were carried out to calibrate the required $p_{rms}$ level more precisely. Various numerical trials (not shown) have revealed that the flow conditions over a blunt-nosed geometry are more sensitive to the $p_{rms}$ level than the sharp cone runs.

It is important to note that the current findings are based on dynamic LESs, which are not intended to encompass all scales present in the turbulent portion of the flow field. Instead, their primary goal is to retain high fidelity of the transitional flow region (as demonstrated by Sousa et al. 2023) and of the larger turbulent length scales, when turbulence occurs. Figure 7 shows how the sensitivity to the grid resolution is more pronounced in the region of turbulent breakdown, especially at the highest Reynolds number considered, $Re_m = 6.40 \cdot 10^6\ \text {m}^{-1}$.

At the intermediate Reynolds number of $4.06 \cdot 10^6\ \text {m}^{-1}$, the simulated results align well with experimental measurements, closely matching the values predicted by the turbulent correlation (3.16a–c) in the case of a sharp cone and only slightly deviating below these values for a blunt cone. For the lowest Reynolds number, $Re_m=2.43 \cdot 10^6\ \text {m}^{-1}$, the simulated values coincide with those predicted by the turbulent heat-flux correlation but fall below the experimentally measured values, which also inherently exhibit some uncertainty.

Two reasons are behind this result: first, the transition over a blunt geometry happens further downstream in the computational domain than in the sharp cone simulations. Due to the present grid arrangement (see figure 1), both wall-normal and azimuthal resolutions decrease in the streamwise direction; second, the higher $Re_m$ increases the grid resolution requirements, shifting the LES results further away from a grid converged value of the fully turbulent heat flux. This is confirmed by the grid sensitivity study in figure 8, where the minimum wall-normal grid spacing in wall units

(3.19)\begin{equation} y^+_{{min}} = \frac{y_{{min}} \rho_w u_{\tau}}{\mu_w} \quad \mathrm{where}\ u_{\tau} = \sqrt{\tau_w/\rho_w}, \end{equation}

is plotted vs the axial distance $X$. This analysis reveals that the computational grid leads to an increase in the minimum wall-normal spacing in plus units in the streamwise direction. The increase is especially pronounced for the highest $Re_m=6.40 \cdot 10^6\ \text {m}^{-1}$. This is a result of the expansion of the computational domain in the streamwise direction approximately following the boundary layer growth, as discussed in Sousa et al. (2023). Figure 8 also shows that the $y^+_{{min}}<1$ condition is met in the laminar region, and only $y^+_{{min}}<2$ in the fully turbulent region for all cases except the higher-Reynolds-number blunt case. This fact is the reason behind the results in figure 7 achieving grid convergence of the predicted the location of transition onset, i.e. where the mean heat flux departs from the laminar correlation values, even for the highest $Re_m$ considered, but not achieving grid convergence of the predicted fully turbulent heat flux.

Figure 8. Minimum local wall-normal coordinate as a function of axial distance for all flow conditions, cone geometries and grid refinement levels considered. Results from coarser to finer grid refinement levels ($2048\times 96\times 48$, $2560\times 120\times 64$ and $3072\times 144\times 80$) are shown with progressively darker shades of the colour associated with each flow condition (see table 1).

The results presented in figures 7 and 8 indicate that, even though not all scales governing the dynamics of the hypersonic boundary layer in question are fully represented by the finest grid considered, the transitional process is accurately simulated with respect to the wall-normal grid resolution. Furthermore, the following discussion on turbulent statistics (§ 3.6) is primarily centred on the sharp conical model, where a larger portion of the turbulent spectrum is captured by the employed LES grid.

Figure 7 shows that the TDIBC correctly models the effect of the porous C/C surface in the flow field. The resulting acoustic energy absorption occurring at the boundary decreases the growth rate of the unstable modes in the hypersonic boundary layer, ultimately inducing a transition delay, here defined as the streamwise shift of the inflection points in the normalized heat-flux curves between impermeable- and permeable-wall cases. The simulations overpredict the extent of the transition delay due to porous walls when compared against the experimental data (only available for blunt geometries): the blunt tip simulations predict a delay of approximately 14 cm at $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$, and 18 cm at $Re_m = 6.40 \cdot 10^6\ \text {m}^{-1}$ against the observed 11 and 16 cm, respectively. Possible sources of uncertainty that can affect the quality of the matching may lie in the application of the ultrasonic bench test impedance measurements (obtained in quiescent air) to a flow with a strong background shear, or to the choice of bypassing of the receptivity process by applying perturbations after the shock.

An increase in $Re_m$ leads to a greater transition delay potential for both simulated and experimentally obtained results. That is related to the increased surface pressure at the associated flow conditions and, therefore, to a higher acoustic energy absorption coefficient ($\beta$) of the C/C surface (Sousa et al. 2019). The results over the sharp cone are less clear because the transition occurs early at $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$ and the space where the porous surface can influence the flow, which starts at $X = 0.182$ m, is reduced at those flow conditions. This fact counterbalances the increased absorption coefficient resulting from a higher surface pressure at the flow conditions related to $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$ when compared with ones related to $Re_m = 2.43 \cdot 10^6\ \text {m}^{-1}$.

Since the results on the dynamics of a transitional hypersonic boundary layer presented and analysed in the current manuscript were obtained via the dynamic QSV-LES procedure, it is necessary to show the distribution of the SFS dissipative terms on the computational domain. Figures 9 and 10 depict the SFS heat flux and the SFS shear stress normalized by the respective values associated with a laminar boundary layer at the same flow conditions at $x = 0.2$ m for both sharp and blunt cone geometries at $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$. Furthermore, figures 9 and 10 use a logarithmic scale spanning 3 decades to highlight all the areas where numerical dissipation is relevant in the simulations.

Figure 9. Distribution of the SFS heat flux (a) and shear stress (b) on the surface of the conical model as well as on a wall-normal slice at the centre of the computational domain. The SFS components were normalized by the laminar values for heat flux and shear stress at $x = 0.2$ m, respectively. The displayed results are related to a sharp-nosed geometry at $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$.

Figure 10. Distribution of the SFS heat flux (a) and shear stress (b) on the surface of the conical model as well as on a wall-normal slice at the centre of the computational domain. The SFS components were normalized by the laminar values for heat flux and shear stress at $x = 0.2$ m, respectively. The displayed results are related to a blunt-nosed geometry at $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$.

For both cases shown, and more evidently in the blunt case, the computational domain can be divided into three distinct regions: one to the left, where the numerical dissipation is inactive; one to the right, where the SFS components are fully active and follow the pattern of turbulent structures; and one in between, where the activation of the SFS dissipation is intermittent. If one compares the location of these regions with the laminar, transitional and fully turbulent areas shown in figure 7, a clear correspondence can be traced. Additionally, in § 3.3, the presence of turbulent spots in the current simulation is discussed. The area where the activation of the sub-filter components spikes in the intermittent region, between the fully laminar and fully turbulent locations, is safely assumed to be connected with the existence of a turbulent spot.

Ultimately, it can be concluded that the numerical dissipation is correctly silenced in regions where the flow is laminar and only large-scale fluctuations are present, allowing for the unhindered growth of the second-mode instability waves. When the flow starts to breakdown and small scales start emerging, the magnitude of the added dissipation ramps up until the flow becomes fully turbulent and a balance is reached. This, in turn, dissipates the turbulent scales that are smaller than the grid resolution, making it possible to conduct these high-order simulations stably and accurately at a reasonable computational cost. Finally, the results in figures 9 and 10 serve as further evidence of the dynamic nature of the QSV-LES technique in transitional hypersonic boundary layers, which was already demonstrated previously by Sousa et al. (2023) in axisymmetric simulations.

3.3. Unsteadiness of transition process

Although the mean heat-flux spatial profiles show a smooth transition from laminar to turbulent flow, the transition process itself is highly unsteady. Figure 11 shows that the standard deviation of the heat flux, $\sigma ({q_w})$, in the transitional region is higher than in the laminar portion, and comparable to the turbulent region. The sharp-nosed cases experience an initial decrease in $\sigma ({q_w})$ in the near-tip region ($X \le 0.25$ m) because of the relaxation of the numerical solution past the forcing region paving the way to a clean modal growth. For that reason, the physically relevant laminar $\sigma ({q_w})$ value for these cases is taken between $X = 0.25$ m and transition. The blunt-nosed simulations do not suffer from the same phenomenon because the entropy layer caused by the curved shock structure makes the perturbations stable in the near-tip region and, ultimately, leads to a small $\sigma ({q_w})$. Furthermore, the curves displayed in figure 11 report similar standard deviation levels in the transitional and turbulent regions across all cases considered but, at the same time, the different stability dynamics present on the laminar flow portion leads to higher $\sigma ({q_w})$ at the near-tip region over sharp cones when compared with their blunt-nosed counterpart. In conclusion, the growth in $\sigma ({q_w})$ is most evident in the flow over blunt-nosed cones, where an approximately 5-fold increase is observed. On the other hand, when a sharp-nosed cone is considered, the standard deviation variation is limited to a twofold increase.

Figure 11. Streamwise profiles of wall-heat flux. Time-averaged results are shown in thicker lines with a shaded region representing one standard deviation. Instantaneous spanwise-averaged values (thin lines) are filtered in the streamwise direction to remove the high wavenumber oscillations. Results are shown over impermeable (——–) and porous (– – –) walls. Laminar and turbulent correlations are shown in dotted lines. See table 1 for the complete colour scheme.

Figure 11 also shows the filtered streamwise distribution of the instantaneous spanwise-averaged heat flux to showcase the possibility that, in some time instants, the local heat flux can be much higher than the overall time-averaged value. This is the result of the presence of fully turbulent spots in the transitional region for both impermeable and porous walls separated by weakly turbulent regions, seen in figure 12 for a 2.5 mm blunt-nosed cone at $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$. The existence of turbulent spots in the transitional region is one of the reasons for an increased standard deviation and, ultimately, this could dictate more stringent requirements for TPS than what would be inferable from the average heat-flux data alone, due to the periodic overshoot to the fully turbulent heat-flux values.

Figure 12. Instantaneous Q-criterion isosurfaces coloured by temperature as well as a numerical Schlieren plane for times: (a) $t_0$; (b) $t_0 + 39.2\ \mathrm {\mu }$s; and (c) $t_0 + 78.4\ \mathrm {\mu }$s for a 2.5 mm blunt-nosed cone at $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$ over impermeable and porous surfaces.

3.4. Frequency content analysis

Figures 13 and 14 show a comparison between the peak wall-pressure disturbance frequencies measured experimentally in the HEG tunnel and the frequency of the maximum pressure oscillation amplitude of the equivalent signal gathered through numerical simulations. Furthermore, the temporal spectrum of the wall-pressure fluctuation is extracted at each grid location on the surface of the cone, and the frequency of maximum spectral energy content, $f_{max}$, is identified. The spatial distribution of $f_{max}$ along the surface of the conical model reveals the presence of streaks where second-mode waves experience a downshift in frequency.

Figure 13. At the centre (b), the frequencies connected to the maximum amplitude of the spanwise-averaged power spectrum density of the wall-pressure oscillation signal over a sharp-nosed cone are compared with the second-mode peak frequencies reported in Wagner et al. (2019). Simulated results are shown in solid lines and a shaded region equivalent to one standard deviation; experimental results are shown in symbols with error bars. The linearly unstable region for each free-stream condition is shown around the location of maximum growth, indicated by dash-dotted lines (– - – - –). At the top (a) and bottom (c), the spatial distribution of the frequencies associated with the highest perturbation magnitude is shown. The red arrow vertical arrow on the left indicates the diameter of the sensitive surface of the printed circuit board (PCB) pressure transducers in comparison with the streaks present in the flow. See table 1 for the complete set of colours and symbols used.

Figure 14. At the centre (b), the frequencies connected to the maximum amplitude of the spanwise-averaged power spectrum density of the wall-pressure oscillation signal over a 2.5 mm blunt-nosed cone are compared with the measured second-mode peak frequencies reported in Wagner et al. (2019) in experiments at the HEG wind tunnel. Simulated results are shown in solid lines and a shaded region equivalent to one standard deviation; experimental results are shown in symbols with their respective error bars. Additionally, the linear unstable region for each free-stream condition is shown around the location of maximum growth, indicated by dash-dotted lines (– - – - –). The masked zone represents the approximate region where the entropy layer stabilizes the second mode. At the top (a) and bottom (c), the spatial distribution of the frequencies associated with the highest perturbation magnitude is shown. A red arrow is used to indicate the relative size of the sensitive surface diameter of the PCB pressure transducers in comparison with the streaks present in the flow. See table 1 for colour scheme.

Figures 13 and 14 also display the frequencies with the maximum growth rate predicted by LST analysis. These are shown with a dot-dashed line (– - – - –) with a shaded region showing the frequency band of positive growth rate between the higher and lower neutral frequencies. The spectrum of a small-amplitude linear second-mode wave packet in a laminar boundary layer is expected to lie between the maximum growth rate line and the upper neutral frequency. In this region, the initially low-amplitude perturbations at a certain frequency would have experienced most of the expected integrated growth and should have the highest magnitude. The LST results are gathered assuming a fully laminar boundary layer over the surface of the conical geometries for each flow condition and serve as an important reference: the deviation from the LST results, in fact, yields valuable insights into the state of the boundary layer flow at various streamwise positions. The LST solver used in the current work is based on a novel numerical approach relying on Laguerre polynomials and its description can be found in Sousa et al. (2023).

The experimental measurements for the lowest $Re_m$ cases, for both sharp (figure 13) and blunt tip cones (figure 14), lie within the expected LST-predicted range which serves as evidence that the boundary layer, at these conditions, remains laminar over the entire surface of the conical model. On the other hand, the experimental data gathered at higher $Re_m$ consistently report lower frequencies of maximum amplitude for a certain streamwise distance from the tip. Most of the measurements for the higher $Re_m$ are made in regions where disturbance amplitudes are large enough to induce nonlinear interactions and mean flow distortions, causing measurable departure from the laminar base flow dynamics predicted by linear theory. In fact, results from the present three-dimensional transitional simulations confirm the frequency shift to lower values observed in the experimentally measured values.

Numerical results reported in figures 13 and 14 show that the most amplified frequencies in the region close to the tip of the simulated conical model are consistent with LST predictions, while progressively shifting to lower values as the flow approaches the transition location, where the mean heat flux departs from the laminar correlation values. In the turbulent region, higher frequencies and standard deviations are observed in comparison with the laminar portion of the flow due to the broadening of the spectra in that region. This will be discussed in § 3.6 but the reader is directed to figure 26 for a visual representation of the discussed phenomenon. The surface distribution of the frequencies corresponding to the maximum pressure spectrum amplitude reveals patches of relatively lower frequencies surrounded by regions of higher values. This indicates the presence of streak-like distributions of the frequency content of the instantaneous surface pressure field.

Although a high variability of the frequency of the most amplified modes’ values is observed in the surface plots, both experimental and spanwise-averaged numerical results strongly support the presence of a shift to lower frequencies. Two reasons are behind this behaviour: first, the width of the low-frequency streaks is smaller than the sensitive diameter of the PCB sensor used in the HEG experiments, and both experimental and numerical results can be considered spanwise averaged; second, the low-frequency streaks carry the bulk of the perturbation energy. This phenomenon can be observed in figures 15 and 16 for sharp- and blunt-nosed cone geometries, respectively. These figures show the relation between the low heat-flux streaks and the lower-frequency regions of the most amplified modes, which happen in spatially similar locations at the surface of the cone for all cases considered. Furthermore, an instantaneous visualization of the pressure oscillation field at the wall indicates the concentration of the high-amplitude oscillations at such locations. In conclusion, a spanwise average of the pressure fluctuation spectrum will always yield an overall downshift of observed frequency content with respect to linear stability predictions.

Figure 15. Surface distribution of averaged wall heat flux (a,c) as well as instantaneous visualization of normalized wall-pressure oscillations (b,d) are shown for a sharp-nosed cone at different free-stream conditions. A red arrow is used to indicate the diameter of the sensitive surface of the PCB pressure transducers used in Wagner et al. (2019), comparing it against the streaks present in the flow. See table 1 for the complete set of colours and symbols used.

Figure 16. Surface distribution of averaged wall heat flux (a,c) and instantaneous visualization of normalized wall-pressure oscillations (b,d) are shown for a 2.5 mm blunt-nosed cone at different free-stream conditions. A red arrow is used to indicate the diameter of the sensitive surface of the PCB pressure transducers used in Wagner et al. (2019), comparing it against the streaks present in the flow. See table 1 for the complete set of colours and symbols used.

The presence of high and low heat-flux streaks in a transitional hypersonic boundary layer was observed previously by Chynoweth et al. (2014) on a flared cone via a temperature sensitive paint technique, and by Gray & Schneider (2020) on a straight slender cone via infrared imaging in Purdue's Boeing/AFOSR Mach-6 Quiet Tunnel (BAM6QT) at low-noise conditions. Direct numerical simulations of a controlled transition scenario at the same flow conditions were performed for both flared (Hader & Fasel 2019) and straight slender (Hader, Leinemann & Fasel 2020) cones. Results for either geometry show that the fundamental resonance breakdown, characterized by the interaction of an axisymmetric primary wave with an oblique secondary wave of the same frequency, is responsible for the formation of such streaks and is the principal mechanism governing transition at the studied flow conditions. Moreover, pseudo-random forcing simulations performed by Hader & Fasel (2018) over a flared cone yielded a similar pattern of streamwise heat-flux streaks and similar dominant frequencies as was observed in the aforementioned experiments. Similarly, in the current simulations, the imposed pseudo-random forcing triggers a wide range of modes, including the most unstable primary and secondary modes. In turn, these experience a higher streamwise growth and induce the formation of the streak pattern observed in figures 15 and 16. In practice, these heat-flux streaks meander in the azimuthal direction (see figure 17) with time scales of approximately 0.2 ms, comparable to overall simulation time. Ultimately, larger simulations times should lead to less pronounced, i.e. more azimuthally homogenous, mean heat-flux distributions. A more in depth discussion on the actual breakdown mechanisms that occur at the experimental conditions studied is left for future work. The strong cooling effects are speculated to be responsible for the large streamwise extent of the simulated streaks and the low azimuthal wavenumbers.

Figure 17. Contours of instantaneous heat flux at the surface of a blunt cone at $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$ for two different instants in time showcase the meandering of the low heat-flux streaks.

Nonetheless, the cause of the shift to low frequencies is explained by spanwise variation of the boundary layer thickness caused by the presence of streamwise vorticity in the flow field. It brings high-temperature fluid from the bulk of the boundary layer closer to the wall, as well as bringing cold fluid from the near-wall region further into the boundary layer. This is phenomenon is observed in the instantaneous temperature field contours shown in figure 18 and is responsible for the observed streak pattern. The region where low heat flux is predominant is connected to a locally thicker boundary layer and, consequently, a lower-frequency perturbation spectrum.

Figure 18. Instantaneous temperature profiles at different distances from the cone tip are shown for a 2.5 mm blunt-nosed cone at $Re_m = 6.40 \cdot 10^6\ \text {m}^{-1}$ (a) and $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$ (b).

In conclusion, the results gathered in figures 13–16 indicate that, for the free-stream conditions studied, the initial second-mode instability promotes the formation of streamwise streaks which, in turn, lead to a downshift in the expected frequency spectrum.

An interesting result is observable for all cases considered but especially for the sharp-nosed cone at $Re_m = 2.43 \cdot 10^6\ \text {m}^{-1}$. Figure 15 shows that the low heat-flux streaks formed in the transitional region of the flow persist well downstream of $X = 0.7$ m, i.e. the streamwise location where the peak in heat-flux numerical measurements occurs. From that point on the flow regime is considered fully turbulent. The continued presence of the low heat-flux streaks in the fully turbulent region indicate that memory from the earlier stages is retained and affects the turbulent flow. Moreover, the contours of instantaneous pressure oscillation at the surface indicate that there is a correlation between the location of the persistent low heat-flux streaks and the second-mode-like pressure oscillations that occur in the turbulent region. Similar wave trains of wall-pressure oscillations were observed in hypersonic turbulent channel flow simulations by Chen & Scalo (2021) at $M_b = 6.0$ over impermeable and porous walls. In that work, it was also reported that walls with relatively high permeability, much higher than what is measured for the porous C/C material considered in this work, can render such waves linearly unstable, further increasing their importance in the flow field dynamics.

3.5. Comparison against wall-pressure oscillation frequency spectra

In this subsection, the time signal of wall-pressure oscillations measured both experimentally and by means of simulations are compared against each other with respect to their frequency content at different streamwise distances from the cone's nose tip for two flow conditions: $Re_m = 1.46 \cdot 10^6\ \text {m}^{-1}$ and $Re_m = 4.0 \cdot 10^6\ \text {m}^{-1}$.

At the lower $Re_m$, the boundary layer over the blunt cone does not become turbulent, as is evidenced by the measurements reported in figure 14 that show recorded frequencies of maximum oscillation amplitude being within the second-mode unstable region predicted by linear stability. Additional evidence can be found in Wagner et al. (2019), where plot 2-7 gathers the normalized heat flux at the flow condition in question and shows measurements close to the expected value for the heat flux of a laminar boundary layer until the end of the conical model. Measurements of a fully laminar flow over the blunt-nosed cone are a valuable opportunity for verification of the adopted numerical set-up. First, it is possible to simulate the instability mechanisms by using axisymmetric simulations since the unstable mode with the highest growth is longitudinal. Second, the region of interest in the frequency domain is the one around the main instability amplitude peak ( $f=[250,275]$ kHz) since both higher and lower frequencies are stable at each streamwise location and, therefore, will be exponentially decaying. Third, in the region of linear instability growth, central focus is given to the unstable frequency region and to its growth rate. This is because these data values are independent of the actual amplitude of the signal, given that its magnitude is low enough to only induce linear perturbations to the laminar base state. Ultimately, this allows for the scaling of data without any loss of generality, as done in figure 19.

Figure 19. Wall-pressure fluctuations of normalized amplitude spectrum (AS) at two locations on the surface of the cone. The experimental conditions correspond to a blunt cone with $R_{{tip}} = 2.5$ mm at $Re_m = 1.46 \cdot 10^6\ \text {m}^{-1}$. Experimental data are shown with symbols, axisymmetric simulation data are shown with blue lines and LST results from Wartemann et al. (2015) are shown in magenta lines. Impermeable walls: ——– and $\bigstar$; porous surfaces: – – – and .

For the blunt-nosed cone case at $Re_m = 1.46 \cdot 10^6\ \text {m}^{-1}$, PCB sensors at two different locations downstream were able to measure the second-mode pressure signal above the noise floor and its frequency content could be recovered for both impermeable and porous walls. An axisymmetric simulation of the same conditions with $n_x \times n_y = 3072 \times 144$ points was performed with both impermeable and impedance boundary conditions, and perturbed via the pseudo-random forcing outlined in § 3.1 with the same amplitude as the three-dimensional runs previously reported. The normalized amplitude spectra (AS) of the wall-pressure fluctuations, obtained via a discrete Fourier transform, for both experimental and simulated data, are shown in figure 19, and good agreement is presented in both the relative amplitude growth between $X = 650$ and $X = 785$ mm and relative attenuation of the signal caused by the presence of the C/C surface.

The frequency of maximum amplitude observed in the axisymmetric simulations displays a shift to higher frequencies of $\approx$8 % in comparison with the experimentally measured spectra. A frequency shift of such order of magnitude is common in comparisons between experimental results and linear stability analyses. A systematic over-prediction of the order of 10 % in the frequencies with maximum amplitude at each spatial location was also observed by Wartemann et al. (2014) when comparing the results obtained from the DLR stability code, NOLOT (Hein et al. 1994; Wartemann, Lüdeke & Sandham 2012), with the same experimental results. Moreover, in the study performed by Wartemann et al. (2014), an assessment of the influence of inaccuracies in the measurements of the free-stream condition was performed by increasing the free-stream $Re_m$ by a factor of 5 %, of the order of the uncertainties present in the HEG shock tunnel experimental runs. Ultimately, it was shown that a 5 % uncertainty on $Re_m$ can lead to variations of the maximum frequency distribution of the order of 10 %. Wartemann et al. also looked into inaccuracies in the wall temperature determination and concluded that a variation of 10 % in the wall temperature leads to an uncertainty of around 2 % in the maximum frequency spatial profile. A similar result was obtained in the current simulations (not shown), as a small decrease in the predicted frequency was obtained by increasing the wall temperature from 300 to 330 K. A different study performed by Alba et al. (2010) at the BAM6QT comparing experimental results and the stability theory predictions by STABL (Johnson & Candler 2005) also reported inaccuracies that ranged from 2 % to 14 %. In an effort to reduce the frequency shift, Wartemann et al. (2019) performed base flow simulations including the nozzle, test chamber and cone model itself and reported a reduction to 5 % in the error, which indicates another possible cause for inaccuracies in the base flow determination.

Moving to a higher $Re_m$ condition, at $Re_m = 4.0 \cdot 10^6\ \text {m}^{-1}$, the flow over the blunt-nosed cone is not fully laminar anymore and the simplifications previously made for the simulation and analysis of the lower $Re_m$ case are not applicable. First, three-dimensional QSV-LES runs at these conditions were necessary and data collected from them were used for the comparison against experimental measurements. Second, a higher range of frequencies is present in the transitional flow, up to 1 MHz as shown in figure 20. Finally, a simple scaling of the data does not suffice in this case and a comparison against the magnitude of the measured and simulated values obtained is presented in figure 20, showing good agreement.

Figure 20. Amplitude spectra of the wall-pressure fluctuations at two locations on the surface of the cone. The experimental conditions correspond to a blunt cone with $R_{{tip}} = 2.5$ mm at $Re_m = 4.0 \cdot 10^6\ \text {m}^{-1}$. Experimental data are shown with symbols and three-dimensional simulation data are shown with lines. Relative amplitude of the numerically measured ASD curve over porous walls is set to be equal to the experimental values at $X = 650$ mm to showcase amplitude evolution. Impermeable walls: ——– and $\bullet$; porous surfaces: – – – and $\circ$.

For the aforementioned flow conditions, the PCB sensors at $X = 650$ mm and $X = 785$ mm are in the transitional region, being in the early transition stage for the run over porous surfaces and in mid-to-late stage in the run over impermeable walls, as shown in figure 20. Although the flow is transitional, the first sensor is able to measure a distinct frequency band over the noise floor for both surfaces, with impermeable-wall results exhibiting lower frequencies than a porous wall and with both peak locations being below what is expected by a linear stability analysis of the corresponding laminar base flow (see figure 14). The different stages of transition to turbulence have led to distinct levels of mean flow departure from the laminar base flow, i.e. thicker boundary layers, and ultimately, have decreased the frequency band of the local instability further for the impermeable-wall cases, rather than for the porous wall. The results from the sensor at $X = 785$ mm for the impermeable wall shows a broad amplitude spectrum, an indication of turbulence onset. Moreover, the amplitude spectrum over the porous walls at this location shows a further peak shift to low frequencies when compared with the shift at $X = 650$ mm, indicating a more advanced stage of transition.

These experimental results are also compared against numerically obtained AS curves in figure 20 but, since the current simulations exhibit transition to turbulence earlier than the experiments, a turbulent AS curve is observed for both locations over an impermeable wall. At $X = 785$ mm, the numerically obtained turbulence AS spectrum over impermeable walls exhibits similar behaviour to the experimental results, with a relatively flat AS curve. Additionally, a good agreement between the amplitude of the experimentally measured and simulated wall-pressure signals when analysed in the frequency domain is obtained. Furthermore, a comparable decrease in maximum amplitude as the flow advects downstream is observed in the numerical and experimental results when a porous surface is considered.

3.6. Path to turbulence over impermeable and porous walls

The evolution of the averaged flow variables and fluctuation intensity profiles as the boundary layer transitions to a fully turbulent state over a sharp-nosed cone is analysed. Similar results are found for the flow over a blunt cone, and although they are not repeated here for the sake of conciseness, they can be found at Sousa et al. (2022).

Figure 21 shows the normalized streamwise velocity and temperature profiles for both impermeable and porous surfaces as a function of the distance from the tip of the cone at $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$ and $Re_m = 2.43 \cdot 10^6\ \text {m}^{-1}$. Initially, the profiles for both quantities are close to laminar. This indicates that the imposed perturbation field does not have enough energy to distort the mean flow. As the disturbances move downstream their amplitude grows, a gradual departure from the laminar profile is observed. The flow over an impermeable surface shows a faster evolution towards a thicker boundary layer and a higher near-wall gradient than the flow over porous walls, showing an effective transition delay.

Figure 21. Reynolds-averaged streamwise velocity (a,c) and temperature (b,d) profiles for the following Reynolds numbers: (a,b) $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$ (green); and (c,d) $Re_m = 2.43 \cdot 10^6\ \text {m}^{-1}$ (grey) over a sharp cone with impermeable (——–) and porous (– – –) walls.

Further insight into the effects of porous walls on the transition path to a fully turbulent state can be gathered by plotting the mean streamwise velocity normalized in wall units. In the study of compressible turbulence, the transformation developed by van Driest (1951)

(3.20)\begin{equation} U^+_{{VD}} = \int_0^{\bar{U}/u_{\tau}} \left(\frac{\bar{\rho}}{\rho_w}\right)^{1/2} \, {\rm d} (\bar{U}/u_{\tau}), \end{equation}

and plotted against the local wall coordinate

(3.21)\begin{equation} y^+= \frac{y \rho_w u_{\tau}}{\mu_w} \quad \mathrm{where}\ u_{\tau} = \sqrt{\tau_w/\rho_w}, \end{equation}

has proven successful collapse of the data onto the incompressible log law for adiabatic wall flow conditions. Increased heat transfer rates, however, lead to a departure from the incompressible log law

(3.22)\begin{equation} U^+= \frac{1}{\kappa} \ln(y^{+{/}*}) + C \quad \mathrm{where} \ \kappa = 0.41, \end{equation}

exhibiting a larger intercept ($C$), decreasing the slope of the viscous sublayer (Trettel & Larsson 2016; Zhang, Duan & Choudhari 2018). This behaviour is also observed in the current simulations where $T_w \approx 0.1\ T_{{{ad}}}$, as shown in figure 22.

Figure 22. Mean streamwise velocity profiles $m$ for Reynolds numbers: (a,b) $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$ at $x = 0.35, 0.45,\ 0.55$ m (green); and (c,d) $Re_m = 2.43 \cdot 10^6\ \text {m}^{-1}$ at $x = 0.55,\ 0.70,\ 0.85$ m (grey) over a sharp-nosed cone, scaled based on transformations by Trettel & Larsson (2016) (a,c) and by van Driest (b,d). Results are shown over impermeable (——–) and porous (– – –) walls. The reference log laws are displayed (dotted lines) for $\kappa = 0.41$ and $C = 5.2$. In symbols (), DNS data by Zhang et al. (2018) for $M = 7.87$ and $T_w/T_{{ad}} \approx 0.5$ over a flat plate.

To achieve the desired collapse in the viscous sublayer, the mean velocity transformation proposed by Trettel & Larsson (2016)

(3.23)\begin{equation} U^+_{{TL}} = \int_0^{\bar{U}/u_{\tau}} \left(\frac{\bar{\rho}}{\rho_w}\right)^{1/2} \left[ 1 + \frac{1}{2} \frac{1}{\bar{\rho}} \frac{\partial \bar{\rho}}{\partial y}y - \frac{1}{\bar{\mu}} \frac{\partial \bar{\mu}}{\partial y}y\right]{\rm d} (\bar{U}/u_{\tau}), \end{equation}

is used and plotted against the semi-local wall coordinate (Morkovin 1962; Huang, Coleman & Bradshaw 1995) given by

(3.24a,b)\begin{equation} y^* =\frac{\bar{\rho}(y)u_{\tau}^*}{\bar{\mu}(y)}y \quad \mathrm{and}\quad u_{\tau}^* = \sqrt{\overline{\tau_w}/\bar{\rho}(y)}. \end{equation}

Figure 22 shows that, by using this transformation, the correct viscous sublayer slope is recovered once fully turbulent conditions are achieved but, although the correct slope of the log layer seems to be recovered, the Trettel & Larsson (TL) transform is not able to collapse the data onto the incompressible log law. Although the flow over impermeable walls reaches fully turbulent conditions earlier in comparison with the flow over porous walls, no significant difference in the turbulent mean flow profiles is observed.

Zhang et al. (2018) performed a comprehensive analysis of DNS data in the Mach-number range $M = [2.5, 14]$ and $T_w/T_{{{ad}}} = [0.18,1.0]$ and reported that, although the transformation by Trettel & Larsson (2016) is able to collapse the log-law intercept of hypersonic channel flows with incompressible data, it leads to an upwards shift when applied to high-speed flat plate boundary layers. An example of such behaviour, also observed in the current dataset, is displayed in data from a DNS simulation with $M = 7.87$ and $T_w/T_{{{ad}}} \approx 0.5$, also plotted in figure 22. It is worth noting that small differences between the boundary layer turbulent profiles are expected for flows over a flat plate and conical geometry, similarly to the ones observed in the skin friction coefficient given by (3.18), for example.

Figure 23 shows profiles of velocity fluctuation intensity in the streamwise, wall-normal and azimuthal directions. As the flow develops towards a turbulent state, which happens faster for the flow over impermeable walls when compared with the one over porous walls, the characteristic peak of streamwise velocity fluctuations forms in the near-wall region. Moreover, the velocity fluctuation intensity in the azimuthal directions increase gradually as transition develops.

Figure 23. Velocity fluctuation intensity profiles for (a) $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$ (green) and (b) $Re_m = 2.43 \cdot 10^6\ \text {m}^{-1}$ (grey) over a sharp cone with impermeable (——–) and porous (– – –) walls scaled by the local boundary layer thickness.

A non-local scaling (3.24a,b) was used in figure 23 because better collapse between different wall cooling and compressibility levels is achieved when the variation in density in the wall-normal direction is accounted for (Duan, Beekman & Martin 2010). However, the collapse is not perfect. Higher peaks of normalized turbulence intensity in the streamwise direction are observed for cooled walls in comparison with adiabatic conditions at the same Mach number (Duan et al. 2010) as well as for higher Mach numbers at adiabatic wall conditions (Duan, Beekman & Martin 2011). At the current flow conditions of $M \approx 7.5$ and $T_w \approx 0.1 T_{{ad}}$ the observed peak value is around $4$, which is higher than the reported $\approx 2.9$ and $\approx 3.2$ for a Mach 5 flow over an adiabatic wall surface and over a cooled wall with $T_w \approx 0.2 T_{{ad}}$, respectively (Duan et al. 2010).

Figure 23 also serves the purpose of assessing if the presence of surface porosity affects the velocity fluctuation intensity profiles. Comparing the results obtained at the furthest downstream locations probed for both $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$ and $Re_m = 2.43 \cdot 10^6\ \text {m}^{-1}$ free-stream conditions, one can conclude that there is not a significant effect of porous walls representative of a C/C material on the flow field, apart from delaying the transition to turbulence. It is noted, however, that the current modelling approach neglects distributed surface roughness effects, which may have an important hydrodynamic effect on the flow.

The contours of the Reynolds shear stress and heat flux are shown in figures 24 and 25, respectively. These represent the intensity of the turbulent transport of momentum and thermal energy in the wall-normal direction. The sign of the Reynolds shear stress contours indicate that high-momentum fluid is transported by turbulent fluctuations towards the wall and low-momentum fluid in the opposite direction, across the boundary layer. This mean transport of momentum leads to a higher velocity gradient near the wall and, ultimately, a higher skin friction coefficient. On the other hand, the turbulence-induced transport leads to the mixing of the highly concentrated hot fluid present in the inner portion of a laminar hypersonic boundary layer with colder fluid present both closer to and further away from the wall. That is observed in the turbulent heat-flux contours as a predominant negative region away from the wall with a simultaneous positive region close to the wall. The mean transport of thermal energy towards the wall is responsible for the increase of the mean temperature in that region and, ultimately, for the increased the heat flux due to transition.

Figure 24. Contours of Reynolds shear stress for $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$ (a) and $Re_m = 2.43 \cdot 10^6\ \text {m}^{-1}$ (b) over a sharp cone with impermeable and porous walls.

Figure 25. Contours of turbulent heat flux for $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$ (a) and $Re_m = 2.43 \cdot 10^6\ \text {m}^{-1}$ (b) over a sharp cone with impermeable and porous walls.

Moreover, a peak in energy transport intensity is observed at the end of the transitional region, which is moved downstream when a porous surface is present, in both Reynolds shear stress and heat-flux contours for the two free-stream conditions considered. Differently than what was observed by Franko & Lele (2013), this peak does not translate to a skin friction coefficient or heat-flux overshoot with respect to turbulent correlations; rather, it is concentrated in the transitional region where the biggest change in mean profiles happen. Exploratory transitional simulations with different wall-to-adiabatic temperature ratios (not shown) show that the overshoot starts to occur when the wall temperature approaches the adiabatic wall temperature. At this point it is not clear, however, if the resolution of the LES is insufficient to capture the overshoot in heat flux when cooled walls are considered or if this is a physical phenomenon. Ultimately, further study needs to be performed before a conclusion is made.

Near the cone nose tip for both $Re_m$ values, the pressure oscillation spectra over the two surface types shown in figure 26 display clear harmonics of a fundamental frequency being slightly attenuated over porous walls in comparison with impermeable walls. Moving further downstream, it is possible to observe that the flow over impermeable walls develops more rapidly into a fully turbulent state, characterized by a broader energy distribution, in comparison with the flow over a porous wall. Near the end of the computational domain the recovered spectra over impermeable and porous walls seem to converge, which is further evidence that the C/C's surface porosity delay's transition to turbulence but does not have a significant effect on the turbulence itself.

Figure 26. Energy spectra of wall-pressure fluctuations for $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$ (a) and $Re_m = 2.43 \cdot 10^6\ \text {m}^{-1}$ (b) over a sharp cone with impermeable (——–) and porous (– – –) walls.

Figure 27 shows the distribution of wall-pressure disturbance energy amongst scales in both time and azimuthal direction at different streamwise locations. The wavelengths associated with the wavenumbers in question are $8^{\circ }/k_\theta$. Near the nose tip, only high-frequency components that span wavenumbers in the range $k_\theta \in [0, 5]$ are present for both $Re_m$ values considered. As the flow moves downstream and further into the transitional region, a low-frequency and low wavenumber component is induced, more rapidly over impermeable walls than over porous walls. Near the end of the computational domain, differences in the distribution of energy amongst scales for the two boundary conditions considered are not noticeable. The results shown in figure 27 also serve as evidence that the spanwise resolution is adequate to solve for the flow conditions considered as 3 decades of wall-pressure disturbance energy are captured with approximately half of the azimuthal modes available.

Figure 27. Contour plots of the energy spectra of wall-pressure fluctuations for different frequencies and azimuthal wavenumbers at $Re_m = 4.06 \cdot 10^6\ {\rm m}^{-1}$ (a) and $Re_m = 2.43 \cdot 10^6\ {\rm m}^{-1}$ (b) over a sharp cone with impermeable and porous walls.

4. Conclusion

Numerical simulations of hypersonic transition delay due to the presence of porous walls were carried out. The simulations were designed to reproduce experiments in the HEG wind tunnel on a 7$^\circ$ half-angle cone at $M_\infty =7.4$, which measured transition delay when porous C/C surfaces were present (Wagner et al. 2013; Wagner 2014; Wagner et al. 2019). The acoustic absorption of such a material was experimentally measured and this information was used to reconstruct the real and imaginary parts of the broadband complex impedance in the frequency domain via low-order models. A broadband TDIBC was then used to introduce the porosity effects in the high-fidelity simulations. The use of this technique allows the application of the impedance boundary at all frequencies concomitantly, which expands the modelling capacity used in previous numerical studies (Egorov et al. 2008; Wang & Zhong 2011; Lukashevich et al. 2012), which were restricted to a frequency-by-frequency analysis. The current approach allows us to assess the influence of a porous surface in a natural transition scenario, from the initial second-mode growth to the full breakdown to turbulence.

Three-dimensional simulations were performed and the boundary layer was excited with grid-independent filtered pseudo-random pressure perturbations. These are advected downstream and amplified by the instability mechanisms inside the boundary layer, leading to nonlinear mode interactions and turbulent breakdown. Throughout the transition process, a noticeable downshift of the unstable frequency band was observed, of approximately 100 kHz for the sharp-tipped conical model and 200 kHz for the blunt-nosed model. This observation suggests a potential influence stemming from nonlinear perturbations causing mean flow distortion. Moreover, the presence of the impedance boundary condition, simulating the porous surface, resulted in a transition delay of a similar magnitude, with a slight over-prediction when compared against experimental data obtained on the blunt-nosed model. In the simulations, a transition delay of around 14 cm was predicted at $Re_m = 4.06 \cdot 10^6\ \text {m}^{-1}$ and 18 cm at $Re_m = 6.40 \cdot 10^6\ \text {m}^{-1}$, in contrast to the observed experimental values of 11 and 16 cm, respectively. Although the presence of the porous surface was critical in delaying transition, it did not yield a significant change in the structure of the near-wall turbulence or the wall-pressure oscillation spectra.

As a final remark, the numerical approach adopted here assumes that the characteristic length scales of the boundary layer instability waves and near-wall turbulence are much larger than the surface pore sizes of the C/C material. This justifies the imposition of a spatially averaged, homogeneous impedance. Since the pores in a C/C surface are of the order of $1\unicode{x2013}100\ \mathrm {\mu }$m, this assumption is acceptable in regions of modal growth, where the wavelengths involved are of the order of $1\unicode{x2013}10$ mm; however, some transitional and turbulent structures exhibit smaller scales, of the order of a fraction of a millimetre, which may be hydrodynamically affected by the physical microscale structure of the surface porosity.