2.1. Pridmore-Brown equation

To derive the underlying inviscid linear stability equation, the PBE, we start with the compressible Euler equations, neglecting viscosity and heat conduction (Spurk & Aksel Reference Spurk and Aksel2019). Accordingly, the model problem considered in this work is intended to shed light on the inviscid acoustic perturbations in the high-Reynolds-number limit (Rienstra & Hirschberg Reference Rienstra and Hirschberg2001; Delfs Reference Delfs2016). These assumptions hold especially for the high-frequency limit, as shown by Aurégan et al. (Reference Aurégan, Starobinski and Pagneux2001) and Brambley (Reference Brambley2011), which analyse the effects of dissipation and shear on the BC at the wall.

The equations are linearised around a mean state, indexed by $0$ , giving the linearised Euler equations:

(2.1a) \begin{align} \frac {\partial \rho '}{\partial t}+\boldsymbol {v}_0\boldsymbol \cdot \boldsymbol \nabla \rho '+\rho _0\boldsymbol \nabla \boldsymbol \cdot \boldsymbol {v}'+\boldsymbol {v}'\boldsymbol \cdot \boldsymbol \nabla \rho _0+\rho '\boldsymbol \nabla \boldsymbol \cdot \boldsymbol {v}_0 &= 0 ,\end{align}

(2.1b) \begin{align} \rho _0\left (\frac {\partial \boldsymbol v'}{\partial t}+\boldsymbol {v}_0\boldsymbol \cdot \boldsymbol \nabla \boldsymbol v'\right )+\boldsymbol \nabla p'+\rho _0\boldsymbol v'\boldsymbol \cdot \boldsymbol \nabla \boldsymbol {v}_0+\rho '\boldsymbol {v}_0\boldsymbol \cdot \boldsymbol \nabla \boldsymbol {v}_0 &= 0 ,\end{align}

(2.1c) \begin{align} p'&= c^2\rho ',\end{align}

where $\boldsymbol {v}$ denotes the velocity vector, $\rho$ the density and $p$ the pressure. The small unsteady perturbations are denoted by a prime. Equation (2.1c ) represents the linearised thermodynamic equation of state $p=p(\rho , s)$ under the assumption of constant entropy, with $c=\sqrt { (\partial p/\partial \rho )_{s}}$ describing the speed of sound.

For the boundary-layer flow, we assume a parallel mean flow of the form $\boldsymbol {v}_0=U_0(y)\cdot \boldsymbol {e}_x$ , where $y$ describes the wall-normal direction. In addition, we assume the mean flow as isothermal, implying a constant mean density $\rho _0$ .

It is noted that the isothermal assumption is a simplification, especially for the considered Mach number range, as thermal effects will influence the results quantitatively. Eigenmodes are particularly influenced by the character and position of the critical layer, which exists for both an isothermal and a non-isothermal boundary layer. Since the position of the critical layer is slightly moved by a non-constant mean temperature profile, a quantitative modification of the isothermal results shown in this paper is to be expected. However, the qualitative observations made in this work will also be retained in the non-isothermal case. The aim of this work is therefore to reveal fundamental effects employing a model problem.

The resulting linearised Euler equations for the parallel shear flow are non-dimensionalised using $\rho _0$ , the far-field velocity $U_{\infty }$ and the boundary-layer thickness $\delta$ . Due to homogeneity of the equations with regard to $x$ , $z$ and $t$ , a normal-mode approach is applied, which for 2-D perturbations reads

(2.2) \begin{equation} q'(x,y,t) =\hat {q}(y)e^{i(\alpha x-\omega t)} \qquad \text {with} \quad q\in \left \{u,v,\rho \right \}. \end{equation}

While in this section only 2-D disturbances are considered, the explicit $z$ dependence of disturbances through the wavenumber $\beta$ is also taken into account in § 3. In the approach (2.2), $\hat {q}(y)$ denotes the complex amplitudes of the perturbations. For temporal stability analysis, we assume $\omega \in \mathbb {C}$ with its real part $\omega _r$ representing the temporal frequency and its imaginary part $\omega _i$ the temporal growth rate. Parameter $\alpha \in \mathbb {R}$ is the non-dimensionalised streamwise wavenumber. Applying (2.2) to the non-dimensionalised linearised Euler equations yields a coupled system of first-order ordinary differential equations (ODEs) for the perturbation amplitudes, which can be converted to a single ODE of second order for the density perturbation amplitude $\hat {\rho }(y)$ :

(2.3) \begin{equation} \frac {{\rm d}^{2}\hat {\rho }}{{\rm d}y^{2}}+\frac {2\alpha }{\omega -U_0\alpha }\,\frac {{\rm d}U_0}{{\rm d}y}\frac {{\rm d}\hat {\rho }}{{\rm d}y}+\left [M^{2}(\omega -U_0\alpha )^{2}-\alpha ^{2}\right ]\hat {\rho }=0\, , \end{equation}

which is the so-called PBE. Here $M=U_{\infty }/c$ is the free-stream Mach number. The relations for the remaining perturbation amplitudes in terms of $\hat {\rho }(y)$ are given in Appendix A.

2.2. Boundary conditions

The eigenvalue problem is given by finding those values $(\alpha$ , $\omega )$ that ensure compatibility of the solution of the PBE (2.3) with the BCs of the problem. Since we aim to investigate instabilities originating from the boundary layer itself, we do not allow energy transfer from the far field to the inside of the boundary layer. Furthermore, an infinite perturbation amplitude in the far field would seem unphysical. This leads to the BC of bounded perturbation amplitude in the far field ( $y \rightarrow \infty$ ) (Mack Reference Mack1984). Also common in stability analysis is the more restrictive requirement for the amplitudes to vanish in the far field, giving

(2.4) \begin{equation} \lim _{y\to \infty }\hat {\rho }(y)=0. \end{equation}

We will see later in this section that for the temporal stability problem with $\omega _i\neq 0$ the condition of bounded far-field amplitudes reduces to (2.4) for the model equations in this paper, which means that both far-field BCs yield the same stability eigenvalues.

While we use the BC (2.4) here, it should be noted that other literature advocates a different BC. Brazier-Smith & Scott (Reference Brazier-Smith and Scott1984), Huerre & Monkewitz (Reference Huerre and Monkewitz1985), Crighton (Reference Crighton1989), Riedinger et al. (Reference Riedinger, Les Dizès and Meunier2010) and Brambley & Gabard (Reference Brambley and Gabard2014), among others, apply causality arguments that lead to unbounded behaviour in the far field. Riedinger et al. (Reference Riedinger, Les Dizès and Meunier2010, p. 258) comment in particular that ‘The difficulty with this condition is that it does not necessarily imply that the solution vanishes at infinity when the modes are neutral or damped ( $\mathrm {Im}\,\omega \lt 0$ )’, while Brambley & Gabard (Reference Brambley and Gabard2014, p. 5562) find ‘leaky surface waves excited by a line source that grow exponentially away from the surface’. However, for temporal instabilities with $\mathrm {Im}\,\omega \gt 0$ , there is a consensus that the behaviour at infinity should be bounded, in agreement with (2.4). Thus, the results presented in the following of this paper could be easily modified to use alternative BCs instead of (2.4).

The second BC for $\hat {\rho }$ is set up on the wall ( $y=0$ ). Since we consider an acoustic metasurface, we apply the acoustic impedance BC, which in the case of a straight wall without mean flow along the wall reads (Malmuth et al. Reference Malmuth, Fedorov, Shalaev, Cole, Hites, Williams and Khokhlov1998; Rienstra & Darau Reference Rienstra and Darau2011)

(2.5) \begin{equation} \hat {p}=-Z\,\hat {v} \quad \text {for }\, y=0 . \end{equation}

Equation (2.5) links the wall-normal velocity perturbation and the pressure perturbation at the wall via the complex surface impedance $Z$ , which depends on the wall and mean flow properties, such as the Mach number, as well as the perturbation frequency and wavenumber (Malmuth et al. Reference Malmuth, Fedorov, Shalaev, Cole, Hites, Williams and Khokhlov1998; Zhao et al. Reference Zhao, Wen, Zhou, Tu and Lei2022). Consequently, (2.5) describes the effect of the surface on the boundary-layer stability. The limit $Z\rightarrow \infty$ describes a rigid wall (Campos & Serrão Reference Campos and Serrão1998). In order to express (2.5) in terms of the density perturbation $\hat {\rho }$ , we substitute $\hat {v}$ by (A2) and $\hat {p}$ by $\hat p= ({1}/{M^2}) \hat {\rho }$ resulting from (2.1c ) by non-dimensionalising, which yields

(2.6) \begin{equation} \left.\frac {{\rm d}\hat {\rho }}{{\rm d}y}\right |_{y=0}=-i\omega \,Y\hat {\rho }(0) \, \end{equation}

with $Y=Z^{-1}$ denoting the acoustic admittance and thus $Y=0$ representing the rigid-wall case. The two BCs (2.4) and (2.6) together with the PBE (2.3) form the eigenvalue problem for $(\alpha , \omega )$ , which we solve for $\alpha \in \mathbb {R}$ and $\omega \in \mathbb {C}$ .

Using symmetry considerations given in Appendix B, it can be shown that in the case of rigid walls $Y=0$ , the eigenvalues $(\alpha ,\omega )$ occur both complexly conjugated and reflected, meaning that, together with an eigenvalue $(\alpha , \omega )$ , also $(\alpha , \omega ^*)$ and $(-\alpha , -\omega )$ occur as further solutions, having the opposite stability behaviour to $(\alpha , \omega )$ . Therefore, for rigid walls, there is coupling between stable and unstable solutions. Introducing a wall admittance $Y(\omega )$ , however, generally breaks these stable–unstable couplings, as the eigenvalues no longer occur in complex conjugate or reflected pairs. This is quite essential in view of the goal of using impedance walls for stabilisation. Exceptions are those impedance models that satisfy $Y^*(\omega )=-Y(\omega ^*)$ and both $Y_r(\omega )=0$ and $Y_i(-\omega )=-Y_i(\omega )$ for $\omega \in \mathbb {R}$ as well as those impedance models that satisfy $Y(-\omega )=-Y(\omega )$ and $Y_r(\omega ))=0$ for $\omega \in \mathbb {R}$ . Terms $Y_r$ and $Y_i$ denote the real and imaginary parts of the admittance. While for the first type of impedances the complex conjugation of the eigenvalues is preserved, the second type preserves the reflection of the eigenvalues, as derived in Appendix B. Note that reflection is also given for constant real $Y$ .

2.3. Eigenvalue problem for boundary-layer stability

In this section, we give the analytical solution of the PBE for the exponential boundary layer in terms of the CHF. This solution allows us to convert the differential eigenvalue problem consisting of the PBE (2.3) and the BCs (2.4) and (2.6) into a single algebraic eigenvalue equation. To model the boundary-layer flow, we assume an exponential shear profile, which, non-dimensionalised with $U_{\infty }$ and $\delta$ , reads

(2.7) \begin{align} U_0(y)=1-{\rm e}^{-y}. \end{align}

It should be emphasised that with the choice of non-dimensionalisation in this paper, $\delta$ does not correspond to the classical definition of the boundary-layer thickness $\delta _{\,99}$ , representing the height at which $99\, \%$ of the free-stream velocity is reached. Instead, at the wall-normal distance of $\delta$ , we have here $U_0 =1-{\rm e}^{-1}\approx 0.63$ .

Inserting this profile into the PBE (2.3) yields

(2.8) \begin{align} \frac {{\rm d}^{2}\hat {\rho }}{{\rm d}y^{2}}+\frac {2\alpha {\rm e}^{-y}}{\omega -\alpha \left (1-{\rm e}^{-y}\right )}\frac {{\rm d}\hat {\rho }}{{\rm d}y}+\left [M^{2}(\omega -\alpha \left (1-{\rm e}^{-y}\right ))^{2}-\alpha ^{2}\right ]\hat {\rho }=0. \end{align}

As derived by Zhang & Oberlack (Reference Zhang and Oberlack2021) and shown in Appendix C, this equation can be reduced to the confluent Heun differential equation (CHE). Therefore, the general solution for the PBE with exponential profile (2.8) can be given in terms of the CHF, reading

(2.9) \begin{align} \hat \rho (y) & =C_{1}\,{\rm e}^{iM\alpha {\rm e}^{-y}+\sqrt {\theta }y} \cdot \mathrm {HeunC}_1\left (q,\alpha _H,\gamma ,\delta _H,\epsilon ,\frac {\alpha }{\alpha -\omega }{\rm e}^{-y}\right )\nonumber \\ & \quad +C_{2}\,{\rm e}^{iM\alpha {\rm e}^{-y}-\sqrt {\theta }y} \cdot \mathrm {HeunC}_2\left (q,\alpha _H,\gamma ,\delta _H,\epsilon ,\frac {\alpha }{\alpha -\omega }{\rm e}^{-y}\right ). \end{align}

Here, $\mathrm {HeunC}_{1/2}(;z)$ denote the two linearly independent basic power series representations of the CHF, as described in more detail in Appendix C. The far-field exponent is given by

(2.10) \begin{align} \sqrt {\theta }=\sqrt {\alpha ^2-M^2\left (\omega -\alpha \right )^2}, \end{align}

describing both the exponential decay and the oscillatory behaviour of the solution in the far field, as is detailed in the next step.

With the exact PBE solution (2.9), we are now in the position to derive an algebraic eigenvalue equation for the eigenvalue problem of the PBE (2.8) with the BCs (2.4) and (2.6). To this end, the BCs are incorporated into the solution (2.9). We start with incorporating the far-field BC $\hat \rho (y\to \infty )=0$ . As $\mathrm {HeunC}(\ldots ,0)=1$ holds (Ronveaux & Arscott Reference Ronveaux and Arscott1995), the solution (2.9) simplifies in the far field to

(2.11) \begin{equation} \lim _{y\to \infty }\left (C_{1}{\rm e}^{\sqrt {\theta }y}+C_{2}{\rm e}^{-\sqrt {\theta }y}\right )=0 , \end{equation}

where the square root of the complex-valued $\theta$ is given by

(2.12) \begin{align} \sqrt {\theta }=\frac {\sqrt {2}}{2}\left (\sqrt {\sqrt {\theta _{r}^{2}+\theta _{i}^{2}}+\theta _{r}}+\operatorname {sign}\left (\theta _{i}\right ) i\sqrt {\sqrt {\theta _{r}^{2}+\theta _{i}^{2}}-\theta _{r}}\right ). \end{align}

Since we choose the positive solution branch of the square root function, (2.11) together with (2.12) requires $C_1=0$ for the temporal stability case $\omega _i \neq 0$ , since this solution branch always grows indefinitely in the far field, giving no physically plausible eigenfunction. Due to the exponential behaviour, $C_1=0$ also follows from the more general far-field BC of bounded amplitudes for $\omega _i\neq 0$ . The remaining solution branch

(2.13) \begin{align} \hat \rho (y) =C_{2}\,{\rm e}^{iM\alpha {\rm e}^{-y}-\sqrt {\theta }y} \cdot \mathrm {HeunC}_2\left (q,\alpha _H,\gamma ,\delta _H,\epsilon ,\frac {\alpha }{\alpha -\omega }{\rm e}^{-y}\right ) \end{align}

is incorporated into the wall BC (2.6). This yields a single determining equation for the eigenvalues $(\alpha , \omega )$ fulfilling both the far-field and wall BCs (2.4) and (2.6). This eigenvalue equation is given by

(2.14) \begin{equation} \left [ \left ( -iM\alpha -\sqrt {\theta } \right ) +i\omega Y\right ]\mathrm {HeunC_2}\left (\,;\frac {\alpha }{\alpha -{\omega }}\right )-\frac {\alpha }{\alpha -{\omega }}\mathrm {HeunC_2}'\left (\,;\frac {\alpha }{\alpha -{\omega }}\right )=0 \,, \end{equation}

solved in this paper for the temporal stability eigenvalues $\omega (\alpha ; M,Y)$ . Note that the prime denotes differentiation of $\mathrm {HeunC}_2(;z)$ with respect to its argument $z$ . Based on the eigenvalue equation (2.14), we examine the boundary-layer stability under the effect of the wall admittance $Y$ both by asymptotic analysis of the eigenvalue equation in § 4 and by numerical solution in § 5. For the evaluation of the CHF in (2.14), we use an implementation in Matlab by Motygin (Reference Motygin2018). Further details of this implementation are given in Appendix C.

3.1. Equivalence transformation for the PBE

We start with formulating the normal-mode approach for 3-D perturbations:

(3.1) \begin{equation} q(x,y,z,t)=\hat {q}(y){\rm e}^{i(\alpha _{\mathrm {3D}} x+\beta z-\omega _{\mathrm {3D}} t)}, \end{equation}

where $\beta$ denotes the spanwise wavenumber in the $z$ direction. With this approach, the PBE for 3-D perturbations can be derived in the same way as before, reading

(3.2) \begin{equation} \frac {{\rm d}^{2}\hat {\rho }}{{\rm d}y^{2}}+\frac {2\alpha _{\mathrm {3D}}}{\omega _{\mathrm {3D}}-U_0\alpha _{\mathrm {3D}}}\,\frac {{\rm d}U_0}{{\rm d}y}\frac {{\rm d}\hat {\rho }}{{\rm d}y}+\left [M_{\mathrm {3D}}^{2}(\omega _{\mathrm {3D}}-U_0\alpha _{\mathrm {3D}})^{2}-k^{2}\right ]\hat {\rho }=0 , \end{equation}

with $k^2=\alpha _{\mathrm {3D}}^2+\beta ^2$ . Comparison of (3.2) with the 2-D PBE (2.3),

(3.3) \begin{equation} \frac {{\rm d}^{2}\hat {\rho }}{{\rm d}y^{2}}+\frac {2\alpha _{\mathrm {2D}}}{\omega _{\mathrm {2D}}-U_0\alpha _{\mathrm {2D}}}\,\frac {{\rm d}U_0}{{\rm d}y}\frac {{\rm d}\hat {\rho }}{{\rm d}y}+\left [M_{\mathrm {2D}}^{2}(\omega _{\mathrm {2D}}-U_0\alpha _{\mathrm {2D}})^{2}-\alpha _{\mathrm {2D}}^{2}\right ]\hat {\rho }=0 , \end{equation}

reveals a comparable structure of the two latter equations, suggesting the existence of an equivalence transformation. The subscripts ‘2D’ and ‘3D’ indicate that the equations describe the two-dimensionally and three-dimensionally evolving perturbations, respectively. Here, $M_{2{\rm D}}$ denotes the Mach number at which the 2-D modes occur, whereas $M_{3{\rm D}}$ is the Mach number for the 3-D modes. This means that equations (3.3) and (3.2) describe two different solution spaces. In order to transfer solutions between these solution spaces, so-called equivalence transformations for the PBE are required, by which the 2-D and 3-D equations become form invariant.

For the 2-D and 3-D representations (3.3) and (3.2) of the PBE to be form invariant, the quantities that hold specifically for the 2-D and 3-D solution space must fulfil the relations

(3.4a) \begin{align} k^2 & =\alpha _{\mathrm {3D}}^2+\beta ^2=\alpha _{\mathrm {2D}}^2,\end{align}

(3.4b) \begin{align} M_{\mathrm {2D}}^{2}(\omega _{\mathrm {2D}}-U_0\alpha _{\mathrm {2D}})^{2} & =M_{\mathrm {3D}}^{2}(\omega _{\mathrm {3D}}-U_0\alpha _{\mathrm {3D}})^{2},\end{align}

(3.4c) \begin{align} \frac {\alpha _{\mathrm {2D}}}{\omega _{\mathrm {2D}}-U_0\alpha _{\mathrm {2D}}} & =\frac {\alpha _{\mathrm {3D}}}{\omega _{\mathrm {3D}}-U_0\alpha _{\mathrm {3D}}}. \end{align}

Combining (3.4b ) and (3.4c ) yields

(3.5) \begin{equation} \frac {M_{\mathrm {2D}}^2}{M_{\mathrm {3D}}^2}=\frac {\alpha _{\mathrm {3D}}^2}{\alpha _{\mathrm {2D}}^2}=\varPsi ^2 , \end{equation}

where $\varPsi =M_{\mathrm {2D}}/M_{\mathrm {3D}}$ is introduced for the Mach number ratio between the 2-D and 3-D solution space. Using $\varPsi$ , the transformation relations (3.4) can be written as

(3.6a) \begin{align} \alpha _{\mathrm {3D}}=\pm \,\varPsi \,\alpha _{\mathrm {2D}}, \qquad \omega _{\mathrm {3D}}=\pm \,\varPsi \,\omega _{\mathrm {2D}} ,\end{align}

(3.6b) \begin{align} \beta =\pm \,\sqrt {1-\varPsi ^2}\,\alpha _{\mathrm {2D}}.\end{align}

Equations (3.6) form an equivalence transformation of the PBE, which means that they convert the 3-D PBE (3.2) into the 2-D PBE (3.3). For temporal stability considerations, which presuppose $\alpha ,\beta \in \mathbb {R}$ , the $\beta$ -relation (3.6b) leads to the restriction $\varPsi \lt 1$ .

The advantage of the equivalence transformation (3.6) is that it allows 3-D eigenvalues $(\alpha _{\mathrm {3D}}, \beta , \omega _{\mathrm {3D}})$ at Mach number $M$ to be derived from 2-D eigenvalues $(\alpha _{\mathrm {2D}}, \omega _{\mathrm {2D}})$ at a lower Mach number $\varPsi \!\cdot \! M$ . The shift of the Mach number becomes clear when inserting transformation (3.6) into the 3-D PBE (3.2), giving

(3.7) \begin{equation} \frac {{\rm d}^{2}\hat {\rho }}{{\rm d}y^{2}}+\frac {2\alpha _{\mathrm {2D}}}{\omega _{\mathrm {2D}}-U_0\alpha _{\mathrm {2D}}}\,\frac {{\rm d}U_0}{{\rm d}y}\frac {{\rm d}\hat {\rho }}{{\rm d}y}+\left [\left (M\cdot \varPsi \right )^{2}(\omega _{\mathrm {2D}}-U_0\alpha _{\mathrm {2D}})^{2}-\alpha _{\mathrm {2D}}^{2}\right ]\hat {\rho }=0. \end{equation}

Furthermore, by means of (3.6), an infinite number of 3-D solutions can be obtained from one 2-D solution set by varying the Mach number ratio $\varPsi$ .

To ensure that the transformation (3.6) of the PBE constitutes an equivalence transformation of the full eigenvalue problem, the BCs also have to be taken into account. For 2-D perturbations the BCs read

(3.8) \begin{align} \lim _{y\to \infty }\hat {\rho }(y)&=0 , \end{align}

(3.9) \begin{align} \left.\frac {{\rm d}\hat {\rho }}{{\rm d}y}\right |_{y=0}&=-i\,\omega _{\mathrm {2D}}\,\hat {\rho }(0)\left.Y(\omega )\right |_{\omega =\omega _{\mathrm {2D}}} ,\end{align}

as already derived in the previous section. To convert the BCs into the 3-D solution space, we insert the transformation (3.6). This leaves the far-field BC (3.8) unchanged. The BC at the wall, on the other hand, is modified to

(3.10) \begin{equation} \left.\frac {{\rm d}\hat {\rho }}{{\rm d}y}\right |_{y=0}=\pm \,i\,\varPsi ^{-1}\,\omega _{\mathrm {3D}}\,\hat {\rho }(0)\left.Y\left (\omega \right )\right |_{\,\omega =\mp \,\varPsi ^{-1}\,\omega _{\mathrm {3D}}} , \end{equation}

where the parameter $\varPsi$ now appears. Only in the case of rigid walls with $Y=0$ is the wall BC equivalent in the 2-D and 3-D space. Therefore, only for rigid walls does (3.6) constitute an equivalence transformation of the eigenvalue problem. Consequently, the following conclusions from the equivalence transformation only apply to rigid-wall boundary layers.

3.2. Conclusions for rigid-wall boundary layers

The equivalence transformation allows us to deduce conclusions regarding the stability behaviour of the 2-D compared with the 3-D modes. We discuss two aspects in the following. Firstly, we deduce a statement regarding the critical Mach numbers for the occurrence of 2-D and of 3-D instabilities. Secondly, we derive a condition under which 2-D modes represent the most unstable flow perturbations.

The term critical Mach number $M_{\mathrm {crit}}$ refers to the Mach number above which a certain instability mode occurs, analogous to the critical Reynolds number in the viscous case. From the transformation for the Mach numbers (3.5) with the constraint $\varPsi \lt 1$ in the temporal stability case, it can now be seen that if such a critical Mach number exists for a certain instability, then the critical 2-D Mach number $M_{\mathrm {crit,\,2D}}$ for the occurrence of the 2-D evolving instability is smaller than the critical 3-D Mach number $M_{\mathrm {crit,\,3D}}$ for the corresponding 3-D instability. This means that 2-D instabilities occur at a lower physical Mach number than 3-D instabilities.

In addition to the latter conclusion about the critical Mach number, it is also possible to use the equivalence transformation to obtain a condition under which among all modes the most unstable one is 2-D evolving. This situation can be described by the relation

(3.11) \begin{equation} \left.\omega _{\mathrm {2D}, i, \mathrm {max}}\right |_M \gt \left.\omega _{\mathrm {3D}, i, \mathrm {max}}\right |_M , \end{equation}

which says that the maximum growth rate $\omega _{\mathrm {2D}, i, \mathrm {max}}$ of all 2-D modes over all wavenumbers at Mach number $M$ should be greater than the maximum growth rate $\omega _{\mathrm {3D}, i, \mathrm {max}}$ of all 3-D modes at the same Mach number. Thus, (3.11) describes the case where in a flow with physical Mach number $M$ , the most unstable mode is 2-D.

Applying transformation (3.6) allows us to express the most unstable 3-D growth rates in terms of the most unstable 2-D growth rates according to

(3.12) \begin{equation} \left.\omega _{\mathrm {3D}, i, \mathrm {max}}\right |_M =\varPsi \,\left.\omega _{\mathrm {2D}, i, \mathrm {max}}\right |_{(M_{\mathrm {2D}}=\varPsi \cdot M)}. \end{equation}

With (3.12), (3.11) can be written completely in the 2-D solution space in the form

(3.13) \begin{equation} \frac {\left.\omega _{\mathrm {2D}, i, \mathrm {max}}\right |_M}{\left.\omega _{\mathrm {2D}, i, \mathrm {max}}\right |_{(\varPsi \cdot M)}}\gt \varPsi , \end{equation}

where $\left.\omega _{\mathrm {2D}, i, \mathrm {max}}\right |_M$ denotes the most unstable 2-D growth rate at Mach number $M$ and $\left.\omega _{\mathrm {2D}, i, \mathrm {max}}\right |_{(\varPsi \cdot M)}$ the most unstable 2-D growth rate at smaller Mach number $\varPsi \cdot M$ . Therefore, if condition (3.13) is fulfilled for all $\varPsi \in (0,1)$ , requirement (3.11) is met that the most unstable mode is 2-D evolving. Equation (3.13) thus provides a condition allowing one to assess on the basis of the 2-D modes alone whether they are dominant. In this case, considering only 2-D perturbations is permissible. It is noted that condition (3.13) is obviously fulfilled if $\left.\omega _{\mathrm {2D}, i, \mathrm {max}}\right |_M$ grows monotonically for $M$ .

It should be recalled, however, that the equivalence transformation and therefore the latter conclusions only hold for the case of boundary layers with rigid walls. As soon as an impedance wall with $Y\neq 0$ is introduced, the transformations no longer represent equivalence relations of the eigenvalue problem. Thus, condition (3.13) is not applicable in the case of $Y\neq 0$ , making the occurrence of more unstable 3-D modes conceivable.

4.1. Leading-order eigenvalue equation

For an asymptotic analysis of the eigenvalues $\omega$ for small wavenumbers $\alpha$ , a power series expansion of $\omega$ in the small parameter $\alpha$ is formulated according to

(4.1) \begin{equation} \omega =a_{1}\alpha +a_{2}\alpha ^{2}+a_{3}\alpha ^{3}+\cdots +a_{n}\alpha ^{n}+O(\alpha ^{n+1}). \end{equation}

The constant coefficients $a_m, m=1,2,\dots ,n$ are to be determined in such a way that (4.1) provides a good approximate solution of the eigenvalue equation (2.14) for small $\alpha$ . To this end, (4.1) is inserted into the eigenvalue equation (2.14) and all terms of the equation are then developed as a power series in $\alpha$ . This requires a series expansion of the CHF in (2.14), as presented in Appendix D.1. By subsequently collecting the coefficients of the powers of $\alpha$ in the resulting eigenvalue equation, algebraic equations for determining the coefficients $a_m$ as functions of $M$ and $Y$ are obtained. Inserting these coefficients $a_m$ back into the asymptotic approach (4.1) provides an eigenvalue approximation for small $\alpha$ . For a more detailed explanation of the procedure for the asymptotic analysis, refer to Appendix D.1.

In the context of this asymptotic analysis, it is assumed that the terms of higher-power-order of $\alpha$ are negligible due to the consideration of small wavenumbers $\alpha$ . We therefore restrict ourselves in the following to the investigation of eigenvalues $\omega$ , which are dominated by the leading-order term of the expansion (4.1), meaning $\omega \approx a_1 \alpha$ . We denote these leading-order eigenvalues by $\omega ^{(1)}$ , i.e.

(4.2) \begin{equation} \omega ^{(1)}= a_1 \alpha. \end{equation}

It should be noted that the neglect of the higher-order terms requires assumptions on the convergence of the power series representation of the Heun function, as discussed in Appendix D.2.

To determine the coefficients $a_1$ of the leading-order eigenvalues $\omega ^{(1)}$ , an algebraic equation is derived according to the procedure of collecting the coefficients of the power of $\alpha$ described in Appendix D.1. This equation for $a_1$ can be written in the form

(4.3) \begin{equation} f_{1}(a_{1};M,Y)=a_{1}\left (a_{1}\sqrt {1-(a_{1}-1)^{2}M^{2}}-iY(a_{1}-1)^{2}\right )=0. \end{equation}

The notation $f_1$ refers to (D4) in Appendix D. It should be noted that for $Y=0$ equation (4.3) reduces to a form equivalent to the equation in Zhang & Oberlack (Reference Zhang and Oberlack2021). Due to the branching behaviour of the root functions in (4.3), the solutions $a_1$ determining the leading-order eigenvalues (4.2) behave in a branching manner depending on the parameters $M$ and $Y$ . We therefore want to examine how the solutions $a_1$ and thus the leading-order eigenvalues $\omega ^{(1)}$ change under variation of $Y$ , with a focus on possible branching which can lead to a change in the number and character of the eigenvalues (§ 4.3). For this purpose, we first classify the different solutions $a_1(Y, M)$ of (4.3) in order to determine the $Y$ values depending on $M$ at which branching occurs (§ 4.2). This root classification is based on a reduction of (4.3) to a quartic function, similar to the work of Rienstra (Reference Rienstra2003), where the asymptotic approximation yields a quartic function that allows a classification of the number and character of the surface modes.

4.2. Root classification of the leading-order term

The aim is to formulate criteria for classifying the roots $a_1$ determining the leading-order eigenvalues $\omega ^{(1)}$ . To do this, we consider the leading-order equation (4.3), where we truncate the outer factor $a_1$ , giving the trivial solution. Converting the latter equation by bringing the term $iY(a_{1}-1)^{2}$ to the right-hand side of the equation and subsequently squaring it, we obtain a quartic equation for $a_1$ in the form

(4.4) \begin{equation} \big(M^2-Y^2\big) a_1^4+\big(-2 M^2+4 Y^2\big) a_1^3+\big(M^2-6 Y^2-1\big) a_1^2+4 Y^2 a_1-Y^2=0. \end{equation}

It should be noted that squaring can lead to additional solutions. Therefore, to filter out incorrect solutions, the solutions of (4.4) have to be subsequently reinserted into (4.3). In order to derive classification criteria, (4.4) is transformed in the next step according to Dickson (Reference Dickson1914) to a reduced form without cubic term. For this, we employ the transformation

(4.5) \begin{equation} a_1=\tilde {a}_1-\frac {b}{4} \qquad \mathrm {with}\quad b=\frac {\left (-2 M^2+4 Y^2\right )}{M^2-Y^2} , \end{equation}

where $b$ results from the coefficient of the cubic term in (4.4) divided by the coefficient of the highest power. Applying (4.5) to (4.4) yields the reduced quartic equation

(4.6) \begin{equation} \tilde {a}_1^4+q\,\tilde {a}_1^2+r\,\tilde {a}_1+s=0. \end{equation}

The coefficients $q$ , $r$ , $s$ of (4.6) are specified as functions of $M$ and $Y$ according to

(4.7a) \begin{align} q & =\frac {-M^4+\left (-2 Y^2-2\right ) M^2+2 Y^2}{2\left (M^2-Y^2\right )^2}, \end{align}

(4.7b) \begin{align} r & =\frac {\left (Y^2-1\right ) M^4+3 M^2 Y^2-2 Y^4}{\left (M^2-Y^2\right )^3}, \end{align}

(4.7c) \begin{align} s & =\frac {M^8+\left (-4 Y^2-4\right ) M^6+20 M^4 Y^2-32 M^2 Y^4+16 Y^6}{16\left (M^2-Y^2\right )^4} , \end{align}

having the singularities $Y=\pm M$ induced by $b$ in (4.5).

For real coefficients $q$ , $r$ , $s$ , which are given according to (4.7) in the case of real admittances $Y\in \mathbb {R}$ , the reduced form (4.6) allows us to formulate criteria to classify its roots $\tilde {a}_1$ based on the discriminant $\varDelta$ (see Rees Reference Rees1922). For the reduced quartic polynomial (4.6), the discriminant can be expressed in the form (Prodanov Reference Prodanov2022)

(4.8) \begin{equation} \varDelta =16 q^4 s-4 q^3 r^2-128 q^2 s^2+144 q r^2 s-27 r^4+256 s^3. \end{equation}

Using the relations (4.7), $\varDelta$ is obtained as a function of the parameters $M$ and $Y$ :

(4.9) \begin{equation} \varDelta =\frac {16 Y^2\left [16 Y^4+\left (-8 M^4-20 M^2+1\right ) Y^2+\left (M^8-3 M^6+3 M^4-M^2\right )\right ]}{(M-Y)^6(M+Y)^6}. \end{equation}

As stated in Rees (Reference Rees1922), the roots $\tilde {a}_1$ of (4.6) can now be classified in terms of number and character by means of the sign of the discriminant $\varDelta$ . As a consequence, the number or character of the roots $\tilde {a}_1$ only changes at the zeros of $\varDelta$ , which thus represent branch points of the solutions $\tilde {a}_1$ . From (4.9), the non-trivial zeros of $\varDelta$ can be found at

(4.10a) \begin{align} &Y_{1,4}=\pm \frac {\sqrt {-2+16 M^4+40 M^2+2 \sqrt {\left (8 M^2+1\right )^3}}}{8} \quad \mathrm {for}\quad M\geqslant 0, \end{align}

(4.10b) \begin{align} &Y_{2,3}=\pm \frac {\sqrt {-2+16 M^4+40 M^2-2 \sqrt {\left (8 M^2+1\right )^3}}}{8} \quad \mathrm {for}\quad M\geqslant 1. \end{align}

The restriction of $Y_{2,3}$ to $M\geqslant 1$ follows from the fact that the classification criteria only apply for real admittances $Y$ , as explained above. The key point now is that, due to (4.5), the classification criteria for the solutions $\tilde {a}_1$ apply analogously to the solutions $a_1$ determining the leading-order eigenvalues $\omega ^{(1)}$ . Consequently, branching of the leading-order eigenvalues $\omega ^{(1)}$ occurs at precisely these values $Y_{j}$ given by (4.10).

Figure 1.

Sign of the discriminant $\varDelta$ as main classification criterion: real roots $Y_j, j=1,\dots ,4$ (black curves), where $\varDelta =0$ , dividing $\varDelta \gt 0$ (white areas) and $\varDelta \lt 0$ (blue areas); singularities at $Y=\pm M$ (dot-dashed grey lines).

In figure 1, the branching points $Y_j$ as a function of $M$ are shown as black lines, along which $\varDelta =0$ . They divide the parameter space $(M, Y)$ into areas $\varDelta \gt 0$ and $\varDelta \lt 0$ . Additionally, the singularities $Y=\pm M$ of the discriminant are plotted as grey dot-dashed lines. We see that in the subsonic case $M\lt 1$ only the two branching points $Y_{1,4}$ occur, while the further branching points $Y_{2,3}$ only arise for $M\gt 1$ according to (4.10). In addition, it can be seen that the branching point lines are symmetric to the $Y=0$ axis, which indicates that the eigenvalues have the same branching behaviour for passive walls $Y\gt 0$ as for active walls $Y\lt 0$ . This is consistent with the result of the symmetry analysis in Appendix B, where we found that the transformation (B13) represents a symmetry of the eigenvalue problem for constant real admittances $Y$ , from which it follows that the eigenvalues for $-Y$ result from the eigenvalues for $Y$ by reflection. This implies that the solutions for $Y$ and $-Y$ have the same branching behaviour.

Regarding the results in this section, it should be emphasised again that squaring the leading-order term to convert it into a fourth-order polynomial has led to an expansion of the solution space. This means that potentially not all solutions of the fourth-order polynomial (4.4), on which the root classification is based, are in fact solutions $a_1$ of the leading-order equation (4.3), implying that less than four asymptotic eigenvalues $\omega ^{(1)}$ occur. Indeed, this is what we observe in the following subsection, where we consider the leading-order eigenvalues $\omega ^{(1)}$ for different Mach numbers $M$ under variation of $Y$ . Nevertheless, the achievement of the classification to exactly predict the occurrence of eigenvalue branch points is of key importance and manifests itself in the results for $\omega ^{(1)}$ in the next subsection.

4.3. Leading-order eigenvalues under admittance variation

Based on the previous asymptotic expansion, we now calculate the eigenvalues for small wavenumbers $\alpha$ . Again, we focus on the leading-order eigenvalues $\omega ^{(1)}$ , given by the coefficients $a_1=\omega ^{(1)}/\alpha$ according to (4.2). Thus, we determine the solutions $a_1$ of the leading-order eigenvalue equation (4.3) for different Mach numbers $M$ under variation of the admittance $Y$ . Here we restrict ourselves to constant real admittances in order to be able to compare the results with the predictions of the root classification in the previous section. In the following, we are particularly interested in whether new instability modes can occur due to the wall admittance $Y$ . According to the previous section, such an emergence of new leading-order modes $\omega ^{(1)}$ is to be expected for certain admittance values, the discussed branching points $Y_j$ .

Figure 2.

Solutions $a_1$ for leading-order eigenvalues $\omega ^{(1)}$ plotted over positive wall admittance $Y$ : imaginary parts $a_{1,i}$ (left), real parts $a_{1,r}$ (right). Symbols $\blacklozenge$ and $\blacksquare$ : branching points $Y_1$ , $Y_2$ and singularity $Y=M$ . Mach numbers (a) $M=0.8$ , (b) $M=2$ and (c) $M=4.2$ .

Figure 2 shows for three different Mach numbers ( $M=0.8$ , $M=2$ , $M=4.2$ ) the solutions $a_1$ plotted as curves over increasing positive real $Y$ , starting from the rigid-wall case $Y=0$ . In order to be able to identify the effects of $Y$ on both the stability and the phase velocity of the modes, we have plotted the imaginary part (left-hand panels) and the real part (right-hand panels) of the solutions $a_1$ separately over $Y$ . Therefore, the left-hand panels provide information about the temporal growth rates $a_{1,i}=\omega ^{(1)}_i/\alpha$ of the leading-order modes, while the right-hand panels show the corresponding phase velocity $a_{1,r}=\omega ^{(1)}_r/\alpha$ . The individual solutions $a_1$ are plotted in different colours to allow assignment of the real parts $a_{1,r}$ to the corresponding imaginary parts $a_{1,i}$ .

To validate the results $a_1$ of the small-wavenumber asymptotics, we also solved the eigenvalue equation numerically to a small wavenumber $\alpha =10^{-3}$ . Here, several admittance values $Y$ were considered for all three Mach numbers $M=0.8$ , $M=2$ , $M=4.2$ . These numerical solutions $ {\omega }/{\alpha }$ are plotted in figure 2 as black crosses.

Figure 2 reveals that the behaviour and the number of asymptotic solutions $a_1$ change at selected admittance values, marked on the abscissa by black squares and black diamonds. The black square indicates the singularity $Y=M$ in the leading-order equation (4.3). Here, an additional solution (green) with vanishing imaginary part $a_{1,i}=0$ occurs, as can be observed for all three Mach-number cases. It should be noted at this point that for the asymptotic expansion (4.1) applied to the eigenvalue equation (2.14), real solutions only arise in the leading order $a_1$ . The solutions of the higher-order terms $a_m, m\geqslant 2$ , however, always have imaginary parts, resulting in the asymptotic eigenvalue (4.1) being complex in total. Consequently, the solutions $a_1$ with vanishing imaginary part $a_{1,i}=0$ only represent neutrally stable eigenvalues in the leading order, referred to in the following as neutrally stable leading-order eigenvalues. The black diamonds in figure 2 mark the branching points, given by (4.10), which were derived from the root classification in the previous section. According to (4.10), there are two positive branching points $Y_1$ , $Y_2$ for $M\gt 1$ (see figure 2 b,c), whereas there is only the one branching point $Y_1$ for $M\lt 1$ (see figure 2 a).

5.1. Numerical solution scheme

Before discussing the results, the numerical method used to solve the eigenvalue equation is briefly explained in advance. As discussed in Appendix C, the exact solution occurring in the eigenvalue equation (2.14) contains the CHF. Since the eigenvalue equation represents the impedance wall BC, it requires the evaluation of the CHF at the wall, i.e. at $y=0$ . For this evaluation, we use a Matlab implementation of the CHF, developed by Motygin (Reference Motygin2018). As this implementation has branch cuts on the real axis, the argument of the CHF, $z=\alpha /(\alpha -\omega )$ in the present work, must have an imaginary part. This is obviously the case if stable or unstable modes with $\omega _i\neq 0$ are considered. Vice versa, no conclusions on neutrally stable modes are made in this work and we focus on stable and unstable modes.

These modes are given as complex zeros $\omega =\omega _r+i\,\omega _i$ of the complex eigenvalue equation for a given set of parameters $\alpha$ , $M$ , $Y$ . Due to the structure of the eigenvalue equation with $\omega$ occurring in both the parameters and the argument of the CHF, it cannot be resolved explicitly for $\omega$ . Therefore, eigenvalues are determined numerically by a sign change algorithm. It uses the fact that the eigenvalues are given at points where in the vicinity both the imaginary and the real parts of the eigenvalue equation change their sign. Consequently, the search domain is discretised first, yielding a grid in $\omega _r$ and $\omega _i$ . Thereafter, sign changes are detected on the grid. However, the precision of the eigenvalues found this way depends on the refinement of the grid. Therefore, to improve the precision of the results, the solutions of the sign change algorithm are subsequently used as initial guesses for a Mueller algorithm. This is a root-finding algorithm similar to the Newton method, but with convergence of almost order $2$ . With this calculation method, we ensure high precision and avoid the occurrence of spurious modes.

In order to verify the results, we compare them in selected cases with eigenvalues from a collocation method, adapted for an isothermal boundary layer. Since inviscid equations include various singularities, collocation methods would fail if applied on these equations. Therefore, we validate against results from the viscous compressible linearised Navier–Stokes equations at very high Reynolds number $Re=2\times 10^7$ . For this purpose, a Chebyshev collocation method is applied to the viscous boundary-layer stability problem, where again the base flow is assumed to be isothermal with the exponential main flow profile (2.7). The BCs are set as follows: in addition to the impedance wall BC (2.6), no-slip as well as vanishing temperature fluctuations are assumed at the wall. In the far field it is demanded that the perturbations vanish according to (3.8). Stretching the Chebyshev domain (Schmid & Henningson Reference Schmid and Henningson2001) allows half of the collocation points to be set below a certain limit in order to resolve the near-wall region, in which larger gradients of the eigenfunctions occur, more accurately. The number of collocation points is set to 250. Further details of the method and an application to non-isothermal boundary-layer flows over highly cooled impedance walls can be found in De Broeck et al. (Reference De Broeck, Görtz, Flint, Gonzalez, Oberlack and Lele2022).

The results from both calculations show very good agreement for the most unstable eigenvalues, deviating only of $\mathcal {O}(10^{-3})$ with respect to $\omega _r$ and $\mathcal {O}(10^{-4})$ with respect to $\omega _i$ . Also for smaller Reynolds numbers, such as $Re=10^5$ , the results agree quite well. Examples for the most unstable eigenvalues resulting from both calculation methods in the case of a rigid wall are given in table 1.

Table 1.

Validation of the most unstable inviscid eigenvalue (EV) from the sign change algorithm with a re-iteration using the Mueller algorithm against the most unstable viscous EV from the collocation method for different parameters at $Y=0$ .

5.2. Results for constant wall admittances

Using the method based on the sign change and Mueller algorithm described above, we examine the effects of constant real wall admittances $Y$ in this section. The results presented hereinafter are structured as follows. First, we investigate the admittance effects on the 2-D modes in supersonic boundary layers (§ 5.2.1). The new instability observed for a sufficiently large admittance value Y is subsequently investigated in more detail with respect to the critical admittance value leading to its emergence (§ 5.2.2). Finally, the unstable growth of 2-D modes and 3-D modes is compared, whereby we consider both the original boundary-layer instability modes and the newly occurring instability (§ 5.2.3).

5.2.1. Damping effect and occurrence of a new instability for 2-D modes

The influence of metasurfaces on the stability and the phase velocity of all modes in high-velocity boundary layers should be fundamentally examined, focusing here first on the 2-D modes. For this purpose, we compare the 2-D eigenmodes for constant admittances $Y\neq 0$ with those in the rigid-wall case $Y=0$ with regard to their temporal growth rate $\omega _i$ and frequency $\omega _r$ .

In figure 3, the imaginary part $\omega _i$ and real part $\omega _r$ of the eigenvalues $\omega$ in the Mach $M=4.2$ boundary layer are plotted separately over the streamwise wavenumber $\alpha$ for three admittance values: $Y=0$ , $Y=0.2$ , $Y=10$ .

Figure 3.

Growth rate $\omega _i$ (left) and frequency $\omega _r$ (right) of the 2-D perturbations plotted over wavenumber $\alpha$ for Mach number M = 4.2 and three different constant wall admittances Y: (a) rigid wall Y = 0, (b) Y = 0.2, (c) Y = 10.

The left-hand panels showing $\omega _i$ over $\alpha$ reveal that in the rigid-wall case $Y=0$ the unstable and stable solutions occur in a complex conjugate manner, as predicted by the symmetry analysis in Appendix B. Further, we observe in the left-hand panel of figure 3(a) the maximum growth of the most unstable mode (plotted in yellow) at $\alpha \approx 1.5$ , i.e. at a wavenumber of $1.5$ times the boundary-layer thickness. In this context, however, it is crucial to mention that scaling the decay rate of the velocity profile, i.e. introducing the profile $U_0 (\frac {1}{c} y )$ , leads to scaled eigenvalues and eigenfunctions $ ( ({1}/{c}\alpha) , ({1}/{c}\omega) , (c^2\hat {\rho }) )$ , provided that constant admittances $Y$ are assumed (see Appendix E). In other words, scaling the results $(\alpha , \omega , \hat {\rho })$ with a constant factor again yields eigenvalues and eigenfunctions of the boundary-layer problem with a rescaled velocity profile. Consequently, the results presented here for constant $Y$ are to be assessed qualitatively in terms of fundamental effects, not quantitatively. In previous experimental studies by Stetson & Kimmel (Reference Stetson and Kimmel1993) or Kimmel et al. (Reference Kimmel, Demetriades and Donaldson1996), the maximum growth rate of the second mode was found at a wavenumber of about half the layer thickness.

If we now look at figure 3(b,c) for $Y\neq 0$ , we see that the eigenvalues are no longer complex conjugate. Accordingly, the constant admittance breaks this coupling of stable and unstable solutions, as already observed in Appendix B. Furthermore, the left-hand panel of figure 3(b) reveals that a stabilisation of all modes could be achieved by $Y=0.2$ .

When increasing the admittance to $Y=10$ , however, the emergence of a new instability (black) at small wavenumbers can be observed in figure 3(c), which confirms the results from the asymptotics. This phenomenon is of central physical importance and is therefore examined in more detail in the following section.

From the right-hand panels for $\omega _r$ over $\alpha$ , it can be seen that in the rigid-wall case $Y=0$ (see figure 3 a) the modes propagate with phase velocities $\omega _r/\alpha \in [0,1]$ as soon as they show pronounced unstable growth rates. This is consistent with the critical layer theory as explained previously. The same applies to the newly occurring instability (black) in the case $Y=10$ , which propagates with $\omega _r/\alpha \approx 1-(1/M)$ (see inset in figure 3 c). Comparing the $\omega _r$ – $\alpha$ plots of the three different $Y$ values reveals that the phase velocities of the modes, already initially occurring for a rigid wall, are only slightly affected by the smaller wall admittance $Y=0.2$ . For high wall admittances, such as $Y=10$ , however, we observe an acceleration of the stabilised modes. As a result, they reach phase velocities slightly above $1+ ({1}/{M})$ for small $\alpha$ , as can be seen exemplarily in the insets of the two lowest modes in the right-hand panel of figure 3(c). Further, we see for the two lowest modes that they fall below the $1- ({1}/{M})$ line for smaller $\alpha$ than in the cases $Y=0$ and $Y=0.2$ . In dimensional representation, $1\pm ({1}/{M})$ corresponds to $U_{\infty }\pm c$ , which means that the modes propagate partially supersonically relative to the free-stream velocity.

To assess the character of the instability modes, we consider the associated eigenfunctions with regard to their zero crossings, as described in Mack (Reference Mack1984). For the eigenfunctions, we could derive the exact representation (2.13) from the PBE solution. Plotting the eigenfunctions of the modes in the boundary layer with $Y=0$ reveals that the dominant instability (yellow) in the rigid-wall case represents a second mode. However, this is only the case for wavenumbers where this dominant mode shows clearly unstable growth. In the range of smaller wavenumbers, on the other hand, where only vanishing small growth rates $\omega _i$ occur, this solution represents a first mode. The additional unstable modes with smaller growth rates that occur at higher wavenumbers in the rigid-wall case are further higher acoustic Mack modes. The type of all these modes, which already exist in the rigid-wall boundary layer, does not change when stabilised by wall admittances $Y\neq 0$ .

Of particular interest is the instability emerging for large admittances, which is found to be a first mode. This finding is consistent with previous studies for viscous boundary layers, observing slight destabilisation of the first mode by porous surfaces (Fedorov et al. Reference Fedorov, Shiplyuk, Maslov, Burov and Malmuth2003). The results in this work, though, explicitly show for the inviscid case that the first mode, which already occurs for $Y=0$ , is indeed stabilised by wall admittances, like the higher modes. However, as the admittance is increased, the stimulating effect of the wall coating on the first mode becomes apparent, leading to the appearance of a new unstable first mode.

It should be noted that for the characterisation of the modes based on the eigenfunctions’ zero crossings, as in Mack (Reference Mack1984), those zero crossings are not taken into account which are caused by the oscillatory behaviour of the far-field solution, as can be seen in the eigenfunction plots in Appendix A. The oscillation of the far-field solution, which is given by (2.11), follows from the fact that the square root of the complex number $\theta$ , defined by (5.1), is itself a complex number. This means that the far-field solution can be decomposed according to

(5.2) \begin{equation} \hat \rho _{ff}(y) =C_{2}{\rm e}^{-\sqrt {\theta }y}=C_{2}{\rm e}^{- ( (\sqrt {\theta } )_r+i (\sqrt {\theta } )_i )y} , \end{equation}

revealing its exponential as well as oscillatory behaviour. Here it is worth noting again that the positive square root branch is selected.

5.2.2. Influence of wall admittance on the new instability

In view of the goal of delaying the boundary-layer transition, the phenomenon of the newly occurring first-mode instability due to certain wall admittances is of great practical importance. Therefore, the behaviour of the growth rate $\omega _i$ of the new instability under increasing $Y$ is examined more closely for different Mach numbers in the following ( figure 4). In this context, the question arises in particular as to the critical admittance value above which the new instability sets in. The results presented in figure 4 were calculated for the small wavenumber $\alpha =0.1$ .

Figure 4.

Growth rate $\omega _i$ of the new instability plotted over positive wall admittance $Y$ for different Mach numbers.

As figure 4 reveals, for all Mach numbers $M\gt 1$ the instability actually only emerges in each case from a critical value $Y_{\mathrm {crit,2D}}$ , which depends on the Mach number. We observe that the value $Y_{\mathrm {crit,2D}}$ decreases when $M$ is reduced, meaning that in boundary layers of smaller Mach number, the first-mode instability already occurs at lower admittances. For $M\lt 1$ , the instability even emerges directly from all admittance $Y\gt 0$ . We use the designation $Y_{\mathrm {crit,2D}}$ as the critical value for the occurrence of the 2-D disturbances considered here.

Furthermore, it can be observed that the unstable growth rates $\omega _i$ grow when $Y$ is increased beyond the respective critical value $Y_{\mathrm {crit,2D}}$ . In the subsonic cases $M\lt 1$ , a maximum growth rate occurs at around $Y=1$ before a saturation rate is reached as $Y$ increases. For all $M\gt 1$ , on the other hand, the growth rate increases monotonically with increasing $Y$ up to the saturation value. Both the saturation growth rate for $Y\rightarrow \infty$ as well as the maximum growth rate in the cases $M\lt 1$ all prove to be the greater the smaller is $M$ . Accordingly, a decrease of $M$ leads to earlier occurrence and stronger growth of the new first-mode instability.

These findings are consistent with the results of the asymptotics in § 4.3, which qualitatively observed the same effects of a reduction of $M$ on the critical admittance value, there given by the branching point $Y_2$ , and on the saturation growth rate $\omega _i$ for large $Y$ (see figures 1 and 2).

5.2.3. Comparison of 2-D and 3-D evolving unstable modes

So far, only 2-D modes have been considered. For a comprehensive understanding of the influence of metasurfaces on the boundary-layer stability, the question of the growth of 3-D modes compared with that of 2-D modes is addressed in this section. In this context, condition (3.13) has been derived for the dominating instabilities in boundary layers with rigid wall ( $Y=0$ ) to be 2-D modes. This condition is checked in the following on the basis of the numerical eigenvalues. However, such a condition for 2-D dominance cannot be formulated for admittances $Y\neq 0$ , as shown in § 3, since for $Y\neq 0$ the 2-D and 3-D modes are not linked via a transformation. Therefore, we also compare the growth of 2-D and 3-D modes in the case $Y\neq 0$ .

Figure 5.

Maximum growth rate $\omega _{\mathrm {2D}, i, \mathrm {max}}$ of all 2-D modes plotted over Mach number $M$ for $Y=0$ .

Figure 6.

Growth rate $\omega _i$ of new instability (brown) for $Y=10$ and dominant second mode (black) for $Y = 0$ plotted over wavenumber $\beta$ .

To check condition (3.13) for dominance of the 2-D instability perturbations in the rigid-wall boundary layer, figure 5 shows the maximum growth rate $\omega _{\mathrm {2D}, i, \mathrm {max}}$ , which occurs under all 2-D modes for arbitrary wavenumber, plotted against the Mach number $M$ . This was done by numerically calculating all 2-D modes over a broad wavenumber spectrum for varying Mach number. As one can see from figure 5, for $M\gt 2.4$ , $\omega _{\mathrm {2D}, i, \mathrm {max}}$ grows monotonically with increasing $M$ and is greater than for all smaller Mach numbers. Therefore, in the considered Mach number range, condition (3.13) is fulfilled for all $M\gt 2.4$ , which implies dominance of the 2-D instabilities.

Since the condition (3.13) does not apply for $Y\neq 0$ , we now focus on the new first-mode instability, occurring for admittances $Y\gt Y_{\mathrm {crit,2D}}$ . The aim is to examine whether this mode is more unstable for 2-D or 3-D propagation. For this purpose, in figure 6 the maximum growth rate $\omega _i$ of the new instability (brown) is plotted against the spanwise wavenumber $\beta$ . The results were obtained for $M=4.2$ , $Y=10$ . Wavenumber $\alpha =0.08$ was chosen as the streamwise wavenumber, as for this the instability shows the maximum growth $\omega _i$ when propagating in two dimensions. For comparison, the growth rate $\omega _i$ of the dominant second mode, occurring in the rigid-wall boundary layer at $M=4.2$ , is also shown (black). Here, $\alpha =1.5$ was chosen, representing the point of maximum growth of the 2-D second mode (see figure 3 a).

As we can see from figure 6, the second-mode instability of the case $Y=0$ is most unstable for $\beta =0$ , i.e. when propagating in two dimensions, which confirms condition (3.13). For the newly occurring first-mode instability, on the other hand, maximum growth occurs at $\beta \neq 0$ , namely at $\beta \approx 0.2$ for the parameters considered here. This implies that stronger growth occurs in the 3-D than in the 2-D direction. This can be attributed to the character of this instability as first mode. In the direction of maximum growth, the first mode can even reach growth rates of approximately the maximum growth rate of the dominant rigid-wall second mode, as observed here for $Y=10$ at $M=4.2$ .

5.3. Results for porous wall admittance

So far we have only investigated the effects in the case of constant admittances $Y$ . We now extend our investigations to the practically relevant case of walls with porous coating.

5.3.1. Porous wall model

For the mathematical description of such porous walls, we use a model proposed by Fedorov et al. (Reference Fedorov, Malmuth, Rasheed and Hornung2001), building on the works of Daniels (Reference Daniels1950), Benade (Reference Benade1968) and Stinson & Champoux (Reference Stinson and Champoux1992). According to this model, the porous wall admittance is given by

(5.3) \begin{equation} Y=-\frac {\phi }{Z_0}\mathrm {tanh}\left (\Lambda h\right ). \end{equation}

Here, $h$ represents the porous-layer thickness and $\phi$ the porosity. The latter is given by the ratio of the pore radius $r$ to the pore spacing $s$ according to $\phi =\pi (r/s )^2$ . Thus, $\phi$ lies in the range $\phi \in [0,\pi /4 )$ , since $r$ cannot exceed the limit $r=s/2$ . These quantities were scaled with the boundary-layer thickness $\delta$ . Parameters $Z_0$ and $\varLambda$ describe the characteristic impedance and the propagation constant, which can be expressed by

(5.4) \begin{equation} Z_0 = \sqrt {\frac {Z_1}{Y_1}}, \qquad \varLambda =\sqrt {Z_1 Y_1}, \quad \text {where}\quad {\rm Re} \left (\varLambda \right )\lt 0. \end{equation}

The quantities $Z_1$ and $Y_1$ are given by

(5.5) \begin{equation} Z_1=i\omega \frac {1}{T_w}\frac {J_0(k_v)}{J_2(k_v)}, \qquad Y_1=-i\omega M^2\left [\gamma +(\gamma -1)\frac {J_2(k_t)}{J_0(k_t)}\right ] , \end{equation}

as stated in Fedorov et al. (Reference Fedorov, Malmuth, Rasheed and Hornung2001). Here, $J_0$ and $J_2$ are Bessel functions of the first kind and $T_w$ is the wall-surface temperature, which is $T_w=1$ for the considered isothermal case. The arguments $k_v$ and $k_t$ of the Bessel functions read

(5.6) \begin{equation} k_v= r\sqrt {\frac {i\omega \rho _w}{\mu _w}Re}, \qquad k_t=k_v\sqrt {Pr}, \end{equation}

with $\rho _w$ and $\mu _w$ as the mean density and viscosity at the wall surface. Further, $Pr$ denotes the Prandtl number. The quantities are scaled with the free-stream parameters. Thus, $Re= U_\infty \rho _\infty \delta /\mu _\infty $ is the free-stream Reynolds number. With regard to the dependency of the porous wall admittance (5.3) on the key parameters $h$ and $\phi$ , it is reasonable that for vanishing porous-layer thickness $h=0$ or porosity $\phi =0$ , (5.3) yields the admittance $Y=0$ of rigid walls.

In the context of the inviscid consideration adopted in this paper, we simplify the admittance model for the high-Reynolds-number limit $ Re \to \infty$ . Using the asymptotic expansions of the Bessel functions $J_{\nu }(z)$ for $z \to \infty$ in case $|\mathrm {arg}(z)|\lt \pi$ , as provided in Bender & Orszag (Reference Bender and Orszag1999), it can be derived that $\lim _{z \to \infty }J_0(z)/J_2(z)=-1$ holds. With this, (5.5) simplifies in the inviscid limit to

(5.7) \begin{equation} Z_1=i\omega \frac {1}{T_w}, \qquad Y_1=-i\omega M^2. \end{equation}

Substituting (5.7) into (5.4) yields for $Z_0$ and $\varLambda$ the form

(5.8) \begin{equation} Z_0 = \frac {1}{M\sqrt {T_w}}, \qquad \varLambda =\mathrm {sign}(\omega _i)\frac {i\omega M}{\sqrt {T_w}} , \end{equation}

ensuring ${\rm Re} (\varLambda )\lt 0$ . Plotting the admittance $Y$ of the viscous model against the admittance $Y$ of our simplified model for the inviscid case reveals that the viscous $Y$ approaches the inviscid limit continuously as $k_v$ is increased, justifying the use of the simplified inviscid model.

We are now interested in the effects of a porous layer on the boundary-layer stability. To this end, we investigate the influence of the porous-layer thickness $h$ and the porosity $\phi$ on the instability modes. First, in § 5.3.2 we investigate the influence on the 2-D modes in the $M=4.2$ boundary layer to examine whether the effects we observed for constant admittances $Y$ in § 5.2.1 also hold for porous walls. Subsequently, in § 5.3.3, we also consider 3-D modes. The aim is to investigate the maximum growth with a focus on the maximum growth rate and the corresponding propagation direction under variation of the layer thickness $h$ .

5.3.2. Effects of porous walls on 2-D modes

To study the effects of porous layers on the 2-D instabilities, figure 7 shows the growth rate $\omega _i$ of the 2-D modes plotted against the wavenumber $\alpha$ for different porosity $\phi$ . The results were obtained once for a very large porous-layer thickness $h = 300$ (figure 7 a) and once for a very small layer thickness $h = 0.5$ (figure 7 b), both at a Mach number of $M = 4.2$ . As reference, the curve for the rigid-wall case $\phi =0$ (black curve) is plotted, which represents the dominant second-mode instability in figure 3(a).

Figure 7.

Growth rate $\omega _i$ of the 2-D modes plotted against $\alpha$ in the $M=4.2$ boundary layer for different wall porosities $\phi$ and two different porous-layer thicknesses $h$ : (a) $h=300$ , (b) $h=0.5$ .

In the following, we first discuss the results for the large layer thickness $h=300$ . Subsequently, we compare them with those for the small layer thickness $h=0.5$ .

Effects for large porous-layer thickness

Figure 7(a) reveals that the same effects can be observed for porous walls with large porous-layer thickness $h$ as for constant wall admittance $Y$ in § 5.2. On the one hand, we see that the dominant second-mode instability, which shows significant growth rates in the range of larger wavenumbers $\alpha$ in the rigid-wall case $\phi =0$ , is damped by the thick porous layer for porosities $\phi \gt 0$ , i.e. the growth rates $\omega _i$ in the larger $\alpha$ spectrum are reduced for $\phi \gt 0$ . Increasing the wall porosity $\phi$ leads to a stronger reduction of the second-mode growth rates, until for sufficiently large porosity $\phi$ quasi-complete stabilisation is achieved in the range of large $\alpha$ , as can be seen here for example for $\phi =0.5$ , where the growth rates for high $\alpha$ are of the order of $10^{-5}$ . This damping effect on the second mode, up to almost complete stabilisation for large admittances, is consistent with the observations for constant $Y$ in § 5.2.

In addition to the damping effect, however, we also observe a destabilising effect of the porous wall. Above certain porosities $\phi$ , here in figure 7(a) for $\phi \geqslant 0.1$ , new pronounced unstable growth rates $\omega _i$ occur in the range of smaller $\alpha$ , just as found before for very large constant admittances $Y$ . This new instability at small $\alpha$ is a first mode, as already observed for constant admittances $Y$ . The corresponding eigenfunction is shown in Appendix A. This first-mode instability becomes more unstable with an increase in porosity $\phi$ up to a certain value. However, for all $\phi$ , its growth rates are well below the maximum growth rate of the second mode in the rigid-wall case.

The aforementioned effects, observed in figure 7(a) for the thick porous layer, are qualitatively consistent with previous literature. Fedorov et al. (Reference Fedorov, Shiplyuk, Maslov, Burov and Malmuth2003), for example, also found for boundary layers with porous coating a strong stabilisation of the second mode, dominant in the range of high frequencies, and a slight destabilisation of the first mode at low frequencies.

Effects for small porous-layer thickness

If we look at figure 7(b), we see that walls with a very small porous-layer thickness $h$ also lead to a damping of the second mode at large wavenumbers $\alpha$ , i.e. there is also a reduction in the second-mode growth rate for porosities $\phi \gt 0$ compared with the rigid-wall case $\phi =0$ (black curve). In comparison with the damping effect resulting from the large layer thickness in figure 7(a), however, the second-mode damping is significantly weaker for the small layer thickness $h=0.5$ here in figure 7(b) and complete stabilisation is not achieved, even for high porosity $\phi$ .

Apart from this damping effect, we also observe the occurrence of a new instability for the small porous-layer thickness $h$ in figure 7(b). However, this destabilisation resulting from the thin porous layerdiffers significantly from the destabilisation resulting from thick layers, both in terms of the wavenumber spectrum $\alpha$ , where the destabilisation takes place, and in terms of the magnitude of the growth rates $\omega _i$ of the excited modes. As can be seen from figure 7(b), for very small porosity (here, for example, for $\phi =0.01$ ), the new instability occurs primarily in an $\alpha$ range around a peak (here around $\alpha \approx 0.8$ ), coinciding with the range of the smallest wavenumbers, where the second mode in the rigid-wall case becomes relevant. With increasing $\phi$ , the wavenumber range of the excited instability expands to larger $\alpha$ . At the same time, its growth rates increase significantly. With regard to the increase in the growth rates for increasing $\phi$ , we observe two things. Firstly, we see that there is a critical value of porosity, here slightly greater than $\phi =0.1$ , above which the excited mode shows a larger maximum growth rate than the damped second mode at its peak. Secondly, for even higher $\phi$ , as here for $\phi =0.5$ , the growth rates of the new instability become considerably more unstable than the growth rates of rigid-wall second mode (black curve) over the entire $\alpha$ range. At the same time, there is no longer a pronounced peak with respect to $\omega _i$ . Instead, the maximum growth rate shifts towards larger wavenumbers. This result of the occurrence of instabilities for very thin porous layers with high porosity, which are significantly more unstable than the dominant mode in the rigid-wall case, is quite central.

The eigenfunction of the destabilised mode given in Appendix A suggests that it is a second mode, which is consistent with it occurring at large wavenumbers $\alpha$ . This implies that porous walls with very small layer thickness $h$ lead to the destabilisation of a different type of mode from that observed for large porous-layer thickness $h$ or in the literature where, as discussed above, the destabilisation manifests itself in the form of the first mode while the second mode is stabilised.

5.3.3. Influence of porous-layer thickness on maximum growth of 3-D modes

As we saw in the previous section, for walls with sufficiently high porosity $\phi$ , the dominant growth is no longer given by the original second mode, but by a mode that is destabilised by the porous wall. Both the wavenumbers where this excited instability is relevant and its character, i.e. whether it is a first or second mode, depend strongly on the porous-layer thickness $h$ , while the porosity $\phi$ mainly determines the strength of the stabilising and destabilising effects. However, the results in the previous section were limited to 2-D modes. We therefore now extend the investigations to 3-D modes. The focus is on the propagation direction of the maximum growth of the excited mode and on its maximum growth rate itself. We want to answer the question of how these change with variation of the layer thickness.

For this purpose, figure 8 shows, on the one hand, the maximum growth rate $\omega _{i,\mathrm {max}}$ , occurring over all wavenumbers $\alpha$ and $\beta$ (see figure 8 a), and, on the other hand, the corresponding direction of propagation $\gamma =\mathrm {arctan}(\beta /\alpha )$ of this maximum growth (see figure 8 b), each plotted over the layer thickness $h$ .

Figure 8.

(a) Maximum growth rate $\omega _{i,\mathrm {max}}$ and (b) corresponding direction of propagation $\gamma$ and wavenumbers $\alpha$ , $\beta$ of the maximum growth plotted against the porous-layer thickness $h$ for $M=4.2$ and wall porosity $\phi =0.5$ .

Additionally, the corresponding wavenumbers $\alpha$ and $\beta$ , at which the maximum growth rate occurs, are shown in figure 8(b). The results were obtained for Mach number $M=4.2$ at a relatively high porosity of $\phi =0.5$ , ensuring that the dominant growth rates are given by the excited mode for layer thicknesses $h\gt 0$ . The curve for $h=0$ correspond to the maximum growth of the second mode in the rigid-wall case.

Influence on the maximum growth rate

If we now consider figure 8(a) for the maximum growth rate $\omega _{i,\mathrm {max}}$ , we observe in the range of thin porous layers $h$ that the maximum growth rate increases steeply with increasing $h$ . This suggests the conclusion that thin porous layers can induce growth rates that are higher than the maximum second-mode growth rates in the rigid-wall case $h=0$ . The occurrence of growth rates greater than that for the rigid wall has already been observed in figure 7(b) for small thickness $h=0.5$ .

Figure 8(a) shows that $\omega _{i,\mathrm {max}}$ increases up to a layer thickness marked by $h_{\mathrm {max}}$ . For $h\gt h_{\mathrm {max}}$ , an increase in $h$ leads to a successive reduction in the maximum growth rate, whereby above a layer thickness $\tilde {h}$ damping is achieved compared with the dominant rigid-wall instability. This effect of reducing the maximum growth rate by thick porous layers has already been seen in figure 7(b), where the destabilised first mode shows significantly lower growth rates than the rigid-wall second mode.

Influence on the direction of maximum growth

To answer the question regarding the influence of the layer thickness $h$ on the direction of propagation of the maximum growth, we consider figure 8(b) showing the propagation angle $\gamma$ of $\omega _{i,\mathrm {max}}$ over $h$ .

Most importantly, the figure reveals that the porous-layer thickness $h$ is a decisive factor in determining whether maximum growth occurs in the 2-D or 3-D direction. For porous layers $h$ above a threshold value $h_{\mathrm {3D}}$ the strongest growth occurs in the 3-D direction. This can be attributed to the fact that in these cases $h\gt h_{\mathrm {3D}}$ the amplified instability can be identified as a first mode. The result that the destabilised first mode is most unstable in the 3-D direction is consistent with the results for constant large admittances in § 5.2 as well as with earlier literature on porous walls, such as Fedorov et al. (Reference Fedorov, Shiplyuk, Maslov, Burov and Malmuth2003). More specifically, figure 8(b) indicates that the 3-D angle of the maximum growth of the first mode remains almost the same for all layer thicknesses $h\gt h_{\mathrm {3D}}$ , only increasing marginally with $h$ .

In comparison with the behaviour for $h\gt h_{\mathrm {3D}}$ , on the other hand, it is observed that in the range of smaller porous-wall-layer thicknesses $h\lt h_{\mathrm {3D}}$ , the dominant instability propagates in the 2-D direction. This 2-D character supports the previous interpretation that the mode induced by a very thin porous layer is a second mode, as discussed in the context of figure 7(b).

For layers $h_{\mathrm {max}}\lt h\lt h_{\mathrm {3D}}$ , the eigenfunction of the destabilised mode can be interpreted as a first mode (see Appendix A), despite its 2-D character. However, it should be noted that the characterisation of the eigenfunction is to be discussed.

Comparing figures 8(a) and 8(b), we observe another effect. Relative to the dominant instability for rigid walls $h=0$ , all destabilised 2-D modes are more unstable, while all destabilised 3-D modes with $h\gt \tilde {h}$ are less unstable than the rigid-wall instability. However, for the parameters considered here, 3-D instabilities with growth rates greater than in the rigid-wall case can be observed in a narrow range of mean layer thicknesses $h_{\mathrm {3D}}\lt h\lt \tilde {h}$ .

Influence on the wavenumbers

Regarding the wavenumbers, figure 8(b) reveals that, apart from the range of very small $h$ , $\alpha$ decreases with $h$ , which implies that the maximum growth shifts to smaller streamwise wavenumbers as $h$ increases. This is consistent with the results in figure 7, where it was found for high porosity $\phi$ that in the case of small $h=0.5$ the strongest instability occurs in the form of the destabilised second mode at high wavenumbers $\alpha$ , while for large $h=300$ the destabilised first mode at small $\alpha$ is the dominant mode. In the range of 3-D growth, $\alpha$ remains almost constant, decreasing only slightly. Regarding the spanwise wavenumber $\beta$ of maximum growth, we observe for $h\gt h_{\mathrm {3D}}$ that it occurs around $\beta \approx 0.2$ given the parameters considered. This is consistent with the results for large constant admittances $Y$ in figure 6.

To summarise, three ranges of the porous-layer thickness can be observed with respect to the destabilisation of porous walls with high porosity. Firstly, the range of small layer thicknesses $h$ , where a 2-D instability is induced at high wavenumbers $\alpha$ , leading to growth rates greater than those for rigid walls; secondly, a narrow range of mean layer thicknesses $h$ , where a 3-D mode with growth rates greater than for rigid walls occurs; thirdly, the range of large layer thicknesses $h$ , with the strongest growth occurring in the 3-D direction at small $\alpha$ , which is damped compared with a rigid wall.