Disturbance energy budget of linear and nonlinear second-mode waves

Abstract

Linear and nonlinear mechanisms governing the growth of second-mode waves are analysed using a newly derived disturbance energy conservation equation that highlights the physical processes responsible for fluctuation energy production, flux-based transport and destruction. Axisymmetric direct numerical simulations (DNS) data from a Mach 6 hypersonic boundary layer, simulated over a $3^\circ$ half-angle sharp cone at zero angle of attack, is used as a reference. A Legendre polynomial-based forcing methodology is used to trigger transition in the DNS over a range of various amplitude levels and different forcing frequency content. Closure of the disturbance energy budgets is demonstrated numerically using the DNS data. The terms responsible for the amplification of the disturbance are identified, and nonlinear attenuation effects are discussed. We show that the interaction between the entropy/velocity fluctuations and the base temperature gradient governs the second-mode growth in the linear regime. Energy production occurs in the critical layer due to non-isentropic processes and accumulates in acoustic form below the relative sonic line through downward transport. At higher forcing amplitudes, nonlinear spectral broadening is observed, with simultaneous thermoviscous diffusion attenuating the disturbance energy growth. This effect is responsible for the non-monotonic streamwise variation of the wall-pressure spectrum. Phase speed and growth rate analyses, informed by linear stability theory (LST), reveal wave steepening effects preceding this nonlinear attenuation effect. The disturbance energy is observed to match the LST predictions at lower forcing amplitudes, deviating, as expected, at higher amplitudes.

1. Introduction

Hypersonic laminar-to-turbulent transition has been the subject of numerous studies for the past 50 years (Lees 1957). A turbulent boundary layer (BL) experiences higher skin friction drag and heat transfer than a laminar one. Hence, it is desirable to predict and control the laminar–turbulent transition process, especially for high-speed vehicles (Reshotko 1984); Kimmel (2003). Apart from the engineering benefit, the physics of transition has continually involved scientific gaps that have been the subject of many experimental (Lee & Chen 2019) and theoretical (Fedorov 2011) studies. The transition phenomenon begins with the natural introduction of instabilities into the BL, either due to free stream disturbances or generated by surface irregularities. These instabilities are then amplified through different mechanisms, leading to the eventual breakdown to turbulence. This process can occur through several pathways – a general overview of which is given in the seminal works of Morkovin, Reshotko & Herbert (1994). For wall-bounded flows with low-amplitude fluctuations, the typical path to transition is via the linear amplification of waves of specific frequencies. These unstable frequencies are sensitive to the physical characteristics of the BL as well as the initial disturbance levels in the flow, which creates uncertainty in the turbulence onset prediction (Schneider 2004). Although the linear regime has been extensively studied, little insight is available on the physics of the nonlinear breakdown phase.

1.1. Transition in high-speed flows

For flows in a low-disturbance environment and over slender geometries, the BL selectively amplifies frequencies (or modes) from a broadband disturbance spectrum. As the disturbance propagates downstream, those unstable modes are progressively amplified until nonlinear interaction between modes is triggered. Transition in low-speed flows (Ma $\lt$ 1) can be attributed to the modal amplification of Tollmien–Schlichting waves, which are longitudinal, vortical modes affected by viscous action at the no-slip wall. At higher Mach numbers (Ma $\gt$ 1), the viscous instability mechanism weakens and the inviscid instabilities are amplified. Lees & Lin (1946) suggest that this is due to the presence of a generalized inflection point $ (\partial /\partial y (\rho _0 \partial u_0/\partial y ) )$ in the BL, which forces transition to be dominated by first modes (Mack 1990), which are typically oblique in nature. Liang et al. (2023) note that the mechanisms sustaining the first Mack modes are concentrated at the inflection point off the wall, in contrast to the Tollmien–Schlichting waves driven by near-wall relative viscous phasing of streamwise and wall-normal fluctuations. Higher Mach numbers (Ma $\gt$ 4) and lower wall-to-adiabatic temperature ratios stabilize the first mode while destabilizing the second modes, which are fundamentally inviscid but not isentropic in nature, as shown in the present manuscript. In all hypersonic BLs, a region exists where the mean flow is supersonic relative to the propagating instability waves, defined by the relative sonic line. This region allows the disturbances to form a standing wave-like pattern and acts as a waveguide trapping a family of acoustic modes (Mack 1975), out of which the second mode is the fundamental harmonic. These two-dimensional waves, first predicted by Mack through a linear stability analysis (Mack 1969), were observed experimentally on slender cones by Demetriades (1974). The experiments of Stetson & Kimmel (1993) on a sharp cone in Mach $6$ flow confirmed the second mode as the primary instability responsible for hypersonic BL transition.

1.2. DNS/LST studies in hypersonic transition

The onset of transition is usually predicted by BL stability analysis. Locally parallel linear stability theory (LST) (Mack 1984) and nonlinear parabolized stability equations (PSE) (Bertolotti, Herbert & Spalart 1992) have received widespread attention for transition prediction aided by the $e^N$ method by Smith (1956). Alba et al. (2010) performed a detailed comparison of the experimental second-mode growths with PSE predictions for $0^\circ$ angle of attack hypersonic flows over a $7^\circ$ half-angle slender cone. Contemporarily, it was shown by Johnson, Candler & Alba (2010) that three-dimensional PSE results for a hypersonic BL are the same as two-dimensional results for canonical geometries, consistent with the fact that the unstable second modes are inherently two-dimensional in nature. The readiness with which the unstable mode frequencies can be identified has made these methods very attractive in the transition community (Reed, Saric & Arnal 1996).

With the advent of considerable computing power, direct numerical simulations (DNS) of high-speed flows have received renewed interest. Accurate prediction of the transition process involves accurate resolution of the spatial and temporal scales in the transitional flow (Zhong & Wang 2012). The advantage of DNS over the linearized stability methods lies in the fact that a complete numerical solution of the governing equations is attempted without any simplifying assumptions. Thus, nonlinear mechanisms can be captured by DNS, as done in the present study. Nonlinear breakdown to turbulence has been theorized to occur through three distinct mechanisms. Firstly, the fundamental resonance (Klebanoff or K-type) occurs when the two-dimensional primary longitudinal mode interacts with two oblique waves of the same frequency, having equal and opposite wave angles (Bake, Meyer & Rist 2002). Secondly, the subharmonic resonance (Herbert or H-type) involves the interaction of primary second modes with two subharmonic oblique modes having frequencies half of that of the primary mode (Kosinov, Maslov & Semionov 1997). Finally, the oblique breakdown occurs due to the interaction of two high amplitude three-dimensional oblique waves (Fasel, Thumm & Bestek 1993). Sivasubramanian & Fasel (2014) performed three-dimensional DNS on a Mach $6$ BL over a $7^\circ$ half-angle sharp cone. The axisymmetric second modes were reported to be the primary instability. A stronger fundamental resonance of the primary modes with the oblique modes was observed, with a weaker subharmonic resonance. Similar behaviour was observed for transition on a $2.5^\circ$ half-angle slender cone (Hader, Leinemann & Fasel 2020) and a flared cone (Hader & Fasel 2019) exposed to Mach $6$ flow.

The analysis of nonlinear mechanisms in transition is usually done from the point of view of mode-to-mode interactions, which is a logical successor to the modal stability theory. The DNS or PSE can accurately simulate these modal interactions, but little has been explored on the physical mechanisms regulating the linear and nonlinear disturbance energy budgets for second-mode waves.

1.3. Disturbance energy budget

