2.1. Computational set-up

The numerical study performed in this work is based on experiments conducted by Wagner et al. (Reference Wagner, Kuhn, Schramm and Hannemann2013) in the HEG. A complete description of such facility can be found in Hannemann et al. (Reference Hannemann, Martinez Schramm, Wagner and Ponchio Camillo2018). 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. (Reference Wagner, Kuhn, Schramm and Hannemann2013) 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 (Reference Taylor and Maccoll1933) inviscid solution with a viscous solution for the boundary layer. The latter is derived by applying the Mangler transformation for bodies of rotation (Mangler Reference Mangler1948) to the compressible boundary layer similarity solution for a perfect gas over a flat plate (Illingworth Reference Illingworth1949), as done by Lees (Reference Lees1956) 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 Reference Kroll and Fassbender2006).

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. (Reference Sousa, Wartemann, Wagner and Scalo2023) 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. (Reference Wagner, Kuhn, Schramm and Hannemann2013).

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 Reference Sousa and Scalo2022a), 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. (Reference Sousa, Wartemann, Wagner and Scalo2023). 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 (Reference Jordan1999) and Nagarajan, Lele & Ferziger (Reference Nagarajan, Lele and Ferziger2007), 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 (Reference Vreman, Geurts and Kuerten1995), who showed those are negligible compared with other SFS terms. Following Nagarajan, Lele & Ferziger (Reference Nagarajan, Lele and Ferziger2003), 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 (Reference Sousa and Scalo2022a). 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 Reference Sousa and Scalo2022b) of the original manuscript (Sousa & Scalo Reference Sousa and Scalo2022a) outlining the QSV methodology. The full description of the spectral modulation operation done via $\tilde {G}_{{qsv}}$ can be found in Sousa & Scalo (Reference Sousa and Scalo2022a). 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 (Reference Sousa and Scalo2022a) 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 Reference Lele1992) 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. Reference Nagarajan, Lele and Ferziger2003). 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. (Reference Sousa, Wartemann, Wagner and Scalo2023). 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. (Reference Sousa, Wartemann, Wagner and Scalo2023), 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 Reference Thirani, Gupta and Scalo2020), 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. Reference Kinsler, Frey, Coppens and Sanders1999). 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. Reference Fedorov, Shiplyuk, Maslov, Burov and Malmuth2003; Fedorov & Malmuth Reference Fedorov and Malmuth2008). 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 (Reference Fung and Ju2004) and the numerical integration strategy in a compressible Navier–Stokes solver by Scalo et al. (Reference Scalo, Bodart and Lele2015), and subsequent work by Lin, Scalo & Hesselink (Reference Lin, Scalo and Hesselink2016) and Douasbin et al. (Reference Douasbin, Scalo, Selle and Poinsot2018). 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 Reference Rienstra2006). 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 (Reference Gustavsen and Semlyen1999), 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. (Reference Wagner, Schramm, Dittert, Sousa, Patel and Carlo2018a) and Sousa et al. (Reference Sousa, Wartemann, Wagner and Scalo2023).

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. (Reference Sousa, Wartemann, Wagner and Scalo2023).

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 Reference Fung and Ju2004), 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. (Reference Scalo, Bodart and Lele2015) 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 (Reference Fung and Ju2004), 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 (Reference Chen and Scalo2021) 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.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. Reference Wagner, Kuhn, Schramm and Hannemann2013) 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. (Reference Wagner, Schülein, Petervari, Hannemann, Ali, Cerminara and Sandham2018b) (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. (Reference Wagner, Schülein, Petervari, Hannemann, Ali, Cerminara and Sandham2018b) 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 Reference White and Corfield2006)

(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 Reference Lees1956)

(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 Reference Franko and Lele2013)

(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 (Reference White and Corfield2006) 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 Reference Softley1969) 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. (Reference Wagner, Wartemann, Dittert and Kütemeyer2019), 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. Reference Sousa, Wartemann, Wagner and Scalo2023) 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. (Reference Sousa, Wartemann, Wagner and Scalo2023). 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. Reference Sousa, Patel, Chapelier, Wartemann, Wagner and Scalo2019). 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. (Reference Sousa, Wartemann, Wagner and Scalo2023) 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. (Reference Wagner, Wartemann, Dittert and Kütemeyer2019). 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. (Reference Wagner, Wartemann, Dittert and Kütemeyer2019) 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. (Reference Sousa, Wartemann, Wagner and Scalo2023).

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. (Reference Wagner, Wartemann, Dittert and Kütemeyer2019), 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. (Reference Wagner, Wartemann, Dittert and Kütemeyer2019), 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. (Reference Chynoweth, Ward, Henderson, Moraru, Greenwood, Abney and Schneider2014) on a flared cone via a temperature sensitive paint technique, and by Gray & Schneider (Reference Gray and Schneider2020) 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 Reference Hader and Fasel2019) and straight slender (Hader, Leinemann & Fasel Reference Hader, Leinemann and Fasel2020) 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 (Reference Hader and Fasel2018) 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 (Reference Chen and Scalo2021) 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. (Reference Wagner, Wartemann, Dittert and Kütemeyer2019), 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. (Reference Wartemann, Wagner, Kuhn, Eggers and Hannemann2015) 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. (Reference Wartemann, Wagner, Giese, Eggers and Hannemann2014) when comparing the results obtained from the DLR stability code, NOLOT (Hein et al. Reference Hein, Bertolotti, Simen, Hanifi and Henningson1994; Wartemann, Lüdeke & Sandham Reference Wartemann, Lüdeke and Sandham2012), with the same experimental results. Moreover, in the study performed by Wartemann et al. (Reference Wartemann, Wagner, Giese, Eggers and Hannemann2014), 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. (Reference Alba, Casper, Beresh and Schneider2010) at the BAM6QT comparing experimental results and the stability theory predictions by STABL (Johnson & Candler Reference Johnson and Candler2005) also reported inaccuracies that ranged from 2 % to 14 %. In an effort to reduce the frequency shift, Wartemann et al. (Reference Wartemann, Wagner, Wagnild, Pinna, Miró Miró, Tanno and Johnson2019) 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. (Reference Sousa, Wartemann, Wagner and Scalo2022).

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 (Reference van Driest1951)

(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 Reference Trettel and Larsson2016; Zhang, Duan & Choudhari Reference Zhang, Duan and Choudhari2018). 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 (Reference Trettel and Larsson2016) (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. (Reference Zhang, Duan and Choudhari2018) 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 (Reference Trettel and Larsson2016)

(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 Reference Morkovin1962; Huang, Coleman & Bradshaw Reference Huang, Coleman and Bradshaw1995) 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. (Reference Zhang, Duan and Choudhari2018) 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 (Reference Trettel and Larsson2016) 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 Reference Duan, Beekman and Martin2010). 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. Reference Duan, Beekman and Martin2010) as well as for higher Mach numbers at adiabatic wall conditions (Duan, Beekman & Martin Reference Duan, Beekman and Martin2011). 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. Reference Duan, Beekman and Martin2010).

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 (Reference Franko and Lele2013), 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.