2.1. Linearized perturbation equations

According to linear stability theory, the primitive flow variables $\boldsymbol {q} = [\rho , u, v, w, T]^\mathrm {T}$ are split into a steady reference state $\bar {\boldsymbol {q}}$, also known as base flow, and a small unsteady perturbation field $\tilde {\boldsymbol {q}}$:

(2.8)\begin{equation} \boldsymbol{q} = \bar{\boldsymbol{q}} + \epsilon \tilde{\boldsymbol{q}} + O(\epsilon^2), \end{equation}

with $\epsilon \ll 1$. Introducing this splitting into the Navier–Stokes system (2.1) and neglecting nonlinear terms (terms of order $O (\epsilon ^2)$) results in the linearized perturbation equations.

In this work, the base flow is assumed to be locally parallel in the streamwise direction, so that $\bar {\boldsymbol {q}} = \bar {\boldsymbol {q}} (y, z)$ at a given $x$ coordinate, and the wall-normal base-flow velocity component ($\bar {v}$) is assumed to be zero. The unsteady perturbations are three-dimensional, i.e. $\tilde {\boldsymbol {q}} = \tilde {\boldsymbol {q}} (x,y,z,t)$. For these assumptions, the resulting linearized disturbance equations are reported in Appendix A.

2.2. Formulation of the two-dimensional local linear stability problem

To allow the study of the instabilities growing in the wake behind an isolated roughness element, the amplitude of the perturbations is considered to be a function of both the wall-normal and the spanwise directions. This assumption leads to a two-dimensional local linear stability theory (2D-LST), for which the ansatz describing modal perturbations can be written as follows (see for example Theofilis Reference Theofilis2003):

(2.9)\begin{equation} \tilde{\boldsymbol{q}} (x,y,z,t) = \hat{\boldsymbol{q}} (y,z) \exp [\mathrm{i} (\alpha x - \omega t)] + \mathrm{c.c.} = 2 \mathrm{Re} \left\{\hat{\boldsymbol{q}} (y,z) \exp [\mathrm{i} (\alpha x - \omega t)] \right\}, \end{equation}

where $\hat {\boldsymbol {q}}$ is the vector containing the two-dimensional amplitude functions, $\alpha$ is the wavenumber along the streamwise direction, $\omega$ is the angular frequency and c.c. denotes the complex conjugate. Both $\alpha$ and $\omega$ can be complex in general. If $\alpha$ is assumed to be real and $\omega$ is complex, then a temporal stability approach is adopted, in which the imaginary part of $\omega$, denoted by $\omega _i$, represents the temporal growth rate of the perturbation. On the other hand, if $\alpha$ is complex and $\omega$ is real, the spatial evolution of the instabilities is considered instead, and the spatial growth rate is given by $-\alpha _i$, where $\alpha _i$ is the imaginary part of $\alpha$.

By introducing (2.9) into the linearized perturbation equations (A1), the 2D-LST equations are obtained. They constitute a partial differential generalized eigenvalue problem which, depending on the approach, can be expressed in matrix form as

(2.10)\begin{equation} \boldsymbol{\mathsf{A}}_\omega \hat{\boldsymbol{q}} = \omega \boldsymbol{\mathsf{B}}_\omega \hat{\boldsymbol{q}} \end{equation}

for a temporal analysis, or as

(2.11)\begin{equation} \boldsymbol{\mathsf{A}}_\alpha \hat{\boldsymbol{q}} = \alpha \boldsymbol{\mathsf{B}}_\alpha \hat{\boldsymbol{q}} + \alpha^2 \boldsymbol{\mathsf{C}}_\alpha \hat{\boldsymbol{q}} \end{equation}

for a spatial framework. The matrices $\boldsymbol{\mathsf{A}}$, $\boldsymbol{\mathsf{B}}$ and $\boldsymbol{\mathsf{C}}$ contain the coefficients of the 2D-LST system of equations that respectively multiply the zeroth, first and second powers of the eigenvalue in each approach. In the spatial case, it should be noted that the problem is nonlinear in the eigenvalue $\alpha$. To convert the system into a linear eigenvalue problem, the matrix companion method (Bridges & Morris Reference Bridges and Morris1984) is employed, which introduces four auxiliary variables into the vector of amplitude functions: $\alpha \hat {u}$, $\alpha \hat {v}$, $\alpha \hat {w}$ and $\alpha \hat {T}$. This procedure significantly increases the computational cost associated with solving the spatial eigenvalue problem, compared to the temporal one.

3.1. Temporal growth-rate decomposition

The substitution of the perturbation ansatz (2.9) into the disturbance energy (3.1) leads to a decomposition of the temporal growth rate into the different terms that contribute to an increase or decrease of the amplitude of a perturbation, which are directly linked to the production and the dissipation of disturbance energy. To derive this decomposition, a temporal stability framework is employed, such that the disturbances are periodic in the streamwise direction (real streamwise wavenumber $\alpha = \alpha _r$) but allowed to grow or decay in time (complex angular frequency $\omega = \omega _r + \mathrm {i} \omega _i$).

It is worth noting that all the integrands in (3.1) involve products of two disturbance quantities. Following Weder et al. (Reference Weder, Gloor and Kleiser2015), the domain size in the streamwise direction is chosen to be equal to a single wavelength, i.e. $2 {\rm \pi}/ \alpha _r$, which allows the streamwise integration of disturbance products to be simplified as follows (see also Weder Reference Weder2012):

(3.7)\begin{equation} \int_0^{2 {\rm \pi}/ \alpha_r} \tilde{f} \tilde{g}\,\textrm{d} x = \frac{4 {\rm \pi}}{\alpha_r} \exp (2 \omega_i t) \langle \hat{f}, \hat{g} \rangle, \end{equation}

where $\tilde {f}$ and $\tilde {g}$ are two generic disturbance variables, and $\langle \hat {f}, \hat {g} \rangle = (\hat {f}^* \hat {g} + \hat {f} \hat {g}^*)/2$, with $(\cdot )^*$ denoting the complex conjugate. In addition, because of streamwise periodicity, all the terms with volume integrals of the form $\int _\varOmega \partial (\cdot ) / \partial x \,\textrm {d} \mathcal {V}$ or with surface integrals of the form $\int _\varGamma (\cdot ) n_x \,\textrm {d} S$ become zero. Using (3.7), the disturbance energy can be expressed in terms of the amplitude function of the primitive variables as

(3.8)\begin{align} \tilde{E} &= \frac{4 {\rm \pi}}{\alpha_r} \exp (2 \omega_i t) \frac{1}{2} \int_{z^-}^{z^+} \int_{y^-}^{y^+} \left[ \bar{\rho} \left( \langle \hat{u}, \hat{u} \rangle + \langle \hat{v}, \hat{v} \rangle + \langle \hat{w}, \hat{w} \rangle \right) + \frac{\bar{T}}{\gamma M^2 \bar{\rho}} \langle \hat{\rho}, \hat{\rho} \rangle \right.\nonumber\\ &\quad + \left. \frac{\bar{\rho}}{\gamma (\gamma - 1) M^2 \bar{T}} \langle \hat{T}, \hat{T} \rangle \right]\, \textrm{d} y\, \textrm{d} z = \frac{4 {\rm \pi}}{\alpha_r} \exp (2 \omega_i t) \hat{E}, \end{align}

where $\hat {E}$ denotes the perturbation energy amplitude function, $z^-$ and $z^+$ are respectively the lower and upper limits of the domain along the spanwise direction, and $y^-$ and $y^+$ are the corresponding limits along the wall-normal direction. Making use of the relation (3.8), the left-hand side of (3.1) becomes $\textrm {d} \tilde {E} / \textrm {d}t = 4 {\rm \pi}/ \alpha _r \exp (2 \omega _i t) 2 \omega _i \hat {E}$. Employing (3.7) in all the remaining terms of the disturbance energy equation (3.1) and noting that the factor $4 {\rm \pi}/ \alpha _r \exp (2 \omega _i t)$ is common to all terms and cancels out, the following decomposition for the temporal growth rate is obtained:

(3.9)\begin{equation} \omega_i = \frac{1}{2 \hat{E}} ( \hat{P} + \hat{D} + \hat{F} ), \end{equation}

where $\hat {P}$ contains the growth-rate contributions of the energy production terms, which can be either positive or negative depending on the particular conditions; $\hat {D}$ encloses the contributions of the energy dissipation terms, which are always negative; and $\hat {F}$ comprises the contributions due to the energy fluxes across the domain boundary $\varGamma$. Here, the following classification of terms is adopted:

(3.10a)\begin{gather} \hat{P} = \hat{P}_{RS} + \hat{P}_s + \hat{P}_Q + \hat{P}_{mom} + \hat{P}_{\boldsymbol{\nabla} p} + \hat{P}_{dil} + \hat{P}_{s,T} + \hat{P}_{con}, \end{gather}

(3.10b)\begin{gather}\hat{D} = \hat{D}_\mu + \hat{D}_k, \end{gather}

(3.10c)\begin{gather}\hat{F} = \hat{F}_p + \hat{F}_\mu + \hat{F}_k + \hat{F}_{d\mu/dT}, \end{gather}

with

(3.11a)\begin{align} \hat{P}_{RS} &={-}\int_{z^-}^{z^+} \int_{y^-}^{y^+} \bar{\rho} \langle \hat{u}, \hat{v} \rangle \dfrac{\partial \bar{u}}{\partial y}\, \textrm{d} y\, \textrm{d} z - \int_{z^-}^{z^+} \int_{y^-}^{y^+} \bar{\rho} \langle \hat{u}, \hat{w} \rangle \dfrac{\partial \bar{u}}{\partial z}\, \textrm{d} y\, \textrm{d} z\nonumber\\ &\quad - \int_{z^-}^{z^+} \int_{y^-}^{y^+} \bar{\rho} \langle \hat{v}, \hat{v} \rangle \dfrac{\partial \bar{v}}{\partial y}\, \textrm{d} y\, \textrm{d} z - \int_{z^-}^{z^+} \int_{y^-}^{y^+} \bar{\rho} \langle \hat{v}, \hat{w} \rangle \dfrac{\partial \bar{v}}{\partial z}\, \textrm{d} y\, \textrm{d} z\nonumber\\ &\quad - \int_{z^-}^{z^+} \int_{y^-}^{y^+} \bar{\rho} \langle \hat{w}, \hat{v} \rangle \dfrac{\partial \bar{w}}{\partial y}\, \textrm{d} y\, \textrm{d} z - \int_{z^-}^{z^+} \int_{y^-}^{y^+} \bar{\rho} \langle \hat{w}, \hat{w} \rangle \dfrac{\partial \bar{w}}{\partial z}\, \textrm{d} y\, \textrm{d} z, \end{align}

(3.11b)\begin{align} \hat{P}_{s} &={-} \int_{z^-}^{z^+} \int_{y^-}^{y^+} \bar{\rho} \langle \hat{v}, \hat{s} \rangle \dfrac{\partial \bar{T}}{\partial y}\, \textrm{d} y\, \textrm{d} z - \int_{z^-}^{z^+} \int_{y^-}^{y^+} \bar{\rho} \langle \hat{w}, \hat{s} \rangle \dfrac{\partial \bar{T}}{\partial z}\, \textrm{d} y\, \textrm{d} z, \end{align}

(3.11c)\begin{align} \hat{P}_{Q} &= \int_{z^-}^{z^+} \int_{y^-}^{y^+} \frac{\langle \hat{Q}, \hat{T} \rangle}{\bar{T}}\, \textrm{d} y\, \textrm{d} z, \end{align}

