1. Introduction
Predicting the stability of hypervelocity boundary layers is a crucial aspect of hypersonic vehicle design. Transition of the boundary layer from laminar to turbulent can increase surface heating loads by a factor of 4–10 (Lau Reference Lau2008; Lin Reference Lin2008; Hollis Reference Hollis2012), thereby increasing the cooling requirements and the weight of the thermal protection system. The location of the transition region is therefore a key design consideration, yet it remains extremely challenging to predict reliably. This uncertainty necessarily leads to a conservative aircraft design, which in turn increases the vehicle weight and power requirements. As a result, any gains in the understanding of boundary layer instability or advances in transition prediction translate into substantial savings in weight and cost.
Because of the high cost of flight tests and the difficulty of making instability measurements during actual flight, most experimental investigations of hypersonic boundary layer stability have been conducted in ground-based facilities. Reviews of the experimental data from such facilities are available from Morkovin (Reference Morkovin1969) and Schneider (Reference Schneider2001, Reference Schneider2004). Ground-based facilities cannot simultaneously replicate all of the flow conditions from real flight. Most hypersonic wind tunnels operate at low stagnation enthalpies in order to achieve long run times, reduce the difficulty of making measurements, or reach low levels of acoustic noise in the free stream. Other facilities, such as shock tunnels and expansion tubes, have been designed to match the stagnation enthalpy of real flight conditions, but this comes at the expense of short test times and high levels of free-stream noise. Past experimental campaigns (Adam & Hornung Reference Adam and Hornung1997; Germain & Hornung Reference Germain and Hornung1997; Rasheed et al. Reference Rasheed, Hornung, Fedorov and Malmuth2002; Fujii & Hornung Reference Fujii and Hornung2003) in such high-enthalpy facilities have focused mainly on non-equilibrium flow effects, both thermal and chemical. However, the temperature distribution within the boundary layer of a high-enthalpy flow is often very different from that of a low-enthalpy flow: at high stagnation enthalpy, the wall temperature can be a small fraction of the free-stream temperature, while for low-enthalpy conditions this is usually not possible. This can cause high- and low-enthalpy flows to have significantly different stability characteristics, even in the absence of thermochemical reactions. In comparing test data from low- and high-enthalpy test facilities, understanding these differences is important.
It is well known that the linear stability of compressible boundary layers is very sensitive to the temperature of the wall. This conclusion was first reached by Lees & Lin (Reference Lees and Lin1946), who generalized Rayleigh’s inflection point theorem (Rayleigh Reference Rayleigh1880) to compressible flows and demonstrated that, in the inviscid limit, an extremum of mean angular momentum
$(\bar{{\it\rho}}\bar{U}^{\prime })^{\prime }=0$
is a necessary condition for instability. Lees & Lin observed that sufficient wall cooling can eliminate this generalized inflection point, thereby stabilizing the boundary layer to inviscid disturbances. Subsequent work by Lees (Reference Lees1947) included the effects of viscosity, which are important because Tollmien–Schlichting waves can be more unstable at finite Reynolds number than in the inviscid limit (Lin Reference Lin1944), though this is usually only true for low Mach numbers (Masad, Nayfeh & Al-Maaitah Reference Masad, Nayfeh and Al-Maaitah1992). Lees showed that, in contrast to subsonic flows, supersonic flows can be completely stabilized at all Reynolds numbers if the level of wall cooling is great enough. This point was further investigated by Van Driest (Reference Van Driest1952), who provided calculations of the critical Reynolds number as a function of wall temperature and Mach number and concluded that the first mode can be stabilized for all Reynolds numbers as long as the Mach number falls in the range
$1<Ma<9$
. The previously mentioned studies of Lees and Van Driest considered only two-dimensional (2D) disturbances; however, oblique Tollmien–Schlichting waves are considerably more unstable than 2D waves (Mack Reference Mack1969). Nevertheless, it has been shown that the three-dimensional (3D) first-mode waves can still be stabilized with sufficient levels of wall cooling (Masad et al.
Reference Masad, Nayfeh and Al-Maaitah1992).
The investigations cited above were later found to be incomplete, since at high Mach numbers additional instability modes collectively known as the ‘Mack modes’ exist and are highly unstable. The calculations of Mack (Reference Mack1969, Reference Mack1975, Reference Mack1984, Reference Mack1987) revealed that these modes are acoustic waves reflecting between the wall and the relative sonic line. Although the Tollmien–Schlichting waves mentioned above can indeed be stabilized by wall cooling, the acoustic modes are destabilized. The most unstable acoustic mode, the ‘second mode’, has been widely studied and numerous factors that influence its stability have been investigated. These factors include wall cooling (Mack Reference Mack1975, Reference Mack1987; Malik Reference Malik1989a ; Masad et al. Reference Masad, Nayfeh and Al-Maaitah1992; Mack Reference Mack, Ashpis, Gatski and Hirsch1993), pressure gradient (Malik Reference Malik1989a ), chemical reactions (Malik & Anderson Reference Malik and Anderson1991; Stuckert & Reed Reference Stuckert and Reed1994; Chang, Vinh & Malik Reference Chang, Vinh and Malik1997; Hudson, Chokani & Candler Reference Hudson, Chokani and Candler1997; Johnson, Seipp & Candler Reference Johnson, Seipp and Candler1998; Malik Reference Malik2003), thermal non-equilibrium (Bertolotti Reference Bertolotti1998), transverse curvature (Malik & Spall Reference Malik and Spall1991), porous walls (Fedorov et al. Reference Fedorov, Malmuth, Rasheed and Hornung2001, Reference Fedorov, Shiplyuk, Maslov, Burov and Malmuth2003; Bres et al. Reference Bres, Inkman, Colonius and Fedorov2013) and many others.
Although many of the studies mentioned above have investigated the effects of wall cooling on boundary layer stability, most of them have considered only relatively modest levels of cooling. Many papers that analyse high-enthalpy flows (Malik & Anderson Reference Malik and Anderson1991; Stuckert & Reed Reference Stuckert and Reed1994; Hudson et al.
Reference Hudson, Chokani and Candler1997; Johnson, Alba & Candler Reference Johnson, Alba and Candler2008) do not consider wall-to-edge temperature ratios
$T_{w}/T_{e}<1$
, and, although a few papers do include such conditions (Chang et al.
Reference Chang, Vinh and Malik1997; Johnson et al.
Reference Johnson, Seipp and Candler1998; Malik Reference Malik2003), they focus on the effects of chemical reactions rather than on the extreme wall cooling. Presumably, cases with
$T_{w}/T_{e}<1$
are not often considered because most experimental facilities operate at low stagnation temperatures, hence the free-stream temperature is small (often less than 100 K), and cooling the wall below the free-stream temperature is not practicable. In contrast, high-enthalpy test facilities can reach free-stream temperatures of 1000–2500 K, and since the wall temperature remains ambient during the short test duration, wall-to-edge temperature ratios of
$T_{w}/T_{e}=0.1{-}0.3$
are typical. Such high levels of wall cooling are also relevant to some real flight scenarios; for example, Malik (Reference Malik2003) analysed the re-entry flight experiment reported by Sherman & Nakamura (Reference Sherman and Nakamura1970), which involved a 22° half-angle blunt cone at Mach 6, and Malik computed a post-shock temperature of approximately 3300 K as compared to a wall temperature of approximately 730 K, leading to
$T_{w}/T_{e}=0.22$
.
As will be seen, very high levels of wall cooling lead to unique features of the disturbance spectrum that are not present for lower levels of wall cooling. In particular, unstable supersonic modes are found that cause second-mode instability over a wider frequency band than for subsonic modes. Unstable supersonic modes have been noted in a few other studies of flat plates, wedges and cones (Mack Reference Mack1987; Chang, Malik & Hussaini Reference Chang, Malik and Hussaini1990; Chang et al. Reference Chang, Vinh and Malik1997), and they have also been encountered and found to be important for flows involving gas injection (Fedorov, Soudakov & Leyva Reference Fedorov, Soudakov and Leyva2014) and porous surfaces (Fedorov et al. Reference Fedorov, Bres, Inkman and Colonius2011; Bres et al. Reference Bres, Inkman, Colonius and Fedorov2013). However, for highly cooled boundary layers, there remains some uncertainty regarding the importance of these modes and the range of conditions for which they exist. In this paper, spatial, locally parallel stability calculations are used to investigate these effects. Air is chosen as the test gas, and the geometry is a flat plate with a sharp leading edge. Flow conditions are modelled after those found in a high-enthalpy shock tunnel. Chemical reactions are excluded in the modelling, which limits the validity of the results to flows having peak temperatures below approximately 2600 K. However, this restriction is not too severe; for instance, the model is still applicable to most experimental conditions in the T5 shock tunnel facility at Caltech (Germain & Hornung Reference Germain and Hornung1997; Parziale, Shepherd & Hornung Reference Parziale, Shepherd and Hornung2013).
2. Numerical methods
2.1. Governing equations
Although this study focuses on flow conditions for which chemical reactions are negligible, the gas temperature is sufficiently high that the exchange of energy between translational and vibrational molecular motion may not be negligible, and the possibility of thermal non-equilibrium must be taken into account. For a chemically inert gas mixture, the transport of mass, momentum, vibrational energy and total energy can be described by the following equations, which incorporate the two-temperature model of Park (Reference Park1990):
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\rho}^{\ast }}{\partial t^{\ast }}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\rho}^{\ast }\boldsymbol{u}^{\ast })=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}^{\ast }\frac{\text{D}\boldsymbol{u}^{\ast }}{\text{D}t^{\ast }}+\boldsymbol{{\rm\nabla}}p^{\ast }=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\tau}^{\ast }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}^{\ast }\frac{\text{D}e_{v}^{\ast }}{\text{D}t^{\ast }}=-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{q}_{v}^{\ast }+Q^{\ast }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}^{\ast }\frac{\text{D}h_{tr}^{\ast }}{\text{D}t^{\ast }}+{\it\rho}^{\ast }\frac{\text{D}e_{v}^{\ast }}{\text{D}t^{\ast }}-\frac{\text{D}p^{\ast }}{\text{D}t^{\ast }}=-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\boldsymbol{q}_{tr}^{\ast }+\boldsymbol{q}_{v}^{\ast })+{\bf\tau}^{\ast }:\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\ast }. & \displaystyle\end{eqnarray}$$
${\bf\tau}^{\ast }$
is the viscous stress tensor. The velocity vector
$\boldsymbol{u}^{\ast }$
has elements
$u^{\ast }$
,
$v^{\ast }$
and
$w^{\ast }$
in the streamwise, wall-normal and spanwise directions, and
$p^{\ast }$
and
${\it\rho}^{\ast }$
are the pressure and mass density. The system is closed with the equation of state: 
$$\begin{eqnarray}p^{\ast }={\it\rho}^{\ast }\mathscr{R}^{\ast }T^{\ast },\end{eqnarray}$$
where
$\mathscr{R}^{\ast }$
is the gas constant of the mixture and
$T^{\ast }$
is the translational temperature. Although the gas is chemically inert, thermal non-equilibrium is included in the model by dividing the flow enthalpy into a translational enthalpy
$h_{tr}^{\ast }$
and a vibrational energy
$e_{v}^{\ast }$
. Likewise, the heat flux vector is split into a vibrational part
$\boldsymbol{q}_{v}^{\ast }$
and a translational part
$\boldsymbol{q}_{tr}^{\ast }$
. The exchange of energy between translational and vibrational modes is taken into account by the source term
$Q^{\ast }$
.
It is convenient to define a vibrational temperature
$T_{v}^{\ast }$
such that, for the
$i$
th species,
$$\begin{eqnarray}\frac{\partial e_{v,i}^{\ast }}{\partial T_{v}^{\ast }}=c_{v,v,i}^{\ast }=\mathscr{R}_{i}^{\ast }\left(\frac{{\it\Theta}_{v,i}^{\ast }}{T_{v}^{\ast }}\right)^{2}\frac{\text{e}^{{\it\Theta}_{v,i}^{\ast }/T_{v}^{\ast }}}{(\text{e}^{{\it\Theta}_{v,i}^{\ast }/T_{v}^{\ast }}-1)^{2}},\end{eqnarray}$$
where
$c_{v,v,i}^{\ast }$
is the vibrational specific heat of the
$i$
th species,
$\mathscr{R}_{i}^{\ast }$
is its gas constant and
${\it\Theta}_{v,i}^{\ast }$
is its vibrational activation temperature, taken to be 3390 K for
$\text{N}_{2}$
and 2270 K for
$\text{O}_{2}$
(Vincenti & Kruger Reference Vincenti and Kruger1967). In all cases the gas mixture is assumed to be air, modelled as 78 %
$\text{N}_{2}$
and 22 %
$\text{O}_{2}$
by mole.
The model used in (2.3) assumes that the vibrational energy states of each molecule satisfy a Boltzmann distribution characterized by the vibrational temperature
$T_{v}^{\ast }$
. Equations (2.1) and (2.3) also assume that all diatomic species share the same vibrational temperature. This model is only accurate if vibration–vibration (V–V) energy exchange between unlike molecules (e.g.
$\text{N}_{2}$
and
$\text{O}_{2}$
) is much faster than vibration–translation (V–T) exchange. Experimental data compiled by Taylor & Bitterman (Reference Taylor and Bitterman1969) reveal that, in the temperature range of interest in this study, 1000–3000 K, V–V transfer between
$\text{N}_{2}$
and
$\text{O}_{2}$
is several orders of magnitude faster than V–T transfer for
$\text{N}_{2}$
, but is slightly slower than V–T transfer for
$\text{O}_{2}$
. This suggests that the use of a single vibrational temperature is a reasonable first approximation, but more accurate results would be obtained if separate vibrational temperatures were used and the V–V exchange process were explicitly included in the model, as was done by Bertolotti (Reference Bertolotti1998). For simplicity, however, a single vibrational temperature is assumed in the remainder of this paper.
With the assumption of a single vibrational temperature, the rate of V–T energy exchange is computed using the Landau–Teller model (Vincenti & Kruger Reference Vincenti and Kruger1967; Park Reference Park1990):
$$\begin{eqnarray}Q^{\ast }={\it\rho}^{\ast }\mathop{\sum }_{i=1}^{n_{sp}}Y_{i}\frac{\mathit{e}_{v,i}^{\ast }(T^{\ast })-\mathit{e}_{v,i}^{\ast }(T_{v}^{\ast })}{{\it\tau}_{i}^{\ast }}.\end{eqnarray}$$
Here
$Y_{i}$
is the species mass fraction and
${\it\tau}_{i}^{\ast }$
is the effective relaxation time for species
$i$
, which accounts for collisions between molecules of type
$i$
and all other types of molecules. For a mixture consisting of
$n_{sp}$
species, the effective relaxation time
${\it\tau}_{i}^{\ast }$
is given by (Millikan & White Reference Millikan and White1963)
$$\begin{eqnarray}{\it\tau}_{i}^{\ast }=\left[\mathop{\sum }_{j=1}^{n_{sp}}\frac{X_{i}}{{\it\tau}_{ij}^{\ast }}\right]^{-1},\end{eqnarray}$$
where
$X_{i}$
is the mole fraction of species
$i$
and
${\it\tau}_{ij}^{\ast }$
is the relaxation time for a dilute species
$i$
in an isothermal bath of molecules
$j$
; its value is computed using the correlations of Millikan & White (Reference Millikan and White1963). With the assumption of a single vibrational temperature, the heat flux vectors can be written in terms of the temperature gradients,
$$\begin{eqnarray}\displaystyle & \boldsymbol{q}_{v}^{\ast }=-k_{v}^{\ast }\boldsymbol{{\rm\nabla}}T_{v}^{\ast }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{q}_{tr}^{\ast }=-k^{\ast }\boldsymbol{{\rm\nabla}}T^{\ast }, & \displaystyle\end{eqnarray}$$
$k^{\ast }$
is the translational thermal conductivity and
$k_{v}^{\ast }$
is its vibrational counterpart. For each species, the thermal conductivity is evaluated using Euken’s relation (Vincenti & Kruger Reference Vincenti and Kruger1967), whereby the total conductivity is split into a translational/rotational part
$k_{i}^{\ast }$
and a vibrational part
$k_{v,i}^{\ast }$
: 
$$\begin{eqnarray}\displaystyle & k_{i}^{\ast }={\textstyle \frac{5}{2}}{\it\mu}_{i}^{\ast }c_{v,tr,i}^{\ast }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & k_{v,i}^{\ast }={\it\mu}_{i}^{\ast }c_{v,v,i}^{\ast }. & \displaystyle\end{eqnarray}$$
${\it\mu}_{i}^{\ast }$
is the dynamic viscosity of species
$i$
and
$c_{v,tr,i}^{\ast }$
is its translational specific heat at constant volume. The thermal conductivity of the mixture is then obtained using Wilke’s mixing rule (Wilke Reference Wilke1950). At low temperatures, the vibrational thermal conductivity vanishes and the Prandtl number reduces to a constant value of 0.737.
Separate viscosity models are needed for each species in order to apply Euken’s relation in (2.7). In this paper, Sutherland models for
$\text{N}_{2}$
and
$\text{O}_{2}$
were deduced by fitting Sutherland’s coefficient to the correlations of Cole & Wakeham (Reference Cole and Wakeham1985), which are valid up to 2100 K for nitrogen and 2500 K for oxygen. The viscosity of the mixture was then calculated using Wilke’s mixing rule. The mixture viscosity computed in this way differs by less than 4 % from the correlations of Kadoya, Matsunaga & Nagashima (Reference Kadoya, Matsunaga and Nagashima1985) over the temperature range 85–2000 K.
Equations (2.1) are made dimensionless using the following variables:
$$\begin{eqnarray}{\it\rho}=\frac{{\it\rho}^{\ast }}{{\it\rho}_{e}^{\ast }},\quad T=\frac{T^{\ast }}{T_{e}^{\ast }},\quad T_{v}=\frac{T_{v}^{\ast }}{T_{e}^{\ast }},\quad \boldsymbol{u}=\frac{\boldsymbol{u}^{\ast }}{U_{e}^{\ast }},\quad p=\frac{p^{\ast }}{{\it\rho}_{e}^{\ast }U_{e}^{\ast 2}},\quad t=\frac{t^{\ast }U_{e}^{\ast }}{L^{\ast }},\end{eqnarray}$$
$$\begin{eqnarray}c_{v,v}=\frac{c_{v,v}^{\ast }}{c_{p,tr}^{\ast }},\quad \boldsymbol{x}=\frac{\boldsymbol{x}^{\ast }}{L^{\ast }},\quad k=\frac{k^{\ast }}{k_{e}^{\ast }},\quad k_{v}^{\ast }=\frac{k_{v}^{\ast }}{k_{e}^{\ast }},\quad {\bf\tau}=\frac{{\bf\tau}^{\ast }}{{\it\mu}_{e}^{\ast }U_{e}^{\ast }/L^{\ast }},\quad Q=\frac{Q^{\ast }{\it\nu}_{e}^{\ast }}{{\it\rho}_{e}^{\ast }c_{p,tr}^{\ast }T_{e}^{\ast }U_{e}^{\ast 2}}.\end{eqnarray}$$
In these equations, subscript
$e$
refers to the edge conditions and
$L^{\ast }$
is a reference length. The vector
$\boldsymbol{x}$
has the components
$x$
,
$y$
and
$z$
corresponding to the streamwise, vertical and spanwise directions. The translational specific heat
$c_{p,tr}^{\ast }$
is a constant, since the gas composition is fixed. The vibrational source term
$Q$
is non-dimensionalized using the volumetric energy density
${\it\rho}_{e}^{\ast }c_{p,tr}^{\ast }T_{e}^{\ast }$
and the time scale
${\it\nu}_{e}^{\ast }/U_{e}^{\ast 2}$
.
After rewriting the governing equations (2.1) in terms of these non-dimensional parameters and subtracting the vibrational energy equation from the total energy equation, the result is
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\rho}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\rho}\boldsymbol{u})=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}\frac{\text{D}\boldsymbol{u}}{\text{D}t}+\boldsymbol{{\rm\nabla}}p=\frac{1}{\mathit{Re}_{L}}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\tau}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}c_{v,v}\frac{\text{D}T_{v}}{\text{D}t}=\frac{1}{\mathit{Re}_{L}{\it\sigma}_{e}}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(k_{v}\boldsymbol{{\rm\nabla}}T_{v})+\mathit{Re}_{L}Q, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}\frac{\text{D}T}{\text{D}t}-Ma^{2}({\it\gamma}_{e}-1)\frac{\text{D}p}{\text{D}t}=\frac{1}{\mathit{Re}_{L}{\it\sigma}_{e}}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(k\boldsymbol{{\rm\nabla}}T)+\frac{Ma^{2}({\it\gamma}_{e}-1)}{\mathit{Re}_{L}}({\bf\tau}:\boldsymbol{{\rm\nabla}}\boldsymbol{u})-\mathit{Re}_{L}Q,\qquad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\it\gamma}_{e}Ma^{2}p={\it\rho}ReT. & \displaystyle\end{eqnarray}$$
$Ma$
is the frozen Mach number in the free stream,
${\it\sigma}_{e}=c_{p,tr}^{\ast }{\it\mu}_{e}^{\ast }/k_{e}^{\ast }$
is the Prandtl number,
${\it\gamma}_{e}$
is the ratio of frozen specific heats, and
$\mathit{Re}_{L}$
is the Reynolds number based on the length scale
$L^{\ast }$
.
2.2. Base flow calculation
Since the boundary layer is in a state of thermal non-equilibrium, the velocity and temperature profiles cannot be computed using a similarity solution, and a suitable non-similar boundary layer solver is used instead. By applying Prandtl’s boundary layer scaling arguments (Schlichting & Gersten Reference Schlichting and Gersten2000) to (2.10) and assuming zero pressure gradient, one arrives at the result
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial x}(\bar{{\it\rho}}\bar{U})+\frac{\partial }{\partial y}(\bar{{\it\rho}}\bar{V})=0,\qquad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \bar{{\it\rho}}\bar{U}\frac{\partial \bar{U}}{\partial x}+\bar{{\it\rho}}\bar{V}\frac{\partial \bar{U}}{\partial y}=\frac{1}{\mathit{Re}_{L}}\frac{\partial }{\partial y}\left(\bar{{\it\mu}}\frac{\partial \bar{U}}{\partial y}\right),\qquad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \bar{{\it\rho}}\left(\bar{U}\frac{\partial \bar{T}}{\partial x}+\bar{V}\frac{\partial \bar{T}}{\partial y}\right)=\frac{1}{\mathit{Re}_{L}{\it\sigma}_{e}}\frac{\partial }{\partial y}\left(\bar{k}\frac{\partial \bar{T}}{\partial y}\right)-\mathit{Re}_{L}\bar{Q}+\frac{({\it\gamma}_{e}-1)Ma^{2}}{\mathit{Re}_{L}}\bar{{\it\mu}}\left(\frac{\partial \bar{U}}{\partial y}\right)^{2},\qquad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \bar{{\it\rho}}\bar{c}_{v,v}\left(\bar{U}\frac{\partial \bar{T}_{v}}{\partial x}+\bar{V}\frac{\partial \bar{T}_{v}}{\partial y}\right)=\frac{1}{\mathit{Re}_{L}{\it\sigma}_{e}}\frac{\partial }{\partial y}\left(\bar{k}_{v}\frac{\partial \bar{T}_{v}}{\partial y}\right)+\mathit{Re}_{L}\bar{Q},\qquad & \displaystyle\end{eqnarray}$$
$T_{w}^{\ast }=300~\text{K}$
, which is selected to facilitate comparisons with experiments in impulse facilities. Throughout this paper, the vibrational and translational temperatures are set equal in the free stream. This is often a good approximation for flows of air over slender bodies in high-enthalpy shock tunnels, but for some flows, especially those involving pure nitrogen, the vibrational temperature can be much greater than the translational temperature in the free stream. The impact of this phenomenon has been investigated by Bertolotti (Reference Bertolotti1998).
It is convenient to rewrite the boundary layer equations (2.11) in terms of the Blasius similarity variable
${\it\eta}$
:
$$\begin{eqnarray}{\it\eta}=y^{\ast }\sqrt{\frac{U_{e}^{\ast }}{{\it\nu}_{e}^{\ast }x^{\ast }}}.\end{eqnarray}$$
Although the flow is not self-similar, the use of similarity coordinates has the advantages of a fixed domain size, a fixed grid and removal of the singularity at the leading edge of the plate. By writing the continuity equation (2.11a
) in terms of the similarity variable, integrating with respect to
${\it\eta}$
and using integration by parts, one finds the relation
$$\begin{eqnarray}\bar{{\it\rho}}\bar{V}=-\sqrt{\frac{{\it\nu}_{e}^{\ast }x^{\ast }}{U_{e}^{\ast }}}\left(2\frac{\partial g}{\partial x^{\ast }}-\frac{1}{2x^{\ast }}[\bar{{\it\rho}}\bar{U}{\it\eta}-2g]\right),\end{eqnarray}$$
where, following the notation of Klunker & McLean (Reference Klunker and McLean1953),
$$\begin{eqnarray}g\equiv \frac{1}{2}\int _{0}^{{\it\eta}}(\bar{{\it\rho}}\bar{U})\,\text{d}{\it\eta}^{\prime }.\end{eqnarray}$$
Using (2.13) to eliminate the vertical velocity, the momentum and energy equations (2.11b )–(2.11d ) can be written in terms of the similarity variable as follows:
$$\begin{eqnarray}\displaystyle & \displaystyle \bar{{\it\rho}}\bar{U}x^{\ast }\frac{\partial \bar{U}}{\partial x^{\ast }}-\frac{\partial }{\partial {\it\eta}}\left(\bar{{\it\mu}}\frac{\partial \bar{U}}{\partial {\it\eta}}\right)-\frac{\partial \bar{U}}{\partial {\it\eta}}\left(g+2x^{\ast }\frac{\partial g}{\partial x^{\ast }}\right)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \bar{{\it\rho}}\bar{U}x^{\ast }\frac{\partial \bar{T}}{\partial x^{\ast }}-\frac{1}{{\it\sigma}_{e}}\frac{\partial }{\partial {\it\eta}}\left(\bar{k}\frac{\partial \bar{T}}{\partial {\it\eta}}\right)-({\it\gamma}_{e}\!-\!1)Ma^{2}\bar{{\it\mu}}\left(\frac{\partial \bar{U}}{\partial {\it\eta}}\right)^{2}\!+\mathit{Re}_{x}\bar{Q}-\frac{\partial \bar{T}}{\partial {\it\eta}}\left(\!g+2x^{\ast }\frac{\partial g}{\partial x^{\ast }}\right)=0, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \bar{{\it\rho}}\bar{U}\bar{c}_{v,v}x^{\ast }\frac{\partial \bar{T}_{v}}{\partial x^{\ast }}-\frac{1}{{\it\sigma}_{e}}\frac{\partial }{\partial {\it\eta}}\left(\bar{k}_{v}\frac{\partial \bar{T}_{v}}{\partial {\it\eta}}\right)-\mathit{Re}_{x}\bar{Q}-\bar{c}_{v,v}\frac{\partial \bar{T}_{v}}{\partial {\it\eta}}\left(g+2x^{\ast }\frac{\partial g}{\partial x^{\ast }}\right)=0. & \displaystyle\end{eqnarray}$$
$x^{\ast }\rightarrow 0$
, equations (2.15) reduce to the self-similar equations of Klunker & McLean (Reference Klunker and McLean1953) along with a decoupled vibrational energy equation: 
$$\begin{eqnarray}\displaystyle & \displaystyle g\frac{\partial \bar{U}}{\partial {\it\eta}}+\frac{\partial }{\partial {\it\eta}}\left(\bar{{\it\mu}}\frac{\partial \bar{U}}{\partial {\it\eta}}\right)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle g\frac{\partial \bar{T}}{\partial {\it\eta}}+\frac{1}{{\it\sigma}_{e}}\frac{\partial }{\partial {\it\eta}}\left(\bar{k}\frac{\partial \bar{T}}{\partial {\it\eta}}\right)+({\it\gamma}_{e}-1)Ma^{2}\bar{{\it\mu}}\left(\frac{\partial \bar{U}}{\partial {\it\eta}}\right)^{2}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle g\bar{c}_{v,v}\frac{\partial \bar{T}_{v}}{\partial {\it\eta}}+\frac{1}{{\it\sigma}_{e}}\frac{\partial }{\partial {\it\eta}}\left(\bar{k}_{v}\frac{\partial \bar{T}_{v}}{\partial {\it\eta}}\right)=0. & \displaystyle\end{eqnarray}$$
$10^{-7}$
.
Examples of temperature and velocity profiles calculated using this method are shown in figure 1. For validation, these results are compared with Navier–Stokes simulations, which were conducted using the STABL stability software package developed at the University of Minnesota (Johnson et al.
Reference Johnson, Seipp and Candler1998; Wright, Candler & Bose Reference Wright, Candler and Bose1998; Johnson Reference Johnson2000); this software package includes a shock-capturing finite volume solver that models non-equilibrium chemically reacting flows. Navier–Stokes simulations shown here and throughout this paper were conducted using a grid of 450 × 450 points, with exponential clustering of points at the wall and near the leading edge of the plate. Convergence was checked on a finer grid of 600 streamwise
$\times$
650 wall-normal points, which led to no significant differences in the flow profiles or stability calculations based on them.
Figure 1. Comparison of boundary layer profiles at
$\mathit{Re}_{x}=2.25\times 10^{6}$
from STABL Navier–Stokes solver (symbols) and boundary layer code (solid lines). Profiles are velocity (▿), translational/rotational temperature (▫) and vibrational temperature (♢). (a) Low enthalpy,
$P_{e}^{\ast }=5~\text{kPa}$
,
$Ma=5$
,
$T_{w}^{\ast }=300~\text{K}$
,
$T_{e}^{\ast }=T_{v,e}^{\ast }=300~\text{K}$
,
$x=0.42~\text{m}$
. (b) High enthalpy,
$P_{e}^{\ast }=20~\text{kPa}$
,
$Ma=5$
,
$T_{w}^{\ast }=300~\text{K}$
,
$T_{e}^{\ast }=T_{v,e}^{\ast }=1500~\text{K}$
,
$x^{\ast }=0.70~\text{m}$
.
Figure 1(a) is a fairly low-enthalpy flow that is representative of a slender body in free flight at Mach 5, with
$T_{e}^{\ast }=T_{v,e}^{\ast }=300~\text{K}$
. At
$\mathit{Re}_{x}=2.25\times 10^{6}$
, the vibrational temperature remains close to its initial value and the translational temperature distribution is that of a self-similar frozen flow. Figure 1(b) is a much higher-enthalpy case that is representative of a shock tunnel, with
$T_{e}^{\ast }=T_{v,e}^{\ast }=1500~\text{K}$
. In this case the rate of vibrational energy transfer is more significant, so that by
$\mathit{Re}_{x}=2.25\times 10^{6}$
the vibrational temperature is fairly close to the equilibrium temperature distribution. As
$\mathit{Re}_{x}$
increases further, the vibrational and translational temperature profiles both approach a single equilibrium temperature profile. For both cases shown in figure 1, good agreement between the boundary layer solver and the Navier–Stokes solver is seen, and similar agreement has been found at other Reynolds numbers.
2.3. Stability calculations
The stability calculations presented in this paper are based on the locally parallel linearized Navier–Stokes equations in which gradients of the mean flow in the streamwise direction and the mean vertical velocity are neglected. According to this scheme, the three velocities (
$u,v,w$
), the pressure
$p$
and the two temperatures
$T$
and
$T_{v}$
are divided into mean and fluctuating components as follows:
$$\begin{eqnarray}\left(\begin{array}{@{}c@{}}u\\ v\\ w\\ p\\ T\\ T_{v}\end{array}\right)=\left(\begin{array}{@{}c@{}}\bar{U}(y)\\ 0\\ \bar{W}(y)\\ \bar{P}\\ \bar{T}(y)\\ \bar{T}_{v}(y)\end{array}\right)+\text{Re}\left\{\left(\begin{array}{@{}c@{}}\hat{u} (y)\\ \hat{v}(y)\\ {\hat{w}}(y)\\ \hat{p}(y)\\ \hat{{\it\theta}}(y)\\ \hat{{\it\theta}}_{v}(y)\end{array}\right)\text{e}^{\text{i}({\it\alpha}x+{\it\beta}z-{\it\omega}t)}\right\}.\end{eqnarray}$$
In this system of equations the overbars signify mean flow variables, hats designate the complex amplitude coefficients of the fluctuations, and
$\text{Re}$
signifies the real part of a complex variable. In all cases a spatial analysis is performed with real frequency
${\it\omega}$
, real spanwise wavenumber
${\it\beta}$
and complex streamwise wavenumber
${\it\alpha}$
. The wavenumbers are non-dimensionalized using the local Blasius boundary layer thickness
${\it\delta}=\sqrt{{\it\nu}_{e}^{\ast }x^{\ast }/U_{e}^{\ast }}$
and the frequency is non-dimensionalized using the time scale
${\it\delta}/U_{e}^{\ast }$
. The fluctuations (2.17) are substituted into the Navier–Stokes equations (2.10) and linearized, and the result is expressed as a system of 10 first-order differential equations of the form:
$$\begin{eqnarray}\frac{\text{d}\hat{\boldsymbol{q}}}{\text{d}y}=\unicode[STIX]{x1D63C}\hat{\boldsymbol{q}},\end{eqnarray}$$
where
$\unicode[STIX]{x1D63C}$
is a 10 × 10 matrix whose coefficients are given in appendix A, and
$\hat{\boldsymbol{q}}$
is the vector of disturbance variables,
$$\begin{eqnarray}\hat{\boldsymbol{q}}=\left(\hat{u} ,\frac{\text{d}\hat{u} }{\text{d}y},\hat{v},\hat{p},\hat{{\it\theta}},\frac{\text{d}\hat{{\it\theta}}}{\text{d}y},{\hat{w}},\frac{\text{d}{\hat{w}}}{\text{d}y},\hat{{\it\theta}}_{v},\frac{\text{d}\hat{{\it\theta}}_{v}}{\text{d}y}\right)^{\text{T}}.\end{eqnarray}$$
The boundary conditions at the wall are taken to be
$$\begin{eqnarray}\hat{u} =\hat{v}={\hat{w}}=\hat{{\it\theta}}=\hat{{\it\theta}}_{v}=0,\quad y=0.\end{eqnarray}$$
At the wall the fluctuations in streamwise, vertical and spanwise velocities
$\hat{u}$
,
$\hat{v}$
and
${\hat{w}}$
are set to zero in accordance with the no-slip condition, and the fluctuations in translational and vibrational temperatures
$\hat{{\it\theta}}$
and
$\hat{{\it\theta}}_{v}$
are set to zero under the assumption that the thermal inertia of the wall is high at the frequencies of interest. Ma & Zhong (Reference Ma and Zhong2003) compared isothermal and adiabatic boundary conditions for the disturbances, and they found that the spatial growth rates differed only slightly, with adiabatic disturbances being less stable. Along with the boundary conditions (2.20), the normalization condition
$\hat{p}(0)=1$
is used.
In the free stream the disturbances are required to remain bounded:
$$\begin{eqnarray}\hat{u} ,\hat{v},{\hat{w}},\hat{p},\hat{{\it\theta}},\hat{{\it\theta}}_{v}<\infty ,\quad y\rightarrow \infty .\end{eqnarray}$$
These boundary conditions are enforced using the method of Mack (Reference Mack, Alder, Fernback and Rotenberg1965) in which the asymptotic behaviour of the 10 fundamental solutions of (2.18) is determined in the free stream and the five exponentially increasing solutions are rejected. The remaining solutions are integrated towards the wall using a fourth-order Runge–Kutta method, and a linear combination of these fundamental solutions is taken to satisfy the normalization condition
$\hat{p}(0)=1$
and all of the boundary conditions except
$\hat{{\it\theta}}(0)=0$
. The remaining boundary condition on the temperature perturbation is satisfied by a complex eigenvalue search, which is carried out using the secant method. During the Runge–Kutta integration from the free stream to the wall, Gram–Schmidt orthonormalization is employed intermittently in order to control the parasitic growth of numerical errors (Conte Reference Conte1966; Garg Reference Garg1980; Davey Reference Davey1983). Orthonormalization is applied whenever the norm of one of the fundamental solutions exceeds 106.
The use of the asymptotic boundary condition (2.21) provides a significant advantage over the Dirichlet approximation that is often used. If Dirichlet boundary conditions are used, then the domain must be much larger than the boundary layer thickness in order to capture the asymptotic behaviour of the eigenfunctions. Because of this large domain, a highly stretched grid is needed to achieve sufficient resolution of the instability within the boundary layer and near the wall. In contrast, if the asymptotic boundary condition is used, then the domain needs only to be slightly larger than the boundary layer thickness and hence the grid points are used more efficiently. As will be seen, this feature is especially useful when simulating hypersonic boundary layers with very cold walls; at these conditions the unstable disturbances can travel supersonically relative to the mean flow and their eigenfunctions are nearly non-decaying in the free stream. This non-decaying behaviour can be captured only by using the asymptotic boundary conditions or a very large domain, and the eigenvalues can be affected if Dirichlet boundary conditions are used with too small a domain.
For all simulations the number of grid points in the boundary layer is 1000. This number has been found to produce grid-independent results for Reynolds numbers up to approximately
$Re=6000$
, where
$Re=U_{e}^{\ast }{\it\delta}/{\it\nu}_{e}^{\ast }=\sqrt{\mathit{Re }_{x}}$
is the Reynolds number based on boundary layer thickness. The size of the domain used depends on the Mach number, and is selected to be a factor of 1.5–5 larger than
${\it\delta}_{99}$
, which is the boundary thickness at which the mean streamwise velocity reaches 99 % of the edge velocity. The results obtained are unchanged if the domain is made larger than this, even when the eigenfunctions are slowly decaying outside of the boundary layer.
3. Results
3.1. Model verification
For verification of the stability analysis, the results of the present method are compared with those of Fedorov & Tumin (Reference Fedorov and Tumin2011) for a perfect gas. Figure 2 overlays the present results (symbols) against figure 10 from that paper (solid lines). The flow is a perfect gas with constant specific heats with
$Ma=4.2$
, Prandtl number of 0.72, Sutherland constant of 110.4 K and an adiabatic wall. Excellent agreement is seen for both growth rates (a) and phase speeds (b), where
$c_{r}$
is the real part of the complex phase speed
$c={\it\omega}/{\it\alpha}$
.
Figure 2. Comparison of present stability analysis (symbols) with the results of Fedorov & Tumin (Reference Fedorov and Tumin2011) (solid lines). Flow is a perfect gas with
$Ma=4.2$
,
$T_{o}^{\ast }=300~\text{K}$
,
$Re=2000$
, adiabatic wall. Waves are 2D with
${\it\beta}=0$
; (a) growth rates; (b) phase speeds.
As a second validation case, a higher-enthalpy flow is selected in which the variation of specific heats is significant and the base flow is not in equilibrium. Figure 3 compares our spatial growth rates with those from the PSE-Chem portion of the STABL software package (Johnson et al.
Reference Johnson, Seipp and Candler1998; Wright et al.
Reference Wright, Candler and Bose1998; Johnson Reference Johnson2000). The PSE-Chem computations were conducted using 300 points in the wall-normal direction, and convergence was checked by doubling the number of points, which affected the growth rates by at most 1 %. In performing the PSE-Chem analysis, the options for PSE marching and non-parallel effects were turned off so as to make as direct a comparison as possible with our results. However, chemical and vibrational non-equilibria were retained in the PSE-Chem calculation, with the chemistry model consisting of five species:
$\text{N}_{2}$
,
$\text{O}_{2}$
, NO, N and O. The free-stream concentrations of these species were set to 78 %
$\text{N}_{2}$
and 22 %
$\text{O}_{2}$
by mole, which matches the present linear stability calculation.
Figure 3. Comparison of spatial growth rates from STABL PSE-Chem (symbols) and present LST analysis (solid lines). (a) Low enthalpy,
$P_{e}^{\ast }=5~\text{kPa}$
,
$Ma=5$
,
$T_{w}^{\ast }=300~\text{K}$
,
$T_{e}^{\ast }=T_{v,e}^{\ast }=300~\text{K}$
. (b) High enthalpy,
$P_{e}^{\ast }=20~\text{kPa}$
,
$Ma=5$
,
$T_{w}^{\ast }=300~\text{K}$
,
$T_{e}^{\ast }=T_{v,e}^{\ast }=1500~\text{K}$
. Only 2D waves of
${\it\beta}=0$
are included.
Figure 4. High-enthalpy comparison between current method and STABL software:
$P_{e}^{\ast }=20~\text{kPa}$
,
$Ma=5$
,
$T_{w}^{\ast }=300~\text{K}$
,
$T_{e}^{\ast }=T_{v,e}^{\ast }=2500~\text{K}$
. (a) Comparison of boundary layer profiles at
$Re=1500$
(
$x^{\ast }=1.2~\text{m}$
) from STABL Navier–Stokes solver (symbols) and boundary layer code (solid lines). Profiles are velocity (▿), translational/rotational temperature (▫) and vibrational temperature (♢). (b) Comparison of spatial growth rates from STABL PSE-Chem (symbols) and present LST analysis (solid lines). Only 2D waves of
${\it\beta}=0$
are included.
Figure 3(a) is a relatively low-enthalpy flow that is representative of the free flight of a slender body. At the edge of the boundary layer both translational and vibrational temperatures are 300 K. Figure 3(b) is a higher-enthalpy flow that is representative of a shock tunnel, with
$T_{e}^{\ast }=T_{v,e}^{\ast }=1500~\text{K}$
. The base flows corresponding to figure 3 are given in figure 1 for
$Re=1500$
. Very close agreement between the two methods is seen, though slightly higher growth rates are predicted by the present model based on conventional linear stability theory (LST). The small differences are greater near the leading edge of the plate and may be caused by the weak shock that develops through the viscous–inviscid interaction: this effect is included in the PSE-Chem analysis, but is excluded from the present LST model. This explanation is consistent with the findings of Chang et al. (Reference Chang, Malik and Hussaini1990), who showed that, when the shock is close to the boundary layer edge, there is a slight stabilization of both the first and second modes. Differences in transport properties also contribute to the discrepancies in growth rates.
This chemically frozen stability analysis might be expected to perform poorly for the very high-enthalpy conditions that are achievable in reflected shock tunnel facilities. For instance, in a uniform volume of air at 20 kPa, the enthalpy of chemical equilibrium air first exceeds that of chemically frozen air by 10 % at approximately 2600 K. This temperature can be significantly exceeded in shock tunnel facilities. However, we have found that, for slender bodies in air, the stability characteristics at such high temperatures can still be predicted reasonably well under the assumption of a chemically inert gas.
Figure 4 shows mean boundary layer profiles and stability calculations for a high-enthalpy boundary layer with an edge temperature of 2500 K and Mach number of
$Ma=5$
. This corresponds to a stagnation enthalpy of approximately
$15~\text{MJ}~\text{kg}^{-1}$
, which is close to the maximum stagnation enthalpy achievable in reflected shock tunnel facilities like the T5 facility at Caltech. For this case, the peak temperature within the boundary layer exceeds 3500 K, which is well above the limit at which dissociation of
$\text{O}_{2}$
becomes significant. For instance, the chemically reacting Navier–Stokes simulations reveal that at
$Re=1500$
the mass fraction of atomic oxygen reaches a peak value of 0.72 % in the interior of the boundary layer. The mean boundary layer profiles (figure 4
a) are in reasonable agreement with the Navier–Stokes simulations from STABL, although the boundary layer code produces noticeably greater temperatures by neglecting chemistry. The difference in peak temperature is large at the high Reynolds number shown, but is much smaller at lower Reynolds numbers, where the chemical reactions have not progressed as far towards equilibrium.
In figure 4(b), the chemically frozen stability predictions from the present method are compared with the chemically reacting results from STABL. The agreement is reasonably good, especially considering that we have neglected chemical reactions and are using different transport property models. These results suggest that, despite the high temperatures in the interior of the boundary layer, the present chemically frozen analysis can be used to study the stability of flows of air over slender bodies in high-enthalpy shock tunnels without incurring serious errors. This finding can be explained in part by the analysis of Fujii & Hornung (Reference Fujii and Hornung2003), which showed that, in air at 1000–3000 K, the absorption of sound waves caused by thermochemical non-equilibrium occurs at frequencies much lower than those of second-mode disturbances in hypervelocity boundary layers. For blunt geometries, however, the extremely high temperatures and low flow velocities that occur behind the bow shock may lead to more significant chemical effects and the frozen analysis may no longer be adequate.
3.2. Mode characteristics and terminology
To demonstrate the mode characteristics and introduce some terminology, figure 5 shows the disturbance spectrum for a high-Mach-number, low-enthalpy flow. Spatial growth rates (a) and phase speeds (b) are shown for
$Ma=4.5$
and
$Re=\sqrt{U_{e}^{\ast }x^{\ast }/{\it\nu}_{e}^{\ast }}=2000$
with a cold free stream (
$T_{e}^{\ast }=65.15~\text{K}$
) and an adiabatic wall. This same case was computed by Ma & Zhong (Reference Ma and Zhong2003), and our results are identical to theirs to within plotting accuracy.
Figure 5. Spatial growth rates (a) and phase speeds (b) for low-enthalpy flow with
$Ma=4.5$
,
$Re=2000$
,
$T_{e}^{\ast }=65.15~\text{K}$
,
$P_{e}^{\ast }=728~\text{Pa}$
, adiabatic wall. Flow conditions modelled after Ma & Zhong (Reference Ma and Zhong2003). Only 2D waves with
${\it\beta}=0$
are shown.
Figure 6. Spatial growth rates (a) and phase speeds (b) for high-enthalpy flow with
$Ma=4.5$
,
$Re=2000$
,
$T_{e}^{\ast }=1500~\text{K}$
,
$T_{w}^{\ast }=300~\text{K}$
,
$P_{e}^{\ast }=10~\text{kPa}$
. Only 2D waves with
${\it\beta}=0$
are shown.
The disturbance spectrum features an infinite sequence of modes that ‘cut in’ as the frequency increases. By ‘cut in’, it is meant that these modes transition from the continuous spectrum to the discrete spectrum. Each mode cuts in with a dimensionless phase speed of
$1+1/Ma$
, which is the speed of fast acoustic waves in the free stream, hence these modes are termed ‘fast modes’ and designated
$F_{1}$
,
$F_{2}$
, etc. in order of increasing cut-in frequency. Each additional mode that cuts in is characterized by an additional zero in the real part of the pressure eigenfunction, with mode
$F_{n}$
having
$n-1$
zeros. There is also one mode that has a phase speed of
$1-1/Ma$
at
${\it\omega}=0$
; this mode is termed the ‘slow mode’ and labelled
$S_{1}$
.
At the locations where the phase speeds of the fast and slow modes cross, there is a possibility of synchronization, which may lead to increased growth rate and instability. At these crossing points, the phase speeds of the slow and fast modes are equal and their eigenfunctions are nearly identical. This synchronism has been analysed by Fedorov & Khokhlov (Reference Fedorov and Khokhlov2001), who proposed a reduced-order model of the dispersion relation in the vicinity of the synchronism region, which captures the topology of the spectrum that is found there. These synchronous interactions between fast and slow modes are frequently called the ‘Mack modes’ and are labelled Mack’s second mode, third mode, etc. in order of increasing frequency. As discussed by Fedorov & Tumin (Reference Fedorov and Tumin2011), the Mack modes are not ‘modes’ in the mathematical sense as are modes
$F_{n}$
and
$S_{1}$
, but since this terminology is in common use it will be used throughout in this paper.
Lastly, we point out that the lowest-frequency unstable region in figure 5, which was labelled the ‘first mode’ by Mack, does not involve a synchronism between fast and slow acoustic modes. This instability is the compressible counterpart to the Tollmien–Schlichting waves that are found in incompressible flows, and its instability characteristics depend differently on Mach number, Reynolds number and wall temperature compared to the acoustic instabilities.
Figure 7. Effect of wall cooling on the spatial growth rate (a) and phase speed (b). For all cases,
$Ma=5$
,
$Re=2000$
,
$T_{w}^{\ast }=300~\text{K}$
,
$P_{e}^{\ast }=10~\text{kPa}$
. Lines with symbols designate an additional mode that appears when
$T_{w}^{\ast }\ll T_{e}^{\ast }$
. Only 2D waves with
${\it\beta}=0$
are shown.
For comparison with the low-enthalpy adiabatic case given in figure 5, the stability diagram for a high-enthalpy flow with a cold wall is shown in figure 6. Here the Mach number is again 4.5 and the Reynolds number is
$Re=2000$
, but the free-stream temperature is now 1500 K and the wall temperature is 300 K. These conditions are representative of a moderate-enthalpy flow in a reflected shock tunnel. By comparing figures 5 and 6, it is seen that the maximum growth rate is more than a factor of 2 greater for the cold-wall case, which is consistent with the well-known fact that the second-mode instability is destabilized by wall cooling. Also the first-mode instability is absent in the case with a cold wall, which again is well-known behaviour for hypersonic flows. An additional difference is that in figure 5 the mode
$S_{1}$
is unstable while in figure 6 the mode
$F_{1}$
is unstable. This behaviour has also been observed by Fedorov & Khokhlov (Reference Fedorov and Khokhlov2001), who demonstrated that the acoustic instabilities are associated with a pair of branch points in the complex plane whose location depends on the Reynolds number, Mach number, frequency, wall temperature and other parameters. As the frequency is increased, either
$F_{1}$
or
$S_{1}$
can be unstable depending on whether the dispersion curve passes above or below these branch points. Further discussion regarding this behaviour is available from Fedorov & Tumin (Reference Fedorov and Tumin2011).
The above-mentioned differences between figures 5 and 6 are well known, but we also observe two new features in figure 6 that have not been widely reported in the literature. The first is an abrupt change in slope of the second-mode unstable region at approximately
${\it\omega}=0.5$
, which leads to a wider range of unstable frequencies for the cold-wall case. This behaviour is also visible in a few of the stability calculations reported by Klentzman & Tumin (Reference Klentzman and Tumin2013), which involved high-enthalpy flows of oxygen. The second new feature is that the mode labelled
$F_{1}$
ceases to exist for frequencies greater than approximately
${\it\omega}=0.62$
, but a new related mode (dashed line) appears at
${\it\omega}=0.54$
. These two modes have nearly the same phase speed over the frequency range
$0.54<{\it\omega}<0.62$
for which they both exist. Note that a similar phenomenon involving the disappearance of one mode and the appearance of another occurs for
$F_{1}$
at
${\it\omega}=0.32$
, but the effect is hardly visible at the scale shown in the figure.
3.3. Shape of dispersion curve
We first investigate the sharp change in slope of the dispersion curve observed in figure 6. Figure 7 shows the growth rates and phase speeds for flow at
$Ma=5$
and
$Re=2000$
. Three cases are shown with free-stream temperatures of 300, 1000 and 1500 K and with the wall temperature fixed at 300 K in all cases. As the ratio
$T_{w}^{\ast }/T_{e}^{\ast }$
decreases, a kink appears in the high-frequency end of the growth rate curve and the width of the unstable region increases. For each case, the frequency at which the kink occurs is marked by a vertical dashed line. From the intersection between these vertical lines and the phase speed diagram below, one can see that the change in slope of the growth rate curve occurs precisely when the dimensionless phase speed falls below
$1-1/Ma$
. That is, the change in shape of the dispersion curve occurs when the unstable modes travel supersonically relative to the free stream. Supersonic modes can exist for all wall temperature conditions, but they only become unstable when the wall is highly cooled, usually (but not necessarily) below the free-stream temperature. This behaviour was also encountered by Chang et al. (Reference Chang, Vinh and Malik1997) in their study of chemically reacting, high-enthalpy flow over a wedge, and the existence of supersonic unstable modes was also briefly mentioned by Mack (Reference Mack1969, Reference Mack1987) for the inviscid case.
The eigenfunctions corresponding to two values of
${\it\omega}$
from figure 7 are shown in figure 8 with
$T_{e}^{\ast }=1500~\text{K}$
. In each case, the locations of the critical layer
$c_{r}=\bar{U}$
and the sonic lines
$c_{r}=\bar{U}\pm \bar{a}$
are indicated by horizontal dashed lines. For
${\it\omega}=0.35$
(just before the kink in the dispersion curve), the phase speed is above
$c_{r}>1-1/Ma$
and only one sonic line exists, but for
${\it\omega}=0.45$
(just after the kink) there are two sonic lines because
$c_{r}<1-1/Ma$
. Between these two lines the disturbance is travelling subsonically relative to the fluid, but outside of these lines the disturbance travels supersonically.
Figure 8. Eigenfunctions corresponding to two points along the dispersion curve in figure 7, for
$Re=2000$
,
$T_{e}^{\ast }=1500~\text{K}$
,
$T_{w}^{\ast }=300~\text{K}$
,
$Ma=5$
: (a)
${\it\omega}=0.35$
, subsonic mode; (b)
${\it\omega}=0.45$
, supersonic mode.
From figure 8 one can see that for
${\it\omega}=0.45$
the eigenfunctions decay very slowly in the free stream. This is highlighted by the contours of the temperature disturbance given in figure 9. In the interior of the boundary layer, the structure is similar for
${\it\omega}=0.35$
and 0.45, and a large peak in temperature is visible at the critical layer. For
${\it\omega}=0.35$
, the waves decay rapidly in the free stream, but for
${\it\omega}=0.45$
the second sonic line acts as a turning point, leading to oscillatory waves that are radiated into the free stream. The eigenfunction behaviour in the free stream is similar for the pressure and velocity fluctuations.
Figure 9. Contours of temperature fluctuations at two points along the dispersion curve in figure 7 for
$Re=2000$
,
$T_{e}^{\ast }=1500~\text{K}$
,
$T_{w}^{\ast }=300~\text{K}$
,
$Ma=5$
: (a)
${\it\omega}=0.35$
, subsonic mode; (b)
${\it\omega}=0.45$
, supersonic mode.
This behaviour can be anticipated from the compressible Rayleigh equation for the pressure (see Lees & Lin Reference Lees and Lin1946; Lees & Reshotko Reference Lees and Reshotko1962; Mack Reference Mack1969)
$$\begin{eqnarray}\frac{\text{d}^{2}\hat{p}}{\text{d}y^{2}}-\left(\frac{2\bar{U}^{\prime }}{\bar{U}-c}-\frac{\bar{T}^{\prime }}{\bar{T}}\right)\frac{\text{d}\hat{p}}{\text{d}y}-{\it\alpha}^{2}\left(1-\frac{Ma^{2}(\bar{U}-c)^{2}}{\bar{T}}\right)\hat{p}=0.\end{eqnarray}$$
In the free stream this equation reduces to
$$\begin{eqnarray}\frac{\text{d}^{2}\hat{p}}{\text{d}y^{2}}-{\it\alpha}^{2}[1-Ma^{2}(1-c)^{2}]\hat{p}=0,\end{eqnarray}$$
which has the solution
$$\begin{eqnarray}\hat{p}=\exp ({\it\alpha}\sqrt{1-Ma^{2}(1-c)^{2}}y).\end{eqnarray}$$
Here the branch of the square root is chosen such that the real part of the exponential’s argument is negative, so the pressure is bounded. Since
$|c_{i}|\ll |c_{r}|$
and
$|{\it\alpha}_{i}|\ll |{\it\alpha}_{r}|$
, the pressure is mainly exponentially decaying if
$c_{r}>1-1/Ma$
with slight oscillations arising from the imaginary parts of
${\it\alpha}$
and
$c$
. Conversely, for
$c_{r}<1-1/Ma$
the solution is mainly oscillatory, with slight damping caused by the imaginary parts of
${\it\alpha}$
and
$c$
. Therefore, slowly decaying eigenfunctions can always be expected for supersonic modes at high Reynolds numbers.
Figure 7 showed that cases which develop a supersonic unstable mode also experience a discontinuity in the dispersion curve. As the frequency increases, the fast mode (labelled
$F_{1}^{+}$
in the figure) disappears and a new mode (labelled
$F_{1}^{-}$
) appears. To further investigate this effect, the real and imaginary parts of the complex phase velocity are plotted for mode
$F_{1}$
in figure 10. The data are exactly the same as in figure 7, except that now the real and imaginary parts of the phase speed are plotted with
${\it\omega}$
as a parameter along the curves. Also shown in the figure are three thick black lines which mark the branch cuts corresponding to the continuous spectra; the two nearly horizontal branches are the acoustic branch cuts, while the vertical branch near
$c_{r}=1$
(which is actually two overlapping branches) corresponds to vorticity and entropy waves. Further discussion about the location of these branch cuts and methods for their computation are available from Balakumar & Malik (Reference Balakumar and Malik1992). It should be noted that the continuous spectrum is a (weak) function of
${\it\omega}$
, and the spectra plotted in figure 10 are evaluated at a single value of
${\it\omega}=0.4$
.
Figure 10. Map of real and imaginary parts of phase speed in the complex plane, for
$Ma=5$
and
$Re=2000$
. Solid black lines correspond to the continuous spectra evaluated at
${\it\omega}=0.4$
. Lines with symbols designate an additional mode that appears when
$T_{w}^{\ast }\ll T_{e}^{\ast }$
.
At
${\it\omega}=0$
, the mode
$F_{1}^{+}$
is a neutral wave with phase speed
$c_{r}=1+1/Ma$
, which originates at the branch point of the fast acoustic branch cut. As indicated by the arrow in figure 10, as the frequency increases, the mode’s phase speed reduces and it becomes damped. As the mode crosses the branch cut at
$c_{r}\sim 1$
, the growth rate suffers a small jump while the real part of
${\it\alpha}$
remains nearly continuous. This phenomenon is discussed further by Fedorov & Khokhlov (Reference Fedorov and Khokhlov2001), who showed that the synchronism between the eigenvalue and the vorticity/entropy waves of the continuous spectrum causes the mode
$F_{1}$
to be especially receptive to entropy spots and vortical free-stream disturbances at this condition. Additional analysis of the receptivity to entropy disturbances is available from Fedorov & Tumin (Reference Fedorov and Tumin2003), who solved the initial value problem of a temperature spot located at the edge of a boundary layer.
As the frequency is increased further, the mode becomes unstable. For the case with
$T_{e}^{\ast }=300~\text{K}$
, the mode becomes stable again before the phase speed falls below
$1-1/Ma$
, so the mode passes below the slow acoustic branch point in the complex plane. This is typical of boundary layers that are adiabatic or have wall temperatures greater than the free-stream temperature. In contrast, for the cases with
$T_{e}^{\ast }=1000$
and 1500 K, the mode remains unstable for
$c_{r}<1-1/Ma$
and therefore passes above the branch point of the slow acoustic waves.
Because the mode with
$T_{e}^{\ast }=1000$
or 1500 K passes above the branch point, it must eventually cross the slow acoustic branch cut. When this happens, the mode coalesces with the branch cut and ceases to be a discrete mode; however, a new mode emerges from the other side of the branch cut. This new mode emerges from the branch cut at a slightly lower frequency than the one at which the coalescing mode disappears. The modes on the top and bottom of the branch cut are distinguished by the labels
$F_{1}^{+}$
and
$F_{1}^{-}$
, similar to the notation used by Fedorov & Tumin (Reference Fedorov and Tumin2011) to describe modes crossing the branch cut at
$c_{r}\sim 1$
.
This crossing of the branch cut results in a synchronism between the instability mode and acoustic disturbances in the free stream. This is similar to the synchronism with vorticity/entropy waves that occurs when the branch cut at
$c_{r}=1$
is crossed (Fedorov & Khokhlov Reference Fedorov and Khokhlov2001; Fedorov et al.
Reference Fedorov, Ryzhov, Soudakov and Utyuzhnikov2013), which has been found to produce increased receptivity to vorticity and entropy spottiness. The synchronism identified here is expected to cause a similar effect, but with enhanced receptivity to acoustic disturbances in the free stream rather than vorticity/entropy disturbances. However, unlike the vorticity/entropy synchronism, which occurs upstream of the lower neutral branch, the synchronism with free-stream acoustic waves takes place downstream of the upper neutral branch. As a result, modes excited in this manner by free-stream acoustic waves are unlikely to experience amplification. Nevertheless, the synchronism can still affect the downstream development of amplified waves, which may be important during the nonlinear stages of transition.
3.4. Effect of Reynolds and Mach numbers
In figure 7, it was shown that the supersonic unstable modes only appear when the wall is sufficiently cold. In this section we investigate also how the Reynolds and Mach numbers influence these supersonic modes. Figure 11(a) compares the spatial growth rates for a high-enthalpy boundary layer at several different Reynolds numbers. For low enough Reynolds numbers, the mode remains subsonic throughout the entire unstable region, but for higher
$Re$
the ‘tail’ on the dispersion curve appears and grows considerably larger as the Reynolds number is further increased.
Figure 11. (a) Effect of Reynolds number on the shape of the dispersion curve for a cold-wall case:
$Ma=5$
,
$T_{e}^{\ast }=1500~\text{K}$
,
$T_{w}^{\ast }=300~\text{K}$
,
$P_{e}^{\ast }=20~\text{kPa}$
. (b) The
$N$
factor curve corresponding to (a). The heavy black line is the complete
$N$
factor curve. The dot-dashed red line includes only values of
${\it\omega}<0.4$
. Only 2D waves with
${\it\beta}=0$
are included.
As can be noted from figure 11(a), at higher Reynolds numbers the supersonic modes contribute significantly to the area under the dispersion curve, which is related to the
$N$
factor distribution. This brings into question whether the supersonic modes significantly influence the
$N$
factors. Figure 11(b) shows the
$N$
factor envelope curves for a cold-wall case with
$Ma=5$
and
$T_{e}^{\ast }=1500~\text{K}$
. Two envelope curves are included in this figure: the dashed red line excludes growth rates having frequencies above
${\it\omega}=0.4$
, while the solid black line is the full
$N$
factor curve. The difference between the two curves then indicates the contribution of supersonic modes. Clearly, the supersonic modes do not contribute much to the
$N$
factor until the Reynolds number exceeds approximately 1500. However, at this point the
$N$
factor is already fairly high, reaching the level of 5–10 at which transition is typically observed experimentally. This suggests that the supersonic modes do not significantly increase the level of amplification, despite their influence on the shape and behaviour of the dispersion curve.
Lastly, the effect of the Mach number on the stability of highly cooled flows is considered. Figure 12 shows dispersion curves for several Mach numbers with
$T_{e}^{\ast }=1500~\text{K}$
and
$T_{w}^{\ast }=300~\text{K}$
. The left panel is for a lower Reynolds number of
$Re=1000$
, while the right one is for
$Re=2000$
. These plots reveal that the supersonic unstable modes first appear and experience the largest growth rate at approximately
$Ma=5$
.
Figure 12. Dispersion curves for cold-wall flows (
$T_{w}^{\ast }/T_{e}^{\ast }=0.2$
) for various Mach numbers: (a)
$Re=1000$
; (b)
$Re=2000$
.
3.5. Vibrational non-equilibrium effects
In this section, the influence of thermal non-equilibrium on the growth of disturbances is examined. Vibrational non-equilibrium primarily enters the stability analysis in two ways. First, the mean profiles of temperature, density and (to a smaller extent) velocity are influenced by the exchange of energy between the translational and vibrational modes. This indirectly affects the disturbances through changes to the mean flow profile. The second effect of vibrational non-equilibrium is the attenuation of sound waves caused by the phase lag between the kinetic energy of molecules and their internal energy modes (Lighthill Reference Lighthill, Batchelor and Davies1956; Clarke & McChesney Reference Clarke and McChesney1964; Vincenti & Kruger Reference Vincenti and Kruger1967; Fujii & Hornung Reference Fujii and Hornung2003).
The first effect of vibrational non-equilibrium, changes to the base flow, can be assessed by conducting a stability analysis for three different base flows: a frozen flow with vibrational energy neglected, a non-equilibrium flow with finite rates of vibrational energy transfer, and an equilibrium flow in which vibrational energy is assumed to be exchanged at an infinite rate. In computing the stability of these three base flows, fully non-equilibrium disturbances are modelled. The resulting stability diagram is shown in figure 13(a). As could be expected, the non-equilibrium solution agrees best with the frozen solution for small Reynolds numbers, and tends to the equilibrium solution for large
$Re$
. The corresponding
$N$
factors shown in figure 13(b) reveal that the equilibrium model lies closest to the non-equilibrium one, but there is a noticeable discrepancy between them.
Figure 13. (a) Stability diagrams for
$Ma=5$
,
$T_{e}^{\ast }=T_{v,e}^{\ast }=1500~\text{K}$
and
$P_{e}^{\ast }=20~\text{kPa}$
. Three different models of the mean flow: frozen, non-equilibrium and equilibrium. Disturbances are non-equilibrium. (b) The
$N$
factor diagram corresponding to (a). Only 2D waves with
${\it\beta}=0$
are included.
Figure 14. (a) Stability diagrams for
$Ma=5$
,
$T_{e}^{\ast }=T_{v,e}^{\ast }=1500~\text{K}$
and
$P_{e}^{\ast }=20~\text{kPa}$
. (a) Three different models of the disturbances: frozen, non-equilibrium and equilibrium. Mean flow is non-equilibrium. (b) The
$N$
factor diagram corresponding to (a). Only 2D waves with
${\it\beta}=0$
are included.
The influence of acoustic absorption on the disturbance growth rates can be investigated in a similar manner. Three different disturbance models are now applied to a single fully non-equilibrium base flow. The three disturbance models are frozen, thermal equilibrium and full thermal non-equilibrium. The growth rates are compared in figure 14(a) and the
$N$
factors in figure 14(b). In this case, the frozen disturbance model is nearly indistinguishable from the non-equilibrium one. This is consistent with the findings of Fujii & Hornung (Reference Fujii and Hornung2003), who showed that, for air at temperatures below approximately 3000 K, the frequency of maximum acoustic absorption is several orders of magnitude lower than the frequency (
${\sim}\text{MHz}$
) of second-mode disturbances in hypervelocity boundary layers.
Similar calculations to figures 13 and 14 have been carried out for a wide range of Mach numbers and wall temperature conditions relevant to shock tunnel experiments. In some cases it was found that a non-equilibrium base flow was needed to fully capture the stability characteristics, but the frozen disturbance model never differed from the non-equilibrium one in any significant way. This finding is independent of the mean pressure of the flow, and results given in terms of the Reynolds number, dimensionless wavenumber and dimensionless frequency apply for any pressure. This is because both the viscous time scale and the time scale of vibrational relaxation are inversely proportional to the pressure, meaning that the effects of vibrational non-equilibrium are always the same at a given Reynolds number. The pressure independence of the results can also be deduced from the dimensionless governing equations (2.10), in which none of the parameters depend on the mean pressure.
3.6. Stability trends
This section summarizes the effects of wall temperature and Mach number on the stability characteristics. Although other researchers have conducted similar investigations, these results are included because they provide a complete picture of the instability behaviour over a much wider range of conditions than is considered in other references. Also, previous work (Mack Reference Mack1969, Reference Mack1984; Malik Reference Malik, Arnal and Michel1989b
; Masad et al.
Reference Masad, Nayfeh and Al-Maaitah1992) has targeted low-enthalpy conditions, whereas these results are relevant to high-enthalpy flows. The results given here are for a flat plate, but the growth rates and frequencies can be related to those for a sharp cone through the Mangler transform (Malik & Spall Reference Malik and Spall1991), by which it is found that the boundary layer is a factor of
$\sqrt{3}$
thinner for a cone than for a flat plate. Therefore, the results presented here can be related to those for a sharp cone by the transformations
$Re\rightarrow Re\sqrt{3}$
,
${\it\alpha}\rightarrow {\it\alpha}\sqrt{3}$
and
${\it\omega}\rightarrow {\it\omega}\sqrt{3}$
.
A sequence of simulations was conducted in which the wall temperature was fixed at 300 K while the free-stream temperature was incrementally raised from 70 to 2000 K. Although wall temperature effects are often spoken of in terms of ‘wall cooling’, we choose here to keep the wall temperature fixed and adjust the free-stream temperature so that comparisons with experiments in both high- and low-enthalpy impulse facilities can be made more easily. However, it should be recognized that, in terms of the stability behaviour, raising the free-stream temperature is nearly equivalent to cooling the wall, with differences being caused only by the temperature dependence of transport properties and specific heats, as well as non-equilibrium effects. Accordingly, we refer to cooling of the wall and heating of the free stream interchangeably. These results are also presented in terms of
$T_{w}^{\ast }/T_{e}^{\ast }$
rather than the customary
$T_{w}^{\ast }/T_{ad}^{\ast }$
because the adiabatic wall temperature depends on one’s assumptions about specific heats, chemistry and transport models, and must in general be simulated. Thus it can be difficult to reproduce results that are reported in terms of
$T_{ad}^{\ast }$
.
Figure 15(a) shows the maximum spatial growth rate as a function of Mach number for several values of the ratio
$T_{w}^{\ast }/T_{e}^{\ast }$
. Two sets of curves are shown: the dashed lines correspond to the first-mode instability and the solid lines to the second mode. For each curve, the spatial growth rate has been maximized over all frequencies
${\it\omega}$
and all spanwise wavenumbers
${\it\beta}$
. These curves are very similar to those of Mack (Reference Mack1984) and Masad et al. (Reference Masad, Nayfeh and Al-Maaitah1992), except that we consider a wider range of wall temperature conditions and hold
$T_{w}^{\ast }/T_{e}^{\ast }$
fixed along each curve rather than
$T_{w}^{\ast }/T_{ad}^{\ast }$
. Note that, in this figure and the ones that follow, the neglect of chemistry may not be valid for some of the high-Mach-number cases, but these results are still included for completeness.
Figure 15. Effect of wall cooling on the maximum spatial growth rates of the first (dashed lines) and second (solid lines) modes for
$Re=1500$
,
$T_{w}^{\ast }=300~\text{K}$
,
$T_{e}^{\ast }=70{-}2000~\text{K}$
,
$P_{e}^{\ast }=10~\text{kPa}$
. (a) Maximum growth rates for first and second modes. (b) Maximum growth rate of second mode normalized by
${\it\delta}_{99}$
.
For the first mode, no instability was found for some of the cases having
$T_{w}^{\ast }/T_{e}^{\ast }<1$
, so these lines are not included on the figure. For the second mode, a monotonic increase in growth rate is observed as the wall is cooled. Additionally, the second mode becomes unstable at lower Mach numbers as the wall is cooled, with instability occurring at Mach numbers as low as 2.5. Although it is commonly stated that the second mode becomes dominant above approximately
$Ma=4$
, this statement is true only when the wall is nearly adiabatic. This result is consistent with that of Mack (Reference Mack, Ashpis, Gatski and Hirsch1993), who also observed second-mode instability at
$Ma=3$
.
The destabilizing influence of wall cooling appears to be caused mainly by thinning of the boundary layer. This is illustrated in figure 15(b) in which the spatial growth rate of the second mode is normalized by
${\it\delta}_{99}$
, the boundary layer thickness at which the streamwise velocity reaches 99 % of the free-stream velocity, rather than the Blasius thickness
${\it\delta}=\sqrt{{\it\nu}_{e}^{\ast }x^{\ast }/U_{e}^{\ast }}$
. The maximum second-mode growth rates normalized in this way experience much less variation with both Mach number and wall temperature than in figure 15(a), indicating that the spatial growth rate is nearly inversely proportional to
${\it\delta}_{99}$
. Although it is well known that the frequency and wavenumber of the second mode scale inversely with the boundary layer thickness (Demetriades Reference Demetriades1977; Stetson et al.
Reference Stetson, Thompson, Donaldson and Siler1983), the correlation of growth rate with boundary layer thickness has been less widely recognized.
The relationship between the frequency of second-mode disturbances and the boundary layer thickness was first examined experimentally by Demetriades (Reference Demetriades1977) and Stetson et al. (Reference Stetson, Thompson, Donaldson and Siler1983), who both made hot-wire measurements of disturbance fluctuations in the same Mach 8 wind tunnel. Their results showed that the non-dimensional frequency
$2f^{\ast }{\it\delta}_{99}/U_{e}^{\ast }$
is nearly constant, decreasing only slightly with Reynolds number and with increased wall cooling. In this relation
$f^{\ast }$
is the dimensional frequency in hertz. Figure 16 shows the variation of dimensionless frequency with Mach number and wall temperature ratio at
$Re=1000$
. For figure 16(a), the plotted frequency is the most amplified one, i.e. the frequency for which the
$N$
factor is largest; this is the frequency that is measured most often in experiments. Figure 16(b) shows the most unstable frequency, that is, the frequency that has the largest spatial growth rate
$-{\it\alpha}_{i}$
. For flows with a cold wall, the most unstable and most amplified frequencies follow similar trends, and the most unstable frequency is slightly smaller. This is consistent with the results of Marineau et al. (Reference Marineau, Moraru, Lewis, Norris, Lafferty, Wagnild and Smith2014), who measured both the most unstable and the most amplified frequencies for flow over a sharp cone. However, for an adiabatic wall, the most unstable and most amplified frequencies exhibit very different variations with Mach number.
Figure 16. Variation of second-mode frequencies with Mach number and wall temperature at
$Re=1000$
. For
$T_{w}^{\ast }=T_{ad}^{\ast }$
, the edge temperature is
$T_{e}^{\ast }=70~\text{K}$
; for all other cases,
$T_{w}^{\ast }=300~\text{K}$
. (a) Most amplified frequency (max
$N$
factor). (b) Most unstable frequency (max
$-{\it\alpha}_{i}$
).
For flows with highly cooled walls, figures 15(b) and 16 reveal that the effect of the wall temperature can be nearly scaled out of the problem by normalizing both the frequency and the growth rate by
${\it\delta}_{99}$
. This length scale is preferable to the Blasius boundary layer thickness
${\it\delta}=\sqrt{{\it\nu}_{e}^{\ast }x^{\ast }/U_{e}^{\ast }}$
, which takes into account only the free-stream conditions, whereas the
${\it\delta}_{99}$
length scale includes the physical structure of the boundary layer and encapsulates the change in thickness that occurs when the wall is cooled. These results suggest that by scaling the growth rate and frequency with
${\it\delta}_{99}$
one can generate a single dispersion curve that, to a good approximation, describes a range of wall temperature conditions. An example of such a dispersion curve is given in figure 17. In the left panel the growth rate and frequency are scaled by the usual Blasius boundary layer thickness
${\it\delta}$
, whereas in the right one they are scaled by
${\it\delta}_{99}$
. The latter scaling nearly collapses the dispersion curves onto one another, though for the case with large
$T_{w}/T_{e}$
the collapse is not as good. For the cases with smaller
$T_{w}/T_{e}$
, however, the collapse is excellent, especially along the lower neutral branch, which contributes most significantly to the
$N$
factor.
Figure 17. Spatial growth rates for several different wall temperature conditions:
$Ma=4.5$
,
$T_{w}^{\ast }=300~\text{K}$
,
$Re=1500$
,
$P_{e}^{\ast }=10~\text{kPa}$
. Only 2D waves (
${\it\beta}=0$
) are included. (a) Normalized using Blasius thickness
${\it\delta}=\sqrt{{\it\nu}_{e}^{\ast }x^{\ast }/U_{e}^{\ast }}$
. (b) Normalized using
${\it\delta}_{99}$
boundary layer thickness.
The collapse of the scaled dispersion curves in figure 17(b) suggests that one might be able to collapse the
$N$
factors in a similar fashion. To achieve this, it is assumed that for a given Mach number the dispersion curves for all wall temperatures are well approximated by a two-parameter function
$\mathscr{F}$
of the form
$$\begin{eqnarray}-{\it\alpha}_{i}^{\ast }{\it\delta}_{99}=\mathscr{F}\left(Re,\frac{{\it\omega}^{\ast }{\it\delta}_{99}}{U_{e}^{\ast }}\right).\end{eqnarray}$$
This functional dependence follows from the collapsed dispersion curves in figure 17(b) and the fact that the dispersion curve is also a function of the Reynolds number. One can show that the following definitions of the
$N$
factor are equivalent:
$$\begin{eqnarray}N(F)=\int _{x_{o}^{\ast }}^{x^{\ast }}-{\it\alpha}_{i}^{\ast }(F,x^{\ast })\,\text{d}x^{\ast }=\int _{Re_{o}}^{Re}-2{\it\alpha}_{i}(F,Re)\,\text{d}Re.\end{eqnarray}$$
Here
$Re_{o}=\sqrt{U_{e}^{\ast }x_{o}^{\ast }/{\it\nu}_{e}^{\ast }}$
is the Reynolds number at which the frequency
$F$
first becomes unstable,
${\it\alpha}_{i}^{\ast }$
is the dimensional spatial growth rate and
${\it\alpha}_{i}$
is the non-dimensional growth rate normalized by the Blasius thickness
${\it\delta}$
. In (3.5),
$F$
is the frequency parameter defined by
$$\begin{eqnarray}F=\frac{{\it\omega}^{\ast }{\it\nu}_{e}^{\ast }}{U_{e}^{\ast 2}}.\end{eqnarray}$$
Making use of these definitions, one can rewrite the dispersion relation (3.4) in the form
$$\begin{eqnarray}{\it\alpha}_{i}\frac{{\it\delta}_{99}}{{\it\delta}}=\mathscr{F}\left(Re,F\frac{{\it\delta}_{99}}{{\it\delta}}\right).\end{eqnarray}$$
By substituting this result into (3.5), one arrives at the functional dependence:
$$\begin{eqnarray}N{\it\delta}_{99}/{\it\delta}=\mathscr{F}(Re,F{\it\delta}_{99}/{\it\delta}).\end{eqnarray}$$
The ratio
${\it\delta}_{99}/{\it\delta}$
is constant for a given Mach number and usually falls between approximately 3 and 20, increasing with both Mach number and wall temperature. The functional form of (3.8) suggests that, if one makes measurements of maximum
$N$
factors
$N_{1}$
and
$N_{2}$
for two different wall temperatures, then the maximum
$N$
factors and the corresponding most amplified frequencies
$F_{1}$
and
$F_{2}$
are related by
$$\begin{eqnarray}N_{2}=N_{1}\frac{{\it\delta}_{99,1}/{\it\delta}_{1}}{{\it\delta}_{99,2}/{\it\delta}_{2}},\quad F_{2}=F_{1}\frac{{\it\delta}_{99,1}/{\it\delta}_{1}}{{\it\delta}_{99,2}/{\it\delta}_{2}}.\end{eqnarray}$$
The effectiveness of this scaling is demonstrated by several examples given in figure 18. Figure 18(a) shows the maximum
$N$
factors for several different wall temperature conditions, as well as individual
$N$
factor curves for the frequency
$F=4\times 10^{-4}$
. The
$N$
factor at
$Re=1500$
is approximately 2.6 times greater for the cold wall than for
$T_{w}/T_{e}=3$
, and the frequency
$F=4\times 10^{-4}$
is amplified at different locations for each wall temperature condition.
Figure 18. Maximum
$N$
factors for several different wall temperature conditions:
$Ma=4.5$
,
$T_{w}^{\ast }=300~\text{K}$
,
$Re=1500$
,
$P_{e}^{\ast }=10~\text{kPa}$
. Only 2D waves (
${\it\beta}=0$
) are included. (a) Maximum
$N$
factors and individual
$N$
factor curves for
$F=4\times 10^{-4}$
. (b) Maximum
$N$
factors scaled by
${\it\delta}_{99}/{\it\delta}$
and individual
$N$
factor curves corresponding to
$F({\it\delta}_{99}/{\it\delta})=2\times 10^{-3}$
.
Figure 18(b) shows the same
$N$
factor curves, but scaled by the ratio
${\it\delta}_{99}/{\it\delta}$
, which significantly reduces the spread between the curves. Also shown in figure 18(b) are the
$N$
factors for the frequencies satisfying
$F({\it\delta}_{99}/{\it\delta})=2\times 10^{-3}$
. All of the frequencies scaled in this way are amplified at nearly the same Reynolds number. For both the frequencies and the
$N$
factors, the collapse is poorer for large
$T_{w}/T_{e}$
. The poorer collapse for large
$T_{w}/T_{e}$
may in part be associated with the much smaller free-stream temperatures used in these cases, since the viscosity varies more rapidly at low temperatures.
This manner of collapsing the
$N$
factors has been applied over a wide range of wall temperatures and Mach numbers. The result is shown in figure 19. Figure 19(a) reports the maximum
$N$
factors at
$Re=1500$
without using the scaling of (3.8). This plot includes results only for 2D waves (
${\it\beta}=0$
) since these are most amplified for the second mode, but it should be noted that 3D first-mode waves can produce larger
$N$
factors for some of the cases with lower Mach numbers and large values of
$T_{w}/T_{e}$
. The maximum
$N$
factors shown in the figure exhibit the same trends as the growth rates from figure 15, namely, maximum amplification at approximately
$Ma=5$
, a systematic increase in growth as the wall is cooled and reduction in the most amplified Mach number as the wall is cooled. This indicates that, although the influences of Mach number and wall temperature have historically been characterized mainly using the maximum spatial growth rates, the
$N$
factors behave similarly.
Figure 19. Effect of Mach number and wall cooling on the maximum
$N$
factors for the second mode. Only 2D waves (
${\it\beta}=0$
) are considered. In all cases,
$Re=1500$
,
$T_{w}^{\ast }=300~\text{K}$
and
$P_{e}^{\ast }=10~\text{kPa}$
. (a) Maximum
$N$
factors. (b) Maximum
$N$
factors scaled by
${\it\delta}_{99}/{\it\delta}$
.
Figure 19(b) attempts to collapse the maximum
$N$
factors for different wall temperatures by scaling the
$N$
factor by
${\it\delta}_{99}/{\it\delta}$
. This scaling is quite successful for the cases with
$0.2\leqslant T_{w}/T_{e}\leqslant 1.5$
, but performs poorly for the cases with large
$T_{w}/T_{e}$
. In spite of the imperfect collapse of the data, the scaling arguments demonstrated here may be useful for estimating how a change in the wall temperature or stagnation temperature might affect the stability characteristics and transition location in an experiment.
The frequency scaling presented in figure 16 and its variation with wall temperature and Mach number are compared with experimental measurements in table 1. For both experiments and computations, the frequency given is the most amplified one; in the computations only 2D waves with
${\it\beta}=0$
are included since the measured frequencies are believed to be second-mode waves. Measurements in the AEDC Tunnel 9 were made by Marineau et al. (Reference Marineau, Moraru, Lewis, Norris, Lafferty, Wagnild and Smith2014) using PCB pressure transducer arrays. Measurements in Caltech’s T5 reflected shock tunnel were made by Parziale (Reference Parziale2013) using focused laser differential interferometry (FLDI). The FLDI technique was also used in the Hypersonic Ludwieg tube at Braunschweig (HLB) (D. Heitmann, unpublished observations). Further details regarding the operation of this facility are given by Estorf, Wolf & Radespiel (Reference Estorf, Wolf and Radespiel2004) and Heitmann, Radespiel & Knauss (Reference Heitmann, Radespiel and Knauss2011). Measurements in the High Enthalpy shock tunnel at Göttingen (HEG) were made by Laurence et al. (Reference Laurence, Wagner, Ozawa, Schramm and Hannemann2014) using schlieren deflectometry and full-field pulsed-laser schlieren photography. Measurements in AEDC Tunnel B were made using hot wires by Stetson et al. (Reference Stetson, Thompson, Donaldson and Siler1983); these data were more recently re-plotted and analysed by Schneider (Reference Schneider2006).
Table 1. Comparison of most amplified frequencies from 2D computation (
${\it\beta}=0$
) and experiment. Here
$Ma$
is the post-shock Mach number; and
$Re$
is the Reynolds number for flow over a cone, equal to
$\sqrt{3}$
times the value for a plate. The streamwise distance is
$s$
and the nose radius is
$R_{n}$
. For AEDC Tunnel B, the wall temperature is nearly adiabatic. For all other cases, the wall temperature is ambient.
All of the measured frequencies agree reasonably well with the simulated ones, although the agreement is slightly poorer for the reflected shock tunnels T5 and HEG, where the free-stream conditions are known with less certainty. The agreement between experiment and computation is also generally better at low Reynolds numbers, where the measured disturbances remain in the linear regime. In all cases except HEG, the cones are nominally sharp, and the reported nose radii (which can be determined from table 1) are small enough that the measurement stations are far downstream of the entropy swallowing length. For the HEG data, the cone is slightly blunted with a nose radius of
$R_{n}=2.5~\text{mm}$
. According to the swallowing length correlations reported by Stetson (Reference Stetson1980), the measurement station is at about twice the swallowing length, but the local boundary layer thickness may still be influenced by the nose bluntness (Stetson et al.
Reference Stetson, Thompson, Donaldson and Siler1984). This may in part explain why the measured frequencies are slightly larger than the sharp-cone predictions. This may also be a factor in the T5 data, where the ratio of the streamwise distance
$s$
to the nose radius
$R_{n}$
is approximately
$s/R_{n}=1260$
and the measured frequencies are again somewhat greater than the sharp-cone prediction. The experiments tabulated here employ an array of different measurement techniques and span a wide range of flow conditions, from low-enthalpy adiabatic flow (AEDC Tunnel B) to high-enthalpy reflected shock tunnels (HEG and T5). Although the Reynolds numbers of the measurements are slightly different from the value
$Re=1000$
used in figure 16, the measured data confirm the major features of the trends shown in the figure.
A final illustration of the effects of Mach number and wall temperature on the stability behaviour is given in figure 20. This plot contains the same data from figure 15, but is organized in contour form to better portray the stability boundaries. The contour levels indicate the spatial growth rate, which is optimized over all values of the frequency
${\it\omega}$
and spanwise wavenumber
${\it\beta}$
. For reference, a line (symbols) corresponding to the adiabatic wall condition is also included. Most practical conditions would fall below this line, as the points above it correspond to a heated wall. There is a clear region at low Mach numbers and small values of
$T_{w}^{\ast }/T_{e}^{\ast }$
for which no modal instabilities are found. Although this region shrinks slightly as the Reynolds number is increased, it is still present even at rather high Reynolds numbers,
$Re>4000$
. Since no modal instabilities are found at these conditions, one might conclude that large-amplitude (nonlinear) disturbances are needed to cause transition when the Mach number is low and the wall is cold. However, recent transient growth calculations have reported elevated levels of non-modal amplification at these conditions (Bitter & Shepherd Reference Bitter and Shepherd2014), which offers an alternative transition mechanism for infinitesimal disturbances.
Figure 20. Map of unstable regions for first and second modes:
$Re=1500$
,
$T_{w}^{\ast }=300~\text{K}$
,
$T_{e}^{\ast }=70{-}2000~\text{K}$
,
$P_{e}^{\ast }=10~\text{kPa}$
. Contour levels are spatial growth rates maximized over all values of
${\it\omega}$
and
${\it\beta}$
for the first mode (dashed) and second mode (solid).
4. Conclusions
The influence of very high levels of wall cooling,
$T_{w}^{\ast }/T_{e}^{\ast }\ll 1$
, on the stability of hypervelocity boundary layers is studied using the locally parallel, linear stability framework. The appearance of supersonic unstable modes at very high levels of wall cooling is investigated. The supersonic unstable modes are most prominent at approximately
$Ma=5$
and are enhanced at high Reynolds numbers. These modes cause the flow to be unstable over a much wider band of frequencies, especially as the Reynolds number is increased. Although the supersonic unstable modes significantly alter the behaviour of the dispersion curve, they do not influence the integrated
$N$
factors until the Reynolds number is quite high. A synchronism between free-stream acoustic waves and the supersonic modes was identified, which has also been observed in boundary layers with gas injection (Fedorov et al.
Reference Fedorov, Soudakov and Leyva2014). Since this synchronism is located downstream of the upper neutral branch, it is unlikely to affect the amplification of waves, but it may still influence the nonlinear stages of their downstream development.
The importance of modelling vibrational non-equilibrium for flows of air over slender bodies is briefly examined for flow conditions relevant to impulse facilities. It is found that the temperature profiles of the base flow are significantly influenced by vibrational energy transfer, which indirectly affects the stability characteristics, and therefore a non-equilibrium base flow calculation is generally needed. In contrast, a vibrationally frozen stability analysis is nearly indistinguishable from a non-equilibrium one, since for air the frequency of second-mode disturbances is too high for significant vibrational energy transfer to take place.
The effects of wall temperature and Mach number on the boundary layer stability characteristics are summarized over a wide range of conditions. Unlike previous investigations of this sort, the current results focus on high-enthalpy flows and conditions that are relevant to impulse facilities. For highly cooled walls, the second mode is found to be unstable for Mach numbers as low as 2.5. The tuning of second-mode frequencies to the boundary layer thickness and its dependence on wall temperature and Mach number are explored and compared with experimental data. Spatial growth rates are also found to scale with the boundary layer thickness
${\it\delta}_{99}$
, and this scaling is extended to show that the
$N$
factor and most amplified frequency follow a simple scaling with
${\it\delta}_{99}$
as well.
Acknowledgements
The authors are grateful to A. Fedorov and A. Tumin for key insights regarding the interpretation of the supersonic modes and helpful discussions of other aspects of the paper. The authors are also grateful to H. Johnson and G. Candler for generously providing use of their STABL stability software package, which was essential for validation at high-enthalpy conditions.
Appendix A
In this appendix the non-zero elements of the 10 × 10 matrix
$\unicode[STIX]{x1D63C}$
from (2.18) are given. The notation of Malik (Reference Malik1990) is followed as closely as possible. For convenience, the following three parameters from that paper are used:
$$\begin{eqnarray}{\it\xi}\equiv {\it\alpha}\bar{U}+{\it\beta}\bar{W}+{\it\omega},\quad {\it\chi}\equiv \left[\frac{Re}{\bar{{\it\mu}}}+\text{i}{\it\gamma}_{e}Ma^{2}{\it\xi}(2+r)\right]^{-1},\quad r\equiv \frac{{\it\lambda}}{{\it\mu}},\end{eqnarray}$$
where
${\it\lambda}$
is the second coefficient of viscosity, which is taken to be
${\it\lambda}=-2{\it\mu}/3$
in accordance with Stokes’ hypothesis of zero bulk viscosity. It should be noted, however, that some of the effects of non-zero bulk viscosity are included through the explicit modelling of vibrational non-equilibrium. With the above definitions in place, the non-zero coefficients of the matrix
$\unicode[STIX]{x1D63C}$
are the following:
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D608}_{12}=1, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{21}=\text{i}{\it\xi}\frac{Re}{\bar{{\it\mu}}\bar{T}}+{\it\alpha}^{2}+{\it\beta}^{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{22}=-\frac{1}{\bar{{\it\mu}}}\frac{\partial {\it\mu}}{\partial T}\bar{T}^{\prime }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{23}=-\text{i}{\it\alpha}\frac{1}{\bar{{\it\mu}}}\frac{\partial {\it\mu}}{\partial T}\bar{T}^{\prime }+\frac{Re}{\bar{{\it\mu}}\bar{T}}\bar{U}^{\prime }-\text{i}{\it\alpha}(1+r)\frac{\bar{T}^{\prime }}{\bar{T}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{24}=\text{i}{\it\alpha}\frac{Re}{\bar{{\it\mu}}}-{\it\alpha}{\it\xi}(1+r){\it\gamma}_{e}Ma^{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{25}=-\frac{1}{\bar{{\it\mu}}}\frac{\partial {\it\mu}}{\partial T}\bar{U}^{\prime \prime }-\frac{1}{\bar{{\it\mu}}}\frac{\partial ^{2}{\it\mu}}{\partial T^{2}}\bar{U}^{\prime }\bar{T}^{\prime }+(1+r)\frac{{\it\alpha}{\it\xi}}{\bar{T}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{26}=-\frac{1}{\bar{{\it\mu}}}\frac{\partial {\it\mu}}{\partial T}\bar{U}^{\prime }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D608}_{31}=-\text{i}{\it\alpha}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{33}=\frac{\bar{T}^{\prime }}{\bar{T}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D608}_{34}=-\text{i}{\it\xi}{\it\gamma}_{e}Ma^{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{35}=\frac{\text{i}{\it\xi}}{\bar{T}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D608}_{37}=-\text{i}{\it\beta}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{41}=-\text{i}{\it\alpha}{\it\chi}\left[(2+r)\frac{\bar{T}^{\prime }}{\bar{T}}+\frac{2}{\bar{{\it\mu}}}\frac{\partial {\it\mu}}{\partial T}\bar{T}^{\prime }\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D608}_{42}=-\text{i}{\it\alpha}{\it\chi}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{43}=-{\it\chi}\left[{\it\alpha}^{2}+{\it\beta}^{2}+\text{i}{\it\xi}\frac{Re}{\bar{{\it\mu}}\bar{T}}-(2+r)\left(\frac{\bar{T}^{\prime \prime }}{\bar{T}}+\frac{\bar{T}^{\prime 2}}{\bar{T}}\frac{1}{\bar{{\it\mu}}}\frac{\partial {\it\mu}}{\partial T}\right)\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{44}=-\text{i}{\it\chi}(2+r){\it\gamma}_{e}Ma^{2}\left[({\it\alpha}\bar{U}^{\prime }+{\it\beta}\bar{W}^{\prime })+{\it\xi}\left(\frac{\bar{T}^{\prime }}{\bar{T}}+\frac{1}{\bar{{\it\mu}}}\frac{\partial {\it\mu}}{\partial T}\bar{T}^{\prime }\right)\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{45}=\text{i}{\it\chi}\left[({\it\alpha}\bar{U}^{\prime }+{\it\beta}\bar{W}^{\prime })\left(\frac{1}{\bar{{\it\mu}}}\frac{\partial {\it\mu}}{\partial T}+\frac{2+r}{\bar{T}}\right)+\frac{{\it\xi}}{\bar{T}}(2+r)\frac{1}{\bar{{\it\mu}}}\frac{\partial {\it\mu}}{\partial T}\bar{T}^{\prime }\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{46}=\frac{\text{i}{\it\xi}{\it\chi}(2+r)}{\bar{T}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{47}=-\text{i}{\it\beta}{\it\chi}\left[(2+r)\frac{\bar{T}^{\prime }}{\bar{T}}+\frac{2}{\bar{{\it\mu}}}\frac{\partial {\it\mu}}{\partial T}\bar{T}^{\prime }\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D608}_{48}=-\text{i}{\it\beta}{\it\chi}, & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D608}_{56}=1, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{62}=-\frac{Ma^{2}({\it\gamma}_{e}-1)\bar{{\it\mu}}{\it\sigma}_{e}}{\bar{k}}2\bar{U}^{\prime }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{63}=-2\text{i}\frac{Ma^{2}({\it\gamma}_{e}-1)\bar{{\it\mu}}{\it\sigma}_{e}}{\bar{k}}({\it\alpha}\bar{U}^{\prime }+{\it\beta}\bar{W}^{\prime })+\frac{Re{\it\sigma}_{e}}{\bar{T}\bar{k}}\bar{T}^{\prime }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{64}=-\text{i}{\it\xi}({\it\gamma}_{e}-1)Ma^{2}\frac{Re{\it\sigma}_{e}}{\bar{k}}+\underline{\frac{Re^{2}{\it\sigma}_{e}}{\bar{k}}\frac{\partial Q}{\partial p}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{65}={\it\alpha}^{2}+{\it\beta}^{2}-\frac{\bar{k}^{\prime \prime }}{\bar{k}}-\frac{Ma^{2}({\it\gamma}_{e}-1){\it\sigma}_{e}}{\bar{k}}\frac{\partial {\it\mu}}{\partial T}(\bar{U}^{\prime 2}+\bar{W}^{\prime 2})+\text{i}{\it\xi}\frac{Re{\it\sigma}_{e}}{\bar{k}\bar{T}}+\underline{\frac{Re^{2}{\it\sigma}_{e}}{\bar{k}}\frac{\partial Q}{\partial T}},\quad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{66}=-2\frac{\bar{k}^{\prime }}{\bar{k}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{68}=-\frac{Ma^{2}({\it\gamma}_{e}-1)\bar{{\it\mu}}{\it\sigma}_{e}}{\bar{k}}2\bar{W}^{\prime }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{69}=\underline{\frac{Re^{2}{\it\sigma}_{e}}{\bar{k}}\frac{\partial Q}{\partial T_{v}}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D608}_{78}=1, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{83}=\frac{Re\bar{W}^{\prime }}{\bar{{\it\mu}}\bar{T}}-\text{i}{\it\beta}\frac{1}{\bar{{\it\mu}}}\frac{\partial {\it\mu}}{\partial T}\bar{T}^{\prime }-\text{i}{\it\beta}(1+r)\frac{\bar{T}^{\prime }}{\bar{T}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{84}=\text{i}{\it\beta}\frac{Re}{\bar{{\it\mu}}}-(1+r){\it\beta}{\it\xi}{\it\gamma}_{e}Ma^{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{85}=\frac{{\it\beta}{\it\xi}(1+r)}{\bar{T}}-\frac{1}{\bar{{\it\mu}}}\frac{\partial {\it\mu}}{\partial T}\bar{W}^{\prime \prime }-\frac{1}{\bar{{\it\mu}}}\frac{\partial ^{2}{\it\mu}}{\partial T^{2}}\bar{T}^{\prime }\bar{W}^{\prime }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{86}=-\frac{1}{\bar{{\it\mu}}}\frac{\partial {\it\mu}}{\partial T}\bar{W}^{\prime }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{87}=\text{i}{\it\xi}\frac{Re}{\bar{{\it\mu}}\bar{T}}+{\it\alpha}^{2}+{\it\beta}^{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{88}=-\frac{1}{\bar{{\it\mu}}}\frac{\partial {\it\mu}}{\partial T}\bar{T}^{\prime }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D608}_{9,10}=\underline{1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{10,3}=\underline{\frac{Re{\it\sigma}_{e}}{\bar{k}_{v}\bar{T}}\bar{c}_{v,v}\bar{T}_{v}^{\prime }}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{10,4}=-\underline{\frac{Re^{2}{\it\sigma}_{e}}{\bar{k}_{v}}\frac{\partial Q}{\partial p}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{10,5}=\underline{-\frac{Re^{2}{\it\sigma}_{e}}{\bar{k}_{v}}\frac{\partial Q}{\partial T}-\frac{\bar{T}_{v}^{\prime \prime }}{\bar{k}_{v}}\frac{\partial k_{v}}{\partial T}-\frac{\bar{T}_{v}^{\prime }}{\bar{k}_{v}}\left(\frac{\partial k_{v}}{\partial T}\right)^{\prime }}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{10,6}=\underline{-\frac{\bar{T}_{v}^{\prime }}{\bar{k}_{v}}\frac{\partial k_{v}}{\partial T}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{10,9}=\underline{{\it\alpha}^{2}+{\it\beta}^{2}+i{\it\xi}\frac{Re{\it\sigma}_{e}}{\bar{k}_{v}\bar{T}}\bar{c}_{v,v}-\frac{Re^{2}{\it\sigma}_{e}}{\bar{k}_{v}}\frac{\partial Q}{\partial T_{v}}-\frac{\bar{T}_{v}^{\prime \prime }}{\bar{k}_{v}}\frac{\partial k_{v}}{\partial T_{v}}-\frac{\bar{T}_{v}^{\prime }}{\bar{k}_{v}}\left(\frac{\partial k_{v}}{\partial T_{v}}\right)^{\prime }}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D608}_{10,10}=\underline{-\frac{\bar{T}_{v}^{\prime }}{\bar{k}_{v}}\frac{\partial k_{v}}{\partial T_{v}}-\frac{k_{v}^{\prime }}{\bar{k}_{v}}}. & \displaystyle\end{eqnarray}$$
