1. Introduction
Boundary-layer transition has always been a critical topic of study in the field of fluid mechanics, especially for the hypersonic flow over aerospace vehicles. Since the turbulent boundary layer can cause higher aerodynamic drag and increased heating on hypersonic vehicles, an accurate prediction of the laminar–turbulent transition location is crucial for vehicle design in the placement of thermal protection. Morkovin (Reference Morkovin1994) first characterised the boundary-layer laminar–turbulent transition into various paths depending on the disturbance amplitude. Figure 1 is adopted from Zhong & Wang (Reference Zhong and Wang2012) to show the various paths of boundary-layer transition. For all paths, the transition is initiated by the initial perturbations generated from ambient forcing disturbances through the receptivity mechanism. Depending on the amplitude of the disturbance, the instability waves can experience a linear modal growth, a transient growth due to non-orthogonal modal (non-modal) interactions, or bypass to turbulence. In order to control and suppress the transition to turbulence, understanding the transition mechanism and its underlying pathways is paramount for controlling turbulence onset.
Path to transition with respect to the disturbance amplitude, adopted from Zhong & Wang (Reference Zhong and Wang2012).
This work focuses primarily on the small-perturbation regime, which closely resembles flight environments. In this regime, the disturbances from the free stream enter the boundary layer and give rise to initial perturbations during the receptivity process. These weak initial disturbances experience linear growth, leading to nonlinear-mode interactions and eventual breakdown. To predict the transition location in path A in figure 1, current prediction procedures include methods such as the
$\textrm{e}^N$
method and the amplitude method by Mack (Reference Mack1977). While the widely applied
$\textrm{e}^N$
method can be easily implemented, the method does not account for broadband effects and the receptivity mechanism at work. The amplitude method, proposed by Mack, provides a more accurate prediction of boundary-layer transition by including receptivity, nonlinear and broadband effects. Integrating across the wavenumber and frequency spectra, this method directly estimates the maximum disturbance amplitude at a given position of the boundary layer with the initial amplitude. Although straightforward in theory, Mack’s method can pose challenges on the implementation without assumptions, such as the infinite frequency spectrum and the availability of receptivity data. Marineau (Reference Marineau2017) addressed these challenges by making practical assumptions, such as approximating the integral across the finite-bandwidth-disturbance spectra as a summation, and with experimental correlations of receptivity coefficients. Given the critical role that receptivity data play in the amplitude method, an accurate evaluation of receptivity coefficients is essential. To address this need, the current study aims to develop and demonstrate a combined computational and theoretical analysis tool capable of extracting receptivity data directly.
In this work, the process of external disturbances entering the flow and exciting instability, defined as the free-stream receptivity by Zhong & Wang (Reference Zhong and Wang2012), is focused on. Extensive studies including theoretical, experimental and numerical approaches have been conducted in the investigation of free-stream receptivity. Goldstein & Hultgren (Reference Goldstein and Hultgren1989), Fedorov & Khokhlov (Reference Fedorov and Khokhlov1991) and Qin & Wu (Reference Qin and Wu2016) have made substantial progress in receptivity studies using the theoretical approach, which includes the application of asymptotic theory alone or in combination with bi-orthogonal decomposition. However, analytical solutions are limited to simple flow geometries and conditions due to the complexity of the analysis. Moreover, the experimental work of Marineau et al. (Reference Marineau, Moraru, Lewis, Norris, Lafferty, Wagnild and Smith2014) investigated a Mach
$10$
boundary-layer flow over sharp and blunted cones using the AEDC Hypervelocity Wind Tunnel 9. Marineau observed, for a large bluntness, no significant second-mode instability was observed at the start of the transition. In this case, the transition might be linked to tunnel noise or non-modal effects. The latter indicates the need for a multimode method in the theoretical analysis for validation.
Recently, significant improvements in numerical simulations, such as direct numerical simulation (DNS), have been found to provide detailed insight into the flow physics and allow for the investigation of complex flow configurations. For hypersonic flows, Zhong (Reference Zhong1998) developed a high-order finite difference (FD) shock-fitting method and later Ma & Zhong (Reference Ma and Zhong2003) applied this method to study the receptivity of free-stream acoustic waves with fixed frequencies for a Mach 4.5 flow over a flat plate. Ma and Zhong identified discrete modes I, II and higher that arise from the fast acoustic spectra, and also observed the presence of the
$1$
st,
$2$
nd and higher modes originating from the slow acoustic spectra. These discrete modes were named the discrete F and S modes accordingly by Fedorov & Tumin (Reference Fedorov and Tumin2011) and this nomenclature will be used throughout the remainder of this chapter. Extending the work, Huang & Zhong (Reference Huang and Zhong2014) studied the broadband effects of a free-stream entropy hotspot on receptivity over a blunt cone and compared the simulated results with Linear Stability Theory (LST). Similar to the DNS result from Ma & Zhong (Reference Ma and Zhong2003), the broadband receptivity results exhibit an oscillating behaviour of the phase speed near the fast acoustic phase speed indicating a significant contribution from the fast acoustic wave in the region near the synchronisation point. Huang and Zhong proposed an approach for obtaining the receptivity coefficient of the second mode by approximating the initial amplitude of the dominant mode with the
$N$
-factor growth rate. This procedure, referred to as the backtracking method in the rest of this study, divides the second-mode amplitude from a dominant location by the LST
$N$
factor. The receptivity coefficient or the initial amplitude is then defined as the ratio between the disturbance amplitude,
$A (s^*, f_n )$
, obtained from the unsteady DNS data, and the
$e^N$
value at a particular streamwise location
$s^*$
and frequency
$f_n$
:
\begin{equation} C_{\textit{rec}}\left (f_n\right )=A_0\left (f_n\right )=\frac {A\left (s^*, f_n\right )}{e^{N\left (s^*, f_n\right )}}. \end{equation}
Here, the non-dimensional amplitude spectral density,
$A (s^*, f_n )$
, is normalised by the free-stream pulse density perturbation. This transcribes to a normalised
$C_{\textit{rec}}$
as well. This definition is similar to the receptivity coefficient used in the work of Marineau (Reference Marineau2017) and can be applied directly in the amplitude method. He & Zhong (Reference He and Zhong2021a
, Reference He and Zhong2021b
, Reference He and Zhong2022) expanded further on the research of broadband receptivity. Recently, He & Zhong (Reference He and Zhong2023) applied their receptivity results to the amplitude method implemented by Marineau (Reference Marineau2017) to improve the accuracy of the transition-prediction method. With the numerical results, He & Zhong (Reference He and Zhong2023) was able to generate the receptivity coefficients for various flow conditions and scenarios. However, the backtracking method only focuses on the contribution of the unstable mode and neglects other modes, such as the continuous- and entropy-layer modes. The bi-orthogonal decomposition, although also relying on the LST assumptions, considers entropy-layer and non-dominant modes that would effectively enhance the understanding of receptivity. Furthermore, the backtracking approach applies the
$\textrm{e}^N$
method directly to the DNS flow field, which includes contributions from both modal disturbances and free-stream forcing. As a result, the computed receptivity coefficient inherently contains non-modal components associated with the initial disturbance pulse. In contrast, the present study employs bi-orthogonal decomposition to isolate the true modal amplitudes, thereby reducing the influence from free-stream forcing.
The bi-orthogonal decomposition has been proposed as an enhancement to asymptotic and numerical methods to solve the receptivity problem by Huang & Zhong (Reference Huang and Zhong2014) and He (Reference He2022). Proven in the works of Salwen & Grosch (Reference Salwen and Grosch1981) and Tumin & Fedorov (Reference Tumin and Fedorov1983), the solutions to both temporal and spatial problems of linear growth in the boundary layer include various discrete and continuous normal modes. Another set of solutions, which are adjoint and orthogonal to the original problem, exists and forms a bi-orthogonal eigenfunction system (BES) with the original solution. The bi-orthogonal decomposition involves applying adjoint fields, which are orthogonal to the original solution, to obtain the amplitude coefficients for each modal solution. In this context, the discrete modes consist of two wave types, mode F and mode S, as defined by Fedorov & Tumin (Reference Fedorov and Tumin2011). In addition, the continuous modes are referred to as the non-decaying entropy, vorticity and fast/slow acoustic waves. These continuous modes originate from the boundary layer and extend to the outer shock layer. Furthermore, the intermodal energy exchange, known as synchronisation, between discrete and continuous modes plays a significant role in the transition. More specifically, the energy exchange between the discrete first mode and second/higher-order modes leads to the well-known second-mode instability in both computational and theoretical studies by Fedorov & Khokhlov (Reference Fedorov and Khokhlov2001) and Ma & Zhong (Reference Ma and Zhong2003). Moreover, Tumin (Reference Tumin2011) also presented a survey of applications to the receptivity problem in computational and experimental work. Saikia, Al Hasnine & Brehm (Reference Saikia, Al Hasnine and Brehm2022) used the BES to examine the amplification of the supersonic discrete mode of a high-enthalpy flow over a flat plate. Hasnine et al. (Reference Al Hasnine, Russo, Browne, Tumin and Brehm2020) studied a particulate-induced disturbance over a plate boundary layer of a high-speed flow. Furthermore, Miselis, Huang & Zhong (Reference Miselis, Huang and Zhong2016) performed the decomposition for an entropy-spot perturbed flow field over a blunt cone. However, there is only limited literature on the application of the method to flow over blunt geometries. In the current study, we aim to deploy the bi-orthogonal decomposition more extensively on a hypersonic flow over a blunt cone to investigate the boundary-layer and entropy-layer disturbances in depth.
Furthermore, input–output, or resolvent, analysis is another popular tool for transition analysis due to its capability of considering non-modal effects. In the input–output approach, an arbitrary forcing, including the nonlinear term from the full nonlinear Navier–Stokes equation or an external forcing on the linearised Navier–Stokes equations acts as the input to the discrete operator and the response of the forcing is studied. To use this tool, the Navier–Stokes equations for the full domain are discretised, using the spectral collocation method by McKeon & Sharma (Reference McKeon and Sharma2010), Jovanović & Bamieh (Reference Jovanović and Bamieh2005), the finite volume method by Dwivedi, Sidharth & Jovanović (Reference Dwivedi, Sidharth and Jovanović2021) or a numerical Jacobian method by extracting the linear operator from known solvers in the works of Cook et al. (Reference Cook, Thome, Brock, Nichols and Candler2018) and Nichols & Candler (Reference Nichols and Candler2019). Although this analysis accounts for the non-modal effects, the resolvent analysis decomposes disturbances into modes with respect to the energy norm presented by Chu (Reference Chu1965), which contains both modal and non-modal content. A recent study by Cook & Nichols (Reference Cook and Nichols2024) has presented a comparison between the LST mode F and a direct response mode using the resolvent analysis, showing that such a response mode contains both mode F in the boundary later and non-modal content in the entropy layer. Input–output analysis highlights the early, non-modal amplification of external disturbances, whereas bi-orthogonal decomposition provides a clear modal decomposition of the flow field and isolates each mode’s amplitude evolution. Thus, the input–output analysis and modal analysis, such as the bi-orthogonal decomposition, can complement each other and be used to better understand the source of receptivity, whether modal or non-modal.
Although the receptivity mechanism has been extensively studied with theoretical, experimental and DNS approaches, direct evaluation of receptivity coefficients at the branch I location for hypersonic flows over blunt cones remains very limited. In addition, the previous backtracking method incurs significant computational costs because it requires simulations to be performed until a second-mode dominant region is reached. To address current limitations, the primary goal of this paper is to present a practical tool employing bi-orthogonal eigenfunction decomposition for directly and accurately evaluating receptivity coefficients at the branch I location. Implementing the high-order finite difference method used by Zou & Zhong (Reference Zou and Zhong2023), the BES including discrete boundary-layer and entropy-layer modes is obtained for a hypersonic flow over a blunt cone. The preliminary decomposed result of the receptivity flow field has shown agreement between the decomposed mode S amplitude and the backtracked mode S amplitude based on the
$\textrm{e}^N$
method. With the initial amplitudes obtained with the bi-orthogonal decomposition, the receptivity coefficients have been computed directly for the first time to enhance the amplitude method in transition prediction. Furthermore, an entropy-layer discrete mode, existing only in flow over blunt cone geometries, is incorporated into the bi-orthogonal decomposition and its contribution to the flow field is explored.
2. Governing equations
2.1. Navier–Stokes equations
While hypersonic flows are typically characterised by high enthalpy and chemical non-equilibrium effects, as discussed in Mortensen & Zhong (Reference Mortensen and Zhong2016), many experimental investigations, including Marineau et al. (Reference Marineau, Moraru, Lewis, Norris, Lafferty, Wagnild and Smith2014, Reference Marineau, Moraru, Lewis, Norris, Lafferty and Johnson2015), are conducted under low-enthalpy conditions. To align with these experiments and improve the accuracy of receptivity data, the current study adopts the low-enthalpy free-stream conditions with the calorically perfect gas assumption. For the compressible boundary-layer flow, all flow variables are non-dimensionalised by their corresponding boundary-layer-edge values, velocities
$\boldsymbol{u}$
by the edge velocity
$u_e$
, pressure
$p$
by the edge dynamic pressure
$\rho _e u_e^2$
, density
$\rho$
density by the edge density
$\rho _e$
. A reference length scale
$L$
is assumed for all lengths and the time variable,
$t$
, is scaled by
$L/u_e$
. Thus, the dimensionless Navier–Stokes equations for a viscous compressible flow of ideal gas in Cartesian coordinates are used:
\begin{align}&\qquad\qquad\qquad\qquad\qquad\quad \frac {\partial \rho }{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{u}) = 0, \end{align}
\begin{align}& \rho \left ( \frac {\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} \right ) = -\boldsymbol{\nabla }p + \frac {1}{\textit{Re}} \boldsymbol{\nabla }\boldsymbol{\cdot }\left [ \lambda (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}) \boldsymbol{I} + \mu \left ( \boldsymbol{\nabla }\boldsymbol{u} + (\boldsymbol{\nabla }\boldsymbol{u})^\textrm{T} \right ) \right ]\!, \end{align}
\begin{align}&\qquad\qquad \rho c_p \left ( \frac {\partial T}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }T \right ) = \frac {1}{{\textit{Re}}\boldsymbol{\cdot } {\textit{Pr}}} \boldsymbol{\nabla }\boldsymbol{\cdot }(\kappa \boldsymbol{\nabla }T) + \frac {1}{{\textit{Re}}} \varPhi , \end{align}
where
$Re$
is the Reynold number and
$Pr$
is the Prandtl number. The dimensionless viscous dissipation function is
\begin{equation} \varPhi = \lambda (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u})^2 + \frac {\mu }{2} \left [ \boldsymbol{\nabla }\boldsymbol{u} + (\boldsymbol{\nabla }\boldsymbol{u})^\textrm{T} \right ]^2. \end{equation}
The ideal gas law is given by
\begin{equation} p = \frac {\rho T}{\gamma M_\infty ^2}. \end{equation}
We also make use of the Stokes’ hypothesis of
$\lambda ^{*} =-2 \mu ^{*} / 3 .$
The viscosity
$\mu$
is given by the Sutherland formula with the reference air temperature set to
$300.55\, \rm K$
. The thermal conductivity
$k$
may also be formulated by the Sutherland formula with the Sutherland temperature
$S_k$
. For the result of this paper, however,
$k$
is computed by assuming a Prandtl number
$Pr = 0.72$
, defined as
$\mu c_{p}/k$
where
$c_p$
is the specific heat at constant pressure and assumed to be constant for a perfect gas.
2.2. Linear stability theory
To study the property of instability waves in the boundary-layer flow, a small perturbation is assumed. From the Navier–Stokes equations, we can derive the linearised Navier–Stokes equations (LNSE) by assuming a small disturbance added to the flow variables following previous works of Malik (Reference Malik1990) and Tumin (Reference Tumin2007). The quasi-parallel assumption, which assumes that the mean variable depends only on the wall-normal direction, i.e. Mean-flow velocities,
$U \approx U(y)$
and
$V = 0$
, is applied. By neglecting the higher-order terms such as the disturbance velocity squared,
$\tilde{u}^2$
, and subtracting the mean-flow equations, we can obtain the non-dimensional linear stability equations. For a flat-plate profile, the assumption of
${\rm d} P / {\rm d} y \rightarrow 0$
causes the pressure profile to be constant,
$P=1/(\gamma M^2)$
. For mean-flow profiles of other geometry, as suggested by Miselis et al. (Reference Miselis, Huang and Zhong2016), the assumption is not applicable and the mean-flow pressure distribution has to be accounted for in the linearised equations. The mean-flow velocity boundary condition at the wall,
$U_{wall}$
, will be the no-slip condition. The mean-flow temperature boundary condition,
$T_{wall}$
, can be adiabatic or isothermal at the wall:
\begin{equation} \begin{gathered} U_{w\textit{all}} = 0\\ \frac {{\rm d} T_{w\textit{all}}}{{\rm d} y} = 0 \; \text{or} \; T_{w\textit{all}} = \text{constant}. \end{gathered} \end{equation}
For the fluctuations velocities at the wall,
$u_{wall}$
,
$v_{wall}$
and
$w_{wall}$
, will still have the no-slip condition. However, the fluctuation temperature boundary condition,
$\theta_{wall}$
, is set to be zero because of the thermal inertia of the wall:
\begin{equation} \begin{gathered} \tilde {u}_{w\textit{all}},\tilde {v}_{w\textit{all}},\tilde {w}_{w\textit{all}} = 0\\ {\tilde{\theta}_{wall}} = 0. \end{gathered} \end{equation}
In the free stream, the flow parameters are governed by the shock using the Rankine–Hugoniot condition in the DNS of Zhong (Reference Zhong1998). The free-stream perturbations, including velocities
$u,v,w$
, pressure
$p$
, and temperature
$\theta$
are defined to be bounded by Salwen & Grosch (Reference Salwen and Grosch1981), allowing the existence of neutral continuous modes,
\begin{equation} y \rightarrow \infty : \big |\tilde {u},\tilde {v},\tilde {w}, \tilde {p}, \tilde {\theta } \big |\lt \infty , \end{equation}
as well as the corresponding spatial derivatives in the streamwise and spanwise directions.
From the set of linear stability equations, a periodic-in-time perturbation is assumed, which leads to a solution in the form of
$\exp (-\textrm{i}\omega t)$
after a Fourier transformation in time is performed. After the transformation, the perturbation variables and their corresponding spatial derivatives are defined in a column vector,
$\boldsymbol{A}$
, with
$16$
components:
\begin{equation} \begin{gathered} \boldsymbol{A}(x, y, z)= (u, \partial u / \partial y, v, p, \theta , \partial \theta / \partial y, w, \partial w / \partial y, \partial u / \partial x, \partial v / \partial x,\\ \partial \theta / \partial x, \partial w / \partial x, \partial u / \partial z, \partial v / \partial z, \partial \theta / \partial z, \partial w / \partial z)^{T}. \end{gathered} \end{equation}
The
$\hat {u}$
notation is dropped and written as
$u$
for simplicity in notation. Vector
$\boldsymbol{A}$
is defined with
$16$
variables such that the LNSE can be written to be linear with respect to the streamwise wavenumber
$\alpha$
following Tumin (Reference Tumin2007). With the
$16$
-variable form, the orthogonality relation will be defined in a later section. In this study, the perturbation is assumed to take the form of a travelling wave, such as the following streamwise velocity disturbance variable,
$\tilde {u}$
:
\begin{equation} \tilde {u}=\hat {u}(y) \exp [\textrm{i}( \alpha x+ \beta z-\omega t)] \end{equation}
where
$\alpha$
and
$\beta$
are the streamwise and spanwise wavenumbers,
$c$
is the wave propagation speed, and
$\omega = \alpha c$
is the wave frequency. The wave shapefunction is denoted by
$\hat {u}(y)$
. In the temporal problem, both wavenumbers,
$\alpha$
and
$\beta$
are assumed real and the complex frequency
$\omega$
is determined. On the other hand, for the spatial problem, the spanwise wavenumber
$\beta$
and the frequency
$\omega$
are real valued and the streamwise wavenumber
$\alpha$
is complex. For discrete modes, the streamwise wavenumber
$\alpha$
is solved as a complex eigenvalue. For continuous modes, however, the streamwise wavenumber
$\alpha$
is determined from the branch cuts in the solution of the Cauchy problem presented in a later section. Substituting the solution form of (2.10) into the LNSE and eliminating
$ \exp [\textrm{i}( \alpha x+ \beta z-\omega t)]$
on both sides of the equation, the system with
$\boldsymbol{A}_{\alpha \beta }$
being the transformed variable of
$\boldsymbol{A}$
in (2.9) can be written as
\begin{equation} \begin{gathered} \frac {\partial }{\partial y}\left (\boldsymbol{L}_{0} \frac {\partial \boldsymbol{A}_{\alpha \beta }}{\partial y}\right )+\boldsymbol{L}_{1} \frac {\partial \boldsymbol{A}_{\alpha \beta }}{\partial y}=\boldsymbol{H}_{1} \boldsymbol{A}_{\alpha \beta }+{\rm{i}} \alpha \boldsymbol{H}_{2} \boldsymbol{A}_{\alpha \beta }+{\rm{i}} \beta \boldsymbol{H}_{3} \boldsymbol{A}_{\alpha \beta } \\ y=0: \quad (\boldsymbol{Z_{\alpha \beta }})_i=0, \quad i=1,3,5,7 \\ y \rightarrow \infty : \quad \left |(\boldsymbol{A_{\alpha \beta }})_i\right |\lt \infty , \quad i=1, \ldots , 16. \end{gathered} \end{equation}
The elements of
$\boldsymbol{L}_0$
,
$\boldsymbol{L}_1$
,
$\boldsymbol{H}_1$
,
$\boldsymbol{H}_2$
and
$\boldsymbol{H}_3$
are given in Appendix A. For implementation, (2.11) can be recast as a first-order system of eight components of the following with the corresponding boundary conditions:
\begin{equation} \begin{gathered} \frac {{\rm d} \boldsymbol{z}_{\alpha \beta }}{{\rm d} y}=\boldsymbol{H} \boldsymbol{z}_{\alpha \beta }\\ y=0: (\boldsymbol{z_{\alpha \beta }})_i=0, \quad i=1,3,5,7 \\ y \rightarrow \infty : \left |(\boldsymbol{z_{\alpha \beta }})_i\right |\lt \infty , \quad i=1, \ldots , 8, \end{gathered} \end{equation}
with
$\boldsymbol{z_{\alpha \beta }}$
defined as
\begin{equation} \boldsymbol{z_{\alpha \beta }} = (u, \partial u / \partial y, v, p, \theta , \partial \theta / \partial y, w, \partial w / \partial y). \end{equation}
Here,
$\boldsymbol{H}$
is an
$8$
by
$8$
coefficient matrix. The elements of
$\boldsymbol{H}$
can be found in Appendix A. The solutions of (2.12) can be classified into two types, discrete and continuous modes, according to the free-stream boundary conditions. From (2.12), a total of eight fundamental solutions can be obtained with each fundamental solution
$\boldsymbol{z_{\alpha \beta }^{(m)}}(y)$
,
$m = 1,{\cdots} ,8$
, where
$m$
is from 1 to 8 for the 8 fundamental solutions, satisfying the following:
\begin{equation} \begin{gathered} \frac {{\rm d} \boldsymbol{z^{(m)}_{\alpha \beta }}}{{\rm d} y}=\boldsymbol{H}\boldsymbol{z^{(m)}_{\alpha \beta }}(y), \quad m =1,\;{\cdots} ,\; 8,\\ y \rightarrow \infty : \boldsymbol{z^{(m)}_{\alpha \beta }}(y) = \boldsymbol{z^{(m)}_{\alpha \beta _0}}\textrm{e}^{\lambda _m y}, \quad m =1,\;{\cdots} ,\; 8, \end{gathered} \end{equation}
where
$ \boldsymbol{z_{\alpha \beta _0}}^{(m)}$
is the
$m$
th fundamental solution in the free stream. The derivation of eigenvalues
$\lambda _m$
and the corresponding asymptotic solutions
$\boldsymbol{z^{(m)}_{\alpha \beta _0}}$
can be found in previous works of Tumin (Reference Tumin2007) and Zou & Zhong (Reference Zou and Zhong2024).
2.2.1. Discrete modes
Four solutions having the asymptotic behaviour of
$\sim \exp {(\lambda _m y)}, \; m = 2,\;4,\;6,\;8$
where
${\textit{Re}}(\lambda _m) \gt 0$
and the other four solutions being
$\sim \exp {(\lambda _m y)}, \; m = 1,\;3,\;5,\;7$
where
${\textit{Re}}(\lambda _m) \lt 0$
. For discrete modes, the boundary condition at the free stream can be formulated as a decaying condition as a subset of the bounded condition in (2.12). These discrete modes include free-stream decaying modes such as mode F, originating from the fast acoustic spectrum, mode S, originating from the slow acoustic spectrum and higher-order modes. Based on the decaying nature of these fundamental solutions in free stream, such boundary conditions can also be written as
\begin{equation} y \rightarrow \infty : (\boldsymbol{z_{\alpha \beta }})_{k}=0, \quad k=1, 3,4,5,7. \end{equation}
These free-stream boundary conditions can replace the free-stream boundary conditions in (2.12) and form a homogeneous problem. For implementation, the homogeneous problem can be formulated as a nonlinear boundary-value problem given in Zou & Zhong (Reference Zou and Zhong2024).
For the case with a bow shock at the end of the domain, a set of shock boundary conditions based on the linearised Rankine–Hugoniot relations can also be implemented following Chang, Malik & Hussaini (Reference Chang, Malik and Hussaini1990) and Knisely & Zhong (Reference Knisely and Zhong2019a ). These are also the jump conditions implemented by Cook & Nichols (Reference Cook and Nichols2022) in their resolvent analysis. However, these linearised Rankine–Hugoniot relations were concluded to have negligible influence on the direct LST result by Knisely & Zhong (Reference Knisely and Zhong2019a ). Thus, this study will impose the decaying conditions as Tumin (Reference Tumin2007).
2.2.2. Continuous spectra
In addition to the discrete modes, the boundedness condition allows the existence of neutral waves oscillating beyond the boundary layer known as continuous modes. According to the work of Tumin, Wang & Zhong (Reference Tumin, Wang and Zhong2007), the continuous acoustic modes contribute mainly to the outside of the boundary layer, apart from the free-stream forcing. Furthermore, for the flow conditions considered in the current study, previous work by He & Zhong (Reference He and Zhong2022) has shown that the phase speed of instability at the branch I location is far from the fast and slow acoustic spectrum. Thus, the current study will focus on the discrete-mode receptivity.
Nevertheless, the continuous mode formulation can be found in the previous work of Zou & Zhong (Reference Zou and Zhong2024) for future interest. The two types of solutions combine to form a complete set of solutions as a superposition of discrete and continuous eigenfunctions accordingly:
\begin{equation} \boldsymbol{A}=\sum _n C_n \boldsymbol{A}_{\alpha _n}(y) \textrm{e}^{{\rm{i}} \alpha _n x}+\sum _j \int _0^{\infty } C_j \boldsymbol{A}_{\alpha _j}(y) \textrm{e}^{{\rm{i}} \alpha _j(k) x} {\rm d} k, \end{equation}
with
$A_{{\alpha}_{n}}$
representing the discrete modes and
$A_{{\alpha}_{j}}$
being the continuous modes, containing a superposition of the fundamental solution. Since any perturbation can be decomposed into discrete and continuous modes, the bi-orthogonal decomposition procedure can be formulated. The numerical implementation to obtain both discrete and continuous modes using the finite difference method will be discussed in § 3.
2.3. Adjoint problem and BES
With the 16-variable direct system of (2.11) defined, the corresponding adjoint system with complex-conjugate adjoint variables
$\boldsymbol{B}_{\alpha \beta }$
can be introduced:
\begin{equation} \begin{gathered} \frac {\partial }{\partial y}\left (\boldsymbol{L}_{0}^{\textrm {T}} \frac {\partial \boldsymbol{B}_{\alpha \beta }}{\partial y}\right )-\boldsymbol{L}_{1}^{\textrm{T}} \frac {\partial \boldsymbol{B}_{\alpha \beta }}{\partial y}=\boldsymbol{H}_{1}^{\textrm{T}} \boldsymbol{B}_{\alpha \beta }+{\rm{i}} \alpha \boldsymbol{H}_{2}^{\textrm{T}} \boldsymbol{B}_{\alpha \beta }+{\rm{i}} \beta \boldsymbol{H}_{3}^{\textrm{T}} \boldsymbol{B}_{\alpha \beta } \\ y=0: \quad (B_{\alpha \beta })_i=0, \quad i=2,4,6,8 \\ y \rightarrow \infty : \quad \left |(B_{\alpha \beta })_i\right |\lt \infty , \quad i=1, \ldots , 16. \end{gathered} \end{equation}
Here, the transpose is the conventional transpose without conjugation since the complex conjugate is introduced in the definition of
$\boldsymbol{B}_{\alpha \beta }$
. The boundary conditions for the adjoint problems are also obtained from the inner product definition with the direct-problem boundary conditions. The adjoint system can also be recast in the form of first-order equations, with the adjoint variable vector
$\boldsymbol{Y}_{\alpha \beta }$
presented in the works of Tumin (Reference Tumin2007). With the dual systems, the orthogonality condition exists between two different modes with eigenvalues
$\alpha$
and
$\alpha'$
\begin{equation} {\rm{i}}\left (\alpha -\alpha ^{\prime }\right ) \int _{0}^{\infty }\left (\boldsymbol{H}_{2} \boldsymbol{A}_{\alpha \beta }, \boldsymbol{B}_{\alpha ^{\prime } \beta }\right ) {\rm d} y=0. \end{equation}
This orthogonality relation can be obtained by an integration by part of the dot product between
$\boldsymbol{B}_{\alpha ^{\prime } \beta }$
and (2.11). Equation (2.18) can be rewritten into this form:
\begin{equation} \left \langle \boldsymbol{H}_{2} \boldsymbol{A}_{\alpha \beta }, \boldsymbol{B}_{\alpha ^{\prime } \beta }\right \rangle =\int _{0}^{\infty }\left (\left (\boldsymbol{H}_{2} \boldsymbol{A}_{\alpha \beta }\right ), \boldsymbol{B}_{\alpha ^{\prime } \beta }\right ) {\rm d} y=\varGamma \Delta _{\alpha \alpha ^{\prime }}. \end{equation}
According to Salwen & Grosch (Reference Salwen and Grosch1981), for the discrete modes, the term
$\Delta _{\alpha \alpha ^{\prime }}$
represents a Kronecker Delta, which equals
$1$
if the eigenvalues of the two modes are the same. This is because the decaying behaviour of the discrete modes and the integral should result in a constant
$\varGamma$
value depending on the normalisation. For the continuous modes, the same term represents a Dirac delta. Tumin (Reference Tumin2007) also derived the following inner product relation between
$ \boldsymbol{A}_{\alpha \beta }$
and
$\boldsymbol{z}_{\alpha \beta }$
, which would be helpful in the numerical implementation:
\begin{equation} \left \langle \boldsymbol{H_2} \boldsymbol{A}_{\alpha \beta }, \boldsymbol{B}_{\alpha \beta }\right \rangle =-\mathrm{i}\left \langle \frac {\partial \boldsymbol{H}}{\partial \alpha } z_{\alpha \beta }, \boldsymbol{Y}_{\alpha \beta }\right \rangle\!. \end{equation}
To verify the bi-orthogonality of the modes, the orthogonality relationship is evaluated numerically in a later section. Subsequently, the coefficients for each modes can be obtained from the DNS flowfield
$Z_{\textit{DNS}}$
using the following relation assuming local parallel flow:
\begin{equation} C_{mode}=\frac {-\mathrm{i}\left \langle \frac {\partial \boldsymbol{H}}{\partial \alpha } z_{\textit{DNS}}, \boldsymbol{Y}_{\alpha \beta }\right \rangle }{\varGamma }. \end{equation}
3. Numerical methods
3.1. Direct numerical simulation
A receptivity simulation is performed with a high-order shock-fitting DNS developed by Zhong (Reference Zhong1998) for a hypersonic flow over a cone on a curvilinear grid. The governing equations in vector form are written as
\begin{equation} \frac {\partial U}{\partial t} \;+\; \frac {\partial F_{\!j}}{\partial x_{\textit{ij}}} \;+\; \frac {\partial G_{\!j}}{\partial x_{\!j}} \;=\;0, \qquad (j=1,2,3), \end{equation}
where
$U$
is the state vector of conserved quantities and
$F_j$
and
$G_j$
are the inviscid and viscous flux vectors, respectively. Here, the
$j$
indices indicate curvilinear coordinates in the streamwise, radial and azimuthal directions about the cone. The conservative vector
$U$
comprises five conservative flow variables for mass, momentum and energy. The shock is treated as a moving boundary and the flow condition behind the shock is determined with the Rankine–Hugoniot relations. The numerical scheme used is a fifth-order upwind scheme for the inviscid fluxes and a sixth-order central scheme for the viscous fluxes. Additional details of the general numerical method used for the DNS and the shock-fitting scheme can be found in the work of Zhong (Reference Zhong1998). This shock-fitting scheme is used for high-order accuracy and the ability to resolve the flow field after the shock. For steady mean flow, a parallel assumption is applied and the mean-flow variables are functions of the wall-normal direction coordinate only. Under this assumption, the flow field data are divided into wall-normal snapshots for each streamwise location and each snapshot at these streamwise sampling locations is analysed.
3.2. Free-stream disturbance model
Following the previous DNS work by He & Zhong (Reference He and Zhong2021a , Reference He and Zhong2022), the perturbed flow field is obtained by imposing a propagating Gaussian distributed free-stream disturbance throughout the domain:
\begin{equation} q(x, y, z, t)=\left |q^{\prime }\right |_{\infty } \exp \left (-\frac {\left (R_c\right )^2}{2 \sigma ^2}\right )+q_{\infty }. \end{equation}
The term
$q$
is the disturbance variable and
$\big |q^{\prime }\big |_{\infty }$
stands for the peak free-stream perturbation. The term
$R_c$
refers to the radial distance from the centre of the pulse to a point in the flow field. The variable
$\sigma$
controls both the spatial width of the pulse as well as the frequency bandwidth of the disturbance since
$R_c$
is defined as a function of
$x-x_{pulse}$
where
$x_{pulse} = x_0 + C_\infty t$
with
$x_0$
to be the starting location and
$C_{\infty}t$
to be the phase speed times time. For the detailed formulation and parameters of the current case, one can refer to case I from the previous study of He & Zhong (Reference He and Zhong2022).
To study the receptivity behaviour to various disturbances, the standard procedure is to decompose the perturbed flow field resulting from each pulse into their spectral frequency components, deployed in various works of Sivasubramanian & Fasel (Reference Sivasubramanian and Fasel2014) and A. Al Hasnine et al. (Reference Al Hasnine, Russo, Tumin and Brehm2023). The unsteady DNS flow field is analysed with a fast Fourier transform (FFT) to obtain the disturbance information in the discrete frequency spectrum. Although the FFT is ideally suited for statistically stationary or harmonically forced signals, the semi-infinite streamwise computational domain in the current study allows the initial disturbance to propagate without significant reflections. Furthermore, the agreement between the spectral components of the DNS flow field at the second-mode dominant location and the discrete mode S from the BES would also validate the approach. The bi-orthogonal decomposition is then applied to the spectral components of the flow field to further analyse the modal contribution.
3.3. High-order finite difference method for discrete and continuous modes
The general high-order finite difference method on a non-uniform grid by Zou & Zhong (Reference Zou and Zhong2024) is applied to obtain discrete and continuous modes. Using the boundary-value-problem approach, the method can be adapted to compute both discrete and continuous modes using different boundary conditions and eliminates the need for intricate Gram–Schmidt orthonormalisation and the associated recovery process in the initial-value-problem approach. The domain in the wall-normal direction is non-uniform enabling a direct discretisation. The distribution of
$N$
number of grid points over the domain
$[a,b]$
will follow Kosloff & Tal-Ezer (Reference Kosloff and Tal-Ezer1993), Zhong & Tatineni (Reference Zhong and Tatineni2003), and Shukla, Tatineni & Zhong (Reference Shukla, Tatineni and Zhong2007),
\begin{equation} y_{\!j}=\frac {b+a}{2}+(b-a) \frac {\sin ^{-1}(-\alpha _{g} \cos (\pi j / N))}{2 \sin ^{-1} \alpha }, \quad j=0, \ldots , N, \end{equation}
where
$a = 0$
and
$b = y_{\infty }$
. In (3.3),
$\alpha _{g}$
refers to the grid stretching parameter and does not relate to the streamwise wavenumber mentioned in other parts of the paper. The grid stretching parameter,
$\alpha _{g}$
, controls how dense the grid points are distributed near the domain edges, with
$0$
being the Chebyshev distribution and
$1$
being the uniform distribution. This parameter is set to
$0.9995$
to ensure stability following the work of Zhong & Tatineni (Reference Zhong and Tatineni2003). With the grid defined, we apply the finite difference method. In order to approximate the solution at each
$y_i$
, the local stencil of
$l$
points,
$ l= n + m + 1$
where
$n$
being the number of points before point
$i$
and
$m$
being the number of points after
$i$
, can be used. Under this convention, the stencil length
$l$
is used to identify the schemes. Based on this grid, the differential equations of (2.12) can be solved with a finite difference method. By taking the derivative of the Lagrange polynomial, the first derivative for all points
$i = 0, {\cdots} ,N$
can be approximated as
\begin{equation} \frac {{\rm d} \boldsymbol {z_{\alpha \beta }}\left (y_{i}\right )}{{\rm d} y_{i}}=\sum _{j=i - n}^{ i + m} \boldsymbol{z_{\alpha \beta }}\left (y_{\!j}\right ) a_{ji}, \quad i=0, \ldots , N, \end{equation}
where the coefficient
$a_{ji}$
corresponds to the
$j$
th coefficient in the first derivative of the Lagrange polynomial for the point
$y_i$
, given as
\begin{equation} a_{ji}=\sum _{p \neq j} \frac {1}{y_{\!j}-y_{p}} \prod _{q \neq (j, p)} \frac {y_i-y_{q}}{y_{\!j}-y_{q}}. \end{equation}
With the first derivative approximated at each
$y_i$
, the matrix
$\boldsymbol{H}$
of (2.12) also varies for each
$y_i$
location from
$i = 0$
to
$N$
. Hence, the approximated (2.12) can be written as
\begin{equation} \frac {{\rm d} \boldsymbol {z_{\alpha \beta }}\left (y_{i}\right )}{{\rm d} y_{i}}=\sum _{j=i - n}^{ i + m} \boldsymbol{z_{\alpha \beta }}\left (y_{\!j}\right ) a_{ji} = \boldsymbol{H}_i\boldsymbol{z_{\alpha \beta }}\left (y_{i}\right ), \quad i=0, {\cdots} , N. \end{equation}
For the implementation of boundary conditions, Zou & Zhong (Reference Zou and Zhong2024) presented a detailed formulation and implementation procedure for both discrete and continuous modes. With the discrete system, the sparse lower-upper decomposition from MATLAB is used to solve the linear system. With the discrete and continuous modes computed, the bi-orthogonal decomposition method to obtain receptivity coefficients is described.
3.4. Bi-orthogonal decomposition for receptivity
The bi-orthogonal decomposition procedure follows previous works of Tumin (Reference Tumin2003, Reference Tumin2007). Let us consider a decomposition of computational data
$\boldsymbol{A}_{D}$
with discrete and continuous modes. The components of
$\boldsymbol{A}_{D}$
are the spectral components of the perturbation variables for a given frequency,
$f$
, arranged in the order defined in (2.9). Equation (2.16) of any perturbation can be expressed as follows, where
$P=2$
represents the number of discrete modes (F and S), and
$Q=5$
represents the number of continuous modes (fast and slow acoustic modes, two vorticity modes and entropy mode). For each continuous branch, the integral in (2.16) can be approximated with a Riemann sum or other numerical integration scheme. The discrete branch parameter
$K_l, \; l = 1, \;{\cdots} ,\; N$
represents the range of
$K$
for
$N$
discrete points. For every discrete branch parameter of
$K$
, the continuous mode coefficients,
$C_j(K)$
, need to be computed and summed to account for the contribution of the entire branch. Thus, following (2.16), the numerical result from DNS can be approximated as
\begin{equation} \boldsymbol{A}_{D}(x,y,f)\approx \sum _n^{P} C_n \boldsymbol{A}_{\alpha \beta _n}(y) e^{i \alpha _n x}+\sum _j^{Q} \sum _l^N C_j(K_l) \boldsymbol{A}_{\alpha \beta _j}(y,K_l) e^{i \alpha _j(K_l) x} \Delta K. \end{equation}
where
$C_n$
corresponds the modal coefficient of the
$n$
th discrete mode and
$\boldsymbol{A}_{\alpha \beta _n}(y)$
being the discrete-mode eigenfunction. The coefficient
$C_j$
refers to the coefficient of the
$j$
th continuous mode and
$\boldsymbol{A}_{\alpha \beta _j}(y)$
is the corresponding wave shapefunction. Let the vector
$\boldsymbol{z}_{D}$
be the first eight components of
$\boldsymbol{A}_{D}$
; the coefficients for each mode can be obtained with (2.21). The modal coefficient,
$C_{m}$
, where
$m = n$
for the
$n$
th discrete mode and
$m = j$
for the
$j$
th continuous mode, is as follows:
\begin{equation} C_{m}=\frac {\left \langle \dfrac {\partial \boldsymbol{H}}{\partial \alpha } \boldsymbol{z}_{D}, \boldsymbol{Y}_{\alpha \beta }^{m}\right \rangle }{\left \langle \dfrac {\partial \boldsymbol{H}}{\partial \alpha } \boldsymbol{z}_{\alpha \beta }^{m}, \boldsymbol{Y}_{\alpha \beta }^{m}\right \rangle }, \end{equation}
with
$\boldsymbol{z}_{\alpha \beta }^{m}$
being the direct eigenfunction of the
$n$
th discrete mode or the
$j$
th continuous mode from (2.16) and
$\boldsymbol{Y}_{\alpha \beta }^{m}$
being the adjoint eigenfunction for the corresponding discrete or continuous mode. Using a non-uniform trapezoidal rule, the integration can be numerically implemented. Note that a consistent normalisation across the modes is necessary as the normalisation is dependent on the normalisation boundary conditions already mentioned.
Utilising the bi-orthogonal decomposition already shown, the individual amplitude of perturbation variables can be obtained for each mode and the result can be used to aid receptivity studies, particularly in the computation of receptivity coefficient. Theoretically, the receptivity coefficient is defined as the ratio of initial perturbation amplitude of the unstable mode,
$A_{{I}}$
, where the subscript I corresponds to the data measured at the branch I location of the neutral stability curve given a frequency, over the free-stream disturbance amplitude,
$A_{\infty }$
, in the works of Marineau (Reference Marineau2017) and He & Zhong (Reference He and Zhong2021a
), given by
\begin{equation} C_{\textit{rec}}=\frac {A_{{I}}}{A_{\infty }}, \end{equation}
where
$A_{{I}}$
is the unstable mode amplitude of a perturbation variable, which can be
$u,v,p,\theta$
, in the vector
$\boldsymbol{A}$
from (2.9), for a specific frequency
$\omega$
. This amplitude
$A_{{I}}$
contains both magnitude and phase angle. Most applications using the initial amplitude, such as the computation of receptivity coefficient, only require the magnitude of
$A_{{I}}$
. In DNS studies, the phase angle is also significant to the simulation result. The pressure component at the wall
$p_{{\textit{wall}}}$
, which is the
$4$
th component of
$\boldsymbol{A}$
at
$y=0$
, is used by He & Zhong (Reference He and Zhong2021a
) as
$A_{{I}}$
in their cases while the free-stream pulse density perturbation of the same frequency is chosen for
$A_{\infty }$
in (3.9) to have an uniform normalisation over various types of free-stream disturbances such as acoustic and entropy pulses that were considered in the literature.
To obtain the initial amplitudes of the unstable mode at the branch I location, the bi-orthogonal decomposition described is applied. Using the
$m$
th-mode amplitude coefficient of
$C_{m}$
from (2.21) and the solution of the corresponding mode,
$\boldsymbol{z}_{\alpha \beta }^{m}$
or
$\boldsymbol{A}_{\alpha \beta }^{m}$
, evaluated at the branch I location, the receptivity coefficient of the
$m$
th discrete or continuous mode for the
$m$
th mode can be formulated as
\begin{equation} C_{\textit{rec}}=\frac {C_{m} \boldsymbol{A}_{\alpha \beta }^{m}(y=0)}{A_{\infty }} =\frac {C_{m} \boldsymbol{z}_{\alpha \beta }^{m}( y=0)}{A_{\infty }}. \end{equation}
Note that the bi-orthogonal eigenfunctions formulated are normalised by the wall pressure. Thus,
$C_{m} \boldsymbol{A}_{\alpha \beta }^{m}(y=0)$
is multiplied by the wall-pressure amplitude obtained by DNS to recover the absolute amplitude at the wall. The aforesaid presents the procedure of obtaining the receptivity coefficient of an unstable mode by bi-orthogonal decomposition. A case of receptivity analysis for a hypersonic flow over a blunt nose cone will be discussed in a later section.
4. Results
4.1. Direct numerical simulation
In this section, the steady and unsteady DNS results for a hypersonic flow over a blunt cone subjected to an acoustic pulse are presented. For the current case, we adopted case I from He (Reference He2022) with the flow conditions presented in table 1.
Free-stream conditions of hypersonic flow over a blunt nose cone following He & Zhong (Reference He and Zhong2022).
With the low freestream total enthalpy,
$h_{0,\infty}$
given in table 1 along with other tunnel conditions, figure 2 shows the mean-flow contours of Mach number, pressure and temperature plotted on the physical coordinates
$X$
and
$Y$
. From figure 2(a), the boundary-layer development can be observed near the wall. Furthermore, the shock layer was converged with the shock as the upper boundary through the Rankine–Hugoniot relations. The mean-flow pressure near the nose region is shown in figure 2(b), normalised by the free-stream value. The pressure contour indicates that the pressure throughout the shock layer is nearly constant downstream at the centreline distance
$X =0.01 \,\rm m$
, which is relatively upstream compared with the full domain length of
$1.5 \,\rm m$
, which is at
$295.28$
nose radius lengths. Moreover, the temperature contour in figure 2(c) depicts an entropy layer being present and extending downstream. This entropy layer is associated with the transition reversal process observed in nose bluntness studies by Lei & Zhong (Reference Lei and Zhong2012) and Balakumar & Chou (Reference Balakumar and Chou2018). Before the reversal, the increasing nose bluntness reduces the local Reynolds number and delays the transition. However, when the nose bluntness exceeds a critical value, this reduction is reversed with no significant contribution from the second-mode instability waves, suggesting other forms of excitation. Traditionally, the existence of entropy-layer instabilities observed by Stetson (Reference Stetson1984) at low frequencies was identified with small amplification rates. Yet, another type of entropy-layer disturbance was found by Dietz & Hein (Reference Dietz and Hein1999) located at a second generalised inflection point outside the boundary-layer edge. Wan, Su & Chen (Reference Wan, Su and Chen2020) further concluded with the comparison between DNS and LST that this entropy-layer disturbance is dominant after the first mode decay and further induces instability when the entropy layer is swallowed by the developing boundary layer. Subsequently, Wan et al. (Reference Wan, Chen, Tu, Xiang, Yuan and Duan2023) identified this entropy-layer instability as an entropy-layer mode. The bi-orthogonal decomposition at the second-mode branch I neutral stability location in a later section also indicates that a higher-order discrete mode other than mode S contributes to the entropy-layer disturbance.
(a) Mach number, (b) normalised pressure and (c) normalised temperature contour near the blunt nose region. Note that the full domain of the cone ends at the central line distance of
$ X = 1.5$
m or
$295.28$
radius lengths.
Subsequently, an unsteady simulation is performed to obtain the perturbed flow field for bi-orthogonal decomposition. This unsteady simulation introduces a free-stream pulse in front of the blunt nose, and the flow field is recorded as the pulse propagates downstream. The finite spherical pulse parameters and schematic set-up are provided in the work of He & Zhong (Reference He and Zhong2022), referred to as case I. The slow acoustic disturbance in the free stream has a peak density perturbation amplitude,
$|\rho ^{\prime }|_{\infty }$
, of
$1 \times 10^{-6}$
as well as a peak pressure perturbation amplitude
$|P^{\prime }|_{\infty }=|\rho ^{\prime }|_{\infty }\gamma$
. The parameter
$\sigma$
for the finite spherical pulse is set to
$1\times 10^{-3}$
. The pulse is set to start the advection at a location
$x_0 = -0.02 \,\rm m$
with a slow acoustic disturbance speed of
$c_{\infty }=u_{\infty }-a_{\infty }$
, where
$u_{\infty}$
is the freestream velocity and
$a_{\infty}$
is the freestream speed of sound. Figure 3 presents the unsteady flow field in the time domain in a downstream region where the second-mode instability is dominant. The end of the finite pulse is shown in figure 3(a) at the left end of the domain. Following the finite pulse, the acoustic modulation inside the boundary layer indicates the existence of an excited boundary-layer mode, resembling a second-mode disturbance. Furthermore, the weaker disturbances radiating from the second-mode instability into the rest of the shock layer indicate the existence of a potential supersonic mode. This supersonic mode, as Knisely & Zhong (Reference Knisely and Zhong2019b
) have suggested, can be excited by the interactions between the unstable mode S, the stable mode F and the slow acoustic spectrum. The phase speed result shown in the next section indicates that such synchronisation might occur at the point where the normal modes intersect with the phase speed of
$c_r =1-(1/M) = 0.898$
. Thus, the influence of this possible supersonic mode can manifest itself as an oscillatory behaviour in the free stream, similar to the continuous modes. Thus, the bi-orthogonal decomposition was performed with acoustic continuous modes with results shown in a later section to further identify the source of free-stream behaviour.
Normalised second-mode pressure disturbance in the downstream region following the propagating pulse’s end towards the domain’s end. (a) The front of second-mode instability. (b) The end of second-mode instability.
A time history of the normalised total surface pressure perturbations at various streamwise positions, before and at the second-mode dominant region along the cone, is presented in figure 4. In the rest of the study, the streamwise coordinate
$s$
, defined as the arclength of the cone surface, is used when the streamwise locations are referred to for consistency with previous studies by He & Zhong (Reference He and Zhong2021a
, Reference He and Zhong2023). The magnitude of the wavepacket sampled at
$s = 0.2547\;{\rm{m}}\;(s/R_n=50.1)$
is relatively small compared with the magnitude of the wavepacket at the second-mode dominant region as expected. This upstream wave is also uncertain in wave shape and might contain a multitude of forcing waves excited by the receptivity mechanism. At
$s=0.5829\;{\rm{m}}\;(s/R_n=114.8)$
, which is right after the branch I location of the unstable mode S, the second-mode instability wavepacket starts to form. At this location, the uncertain forcing originated from the receptivity is still present and had a greater magnitude than the second-mode instability wavepacket. This forcing eventually decays as the packet propagates downstream to the unstable second-mode region, where the more defined modal wavepacket of the second-mode instability dominates. In addition, in the unstable region, a secondary wave packet is also present at
$s=1.0104\;{\rm{m}}\;(s/R_n=200)$
for the slow acoustic disturbance case. This structure includes a clearly amplified second mode as well as a nonlinear interaction represented by the valley between the two peaks. Previous work of Sivasubramanian & Fasel (Reference Sivasubramanian and Fasel2014) has suggested that this nonlinear behaviour is a result of the generation of higher and lower harmonics.
Time history of normalised surface pressure perturbations at different locations along the streamwise direction.
With the temporal disturbance data of the unsteady flow field depicted in figure 4, the time series data at various streamwise locations are then transformed by using FFT to obtain the spectral components for the bi-orthogonal decomposition. Utilising the FFT perturbation data at a particular frequency, a projection onto the discrete and continuous modes can be performed, enabling a comprehensive analysis of the flow field’s modal composition.
4.2. Discrete mode F and S decomposition for receptivity
In this section, the spectral components of the unsteady DNS flow field are decomposed by the discrete modes F and S. According to earlier LST results from Zou & Zhong (Reference Zou and Zhong2024), the unstable mode for this case is the discrete mode S, emerging from the slow acoustic line. The current study will decompose the unsteady flow field with discrete modes F and S to further validate this behaviour.
Previously in the works of Huang & Zhong (Reference Huang and Zhong2014) and He & Zhong (Reference He and Zhong2022), the disturbance behaviour was only studied through the combination of LST and the unsteady DNS result, without decomposing the signal into boundary-layer modes. With the bi-orthogonal decomposition, the unstable mode behaviour can be isolated in a mix of uncertain forcing and modal wavepacket for further analysis. Prior to presenting the discrete mode F and S projections on the unsteady DNS flow field, the discrete mode F and S eigenfunctions should be presented. Figure 5 shows the eigenfunctions for discrete modes F and S at the
$200 \,\rm kHz$
branch I location of
$s=0.5772\;{\rm{m}}\;(s/R=113.6)$
. Here,
$y/L$
denotes the wall-normal coordinate normalised by the length from the wall to the shock. From the eigenfunctions, the two discrete modes have perturbations that decay around the wall-normal position of
$y/L = 0.1$
, which is the location of the boundary-layer edge at this location. This decaying behaviour aligns with the general characteristic of a discrete boundary-layer mode. Furthermore, the discrete mode S is manifested with one peak inside the boundary layer as trapped acoustic waves, while mode F has multiple peaks near the end of the boundary layer. This spatial structure is not fixed and can vary with the streamwise location. Particularly for the results shown in figure 5, the phase speeds of mode F and the entropy-layer mode, which we will discuss later, have intersected each other and energy exchange may alter the mode shapes. In addition, the shock boundary conditions from the work of Knisely & Zhong (Reference Knisely and Zhong2019a
) were also implemented and no significant effect was observed in the resulting eigenfunctions.
Discrete mode F and S eigenfunctions with the decaying and shock boundary conditions at
$s=0.5772\;\textrm{m}\;(s/R=113.6)$
. Real part of (a) velocity,
$u_r$
. (b) Pressure,
$p_r$
. (c) Temperature,
$\theta_r$
.
Subsequently, with the discrete modes obtained, the bi-orthogonal decomposition can be performed by projecting the discrete-mode eigenfunctions to the spectral components of the DNS flow field with (3.8). Figure 6 presents various unsteady DNS flow fields and the corresponding discrete-mode projections along the streamwise direction. Overall, the discrete-mode projections indicate that both discrete modes F and S contribute to the disturbance inside the boundary layer upstream, while mode S becomes more dominant travelling downstream. The disturbance of figure 6(a–c) also indicates the existence of a high-amplitude entropy-layer perturbation, from about
$y/L=0.09$
to
$0.2$
, near the branch I region, visualised by a dominant peak residing in the entropy layer and extending into the free stream beyond the boundary layer. According to Tumin (Reference Tumin2007), the BES forms a complete set of solutions in the boundary layer. However, the current DNS result beyond the boundary layer also contains peaks attributed to the free-stream forcing from the shock layer, depicted in figure 4, and a possible discrete entropy-layer mode observed by Wan et al. (Reference Wan, Su and Chen2020), which will be discussed in a later section. Figure 6(d–i) further suggest that this perturbation in the entropy layer emerges closer to the wall and decays downstream as a result of the entropy layer being swallowed by the boundary layer along the streamwise direction. This observation aligns with previous findings in the work of Wan et al. (Reference Wan, Su and Chen2020) that the mode S instability could also be induced by an entropy-layer instability mode, referred to as the discrete entropy-layer mode in a later section. Investigating further into the source of such a disturbance, we found a discrete entropy-layer mode near the branch I location in addition to the discrete modes F and S. This additional discrete entropy-layer mode will be discussed in the next section. Eventually, in the second-mode dominant region shown in figure 6( j–l), the decomposed result indicates that the DNS flow field is composed solely of the discrete mode S, with the discrete mode F and entropy-layer disturbance amplitudes being negligible. Note that the reconstructed signal has a slightly higher amplitude than the FFT result due to parasitic errors of the numerical integration.
Discrete mode F and S projections,
$C_{F}A_{F}$
and
$C_{S}A_{S}$
, plotted against the unsteady DNS result at various streamwise locations. (a) Velocity at
$s/R=113.6$
. (b) Pressure at
$s/R=113.6$
. (c) Temperature at
$s/R=113.6$
. (d) Velocity at
$s/R=140.6$
. (e) Pressure at
$s/R=140.6$
. ( f) Temperature at
$s/R=140.6$
. (g) Velocity at
$s/R=157.2$
. (h) Pressure at
$s/R=157.2$
. (i) Temperature at
$s/R=157.2$
. ( j) Velocity at
$s/R =175.5$
. (k) Pressure at
$s/R =175.5$
. (l) Temperature at
$s/R =175.5$
.
Mode F and S amplitudes at various streamwise locations compared with
$\textrm{e}^N$
and DNS results. (a) Modal amplitudes. (b) Mode S amplitudes at branch I.
To further validate the accuracy of the decomposition, the modal amplitudes of the discrete-mode projections onto the unsteady DNS are also examined. Using the bi-orthogonal eigenfunctions of the discrete modes, the orthogonality condition in (3.8) is applied to obtain the modal coefficient for that particular mode. Subsequently, the modal coefficient is multiplied by the DNS surface amplitude,
$|A|_{\textit{DNS}}$
, to recover the modal amplitude. The surface pressure amplitude of
$C_s A_s$
is plotted against both the unsteady DNS and LST results in figure 7. Here
$C_s$
is the modal coefficient from the bi-orthogonal decomposition, used solely to only account for the mode S amplitude and similarly for the other modes. In the figure, two
$\textrm{e}^N$
approaches are included: the conventional
$\textrm{e}^N$
downstream predictions,
$|A|_{e^N}$
, using the decomposed initial amplitude, and the backtracked amplitudes,
$|A|_{e^N}$
(backtracked), using the method of He & Zhong (Reference He and Zhong2023). The initial amplitude from backtracking can be obtained as
$|A_0|=|A(s^*,f_n )/e^{N(s^*,f_n )} |$
where
$N=\int ^{s^*}_{s_0} - \alpha _i {\rm d}s^{\prime }$
,
$s^*$
is a downstream second-mode dominant location, and
$s_0$
refers to the branch I location. The branch I neutral stability point is located at
$s_0 = 0.5772\;{\rm{m}}\;(s/R=113.6)$
for the
$200 \,\rm kHz$
spectral component. The red dashed line represents the downstream amplitudes computed by the direct
$\textrm{e}^N$
method with the initial amplitude extracted by the decomposition. In comparison, the blue dashed line presents the backtracked amplitudes from a second-mode dominant location of
$s^* = 1.0104\;{\rm{m}}\;(s/R=200)$
. Nevertheless, the black dots represent the surface pressure amplitude obtained directly from the unsteady DNS. From the comparison, the pure DNS result displays a much higher surface pressure amplitude than both the
$\textrm{e}^N$
and the decomposed amplitudes near the branch I location. This higher amplitude has two possible causes: the presence of multiple prevalent modes other than mode S and the non-modal excitation as a result of multimode interaction. From the decomposition presented later, the superposition of the current discrete modes aligns closely with the DNS profile but with possible non-modal discrepancies. Moreover, the entropy-layer mode results in the following section also suggest a possible destabilising effect as the entropy layer is swallowed. Subsequently, the modal amplitudes of three different discrete modes, including mode F, mode S and the entropy-layer mode, are shown with solid lines. The decomposed mode S amplitudes align relatively well with the two
$\textrm{e}^N$
predictions in the second-mode dominant region, up to the branch I location. The alignment of the results suggests that only the initial amplitude is required to predict downstream modal amplitudes. Moreover, the modal amplitude becomes more oscillatory during the synchronisation between discrete modes, which occurs before the unstable mode branch I location. This result is expected because the bi-orthogonality relation of (3.8) depends on the quasi-parallel assumption and does not account for the non-parallel effects that cause the synchronisation. While synchronisation between these modes can enhance energy exchange and influence the amplification of instabilities, it is not the cause of the second-mode instability itself according to previous literature by Fedorov & Tumin (Reference Fedorov and Tumin2011). Nevertheless, the current bi-orthogonal decomposition method is capable of recovering the modal amplitude at the branch I location for receptivity analysis. Furthermore, the discrete mode F amplitudes are also plotted to show the decreasing trend upstream, which broadly aligns with the pure DNS result without oscillations. Moreover, the discrete mode F has been identified as computationally challenging to resolve in the synchronisation region from approximately
$s/R=95$
to
$110$
, since modes F and S synchronise and the eigenpairs resemble each other when the LST system is solved. A possible solution to resolve this phenomenon can be the incorporation of weakly non-parallel terms in the formulation suggested by Tumin (Reference Tumin2011) as well as using the parabolised stability equations (PSE). Since the focus of the current research is on the direct evaluation of receptivity data, the LST assumption is sufficient to extract the unstable mode with the BES. In addition, a newly discovered entropy-layer discrete mode has been identified in the synchronisation region. The presence of this entropy-layer mode explains the higher DNS amplitudes in this region. This mode contributes to the total disturbance and eventually also plays a role in the second-mode instability, as discussed in a later section.
Decomposed mode F and S (a) phase speed,
$c_r$
and (b) normalised growth rates,
$- \alpha_i$
, compared with DNS and LST results.
Moreover, the phase speed and growth rate for the discrete mode S are computed by applying the modal amplitudes obtained from the decomposition. Figure 8 compares the DNS phase speed and growth rate of He & Zhong (Reference He and Zhong2022) with the results calculated from the bi-orthogonal decomposition. Here, the growth rate is normalised by the reference length as per previous studies by He & Zhong (Reference He and Zhong2022) and Zou & Zhong (Reference Zou and Zhong2024). Near the branch I location, the DNS result of He and Zhong (Reference He and Zhong2022) exhibits high oscillations due to multimode effects. Applying the bi-orthogonal decomposition, the mode S amplitudes can be extracted, and the oscillations of the resultant growth rate near the branch I location are reduced. Note that near the branch I location, the numerical integration in the decomposition is performed up to the boundary-layer edge for the discrete modes to avoid additional numerical error. In addition, the phase speed and growth rates for discrete mode F are also recovered with reduced oscillations prior to the synchronisation region, where mode F is dominant. In the mode S dominant region, although with larger oscillations than the dominant mode, the phase speed and the growth rate of the non-dominant mode F can also be calculated to reveal the overall trend. It is worth noting that the non-dominant mode shows larger oscillations than the dominant mode. This occurs because, in the region where mode S dominates, mode F’s amplitude is very small. With small amplitudes, the influence of numerical error in the decomposition is enlarged. Nevertheless, the resultant trend is broadly consistent with the LST result and is not demonstrated from the DNS result alone. Furthermore, oscillations still persist in the decomposed results towards the synchronisation region, due primarily to the mathematical singularity that exists at the location where the phase speeds of modes F and S match. However, at branch I, which is downstream of the synchronisation region, the eigenfunctions and their corresponding adjoint modes remain distinct, allowing the bi-orthogonal decomposition to resolve the modal contributions.
Modes F and S surface pressure amplitude normalised by the freestream pressure at various streamwise locations compared with unsteady DNS results (log scale).
While the previous results focused on the 200 kHz decomposition, a broader frequency spectrum must be analysed to capture the full disturbance behaviour. Although He & Zhong (Reference He and Zhong2021a
) examined the broadband behaviour of the unsteady flow field, the theoretical analysis of bi-orthogonal decomposition has not been incorporated. Figure 9 illustrates the combined surface pressure amplitude,
${\rm d}p$
, of the discrete modes F and S across a
$160$
to
$260 \,\rm kHz$
band at various streamwise locations, normalised by the free-stream pressure in log scale. In general, the surface pressure amplitude from the unsteady DNS can be recovered with the combination of discrete modes S and F. Concurring with the LST
$N$
-factor result shown in previous work by He (Reference He2022), the pressure amplitude of the unstable-mode frequency band increases downstream while the unstable frequency band shifts to lower frequencies. Nevertheless, it is important to acknowledge that free-stream forcing and entropy-layer disturbances are also present in the flow. However, these disturbances reside mainly beyond the boundary layer and have a minor impact on the surface amplitude.
To highlight modal dominance at different frequencies, the individual modal surface pressure amplitudes are plotted in figures 10 and 11. Each figure includes a linear scale plot that focuses on the second-mode amplification and a zoomed plot that reveals additional details before the unstable region for lower frequencies. From figure 10, the discrete mode F is dominant before
$f=180\,\rm kHz$
for the cases from
$s=0.636\;{\rm{m}}\;(s/R=125)$
to
$0.799\;{\rm{m}}\;(s/R=157)$
. For locations
$s=0.844\;{\rm{m}} \;(s/R=166)$
and
$s=0.891\;{\rm{m}}\;(s/R=175)$
, the second-mode regions begin before
$f = 175 \,\rm kHz$
. For all locations, the transition of modal dominance from mode F to mode S is suggested before the unstable region. Notably, near the end of the unstable second-mode region, a mode F contribution remains visible, indicating intermodal energy exchange consistent with the DNS result in the streamwise direction from figure 8. Furthermore, figure 11 presents the decomposed mode S amplitudes, and a clear dominance can be observed in the unstable frequency band. The intermodal exchange between mode F and S is manifested from
$f=175\,\rm$
to
$200 \,\rm kHz$
for the selected streamwise locations. From both mode F and S results shown, the physical insight from the decomposed result agrees with previous unsteady DNS findings by He & Zhong (Reference He and Zhong2022) that oscillations exist before and after the unstable region. The current result also confirms that a discrete mode F contributes to a part of the DNS oscillations.
Mode F normalised surface pressure amplitudes at various streamwise locations compared with unsteady DNS results. (a) Mode F. (b) Mode F (zoomed).
Mode S normalised surface pressure amplitudes at various streamwise locations compared with unsteady DNS results. (a) Mode S. (b) Mode S (zoomed).
Moreover, the spectral receptivity coefficient for each of the disturbances is calculated from (3.10) using decomposed surface pressure perturbation data. Following (3.10), the decomposed spectral disturbance data are normalised by the amplitude spectra of the incident pulse. In previous work by He & Zhong (Reference He and Zhong2022), the receptivity coefficient obtained via the backtracking method was shown to be an indirect measure that exhibits oscillations, suggesting free-stream forcing and multimode contributions. To address this issue, we now calculate the receptivity coefficient directly from the decomposed mode S amplitude at the branch I neutral stable location. In figure 12, the decomposed receptivity coefficients are broadly in agreement with the DNS backtracking result but with smaller oscillation amplitudes. Since the backtracking method scales the DNS amplitude directly using the
$\textrm{e}^N$
method from the second-mode region, the large oscillations in the backtracking result stem from the DNS flow field that contains the free-stream forcing penetrating the boundary layer. The bi-orthogonal decomposition is shown to reduce such an influence. Although the decomposed amplitudes at branch I exceed those of the backtracking method in the most amplified frequency band, the overall trend of the receptivity coefficients aligns with the previous result from He (Reference He2022), increasing toward the peak near
$235\,\rm kHz$
. He & Zhong (Reference He and Zhong2023) further suggested that intermodal interactions in the backtracking approach lower the receptivity coefficient, contributing to their prediction of subsequent downstream transition location compared with experiments. In contrast, the bi-orthogonal decomposition appears to exclude these intermodal effects, yielding higher amplitudes. These increased receptivity coefficients, in theory, would yield an earlier transition location. Hence, an important future step will be the application of such coefficients into the amplitude method to validate the predicted transition location. The expansion of a broader frequency range would also be considered to completely capture the receptivity spectrum.
Receptivity coefficients,
$c_{rec}$
, of discrete mode S plotted against results from He & Zhong (Reference He and Zhong2023).
4.3. Flow-field recovery with higher discrete and continuous modes
4.3.1. Discrete entropy-layer mode
In addition to the analysis of discrete modes F and S, the flow physics near the branch I location is less understood due to the lack of a multimode decomposition. In these synchronisation regions, various discrete and continuous normal modes interact and participate in intermodal energy exchanges. Specifically, existing literature including Fedorov (Reference Fedorov2003) and He & Zhong (Reference He and Zhong2022) attributed the synchronisation between the discrete modes F and S to giving rise to the second-mode instability. However, a recent study by Wan et al. (Reference Wan, Su and Chen2020) suggests the possibility of second-mode excitation due to the entropy-layer instabilities. Thus, the bi-orthogonal decomposition of the unsteady flow field with three discrete modes – F, S and entropy-layer mode – is performed with the goal to better understand the modal contributions.
In addition to the discrete modes F and S, an entropy-layer mode has been identified near the synchronisation region. The phase speed and growth rate of this discrete entropy-layer mode are plotted over streamwise locations in figure 13. A grid convergence study was also performed with
$N=300$
and
$N=400$
points. The two results aligned well with small discrepancies on the upstream growth rate since the growth rate is normalised by the reference length that amplified the difference. From the phase speed plot, a discrete mode independent of modes F and S exists near the synchronisation region. This mode emerges upstream between the two conventional discrete modes and approaches the phase speed of an entropy/vorticity mode (
$c_r=1$
). Similarly, the growth rate shows a trend approaching a neutral wave similar to the vorticity and entropy continuous modes. Although such an entropy-layer mode has been shown to exist before synchronisation of discrete modes F and S in previous studies by Chen et al. (Reference Chen, Tu, Wan, Su, Yuan and Chen2021) and Wan et al. (Reference Wan, Chen, Tu, Xiang, Yuan and Duan2023), the current numerical method is capable of resolving the entropy-layer mode downstream into the continuous entropy/vorticity spectra, indicating the possibility of coexistence. This behaviour also aligns with previous LST results in the work of He (Reference He and Zhong2022) suggesting that reduced bluntness would lead to an earlier second-mode amplification due to a thinner entropy layer. From a modal perspective, as the entropy layer merges into the boundary layer, the discrete mode F experiences a reduction in phase speed. This shift facilitates phase speed alignment with entropy-layer modes, enabling energy exchange between the entropy- and boundary-layer modes. Such synchronisation can enhance intermodal coupling and may further promote the destabilisation of the second mode.
(a) Phase speed and (b) growth rate of the entropy-layer mode (E) merges into the continuous spectra, plotted against discrete modes F and S.
Real and imaginary parts of the entropy-layer mode eigenfunction at
$200\,\rm kHz$
at the branch I location (
$s/R =113.6$
). (a) Velocity. (b) Pressure. (c) Temperature.
Subsequently, figure 14 presents the eigenfunctions of the discrete entropy-layer mode at
$200\,\rm kHz$
at the branch I location of the unstable second mode. From the eigenfunction, the discrete entropy-layer mode demonstrates an extended perturbation into the entropy layer around
$y/L= 0.15$
with magnitudes larger than those inside the boundary layer. Furthermore, this discrete mode also aligns with the behaviour of a discrete mode, in which the perturbation decays in the free stream after the peaks in the entropy layer. Moreover, the magnitude of the discrete entropy-layer mode is relatively large compared with the unsteady DNS perturbation at the same location shown in figure 6(a). Thus, a small contribution of this entropy-layer mode would have a large effect on the unsteady flow field. Furthermore, the large amplitude of the entropy-layer mode might trigger non-modal responses in the flow field as the entropy layer is swallowed (Wan et al. Reference Wan, Chen, Tu, Xiang, Yuan and Duan2023). This non-modal effect may manifest as discrepancies between the current DNS result and the reconstructed modal signal and can be further captured by the input–output analysis conducted by Cook & Nichols (Reference Cook and Nichols2022, Reference Cook and Nichols2024).
Streamwise velocity component of the reconstructed signal over various streamwise locations plotted against the
$200\,\rm kHz$
unsteady FFT flow field. (a)
$s/R=113.6$
, (b)
$s/R=125.2$
, (c)
$s/R=140.6$
, (d)
$s/R=113.6$
, (e)
$s/R=125.2$
, ( f)
$s/R=140.6$
.
The inclusion of this entropy-layer mode in the bi-orthogonal decomposition over various streamwise locations from the branch I location of the unstable mode is also performed. With the discrete modes F, S and the entropy-layer mode, figure 15 presents the decomposed result of the unsteady flow field near the synchronisation region at
$200\,\rm kHz$
. Both the modal amplitudes and the combined signal to reconstruct the DNS flow field are shown. From the modal amplitudes, both mode F and mode S can be observed inside the boundary layer while the discrete entropy-layer mode contribution,
$C_{E}A_{E}$
, is responsible for the second peak of the unsteady DNS flow field. The observation aligns with previous finding that both the discrete mode S and the entropy-layer mode are present near the branch I location of the unstable mode. Similar to Wan et al. (Reference Wan, Chen, Tu, Xiang, Yuan and Duan2023), the LST assumption in the formulation would cause a discrepancy between the modal amplitude and the DNS flow field. The PSE can be implemented to better account for the non-parallel effect, such that the amplitude peaks are more aligned Wan et al. (Reference Wan, Chen, Tu, Xiang, Yuan and Duan2023). In addition, non-modal effects discussed by Cook et al. (Reference Cook, Thome, Brock, Nichols and Candler2018), Nichols & Candler (Reference Nichols and Candler2019) and Kumar et al. (Reference Kumar, Sarath, Sasidharan and Gaitonde2025) could also contribute to oscillations in the unsteady DNS. However, direct comparisons between the input–output modes and the wall-normal DNS profile have not been presented in the literature. The existence of the broadband free-stream forcing should also be noted for a possible source of the amplitude mismatch beyond the entropy layer into the shock layer. Further investigation downstream at
$s/R=140.6$
also indicates that the entropy-layer disturbance decays as the entropy layer is swallowed, and the unstable boundary-layer mode is dominant. This finding can also be interpreted as an energy transfer between the entropy-layer mode and the unstable mode S when the entropy layer is swallowed. Beyond the second peak of the DNS flow field, wave structures demonstrate the existence of continuous modes and non-modal free-stream forcing that require additional investigation.
To further assess the robustness and physical relevance of the entropy-layer mode, linear stability analysis was repeated for a blunter nose configuration with an increased nose radius of
$9.525\,\rm mm$
. The mean flow conditions are similar to the previous case with a Mach number
$M=9.79$
, and a stagnation enthalpy of
$1.07$
. For the detailed flow conditions, one can refer to the work of He & Zhong (Reference He and Zhong2021a
). The entropy-layer mode was successfully recovered in the blunter case, exhibiting similar structural features but with enhanced amplitude near the boundary-layer edge and deeper spatial penetration. Figure 16 plots the eigenfunctions of the discrete modes at the same dimensional streamwise location of
$s=0.5772\,\rm m$
. Interestingly, while the entropy layer is broader in the blunter case, the entropy-layer mode exhibits a more confined wall-normal structure, with its energy concentrated closer to the boundary-layer edge. Notably, the mode also displays a higher amplitude relative to the baseline case. This suggests stronger trapping and sharper thermodynamic gradients, which may enhance the local energy density and strengthen the mode’s interaction with boundary-layer modes.
Eigenfunction of the entropy-layer mode at
$200\,\rm kHz$
for a cone with
$9.525\,\rm mm$
nose radius bluntness. (a) Velocity. (b) Pressure. (c) Temperature.
(a) Phase speed and (b) growth rate of the entropy-layer mode merges into the continuous spectra, plotted against discrete modes F and S for a nose radius of
$9.525 \,\rm mm$
bluntness.
In addition, to validate the entropy-layer mode with prior literature, the phase speed and growth rate of the entropy-layer mode at the frequency of
$150\,\rm kHz$
are plotted. Figure 17 plots the phase speed and growth rate of the entropy-layer mode at
$150\,\rm kHz$
, showing the mode that approaches the continuous spectra compared with the discrete-mode LST results reported by He & Zhong (Reference He and Zhong2021a
). From the phase-speed plot, the entropy-layer mode meets the discrete mode F at a more upstream location than the
$200\,\rm kHz$
case. This is expected since the frequency is reduced. Nevertheless, these trends confirm the existence of the entropy-layer mode and suggest possible synchronisation with the boundary-layer modes. Future work applying bi-orthogonal decomposition to unsteady DNS data for the blunter nose configuration will help provide additional physical insight into the entropy-layer mode’s role and its contribution to receptivity.
5. Summary and conclusions
The primary goal of this study was established as directly obtaining receptivity coefficients at the branch I location to enhance the accuracy of the current state-of-the-art amplitude transition-prediction method for realistic systems in atmospheric flight conditions. The previous
$\textrm{e}^N$
backtracking method to compute receptivity coefficients requires conducting DNS beyond the branch I location and inherits the oscillating behaviour attributed to the free-stream forcing. To overcome the current limitation of the lack of accurate receptivity data, a bi-orthogonal decomposition of unsteady DNS receptivity flow field data obtained by He & Zhong (Reference He and Zhong2022) was successfully performed. Using this decomposition, the receptivity data of a hypersonic flow over a blunt cone were extracted from the branch I location. The multimodal analysis for a hypersonic flow over a cone poses unique challenges, particularly regarding flow field reconstruction, due to weakly non-parallel effects and entropy-layer swallowing phenomena. Previous studies of Miselis et al. (Reference Miselis, Huang and Zhong2016), Saikia et al. (Reference Saikia, Al Hasnine and Brehm2022) and Wan et al. (Reference Wan, Chen, Tu, Xiang, Yuan and Duan2023) have addressed only specific aspects of this broader problem. In contrast, current research involves using the high-order finite difference method to systematically compute discrete and continuous modes, including an entropy-layer discrete mode, and applying the orthogonality relation to Fourier-transformed flow field data at a specified frequency. By projecting the perturbation variables onto the bi-orthogonal eigenfunctions, the modal perturbation amplitudes are extracted, enabling a clear separation of the different flow modes.
Remarkably, from the preliminary results of the mode S decomposition along the streamwise direction, an overall trend in the discrete mode S amplitude is observed that is in good agreement with previous results starting from the branch I location. From the phase speed and growth rate comparisons, the bi-orthogonal decomposition yields results with reduced oscillations near the branch I neutral location relative to earlier DNS studies by He & Zhong (Reference He and Zhong2021a , Reference He and Zhong2022 and Reference He and Zhong2023). These dominant mode trends also show better consistency with the LST predictions. Moreover, the general non-dominant mode characteristics can also be isolated from the DNS. Subsequently, the study computes the discrete mode S receptivity coefficients for a band of frequencies. These coefficients not only align with trends from previous studies, but also exhibit a reduction in oscillation due to free-stream forcing. This is a promising result, indicating that the bi-orthogonal decomposition can enhance the robustness of receptivity coefficient calculations. Further extending the analysis to various types of free-stream disturbances within the amplitude method would be critical for gaining insight and better explaining experimental phenomena.
Furthermore, the decomposed results at a downstream location are aligned with the observation that the discrete mode S is dominant, while the decomposition of the flow field near the branch I neutral location indicates multimode contributions with perturbations in the entropy layer. In addition, near the synchronisation region, the existence of a discrete entropy-layer mode is demonstrated. This mode is obtained with the very high-order finite difference scheme proposed by Zou & Zhong (Reference Zou and Zhong2023). Although such an entropy-layer mode has been shown to exist before the emergence of discrete mode S in a previous study by Wan et al. (Reference Wan, Chen, Tu, Xiang, Yuan and Duan2023), the current result indicates that the disturbance of the entropy-layer mode coexists with the unstable mode. Near the branch I location, the entropy-layer perturbation attributed to this mode is larger in magnitude than the unstable mode S. As the entropy layer is swallowed, the entropy-layer mode decays as the unstable second mode develops. It is suspected that the entropy-layer mode, in addition to the synchronisation of modes F and S, may contribute to the destabilisation of the unstable second mode. Further investigation of the disturbance energy with methods such as the input–output analysis would be a next step in analysing the non-modal interaction between the entropy-layer mode and the unstable second mode. Moreover, the role of this entropy-layer mode can be investigated further for various flow scenarios such as supercritical flow and non-ideal fluids where the unstable mode differs from Mack’s second mode, suggested by Ren, Marxen & Pecnik (Reference Ren, Marxen and Pecnik2019). In addition, since the branch I neutral stable location is relatively downstream from the start location of the discrete modes, the process of how fast and slow acoustic waves transform into discrete modes are not focused on in this study. Thus, the decomposition using continuous spectra is not presented. Such continuous modes may have a more significant contribution to the flow field upstream, near the leading-edge region.
This paper presents a combined methodology in theoretical and computational analysis in an attempt to gain deeper physical insight into the development of second-mode instability in hypersonic boundary-layer flows over a blunt body through direct modal analysis. Specifically, it is the first time that the bi-orthogonal eigenfunction decomposition method, incorporating both boundary-layer and entropy-layer modes, has been applied directly at the branch I location for hypersonic boundary-layer flows over a blunt cone. This approach limits the effect of the inherited forcing oscillations on the receptivity coefficients, and reduces the need for conducting DNS beyond the branch I location. Overall, current validation confirms the decomposition’s reliability, aligning closely with previous theoretical and numerical predictions. As the amplitude method described by Mack (Reference Mack1977) critically depends on accurately determining the initial disturbance amplitude, the receptivity coefficient derived from the modal decomposition of the initial perturbation is anticipated to significantly enhance the accuracy of laminar-to-turbulent transition predictions. Consequently, the immediate future application of this research will involve integrating these receptivity results into Mack’s amplitude method to more precisely predict boundary-layer transition.
Funding
The first and third authors would like to acknowledge the support by the Office of Naval Research (ONR) grant no. N000142412371, monitored by Dr Marineau. The second author was supported partially by the Air Force Office of Scientific Research (AFOSR) under AFOSR grant no. FA9550-19-1-0206 and by the ONR under grant no. N00014-17-1-2343 during his time of contribution to this research at UCLA. Primary computational resources used for the hypersonic unsteady flow are provided by ACCESS through the San Diego Super-Computer Center (SDSC). The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the AFOSR, ONR, XSEDE or the US Government.
Declaration of interest
The authors report no conflict of interest.
Appendix A. The matrix elements
The matrix elements of
$\boldsymbol{H}$
in (2.12) are presented here. Similar to Tumin (Reference Tumin2007) and Miselis et al. (Reference Miselis, Huang and Zhong2016), we define
$r=2(\epsilon +2) / 3$
and
$m=2(\epsilon -1) / 3$
where
$\epsilon$
is the ratio of bulk viscosity to dynamic viscosity and equals to
$0$
following Stoke’s hypothesis. For entries not specified, the coefficients are equal to zero. We first define
\begin{align} \chi &=\left [\frac {\textit{Re}}{\mu }-\frac {i \hat {\omega } r}{P}\right ]^{-1} \nonumber\\ \hat {\omega }&=\omega -\alpha U. \end{align}
For
$\boldsymbol{H}$
, we have
\begin{align} & H_0^{1,2}=1, \end{align}
\begin{align} & H_0^{2,1}=\alpha ^2+\beta ^2-i \hat {\omega } P \gamma M^2 \frac {\textit{Re}}{\mu T}, \end{align}
\begin{align} & H_0^{2,2}=-\frac {1}{\mu } \frac {\partial \mu }{\partial y}, \end{align}
\begin{align} & H_0^{2,3}=-i \alpha (m+1) \frac {1}{T} \frac {\partial T}{\partial y}-i \alpha \frac {1}{\mu } \frac {\partial \mu }{\partial y}+{\textit{Re}} P \gamma M^2 \frac {1}{\mu T} \frac {\partial U}{\partial y}+i \alpha (m+1) \frac {1}{P} \frac {\partial P}{\partial y}, \end{align}
\begin{align} & H_0^{2,4}=i \alpha \frac {\textit{Re}}{\mu }+\alpha \hat {\omega }(m+1) \frac {1}{P}, \end{align}
\begin{align} & H_0^{2,5}=-\alpha \hat {\omega }(m+1) \frac {1}{T}-\frac {1}{\mu }\left (\frac {\partial T}{\partial y} \frac {\partial U}{\partial y} \frac {\partial ^2 \mu }{\partial T^2}+\frac {\partial \mu }{\partial T} \frac {\partial ^2 U}{\partial y^2}\right)\!, \end{align}
\begin{align} & H_0^{2,6}=-\frac {1}{\mu } \frac {\partial \mu }{\partial T} \frac {\partial U}{\partial y}, \end{align}
\begin{align} & H_0^{3,1}=-i \alpha, \end{align}
\begin{align} & H_0^{3,3}=\frac {1}{T} \frac {\partial T}{\partial y}-\frac {1}{P} \frac {\partial P}{\partial y}, \end{align}
\begin{align} & H_0^{3,4}=\frac {i \hat {\omega }}{P}, \end{align}
\begin{align} & H_0^{3,5}=-\frac {i \hat {\omega }}{T}, \end{align}
\begin{align} & H_0^{3,7}=-i \beta, \end{align}
\begin{align} & H_0^{4,1}=-i \chi \alpha \left (r \frac {1}{T} \frac {\partial T}{\partial y}+\frac {2}{\mu } \frac {\partial \mu }{\partial y}-r \frac {1}{P} \frac {\partial P}{\partial y}\right )\!, \end{align}
\begin{align} & H_0^{4,2}=-i \chi \alpha, \end{align}
\begin{align} & H_0^{4,3}=\chi \left (-\alpha ^2-\beta ^2+i \hat {\omega } P \gamma M^2 \frac {R e}{\mu T}+r \frac {1}{T} \frac {\partial ^2 T}{\partial y^2}+r \frac {1}{\mu T} \frac {\partial \mu }{\partial y} \frac {\partial T}{\partial y}-r \frac {1}{\mu P} \frac {\partial \mu }{\partial y} \frac {\partial P}{\partial y}\right .\nonumber \\ & \left .\qquad\,\, +\,2 r \frac {1}{P^2}\left [\frac {\partial P}{\partial y}\right ]^2-2 r \frac {1}{P T} \frac {\partial P}{\partial y} \frac {\partial T}{\partial y} -r \frac {1}{P} \frac {\partial ^2 P}{\partial y^2}\right )\!, \end{align}
\begin{align} & H_0^{4,4}=-i \chi r \frac {1}{P}\left (\alpha \frac {\partial U}{\partial y}-\hat {\omega } \frac {1}{T} \frac {\partial T}{\partial y}-\hat {\omega } \frac {1}{\mu } \frac {\partial \mu }{\partial y}+2 \hat {\omega } \frac {1}{P T} \frac {\partial P}{\partial y}\right )\!, \end{align}
\begin{align} & H_0^{4,5}=i \chi \left (\alpha r \frac {1}{T} \frac {\partial U}{\partial y}+\alpha \frac {1}{\mu } \frac {\partial \mu }{\partial T} \frac {\partial U}{\partial y}-r \hat {\omega } \frac {1}{\mu T} \frac {\partial \mu }{\partial y}+r \hat {\omega } \frac {1}{P T} \frac {\partial P}{\partial y}\right )\!, \end{align}
\begin{align} & H_0^{4,6}=-i \chi r \hat {\omega } \frac {1}{T}, \end{align}
\begin{align} & H_0^{4,7}=-i \chi \beta \left (r \frac {1}{T} \frac {\partial T}{\partial y}+2 \frac {1}{\mu } \frac {\partial \mu }{\partial y}-r \frac {1}{P} \frac {\partial P}{\partial y}\right )\!, \end{align}
\begin{align} & H_0^{4,8}=-i \chi \beta, \end{align}
\begin{align} & H_0^{5,6}=1, \end{align}
\begin{align} & H_0^{6,2}=-2(\gamma -1) M^2 {\textit{Pr}} \frac {\partial U}{\partial y}, \end{align}
\begin{align} & H_0^{6,3}=-2 i \alpha (\gamma -1) M^2 {\textit{Pr}} \frac {\partial U}{\partial y}+P\operatorname {\gamma } M^2 {\textit{Pr}} \frac {{\textit{Re}}}{\mu T} \frac {\partial T}{\partial y}-(\gamma -1) M^2 {\textit{Pr}} \frac {{\textit{Re}}}{\mu } \frac {\partial P}{\partial y}, \end{align}
\begin{align} & H_0^{6,4}=i \hat {\omega }(\gamma -1) M^2 {\textit{Pr}} \frac {{\textit{Re}}}{\mu }, \end{align}
\begin{align} & H_0^{6,5}=\alpha ^2+\beta ^2-i \hat {\omega } {\textit{Pr}} M^2 {\textit{Pr}} \frac {{\textit{Re}}}{\mu T}-(\gamma -1) M^2 {\textit{Pr}} \frac {1}{\mu } \frac {\partial \mu }{\partial T}\left (\frac {\partial U}{\partial y}\right )^2-\frac {1}{\mu } \frac {\partial ^2 \mu }{\partial y^2}, \end{align}
\begin{align} & H_0^{6,6}=-\frac {2}{\mu } \frac {\partial \mu }{\partial y}, \end{align}
\begin{align} & H_0^{7,8}=1, \end{align}
\begin{align} & H_0^{8,3}=-i \beta (m+1) \frac {1}{T} \frac {\partial T}{\partial y}-i \beta \frac {1}{\mu } \partial \mu \partial y+i \beta (m+1) \frac {1}{P} \frac {\partial P}{\partial y}, \end{align}
\begin{align} & H_0^{8,4}=\beta \hat {\omega }(m+1) \frac {1}{P}+i \beta \frac {\textit{Re}}{\mu }, \end{align}
\begin{align} & H_0^{8,5}=-\beta \hat {\omega }(m+1) \frac {1}{T}, \end{align}
\begin{align} & H_0^{8,7}=\alpha ^2+\beta ^2-i \hat {\omega } P \gamma M^2 \frac {\textit{Re}}{\mu T}, \end{align}
\begin{align} & H_0^{8,8}=-\frac {1}{\mu } \frac {\partial \mu }{\partial y}. \end{align}