(3.11d)\begin{align} \hat{P}_{mom} &={-} \int_{z^-}^{z^+} \int_{y^-}^{y^+} \langle \hat{\rho}, \hat{u} \rangle \bar{v} \dfrac{\partial \bar{u}}{\partial y}\, \textrm{d} y\, \textrm{d} z - \int_{z^-}^{z^+} \int_{y^-}^{y^+} \langle \hat{\rho}, \hat{u} \rangle \bar{w} \dfrac{\partial \bar{u}}{\partial z}\, \textrm{d} y\, \textrm{d} z\nonumber\\ &\quad - \int_{z^-}^{z^+} \int_{y^-}^{y^+} \langle \hat{\rho}, \hat{v} \rangle \bar{v} \dfrac{\partial \bar{v}}{\partial y}\, \textrm{d} y\, \textrm{d} z - \int_{z^-}^{z^+} \int_{y^-}^{y^+} \langle \hat{\rho}, \hat{v} \rangle \bar{w} \dfrac{\partial \bar{v}}{\partial z}\, \textrm{d} y\, \textrm{d} z\nonumber\\ &\quad - \int_{z^-}^{z^+} \int_{y^-}^{y^+} \langle \hat{\rho}, \hat{w} \rangle \bar{v} \dfrac{\partial \bar{w}}{\partial y}\, \textrm{d} y\, \textrm{d} z - \int_{z^-}^{z^+} \int_{y^-}^{y^+} \langle \hat{\rho}, \hat{w} \rangle \bar{w} \dfrac{\partial \bar{w}}{\partial z}\, \textrm{d} y\, \textrm{d} z, \end{align}

(3.11e)\begin{align} \hat{P}_{\boldsymbol{\nabla} p} &={-} \int_{z^-}^{z^+} \int_{y^-}^{y^+} \frac{\langle \hat{\rho}, \hat{v} \rangle}{\bar{\rho}} \dfrac{\partial \bar{p}}{\partial y}\, \textrm{d} y\, \textrm{d} z - \int_{z^-}^{z^+} \int_{y^-}^{y^+} \frac{\langle \hat{\rho}, \hat{w} \rangle}{\bar{\rho}} \dfrac{\partial \bar{p}}{\partial z}\, \textrm{d} y\, \textrm{d} z, \end{align}

(3.11f)\begin{align} \hat{P}_{dil} &={-} \int_{z^-}^{z^+} \int_{y^-}^{y^+} \frac{\langle \hat{p}, \hat{T} \rangle}{\bar{T}} \dfrac{\partial \bar{v}}{\partial y}\, \textrm{d} y\, \textrm{d} z - \int_{z^-}^{z^+} \int_{y^-}^{y^+} \frac{\langle \hat{p}, \hat{T} \rangle}{\bar{T}} \dfrac{\partial \bar{w}}{\partial z}\, \textrm{d} y\, \textrm{d} z\nonumber\\ &\quad - \int_{z^-}^{z^+} \int_{y^-}^{y^+} \frac{\langle \hat{\rho}, \hat{\rho} \rangle}{\bar{\rho}^2} \bar{p} \dfrac{\partial \bar{v}}{\partial y}\, \textrm{d} y\, \textrm{d} z - \int_{z^-}^{z^+} \int_{y^-}^{y^+} \frac{\langle \hat{\rho}, \hat{\rho} \rangle}{\bar{\rho}^2} \bar{p} \dfrac{\partial \bar{w}}{\partial z}\, \textrm{d} y\, \textrm{d} z, \end{align}

(3.11g)\begin{align} \hat{P}_{s,T} &={-} \frac{1}{\gamma \left( \gamma - 1 \right) M^2} \int_{z^-}^{z^+} \int_{y^-}^{y^+} \frac{\langle \hat{\rho}, \hat{T} \rangle}{\bar{T}} \bar{v} \dfrac{\partial \bar{T}}{\partial y}\, \textrm{d} y\, \textrm{d} z\nonumber\\ &\quad - \frac{1}{\gamma \left( \gamma - 1 \right) M^2} \int_{z^-}^{z^+} \int_{y^-}^{y^+} \frac{\langle \hat{\rho}, \hat{T} \rangle}{\bar{T}} \bar{w} \dfrac{\partial \bar{T}}{\partial z}\, \textrm{d} y\, \textrm{d} z, \end{align}

(3.11h)\begin{align} \hat{P}_{con} &={-} \int_{z^-}^{z^+} \int_{y^-}^{y^+} \left[ \frac{\bar{v} \bar{p}}{\bar{\rho}^2} \langle \hat{\rho}, \partial \hat{\rho} / \partial y \rangle + \frac{\bar{v} \bar{p}}{(\gamma - 1) \bar{T}^2} \langle \hat{T}, \partial \hat{T} / \partial y \rangle \right. \nonumber\\ &\quad \left. +\, \bar{v} \bar{\rho} \left( \langle \hat{u}, \partial \hat{u} / \partial y \rangle + \langle \hat{v}, \partial \hat{v} / \partial y \rangle + \langle \hat{w}, \partial \hat{w} / \partial y \rangle \right) \vphantom{\frac{\bar{v} \bar{p}}{\bar{\rho}^2}}\right]\, \textrm{d} y\, \textrm{d} z\nonumber\\ &\quad - \int_{z^-}^{z^+} \int_{y^-}^{y^+} \left[ \frac{\bar{w} \bar{p}}{\bar{\rho}^2} \langle \hat{\rho}, \partial \hat{\rho} / \partial z \rangle + \frac{\bar{w} \bar{p}}{(\gamma - 1) \bar{T}^2} \langle \hat{T}, \partial \hat{T} / \partial z \rangle \right. \nonumber\\ &\quad \left. +\, \bar{w} \bar{\rho} \left( \langle \hat{u}, \partial \hat{u} / \partial z \rangle + \langle \hat{v}, \partial \hat{v} / \partial z \rangle + \langle \hat{w}, \partial \hat{w} / \partial z \rangle \right) \vphantom{\frac{\bar{v} \bar{p}}{\bar{\rho}^2}}\right]\, \textrm{d} y\, \textrm{d} z, \end{align}

(3.11i)\begin{align} \hat{D}_{\mu} &={-} \int_{z^-}^{z^+} \int_{y^-}^{y^+} \left( \langle \hat{\tau}_{xx}, \mathrm{i} \alpha \hat{u} \rangle + \langle \hat{\tau}_{yy}, \partial \hat{v} / \partial y \rangle + \langle \hat{\tau}_{zz}, \partial \hat{w} / \partial z \rangle + \langle \hat{\tau}_{xy}, \partial \hat{u} / \partial y \rangle \right.\nonumber\\ &\quad \left. +\, \langle \hat{\tau}_{xy}, \mathrm{i} \alpha \hat{v} \rangle + \langle \hat{\tau}_{xz}, \partial \hat{u} / \partial z \rangle + \langle \hat{\tau}_{xz}, \mathrm{i} \alpha \hat{w} \rangle + \langle \hat{\tau}_{yz}, \partial \hat{v} / \partial z \rangle \right.\nonumber\\ &\quad \left. +\, \langle \hat{\tau}_{yz}, \partial \hat{w} / \partial y \rangle \right)\, \textrm{d} y\, \textrm{d} z, \end{align}

(3.11j)\begin{align} \hat{D}_{k} &= \int_{z^-}^{z^+} \int_{y^-}^{y^+} \frac{\langle \hat{q}_x, \mathrm{i} \alpha \hat{T} \rangle + \langle \hat{q}_y, \partial \hat{T} / \partial y \rangle + \langle \hat{q}_z, \partial \hat{T} / \partial z \rangle}{\bar{T}}\, \textrm{d} y\, \textrm{d} z, \end{align}

(3.11k)\begin{align} \hat{F}_{p} &={-}\int_{z^-}^{z^+} \left. \langle \hat{p}, \hat{v} \rangle \right|_{y^-}^{y^+} \,\textrm{d} z -\int_{y^-}^{y^+} \left. \langle \hat{p}, \hat{w} \rangle \right|_{z^-}^{z^+}\, \textrm{d} y, \end{align}

(3.11l)\begin{align} \hat{F}_{\mu} &= \int_{z^-}^{z^+} \left. \langle \hat{\tau}_{yy}, \hat{v} \rangle \right|_{y^-}^{y^+} \,\textrm{d} z + \int_{y^-}^{y^+} \left. \langle \hat{\tau}_{zz}, \hat{w} \rangle \right|_{z^-}^{z^+}\, \textrm{d} y + \int_{z^-}^{z^+} \left. \langle \hat{\tau}_{xy}, \hat{u} \rangle \right|_{y^-}^{y^+} \,\textrm{d} z\nonumber\\ &\quad +\, \int_{y^-}^{y^+} \left. \langle \hat{\tau}_{xz}, \hat{u} \rangle \right|_{z^-}^{z^+}\, \textrm{d} y + \int_{z^-}^{z^+} \left. \langle \hat{\tau}_{yz}, \hat{w} \rangle \right|_{y^-}^{y^+} \,\textrm{d} z + \int_{y^-}^{y^+} \left. \langle \hat{\tau}_{yz}, \hat{v} \rangle \right|_{z^-}^{z^+}\, \textrm{d} y, \end{align}

(3.11m)\begin{align} \hat{F}_{k} &={-}\int_{z^-}^{z^+} \left. \frac{\langle \hat{q}_y, \hat{T} \rangle}{\bar{T}} \right|_{y^-}^{y^+} \,\textrm{d} z -\int_{y^-}^{y^+} \left. \frac{\langle \hat{q}_z, \hat{T} \rangle}{\bar{T}} \right|_{z^-}^{z^+}\, \textrm{d} y, \end{align}

(3.11n)\begin{align} \hat{F}_{d\mu/dT} &= \frac{1}{Re} \int_{z^-}^{z^+} \left. \langle \hat{u}, \hat{T} \rangle \dfrac{\textrm{d} \bar{\mu}}{\textrm{d} \bar{T}} \dfrac{\partial \bar{u}}{\partial y} \right|_{y^-}^{y^+} \,\textrm{d} z + \frac{1}{Re} \int_{y^-}^{y^+} \left. \langle \hat{u}, \hat{T} \rangle \dfrac{\textrm{d} \bar{\mu}}{\textrm{d} \bar{T}} \dfrac{\partial \bar{u}}{\partial z} \right|_{z^-}^{z^+}\, \textrm{d} y\nonumber\\ &\quad + \frac{1}{Re} \int_{y^-}^{y^+} \left. \langle \hat{v}, \hat{T} \rangle \dfrac{\textrm{d} \bar{\mu}}{\textrm{d} \bar{T}} \dfrac{\partial \bar{v}}{\partial z} \right|_{z^-}^{z^+} \, \textrm{d} y + \frac{1}{Re} \int_{z^-}^{z^+} \left. \langle \hat{w}, \hat{T} \rangle \dfrac{\textrm{d} \bar{\mu}}{\textrm{d} \bar{T}} \dfrac{\partial \bar{v}}{\partial z} \right|_{y^-}^{y^+} \,\textrm{d} z\nonumber\\ &\quad + \frac{1}{Re} \int_{y^-}^{y^+} \left. \langle \hat{v}, \hat{T} \rangle \dfrac{\textrm{d} \bar{\mu}}{\textrm{d} \bar{T}} \dfrac{\partial \bar{w}}{\partial y} \right|_{z^-}^{z^+}\, \textrm{d} y + \frac{1}{Re} \int_{z^-}^{z^+} \left. \langle \hat{w}, \hat{T} \rangle \dfrac{\textrm{d} \bar{\mu}}{\textrm{d} \bar{T}} \dfrac{\partial \bar{w}}{\partial y} \right|_{y^-}^{y^+} \,\textrm{d} z\nonumber\\ &\quad + \frac{2}{Re} \int_{z^-}^{z^+} \left. \langle \hat{v}, \hat{T} \rangle \dfrac{\textrm{d} \bar{\mu}}{\textrm{d} \bar{T}} \dfrac{\partial \bar{v}}{\partial y} \right|_{y^-}^{y^+} \,\textrm{d} z + \frac{2}{Re} \int_{y^-}^{y^+} \left. \langle \hat{w}, \hat{T} \rangle \dfrac{\textrm{d} \bar{\mu}}{\textrm{d} \bar{T}} \dfrac{\partial \bar{w}}{\partial z} \right|_{z^-}^{z^+}\, \textrm{d} y\nonumber\\ &\quad + \frac{1}{Re} \int_{z^-}^{z^+} \left. \langle \hat{v}, \hat{T} \rangle \dfrac{\textrm{d} \bar{\lambda}}{\textrm{d} \bar{T}} \dfrac{\partial \bar{v}}{\partial y} \right|_{y^-}^{y^+} \,\textrm{d} z + \frac{1}{Re} \int_{y^-}^{y^+} \left. \langle \hat{w}, \hat{T} \rangle \dfrac{\textrm{d} \bar{\lambda}}{\textrm{d} \bar{T}} \dfrac{\partial \bar{v}}{\partial y} \right|_{z^-}^{z^+}\, \textrm{d} y\nonumber\\ &\quad + \frac{1}{Re} \int_{z^-}^{z^+} \left. \langle \hat{v}, \hat{T} \rangle \dfrac{\textrm{d} \bar{\lambda}}{\textrm{d} \bar{T}} \dfrac{\partial \bar{w}}{\partial z} \right|_{y^-}^{y^+} \,\textrm{d} z + \frac{1}{Re} \int_{y^-}^{y^+} \left. \langle \hat{w}, \hat{T} \rangle \dfrac{\textrm{d} \bar{\lambda}}{\textrm{d} \bar{T}} \dfrac{\partial \bar{w}}{\partial z} \right|_{z^-}^{z^+}\, \textrm{d} y. \end{align}