A different route to analysing the stability of a fluid flow system is through disturbance energy norms. This involves constructing a suitable inner product in the vector space $\mathcal {V}$ , characterizing the disturbance state $ (\mathbf {v}=(p',u'_i,s')^T,\,\mathbf {v}\in \mathcal {V} )$ . The norm so defined offers a metric of the disturbance energy level. The mathematical description of the energy norm has been a topic of much debate, and the general consensus is that it should describe an instability’s acoustic, kinetic and thermal energies. Giauque, Nicoud & Brear (2007) listed a set of mathematical constraints to be respected by an energy norm: (i) the norm is required to be quadratic and positive definite, i.e., it should increase with fluctuation amplitudes; (ii) the energy norm characterizing a viscous or inviscid flow system should monotonically decrease or stay constant, respectively, in the absence of any external sources/sinks or production mechanisms.

The classic conservation equation of acoustic energy given by Rayleigh (1894) started the discussion on disturbance energy characterization in the context of combustion instability. Cantrell & Hart (1964) spearheaded one of the earliest efforts to include mean flow effects in the disturbance energy, assuming an irrotational and homentropic base flow, which was later extended by Morfey (1971) to non-isentropic flows. Chu (1965) critically analysed Rayleigh’s acoustic energy criterion and advocated the inclusion of thermal fluctuation terms describing entropy transport effects. These modifications were finally compiled by Myers (1991) to obtain a second-order equation for disturbance energy constructed from the evolution equations of first-order perturbations. This is termed as the first-order energy corollary. Myers later extended his work to an exact energy corollary, valid for any perturbation amplitude. Later, Giauque et al. (2006) augmented Myers’ norm by including thermochemical non-equilibrium effects. Widely accepted as the most complete description of disturbance energy evolution, Myers’ work serves as a constant guide for further studies on disturbance energy corollaries and inspires the current work. However, in the context of thermoacoustic instabilities, Chu’s energy norm has been deemed more suitable to obtain a stability criterion (Nicoud & Poinsot 2005), which was later supported by George & Sujith (2011). Myers’ norm demonstrated spurious growth in the absence of external sources, while Chu’s norm did not. The absence of spurious growth can be attributed to the self-adjointness of the linearized operator describing the system under investigation with respect to the energy norm given by Chu (1965). George & Sujith (2011) shows that including mean flow velocity terms in Myers’ energy norm introduces non-orthogonality between the acoustic and entropy modes, leading to transient growth. Gupta, Lodato & Scalo (2017) suggested modifying the first-order acoustic energy by including a correction term accounting for second-order isentropic fluctuations. No transient growth was found in high-fidelity simulations of travelling nonlinear thermoacoustic acoustic waves (Gupta & Scalo 2018) using the asserted approach.

Relevant to hypersonic BL transition, Kuehl (2018) used the first-order isentropic acoustic energy corollary to study the thermoacoustic interaction in second-mode instabilities. Thermoacoustic Reynolds stress terms were identified as the potential sources feeding the second modes. Shortly after, Chen, Wang & Fu (2022) used Chu’s energy norm definition to study hypersonic transition, including thermochemical effects. The stabilizing effect due to viscosity was pointed out to be dominant for the first modes while found to be negligible for the second modes. For transition in hypersonic flows, the mean flow is expected to have a profound effect on the disturbance energy growth. Thus, Myers’ energy norm definition, which includes the disturbance’s thermal, acoustic and kinetic energies along with the mean flow energy, offers a more complete description in this context. Instead of analysing the stability of perturbation growth via a modal energy norm, the current study uses the existing energy norm definitions as templates to describe the perturbation energy conservation in a non-modal form. Quantifying the spatiotemporal evolution of perturbation energy allows one to capture the nonlinear growth of second modes and identify the sources/sinks that destabilize/stabilize them. Such a framework can be translated to analyse the linear and nonlinear dynamics of other instability types involved in turbulence breakdown.

1.4. Manuscript outline

This work focuses on capturing the spatiotemporal evolution of second-mode instabilities and exploring the underlying dynamics. In the linear regime, high-fidelity DNS results have been compared with quasiparallel LST results. The derived perturbation energy corollary and the DNS results are then used to analyse the nonlinear regime. This study aims to close the energy budget of the second-mode instabilities and provide insight into the dynamics responsible for the production or destruction of instability energy.

Section 2 documents the flow under inspection and the governing equations that are simulated using DNS. The DNS set-up with the adopted noise generation approach is developed in § 3. In § 4, the second-order perturbation energy corollary is derived. The results of this study are broken into three sections. In § 5, the sensitivity of the second modes to different forcing parameters is examined. Section 6 presents the perturbation energy budget closure and discusses the production/dissipation mechanisms responsible for the amplification/attenuation of second modes. In § 7, the perturbation energy is presented in a non-dimensional form and a scaling analysis is performed. Finally, § 8 summarizes the findings.

2. Problem formulation

2.1. Flow conditions

The chosen canonical geometry is a sharp-tipped $3^\circ$ half-angle cone. The cone is exposed to a Mach $6$ flow at a $0^{\circ }$ angle of attack. The free stream parameters, denoted as $()_\infty$ , are tabulated in table 1. The flow parameters and the geometry are chosen based on the experimental runs at Purdue’s BAM6QT facility (Miller et al. 2023). The cone used experimentally has a length of $L=1.2\,\textrm{m}$ , allowing sufficient room to visualize the transition process. The same geometrical configuration has been maintained in the DNS runs performed in this work. The cone surface is modelled as an isothermal wall maintained at $T_w=300\,\textrm{K}$ . Air is assumed to be calorically perfect ( $C_p= \gamma R / (\gamma -1)$ , $C_v = C_p/\gamma$ ), with specific gas constant $R=287.113 \,J\,kg^{-1}\,K^{-1}$ and specific heat ratio $\gamma =1.4$ . The Prandtl number is assumed to stay constant at $Pr=0.707$ . Viscosity is modelled using Sutherland’s law with the low-temperature correction (Mack 1965).

Table 1. Free stream flow parameters used in the current work.

2.2. Governing equations

Throughout this manuscript, the two-dimensional BL coordinate system attached to the cone surface will be $(x,y)$ , denoting the streamwise and wall-normal coordinate, respectively. The Cartesian coordinate system is $(X,Y,Z)$ , with $X$ being the axis of symmetry of the cone. The azimuthal direction, wrapped around the cone’s surface, will be denoted as $\theta$ . The conservation equations in index notation read

(2.1a) \begin{gather} \frac {\partial \rho }{\partial t} + \frac {\partial }{\partial x_j} \left ( \rho u_j \right ) = 0, \end{gather}

(2.1b) \begin{gather} \frac {\partial \left ( \rho u_i \right )}{\partial t} + \frac {\partial }{\partial x_j} \left ( \rho u_i u_j + p\delta _{ij} - \tau _{ji} \right ) = 0, \end{gather}

(2.1c) \begin{gather} \frac {\partial \left ( \rho E \right )}{\partial t} + \frac {\partial \left ( \rho u_j E \right )}{\partial x_j} = \frac {\partial }{\partial x_j} \left ( \frac {\mu C_p}{Pr} \frac {\partial T}{\partial x_j} \right ) -\frac {\partial \left (u_i p\delta _{ij}\right )}{\partial x_j} + \frac {\partial \left (u_i \tau _{ji} \right )}{\partial x_j} , \end{gather}

(2.1d) \begin{gather} p = \rho R T, \end{gather}

where $(x_1,x_2,x_3) = (X,Y,Z)$ denotes the Cartesian system, which is not to be confused with the two-dimensional, BL-attached system denoted as $(x,y)$ . The effect of gravity is ignored, and other source terms are not considered in the governing equation. The total energy (per unit mass) $E$ is the sum of specific internal energy $e$ and specific kinetic energy (KE), $E=e + V^2/2$ , where $V$ is the norm of the velocity vector. The viscous stress tensor $\tau _{ji}$ is given as $\mu ( \partial u_j/\partial x_i + \partial u_i/\partial x_j ) + \lambda (\partial u_k/\partial x_k) \delta _{ij}$ , where $\lambda =-2/3\mu$ is used with the assumption that bulk viscosity $\mu _b=0$ (Stokes’ hypothesis) and $\delta _{ij}$ is the Kronecker delta.

2.3. Methodology

This work intends to focus on the hypersonic BL transition governed by second-mode waves only. Other transition mechanisms, like cross-flow instabilities, oblique breakdown, etc., have not been considered. Therefore, a fully resolved axisymmetric simulation is sufficient to capture the relevant physics. An artificial perturbation in the wall-normal velocity is used to initiate second-mode waves at different amplitude levels, ranging from linear to strongly nonlinear (while still axisymmetric). The stability of the base flow is also analysed using a spectral linear stability solver (Sousa et al. 2023) and compared with the DNS.

2.3.1. Direct numerical simulation: finite-difference Navier–Stokes solver, CFDSU

The fully compressible Navier–Stokes equations (2.1) are solved using a finite differencing code, CFDSU, originally written by Nagarajan, Lele & Ferziger (2003). It is a Fortran-based code that is continuously being developed at Purdue University. The code CFDSU has been used previously for various types of flow applications, which include large-eddy simulation modelling (Chapelier, Wasistho & Scalo 2018); Chen & Scalo (2021), vortex dynamics (Chapelier, Wasistho & Scalo 2019) and hypersonic transitional flow (Sousa et al. 2019, 2024).

The governing equations are solved on a structured curvilinear grid in a staggered arrangement. Thermodynamic variables, including density $\rho$ , temperature $T$ and pressure $p$ , are stored at the cell centres, while the other conservative variables involving velocities are stored at cell faces. Temporal integration has been performed using a fourth-order accurate, six-stage Runge–Kutta scheme (Allampalli et al. 2009). Spatial discretization is done using the compact sixth-order scheme (C6) by Visbal & Gaitonde (2002). This numerical strategy is expected to achieve spectral-like numerical accuracy with less numerical dissipation. This is desirable for resolving the relevant length and time scales involved in a hypersonic BL transition problem. A grid transformation is done to map the physical domain in Cartesian coordinates to a cylindrical coordinate system, where the governing equations are numerically solved. The curved grid in the cylindrical coordinate system (physical space) is further mapped to an equispaced orthogonal grid (computational space) using a curvilinear grid transformation suitable for applying the structured spatial scheme. The relevant equations for grid transformations involved are given in Nagarajan, Lele & Ferziger (2007).

2.3.2. Linear stability theory: Laguerre–Galerkin spectral linear stability solver

The stability of the BL flow has been analysed using a spectral LST solver, employing the Laguerre functions as the orthogonal basis. The linearized perturbation equations, obtained by linearizing the governing equations (2.1) about a steady base flow, can be posed as an eigenvalue problem. This requires applying the normal mode approximation and the locally parallel flow approximation, which imparts two types of eigenvalue problems: temporal and spatial. The spatial eigenvalue problem has been previously solved using the Chebyshev spectral collocation method (Malik 1990). Recently, Sousa et al. (2023) showed the suitability of the Laguerre–Galerkin spectral method for open BL stability problems. The Laguerre functions are defined on the positive real axis $[0,\infty )$ , and their natural property of decaying exponentially away from zero allows the boundary conditions of open shear flows to be automatically satisfied without a hard Dirichlet imposition. In the present work, the methodology outlined in Sousa et al. (2023) has been used to develop an in-house stability solver to solve the spatial eigenvalue problem.

3. Computational set-up

3.1. Computational domain and grid generation

Axisymmetric simulations are performed in a rectangular domain (see figure 1) attached to the cone surface. Low-order precursor runs, including the shock, are used to initialize the high-order DNS runs, which focus on the BL below the shock. The methodology adopted for base flow generation is described in appendix A. The BL edge quantities, denoted by the subscript “e”, are also reported in appendix A. The computational domain for the high-order DNS extends $x=0.12 - 1.2\,\textrm{m}$ in the streamwise direction and up to $y= 20\,\textrm{mm}$ in the wall-normal direction. Sponge layers are positioned at the inlet and outlet walls on the left- and right-hand boundaries, spanning approximately $8\%$ of the domain length in $x$ . The top wall sponge layer spans $20\%$ of the domain height. The role of the outlet sponge layers on the right-hand side and at the top is to ramp down the outgoing flow gradually to zero, thus preventing the reflection of spurious waves and suppressing numerical oscillations in the flow.

Figure 1. Schematic of the $\psi _c=3^\circ$ half-angle cone of length $L=1.2\,\textrm{m}$ . The computational domain is shown, which starts with an offset of $x_0=0.12\,\textrm{m}$ ( $10\%$ of $L$ ) from the cone tip and has a height of $20\,\textrm{mm}$ in $y$ . The different wall forcing strategies (detailed later in table 3 of § 5) are applied at the location $x_d=0.282\,\textrm{m}$ (at $15\,\%$ of the domain length) across a strip of width $w_d=0.054\,\textrm{m}$ ( $5\%$ of the domain length).

The computational domain has been discretized uniformly in the streamwise direction $x$ , while a stretched grid is used in the wall-normal direction $y$ . A half-tangent hyperbolic stretching is employed in $y$ , such that the coarsest grid has more than $\sim 50$ points within the BL (see figure 1), allowing for accurate resolution of the acoustic-wave-induced Stokes’ BL. Simulations have been carried out on three grid levels ( $6144 \times 256$ , $9216 \times 384$ and $12\,288 \times 512$ ). Grid sensitivity analysis shows the intermediate grid $9216 \times 384$ to be appropriate for the scope of this study (see § 5.1).

The flow is sustained at the inlet via a sponge layer. Homogeneous Neumann conditions are used for all outlet flow quantities. No-slip, no-penetration conditions are imposed on the bottom wall, i.e., the cone surface. Velocity and temperature fluctuations are set to zero at the wall. After a statistically stationary state is reached, the forcing pulse is introduced, and calculations are run for one flow-through time of $1.3\,\textrm{ms}$ .

3.2. Wall forcing strategies

Transition modelling through DNS involves imposing a randomized disturbance as an initial condition to trigger transition. The spectral content of this initial excitation significantly dictates the transition process, especially the onset of the nonlinear phase. Fedorov & Khokhlov (2002) found the hypersonic BL to be more receptive to wall vibrations and blowing/suction type forcing than temperature fluctuations-based forcing. Malik, Lin & Sengupta (1999) demonstrated the transition process to be independent of which state variables the disturbance is forced on, as long as the signal characteristics stay consistent. Thus, one should be careful when introducing numerical noise in transition studies to make a fair comparison with flight or wind tunnel data.

The initial disturbance in the axisymmetric DNS runs is introduced in a strip of width $w_d$ centred at $x_d$ (see figure 1). The noise is introduced at the cone surface as wall-normal velocity fluctuations. Two separate strategies have been used, as described in the following sections.

3.2.1. Blowing/suction wall forcing

A blowing/suction type wall excitation was first proposed by Gaster & Grant (1975) for low-speed transition, which was used later on by Sivasubramanian & Fasel (2014) and Sousa et al.2023 for high-speed transition simulations. The resulting wall-normal velocity fluctuations read

(3.1) \begin{equation} v'_w (x, y=0, t) = A_0 \cos ^3 \left ( \pi \xi \right ) \sin \left (2 \pi f_d t \right ), \quad \forall t \in [0,1/f_d], \end{equation}

where $A_0$ is the maximum amplitude of the pulse and $\xi = ( x - x_d )/w_d$ , such that $\xi \in [ -1/2,1/2 ]$ , where $x_d$ and $w_d$ are portrayed in figure 1.

A finite pulse triggers a wide spectrum of frequencies with no direct control over the frequency spectrum. In the following section, a second forcing method is developed where pseudorandom noise is projected on a spectral basis and filtered, guaranteeing grid independence and control over the spectral distribution.

3.2.2. Legendre-space filtered wall forcing

A univariate pseudorandom noise field $\phi (x_i)$ is sampled on a uniformly spaced grid $x_i$ in the streamwise direction and then projected onto the Legendre spectral space. Therein, the high-wavenumber modes that are unresolvable on the coarsest grid considered in this work are removed (see figure 2). The thus-obtained signal comprising a collection of Legendre modes is transformed back to the physical space with the same analytical expression regardless of grid refinement.

Figure 2. Schematic illustration of the Legendre-space filtered forcing strategy. (a) Sharp filtering of the sampled noise in the Legendre space. The horizontal axis denotes the Legendre mode number $n$ . (b) The effect of the filtering process in the physical domain. The initially sampled noise is shown in black, while the analytical noise obtained upon filtering is displayed in red. The black circles denote the Legendre–Gauss nodes where the discrete noise is initially sampled.

The aim of the modulation operation is to obtain an analytically filtered noise $\bar {\phi }(x)$ that can be evaluated at any point in space and stays invariant with grid refinement. The bar signifies the modulated function in the Legendre space throughout the following formulation. The Legendre-space modulation operator $\mathcal {M}_L$ is given as

(3.2) \begin{equation} \bar {\phi }(x)= \mathcal {M}_L\left (\sigma ^L_n, \phi \right ) = \mathcal {L}^{-1} \left [ \sigma ^L_n \mathcal {L}(\phi (x_i)) \right ], \quad \sigma ^L_n = \begin{cases} 1,\quad \forall n \leq N_{cL}\\ 0,\quad \text {otherwise} \end{cases}, \end{equation}

which will be evaluated at the DNS grid points $x=x_i^g$ for different grid refinement levels. These can be different from the sampled grid points $x_i$ . Here $\sigma ^L_n$ denotes a sharp Legendre-space filter, removing the noise content above a user-defined cutoff mode number $N_{cL}$ . Here $\mathcal {L}$ and $\mathcal {L}^{-1}$ denote the forward and backward discrete Legendre transform operations, respectively,

(3.3) \begin{equation} \tilde {\phi }_n = \mathcal {L}\left [\phi (x_i)\right ] = \frac {1}{\gamma _n}\sum _{i=0}^{N_x} \phi (x_i)L_n(x_i)w_i, \quad \bar {\phi }(x) = \mathcal {L}^{-1}\left [\tilde {\phi }_n\right ] = \sum _{n=0}^{N_{cL}} \tilde {\phi }_n L_n(x), \end{equation}

where, $\gamma _n= ({2}/{2n+1})$ , and $w_i={2}/{ ((1-x_i^2 ) ( L'_{N_x+1}(x_i) )^2)}$ denotes the nodal weights for Legendre–Gauss quadrature (Shen, Tang & Wang 2011). While $N_x$ varies depending on the grid size, $N_{cL}$ is kept constant to achieve grid independence.

Modulation can be explored in other spectral spaces, like Chebyshev, Hermite, etc. However, in this work, a Legendre-space modulation has been chosen due to the more uniform distribution of Legendre–Gauss quadrature nodes over the excitation region. The Chebyshev–Gauss quadrature nodes tend to be more aggressively concentrated at the ends.

Modulation in time is achieved via a series of $m_f$ discrete frequencies that are intended to be introduced into the BL. For each frequency, a different noise field $\bar {\phi }$ is created. The final wall-normal velocity noise signal is a summation of all the $m_f$ reconstructed fields and has the form

(3.4) \begin{equation} v'_w (x, y=0, t) = A_0 \sum _{m=0}^{m_f} c(f_m) \cdot \Re \left ( \bar {\phi }(x, f_m) e^{-2\pi i f_m t} \right ), \end{equation}

where $c(f_m)$ is a frequency-dependent amplitude modulation that can be used to decide the colour of the noise, for e.g., pink, white, brown, etc. Pink noise, which is generally accepted to mimic ambient noises, has been used in this work. Here ‘ $A_0$ ’ is the amplitude parameter that indicates the strength of the introduced signal and is used for all the analyses in this work. The forcing has been applied over the time span $t \in [0,n_{c}/f_0 ]$ in the simulations, where $f_0$ is the lowest frequency and $n_{c}$ represents the number of periodic cycles in time.

4. Theoretical disturbance energy derivation

Transition to turbulence under low free stream disturbance amplitudes and on slender canonical geometries entails an initial modal amplification phase, followed by nonlinear growth, saturation and breakdown phases. In this paper, we also consider the scenario where perturbations are initialized directly with high amplitudes to explore the nonlinearities in second-mode growth through an energy budget analysis.

The derivation of a new disturbance energy equation below is discussed by relying on the energy corollary definitions given by Myers (1991) and Chu (1965). For the sake of generality, the derivation uses the governing equations in the Cartesian coordinates $(x_j=X,Y,Z)$ . The derived equations can be directly applied to the BL attached system $(x,y)$ used in this work (see figure 1), invoking the parallel flow assumption as is typically done in an LST analysis. The conservation equations are first recast with the entropy evolution equation instead of total energy $E$ , yielding,

(4.1a) \begin{equation} \frac {\partial \rho }{\partial t} + \frac {\partial }{\partial x_j} \left ( \rho u_j \right ) = 0, \end{equation}

(4.1b) \begin{equation} \frac {\partial }{\partial t} \left ( \rho u_i \right ) + \frac {\partial }{\partial x_j} \left ( \rho u_i u_j + p\delta _{ij} - \tau _{ji} \right ) = 0, \end{equation}

(4.1c) \begin{equation} \rho T \left ( \frac {\partial s}{\partial t} + u_j \frac {\partial s}{\partial x_j} \right ) = \frac {\partial }{\partial x_j} \left ( \frac {\mu C_p}{Pr} \frac {\partial T}{\partial x_j} \right ) + \tau _{ji} \frac {\partial u_i}{\partial x_j}, \end{equation}

(4.1d) \begin{equation} p = \rho R T, \end{equation}

where the ideal gas law has been included for completeness. The viscous stress tensor $\tau _{ji}$ is defined previously in § 2.2. The flow field variables are then decomposed into a stationary base state, denoted as $(.)_0$ , and a fluctuating component about the base state representing the disturbance, denoted as $(.)'$ :

(4.5) \begin{equation} \left [ u_{i},\rho ,T,p,s\right ]^T = \left [ u_{0i},\rho _0,T_0,p_0,s_0\right ]^T + \left [ u'_{\!i},\rho ',T',p',s'\right ]^T. \end{equation}

The fluctuation terms are not expanded asymptotically, as was done by Myers (1991) since the resulting functional form of the fluctuations, divided based on their orders of magnitude, would not offer physical insight into the nonlinear interactions between different physical quantities. Substituting the decomposition (4.5) into (4.1) and subtracting the base governing equations (4.1) $_0$ , one obtains the full, nonlinear equations governing the fluctuation field:

(4.6a) \begin{equation} \frac {\partial \rho '}{\partial t} + \frac {\partial }{\partial x_j} \left ( \rho _0 u'_{\!\!j} + \rho ' u_{0j} + \underbrace {\rho ' u'_{\!\!j}}_{\textit {O}^{(2)}} \right ) = 0, \end{equation}

(4.6b) \begin{equation} \begin{aligned} &\frac {\partial }{\partial t} \left ( \rho _0 u'_{\!i} + \rho ' u_{0i} + \underbrace {\rho ' u'_{\!i}}_{\textit {O}^{(2)}} \right ) + \frac {\partial }{\partial x_j} \left ( \rho _0 u_{0i}u'_{\!\!j} + \rho _0 u'_{\!i} u_{0j} + \rho ' u_{0i}u_{0j} \right. \\ &\left. + \underbrace {\rho _0 u'_{\!i} u'_{\!\!j} + \rho ' u_{0i}u'_{\!\!j} + \rho ' u'_{\!i} u_{0j} }_{\textit {O}^{(2)}} + \underbrace {\rho ' u'_{\!i} u'_{\!\!j}}_{\textit {O}^{(3)}} + p' \delta _{ij} - \tau'_{\!\!\!ji} \right ) = 0, \end{aligned} \end{equation}

(4.6c) \begin{equation} \begin{aligned} \rho _0 T_0 \left ( \frac {\partial s'}{\partial t} + u_{0j} \frac {\partial s'}{\partial x_j} + u'_{\!\!j} \frac {\partial s_0}{\partial x_j} \right ) &+ \left ( \rho ' T_0 + \rho _0 T' \right ) u_{0j} \frac {\partial s_0}{\partial x_j} \\ + \underbrace {\left ( \rho ' T_0 + \rho _0 T' \right ) \left ( \frac {\partial s'}{\partial t} + u_{0j} \frac {\partial s'}{\partial x_j} + u'_{\!\!j} \frac {\partial s_0}{\partial x_j} \right )}_{\textit {O}^{(2)}} &+ \underbrace {\rho _0 T_0 u'_{\!\!j} \frac {\partial s'}{\partial x_j}}_{\textit {O}^{(2)}} + \underbrace {\rho ' T' u_{0j}\frac {\partial s_0}{\partial x_j}}_{\textit {O}^{(2)}} \\ + \underbrace {\rho ' T' \frac {\partial s'}{\partial t}}_{\textit {O}^{(3)}} + \underbrace {\left ( \rho ' T_0 + \rho _0 T' \right ) u'_{\!\!j} \frac {\partial s'}{\partial x_j}}_{\textit {O}^{(3)}} &+ \underbrace {\rho ' T' \left ( u_{0j}\frac {\partial s'}{\partial x_j} + u'_{\!\!j} \frac {\partial s_0}{\partial x_j} \right )}_{\textit {O}^{(3)}} \\ &+ \underbrace {\rho ' T' u'_{\!\!j} \frac {\partial s'}{\partial x_j}}_{\textit {O}^{(4)}}\\ = \frac {\partial }{\partial x_j} \left ( \frac {\mu _0 C_p}{Pr} \frac {\partial T'}{\partial x_j} + \frac {\mu ' C_p}{Pr} \frac {\partial T_0}{\partial x_j} \right ) + \tau'_{\!\!\!ji}\frac {\partial u_{0i}}{\partial x_j} &+ \tau _{0ji} \frac {\partial u'_{\!i}}{\partial x_j} + \underbrace {\tau'_{\!\!\!ji} \frac {\partial u'_{\!i}}{\partial x_j}}_{\textit {O}^{(2)}} + \underbrace {\frac {\mu ' C_p}{Pr} \frac {\partial T'}{\partial x_j}}_{\textit {O}^{(2)}}, \end{aligned} \end{equation}

(4.6d) \begin{equation} p' = \rho ' RT_0 + \rho _0 RT' + \underbrace {\rho ' RT'}_{\textit {O}^{(2)}}. \end{equation}

Higher orders of nonlinearity in the above equations are denoted with underbraced text below the respective terms. Equation (4.6b ) denotes three separate equations, referred hereafter as (4.6b ) $_i$ whenever necessary. Obtaining the equation for disturbance energy evolution further requires the derivation of two corollaries from the above system of equations (4.6).

The first corollary is derived from a combination of conservation of mass and momentum and comprises a set of three equations, which follows from (4.6b ) $_i - u_{0i}\times$ (4.6a ) $ - u'_{\!i}\times$ (4.1a ) $_0$ , yielding

(4.10) \begin{equation} \begin{aligned} \rho _0 \frac {\partial u'_{\!i}}{\partial t} &+ \frac {\partial }{\partial t} \underbrace {\left ( \rho ' u'_{\!i} \right )}_{\textit {O}^{(2)}} + \left ( \rho _0 u'_{\!\!j} + \rho ' u_{0j} + \underbrace {\rho ' u'_{\!\!j}}_{\textit {O}^{(2)}} \right ) \frac {\partial u_{0i}}{\partial x_j} + \rho _0 u_{0j} \frac {\partial u'_{\!i}}{\partial x_j} \\ &+ \frac {\partial }{\partial x_j} \left ( \underbrace {\rho _0 u'_{\!i} u'_{\!\!j} + \rho ' u'_{\!i} u_{0j}}_{\textit {O}^{(2)}} + \underbrace {\rho ' u'_{\!i} u'_{\!\!j}}_{\textit {O}^{(3)}} + p'\delta _{ij} - \tau '_{ji} \right ) = 0. \end{aligned} \end{equation}

The second step is casting the $\partial \rho '/\partial t$ term in (4.6a ) in terms of pressure and entropy derivatives. This involves using the total differential ${\textrm{d}}\rho = ( \partial \rho/\partial p)_s {\textrm{d}}p + ( \partial \rho/\partial s)_p {\textrm{d}}s$ , which yields the thermodynamic relation $a^2 \textrm{d}\rho = \textrm{d}p - \rho \beta T a^2(C_p)^{-1} \textrm{d}s$ . The ideal gas assumption is further invoked, i.e., $\beta T=1$ . This yields the equation

(4.11) \begin{equation} \begin{aligned} \frac {\partial p'}{\partial t} &+ a_0^2 \frac {\partial }{\partial x_j} \left ( \rho _0 u'_{\!\!j} + \rho ' u_{0j} \right ) + \underbrace {a_0^2 \frac {\partial }{\partial x_j} \left ( \rho ' u'_{\!\!j} \right )}_{\textit {O}^{(2)}} + \underbrace {\gamma RT' \frac {\partial }{\partial x_j}\left ( \rho ' u'_{\!\!j} \right )}_{\textit {O}^{(3)}} \\ &+ \left ( \gamma -1 \right ) \left [ \rho _0T_0\left ( u_{0j}\frac {\partial s'}{\partial x_j} + u'_{\!\!j}\frac {\partial s_0}{\partial x_j} \right ) + \left ( \rho 'T_0 + \rho _0 T' \right )u_{0j}\frac {\partial s_0}{\partial x_j} \right ] \\ &- \left ( \gamma -1 \right ) \left [ \frac {\partial }{\partial x_j}\left ( \frac {\mu _0 C_p}{Pr}\frac {\partial T'}{\partial x_j} \right ) + \frac {\partial }{\partial x_j}\left ( \frac {\mu ' C_p}{Pr}\frac {\partial T_0}{\partial x_j} \right ) + \tau'_{\!\!\!ji}\frac {\partial u_{0i}}{\partial x_j} + \tau _{0ji}\frac {\partial u'_{\!i}}{\partial x_j} \right ] \\ &+ \underbrace {\left ( \gamma -1 \right ) \left [ \left ( \rho 'T_0+\rho _0T' \right )\left ( u_{0j}\frac {\partial s'}{\partial x_j} + u'_{\!\!j}\frac {\partial s_0}{\partial x_j} \right ) + \rho _0T_0u'_{\!\!j}\frac {\partial s'}{\partial x_j} + \rho 'T'u_{0j}\frac {\partial s_0}{\partial x_j} \right ]}_{\textit {O}^{(2)}} \\ &- \left ( \gamma -1 \right ) \left [ \underbrace {\tau'_{\!\!\!ji}\frac {\partial u'_{\!i}}{\partial x_j} + \frac {\partial }{\partial x_j}\left ( \frac {\mu 'C_p}{Pr}\frac {\partial T'}{\partial x_j} \right )}_{\textit {O}^{(2)}} + \underbrace {\rho 'T'U_j'\frac {\partial s'}{\partial x_j}}_{\textit {O}^{(4)}} \right ]\\ &+ \left ( \gamma -1 \right ) \left [ \underbrace {\left ( \rho 'T_0+\rho _0T' \right )u'_{\!\!j}\frac {\partial s'}{\partial x_j} + \rho 'T'\left ( u_{0j}\frac {\partial s'}{\partial x_j}+ u'_{\!\!j}\frac {\partial s_0}{\partial x_j} \right )}_{\textit {O}^{(3)}} \right ] = 0. \end{aligned} \end{equation}

The above set of equations govern the evolution of the disturbance field, from which a disturbance energy conservation equation can be recovered as a corollary. In the following section, we discuss well-established energy definitions from classical mechanics and thermodynamics that support the functional form of the disturbance energy (4.16) chosen in this work, acknowledging that this choice is not uniquely determined.

4.1. Disturbance energy definition

Transition to turbulence is caused by the gradual amplification of disturbances, which represent fluctuations around the steady, laminar base state. Such fluctuating fields evolve in space and time by drawing energy from (and, in some cases, distorting) the base state. This section presents a suitable definition of the disturbance energy, which should account for the increase or decrease of the total instantaneous energy of a fluid element due to the presence of a disturbance wave packet. The total energy of a fluid element constitutes the KE due to motion and potential energy(PE) due to its configuration. The enthalpy, which represents the work and heat interactions across the element’s boundaries, is used to obtain the PE.

The disturbance energy, thus defined, can then be recovered from the perturbation governing equations (4.6) through a suitable linear combination, leading to a disturbance energy conservation equation.

4.1.1. Kinetic energy

For a rigid body of fixed mass $m$ , acted upon by a force $F_i=m \textrm{d}u_i/\textrm{d}t$ displacing it by $\textrm{d}x_i=u_i \textrm{d}t$ , the total work done is $\int \textrm{d}W=\int mdu_i/\textrm{d}t \cdot u_i \textrm{d}t=\int _0^{u_i} mu_i \,\textrm{d}u_i$ , which gets stored as KE.

The KE (per unit volume) of the perturbed flow can be similarly defined on a fluid element as the work done by the force (per unit volume) $\textrm{d}(\rho u_i)/\textrm{d}t$ to displace such element by $\textrm{d}x_i=u_i \textrm{d}t$ . Thus, the total KE of the flow (per unit volume) can be given as

(4.12) \begin{equation} \begin{aligned} E_k &= \frac {1}{2} \rho u_i^2 = \frac {1}{2} \left ( \rho _0+\rho ' \right ) \left ( u_{0i} + u'_{\!i} \right )^2 \\ &= \underbrace {\frac {1}{2} \rho _0 u_{0i}^2}_{\textit {O}^{(0)}} + \underbrace {\rho _0 u_{0i} u'_{\!i} + \frac {1}{2} \rho ' u_{0i}^2}_{\textit {O}^{(1)}} + \underbrace {\frac {1}{2} \rho _0 u_i^{\prime 2} + \rho ' u_{0i} u'_{\!i}}_{\textit {O}^{(2)}} + \underbrace {\frac {1}{2} \rho ' u_i^{\prime 2}}_{\textit {O}^{(3)}}. \end{aligned} \end{equation}

The zeroth-order term represents the base flow KE, while the following two terms constitute the first-order fluctuation energy, comprising the perturbation momentum carried by the non-zero mean flow. Finally, the second-order terms obtained here form the KE component of the disturbance energy expression (4.16).

4.1.2. Acoustic and thermal potential energy

Due to the presence of a disturbance, the mechanical work done on a fluid element to compress it gets stored as its PE, referred to as the acoustic PE in the literature (Kinsler et al. 2000). In a constant-volume fluid element, the compression/expansion work can be quantified by the enthalpy fluctuation when no heat exchange occurs (isentropic). The second-order acoustic PE (per unit volume) then reads

(4.13) \begin{equation} E_a= \int (\rho ^{\prime} dh^{\prime})_s = \int \left (\rho ^{\prime} \frac {\partial h}{\partial p}\right )_s \textrm{d}{p}^{\prime} = \int _0^{p^\prime} \frac {p^\prime}{\gamma p_0} \,\textrm{d}{p^\prime} = \frac {1}{2} \frac {p^{\prime 2}}{\gamma p_0}, \end{equation}

where the expression $ (\rho^{\prime}/\rho _0=p^{\prime}/\gamma p_0 - s^{\prime}/C_p )$ is used to replace $\rho '$ . The above equation indicates that if a control element is compressed in the presence of a propagating disturbance, its acoustic PE increases while it decreases during an expansion. This can be alternatively derived considering a Lagrangian control-mass element, whose volume changes as it moves with the flow. The isentropic $pdV$ work on such an element reads

(4.14) \begin{equation} {E_a} = - \int _{V_0}^V \frac {p^{\prime}}{V} \,\textrm{d}V = \int _0^{\rho ^{\prime}} \frac {p^{\prime}}{\rho } \,\textrm{d}\rho ^{\prime} = \int _0^{p^{\prime}} \frac {p^{\prime}}{\rho } \,\frac {\textrm{d} p^{\prime}}{c^2} = \frac {1}{2}\frac {p^{\prime 2}}{\gamma p_0}, \end{equation}

which is equivalent to (4.13). The negative sign above indicates a PE increase (work done on the element) during compression caused by a positive acoustic pressure $p^{\prime}$ (Kinsler et al. 2000).

Similarly, the heat added to a fluid element gets stored as the thermal PE. On a constant volume element, this heat addition can be quantified by the enthalpy fluctuation when no compression/expansion work is done to cause pressure change (isobaric). This gives the second-order thermal PE as

(4.15) \begin{equation} E_t = -\int (\rho^{\prime} dh^{\prime})_p = -\int \left (\rho^{\prime} \frac {\partial h}{\partial s}\right )_p \textrm{d}s^{\prime} = \int _0^{s^{\prime}} \frac {\rho _0 T_0}{C_p} s^{\prime} \,\textrm{d}s^{\prime} = \frac {\rho _0 T_0}{2C_p}s^{\prime 2}. \end{equation}

The negative sign is needed to associate a density decrease ( $\rho^{\prime}\lt 0$ ) at constant pressure with a heat gain, which indicates an increase in thermal PE. Heat gain at constant pressure occurs due to the influx of high-temperature ( $T^{\prime}\gt 0$ ) fluid into the control volume.

4.1.3. Deriving the energy equation

Finally, combining all of the aforementioned components that make up the total energy of a perturbation to the second order, yields

(4.16) \begin{equation} E_2 = E_t + E_k + E_a = \underbrace {\frac {1}{2} \frac {\rho _0 T_0}{C_p}s^{\prime 2}}_{i} + \underbrace {\rho^{\prime} u_{0i}u_i^{\prime}}_{ii} + \underbrace {\frac {1}{2} \rho _0 u_i^{\prime 2}}_{iii}+ \underbrace {\frac {1}{2}\frac {p^{\prime 2}}{\gamma p_0}}_{iv}. \end{equation}

Throughout this manuscript, each term in the above expression is referred to by the numerals given below. Correctly deriving the disturbance energy equation entails choosing an appropriate combination of (4.6c ), (4.10) and (4.11) to obtain $E_2$ , which reads

(4.17)

The result of this combination contains perturbation terms up to the fifth order of nonlinearity because of (4.6c ) and (4.11). By truncating up to the second order, the desired perturbation energy conservation equation can be recovered. Contrary to Myers’ naming convention, it will be denoted as the second-order energy equation hereafter.

The multipliers in (4.17) are required to achieve the specific perturbation energy $(E_2)$ chosen in the current work (4.16). While its expression is justified by the considerations in §§ 4.1.1 and 4.1.2, we note that one might choose to use a different expression for the perturbation energy. For example, George & Sujith (2012) offers an expression for a family of allowable perturbation norms. The question of the uniqueness of the perturbation energy of a fluctuating thermofluid-dynamic field remains a topic for continued academic investigation.

4.2. Second-order disturbance energy equation

The above-specified combination (4.17) to the second order yields a conservation equation of the disturbance energy having the form

(4.18) \begin{equation} \frac {\partial E_2}{\partial t} + \frac {\partial I_{2j}}{\partial x_j} = D_2, \end{equation}

where the second-order disturbance energy, $E_2$ , is given in (4.16). Term $i$ in (4.16) denotes the thermal part of the disturbances, as entropy fluctuations represent the influx/outflux of thermal energy due to the propagating disturbance. Terms $ii$ and $iii$ are second-order components of the total KE, which, for analysis purposes in this work, will be referred to as the disturbance energy carried by the mean flow and the disturbance KE, respectively. Finally, term $iv$ quantifies the energy content due to compression/rarefaction effects. However, the term $\rho^{\prime}u_{0i}u_i^{\prime}$ is not strictly positive definite, which does not allow (4.16) to be considered as a norm. Instead, the energy definition is used to develop a perturbation energy budget device and analyse the sources of perturbation growth or decay. A discussion regarding the positive definiteness of this term and the reasoning behind its inclusion in the perturbation energy definition is presented in appendix B.1.

The energy flux term $I_{2j}$ is given as

(4.19) \begin{equation} I_{2j} = \underbrace {\frac {1}{2} \frac {\rho _0 T_0}{C_p}s^{\prime 2} u_{0j}}_{i} + \underbrace {\left ( \rho _0 u_{0i}u_i^{\prime} + p^{\prime} \right ) \left ( u_j^{\prime} + \frac {\rho^{\prime} u_{0j}}{\rho _0} \right )}_{ii} + \underbrace {\frac {1}{2} \rho _0 u_i^{\prime 2} u_{0j}}_{iii}, \end{equation}

where terms $i,\,ii,\,iii$ represent the convective fluxes of the thermal, acoustic and kinetic disturbance energies, respectively. Finally, the source/sink term $D_2$ comes out to be

(4.20) \begin{equation} \begin{aligned} &D_2 = - \underbrace {\frac {p^{\prime} u_{0j}}{C_p} \frac {\partial s^{\prime}}{\partial x_j}}_{I} + \underbrace {\left ( \frac {\gamma -1}{\gamma }\frac {p^{\prime}}{p_0} + \frac {s^{\prime}}{C_p} \right ) \frac {\partial }{\partial x_j} \left ( \frac {\mu _0 C_p}{Pr} \frac {\partial T^{\prime}}{\partial x_j} + \frac {\mu ^{\prime} C_p}{Pr} \frac {\partial T_0}{\partial x_j} \right )}_{III} \\ & - \underbrace {\rho^{\prime} u_{0i}\left ( u_j^{\prime} + \frac {\rho^{\prime} u_{0j}}{\rho _0} \right )\frac {\partial u_{0i}}{\partial x_j}}_{II} - \underbrace {\left ( \frac {\gamma -1}{\gamma }\frac {p^{\prime}}{p_0} + \frac {s^{\prime}}{C_p} \right )\left ( \rho _0 T_0 u_j^{\prime} + \rho^{\prime} T_0 u_{0j} + \rho _0 T^{\prime} u_{0j} \right )\frac {\partial s_0}{\partial x_j}}_{IV}\\ & + \underbrace {\left ( u_i^{\prime}+\frac {\rho^{\prime} u_{0i}}{\rho _0} \right ) \frac {\partial \tau _{ji}^{\prime}}{\partial x_j}}_{V} - \underbrace {\frac {p^{\prime}}{\rho _0}\left ( u_j^{\prime} + \frac {\rho^{\prime} u_{0j}}{\rho _0} \right ) \frac {\partial \rho _0}{\partial x_j}}_{VI} + \underbrace {\rho _0 u_{0i}u_j^{\prime} \frac {\partial u_i^{\prime}}{\partial x_j}}_{VII} + \underbrace {s^{\prime 2} \frac {\partial }{\partial x_j} \left ( \frac {\rho _0 T_0 u_{0j}}{C_p} \right )}_{VIII} \\ & + \underbrace {\left ( \frac {\gamma -1}{\gamma } \frac {p^{\prime}}{p_0} + \frac {s^{\prime}}{C_p} \right ) \left ( \tau ^{\prime}_{ji}\frac {\partial u_{0i}}{\partial x_j} + \tau _{0ji} \frac {\partial u_i^{\prime}}{\partial x_j} \right )}_{IX}, \end{aligned} \end{equation}

where, $\tau ^{\prime}_{ji} = \mu _0 ( \partial u_j^{\prime}/\partial x_i + \partial u_i^{\prime}/\partial x_j ) - 2/3\mu _0 (\partial u_k^{\prime}/\partial x_k) \delta _{ij}$ . The terms in $D_2$ are ordered based on their relative importance. A brief explanation of the terms is given here, while a detailed breakdown of the terms will be given in § 6. Terms I and IV indicate an interaction between thermal and acoustic fluctuations, while term VIII only represents the effects of entropy fluctuations. Terms V and IX represent the effects of viscous diffusion, while term III represents thermal diffusion effects. Finally, the terms II, VI and VII define the source/sink effects influenced by the KE of the disturbance and the energy carried by the mean flow. It can be seen that while we eventually use the same disturbance energy definition as Myers’, the flux and the source terms (right-hand side) of our proposed budget equation are functionally different. The current disturbance energy equation has been compared with Myers’ first-order corollary in appendix B.

The objective of this manuscript is, in fact, to find an appropriate functional form of the overall disturbance energy conservation equation to identify the terms responsible for the production/dissipation of disturbance energy during the streamwise propagation of the second-mode waves. The terms in (4.20) unveil the physical mechanisms controlling the second-mode-dominated transition in a hypersonic BL.

5. Sensitivity of second-mode dynamics to forcing type and amplitude

In this work, we drive second-mode waves to nonlinear amplitudes and examine the resulting effects on their perturbation energy budgets. This section sets up the parametric study performed based on the tunable forcing parameters and briefly highlights the separation of linear and nonlinear second mode dynamics. Four amplitude levels, spanning three decades, are used (see table 2). The lowest amplitude chosen shows a linear behaviour, while the upper limit has been chosen based on the maximum achievable DNS resolution. An amplitude of $A_0=10\,\textrm{m}\,\textrm{s}^-{^1}$ leads to the generation of small-scale shocks, yielding numerical resolution issues. Therefore, the four gradually increasing amplitudes chosen for this study will be denoted as linear $ (10^{-3} )$ , weakly nonlinear $ (10^{-2} )$ , nonlinear $ (10^{-1} )$ and strongly nonlinear $ (10^0 )$ throughout the manuscript, as given in table 2.

Table 2. Amplitudes of the wall-normal velocity fluctuations imposed at the forcing region. The velocity amplitudes, non-dimensionalized by the free stream velocity ( $A_0/U_\infty$ ), are also shown. A gradual shift from linear to nonlinear fluctuation levels is intended.

Figure 3. Numerical schlieren $|\nabla \rho |$ showing rope-wave structures (Demetriades 1974), obtained from the DNS runs with Legendre forcing of amplitude $A_0=0.1\,\textrm{m}\, \textrm{s}^-{^1}$ at times $t=0.365\,\textrm{ms}$ (a) and $t=0.805\,\textrm{ms}$ (b).

Figure 4. The DNS schlieren showing the second-mode disturbance at times $t=0.365\,\textrm{ms}$ (a) and $t=0.805\,\textrm{ms}$ (b), with Legendre forcing amplitude $A_0=1\,\textrm{m}\, \textrm{s}^-{^1}$ .

A sign of nonlinearity is often attributed to the occurrence of rope-like structures (Demetriades 1974); Stetson & Kimmel (1993). Figure 4 shows the presence of similar structures for the higher amplitude forcings. These structures occur at the critical layer where the magnitude of the density gradient becomes high enough to trigger nonlinear action, resulting in wave distortion and steepening effects.

The different forcing types used to trigger transition are given in table 3. Two parameters are varied in the Legendre spectral space-based forcing – the spectral content of the disturbance and the time period over which the forcing is applied. This is done since nonlinear action is not only dependent on the initial forcing amplitude but also on the initial forcing characteristics, which are intended to be investigated in this work. The following types are considered:

Table 3. Forcing types considered in this work (see § 3.2 for detailed description). All forcings are applied over the streamwise region $x \in [0.255,0.309]$ m (see figure 1). $n_c$ denotes the number of cycles over which the forcing is applied.

1. (i) BL5– a broad-band Legendre forcing, applied with discrete frequencies ranging from $60-200\,\textrm{kHz}$ at intervals of $10\,\textrm{kHz}$ , applied over five full cycles of the lowest frequency $60\,\textrm{kHz}$ ;

2. (ii) L5– a monotone Legendre forcing, applied over five cycles of the frequency $120\,\textrm{kHz}$ ;

3. (iii) L1– a monotone Legendre forcing, applied over one cycle of the frequency $120\,\textrm{kHz}$ ;

4. (iv) C1– a blowing/suction pulse, applied over onecycle of the frequency $120\,\textrm{kHz}$ .

Here BL5 has been used as the primary forcing type for most analyses in this work since it is potentially the most nonlinear forcing type.

5.1. Power spectra of wall pressure fluctuations

The spatiotemporal growth of the second-mode waves has been tracked using the power spectral density (PSD) of the wall pressure fluctuations. Figure 5(b) shows the wall pressure fluctuations at three streamwise locations, while figure 5(a ) shows the PSD at the respective locations. The power spectra are normalized by the base wall pressure at the corresponding locations, i.e., $\bar {p}=p_0(x,y=0)$ . The PSD results reveal the frequency content of the instabilities responsible for triggering transition. Maximum spectral content is concentrated at frequencies of $110-150\,\textrm{kHz}$ , which is the typical range expected from the second-mode waves for a Mach $6$ sharp cone. The presence of higher-order modes becomes prominent with increasing amplitude levels, which is considered to be an indication of nonlinear dynamics (Stetson 1988).

Figure 5. Spatiotemporal evolution of the second-mode disturbance (b) and its power spectra (a), for increasing BL5 forcing amplitudes, $A_0$ (row-wise) (described in table 3). The PSDs are shown at streamwise locations $x=55\,\textrm{cm}$ (), $75\,\textrm{cm}$ () & $95\,\textrm{cm}$ (). Three increasing grid refinement levels are used: $6144\times 256$ (dotted), $9216\times 384$ (dashed) and $12\,288\times 512$ (solid).

Second-mode growth for the relatively lower forcing amplitudes, $A_0=10^{-3}$ and $10^{-2}\,\textrm{m}\,\textrm{s}^-{^1}$ , are observed to yield a similar trend, consistent with the linear dynamics. The wave packets widen and grow as they propagate downstream, with the frequency of the primary unstable mode gradually decreasing. The frequencies are independent of the disturbance amplitude levels close to the forcing. Farther downstream, the frequency shift and wave packet broadening are more significant for higher disturbance levels. Both the wave packet and the corresponding power spectra evolution for the highest amplitude (figure 5 a iv, b iv) show plateauing of the second-mode growth, which is not observed in the linear cases.

Table 4. The three grid levels considered for the grid sensitivity analysis. Numbers of grid points in $x$ and $y$ are tabulated as $N_x$ and $N_y$ . The uniform grid spacing in $x$ and the minimum spacing of the stretched grid along $y$ are shown as $\triangle x$ and $\triangle y_{wall}$ , respectively.

5.1.1. Grid sensitivity

The results of grid sensitivity analysis with the forcing type BL5 are also included in figure 5. Three grid refinement levels are considered in this work, the details of which are given in table 4. Results for the grids $9216\times 384$ (dashed) and $12\,288\times 512$ (solid) are in good agreement, as the coarse grid $6144\times 256$ (dotted) differs slightly. The DNS results from the grid $9216\times 384$ have been used for all succeeding analyses.

5.1.2. Behaviour of power spectra at higher amplitudes

Figure 6. Power spectra evolution in the streamwise direction is shown for DNS runs on the grid $9216\times 384$ with forcing type BL5 (described in table 3). The PSDs are shown for streamwise locations in the range x = 35 () to 95 cm (), using a sequential colour map (arrows designate increasing $x$ ). The forcing amplitude $A_0$ increases row-wise. A temporary attenuation effect in the PSD growth, observed at higher fluctuation levels, is highlighted with dashed boxes.

Figure 6 shows a more detailed streamwise evolution of the DNS power spectra. The PSDs are shown at uniformly spaced locations in the range $x=35-95\,\textrm{cm}$ at intervals of $5\,\textrm{cm}$ . The figure has been broken down into three rows to avoid plot overcrowding. The PSD evolution of the linear amplitude case is omitted since it is qualitatively similar to the weakly nonlinear amplitude. The spatial PSD evolution can be divided into two distinct stages. In the first stage, immediately downstream of the forcing region between $35$ and $50\,\textrm{cm}$ , the spectral content at all frequencies increases monotonically without any frequency shift. In the next stage, after $50\,\textrm{cm}$ , the second modes start to amplify with a steady frequency shift. The forcing region having ended at $x=30.9\,\textrm{cm}$ , the unstable modes have not modulated with the BL height in the first phase and cannot sustain standing wave resonance yet. Sivasubramanian & Fasel (2014) points out that just downstream of a localized forcing, the wave packets are predominantly three-dimensional, and the axisymmetric second modes do not start to destabilize readily. The steady frequency shift and amplification observed in the second stage is a characteristic of second-mode evolution and is known to occur due to the increasing BL height.

However, the PSD evolution deviates from this trend at higher forcing amplitudes as a temporary nonlinear attenuation is observed, which has been marked in figure 6 with dashed rectangles. The first occurrence of this for $A_0=1\,\textrm{m}\,\textrm{s}^-{^1}$ is at $x=55\,\textrm{cm}$ after the forcing mechanism stabilizes and continues to occur periodically downstream. For the nonlinear amplitude $A_0=0.1\,\textrm{m}\,\textrm{s}^-{^1}$ , this behaviour is only seen temporarily at $x\approx 0.6\,\textrm{m}$ . This effect was observed in the experiments conducted by Casper et al. (2016) at Sandia’s HWT-8 wind tunnel facility on a $7^\circ$ cone. Schlieren imaging showed signs of intermittent turbulence in the BL, which were reasoned to be the culprit for the abrupt attenuation. Chynoweth et al. (2019) also observed a similar effect on a Mach 6 flared cone as second-mode pressure fluctuations were seen to drop with simultaneous spectrum broadening. This was similarly concluded to be a sign of intermittency and ensuing breakdown. It is to be noted that these attenuations of the power spectrum are also observed in the axisymmetric DNS of the current work, where three-dimensional turbulence/intermittency effects are not accounted for, but nonlinear mechanisms are captured. Succeeding the PSD dip, the power spectra amplification of the primary mode starts to saturate, while the higher-frequency modes are seen to amplify more aggressively. In § 6, the perturbation energy corollary is used to investigate this effect further. However, before moving to analyse the nonlinear behaviour, the linear behaviour is first validated against locally parallel LST predictions.

5.2. Linear and weakly nonlinear dynamics

5.2.1. DNS versus LST pressure spectra

Linear stability theory serves two purposes in this work. Firstly, the DNS results for the linear forcing cases are validated against LST results. Secondly, any deviation of the DNS results from LST predictions helps identify the onset of nonlinear effects. Instead of comparing the $N$ -factors obtained from LST, they are used to predict the pressure spectra for direct comparison with DNS, as shown in figure 7. Using the phasor representation of a harmonic wave $e^{i (\alpha x - \omega t )}$ , the $N$ -factor for the spatial stability theory can be given as

(5.1) \begin{equation} N=\ln \left (\frac {\hat {p}(\omega ,x)}{\hat {p}_0(\omega ,x_0)}\right )=\int _{x_0}^x -\alpha _i (\omega ,x) \,\textrm{d}x, \end{equation}

where $\hat {p}_0(\omega ,x_0)$ denotes the Fourier spectrum of the wall pressure fluctuations at an arbitrarily chosen location $x_0$ . For the results shown, the DNS spectrum at $x_0=0.6\,\textrm{m}$ was chosen, which is a location where the applied forcing is deemed to be settled.

Figure 7. Comparison of DNS (solid lines) and LST predictions (lines with markers) of the power spectra of wall pressure fluctuations for frequencies ranging from $110-130\,\textrm{kHz}$ at intervals of $5\,\textrm{kHz}$ . This frequency span covers the typical range of second-mode frequencies for the base flow considered in this work.

In the comparisons shown in figure 7, only the linear ( $A_0=10^{-3}\,\textrm{m}\,\textrm{s}^-{^1}$ ) and the weakly nonlinear ( $A_0=10^{-2}\,\textrm{m}\, \textrm{s}^-{^1}$ ) cases are shown since the highly nonlinear cases deviate significantly from the linear theory predictions. The figure shows the comparison for forcing types BL5, L1 and C1. For type BL5, the linear amplitude forcing matched well with LST predictions, while the weakly nonlinear amplitude shows deviations for frequencies $110\,\textrm{kHz}$ and $115\,\textrm{kHz}$ (highlighted in figure 7 biii). On the other hand, forcing type L1 showed no such deviation. As given in table 3, the forcing was applied for one cycle of $120\,\textrm{kHz}$ frequency for type L1, while it was done for a range of frequencies over five cycles for type BL5. This shows that the onset of nonlinear effects is dependent not only on the forcing amplitude but also on the frequency spectrum and the time period over which the forcing is applied. The blowing/suction forcing (case C1) followed a similar trend as type L1. This suggests that the transition induced through forcing types L1 and C1 are predominantly in the linear regime. Additionally, figure 7 shows that the spectral energy for the blowing/suction pulse is much lower than that for forcing type L1, even though the forcing time periods, frequency content and amplitudes are identical for both. This higher receptivity of the harmonic forcing over the blowing/suction type is due to the richer spatial features in the former compared with the smoother cosine-cubed function used in the latter (3.1). Achieving a similar fluctuation level through the blowing/suction pulse requires the imposition of exceedingly high amplitudes, leading to significant mean flow distortions.

5.2.2. Phase speeds and growth rates of unstable modes: wave steepening

The pressure spectra deviation observed at modes $110\,\textrm{kHz}$ in the previous section is now investigated further. Comparison is made with the behaviour at $120\,\textrm{kHz}$ , which follows the linear theory predictions, as shown in figure 7. The DNS and LST predictions of the growth rates and phase speeds of the two frequencies are evaluated to aid this analysis. Using the phasor representation of wall pressure fluctuations and neglecting the base pressure gradient (Balakumar & Kegerise 2015), it can be seen that

(5.2) \begin{equation} \alpha (\omega ,x)=\alpha _r+j\alpha _i = \frac {1}{j\hat {p}(\omega ,x)}\frac {\partial \hat {p}(\omega ,x)}{\partial x}, \quad c_r(\omega ,x)=\textrm {Re} \left \{\frac {\omega }{\alpha (\omega ,x)}\right \}, \end{equation}

where, $\alpha$ denotes the complex wavenumber in the streamwise direction and $c_r$ denotes the streamwise phase speed. The Fourier spectra, available at a set of discrete frequencies, are interpolated using a piecewise cubic Hermite interpolating polynomial to obtain the pressure spectra at any given frequency. The spatial derivative of the spectra, $\partial \hat {p}/\partial x$ , is then calculated using the same C6 scheme as in DNS.

Figure 8. Phase speeds (left-hand vertical axis, black) and growth rates (right-hand vertical axis, red) for $f=110$ (a) and $120\,\textrm{kHz}$ (b) are shown. Two forcing types, L1 and BL5 (see table 3) are compared, with initial amplitudes $A_0=10^{-3}\,\textrm{m}\,\textrm{s}^-{^1}$ (dashed) and $A_0=10^{-2}\,\textrm{m} \, \textrm{s}^-{^1}$ (solid). The linear theory (LST) predictions are shown as circles. Mode F (S) denotes the fast (slow) acoustic modes, which propagate supersonically (subsonically) relative to the BL edge velocity $U_e$ .

The comparison of growth rates and phase speeds for forcing cases BL5 and L1 have been shown in figure 8. Non-dimensional phase speeds of the fast (Mode F) and the slow acoustic modes (Mode S) are also shown, which are $1+1/M_e$ and $1-1/M_e$ , respectively. The second-mode instabilities are seen to be propagating subsonically relative to the BL edge convective velocity $U_e$ at all times. The DNS growth rates and phase speeds for the linear amplitude agree with the LST predictions for both forcing types. The pressure spectra at 110 kHz are susceptible to low signal-to-noise ratios (SNR), as the spectrum at 110 kHz is typically six orders of magnitude lower than the peak amplitude in the upstream region (until $x\lt 0.7\,\textrm{m}$ ). This is why the phase speed and growth rate extraction from DNS results via a derivative calculation (5.2) are susceptible to the irregularities observed in figure 8. The SNR at 120 kHz is relatively higher, resulting in a smoother prediction. Farther downstream, as the spectra content at the respective frequencies increases, the behaviour of the power spectra derivatives is numerically more stable.

The phase speed for the weakly nonlinear amplitude imposed with forcing type BL5 peaks at the same downstream location ( $x \approx 0.8-0.9\,\textrm{m}$ ), where the pressure amplitude dropped (denoted as nonlinear deviation in figure 7). This indicates the pressure drop to be a consequence of the sudden decrease in growth rate as the phase speed increases (figure 8). An increase/decrease in the phase speed of a mode is an indication of wave steepening. As the phase speed of a mode increases or decreases, it leads to dispersion of the instability wave packet and creates wave steepening effects. This effect is observed to be significant for the more nonlinear forcing type BL5 compared with type L1 with the same amplitude. Thus, the initial forcing characteristics play a major role in initiating the steepening effect. This steepening phenomenon is a precursor to nonlinear breakdown and is an indication of the onset of nonlinearities in the transition process.

5.2.3. Streamwise evolution of growth rates: signs of nonlinearities

Figure 9. Spatial evolution of the modal growth rates $-\alpha _i$ of the disturbance predicted using DNS (contour) and LST (overlaying white isolines). The four forcing types (see table 3) are shown in each row with the same amplitude level $A_0=10^{-2}\,\textrm{m}\, \textrm{s}^-{^1}$ (weakly nonlinear). The development of nonlinear effects can be seen at the localized hot spots as the growth rate spectrum broadens (highlighted with red arrow), which is visible for all the forcing types except the blowing/suction pulse (C1).

The growth rate, $-\alpha _i (\omega ,x)$ , is now presented as a neutral stability curve in figure 9, with the LST predictions overlaid on it. These plots show the streamwise evolution of growth rates at different frequencies and help identify the nonlinear deviations over a broader spectral range. As in figure 8, spurious growth rates appear at low frequencies until $x\lt 0.7\,\textrm{m}$ due to the low SNR of the pressure spectra. A $\tanh$ -smoothening has been enforced along $x$ to remove the spurious artefacts while taking care that the spectral broadening later downstream is not affected. The smoothening is only done for frequencies below the neutral curve.

Local hot spots of the growth rate can be identified (lower right-hand side), after which the disturbance growth is observed to saturate. All four forcing types for the weakly nonlinear forcing ( $A_0=10^{-2}\,\textrm{m}\, \textrm{s}^-{^1}$ ) are now presented together. The blowing/suction pulse (case C1) shows no evidence of nonlinear action, as the DNS and LST neutral curves coincide with each other. Case L1 $(n_{c}=1)$ showed signs of nonlinearity at $x\approx 0.9\,\textrm{m}$ , whereas case L5 $(n_{c}=5)$ behaves similarly to case BL5 $(n_{c}=5)$ , with the spectral broadening occurring earlier at $x\approx 0.8\,\textrm{m}$ . A broad-band forcing frequency content promotes mode-to-mode interactions, while a higher forcing time increases the second-mode amplification rate, triggering nonlinear mechanisms. Due to this nonlinear action, frequencies below the LST lower neutral branch are observed to amplify farther downstream. The presence of such lower frequency content has been observed by Chynoweth et al. (2019) on a flared cone and compared with nonlinear PSE calculations using finite bandwidth initial disturbances. The DNS results in the current work employ a similar broad-band forcing, demonstrating that apart from the spectral bandwidth of the initial disturbance, the forcing time also causes an aggravated broadening effect.

The preceding sections show that the DNS results agree well with the LST predictions for lower amplitude forcing. The DNS captures the full nonlinear dynamics of the BL and is used to identify deviations from the linear theory predictions at the weakly nonlinear forcing. At such forcing levels, the linear regime persists much longer than the shorter nonlinear regime. The streamwise evolution of the modal growth rates is used to identify the linear-to-nonlinear regime transition, as indications of wave steepening and spectral broadening are captured. At higher amplitudes, the growth of the power spectra of wall pressure fluctuations is not monotonic. It is observed to attenuate periodically in between successive growth phases. This is believed to be an indication of nonlinear effects and is analysed in the following section using the disturbance energy corollary.

6. Disturbance energy budget analysis

The perturbation energy corollary is a conservation equation describing the spatiotemporal evolution of fluctuations in a flow. The DNS results have been used to close the energy budgets, which involves equating the energy sources on the right-hand side of the conservation equation (4.18) with the rate of change of energy on the left-hand side. The closure is attempted for the second-order energy equation only. Additionally, as shown in appendix D, the possibility of spurious growth in the perturbation energy definition has also been inspected.

6.1. Closure of the energy budget

For closing the perturbation energy budget using DNS data, (4.18) is integrated over an axisymmetric observation window to obtain the energy evolution over time,

(6.1) \begin{equation} \frac {\partial }{\partial t}\left ( \int _{\Omega _V} E_2 \,d\Omega _V \right ) + \int _{\Omega _V} \frac {\partial I_{2j}}{\partial x_j} \,d\Omega _V = \int _{\Omega _V} D_2 \,d\Omega _V, \quad d\Omega _V=\,\textrm{d}x\,\textrm{d}y\,Y\textrm{d}\theta , \end{equation}

where, $\Omega _V$ denotes the axisymmetric control volume, with $x$ and $y$ being the streamwise and wall-normal cone-bounded coordinates, respectively. Here $Y$ is the Cartesian coordinate (see figure 1) and $\theta$ denote the azimuthal angle. Supplemental information regarding the calculation of entropy gradients and viscous perturbations required in the above calculation is given in appendix C.

Table 5. Spans of the observation window $\Omega _V$ in the streamwise $x$ and wall-normal direction $y$ , which is used for evaluation of the energy budget integral in (6.1).

Figure 10. Closure of the disturbance energy budget (6.1) is demonstrated by equating the rate of change of energy ${\partial E_2}/{\partial t}$ (solid line) with the source term $D_2$ (line with markers). The control volumes (see table 5) used for evaluating the energy budgets have been indicated as red boxes in the numerical schlieren images, which also show the second-mode waves. The row below each schlieren image shows the time evolution of the disturbance energy for increasing forcing (BL5) amplitudes.

6.1.1. Control volume for energy budget calculations

The integral limits $\Omega _{V,x}$ and $\Omega _{V,y}$ in the streamwise and wall-normal directions, respectively, define the control volume or observation window $\Omega _V$ . Keeping the wall-normal span constant at $\triangle y=7\,\textrm{mm}$ , three different streamwise spans have been used (given in table 5). All domains start at $x=0.33\,\textrm{m}$ (downstream of the forcing region) and at $y=0$ (cone surface). The three $\Omega _V$ roughly capture the locations where disturbances of different amplitude levels start to deviate from linear behaviour. The axisymmetricity of the control volume is accounted for analytically using a constant azimuthal span. The time span for the integral has been adequately chosen for each control volume $\Omega _V$ such that the second-mode instabilities do not cross the boundaries of the $\Omega _V$ . This eliminates the energy flux term $I_{2j}$ at the side and the top boundaries by Gauss’ divergence theorem. The energy flux at the cone surface evaluates to zero due to the no-slip boundary condition for the base flow and homogeneous boundary conditions for the velocity fluctuations. This allows the disturbance power $\frac {\partial E_2}{\partial t}$ to be equated with the source terms $D_2$ directly. Note that $E_2$ is energy density; hence, the integral evaluated in (6.1) has the SI units of power and has not been specified in the following figures. For all analyses in this section, DNS results of the broad-band Legendre forcing (type BL5) have been used.

6.1.2. Behaviour of the energy budget

Figure 10 compares the time evolution of the rate of change of disturbance energy (solid lines) and the energy sources (markers). In each of the panels of figure 10, the DNS schlieren is shown in the top row, designating the second-mode position at the end of the temporal span used for the energy budget calculations shown in the following row. As can be observed in the first observation window $\Omega _{V_1}$ , the source term for the strongly nonlinear forcing starts to deviate, whereas the energy budgets for the three lower amplitudes close well (figure 10 a). For the nonlinear amplitude, the budget starts to diverge similarly in $\Omega _{V_2}$ (shown in figure 10 b). The larger control volume $(\Omega _{V_3})$ shows the disturbance energy for the strongly nonlinear forcing ( $A_0=1\,\textrm{m}\, \textrm{s}^-{^1}$ ) to show temporary attenuations periodically between successive growth phases (figure 10 c), while the energy growth persists following a period of attenuation for the nonlinear forcing ( $A_0=10^{-1}\,\textrm{m}\, \textrm{s}^-{^1}$ ).

For the lower amplitude forcings, $A_0=10^{-3}$ and $10^{-2}\,\textrm{m}\, \textrm{s}^-{^1}$ , the disturbance power increases exponentially, closely matching the source. In contrast, the growth is not monotonically increasing at higher amplitudes, as the disturbance energy is observed to saturate. This deviation from the monotonic growth overlaps with the periodic dips seen in the pressure power spectra at higher amplitudes (figure 6). This behaviour of the disturbance energy occurs periodically for the strongly nonlinear case while only occurring temporarily for the nonlinear case. Power spectra in figure 6 confirm the same behaviour, suggesting a common mechanism responsible for this effect. The axisymmetric nature of the simulations dismisses turbulence as a dissipation mechanism. This indicates the presence of a saturation mechanism that causes the disturbance energy growth to show repeating cycles of amplification and attenuation at higher forcing amplitudes.

6.2. Term-by-term investigation of disturbance energy equation

To pinpoint the reason for the simultaneous dip in the energy budget and the power spectra, a term-by-term breakdown of the energy budget is done. The energy $(E_2)$ and the source $(D_2)$ term variations along the wall-normal direction are shown in figure 11, along with the pressure field demonstrating the instability position at a specific time. Furthermore, figure 12 shows the time variation of the $E_2$ (figure 12 a) and $D_2$ (figure 12 b) terms, evaluated in the larger control volume $\Omega _{V_3}$ . The description of the terms and the corresponding colour codes used in the figures are tabulated in table 6 and table 7. These two figures are discussed in the following two sections.

Figure 11. Behaviour of individual terms in the disturbance energy $E_2$ and source terms $D_2$ along the wall-normal direction $y$ (colour codes given in table 6 and 7). Initial forcing (BL5) amplitude levels are increasing row-wise. Pressure contours expose a double-deck structure of the instabilities, with the acoustic part trapped below the relative sonic line and the thermal part localized at the critical layer.

Table 6. Colour codes of the disturbance energy $(E_2)$ terms shown in figures 11 and 12.

Table 7. Colour codes of the energy source terms $(D_2)$ shown in figures 11 and 12. The $[+]$ signs denote the production terms that add to the disturbance energy, while $[-]$ signs denote the dissipation terms.

6.2.1. Analysis of second-order disturbance energy, $E_2$

The disturbance energy comprises the thermal PE (term $i$ ), mean flow energy (term $ii$ ), KE (term $iii$ ) and acoustic PE due to compression/rarefaction (term $iv$ ) (Chu 1965). Thermal energy is observed to dominate the energy content in all cases. Terms $iii$ and $iv$ , which constitute the total acoustic energy in the absence of mean flow, increase monotonically in time as the second modes propagate downstream. The behaviour of term $ii$ , which signifies the acoustic energy carried by the mean flow, is irregular.

Five lines that characterize the BL physics (Knisely & Zhong 2019) are marked in figure 11, which are:

1. (i) the BL edge $ (\delta _{BL} )$ ;

2. (ii) the generalized inflection line– where $\frac {\partial }{\partial y} \left(\rho \frac {\partial u_0}{\partial y} \right)=0$ ;

3. (iii) the critical layer– where relative Mach number $\overline {M}(y)=\frac {u_0(y)-c_{r}}{a(y)}=0$ ;

4. (iv) the relative sonic line– where relative Mach number $\overline {M}(y)=-1$ ;

5. (v) the local sonic line– where local Mach number $M(y)=\frac {u_0(y)}{a(y)}=1$ .

The critical layer and the relative sonic line are calculated using the second mode phase speed $(c_r)$ evaluated earlier in § 5.2.2.

The acoustic PE (term $iv$ ) dominates the disturbance energy content $E_2$ near the wall as the pressure fluctuations are maximum at the wall. In the neighbourhood of the relative sonic line, the KE and mean flow energy dominate the disturbance energy content. The relative sonic line marks the region below which the instability propagates supersonically relative to the mean flow (Fedorov 2011). The acoustic and thermal instabilities can only interact with each other above the relative sonic line when their respective phase speeds synchronize. The acoustic instabilities cannot interact with the base flow above the relative sonic line and consequently propagate as trapped waves beneath it. Farther above, the thermal energy dominates the acoustic energy, as the temperature and density fluctuations reach a maximum near the critical layer $(\approx 0.8\delta _{BL})$ . In contrast to all the other energies, the mean flow acoustic energy can be seen to change signs and reach a minimum at the inflection line (see figure 11). This negative value of the mean flow energy occurs due to its definition, as the density and velocity fluctuations are out of phase with each other at the generalized inflection line. Appendix B.1 explores the positive definiteness of this term in the present context.

Figure 12. Time evolution of individual terms in the disturbance energy $E_2$ (a) and source terms $D_2$ (b), for initial forcing amplitudes $A_0=10^{-2},\,10^{-1},\,1\,\textrm{m}\, \textrm{s}^-{^1}$ . The colour codes are given in tables 6 and 7. In panel (b), the dotted lines indicate the dissipation terms $(-)$ , while the solid lines denote production terms $(+)$ . The balance between these terms dictates the growth/attenuation of the instabilities.

6.2.2. Analysis of second-order disturbance sources, $D_2$

The interplay between the thermal and acoustic fluctuations is now analysed using the source terms (4.20). The dissipation terms in $D_2$ are term III due to heat diffusion and term V due to viscous action. The second modes are known to be inviscid-type instabilities (Mack 1969). Accordingly, term V is observed to not play a major role until farther downstream when the fluctuation levels become significant. On the other hand, thermal diffusion is active even at lower amplitudes. The strongly nonlinear forcing (figure 12 a iii, b iii) shows that the valleys in the total disturbance source $(D_2)$ correspond with the peaks in the viscous (V) and thermal diffusion (III) terms. This suggests that at sufficiently high fluctuation levels, the diffusive effects $(-)$ overpower the energy production $(+)$ periodically in between intervals of growth. These diffusion effects are seen to be significant at the wall and the critical layer (figure 11).

Figure 12 also demonstrates that the thermal diffusion (term III) increases in magnitude or becomes more negative in response to the growth of term IV. At lower (linear) amplitudes, term IV is always higher than term III, as the magnitude of both increases exponentially. Whereas, at higher (nonlinear) amplitudes, term III $(-)$ periodically surpasses term IV $(+)$ , balancing the total disturbance energy production.

Term IV represents the energy production occurring due to the entropy fluctuations advecting through a base entropy gradient. For a hypersonic BL with zero pressure gradient, the base entropy gradient is proportional to the base temperature gradient, i.e., $\partial s_0/\partial y \approx C_p(T_0)^{-1}\partial T_0/\partial y$ . Term IV is dominated primarily by the term $-p_0 s^{\prime} (R C_p )^{-1} v^{\prime} \partial s_0/\partial y\,(+)$ , with a relatively lower contribution from $-p^{\prime} v^{\prime} (C_p)^{-1} \partial s_0/\partial y\,(+)$ . The importance of term IV suggests that the interaction of the fluctuations with the base temperature gradient of the BL is one of the primary mechanisms of second-mode growth. On the other hand, thermal diffusion (term III) is dominated by $s^{\prime}(Pr)^{-1} \partial ^2 T^{\prime}/\partial y^2\,(-)$ near the BL edge, with a lower contribution from $R\mu _0 p^{\prime} (p_0Pr)^{-1} \partial ^2 T^{\prime}/\partial y^2\,(-)$ at the wall. This suggests that the heat generated by the advecting perturbation near the critical layer is diffused away from it. The direction of this diffusion is dictated by the fluctuating temperature field, as is evident from the presence of the term $\partial ^2 T^{\prime}/\partial y^2$ . This energy production near the BL edge can interact with the compression/rarefaction waves below, creating a feedback mechanism that promotes the growth of the acoustic-type second-mode waves.

This interplay is evident from the behaviour of terms I and II (figure 12), which are seen to balance each other irrespective of the amplitude levels. Term I facilitates the acoustic–thermal energy coupling and is dominated by the term $-p^{\prime}u_0 (C_p)^{-1} \partial s^{\prime}/\partial x\,(-)$ . This feeds the terms $-\rho^{\prime} u_0 v^{\prime} \partial u_0/\partial y\,(+)$ and $-(\rho^{\prime} u_0)^2 (\rho _0)^{-1} \partial u_0/\partial x\,(+)$ in term II, causing the disturbance KE (term $iii$ in $E_2$ ) to increase. This is verified from figure 11 as both terms $iii$ of $E_2$ and II of $D_2$ increase at the neighbourhood of the relative sonic line. The disturbance energy is further augmented by term VII, which is observed to only activate when the fluctuation levels are high. This is the only term exhibiting disturbance energy production below the relative sonic line and is dominated by the term $\rho _0 u_0 v^{\prime} \partial u^{\prime}/\partial y\,(+)$ , suggesting the presence of vorticity effects at high fluctuation levels. The viscous diffusion at the wall (term V) is observed to grow in response to the evolution of term VII.

To summarise, the primary source for second-mode energy growth is identified to be the entropy/velocity fluctuations interacting with the base entropy (temperature) gradient. At higher amplitudes, the disturbance energy growth encounters periodic attenuation due to thermoviscous diffusion effects that counteract the production mechanisms. The balance between the terms V $(-)$ and VII $(+)$ near the wall and of terms III $(-)$ and IV $(+)$ near the critical layer is collectively responsible for this effect that leads to saturation of the disturbance energy growth. The periodic dips in the power spectra of the wall pressure fluctuations, discussed in § 5.1.2, can be attributed to this attenuation effect. Whereas at lower amplitudes, the contributions from the production terms IV and VII are always higher than the dissipation effects, facilitating the monotonic growth of the second-mode waves in the linear regime.

7. Non-dimensionalization and scaling analysis

In this section, the disturbance energy saturation is further examined by non-dimensionalizing the disturbance energy. This promotes a scaling analysis of the second mode growth in time, allowing a comparison with the effects seen in the modal growth rates of § 5.2.3.

7.1. Non-dimensional disturbance energy

So far, the disturbance energies and sources have been presented dimensionally. A suitable factor to obtain a non-dimensional measure of the perturbation energy is now introduced. The factor used is the root-mean-square (r.m.s.) of the acoustic energy introduced into the domain at the forcing region. Scaling the second-order disturbance energy by this factor quantifies the extent to which the input acoustic energy is magnified as the second mode amplifies. The instantaneous acoustic power $\mathcal {P}_a (t)$ introduced at the forcing region and the corresponding r.m.s. acoustic power $\mathcal {P}_{rms}$ are defined as

(7.1) \begin{equation} \mathcal {P}_a(t) = \int _{\mathcal {S}} p^{\prime}v^{\prime} \, dS, \quad \mathcal {P}_{rms} = \sqrt {\frac {1}{t_f}\int _0^{t_f} \mathcal {P}_a^2(t) \,\textrm{d}t}, \end{equation}

where $\mathcal {S}$ denotes the forcing region on the cone surface, and $t_f$ denotes the total time over which the forcing is applied. Here $\mathcal {P}_{rms}$ is expected to be $\approx k\rho _w c_w \mathcal {S}A_0^2$ , where subscript ‘ $w$ ’ indicates the wall values and $k$ is a constant of proportionality. The r.m.s. power defined above yields the r.m.s. acoustic energy fed into the flow at the disturbance region as $E_{rms}=\mathcal {P}_{rms}t_f$ . Note that $E_{rms}$ is heavily dependent on the forcing type and scales by the forcing amplitude squared.

7.2. Linear stability theory prediction of disturbance energy time evolution

The LST has been used to estimate the time evolution of disturbance energy in a fixed control volume. As detailed in § 5.2.1, the $N$ -factors from spatial LST can be used to predict the wall pressure spectra evolution along $x$ . From (5.1), it simply follows that PSD $(\omega ,x)=$ PSD $(\omega ,x_0) e^{2N(\omega ,x)}$ , where the power spectrum at an initial location $x_0$ is informed from the DNS data. Thus, LST can directly estimate the evolution of the acoustic part of the second-order energy $E_2$ , which is represented by $p^{\prime 2}/2\gamma p_0$ (term $iv$ ). Additionally, figure 12(a i, b i) shows that all the energy and source terms increase or decrease exponentially for the lower (linear) forcing amplitudes. This qualifies the assumption that in the linear regime, the disturbance energy also follows an exponential trend $(e^{2N})$ similar to the power spectra.

Owing to the second-mode waves having nearly uniform phase speeds across different frequencies (see figure 8), the group velocity of the wave packet can be assumed constant at $v_g=c_r\approx 0.92U_e$ . Starting from an arbitrarily chosen position $(x_0)$ of the wave packet at time $t_0$ , one can obtain positions of the wave packet at other times $t$ as $x=x_0+v_g(t-t_0)$ . Here, the initial chosen position $x_0$ is taken as the streamwise location of the maximum second-order energy $ (E_2 )$ at time $t_0$ . This approach allows one to conveniently use the spatial LST $N$ -factors instead of resorting to temporal LST.

Figure 13. Time evolution of the scaled disturbance energy for different forcing amplitudes $A_0$ . The solid lines show the scaled energy $E_2/E_{rms}$ , and the markers show the disturbance sources $D_2/E_{rms}$ . The dotted line denotes the disturbance energy growth predicted by LST. Increasing amplitudes are shown with an offset ( $n$ ) along the vertical axis. The values of $n$ for each line are shown on the lower left-hand side.

7.3. Scaling analysis results

Figure 13 compares the time evolution of the scaled disturbance energy for all four forcing amplitudes initiated by forcing type BL5. Non-dimensional time, $\tau ^*=t U_\infty /L$ , is used for the horizontal axis, where $U_\infty$ is the free stream velocity and $L$ is the length of the cone. For all succeeding analyses, $\Omega _{V_3}$ (see table 5) has been used as the fixed control volume over which the DNS disturbance energy data are numerically integrated using (4.16). This figure clearly identifies the disturbance energy saturation as deviations from the LST trend. The disturbance energy evolution can be broken down into two distinct phases, which coincide with the phases of power spectra evolution discussed earlier in 5.1.2. The monotonic amplification after the forcing region constitutes the first phase, and the modal amplification of the second mode away from the forcing region constitutes the second. During the first phase $(\tau ^*\lt 0.3)$ , the amplification rate of the disturbance energy is the same for all amplitudes and is steeper compared with the second phase. The power spectra evolution in figure 6(a) confirms this amplitude-independent growth. During the second phase, the amplification rate is lower and amplitude-dependent. During this phase, the second mode resonates below the relative sonic line, and the frequency starts to decrease with increasing BL height downstream. The linear and the weakly nonlinear forcing follow the LST trend closely, while at higher amplitudes, the energy is observed to plateau.

Figure 14. A qualitative analysis of the disturbance energy saturation phenomena: (a) branching of the scaled disturbance energy $(E_2/A_0^2)$ from linear behaviour at a higher forcing amplitude of $A_0=10^{-1}\,\textrm{m}\, \textrm{s}^-{^1}$ , starting at $t\gt 0.4\,\textrm{ms}$ ; (b) schlieren showing the position of the second mode at $t=0.4$ ms as the leading wavefront enters $x\sim 0.6\,\textrm{m}$ ; (c) modal growth rates extracted from DNS data (see § 5.2.3) shows nonlinear spectral broadening to occur $x\sim 0.6\,\textrm{m}$ onwards.

Using the methodology presented in § 5.2.3, DNS predictions of the modal growth rates $(-\alpha _i)$ offer an alternative approach for analysing the saturation phenomena. Figure 14 presents the results for the nonlinear forcing case ( $A_0=10^{-1}\,\textrm{m}\, \textrm{s}^-{^1}$ ). Figure 14(a) presents a simpler scaling of the disturbance energy using $A_0^2$ . This allows the time evolution to be presented on a linear scale, magnifying the branching from linear behaviour. Scaled disturbance energy shows the attenuation by thermoviscous mechanisms to trigger at $t=0.4\,\textrm{ms}$ (highlighted box in figure 14 a) when the leading wavefront of the instability enters $0.6\,\textrm{m}$ (figure 14 b). Figure 14(c) shows the DNS growth rate to deviate from the LST predictions $0.6\,\textrm{m}$ onwards, confirming the thermoviscous attenuation effects to occur simultaneously with the nonlinear spectral broadening. The hot spots in figure 14(c), which are also shown earlier in figure 9, serve as an indication of the transition from linear to nonlinear behaviour of the instabilities.

7.4. Nonlinearities at different noise types

Figure 15. Time evolution of the scaled disturbance energy $E_2/E_{rms}$ for different forcing types (detailed in table 3) with different forcing amplitudes $A_0$ .

The evolution of the scaled disturbance energy is now examined for the four forcing types BL5, L5, L1 and C1. The difference between the forcing types (given in table 3) is repeated here for the reader’s convenience. Forcing type BL5 imposes Legendre-space modulated noise for discrete frequencies in the range $60-200\,\textrm{kHz}$ . The forcing is imposed for a time period of fivecycles of the lowest frequency $60\,\textrm{kHz}$ . Forcing type L5 imposes a Legendre-modulated field at a single frequency of $120\,\textrm{kHz}$ for a period of five cycles. Type L1 imposes the same frequency as L5 but for a single cycle. Finally, the blowing/suction pulse (type C1) has the same temporal behaviour as type L1 but with a $\cos ^3$ -spatial modulation instead.

For identical amplitudes, the forcing types BL5 to L1 signify decreasing levels of intermodal interactions, while C1 is used to study the effect of spatial modulation. Figure 15 compares the temporal evolution of the scaled disturbance energy for all forcing types. The four amplitude levels, $A_0=10^{-3},\,10^{-2},\,10^{-1}$ and $1\,\textrm{m}\, \textrm{s}^-{^1}$ , are shown in separate figure panels. The LST predictions are informed by the DNS power spectra of the BL5 forcing. The effects of disturbance energy saturation are evident for type BL5 at $A_0=0.1\,\textrm{m}\, \textrm{s}^-{^1}$ , while type L1 shows saturation effects only at the highest forcing amplitude, $A_0=1\,\textrm{m}\, \textrm{s}^-{^1}$ . For lower amplitudes, both BL5 and L1 follow the energy growth rates predicted by LST. Type L5, where the forcing is applied for a longer time, follows the trend of type L1 until it shows an abrupt increase in disturbance energy. This abrupt increase is observed consistently at all the forcing amplitudes. Further investigation is required to uncover the mechanism causing this effect and is not pursued in this work.

The blowing/suction type (C1) shows a very different trend than the Legendre-space modulated forcing. For all amplitudes, the disturbance energy for C1 forcing has much lower values compared with the energy introduced at the forcing region $(E_{rms})$ . A major portion of the introduced energy for C1 forcing is convected upwards and exits the BL through the top wall of the computational domain, leading to its linear behaviour at all amplitudes. This highlights the importance of the spatial characteristics of the forcing, which are often overlooked compared with the temporal characteristics. A Legendre space-based forcing triggers richer spatial features in the hypersonic BL that are capable of activating nonlinear mechanisms later downstream, which is not observed with the blowing/suction type forcing.

8. Conclusion

Time-accurate, axisymmetric DNS is performed to investigate the transition to turbulence in a hypersonic BL over a sharp, circular cone. A $3^\circ$ half-angle cone is considered, which is exposed to Mach 6 flow at a $0^\circ$ angle of attack. The linear and nonlinear phases of transition are, respectively, analysed using LST and a newly derived second-order disturbance energy equation that highlights the physical mechanisms causing flux-based transport, production and destruction of fluctuation energy in the BL. A novel noise generation methodology to trigger transition in the DNS runs has been adopted, demonstrating good grid-independent properties. The tunable parameters in the noise generation procedure are varied to perform a parametric study over a range of initial amplitude levels and frequency content. The amplitude levels span three decades, yielding gradually increasing nonlinear behaviour of the second-mode waves, which are the primary instabilities involved in hypersonic transition. The DNS results at the lower (linear) amplitude levels are validated against LST predictions. Phase speeds and growth rates of the second modes are calculated, demonstrating wave steepening effects at the intermediate (weakly nonlinear) amplitude forcing. For the higher (nonlinear) amplitude forcing, power spectra of the wall pressure fluctuations show a temporary attenuation effect, occurring periodically in between phases of steady growth. The derived energy conservation equation has been used to study this effect.

The disturbance energy conservation equation traces the source terms responsible for the growth/attenuation of second modes in the absence of external sources. The DNS results are used to show the closure of the energy budgets for different forcing amplitudes. The derived equation demonstrates the interaction between the acoustic and thermal part of the instabilities, revealing a two-deck structure in the BL. The acoustic PE arising from compression/rarefaction effects is concentrated below the relative sonic line, while the thermal part of the disturbance energy is dominant at the critical layer near the BL edge. The acoustic-type second modes, trapped near the wall, are seen to be destabilized through a feedback mechanism with the thermal fluctuations above. The primary reason for this destabilization is identified as the interaction of the advecting entropy fluctuations with the base entropy gradient in the BL. Without a pressure gradient, the base entropy gradient is equivalent to the base temperature gradient. The entropy fluctuations at the critical layer lead to disturbance energy production, which diffuses and feeds the acoustic potential and kinetic energies below the relative sonic line. At higher (nonlinear) amplitudes, the thermoviscous dissipation effects amplify in response to the energy production mechanisms and periodically overpower them. This is in contrast to the lower (linear) amplitude levels when the energy production always dominates the dissipation mechanisms, causing a monotonic growth of the second modes.

As the second modes amplify to sufficiently high fluctuation levels, thermoviscous effects appear to be significant at the critical layer and near the wall, leading to the above-mentioned temporary attenuation in the wall pressure spectra. This phenomenon was observed in the experimental investigations of Casper et al. (2016) and Chynoweth et al. (2019) and was argued to be an indication of intermittency and a precursor to turbulence breakdown. In the current work, this effect is explored further by a scaling analysis of the disturbance energy. For higher amplitude forcing, the scaled disturbance energy is observed to plateau and deviate from the linear theory predictions, whereas the lower amplitude forcing follows the LST-predicted trend. The scaling analysis exposes the energy level at which nonlinear mechanisms start dominating the instability behaviour, providing a threshold around which the disturbance energy saturates or shows non-monotonic growth.

Thus, the present work shows that the derived disturbance energy equation can be used to describe the spatiotemporal evolution of an instability wave and to identify the mechanisms causing the instability growth. Moreover, the equation also demonstrates the intricate coupling between the different energy terms, which makes isolating the effect of the individual source terms challenging. We note that the expression of the proposed perturbation energy $E_2$ matches Myers’ energy norm. However, the associated fluxes $I_2$ and source terms $D_2$ derived via linear combination of the governing equations for the thermofluid-dynamic fluctuations exhibit a different functional form than Myers’, that is more amenable to a mechanistic interpretation of second-mode dynamics.

Three-dimensional mechanisms in the transition process are not explored owing to the axisymmetric nature of the DNS in this study. In future works, the disturbance energy equation derived here can be applied with no change to analyse three-dimensional breakdown mechanisms in a hypersonic BL, where the interaction between two-dimensional primary instabilities and secondary or oblique modes can be explored.