The production term $\hat {P}_{RS}$ represents the work done by the disturbance Reynolds stresses against the velocity shear layers present in the base flow. The term $\hat {P}_{s}$ describes the energy produced due to the transport of entropy spottiness across the base-flow temperature shear layers. The next term, $\hat {P}_{Q}$, refers to the production of disturbance energy due to the heat source $\hat {Q}$, which is generated by transport phenomena, including the variation of the transport properties with temperature. The contribution denoted by $\hat {P}_{mom}$ comprises the energy produced as a result of the convection of disturbance momentum across the base-flow velocity shear layers. The presence of a non-negligible pressure gradient in the base flow brings an additional energy production term, labelled as $\hat {P}_{\boldsymbol {\nabla } p}$, which can be interpreted as a form of pressure work. Another source of energy production arises from the interaction between pressure and the divergence of the velocity field, which adopts the shape of a pressure dilatation term, $\hat {P}_{dil}$. The term $\hat {P}_{s,T}$ results from the material derivative of the internal energy fluctuation. The last contribution to the disturbance energy production is comprised inside the term $\hat {P}_{con}$, which represents the remaining part of the convective derivative of the total disturbance energy. The dissipation of perturbation energy due to friction and thermal conduction is contained in the dissipation terms $\hat {D}_{\mu }$ and $\hat {D}_k$, respectively. With regard to the disturbance energy flux, the term $\hat {F}_p$ describes the rate of work done on the domain due to the disturbance pressure force. The terms $\hat {F}_\mu$ and $\hat {F}_k$ respectively embody the net flux of mechanical and thermal energy convected across the domain boundary owing to friction and thermal conduction. Finally, the disturbance energy flux contribution given by $\hat {F}_{d \mu / dT}$ appears as a result of the dependence of viscosity on temperature.

In the following sections, the superscript $d$ employed to denote dimensional quantities is dropped for convenience. Values of quantities presented in dimensional form are unequivocally identified by the accompanying SI units.

4.1. Details of the base-flow simulations

Two different roughness geometries are considered in this work, which are depicted in figure 1. The first geometry is a sharp-edged cuboidal element with a square planform shape. The second one consists in a three-dimensional sharp-edged ramp geometry, which is inspired by popular passive flow control devices employed to promote boundary layer transition (see for instance Stouffer et al. Reference Stouffer, Baker, Capriotti and Northam1993).

Figure 1.

Roughness geometries and computational domain (dashed line). Not to scale.

The geometrical configuration of the computational setup under consideration is also represented in figure 1. A roughness element is mounted on top of a flat plate at a given distance behind the leading edge and centred along the spanwise direction at the streamwise symmetry plane ($z = 0$). The leading edge of the flat plate is assumed to be infinitely sharp. The computational domain is located inside the shock layer established between the surface and the weak shock wave induced at the flat plate leading edge. This approach reduces the computational effort necessary to obtain the base-flow solution, and has already been employed in similar studies in the literature (see for example De Tullio et al. Reference De Tullio, Paredes, Sandham and Theofilis2013; Di Giovanni & Stemmer Reference Di Giovanni and Stemmer2018). The top boundary has an angle to prevent roughness-induced compression waves and weak shock waves from impinging on it, thus avoiding potential reflections back into the boundary layer. In addition, the spanwise symmetry of the geometrical configuration is exploited by considering only half of the roughness element.

The dimensions of the computational domain are the same for both roughness geometries and are summarized in table 1, which lists the streamwise coordinate and height of the inlet ($x_{in}, y_{in}$) and outlet ($x_{out}, y_{out}$) planes of the domain, as well as the spanwise domain size ($z_\infty$), which is constant. Both roughness elements have a height of $h = 0.4\ \textrm {mm}$, and their leading edge is located at a streamwise distance of $x_h = 60\ \textrm {mm}$ from the flat plate leading edge. Their length is determined according to the ramp angle $\theta$ as $d = h/ \tan (\theta )$, with $\theta = 10^{\circ }$. The roughness Reynolds number for this configuration is $Re_h = u_h h / \nu _h = 330$, where $u_h$ and $\nu _h$ are respectively the streamwise velocity and kinematic viscosity of the flow at a height equal to $h$ and at a streamwise location of $x_h$ in a self-similar boundary-layer profile. The corresponding ratio of roughness height to the local boundary-layer thickness is $h/\delta _{99} = 0.59$, where $\delta _{99}$ is estimated by applying the total enthalpy criterion ($h_0/h_{0,\infty } = 0.995$) in a self-similar boundary-layer profile at $x_h$. The Reynolds number based on the displacement thickness of the boundary layer at the domain inflow is $Re_{\delta _1} = 14\,847$.

Table 1.

Dimensions of the computational domain. The origin is located at the flat plate leading edge.

The free-stream values considered for the numerical solution are summarized in table 2. They correspond to the high-Reynolds test conditions of the von Karman Institute H3 hypersonic wind tunnel (Tirtey Reference Tirtey2009). Because of the low free-stream static temperature (total temperature $T_0 = 500\ \textrm {K}$), no significant excitation of the vibrational molecular energy mode is expected in the flow, and therefore the assumption of a calorically perfect gas is a good approximation. The flat plate wall is assumed to be isothermal with a wall temperature of $T_w = 300\ \textrm {K}$. This value corresponds to a cold wall condition, since the adiabatic wall temperature for this flow is estimated to be $T_{ad} \approx 430\ \textrm {K}$.

Table 2.

Free-stream conditions used in the base-flow computations.

The numerical solution of the base flow is steady and laminar, and is carried out using the CFD package CFD++^® (Chakravarthy et al. Reference Chakravarthy, Peroomian, Goldberg and Palaniswamy1998). The spatial discretization is based on a second-order upwind finite volume scheme that features a limited total variation diminishing flux interpolation together with the Harten–Lax–van Leer–Contact (HLLC) Riemann solver to minimize numerical oscillations in the vicinity of large flow gradients. For time integration, a point-implicit scheme (backward Euler) with multi-grid acceleration and residual smoothing is employed to enhance convergence to the steady state.

With respect to the base-flow boundary conditions, the primitive flow variables are prescribed at the inlet boundary and the top boundary of the domain, using values obtained from a self-similar boundary layer computation. At the centre ($z = 0$) and at the far-field ($z = z_\infty$) streamwise planes, symmetry conditions are specified. At the wall, a no-slip isothermal condition is enforced. Finally, a supersonic outflow boundary condition is specified at the outlet boundary, in which all the primitive flow variables are extrapolated from the interior of the domain.

The computational grids used to obtain the base-flow solutions are block-structured and consist of hexahedral cells. The reference grid resolution employed for each case in every direction is summarized in table 3, resulting in a total of 65 million cells for the cuboid element and 76 million cells for the ramp geometry. A coarser grid resolution used for testing grid convergence is also reported, which corresponds to a reduction of 1/4 in the number of points along each direction. As shown in figure 2, good agreement is achieved between the results obtained for the two grid resolutions.

Figure 2.

Convergence of results with respect to the base-flow grid resolution: temperature profiles at $x = 0.1\ \textrm {m}$ and at the roughness symmetry plane ($z = 0$) for $(a)$ the cuboid and $(b)$ the ramp roughness geometries; $(c)$ temporal growth rate as a function of the streamwise direction for a sinuous wake instability of wavelength $\lambda _x = 0.325\ \textrm {cm}$ (cuboid) and $\lambda _x = 0.86\ \textrm {cm}$ (ramp); shape of the recirculation bubbles induced by $(d)$ the cuboid and $(e)$ the ramp element at the roughness symmetry plane, represented as isolines of $\bar {u} = 0$.

Table 3.

Grid resolution employed for the base-flow computations. The quantities $N_x$, $N_y$ and $N_z$ denote the number of grid points along each spatial direction.

4.2. Discretization and solution of the eigenvalue problem

The numerical solution of the stability problem is performed using the von Karman Institute Extensible Stability & Transition Analysis (VESTA) toolkit, originally developed by Pinna (Reference Pinna2013). The terms that constitute the matrices of the stability equations are automatically derived and implemented with a module based on the computer algebra system Maxima. The resulting two-dimensional partial differential eigenvalue problem is assembled and discretized in both spatial directions by means of the high-order finite difference scheme with uniform error developed by Hermanns & Hernández (Reference Hermanns and Hernández2008). Similarly to the Chebyshev collocation method, this technique is based on a piecewise Lagrange interpolation constructed on a collocation grid. For a given desired polynomial order $q_p$, the coordinates of the collocation points are optimized to yield a uniform interpolation error. In this work, a polynomial order of $q_p = 8$ is employed, which is deemed to provide an appropriate trade-off between accuracy and computational cost for the problem under consideration (see for instance De Tullio et al. Reference De Tullio, Paredes, Sandham and Theofilis2013; Theiss et al. Reference Theiss, Hein, Ali and Radespiel2016).

To transform the collocation grid into the desired physical domain of study, two different mappings are used. In the wall-normal direction, the transformation introduced by Malik (Reference Malik1990) is applied, which allows half the grid points to be placed below a given coordinate. Denoting the wall-normal coordinate in the collocation grid by $\eta \in [-1,1]$, the wall-normal coordinate in the physical space is obtained as

(4.1)\begin{equation} y = \frac{y_i y_{max} (1 + \eta)}{y_{max} - \eta (y_{max} - 2 y_i)}, \end{equation}

where $y_{max}$ is the maximum wall-normal coordinate (location of the wall-normal far-field boundary) and $y_i$ denotes the coordinate at which the number of grid points is split into two halves. In the spanwise direction, a biquadratic mapping is considered (Esposito Reference Esposito2016; Groot et al. Reference Groot, Serpieri, Pinna and Kotsonis2018), given by

(4.2)\begin{equation} \left. \begin{gathered} z = \frac{a \zeta^2 + b \zeta + c}{d \zeta^2 + f},\quad a = 2 z_{max} (z_{i2} - 3 z_{i1}) - 2 z_{min} (3 z_{i2} - z_{i1}),\\ b ={-}3 (z_{i1} - z_{i2}) (z_{max} - z_{min}),\quad c = z_{i1} (3 z_{max} - z_{min}) + z_{i2} (z_{max} - 3 z_{min}),\\ d = 8 (z_{i2} - z_{i1}) - 4 (z_{max} - z_{min}),\quad f = 2 (z_{i1} - z_{i2}) + 4 (z_{max} - z_{min}), \end{gathered} \right\} \end{equation}

where $\zeta \in [-1,1]$ denotes the spanwise coordinate in the collocation grid, $z_{min}$ and $z_{max}$ are respectively the minimum and maximum spanwise coordinates of the physical domain (location of the spanwise far-field boundaries) and $z_{i,1}$, $z_{i,2}$ are the spanwise locations which divide the domain into three different regions. This transformation is a generalization of (4.1) that concentrates one-third of the grid points in each of the three regions demarcated by $z_{min} < z_{i1} < z_{i2} < z_{max}$. With these mapping relations, the grid used for solving the stability problem can be clustered towards the boundary layer and towards the roughness centreline, which is the region containing the streaks induced by the roughness element and therefore is where the strongest base-flow gradients are encountered. Even though this process transforms the stability grid into the physical domain, the mapped grid defined by the finite-difference-scheme point distribution does not coincide with the mesh employed to obtain the base-flow solution. As a consequence, before building the eigenvalue problem the base-flow data are interpolated on the stability grid by means of a cubic spline interpolation in each spatial direction.

The boundary conditions imposed on the disturbance variables are as follows. At the wall, the no-slip condition is enforced by setting the velocity perturbations to zero with a homogeneous Dirichlet condition. The same is applied for the temperature disturbance, whereas the pressure fluctuation is determined by means of a compatibility condition satisfying the wall-normal momentum equation at the wall. In the wall-normal far-field boundary, the perturbations are forced to decay by the imposition of a Dirichlet boundary condition, with the exception of $\tilde {p}$, which once again satisfies a compatibility condition. Periodic boundary conditions are specified for all perturbation variables in the spanwise direction. The value of each disturbance quantity and its first derivative normal to the boundary are constrained to be equal at the two spanwise boundaries.

For the problem under consideration, the necessary size of the eigenvalue problem to be solved results in a significant computational cost. The implicitly restarted Arnoldi method (Lehoucq & Sorensen Reference Lehoucq and Sorensen1996) is employed together with a shift–invert transformation to allow for the computation of a small number of eigenvalues and eigenvectors in a specific region of interest in the stability spectrum. Once a particular eigenvalue–eigenvector pair is known, an iterative algorithm is also used to track the evolution of the mode along different streamwise locations, wavenumbers or frequencies. This algorithm is based on a generalized Rayleigh quotient iteration (see for example Freitag & Spence Reference Freitag and Spence2007).

4.3. Computation of the temporal growth-rate decomposition terms

The computation of the terms emerging from the temporal growth-rate decomposition (3.9) requires integration along the wall-normal and spanwise directions. Here, the integration is performed using the inverse of part of the differentiation matrices of the discretization, following the methodology of Trefethen (Reference Trefethen2000). This procedure ensures consistency between the order of accuracy of integration and that of differentiation. The results obtained using this method were verified against the results obtained by means of a simple trapezoidal integration rule.

6.1. Temporal stability spectrum

Temporal 2D-LST computations were carried out along the wake induced behind each roughness element for different streamwise wavelengths ($\lambda _x = 2 {\rm \pi}/ \alpha _r$). Figures 7 and 8 show the most relevant portion of the temporal stability spectrum for each roughness shape, obtained at a streamwise location of $x = 0.1\ \textrm {m}$ and for three different streamwise wavelengths. The results for three different grid resolutions are included to illustrate grid convergence of discrete modes. For both roughness geometries, two dominant converged discrete instability modes can be identified, labelled as VCOS1 and SIN1. These modes respectively correspond to the most unstable varicose and sinuous instabilities developing in the roughness wake. The three different streamwise wavelengths chosen for each geometry cover the range of largest temporal amplification for the modes VCOS1 and SIN1, presented in figure 9 for $x = 0.1\ \textrm {m}$. Their values were selected in order to provide an adequate description of the topology of the spectrum and of the associated amplitude functions at the most relevant wavelengths for each instability. As shown in figure 9$(a)$, the varicose perturbation related to the cuboidal element features a higher growth rate than the sinuous one for most of the studied wavelengths. For the ramp-shaped roughness, however, the varicose mode dominates only for the smaller wavelengths ($\lambda _x < 0.4\ \textrm {cm}$), while the sinuous disturbance becomes the leading instability for higher ones.

Figure 7.

Most relevant region of the temporal 2D-LST spectrum for the cuboid roughness geometry at $x = 0.1\ \textrm {m}$ for three different streamwise wavelengths: $(a)$ $\lambda _x = 0.325\ \textrm {cm}$; $(b)$ $\lambda _x = 0.5\ \textrm {cm}$; $(c)$ $\lambda _x = 1\ \textrm {cm}$.

Figure 8.

Most relevant region of the temporal 2D-LST spectrum for the ramp roughness geometry at $x = 0.1\ \textrm {m}$ for three different streamwise wavelengths: $(a)$ $\lambda _x = 0.35\ \textrm {cm}$; $(b)$ $\lambda _x = 0.56\ \textrm {cm}$; $(c)$ $\lambda _x = 0.86\ \textrm {cm}$.

Figure 9.

Temporal growth-rate evolution of the most unstable wake modes for each roughness geometry as a function of the streamwise wavenumber at $x = 0.1\ \textrm {m}$: $(a)$ cuboid; $(b)$ ramp.

The two-dimensional streamwise velocity amplitude function associated to each of these disturbances is displayed in figure 10 for the cuboid element and in figure 11 for the ramp. Both instabilities mainly develop within the high-shear layer surrounding the low-velocity streak that characterizes the wake flow structure. The sinuous mode can be distinguished by its antisymmetric amplitude function with respect to the roughness mid-plane, which translates into an amplitude magnitude of zero at $z = 0$, while the varicose mode is characterized by a symmetric eigenfunction. For both roughness geometries, the biggest amplitude magnitude of the sinuous mode is concentrated in the upper regions of the low-velocity streak that are located right above the core of each counter-rotating vortex, as indicated by the streamwise vorticity contours in figure 5. According to figure 4, these regions are associated with the highest shear magnitude within the central streak. A second peak in amplitude is found at the outward side of each vortex, which coincides with the interface between the high-shear regions near the wall and the sides of the low-velocity streak, where a sudden decrease in shear magnitude is encountered. Similarly, the varicose instability features the main region of amplification in the upper central region of the streak, located above and between the vortex pair and strongly correlating with the region of high shear magnitude surrounding the central streak. Secondary amplitude peaks are once again visible at each side of the vortex cores, associated as well to the strong shear gradients established at the sides of the low-velocity streak. These observations are in good agreement with the disturbance shapes reported by other authors in similar configurations, employing both DNS and linear stability theory (see for instance Groskopf et al. Reference Groskopf, Kloker and Marxen2010a; Choudhari et al. Reference Choudhari, Li, Chang, Norris and Edwards2013; De Tullio et al. Reference De Tullio, Paredes, Sandham and Theofilis2013; De Tullio & Sandham Reference De Tullio and Sandham2015). It is also worth mentioning that, as described before, the ramp-induced central streak is taller, leading to wake instabilities that develop further away from the flat plate wall in comparison with the cuboidal case.

Figure 10.

Normalized magnitude of the streamwise velocity amplitude function for the most unstable varicose and sinuous modes (VCOS1 and SIN1) developing behind the cuboidal roughness element at $x = 0.1\ \textrm {m}$, for three different streamwise wavelengths: $(a{,}d)$ $\lambda _x = 0.325\ \textrm {cm}$; $(b{,}e)$ $\lambda _x = 0.5\ \textrm {cm}$; $(c{,}f)$ $\lambda _x = 1\ \textrm {cm}$; $(a\text {--}c)$ VCOS1 mode; $(d\text {--}f)$ SIN1 mode.

Figure 11.

Normalized magnitude of the streamwise velocity amplitude function for the most unstable varicose and sinuous modes (VCOS1 and SIN1) developing behind the ramp roughness element at $x = 0.1\ \textrm {m}$, for three different streamwise wavelengths: $(a{,}d)$ $\lambda _x = 0.35\ \textrm {cm}$; $(b{,}e)$ $\lambda _x = 0.56\ \textrm {cm}$; $(c{,}f)$ $\lambda _x = 0.86\ \textrm {cm}$; $(a\text {--}c)$ VCOS1 mode; $(d\text {--}f)$ SIN1 mode.

For the cuboid configuration at $\lambda _x = 0.325\ \textrm {cm}$, besides evolving in the three-dimensional shear layer, the varicose perturbation also grows significantly inside the low-velocity streak in a region close to the flat-plate wall. The amplitude function of this mode shows strong similarities with the mode VC (varicose-central) reported by De Tullio & Sandham (Reference De Tullio and Sandham2015) through DNS computations in a very similar set-up. This difference in the region of development of the varicose mode between the results obtained at $\lambda _x = 0.325\ \textrm {cm}$ (figure 10$a$) and those obtained for $\lambda _x = 0.5\ \textrm {cm}$ or $\lambda _x = 1\ \textrm {cm}$ (figure 10$b{,}c$) correlates with the two distinct peaks in growth rate that can be observed in figure 9$(a)$. This behaviour, already observed by Groskopf & Kloker (Reference Groskopf and Kloker2016), can be attributed to two different manifestations of the varicose-mode instability, depending on whether the instability modes evolving in the flat-plate boundary layer, which also interact with the roughness wake, behave as Mack's first- or second-mode instabilities, as discussed further below. Therefore, as introduced by Groskopf & Kloker (Reference Groskopf and Kloker2016), a distinction can also be made between varicose first- and varicose second-mode disturbances in this case.

Alongside the varicose and sinuous perturbations, several other unstable modes are found in the obtained spectra. Most of these modes are distributed in an ordered fashion along the complex plane. In fact, such instabilities belong to the family of Mack modes (first or second modes) developing in the flat-plate boundary layer, which are modulated by the presence of the roughness element. These disturbances are part of a continuous branch consisting of Mack modes representing the infinite spanwise wavenumbers ($\beta$) that are part of the 2D-LST spectrum. Only those spanwise wavenumbers that can be resolved by the discretization of the eigenvalue problem along the spanwise direction (domain size and resolution) are retrieved by the current numerical solution. The shape of this branch is different depending on the streamwise wavelength considered. The initial mode of the family is always a mode with spanwise wavenumber $\beta = 0$, denoted by M2D, followed by oblique modes with a progressively increasing spanwise wavenumber. Depending on the conditions analysed, the shape of the branch in the 2D-LST spectrum tells whether it manifests itself as Mack's first or second mode. In the case of Mack's first mode, the most unstable mode in the family is oblique for the flow conditions under study. As a result, the growth rate of the modes along the branch first increases progressively until reaching the spanwise wavenumber for maximum growth and then progressively decreases for higher values of $\beta$. This behaviour leads to the hook-shaped branch that can be observed in figures 7$(b{,}c)$ and 8$(b{,}c)$. For conditions in which the branch behaves as Mack's second mode, the most unstable mode in the family is the two-dimensional one (M2D), and the growth rate of the oblique modes progressively decreases as $\beta$ increases. The shape of the branch in this case becomes a diagonal line, similar to the one shown in figure 7$(a)$. Figure 8$(a)$ corresponds to conditions in which the branch shape is transitioning from first to second mode.

To illustrate the topology of the instabilities along the Mack-mode branch, figure 12 shows the streamwise velocity amplitude function of the two-dimensional Mack mode (M2D) and of three oblique Mack modes for the cuboid geometry at $\lambda _x = 0.5\ \textrm {cm}$ and for the ramp element at $\lambda _x = 0.56\ \textrm {cm}$. The oblique modes are respectively denoted by MO3, MO4 and MO7 according to their position along the branch (see figures 7$b$ and 8$b$). As can be observed, the amplitude function of these modes is highly modulated by the roughness wake. Far away from the wake along the $z$ direction (approximately $z > 0.005\ \textrm {m}$), the amplitude of mode M2D becomes approximately constant with respect to the spanwise direction, whereas in the wall-normal direction it adopts the shape of the classical Mack instability mode as predicted by one-dimensional local linear stability theory (LST) in a smooth flat-plate boundary layer. Therefore, the eigenfunction of mode M2D consists of a blend between an eigenfunction that evolves exclusively inside the roughness wake and the eigenfunction of a classical two-dimensional Mack mode.

Figure 12.

Normalized magnitude of the streamwise velocity amplitude function for different Mack modes developing behind each roughness geometry at $x = 0.1\ \textrm {m}$: $(a{,}b)$ two-dimensional Mack mode (M2D); $(c{,}d)$ oblique Mack mode MO3; $(e{,}f)$ oblique Mack mode MO4; $(g{,}h)$ oblique Mack mode MO7; $(a{,}c{,}e{,}g)$ $\lambda _x = 0.5\ \textrm {cm}$; $(b{,}d{,}f{,}h)$ $\lambda _x = 0.56\ \textrm {cm}$.

The oblique modes MO3, MO4 and MO7 adopt a periodic spanwise evolution outside of the roughness wake. The spanwise wavelength ($\lambda _z = 2 {\rm \pi}/ \beta$) associated to this periodicity is different for each mode. According to the periodic boundary conditions enforced in the eigenvalue problem and the spanwise symmetry of the base flow, the values of $\lambda _z$ that can be resolved by the discretization are constrained by the distance between the edge of the roughness wake (in this case located approximately at $z = {\pm }0.002\ \textrm {m}$ for the cuboid geometry and at $z = \pm 0.001\ \textrm {m}$ for the ramp; see figure 4$b{,}e$) and the spanwise boundary of the domain. For instance, mode MO3 adopts a spanwise wavelength that is precisely equal to the distance between the edge of the roughness wake and the spanwise boundary. Therefore, it contains one spanwise period within each half of the span of the computational domain. Similarly, the next mode (MO4) adopts a spanwise wavelength that is four-fifths of the distance between the wake edge and the spanwise boundary, thus allowing one and one-fourth of a period inside each half of the domain span. Finally, mode MO7, which is the most unstable Mack mode found in the spectrum at the current conditions for both obstacle shapes, contains two spanwise periods within each half of the domain width, with a spanwise wavelength equal to half of the allowed distance. As in the case of modes MO3 and MO4, each of the oblique modes resolved along the Mack-mode branch differs in one-fourth of a period with respect to the previous one. For example, although not shown here, mode MO5 contains one and a half periods, mode MO6 contains one and three-fourths of a period, and so on.

The remaining discrete unstable modes found in the spectra are also sinuous and varicose deformations of the low-velocity streak, presenting different peak–valley regions of amplitude development located around the three-dimensional shear layer. For the cases investigated in this work, these modes are found to grow at a smaller rate than the leading wake instabilities (VCOS1 and SIN1). The shape of the amplitude function for this kind of disturbance is discussed by Padilla Montero & Pinna (Reference Padilla Montero and Pinna2020). The group of modes distributed along the real axis in the left part of the spectrum of figures 7$(c)$ and 8$(c)$ belong to the slow acoustic continuous branch.

6.1.1. Convergence of the growth rate with respect to grid resolution

A grid convergence analysis on the growth rate of modes SIN1, VCOS1, M2D and MO7 for each roughness geometry is presented in figure 13 for a streamwise station at $x = 0.1\ \textrm {m}$. This figure shows contour plots of the relative error in the temporal growth rate of each mode as a function of the grid resolution in each direction. The relative error ($\varepsilon _r$) for each mode is evaluated with respect to the growth rate obtained at the finest grid resolution, which corresponds to $N_z = N_y = 201$ points, so that $\varepsilon _{r} = |\omega _i - \omega _{i,201\times 201}| / |\omega _{i,201\times 201}|$. The results indicate that all instabilities are more sensitive to the spanwise grid resolution ($N_z$), owing to the fact that the amplitude functions undergo changes over a longer distance along the spanwise direction than along the wall-normal direction. For the computational domain employed for each roughness geometry, a larger grid resolution is required in the spanwise direction to achieve the same degree of convergence as in the wall-normal direction.

Figure 13.

Contours of relative error in the temporal growth rate of different unstable modes for each roughness geometry at $x = 0.1\ \textrm {m}$ as a function of the grid resolution in each direction: $(a\text {--}d)$ cuboid, $\lambda _x = 0.5\ \textrm {cm}$; $(e\text {--}h)$ ramp, $\lambda _x = 0.56\ \textrm {cm}$; $(a{,}e)$ mode SIN1; $(b{,}f)$ mode VCOS1; $(c{,}g)$ mode M2D; $(d{,}h)$ mode MO7. The relative error is evaluated with respect to the growth rate obtained with the finest grid resolution ($N_z = N_y = 201$) for each mode: $\varepsilon _{r} = |\omega _i - \omega _{i,201\times 201}| / |\omega _{i,201\times 201}|$. The star symbol indicates the position of the reference grid resolution employed for the results presented in this work ($N_z = 181$, $N_y = 141$ points).

Besides ensuring a fine grid resolution in the roughness wake, the discretization of oblique Mack modes modulated by the roughness wake also requires enough spanwise grid points within the boundary layer at the sides of the roughness element to resolve the associated desired spanwise wavenumbers. In this study, $N_z = 181$ is chosen as a reference spanwise resolution, since it provides a good degree of convergence for the wake instabilities as well as for a substantial number of oblique Mack modes, making it possible to describe the most relevant region of the continuous Mack-mode branch. Along the wall-normal direction, a reference resolution of $N_y = 141$ points is selected, as it provides an appropriate trade-off between convergence and computational cost for the different cases investigated.

6.1.2. Effect of the domain size on the stability spectrum

To emphasize the continuous nature of the Mack-mode branch obtained in the stability spectra, figure 14 illustrates the effect of changing the spanwise size of the computational domain employed for the stability analysis ($z_{max}$) on the results obtained for the cuboid configuration at $\lambda _x = 0.5\ \textrm {cm}$ and for the ramp geometry at $\lambda _x = 0.56\ \textrm {cm}$. Since the regions of development of the wake instability modes in these conditions are contained within the roughness wake, the effect of changing $z_{max}$ is negligible for these modes. On the other hand, a change in the spanwise domain size modifies the spanwise wavenumbers that can be resolved by the discretization. As a result, for each particular value of $z_{max}$, different oblique Mack modes are obtained in the spectrum. Nevertheless, the shape of the continuous Mack-mode branch must remain the same independently of the spanwise domain size considered. In other words, different values of $z_{max}$ resolve different sets of modes from the infinite number that compose the branch. This behaviour is demonstrated by the spectra presented in figure 14.

Figure 14.

Most relevant region of the temporal stability spectrum for each roughness shape at $x = 0.1\ \textrm {m}$ for three different spanwise domain sizes: $(a)$ cuboid, $\lambda _x = 0.5\ \textrm {cm}$; $(b)$ ramp, $\lambda _x = 0.56\ \textrm {cm}$.

Figure 15 represents the sensitivity of the spectrum to the wall-normal size of the computational domain employed for the stability calculations ($y_{max}$), illustrated for three different wall-normal domain sizes at $x = 0.1\ \textrm {m}$. As can be observed, the influence of the upper boundary of the domain on the discrete modes of the stability spectrum is negligible. The only sensitive modes are those located at the real axis in figure 15$(b)$, which correspond to the slow acoustic continuous branch, and therefore their eigenfunctions do not decay to zero towards the upper boundary of the domain. The value $y_{max} = 300l = 46.50h$ is the actual wall-normal size employed in the results shown throughout this work for $x = 0.1\ \textrm {m}$.

Figure 15.

Most relevant region of the temporal stability spectrum at $x = 0.1\ \textrm {m}$ for three different wall-normal domain sizes ($y_{max}$): $(a)$ cuboid, $\lambda _x = 0.5\ \textrm {cm}$; $(b)$ ramp, $\lambda _x = 0.56\ \textrm {cm}$.

6.2. Streamwise evolution of the growth rate of wake instabilities

The streamwise evolution of the two leading instabilities along the wake induced by each roughness geometry is assessed next. Figure 16 displays the evolution of the temporal growth rate of the SIN1 and VCOS1 modes for five different wavelengths that span the range of amplification shown in figure 9. The first streamwise location considered for all cases corresponds to the end of the downstream recirculation bubble induced at the trailing edge of each roughness element (see figures 2 and 3), located at $x_{b,TE} = 0.0667\ \textrm {m}$ and $x_{b,TE} = 0.0633\ \textrm {m}$ for the cuboid and ramp elements, respectively. The growth rate of both wake instabilities increases very rapidly immediately behind the trailing edge recirculation bubble, reaching its maximum in a very short distance downstream, after which it starts to decrease progressively. The distance at which the maximum growth rate is attained depends on the streamwise wavelength. In the vicinity of the cuboidal roughness trailing edge, the magnitude of the growth rate is significantly higher for the varicose mode than for the sinuous one for most of the streamwise wavelengths considered. By contrast, near the ramp roughness element the sinuous instability reaches higher amplification rates than the varicose one for $\lambda _x \geq 0.56\ \textrm {cm}$. Specifically, the maximum growth rate attained by mode SIN1 at the illustrated wavelengths is $\omega _i \approx 7.2 \times 10^4\ \textrm {rad}\,\textrm {s}^{-1}$ for the ramp configuration and $\omega _i \approx 4.3 \times 10^4\ \textrm {rad}\,\textrm {s}^{-1}$ for the cuboid configuration. In the case of the varicose disturbance, however, mode VCOS1 reaches a maximum of $\omega _i \approx 7.2 \times 10^4\ \textrm {rad}\,\textrm {s}^{-1}$ for the cuboid geometry and $\omega _i \approx 5.5 \times 10^4\ \textrm {rad}\,\textrm {s}^{-1}$ for the ramp case. This constitutes a 67 % increase in the peak growth rate achieved by the ramp-induced sinuous disturbance with respect to that excited by the cuboid, and a 31 % increase in the maximum growth rate attained by the cuboid-induced varicose instability with respect to the ramp counterpart. Of particular interest is the evolution of the varicose mode for the cuboid geometry at $\lambda _x = 0.28\ \textrm {cm}$. As shown in figure 9$(a)$, this wavelength is located within the region of development of the varicose ‘second’ mode. This range of wavelengths exhibits the highest temporal amplification for the varicose disturbance evolving in the cuboid wake.

Figure 16.

Temporal growth-rate evolution of the most unstable wake modes for each roughness geometry as a function of the streamwise coordinate for five different wavelengths: $(a{,}c)$ SIN1 mode; $(b{,}d)$ VCOS1 mode; $(a{,}b)$ cuboid; $(c{,}d)$ ramp.

The rate of decay in growth rate after the location of maximum amplification is also found to be significantly dependent on the streamwise wavelength considered. For the cuboid configuration, the growth rate of both instabilities for $\lambda _x \leq 0.325\ \textrm {cm}$ decreases continuously until the end of the domain at very similar rates. For these wavelengths, the varicose mode is the leading instability along all of the domain investigated. On the other hand, for $\lambda _x \geq 0.5\ \textrm {cm}$, the amplification rate approaches a nearly constant value after a given distance downstream. This distance is found to decrease as the streamwise wavelength increases, and it is found to be significantly shorter for the sinuous mode than for the varicose one. This makes the sinuous disturbance the dominant wake instability for approximately $x > 0.2\ \textrm {m}$ at this range of wavelengths.

Regarding the ramp roughness wake, the streamwise evolution of the sinuous perturbation for $\lambda _x \geq 0.56\ \textrm {cm}$ undergoes first a rapid decay followed by a slower monotonic decrease downstream, in all cases with very similar slopes. For $\lambda _x = 0.35\ \textrm {cm}$ and $\lambda _x = 0.42\ \textrm {cm}$, however, a deviation in this growth-rate behaviour is encountered starting respectively at $x = 0.15\ \textrm {m}$ and at $x = 0.18\ \textrm {m}$, where the disturbance shows an important destabilization. For the varicose perturbation developing in the wake behind the ramp geometry, a rapid decay is also found in the first portion of the streamwise domain under consideration for all the studied wavelengths. However, in this case the rate of decrease is more pronounced than for the sinuous perturbation and extends over a longer streamwise distance. This strong stabilization brings the growth rate of mode VCOS1 below $10^4\ \textrm {rad}\,\textrm {s}^{-1}$. Fluctuations in the growth rate of this mode are also found, in this case for $\lambda _x = 0.35\ \textrm {cm}$ and $\lambda _x = 0.56\ \textrm {cm}$. These observations, together with the results for the cuboidal geometry described before, suggest that there is a mechanism in the roughness wake which, for a certain range of streamwise wavelengths, energizes the wake instabilities, preventing their growth rate from decaying monotonically downstream as would be expected from the continuous dissipation of the base-flow quantities displayed in figure 6. A plausible explanation for this behaviour would be a continuous synchronization between the wake instability modes and the boundary-layer modes that interact with the roughness wake. Further discussion of this phenomenon is provided in §§ 6.3 and 6.4.

For the ramp-induced varicose mode at $\lambda _x = 0.56\ \textrm {cm}$, the rise in the amplification rate that can be observed towards the end of the domain is very similar to the growth-rate evolution inherent to Mack's second-mode instability. Tracking the streamwise evolution of this wake mode shows that it progressively approaches the continuous branch associated to the Mack-mode family and becomes absorbed by it after a given distance downstream. At this location, the amplitude function of mode VCOS1 strongly resembles those of the other modes comprised in this branch, and it becomes difficult to distinguish its evolution from that of the other Mack-mode instabilities by looking at the computed 2D-LST spectra. The curve shown in figure 16$(d)$ corresponds to the evolution obtained by means of a Rayleigh quotient iteration. Two possible scenarios can be devised to explain the rise in growth rate obtained for $x > 0.33\ \textrm {m}$ at $\lambda _x = 0.56\ \textrm {cm}$. On the one hand, as the varicose mode becomes close to the Mack-mode branch, the mode tracked by the algorithm could possibly switch to a different one, in this case corresponding to a particular Mack-mode disturbance. On the other hand, since the varicose-mode disturbance becomes absorbed by the Mack-mode branch in these conditions, another possibility is that the behaviour of the varicose mode is so strongly affected by an interaction with the boundary-layer instability modes that its evolution becomes very similar to that of Mack's second mode. The results of the temporal growth-rate decomposition presented in § 6.3 give more support to the second scenario.

6.2.1. Integrated amplification factors obtained by means of the Gaster transformation

For the roughness Reynolds number considered in this work, the instabilities growing in the roughness wake have a purely convective nature (De Tullio & Sandham Reference De Tullio and Sandham2015; Groskopf & Kloker Reference Groskopf and Kloker2016). Therefore, it is appropriate to track the evolution of their growth in space in order to assess which is the dominant disturbance excited by each geometry in terms of the integrated amplification factor ($N$-factor), given by $N = -\int \alpha _i\,\textrm {d} x$. To avoid the significant computational cost associated with solving the spatial eigenvalue problem (2.11), the relation between the temporal and the spatial growth rates given by Gaster (Reference Gaster1962) is employed. This relation, commonly known as the Gaster transformation, can be expressed as follows:

(6.1)\begin{equation} \alpha_i |_{\omega_i = 0} ={-} \frac{\omega_i |_{\alpha_i = 0}}{c_g} + O(\omega_i |_{\alpha_i = 0}^2), \end{equation}

where $c_g = \partial \omega _r / \partial \alpha _r$ denotes the group velocity of the disturbance. Although, in general, (6.1) is only accurate for convective instabilities with small temporal growth rates, it has been successfully applied by previous authors to the analysis of roughness-induced instabilities at high speed (see Groskopf et al. Reference Groskopf, Kloker, Stephani, Marxen and Iaccarino2010b; Groskopf & Kloker Reference Groskopf and Kloker2016). To verify its applicability to the temporal results presented here, figure 17 shows a comparison between spatial growth rates obtained by means of the Gaster transformation and spatial growth rates computed by solving the spatial eigenvalue problem for the leading wake instabilities. Good overall agreement is found for both roughness geometries, which renders the use of (6.1) a legitimate assumption for the purpose of this analysis. The biggest discrepancies encountered are for the cuboid varicose mode at $\lambda _x = 0.28\ \textrm {cm}$, corresponding to its particular region of development associated to small streamwise wavenumbers (high frequencies), as shown in figure 9$(a)$.

Figure 17.

Spatial growth-rate evolution of the most unstable wake modes obtained by means of the Gaster transformation (6.1) (lines) and by solving the spatial eigenvalue problem (symbols): $(a)$ cuboid; $(b)$ ramp.

Temporal stability computations make it possible to construct a stability diagram of the form $\omega _i = \omega _i(\alpha _r,x)$, or $\omega _i = \omega _i(\lambda _x,x)$. By means of the Gaster transformation, the function $\alpha _i = \alpha _i(\lambda _x,x)$ can then be obtained. However, a meaningful evaluation of the integrated amplification factors requires a surface of the type $\alpha _i = \alpha _i(f,x)$, where $f = \omega _r/ (2 {\rm \pi})$ is the frequency of the disturbance, so that curves of $\alpha _i$ as a function of $x$ for a constant frequency can be integrated to obtain the evolution of $N$ associated to that frequency. Since $\omega _r$ is a solution of the temporal eigenvalue problem, it cannot be fixed constant a priori from the temporal results. Therefore, to obtain $\alpha _i = \alpha _i(f,x)$, an interpolation of the stability diagram $\alpha _i = \alpha _i(\lambda _x,x)$ must be employed. From the temporal stability calculations, the resulting frequency of the disturbances is found to be approximately proportional to $\alpha _r$. In addition, in the unstable region of the stability diagram, the range of frequencies spanned by each of the $\alpha _i = \alpha _i(x)$ curves for a constant $\lambda _x$ is very small in comparison with the range of frequencies over which the instabilities grow. These conditions favour the accuracy of the interpolation process. In this study, this procedure is performed by means of a cubic Hermite interpolation.

The spatial stability diagrams ($\alpha _i = \alpha _i(f,x)$) obtained by means of the Gaster transformation for modes SIN1 and VCOS1 and for each of the roughness geometries under study are displayed in figure 18$(c$–$f)$. The shaded contours represent the spatial growth rate, and the contour lines superimposed illustrate the associated $N$-factor evolution. For comparison purposes, the diagrams shown in figure 18$(a{,}b)$ correspond to the instabilities developing in the smooth flat-plate boundary layer, which are obtained by means of LST for a given spanwise wavenumber. Figure 18$(a)$ depicts contours for the two-dimensional Mack-mode instability ($\beta = 0$), which is fully governed by high frequencies corresponding to second mode behaviour, whereas figure 18$(b)$ shows contours for an oblique Mack-mode instability with a spanwise wavenumber $\lambda _z = 0.4\ \textrm {cm}$, for which both Mack's first- and second-mode regions of development can be observed, respectively associated to the low- and high-frequency ranges. The integration of $N$ is started at the first streamwise station at which the modes become unstable. In the case of the wake instabilities, it is found that for certain frequencies the disturbances are already unstable at the first streamwise location analysed. Under these circumstances, a value of $N = 0$ is assumed at that location. This is a common procedure employed by other authors in similar analyses; see for instance Choudhari et al. (Reference Choudhari, Li, Chang, Norris and Edwards2013), De Tullio et al. (Reference De Tullio, Paredes, Sandham and Theofilis2013) and Theiss et al. (Reference Theiss, Hein, Ali and Radespiel2016). In addition, figure 19$(a)$ presents the streamwise evolution of the $N$-factor envelope ($N_{max}$) curves for each of the wake modes investigated. Here, $N_{max}$ denotes the maximum $N$-factor obtained at each streamwise position for all the different frequencies analysed. The particular frequency associated with $N_{max}$ at each location is shown in figure 19$(b)$.

Figure 18.

Shaded contours of spatial growth rate ($-\alpha _i$) and contour lines of $N$-factor as a function of frequency and streamwise position for $(a)$ two-dimensional Mack mode evolving in the smooth flat-plate boundary layer; $(b)$ oblique Mack mode ($\lambda _z = 0.4\ \textrm {cm}$) evolving in the smooth flat-plate boundary layer; $(c)$ cuboid, SIN1 mode; $(d)$ cuboid, VCOS1 mode; $(e)$ ramp, SIN1 mode; $(f)$ ramp, VCOS1 mode. Results for $(c$–$f)$ were obtained by means of the Gaster transformation (6.1).

Figure 19.

$(a)$ Streamwise evolution of the $N$-factor envelope ($N_{max}$) of the most unstable wake modes for each roughness geometry, obtained by employing the Gaster-transformed spatial growth rates; $(b)$ frequency associated to $N_{max}$ at each streamwise location for each of the curves shown in $(a)$.

These results complete the observations introduced before based on the temporal growth-rate evolution. The $N$-factor envelopes confirm that the cuboid geometry favours the development of a varicose wake instability over a sinuous one, whereas the ramp configuration excites a sinuous disturbance more effectively than a varicose one. Close to the roughness element, the cuboidal geometry is more effective in promoting the linear growth of the dominant wake instability than the ramp-shaped one. However, the ramp element leads to a higher disturbance amplitude further downstream ($x > 0.2\ \textrm {m}$).

As shown in figure 18$(d)$, two distinct regions of development for the cuboid-induced varicose instability can be identified. The frequency range for each of these regions is found to correlate well with that of Mack's first- and second-mode instabilities developing in the smooth flat plate. Although the biggest growth rates for this disturbance are located within the high-frequency band, the decay in this range is more accentuated than for the lower frequencies. As a result, the frequencies that govern the $N$-factor envelope for this mode are at first located in the upper range ($f \approx 280\ \textrm {kHz}$) and further downstream in the lower band ($f \approx 160\ \textrm {kHz}$), as shown in figure 19$(b)$. This leads to the abrupt change in the slope of the envelope curve at $x \approx 0.28\ \textrm {m}$. In accordance with the temporal stability results described above, this is the only instability for which two such different regions of excitation can be clearly observed in the stability diagram.

For the ramp-induced sinuous disturbance (see figure 18$e$), the region of maximum growth is concentrated in the range of $f \approx 100\ \textrm {kHz}$ to $f \approx 160\ \textrm {kHz}$. This frequency band is associated with Mack's first-mode instabilities for the entire domain length. For higher frequencies ($f \approx 230\ \textrm {kHz}$ to $f \approx 270\ \textrm {kHz}$), a small enhancement of the spatial growth rates is observed starting at $x \approx 0.15\ \textrm {m}$. The same phenomenon is visible in the temporal growth rate curves shown in figure 16$(c)$ for $\lambda _x = 0.35\ \textrm {cm}$ and $\lambda _x = 0.42\ \textrm {cm}$. This behaviour correlates with the region of the stability diagram that overlaps with the development of Mack's second-mode instability in the flat-plate boundary layer. This observation further supports the argument of an interaction between the wake modes and the boundary-layer disturbances. In this case, however, this interaction does not have an impact on the $N$-factor envelope of the wake instability.

Although less relevant in terms of $N_{max}$, the cuboid-induced sinuous disturbance and the ramp-induced varicose perturbation also show an influence from the boundary-layer instabilities. As in the ramp configuration, the sinuous mode associated to the cuboid geometry also features its maximum $N$ in the low-frequency band. For $x > 0.2\ \textrm {m}$, its spatial growth becomes approximately constant with respect to the streamwise position for $f < 130\ \textrm {kHz}$ (see figure 18$c$). This feature can also be observed in the temporal growth curves shown in figure 16$(a)$ for $\lambda _x = 0.7\ \textrm {cm}$ and $\lambda _x = 1\ \textrm {cm}$. This reduction in the decay of the growth rate of the instability is attributed to an interaction between the sinuous disturbance and Mack's first-mode instability. It has a significant impact on the evolution of the associated $N$-factor envelope, which shows a clear increase in slope in this region. Mode VCOS1 excited by the ramp roughness element is found to be rapidly stabilized at high frequencies, as shown in figure 18$(d)$. However, its growth rate is found to be influenced by both the first and second boundary-layer disturbances in the range $f \approx 100\ \textrm {kHz}$ to $f \approx 200\ \textrm {kHz}$. The interaction with Mack's first mode is noticeable starting at $x \approx 0.18\ \textrm {m}$, and the interaction with Mack's second mode is noticeable starting at $x \approx 0.33\ \textrm {m}$ (see also figure 16$d$).

6.3. Temporal growth-rate decomposition of the wake instabilities

This section describes the application of the temporal growth-rate decomposition introduced in § 3.1 to the evolution of the leading wake instabilities (modes VCOS1 and SIN1) excited by each roughness shape. Although the decomposition given by (3.9) consists of a large number of terms, only a small number of them are relevant for the physical interpretation of the problem studied. To illustrate which are the most relevant contributions, figures 20 and 21 present the decomposition obtained at $x = 0.1\ \textrm {m}$ for the streamwise wavelengths considered in figures 7 and 8 for the cuboid and ramp configurations. The terms denoted by $\hat {P}_{s,2}$ and $\hat {P}_{s,3}$ respectively refer to the wall-normal and the spanwise contributions that constitute $\hat {P}_s$, i.e.

(6.2a,b)\begin{equation} \hat{P}_{s,2} ={-} \int_{z^-}^{z^+} \int_{y^-}^{y^+} \bar{\rho} \langle \hat{v}, \hat{s} \rangle \dfrac{\partial \bar{T}}{\partial y}\, \textrm{d} y\, \textrm{d} z,\quad \hat{P}_{s,3} ={-} \int_{z^-}^{z^+} \int_{y^-}^{y^+} \bar{\rho} \langle \hat{w}, \hat{s} \rangle \dfrac{\partial \bar{T}}{\partial z}\, \textrm{d} y\, \textrm{d} z, \end{equation}

and the terms $\hat {P}_{RS,12}$ and $\hat {P}_{RS,13}$ are the parts of the total Reynolds stress production term $\hat {P}_{RS}$ that respectively correspond to the stress components given by $\bar {\rho } \tilde {u} \tilde {v}$ and by $\bar {\rho } \tilde {u} \tilde {w}$, i.e.

(6.3a,b)\begin{align} \hat{P}_{RS,12} ={-}\int_{z^-}^{z^+} \int_{y^-}^{y^+} \bar{\rho} \langle \hat{u}, \hat{v} \rangle \dfrac{\partial \bar{u}}{\partial y}\, \textrm{d} y\, \textrm{d} z,\quad \hat{P}_{RS,13} ={-} \int_{z^-}^{z^+} \int_{y^-}^{y^+} \bar{\rho} \langle \hat{u}, \hat{w} \rangle \dfrac{\partial \bar{u}}{\partial z}\, \textrm{d} y\, \textrm{d} z. \end{align}

Figure 20.

Values of the most important contributions to the temporal growth-rate decomposition (3.9) for modes SIN1 and VCOS1 at $x = 0.1\ \textrm {m}$ in the cuboid wake: $(a)$ $\lambda _x = 0.325\ \textrm {cm}$; $(b)$ $\lambda _x = 0.5\ \textrm {cm}$; $(c)$ $\lambda _x = 1\ \textrm {cm}$. Each contribution is divided by the factor $2 \hat {E}$ (see (3.9)).

Figure 21.

Values of the most important contributions to the temporal growth-rate decomposition (3.9) for modes SIN1 and VCOS1 at $x = 0.1\ \textrm {m}$ in the ramp wake: $(a)$ $\lambda _x = 0.35\ \textrm {cm}$; $(b)$ $\lambda _x = 0.56\ \textrm {cm}$; $(c)$ $\lambda _x = 0.86\ \textrm {cm}$. Each contribution is divided by the factor $2 \hat {E}$ (see (3.9)).

Only those terms which have a significant magnitude are included in figures 20 and 21; the contributions of the rest are found to be unimportant (of the order of $\hat {P}_{Q}$ or below) in this case. The biggest contributions to the disturbance kinetic energy for both wake modes come from the production terms due to the relevant disturbance Reynolds stresses ($\hat {P}_{RS,12}$ and $\hat {P}_{RS,13}$), whereas the most important contributions to the potential energy of the disturbance are due to the transport of the entropy perturbation ($\hat {P}_{s,2}$ and $\hat {P}_{s,3}$). The production term due to the disturbance heat source ($\hat {P}_{Q}$) is in all cases one order of magnitude smaller than the $\hat {P}_s$ and the $\hat {P}_{RS}$ production terms, which is in good agreement with the findings of Chu (Reference Chu1965). It is worth emphasizing that, to the best of the authors’ knowledge, previous analyses found in the literature for isolated roughness elements in high-speed flows have only examined the evolution of the production of kinetic energy due to the disturbance Reynolds stresses. However, the current results show that the energy production terms associated to the transport of entropy fluctuations constitute an important contribution to the potential energy of the disturbance. Therefore, although their relative magnitude with respect to the production of kinetic energy due to the disturbance Reynolds stresses cannot be concluded from this formulation of the perturbation energy (see the reasons stated at the beginning of § 3), these results suggest that $\hat {P}_s$ might also play a significant role in high-speed flows. It is important to note that, within the definition of disturbance energy given by Chu (Reference Chu1965), the largest destabilizing (positive) contributions to the temporal growth rate of the wake modes are $\hat {P}_{s,2}$ and $\hat {P}_{s,3}$, while the largest stabilizing (negative) contribution is the dissipation of disturbance energy due to thermal conduction ($\hat {D}_k$). Nevertheless, since the terms $\hat {P}_{s,2}$, $\hat {P}_{s,3}$ and $\hat {D}_k$ originate from a different governing equation than the terms $\hat {P}_{RS,12}$, $\hat {P}_{RS,13}$ and $\hat {D}_\mu$, it is not possible to conclude which contributions are actually larger in practice, since different choices of the multipliers that constitute the disturbance energy formulation alter the relative magnitudes of these sets of terms.

For the cuboid geometry, the decomposition reveals a difference in the source of disturbance energy for the varicose mode at $x = 0.1\ \textrm {m}$ and $\lambda _x = 0.325\ \textrm {cm}$ with respect to the other two wavelengths. As can be noticed, at this wavelength, the growth rate of mode VCOS1 is dominated by the wall-normal production terms ($\hat {P}_{s,2}$ and $\hat {P}_{RS,12}$), whereas for $\lambda _x = 0.5\ \textrm {cm}$ and $\lambda _x = 1\ \textrm {cm}$ it is governed by the spanwise production contributions ($\hat {P}_{s,3}$ and $\hat {P}_{RS,13}$). This reflects once again the occurrence of two different manifestations of the varicose disturbance depending on the streamwise wavelength (or frequency in the spatial case), as discussed previously.

6.3.1. Streamwise evolution of the temporal growth-rate decomposition

With the aim of investigating the particular features found in the growth-rate evolution of the SIN1 and VCOS1 modes, the streamwise evolution of the temporal growth-rate decomposition for each of the two disturbances is presented in figures 22 and 23 for the cuboid configuration and in figures 24 and 25 for the ramp. For the cuboidal geometry, it is observed that immediately behind the trailing-edge recirculation bubble, the energy of both instabilities is governed by the wall-normal production terms $\hat {P}_{s,2}$ and $\hat {P}_{RS,12}$. This implies that in this region the modes extract their energy from the wall-normal base-flow gradients. A short distance further downstream, depending on the streamwise wavelength, the wall-normal production terms rapidly decrease and the disturbance energy is transferred to the spanwise production contributions $\hat {P}_{s,3}$ and $\hat {P}_{RS,13}$, bringing the wall-normal terms to a negative (stabilizing) energy contribution for most of the studied wavelengths. After this first exchange of energy from the wall-normal to the spanwise contributions, the leading term $\hat {P}_{s,3}$ reaches a maximum and begins to decrease at an approximately constant rate. In parallel, the wall-normal terms continue to decrease, which enhances their stabilizing effect. However, after a given distance downstream, the wall-normal production terms start to rise again in most of the cases. This rise is accompanied by an acceleration in the rate of decrease of the spanwise terms, so that there is reciprocal evolution between $\hat {P}_{s,2}$ and $\hat {P}_{s,3}$ and between $\hat {P}_{RS,12}$ and $\hat {P}_{RS,13}$. This implies that another energy transfer begins to take place along the roughness wake, in this case from spanwise to wall-normal production terms, at a slower pace than the first one. This second exchange can be observed for the sinuous mode for all three wavelengths: for $\lambda _x = 0.325\ \textrm {cm}$ and $\lambda _x = 0.5\ \textrm {cm}$, it starts at $x \approx 0.18\ \textrm {m}$ and continues until the end of the domain; for $\lambda _x = 1\ \textrm {cm}$, it starts immediately after the first exchange ($x \approx 0.13\ \textrm {m}$). In the case of the varicose mode, it can also be observed for $\lambda _x = 0.5\ \textrm {cm}$ (starting at $x \approx 0.35\ \textrm {m}$) and for $\lambda _x = 1\ \textrm {cm}$ (starting at $x \approx 0.14\ \textrm {m}$) within the studied domain length.

Figure 22.

Streamwise evolution of the main production $(a\text {--}c)$ and dissipation $(d\text {--}f)$ terms for mode SIN1 developing in the cuboidal roughness wake at three different wavelengths: $(a{,}d)$ $\lambda _x = 0.325\ \textrm {cm}$; $(b{,}e)$ $\lambda _x = 0.5\ \textrm {cm}$; $(c{,}f)$ $\lambda _x = 1\ \textrm {cm}$.

Figure 23.

Streamwise evolution of the main production $(a\text {--}c)$ and dissipation $(d\text {--}f)$ terms for mode VCOS1 developing in the cuboidal roughness wake at three different wavelengths: $(a{,}d)$ $\lambda _x = 0.325\ \textrm {cm}$; $(b{,}e)$ $\lambda _x = 0.5\ \textrm {cm}$; $(c{,}f)$ $\lambda _x = 1\ \textrm {cm}$.

Figure 24.

Streamwise evolution of the main production $(a\text {--}c)$ and dissipation $(d\text {--}f)$ terms for mode SIN1 corresponding to the ramp geometry at three different wavelengths: $(a{,}d)$ $\lambda _x = 0.35\ \textrm {cm}$; $(b{,}e)$ $\lambda _x = 0.56\ \textrm {cm}$; $(c{,}f)$ $\lambda _x = 0.86\ \textrm {cm}$.

Figure 25.

Streamwise evolution of the main production $(a\text {--}c)$ and dissipation $(d\text {--}f)$ terms for mode VCOS1 corresponding to the ramp geometry at three different wavelengths: $(a{,}d)$ $\lambda _x = 0.35\ \textrm {cm}$; $(b{,}e)$ $\lambda _x = 0.56\ \textrm {cm}$; $(c{,}f)$ $\lambda _x = 0.86\ \textrm {cm}$.

From the energy signature of both cuboid-induced modes at $\lambda _x = 1\ \textrm {cm}$, it is clear that this second energy transfer mechanism brings the spanwise production terms asymptotically to approximately zero. Hence, it does not lead to stabilizing contributions associated to the work done through the spanwise base-flow gradients. This implies that the wake instabilities in this condition become governed only by the energy production processes linked to the wall-normal base-flow gradients. This is an inherent characteristic of the energy budget of boundary-layer modes that develop in a smooth flat-plate boundary layer in which no spanwise base-flow gradients exist. In accordance with this observation, it has been found that the occurrence of the second energy exchange correlates with streamwise stations at which the amplitude of the eigenfunctions associated with the wake instabilities VCOS1 and SIN1 begins to increase also in the boundary layer outside of the roughness wake, with amplitude distributions that are equivalent to the Mack-mode instabilities shown in figure 12. This provides evidence that there is an interaction between the wake instabilities and the boundary-layer modes which influences the growth-rate evolution of the wake disturbances. As described in § 6.4, this interaction can be interpreted as a continuous synchronization between the wake modes and the boundary-layer modes.

For the ramp-induced instabilities, the disturbance energy near the roughness trailing edge is also governed by the wall-normal energy production quantities ($\hat {P}_{s,2}$ and $\hat {P}_{RS,12}$). A clear difference between the energy signature of mode SIN1 for both roughness geometries is that, for the ramp element, the first transfer of energy from the wall-normal to the spanwise production terms takes place over a longer streamwise distance. This is particularly evident for the cases with $\lambda _x = 0.56\ \textrm {cm}$ and $\lambda _x = 0.86\ \textrm {cm}$, for which the spanwise term $\hat {P}_{s,3}$ does not become the dominant contribution until approximately halfway into the domain length. Therefore, the wall-normal base-flow gradients remain responsible for the excitation of the leading disturbance in the ramp-induced wake for a long distance downstream. Additionally, for $\lambda _x = 0.35\ \textrm {cm}$, the disturbance energy evolution confirms that the sinuous mode is strongly influenced by the presence of Mack-mode instabilities starting at $x \approx 0.15\ \textrm {m}$, as evidenced by the sharp rise in the term $\hat {P}_{s,2}$ and the sharp decrease in the term $\hat {P}_{s,3}$. This explains the rise in growth rate that can be observed in figure 16$(c)$ for this wavelength. In this case, however, the effect of the synchronization does not last until the end of the domain. Downstream of $x \approx 0.2\ \textrm {m}$, a third energy exchange takes place and the mode recovers the behaviour of a pure wake instability, eventually becoming governed by the spanwise base-flow gradients in a similar way to the other two cases, for which no sign of a synchronization is present.

Significant differences are also found between the energy balances of the varicose mode in the two roughness configurations. For the ramp geometry, the first energy transfer between wall-normal and spanwise terms occurs at a much shorter distance behind the roughness trailing edge. Similarly, the second energy exchange, which leads to a transformation of the instability towards a hybrid wake–boundary-layer disturbance, takes place further upstream and is encountered within the domain length for all three streamwise wavelengths analysed. For $\lambda _x = 0.35\ \textrm {cm}$, it begins to occur at $x \approx 0.17\ \textrm {m}$, for $\lambda _x = 0.56\ \textrm {cm}$ at $x \approx 0.12\ \textrm {m}$ and for $\lambda _x = 0.86\ \textrm {cm}$ at $x \approx 0.09\ \textrm {m}$. This explains the fluctuations in the growth rate encountered for $\lambda _x = 0.35\ \textrm {cm}$ and $\lambda _x = 0.56\ \textrm {cm}$, as shown in figure 16$(d)$. Furthermore, since for $\lambda _x = 0.86\ \textrm {cm}$ the synchronization with the boundary-layer modes begins upstream of $x = 0.1\ \textrm {m}$, the amplitude function of the varicose mode at these conditions (displayed in figure 11$c$) already shows a non-zero amplitude in the boundary layer at the sides of the roughness wake.

The contribution of the disturbance energy dissipation terms follows a similar trend for all the different cases analysed for each roughness geometry. The magnitude of the dissipative energy terms is maximum immediately behind the roughness trailing edge separation region and undergoes a very rapid decrease in a very short distance downstream. This decrease defines the streamwise location of maximum growth rate for each streamwise wavelength. Further downstream, both terms progressively decrease in magnitude. Their evolution is also found to be sensitive to interaction with boundary-layer disturbances.

Finally, it is important to note that no significant differences are observed between the disturbance energy extraction mechanisms of the two types of wake instabilities, only between the streamwise ranges over which the energy exchanges occur. The same finding is reported by Theiss et al. (Reference Theiss, Hein, Ali and Radespiel2016) for the streamwise evolution of the Reynolds stress production terms $\hat {P}_{RS,12}$ and $\hat {P}_{RS,13}$ near the trailing edge of different roughness geometries mounted on the forebody of a reentry capsule. For all the cases analysed by Theiss et al. (Reference Theiss, Hein, Ali and Radespiel2016), the wall-normal contribution ($\hat {P}_{RS,12}$) was also found to be larger than the spanwise one for both wake instabilities in the immediate vicinity of the roughness elements. In addition, the current results show that the streamwise evolution of the growth rate of the varicose and sinuous instabilities is not uniquely governed by the energy production terms associated to the base-flow gradients along a single spatial direction. Rather, the production terms linked to the gradients along $y$ are the leading contributions for certain streamwise regions, whereas the terms related to the gradients along $z$ are the most important contributions for other regions. It is observed that base-flow gradients along both the wall-normal and spanwise directions play an important role in the evolution of both varicose and sinuous disturbances.

6.4. Synchronism between wake instabilities and boundary-layer modes

This last section is devoted to illustrating a synchronization process between the wake-induced disturbances and the boundary-layer instabilities. Figure 26 depicts the magnitude of the streamwise velocity amplitude function of the SIN1 and VCOS1 instabilities excited by the cuboid roughness element at four different streamwise locations for $\lambda _x = 0.5\ \textrm {cm}$. As can be observed, after a given streamwise distance, the amplitude of the wake modes begins to increase in the regions located at the undisturbed boundary layer outside of the roughness wake, in a very similar way to the oblique Mack-mode instabilities shown in figure 12 (the same behaviour is observed for the ramp geometry, not shown here). This provides evidence of an interaction between the wake instability modes and Mack's mode family. For the same configuration, figure 27 shows the streamwise evolution of the phase speed ($c_{ph} = \omega _r/\alpha _r$) of the leading sinuous and varicose modes as well as of different instabilities located along the Mack-mode branch, namely, modes M2D, MO3, MO4 and MO7, already shown in figure 12. The phase speed of mode SIN1 begins to match the phase speed of the Mack-mode instabilities starting at $x \approx 0.2\ \textrm {m}$. This coincides with the location at which the eigenfunction of the sinuous disturbance starts to develop amplitude in the boundary layer at the sides of the roughness wake. The same is encountered for the varicose mode, for which the match in phase speeds begins significantly further downstream, approximately at $x = 0.35\ \textrm {m}$, once again in good agreement with the location of the change in its amplitude function according to figure 26. Since the Mack-mode instability manifests itself as a continuous branch in the 2D-LST spectrum, this means that there is a continuous synchronization between the wake instability modes and the Mack modes. Such a synchronization mechanism was suggested by De Tullio & Sandham (Reference De Tullio and Sandham2015) for the excitation of varicose modes, based on DNS results for a very similar cuboidal roughness configuration. In this work, this phenomenon is observed for both sinuous and varicose disturbances and for the two different roughness geometries investigated, bringing further weight to the importance of this mechanism in roughness-induced transition.

Figure 26.

Normalized magnitude of the streamwise velocity amplitude function of modes SIN1 $(a{,}c{,}e{,}g)$ and VCOS1 $(b{,}d{,}f{,}h)$ for the cuboid geometry at $\lambda _x = 0.5\ \textrm {cm}$ and at four different streamwise locations: $(a{,}b)$ $x = 0.15\ \textrm {m}$; $(c{,}d)$ $x = 0.2\ \textrm {m}$; $(e{,}f)$ $x = 0.27\ \textrm {m}$; $(g{,}h)$ $x = 0.38\ \textrm {m}$. Contour levels range between 0 and $0.5 |\hat {u}| / \max (|\hat {u}|)$ to favour the visualization of structures developing at the sides of the roughness wake.

Figure 27.

Streamwise evolution of the non-dimensional phase speed of different instability modes developing behind the cuboid roughness element for $\lambda _x = 0.5\ \textrm {cm}$. Solid and dashed lines denote the leading wake instability modes (VCOS1 and SIN1), and dash–dot lines with symbols correspond to different Mack instability modes (M2D, MO3, MO4 and MO7).

Under the synchronization process, the wake instabilities progressively transform into hybrid modes that develop in the roughness wake as well as at the surrounding boundary layer, in the same way as the instabilities belonging to the Mack-mode family that are modulated by the roughness wake. This transformation is reflected in the energy decomposition of the sinuous and varicose modes. Recalling the streamwise evolution of the temporal growth-rate decomposition terms (figures 22 and 23), it can be noticed that there is a close relationship between the location at which the synchronization starts and the position at which the leading production terms change their behaviour. For $\lambda _x = 0.5\ \textrm {cm}$, the location at which the term $\hat {P}_{s,2}$ starts to grow after the first energy exchange is approximately $x = 0.2\ \textrm {m}$ for mode SIN1 and $x = 0.35\ \textrm {m}$ for mode VCOS1, which precisely coincide with the locations at which their respective phase speeds match those of the Mack-mode branch. Downstream of these streamwise stations, the energy decomposition of modes SIN1 and VCOS1 is expected to become closer to that of the Mack-mode disturbances.

The temporal growth-rate budget for mode M2D is presented in figure 28. For comparison purposes, the decomposition based on LST for the two-dimensional Mack mode ($\beta = 0$) developing in a smooth flat-plate boundary layer is also included. In this case, the decomposition is evaluated by setting all base-flow spanwise derivatives to zero, so that only wall-normal production terms remain. As can be observed, the energy of the 2D-LST Mack mode is governed by the wall-normal production terms. Although its amplitude function is strongly modulated by the presence of the roughness, the energy signature of the 2D-LST Mack mode is still very close to the LST one, for which no spanwise production terms exist. This implies that the same is also true for its growth-rate evolution. These results confirm that, when the synchronization starts, the disturbance energy decomposition of modes SIN1 and VCOS1 begins to approach the energy signature of the Mack-mode instabilities, namely, spanwise terms ($\hat {P}_{s,3}$ and $\hat {P}_{RS,13}$) decreasing towards zero and wall-normal terms ($\hat {P}_{s,2}$ and $\hat {P}_{RS,12}$) rising. As the wake perturbations begin to develop amplitude in the boundary layer outside of the roughness wake, there is a redistribution of the disturbance energy from the spanwise components to the wall-normal ones.

Figure 28.

Streamwise evolution of the temporal growth-rate decomposition for the two-dimensional Mack-mode instability (M2D) found in the cuboidal roughness element wake for $\lambda _x = 0.5\ \textrm {cm}$: $(a)$ production terms; $(b)$ dissipation terms. Results for the two-dimensional Mack mode developing in the corresponding smooth flat-plate boundary layer, obtained by means of LST, are also included for comparison.

This is one mechanism by which the wake instabilities change their properties along the wake, and which has a direct impact on their amplification rates. It is therefore important to account for this phenomenon in the modelling of the linear stages of roughness-induced transition in conditions where wake and boundary-layer instabilities coexist.