1. Introduction
When high-speed transitional boundary layers encounter surface roughness, the resulting interaction is difficult to anticipate. Depending on the location and shape of the roughness, and the state of the boundary layer, transition to turbulence may be promoted or delayed. Even in the latter case, a longer sustained laminar flow may be a result of taming a particular route to turbulence, but other instability waves may become more amplified. Nonetheless, the potential thermal and mechanical benefits of a sustained laminar state in flight are such that roughness has previously been explored as an effective strategy to delay transition in high-speed boundary layers (Fujii Reference Fujii2006; Marxen, Iaccarino & Shaqfeh Reference Marxen, Iaccarino and Shaqfeh2010; Riley, McNamara & Johnson Reference Riley, McNamara and Johnson2014; Fong et al. Reference Fong, Wang, Huang, Zhong, McKiernan, Fisher and Schneider2015; Zhao, Dong & Yang Reference Zhao, Dong and Yang2019). In the present effort, we examine the influence of roughness location and shape on laminar-to-turbulence transition in a Mach 4.5 zero-pressure-gradient boundary layer using direct numerical simulation (DNS), and explain the observed shifts in transition location. We subsequently perform an optimization of the roughness parameters to achieve the longest possible delay of breakdown to turbulence within our simulation domain.
1.1. Stability of high-speed boundary layers
The work by Lees & Lin (Reference Lees and Lin1946) extended the Rayleigh inflection-point criterion to compressible flows. At the generalized inflection point, 
$\partial _y ( \bar {\rho } \partial _y \bar {u} )$ vanishes, where 
$\bar {\rho }$ and 
$\bar {u}$ are the base-state density and streamwise velocity and 
$\partial _y \equiv \partial / \partial y$ is the derivative in the wall-normal direction. The key implication is that a compressible, zero-pressure-gradient boundary layer can be inviscidly unstable, unlike the incompressible counterpart which has a viscous Tollmien–Schlichting instability. Mack (Reference Mack1969, Reference Mack1984) performed extensive linear stability computations of the compressible boundary layer and identified an infinite sequence of inflectional instabilities at high Mach number, the first two of which are now known as the first and second Mack modes. This spectrum of discrete, unstable modes was also identified by Smith & Brown (Reference Smith and Brown1990) and Cowley & Hall (Reference Cowley and Hall1990) using asymptotic analysis in the limit of infinite Mach number. The higher-order modes are reported to be unstable over relatively small ranges of high frequencies.
Mack's first mode reaches its maximum energy near the generalized inflection point of the boundary-layer profile, or the local maximum of 
$\bar {\rho } \partial _y \bar {u}$. Mack's higher modes, which only exist at high Mach number, are rooted in acoustic waves that are trapped inside the boundary layer near the wall in a region where the phase speed of the wave 
$c$ is locally supersonic relative to the mean flow 
$\bar {u}$, or 
$c - \bar {u} \geq a$ where 
$a$ is the local speed of sound. While most of these modes are inviscid, the analysis by Smith (Reference Smith1989) showed that the first modes are of a viscous-inviscid kind; directed outside of the local wave-Mach-cone direction, i.e. mode angles higher than 
$\tan ^{-1}(\sqrt {M^2_{\infty } - 1})$ where 
$M_\infty$ is the free stream Mach number, these modes exhibit a triple-deck structure.
Mack's modes are potential precursors of transition in compressible boundary layers, and both the first- and second-mode waves have been observed in various experiments (see e.g. Kendall Reference Kendall1975; Lysenko & Maslov Reference Lysenko and Maslov1984; Stetson & Kimmel Reference Stetson and Kimmel1992; Casper et al. Reference Casper, Beresh, Henfling, Spillers, Pruett and Schneider2016; Kegerise & Rufer Reference Kegerise and Rufer2016; Laurence, Wagner & Hannemann Reference Laurence, Wagner and Hannemann2016; Zhu et al. Reference Zhu, Chen, Wu, Chen, Lee and Gad-el Hak2018; Liu et al. Reference Liu, Yi, Xu, Shi, Ouyang and Xiong2019). At high-subsonic and moderate-supersonic speeds, boundary-layer transition in low-disturbance environments occurs as a result of excitation and amplification of instabilities that resemble Mack's first-mode waves. As the Mach number increases, the generalized inflection point moves to the outer region of the boundary layer and the growth rate of first-mode instabilities becomes smaller than the second-mode instabilities. The latter dominate transition in high-Mach-number supersonic and hypersonic boundary layers. According to inviscid linear stability (Mack Reference Mack1984), the growth rate of the second mode exceeds that of the first mode, for an adiabatic flat-plate boundary layer, at 
$M_\infty \approx 4$, where 
$M_\infty$ is the free stream Mach number. For cooled boundary layers, the second mode could become the dominant instability at even lower Mach numbers.
More recent work has pointed out that extending the terminology of first and second model to viscous flows may not be pertinent (Tumin Reference Tumin2007; Fedorov Reference Fedorov2011; Fedorov & Tumin Reference Fedorov and Tumin2011). Instead, the discrete modes were distinguished as follows. The slow mode, or mode S, has a phase speed that approaches the slow acoustic wave 
$c/u_\infty = 1 - 1 /M_{\infty }$ in the limit 
$\alpha \delta \ll 1$, where 
$\alpha$ is the streamwise wavenumber of the wave and 
$\delta$ is the local boundary-layer thickness. This limit is achieved, for example, near the leading edge when 
$\alpha$ is finite and 
$\delta \to 0$. In the same limit, the fast mode, or mode F, has a phase speed that approaches the fast acoustic wave 
$c/u_\infty = 1 + 1 /M_{\infty }$. Fedorov & Tumin (Reference Fedorov and Tumin2011) argued that, in terms of the spatial stability of an adiabatic wall at a finite Reynolds number, there only exists one unstable discrete mode, mode S, which exhibits features of the inviscid first-mode or second-mode instabilities depending on frequency and Reynolds number. Mode F can also become unstable, for example, in a cold-wall boundary layer. As the Reynolds number increases downstream of the leading edge, the phase speed of mode S increases and that of mode F decreases until they synchronize. The synchronization point is an important modulator for many stability features of high-speed boundary layers (Fedorov Reference Fedorov2011; Fedorov & Tumin Reference Fedorov and Tumin2011; Fong et al. Reference Fong, Wang, Huang, Zhong, McKiernan, Fisher and Schneider2015; Zhao et al. Reference Zhao, Wen, Tian, Long and Yuan2018; Park & Zaki Reference Park and Zaki2019; Dong & Zhao Reference Dong and Zhao2021; Jahanbakhshi & Zaki Reference Jahanbakhshi and Zaki2021). For example, Park & Zaki (Reference Park and Zaki2019) examined the sensitivity of the linear stability of high-speed boundary layers to the distortions of the base velocity and temperature profiles; they showed that the sensitivities of modes S and F to the distortions increase with Reynolds number, but near the synchronization point, there is a sudden drop and jump, respectively, in the sensitivities. Despite the arguments by Fedorov and Tumin, the terminology ‘first mode’ and ‘second mode’ remain widely adopted by the high-speed boundary-layer research community, and is therefore retained herein.
1.2. Effects of roughness on transition to turbulence
Introducing isolated or distributed roughness elements can promote breakdown to turbulence, relative to a smooth surface, by increasing the amplification rate of existing instabilities or spurring a new one. Examples include the generation of wakes and unstable shear layers downstream of tall roughness elements (Ergin & White Reference Ergin and White2006), and formation of streamwise vorticity behind shorter roughness elements that can initiate stationary cross-flow instabilities in three-dimensional (3-D) boundary layers (Radeztsky, Reibert & Saric Reference Radeztsky, Reibert and Saric1999). Our interest is, however, in carefully designed roughness elements that can delay transition.
Depending on the shape and location, roughness contributes to (i) receptivity (Dong, Liu & Wu Reference Dong, Liu and Wu2020; Liu, Dong & Wu Reference Liu, Dong and Wu2020) and/or (ii) local scattering of perturbations, by roughness-induced mean-flow distortion, non-homogeneous forcing and non-parallel effects (Wu & Dong Reference Wu and Dong2016; Xu et al. Reference Xu, Sherwin, Hall and Wu2016; Dong & Zhao Reference Dong and Zhao2021). When the roughness is small compared to the local boundary-layer height, these effects can be studied within the framework of triple-deck theory (Stewartson Reference Stewartson1969; Smith Reference Smith1973; Wu & Dong Reference Wu and Dong2016; Dong & Zhao Reference Dong and Zhao2021). Using a large-Reynolds-number asymptotic analysis, Dong et al. (Reference Dong, Liu and Wu2020) showed that the distortion of a small free stream acoustic wave by the curved wall of an isolated surface elements of height 
$h \ll \delta$ contributes to receptivity, and the amplitude of the resulting eigenmode scales with 
${O} ( h / \delta )$. In addition, the interactions between the roughness-induced mean-flow distortion and the acoustic wave leads to receptivity that scales as 
${O} ( [h/\delta ] Re_\delta ^{-1/3} )$, where 
$Re_\delta$ is the Reynolds number based on 
$\delta$. In another study, Liu et al. (Reference Liu, Dong and Wu2020) performed DNS and asymptotic analysis at moderate and large Reynolds numbers, respectively, to show that the amplitude of the excited viscous Mack first mode for the strong receptivity regime scales as 
${O} ( [h/\delta ] Re_x^{1/4} )$, where 
$Re_x$ is the Reynolds number based on the distance from the leading edge.
Two notable early experimental studies that reported transition delay on a roughed wall are the works by James (Reference James1959) and Holloway & Sterrett (Reference Holloway and Sterrett1964). The former examined two-dimensional (2-D) roughness elements in free flight tests with 
$2.8 < {Ma}_{\infty } < 7$, and the latter tested spherical roughness elements mounted on a flat-plate in a Mach 6 wind tunnel. More recent experiments in which roughness-induced delay of transition was observed were conducted by Fujii (Reference Fujii2006), Bountin et al. (Reference Bountin, Chimitov, Maslov, Novikov, Egorov, Fedorov and Utyuzhnikov2013) and Fong et al. (Reference Fong, Wang, Huang, Zhong, McKiernan, Fisher and Schneider2015). The measurements by Fujii (Reference Fujii2006) were on a Mach 7, 5 degree half-angle sharp cone on which either a wavy 2-D or spherical roughness elements were mounted. He found that, at high stagnation temperature and pressure conditions, transition was delayed when the wavelength of wavy-wall roughness is similar to the unstable second mode (approximately 
$2 \delta$); spherical roughness elements had little effect on the location of transition to turbulence. At low stagnation temperature and pressure conditions, however, both the wavy wall and spherical roughness elements promoted transition to turbulence relative to a smooth wall. Bountin et al. (Reference Bountin, Chimitov, Maslov, Novikov, Egorov, Fedorov and Utyuzhnikov2013) examined the effects of a wavy wall on stability of a Mach 6 boundary layer. They observed flow over the shallow-grooved plate was stabilized in a high-frequency band and destabilized at low frequencies, emphasizing that roughness must be carefully selected depending on the flow regime taking into account potential environmental disturbance spectra that may force the boundary layer. Fong et al. (Reference Fong, Wang, Huang, Zhong, McKiernan, Fisher and Schneider2015) studied a Mach 6 flow over a flared cone with initial half-angle of 2 degrees with six equi-spaced, 2-D elliptical roughness elements. They reported that second-mode instabilities can be damped if the roughness is placed downstream of the synchronization point of the fast and slow second modes. However, it is noteworthy that in this experiment, the roughness had an unanticipated effect of promoting the amplification of the first-mode instability.
Computational and theoretical studies have also examined the influence of roughness on high-speed boundary-layer stability (Duan, Wang & Zhong Reference Duan, Wang and Zhong2010; Marxen et al. Reference Marxen, Iaccarino and Shaqfeh2010; Riley et al. Reference Riley, McNamara and Johnson2014; Groskopf & Kloker Reference Groskopf and Kloker2016; Zhao et al. Reference Zhao, Dong and Yang2019; Dong & Zhao Reference Dong and Zhao2021; Haley & Zhong Reference Haley and Zhong2023). We first consider the impact of 2-D modifications of the surface. Using DNS, Duan et al. (Reference Duan, Wang and Zhong2010) studied a Mach 5.92 boundary layer on a flat plate with an isolated 2-D elliptic bump and drew similar conclusions as Fong et al. (Reference Fong, Wang, Huang, Zhong, McKiernan, Fisher and Schneider2015). The DNS by Marxen et al. (Reference Marxen, Iaccarino and Shaqfeh2010) examined the impact of 2-D hyperbolic-shaped isolated roughness on a flat-plate boundary layer at free stream Mach 4.8. Depending on modal frequency, the roughness element can either amplify or damp the disturbance waves. Riley et al. (Reference Riley, McNamara and Johnson2014) investigated the effect of 2-D compliant panels (convex or concave panel buckling) on boundary-layer stability for Mach 4 flow over a wedge. They used linear stability theory and the parabolized stability equations, and showed that placing panels near the leading edge of the wedge promotes naturally occurring high-frequency disturbances. However, placement near the trailing edge enhanced stability. Most recently, Zhao et al. (Reference Zhao, Dong and Yang2019) used the harmonic linearized Navier–Stokes equations to evaluate the effect of a 2-D hump or indentation on a flat wall on boundary-layer stability. They confirmed that the synchronization points of the instability waves are critical modulators of the impact of roughness. The earlier referenced study by Dong & Zhao (Reference Dong and Zhao2021), which developed large-Reynolds-number asymptotic theory for the impact of localized roughness on first and second Mack modes, attributed the dominant roughness effect to the interaction of the oncoming perturbation with the mean flow distortion in the main layer and the inhomogeneous forcing from the curved wall. Most recently, Haley & Zhong (Reference Haley and Zhong2023) investigated how a single roughness strip and an array of six sequential strips influence stability of second modes on a Mach 8 straight blunt cone. For both cases, they observed stabilization of high-frequency and destabilization of low-frequency modes. As for 3-D roughness elements, Groskopf & Kloker (Reference Groskopf and Kloker2016) studied their impact on a Mach 4.8 boundary layer over a flat plate. The roughness parameters included the spanwise width to streamwise length ratio, height and skewing angle with respect to the oncoming flow. The authors observed that, in the wake of obliquely placed elements, strong low-speed streaks are generated due to the induced cross-flow. Oblique placement was also associated with larger amplification of instabilities.
1.3. Objectives
The above discussion highlights both the wealth of discoveries from experimental, theoretical and numerical studies of roughness in transitional high-speed flows, and also the challenge in anticipating the impact of roughness on transition. The outcome depends on the roughness parameters, including shape and location, and on the details of the flow configuration. A choice that does not guarantee robust transition delay can, therefore, lead to an undesirable, potentially catastrophic outcome especially given the uncertainty in the environmental conditions relevant to high-speed flights. In the present work, we examine the key roughness parameters that impact transition to turbulence in a flat-plate boundary layer at Mach 4.5. We then optimize these parameters to achieve maximum transition delay within our computational domain. The inflow condition in our simulations is the nonlinearly most dangerous disturbance at the prescribed level of inlet energy, which was previously computed for the same configuration (Jahanbakhshi & Zaki Reference Jahanbakhshi and Zaki2019). The delay of the associated transition mechanism using surface roughness is analysed in detail.
This paper is organized as follows. Computational aspects, inflow condition and the roughness geometry are introduced in § 2, while validation of the computational model is provided in Appendix A. The results and discussions are reported in § 3, and are followed by a summary in § 4.
2. Computational framework
Direct numerical simulations are performed to study the influence of surface roughness on transitional high-speed boundary layers, and specifically the capacity to delay breakdown to turbulence. The flow satisfies the compressible Navier–Stokes equations, and a sample flow configuration is provided in figure 1. The contours show the streamwise velocity of a boundary layer over a protruding roughness. Throughout this work, flow variables are non-dimensionalized using the free stream velocity 
$\tilde {u}_\infty$, density 
$\tilde {\rho }_\infty$, temperature 
$\tilde {T}_\infty$ and viscosity 
$\tilde {\mu }_\infty$, where 
$\tilde {\bullet }$ represents dimensional quantities. Lengths are normalized using the Blasius scale at the inflow, 
$\tilde {l} = \sqrt {\tilde {\mu }_\infty \tilde {x}_0 / \tilde {\rho }_\infty \tilde {u}_\infty }$, where 
$\tilde {x}_0$ is the distance of the inflow plane from the virtual boundary-layer origin. In terms of these reference scales, we define the Reynolds number 
$Re_l \equiv \tilde {\rho }_\infty \tilde {u}_\infty \tilde {l} / \tilde {\mu }_\infty$, which is equivalent to 
$\sqrt {Re_{x_0}} = \sqrt { \tilde {\rho }_\infty \tilde {u}_\infty \tilde {x}_0 / \tilde {\mu }_\infty }$. Both values are equal to the non-dimensional location of the inflow plane, 
$x_0 = \tilde {x}_0 / \tilde {l} = \sqrt { Re_{x_0} } = Re_l$, while the Reynolds number based on the streamwise distance is given by 
$Re_x = Re_l x$.
Schematic of transitional boundary-layer flow over a flat plate with an isolated, protruding roughness element. Contours show the streamwise velocity 
$u$ normalized by its free stream value, 
$u_\infty$.
In non-dimensional form, the compressible Navier–Stokes equations are
\begin{equation} \frac{\partial \boldsymbol{q} }{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} ( \,\boldsymbol{f}_{i} - \boldsymbol{f}_{v} ) = 0,\end{equation}
where 
$\boldsymbol {q} = [ \rho\, \rho \boldsymbol {u} \ E ]^{\top }$ is the state vector, 
$\boldsymbol {f}_{i} = [ \rho \boldsymbol {u} \ \rho \boldsymbol {u} \boldsymbol {u} + p \boldsymbol{\mathsf{I}} \ \boldsymbol {u} ( E + p ) ]^{\top }$ represents inviscid fluxes, 
$\boldsymbol {f}_{v} = [ 0 \ \boldsymbol {\tau } \ (\boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\tau } - \boldsymbol {\theta })]^{\top }$ are the viscous fluxes and 
$(\bullet )^{\top }$ denotes the transpose. In the equations, 
$\rho$ is the density, 
$\boldsymbol {u} = [ u \ v \ w ]^{\top }$ is the 3-D velocity vector, 
$p$ is the pressure, 
$E = \rho e + 0.5 \rho \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {u}$ is the total energy, 
$e$ is the specific internal energy, 
$\boldsymbol {\tau }$ is the viscous stress tensor, 
$\boldsymbol {\theta }$ is the heat-flux vector and 
$\boldsymbol{\mathsf{I}}$ is the unit tensor. The system of equations is closed by assuming a calorically perfect gas whose thermodynamic properties are related by
\begin{equation} p = ( \gamma - 1 ) \rho e \quad \textrm{and} \quad T = \gamma ( \gamma - 1 ) M_\infty^2 e, \end{equation}
where 
$T$ is the temperature, 
$\gamma$ is the ratio of specific heats, 
$M_\infty = \tilde {u}_\infty [ (\gamma -1) \tilde {c}_p \tilde {T}_\infty ]^{-1/2}$ is the free stream Mach number and 
$\tilde {c}_p$ is the specific heat at constant pressure. The viscous stress and heat flux are modelled as
\begin{align} \boldsymbol{\tau} = \frac{\mu}{Re_l} \left[ \boldsymbol{\nabla} \boldsymbol{u} + ( \boldsymbol{\nabla} \boldsymbol{u})^{\top} + \left( \frac{\mu_b}{\mu} - \frac{2}{3} \right) ( \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} ) \boldsymbol{\mathsf{I}} \right] \quad \textrm{and} \quad \boldsymbol{\theta} ={-} \frac{ \mu \boldsymbol{\nabla} T }{(\gamma-1) M_\infty^2 Re_l Pr}, \end{align}
respectively, where 
$\mu$ is the dynamic shear viscosity, 
$\mu _b$ is the bulk viscosity and 
$Pr$ is the Prandtl number.
The DNS adopt a finite-difference discretization of the governing Navier–Stokes equations on a structured Cartesian grid. A sixth-order-accurate central difference scheme in the split form by Ducros et al. (Reference Ducros, Laporte, Soulères, Guinot, Moinat and Caruelle2000) is adopted for the inviscid fluxes, and is replaced by a fifth-order-accurate weighted essentially non-oscillatory scheme with Roe flux splitting near shocks. The viscous fluxes are computed using a conservative discretization that has the resolution characteristics of a sixth-order scheme. Time is advanced by a fourth-order accurate Runge–Kutta method. To simulate complex geometries, a cut-stencil method (Greene et al. Reference Greene, Eldredge, Zhong and Kim2016) is implemented that changes the discretization of the governing equations near the body and applies the boundary conditions just at the interface between the solid and the fluid. The method generates precise locations for the body. Validation of the original algorithm for Cartesian geometries is available in the literature (see e.g. Larsson & Lele Reference Larsson and Lele2009; Johnsen et al. Reference Johnsen2010; Kawai & Larsson Reference Kawai and Larsson2012; Volpiani, Bernardini & Larsson Reference Volpiani, Bernardini and Larsson2018), and the newly implemented cut-stencil method is discussed in more detail in Appendix A.
The operating gas is air for which the Prandtl number is 
$Pr = 0.72$ and the ratio of specific heats is 
$\gamma = 1.4$. The Mach number in the free stream is 
$M_{\infty } = 4.5$. Sutherland's law (Sutherland Reference Sutherland1893) models the temperature dependence of the dynamic viscosity, and Stokes’ hypothesis relates the dynamic and bulk viscosity coefficients. The streamwise position of the inflow plane is 
$x_0 = \sqrt {Re_{x_0}} = 1800$, which was selected based on the transition Reynolds numbers in high-altitude flight tests being 
$\sqrt {Re_{x_{tr}}} > 2000$ for 
$M_\infty > 4$ (Harvey Reference Harvey1978; Schneider Reference Schneider1999). At the inlet plane, the Blasius base-state is prescribed along with a superposition of linear instability waves. The amplitudes and relative phases of all the modes were optimized in an earlier study (Jahanbakhshi & Zaki Reference Jahanbakhshi and Zaki2019) such that they lead to the earliest possible breakdown to turbulence on an adiabatic flat plate (see also § 2.1). Periodic boundary conditions are enforced in the homogeneous spanwise direction, convective outflow is prescribed at the right and top boundaries, and the bottom boundary is a no-slip adiabatic wall.
The extents of the computational domain in the streamwise, wall-normal and spanwise directions are respectively 
$L_x = 2984$, 
$L_y = 204$ and 
$L_z = 150$. The domain is discretized using a uniform grid in the horizontal directions with 
$N_x = 2985$ and 
$N_z = 151$ points. A hyperbolic tangent stretching of the grid with 
$N_y = 189$ points is used to discretize the wall-normal direction. The wall-normal grid spacing that was adopted throughout this work is reported in figure 2(a), as well as a finer grid that was used to verify grid independence. At the inlet plane, the main grid has 54 points within the boundary layer. Figure 2(b,c) provides evidence of grid independence by comparing the main and finer grids. The figures report the skin-friction coefficient,
\begin{equation} C_f \equiv \frac{ \tau_{wall}}{ 0.5 \rho_{\infty} u_{\infty}^2 }, \end{equation}and the disturbance energy,
\begin{equation} \mathcal{E} \equiv \frac{1}{2} \int_{0}^{L_y} \left( \bar{\rho} \{ \overline{\boldsymbol{u}^{\prime} \boldsymbol{\cdot} \boldsymbol{u}^{\prime}} \} + \frac{\bar{\rho} \bar{T}}{(\gamma-1) \gamma M_\infty^2} \left\{ \frac{\overline{\rho^{\prime 2}}}{\bar{\rho}^2} + \frac{\overline{T^{\prime 2}}}{\bar{T}^2} \right\} \right) {{\rm d} y}, \end{equation}
as a function of 
$\sqrt {Re_x}$. Over-line denotes averaging in time and in the homogeneous spanwise direction, and the primed variables are the perturbations with respect to this average. As demonstrated in figure 2, a finer resolution does not yield any perceptible changes in the skin friction, the disturbance energy or the mean streamwise-velocity contours (we also confirmed grid-independence of the spanwise and time-averaged variance). Additionally, we verified that the instability modes prescribed at the inlet and their nonlinear interactions, leading to transition, are all fully resolved, and that our predictions of transition are grid independent.
(a) Wall-normal grid size and (b–d) sample flow over the optimal roughness for delaying transition. (b) Downstream behaviour of the skin-friction coefficient. (c) Downstream evolution of the disturbance energy. (d) Contours of averaged streamwise velocity in the vicinity of the roughness; colours are from the herein adopted resolution 
$(N_x, N_y, N_z) = (2985, 189, 151)$, and dashed lines are from the finer grid 
$(N_x, N_y, N_z) = (3978, 279, 201)$.
2.1. Inflow instability modes
The inflow disturbance is synthesized as a superposition of linear stability eigenmodes of the local boundary-layer profile, which span the relevant range of frequencies 
$\omega$ and spanwise wavenumbers 
$\beta _z$:
\begin{equation} \boldsymbol{q}^{\prime}(x_0, y, z, t) = \sum_{\omega, \beta_z} \Re \{ \breve{\boldsymbol{q}}_{\langle \omega , \beta_z \rangle }(y) \exp( \alpha x_0 - {\rm i}( \beta_z z + \omega t) )\}. \end{equation}
The eigenmodes at each 
$\langle \omega, \beta _z\rangle$ pair are obtained by substituting the ansatz 
$\breve {\boldsymbol {q}}(y) \exp ( \alpha x - {\rm i} ( \beta _z z + \omega t) )$ in the linear perturbation equations, and the resulting spatial eigenvalue problem is solved for the spectrum of eigenfunctions 
$\breve {\boldsymbol {q}}(y)$ and associated complex eigenvalues 
$\alpha = \alpha _r + i \alpha _i$. Only the linearly most unstable wave, which in the present study is the slow mode, was retained at each 
$\langle \omega, \beta _z\rangle$ pair. The linear-stability operators and solution algorithm are standard (see e.g. Park & Zaki (Reference Park and Zaki2019), for details).
For a prescribed level of the total disturbance energy, Jahanbakhshi & Zaki (Reference Jahanbakhshi and Zaki2019) optimized the amplitudes and phases of the inflow instability waves to achieve the earliest possible transition to turbulence on a flat plate. The spectral makeup of the resulting nonlinearly most dangerous inflow disturbance is provided in figure 3(a), as a function of the normalized frequency and integer spanwise wavenumber:
\begin{equation} F \equiv \frac{\omega}{\sqrt{Re_{x_0}}} \times 10^6 \quad \textrm{and} \quad k_{z} \equiv \frac{\beta_{z} L_z}{2 {\rm \pi}}. \end{equation}
The modal energy, 
$\mathcal {E}_{\langle F , k_z \rangle }$, is defined as
\begin{equation} \mathcal{E}_{\langle F , k_z \rangle } = \frac{1}{2} \int_{0}^{L_y} \left( \bar{\rho} \{ \boldsymbol{\hat{u}}^* \boldsymbol{\hat{u}} \} + \frac{\bar{\rho} \bar{T}}{(\gamma-1) \gamma M_\infty^2} \left\{ \frac{ \hat{\rho}^* \hat{\rho} }{\bar{\rho}^2} + \frac{ \hat{T}^* \hat{T} }{\bar{T}^2} \right\} \right)_{\langle F , k_z \rangle } \,{{\rm d} y} , \end{equation}
where hatted variables are the Fourier coefficients in frequency-spanwise wavenumber space and star denotes the complex-conjugate transpose. At the inlet plane, the total disturbance energy is 
$\sum _{F,k_z} \mathcal {E}_{\langle F , k_z \rangle } = 2 \times 10^{-5}$.
The inflow instability waves. (a) The spectral energy. (b–g) Profiles of the important modes: (b,c) mode 
$\langle 110 , 0 \rangle$, (d,e) mode 
$\langle 100 , 0 \rangle$ and ( f,g) mode 
$\langle 40 , 3 \rangle$. Magnitude (thick); real part (dash–dotted); imaginary (dotted) part of the mode shapes. The horizontal thin lines (top to bottom) mark the 99 % boundary-layer thickness and the relative sonic line 
$c - \bar {u} = a$.
Figure 3(a) shows that the majority of the total energy, approximately 95 %, is assigned to only three inlet waves: 
$\langle 110 , 0 \rangle$; 
$\langle 100 , 0 \rangle$ and 
$\langle 40 , 3 \rangle$. The parameters for these three waves are provided in table 1; our simulations resolve each wavelength with 29 to 75 grid points in the streamwise direction and 50 to 151 points in the span. The streamwise velocity and pressure of the associated eigenfunctions are plotted in figures 3(b)–3(g). The pressure profiles in figures 3(c) and 3(e), for modes 
$\langle 110 , 0 \rangle$ and 
$\langle 100 , 0 \rangle$, have a single zero crossing in the wall-normal direction, and are second-mode instabilities. Jahanbakhshi & Zaki (Reference Jahanbakhshi and Zaki2021) decomposed the momentum-density vector of these waves into acoustic, entropic and solenoidal components, and confirmed that both 
$\langle 110 , 0 \rangle$ and 
$\langle 100 , 0 \rangle$ have acoustic characteristics, in agreement with the interpretation by Mack (Reference Mack1984). The associated pressure perturbations reflect back and forth between the wall and the sonic line of the relative flow. Mode 
$\langle 40 , 3 \rangle$ with no zero crossings is a first-mode vortical instability. Each of these three modes plays an important role in causing the earliest transition location at the present total level of energy. For example, removing one of the modes and redistributing its energy to the other instability waves leads to a downstream shift in transition onset.
Parameters of the inlet instability waves. Definitions are provided in § 2.1.
Figures 3(b), 3(d) and 3( f) show that the second-mode instabilities reach their maximum values below the relative sonic line close to the wall, while the first-mode instability is appreciable beyond the relative sonic line, close to boundary-layer edge. The difference in the wall-normal dependence of the two classes of instabilities makes the second-mode waves more susceptible to control strategies that are applied at the wall, e.g. short roughness elements or wall heating/cooling.
The synchronization of the fast and slow modes has a significant impact on the instability waves. Figure 4 reports the spatial growth rates, streamwise wavenumbers and phase speeds of modes F and S at 
$\langle F , k_z \rangle = \langle 110 , 0 \rangle$ and 
$\langle 100 , 0 \rangle$. These results are obtain using parallel, spatial, linear-stability theory and confirm typical characteristics of slow and fast modes for a high-Mach-number flow over an adiabatic flat-plate (Fedorov & Tumin Reference Fedorov and Tumin2011). The growth rates of modes F, in figures 4(a) and 4(b), show a discontinuity at the location where the phase speed approaches unity, which corresponds to crossing the continuous branch of entropy and vorticity modes. Figures 4 (c) and 4(d) show that at low Reynolds number, modes S and F follow the slow and fast acoustic branches, 
$c = 1 \mp 1/M_{\infty }$, of the spectrum. Figures 4(e) and 4( f) show the evolution of the phase speeds, which approach one another with Reynolds number until they synchronize. The Reynolds numbers at synchronization are 
$\sqrt {Re_x} = 1811$ for mode 
$\langle F , k_z \rangle = \langle 110 , 0 \rangle$ and 
$\sqrt {Re_x} = 1988$ for mode 
$\langle 100 , 0 \rangle$.
(a,b) Spatial growth rate, 
$\alpha _r$, (c,d), streamwise wavenumber, 
$\alpha _i$, and (e, f) phase speed, 
$\omega / \alpha _i$, of modes S (solid line) and F (dash–dotted lines) corresponding to the instabilities (a,c,e) 
$\langle 110 , 0 \rangle$ and (b,d, f) 
$\langle 100 , 0 \rangle$. The red dots indicate the synchronization point, the dotted lines mark the start of the computational domain in our DNS and the dashed lines represent the phase speed of acoustic waves. Reproduction of figure 5 by Jahanbakhshi & Zaki (Reference Jahanbakhshi and Zaki2021).
2.2. Geometry of roughness elements
The wall geometry 
$\hat {h}_r(x)$ is defined over three streamwise segments: 
$x < X_0$, 
$X_0 \leq x \leq X_N$ and 
$x > X_N$. In the first and last segments, the wall is flat, and therefore 
$\hat {h}_r = 0$. In the range 
$X_0 \leq x \leq X_N$, a 2-D roughness element is defined,
\begin{equation} \hat{h}_r = \frac{H_r}{2} [ \tanh ( L_r [ x - X_r + W_r ] ) - \tanh ( L_r [ x - X_r - W_r ] ) ] \sin \left( n{\rm \pi} \frac{ x - X_0}{X_N - X_0} \right), \end{equation}
where 
$H_r \geq 0$, 
$W_r \geq 0$, 
$L_r \geq 0$, 
$X_r = ({X_0 + X_N})/{2}$ and 
$n \in \mathbb {Z}$ are the geometrical parameters that respectively determine the maximum height, streamwise extent, abruptness, centre (location of maximum height) and streamwise integer wavenumber of the roughness. Note that the surface function is not 
$C^2$ continuous at 
$x = X_0$ and 
$X_N$, which can trigger numerical instability. To obtain a surface function that is 
$C^2$ at all 
$x$, we reconstruct the above surface topography using a fourth-order spline in which the first and final knots are at 
$x = X_0$ and 
$x = X_N$. This spline reconstruction is given by 
$\boldsymbol {h}_{r} = \boldsymbol{\mathsf{B}}_{4} [ ( \boldsymbol{\mathsf{B}}^{\top }_4 \boldsymbol{\mathsf{B}}_4 )^{-1} \boldsymbol{\mathsf{B}}^{\top }_4 \hat {\boldsymbol {h}}_{r} ]$, where 
$\boldsymbol{\mathsf{B}}_4 \in \mathbb {R}^{N_x \times N_c}$ is a matrix whose columns are B-splines of order 4 and 
$N_c$ is the number of control points. In a discretized domain, 
$\boldsymbol {h}_{r} \equiv h_r (x)$ and 
$\hat {\boldsymbol {h}}_{r} \equiv \hat {h}_r (x)$ are vector quantities in 
$\mathbb {R}^{N_x \times 1}$ space. The size of the control points, 
$N_c$, is set such that the residuals function 
$(\boldsymbol {h}_{r} - \hat {\boldsymbol {h}}_{r})^{\top }(\boldsymbol {h}_{r} - \hat {\boldsymbol {h}}_{r}) < 10^{-5}$.
The geometrical parameters of eight roughness elements are summarized in table 2. These surface topographies will be adopted in numerical simulations to examine the effect on transition, and the results will be discussed in § 3.1. The eight configurations involve two roughness locations and four wavenumbers: the designation ‘X1’ references cases with roughness elements positioned between 
$X_0 = 2000$ and 
$X_N = 2250$, while ‘X2’ has 
$X_0 = 2250$ and 
$X_N = 2500$. These locations were informed by the synchronization locations of the two most energetic inflow second-mode instabilities. The position X1 is post-synchronization of mode 
$\langle 110 , 0 \rangle$ and pre-synchronization of mode 
$\langle 100 , 0 \rangle$, whereas X2 is located post-synchronization of both second modes. The streamwise extent of these roughness elements, 
$X_N - X_0 = 250$, is much longer than the streamwise wavelength of the three dominant inlet instability waves (see table 1). The designations ‘p1’ and ‘p3’ reference roughness elements with positive values of 
$n = 1$ and 
$n = 3$, respectively, while ‘m1’ and ‘m3’ refer to negative values of 
$n$. Figure 5 shows the dependence of the surface geometry on the four values of 
$n$ at the upstream location. As this figure shows, the main feature of the surfaces defined by 
$n = 1$ and 
$n = -3$ is their protrusion above the flat wall, whereas roughness elements with 
$n = -1$ and 
$n = 3$ are primarily indentations, or cratering of the surface.
Geometrical parameters for examining effects of location and shape of the roughness on transition to turbulence. The boundary-layer thickness in the reference, flat-plate case is 
$\delta _{99} = 14.7$ at X1 and 
$\delta _{99} = 15.6$ at X2.
Shape of the roughness elements positioned between 
$X_0 = 2250$ and 
$X_N = 2500$: (a) 
$n =\pm 1$ and (b) 
$n =\pm 3$. The remaining roughness parameters are reported in table 2.
3. Results and discussion
In this section, we examine of the effects of roughness parameters on the location where the flow transitions to a turbulent state. Specifically, we examine the influence of the roughness location, streamwise integer wavenumber, maximum height, width and abruptness, 
$\{X_r, n, H_r, W_r, L_r\}$ in (2.9). The first two geometrical parameters are examined in § 3.1, while the remaining three are optimized in § 3.2 to achieve maximum transition delay in our computational domain.
3.1. Effects of location and shape of the roughness on transition
While a general surface topography that delays transition can be sought by optimization, we consider localized roughness since the associated impact on the flow is less ambiguous to analyse and due to practical considerations. The optimization can be performed for all the roughness parameters. However, we will first demonstrate the impact of the roughness location and general shape, 
$\{X_r, n\}$, on transition. This initial set of simulations will serve two roles. First, the simulations will be analysed to determine the key impact of roughness on transition dynamics. Second, based on these simulations, we will select the initial location and shape of the roughness whose parameters 
$\{H_r, W_r, L_r\}$ we will further optimize.
For each of the roughness configurations listed in table 2, simulations of transition were performed when the inflow disturbance is the superposition of instability waves summarized in § 2.1. The skin-friction curves from these simulations are reported in figure 6. Figure 6(a) highlights the cases in which the roughness can mainly be considered as a protrusion, whereas figure 6(b) shows the indentation cases. It is evident that both the location and shape of the roughness appreciably affect the evolution of the incoming boundary-layer disturbances as manifested by the observed shift in the location of transition to turbulence. Figure 6(a) shows that the examined protrusions can more effectively delay transition when they are placed at X2, i.e. post-synchronization of the second modes. The present simulations are therefore consistent with previous studies that noted the importance of the position of roughness relative to the synchronization point of the oncoming disturbance waves (Fong et al. Reference Fong, Wang, Huang, Zhong, McKiernan, Fisher and Schneider2015; Zhao et al. Reference Zhao, Wen, Tian, Long and Yuan2018; Dong & Zhao Reference Dong and Zhao2021). The increase in the streamwise wavenumber from one to three further delays transition, at either locations of the roughness. Compared to 
$n = 1$, the roughness with 
$n= - 3$ is a more slender protrusion into the flow that has a relatively larger impact on the boundary-layer thickness due to the cratering that precedes the peak. The results in figure 6(a) are specific to roughness elements whose primary feature is a protrusion. When the primary feature is an indentation, or cratering, transition is similarly delayed relative to the flat plate. However, some of the trends are reversed, as shown in figure 6(b): the most appreciable delay in transition takes place when the roughness location is upstream, at X1, and the streamwise wavenumber 
$n = -1$ is more effective compared to 
$n = 3$.
Skin-friction coefficients for cases where the primary feature of the roughness is (a) protrusion and (b) indentation. The geometric parameters of each roughness are provided in table 2. The results for transition on a flat plate are included for reference.
The shift in transition location reported in figure 6 should be viewed in the context of the effect of each roughness element on the near-wall flow. In hypersonic boundary layers, an important modulator of the amplification of second-mode instabilities is the thickness of the near-wall region where the instability phase speed is supersonic relative to the flow. In this region, the pressure waves reflect at the wall and at the relative sonic line (
$c - \bar {u} = a$) where the waves change from compression to expansion and vice versa (Morkovin Reference Morkovin1987). These waves are typically phase-tuned with the harmonic vorticity and temperature waves that are travelling along the relative sonic line. If the thickness of the relative supersonic region changes due to surface modifications, e.g. roughness elements or localized cooling/heating, the growth rate of the second-mode instabilities is altered. In the case of roughness elements, parameters such as location relative to the synchronization point, streamwise wavenumber, height compared to the relative sonic line and abruptness/width relative to wavelength of instabilities all affect the outcome.
When the primary feature of the roughness is a protrusion, modification of the relative supersonic region pre-synchronization can have an appreciable net destabilizing effect. This trend is reversed when the modification to the relative supersonic region takes place post-synchronization; the net effect on the instability waves becomes stabilizing. This trend is captured in the growth rate of mode 
$\langle 100 , 0 \rangle$ for which the locations X1 and X2 are pre and post its synchronization location (
$\sqrt {Re_x} = 1988$). Figure 7(a,b) quantify this behaviour, where the growth rate of mode 
$\langle 100 , 0 \rangle$ is reported for protrusions at X1p1 and at X2p1 (the figures also show the growth rate of mode 
$\langle 110 , 0 \rangle$). The contours are iso-levels of the instantaneous streamwise velocity. Upstream of the roughness, there is a zone where the iso-lines diverge away from the wall, specifically where the relative supersonic region thickens. Atop the roughness, the iso-lines converge and thus the relative-supersonic region thins. Finally, within a third zone downstream of the roughness, the iso-lines recover their natural boundary-layer spreading. The wave is initially destabilized before it reaches the roughness, followed by a stabilization region on top of the roughness and finally a re-adjustment zone across which the instability essentially recovers its amplification rate for an undisturbed boundary layer. When the protrusion is positioned upstream of the synchronization point of mode 
$\langle 100 , 0 \rangle$, the net effect was a higher modal energy of this wave. Placement of the protruding roughness downstream of the synchronization point is shown in figure 7(b). The final energy attained by mode 
$\langle 100 , 0 \rangle$ in this case is lower than the upstream placement of this roughness element. The above dependence of the amplitude of mode 
$\langle 100 , 0 \rangle$ on the roughness location may seem at odds with transition being delayed relative to the flat plate for both roughness locations (figure 6). For the explanation, it is important to recall that both roughness configurations are downstream of synchronization of the other dominant second-mode wave, namely 
$\langle 110 , 0 \rangle$. As a result, while this mode also exhibits the three-stage stabilization/destabilization/readjustment behaviour across the roughness, it is ultimately stabilized in both cases (figures 7a and 7b). Comparison of the DNS results with computations of modal evolution using linearized Navier–Stokes equations (not shown here) confirmed that the destabilization zones of modes 
$\langle 100 , 0 \rangle$ and 
$\langle 110 , 0 \rangle$ are primarily caused by the roughness-induced base-flow distortion, while the stabilization and re-adjustment zones are strongly influenced by the interaction of the instability waves with the finite slope of the roughness geometry.
Contours of instantaneous streamwise velocity for (a) X1p1, (b) X2p1, (c) X1m1 and (d) X2m1. Right axis is the energy of modes 
$\langle 110 , 0 \rangle$ and 
$\langle 100 , 0 \rangle$ for each case (dotted blue and red lines) compared to the reference flat-plate case (dash–dotted black lines). The relative sonic lines of modes 
$\langle 100 , 0 \rangle$ and 
$\langle 110 , 0 \rangle$ are marked by dark and light green dashed lines, and are mostly overlapping.
When the main roughness feature is an indentation (figure 7c,d), the relative supersonic region initially thins, then thickens and finally thins again, with an associated stabilization/destabilization/stabilization and re-adjustment of the second-mode waves. Roughness location at either X1 or X2 leads to a net stabilization for both second-mode waves. For mode 
$\langle 100 , 0 \rangle$, this outcome is noteworthy because the synchronization location of this mode is downstream of X1 and upstream of X2. According to Dong & Zhao (Reference Dong and Zhao2021) (see their figures 16 and 18), indentations have a similar albeit weaker scattering effect as protrusions, and hence one expects that the instability wave is destabilized when an indentation is introduced pre-synchronization and stabilized when the indentation is downstream of synchronization. In contrast to that work, in our nonlinear simulations where the roughness is relatively short in the streamwise direction, mode 
$\langle 100 , 0 \rangle$ is stabilized in both configurations.
Similar to earlier works, the herein considered roughness parameters were guided by knowledge of (a) the importance of the synchronization point, (b) the sensitivity of the modal amplification to distortions in the base flow and (c) the non-parallel effects induced by the roughness. The roughness parameters were not, however, optimized to guarantee a particular outcome. We next consider such optimization; specifically, we seek an optimal roughness geometry that can lead to sustained laminar flow throughout the entire computational domain in our simulations.
3.2. Optimal protruding roughness for transition delay
In this section, we present the nonlinear optimization of the roughness height, width and abruptness, and examine the flow field associated with the optimal roughness.
3.2.1. Ensemble-variational optimization of the roughness
The base design that is adopted as a starting point of the optimization has streamwise wavenumber 
$n = 1$, which corresponds to a simple protrusion. This choice is primarily motivated by its simplicity to aid the interpretation of the impact on the flow and ease of manufacturing relative to cratering of the surface. Figure 6(a) showed that protrusions at X2 are more effective in delaying transition. Therefore, 
$X_0 = 2250$ and 
$X_N = 2500$ are selected as the start and end positions of the roughness element.
The optimization is performed using an ensemble-variational (EnVar) approach. Starting from the initial design, here a roughness with 
$\{H_r, W_r, L_r\}=\{0.22\delta _{x_0}, 5.37\delta _{x_0}, 0.81 \delta ^{-1}_{x_0} \}$ where 
$\delta _{x_0} = 13.5$ is the boundary-layer thickness at the inlet plane, EnVar updates the estimate of the control vector, 
$\boldsymbol {c} = [ H_r \ W_r \ L_r ]^{\top }$, at the end of each iteration using the gradient of the cost function; the gradient is evaluated from the outcomes of an ensemble of possible solutions. The cost function is the integrated skin-friction coefficient along the plate which, once minimized, ensures the farthest possible downstream location of transition to turbulence. The iterative optimization is halted once the identified roughness can maintain a laminar state throughout the entire computational domain. It should be noted that the solution to the above nonlinear optimization is not unique and, similar to other gradient-based methods, the reported results are only guaranteed to be a local optimum for a given choice of the initial guess. However, as long as the discovered roughness design accomplishes the objective of the optimization, it is deemed successful.
The choice of the EnVar technique is justified because of the low-dimensional nature of the control vector. In addition, unlike adjoint methods which may place limits on the time horizon of the flow solution (Zaki & Wang Reference Zaki and Wang2021), EnVar is perfectly suited for long-time integration and evaluation of statistical cost functions (Mons, Wang & Zaki Reference Mons, Wang and Zaki2019; Mons, Du & Zaki Reference Mons, Du and Zaki2021). EnVar has successfully been adopted in high-speed boundary layers for assimilation of measurements into flow simulations (Buchta & Zaki Reference Buchta and Zaki2021; Buchta, Laurence & Zaki Reference Buchta, Laurence and Zaki2022) and optimization (Jahanbakhshi & Zaki Reference Jahanbakhshi and Zaki2021). The reader is referred to those studies for the details of the algorithm.
The optimization procedure seeks the roughness parameters that minimize the cost function:
\begin{equation} \mathcal{J} = \underbrace{ \tfrac{1}{2} \Vert \boldsymbol{c} - \boldsymbol{c}^{(e)} \Vert_{\boldsymbol{\mathsf{B}}^{{-}1}}^2 }_{\mathcal{J}_p} + \underbrace{ \tfrac{1}{2} \Vert \mathcal{G} ( \boldsymbol{q} ) \Vert_{\boldsymbol{\mathsf{R}}^{{-}1}}^2 }_{ \mathcal{J}_o }. \end{equation}
The first term 
$\mathcal {J}_p$ is a regularization which ensures that the optimal 
$\boldsymbol {c}$ at the end of each iteration does not deviate appreciably from the previous estimate 
$\boldsymbol {c}^{(e)}$. The second term, 
$\mathcal {J}_o$, is the objective function that is formulated by defining 
$\mathcal {G} ( \boldsymbol {q} ) = C_f ({ {{\rm d}\kern0.06em x} / L_x })^{{1}/{2}}$, where the skin-friction coefficient is 
$C_f = \tau _{wall} / ( \frac {1}{2} \rho _{\infty } u_{\infty }^2 )$. In (3.1), 
$\boldsymbol{\mathsf{B}}$ and 
$\boldsymbol{\mathsf{R}}$ are the co-variance matrices of the prior term and the observed skin friction, respectively. To ensure that the maximum slope of the predicted roughness at the end of each iteration of the EnVar algorithm is appropriately resolved by the grid, a constraint is introduced in the optimization procedure. Specifically, we require that 
$| d \boldsymbol {h}_r / {{\rm d}\kern0.06em x} | \leq s$, where 
$s$ is a pre-determined limit on the slope of the roughness which is chosen based on the grid resolution, 
$s \sim {O}({\rm \Delta} y_{{near wall}} / {\rm \Delta} x)$.
The skin friction and normalized cost function, 
$C_f$ and 
$\mathcal {J}_o/\mathcal {J}_{o,0}$, associated with the mean control vector 
$\boldsymbol {c}$ at the end of each EnVar iteration are reported in figure 8(a,b). These plots demonstrate that the location, where the boundary layer transitions to a turbulent state, shifts downstream after consecutive iteration. In other words, the optimization procedure is effective at updating the roughness parameters in a manner to reduce the cost function (3.1) and delay breakdown to turbulence. At the end of iteration # 4, the boundary layer is laminar throughout the computational domain.
The outcomes of EnVar iterations. (a) Skin friction of the reference flow over a flat plate and the predictions at consecutive iterations versus 
$\sqrt {Re_x}$. (b) Cost function of each iteration normalized by the initial guess. (c) The geometry of the optimal roughness, after the fourth EnVar iteration, versus streamwise distance 
$x$.
The height and slope of the optimal roughness after the fourth EnVar iteration are depicted in figure 8(c). The associated geometric parameters are 
$\{H_r, W_r, L_r\} = \{0.31 \delta _{x_0}, 4.84 \delta _{x_0}, 0.89 \delta ^{-1}_{x_0} \}$. Therefore, compared to the initial guess, the optimal roughness is taller, more slender and more abrupt. Despite the larger height, the optimal roughness is still below the relative sonic line inside the boundary layer. As for the width 
$2 W_r$, it is approximately 4.1 times the streamwise wavelength of mode 
$\langle 100, 0 \rangle$ and 4.5 times the wavelength of mode 
$\langle 110, 0 \rangle$ (see table 1). To provide an appreciation for the physical size of the roughness elements, we can convert the current roughness parameters to dimensional quantities. For this purpose, we adopt the reference scales from the experiments by Kendall (Reference Kendall1975) which examined the stability of a flat-plate boundary layer at the same free stream Mach number and temperature. The optimal dimensional parameters of our roughness then become 
$\{H_r, W_r, L_r\} = \{1.1\ \textrm {mm}, 16.3 \ \textrm {mm}, 0.264\ \textrm {mm}^{-1} \}$, which are physically and practically relevant. For example, in the experiment by Fong et al. (Reference Fong, Wang, Huang, Zhong, McKiernan, Fisher and Schneider2015) for a Mach 6 boundary layer, the height and width of the roughness were set to 
$H_r = 0.665 \ \textrm {mm}$ and 
$W_r = 2.66 \ \textrm {mm}$, respectively; Bountin et al. (Reference Bountin, Chimitov, Maslov, Novikov, Egorov, Fedorov and Utyuzhnikov2013) set the roughness height and width in their experiment for a Mach 6 boundary layer to 
$H_r = 1.8 \ \textrm {mm}$ and 
$W_r = 12 \ \textrm {mm}$; Fujii (Reference Fujii2006) Mach 7 experiments included roughness height in the range 
$0.5 \ \textrm {mm}$ to 
$0.9 \ \textrm {mm}$ and roughness width in the range 
$0.5\ \textrm {mm}$ to 
$4 \ \textrm {mm}$.
3.2.2. Effect of optimal roughness on transition mechanism
Two instantaneous flow fields are contrasted in figure 9: the structures are iso-surfaces of the second invariant of the velocity-gradient tensor and are coloured by the spanwise velocity perturbations. Figure 9(a) is the reference simulation over a flat plate, while figure 9(b) shows the outcome of the optimization procedure. The visual difference highlights the efficacy of the optimized roughness in suppressing the transition precursors and delaying the onset of turbulence.
Isosurfaces of 
$Q$-criteria for (a) flow over a flat plate and (b) the optimal roughness.
The impact of the optimized roughness element on the downstream development of key instability waves and nonlinear energy exchanges is examined in figure 10. In figure 10(a,b), we reproduce the shape of the optimized roughness and mark the synchronization locations of the slow and fast modes for waves 
$\langle 110, 0 \rangle$ in figure 10(a) and 
$\langle 100, 0 \rangle$ in figure 10(b); the roughness location is downstream of both synchronization points. Figure 10(c–h) compare the downstream development of the spectral energy, 
$\mathcal {E}_{\langle F , k_z \rangle }$, over a flat plate (grey lines) and in the presence of the optimal roughness (black lines). Figure 10(i–l) report the nonlinear energy transfer among wave triads 
$\{ \langle F , k_{z} \rangle, \langle F_1 , k_{z,1} \rangle, \langle F_2 , k_{z,2} \rangle \}$, which is evaluated using
\begin{equation} \mathcal{I}_{\langle F , k_{z} \rangle } = \sum_{{\pm} F , \pm k_{z}} \int_{0}^{L_y} \vert \boldsymbol{\hat{\mathcal{A}}}^*_{\langle F_1 , k_{z,1} \rangle } \boldsymbol{\hat{\mathcal{B}}}_{\langle F_2 , k_{z,2} \rangle } + \boldsymbol{\hat{\mathcal{A}}}^*_{\langle -F_2 , -k_{z,2} \rangle } \boldsymbol{\hat{\mathcal{B}}}_{\langle -F_1 , -k_{z,1} \rangle } \vert \,{{\rm d} y},\end{equation}
where 
$F = F_2 - F_1$ and 
$k_z = k_{z,2} - k_{z,1}$ (Cheung & Zaki Reference Cheung and Zaki2010; Jahanbakhshi & Zaki Reference Jahanbakhshi and Zaki2019). The symbols 
$\boldsymbol {\hat {\mathcal {A}}}$ and 
$\boldsymbol {\hat {\mathcal {B}}}$ are the Fourier coefficients of 
$\boldsymbol {\mathcal {A}} = [ \rho u^{\prime \prime } u^{\prime \prime } \ \rho u^{\prime \prime } v^{\prime \prime } \ \rho u^{\prime \prime } w^{\prime \prime } ]^{\top }$ and 
$\boldsymbol {\mathcal {B}} = [ u^{\prime \prime } v^{\prime \prime } \ w^{\prime \prime } ]^{\top }$, where the double prime denotes fluctuations with respect to the Favre (density-weighted) average. The mathematical definition of 
$\mathcal {I}_{\langle F , k_z \rangle }$ does not depend on the ordering of the modes within the triad, and hence this expression only measures the energy exchange and not the direction of the transfer. In other words, the designation 
$\langle F_1 , k_{z,1} \rangle + \langle F_{2} , k_{z,2} \rangle \Rightarrow \langle F , k_z \rangle$ in figure 10(i–l) merely identifies the triad and bears no directional significance, although the outcome of the exchange can be gleaned, with caution, by simultaneously considering the spectra of the individual waves (figure 10c–f).
(a,b) The optimal protruding roughness at the end of iteration # 4, with ‘x’ marking the synchronization Reynolds numbers of modes 
$\langle 110,0 \rangle$ and 
$\langle 100,0 \rangle$, respectively. (c–h) Spectral energy, 
$\mathcal {E}_{\langle F , k_z \rangle }$, for selected instability modes versus streamwise coordinate. (i–l) Modal nonlinear energy-transfer coefficient, computed for key triad interactions. Grey lines are reference curves for the flat-plate case; black lines are the results of the boundary layer over the optimal roughness.
To establish the background against which the influence of roughness is examined, we summarize the transition mechanism in the flat-plate case with the aid of the spectra in figure 10(c–h). The three most energetic inflow modes, 
$\{\langle 110 , 0 \rangle, \langle 100 , 0 \rangle, \langle 40 , 3 \rangle \}$, initially amplify and then spur other waves via nonlinear interactions. The pair of second modes participate in 
$\mathcal {I}^{(1)}$ (figure 10i) and generate 
$\langle 10 , 0 \rangle$. Another triad, 
$\mathcal {I}^{(2)}$ in figure 10( j), activates mode 
$\langle 70 , 3 \rangle$ which is visible as oblique structures in figure 9(a) near 
$\sqrt {{Re}_x} \approx 2300$. The last two interactions, 
$\mathcal {I}^{(3)}$ and 
$\mathcal {I}^{(4)}$ in figure 10(k–l), involve both the inflow and newly formed waves, and spur the formation of mode 
$\langle 30 , 3 \rangle$. This mode amplifies faster than any other wave, forms the 
$\varLambda$-shaped structures in figure 9(a) in the range 
$2400 < \sqrt {{Re}_x} < 2550$ and is the ultimate cause of breakdown to turbulence.
Comparing the reference case with the optimal roughness in figure 10(c–l), both modes 
$\langle 110,0\rangle$ and 
$\langle 100,0\rangle$ are attenuated in the latter configuration as the roughness element is approached (see figure 10c,d). The effect is more pronounced for mode 
$\langle 100,0\rangle$ whose synchronization point is much closer to the roughness. The subsequent nonlinear interactions that involve these two instabilities, especially interactions 
$\mathcal {I}^{(1)}$ and 
$\mathcal {I}^{(3)}$ involving mode 
$\langle 100,0\rangle$, are also mitigated. The outcome is an appreciable weakening of the nonlinearly generated modes that cause breakdown to turbulence.
Another observation from figure 10 is that the first-mode wave 
$\langle 40,3\rangle$ is negligibly destabilized as a result of the introduced roughness. This trend is different from previous efforts (see e.g. Bountin et al. Reference Bountin, Chimitov, Maslov, Novikov, Egorov, Fedorov and Utyuzhnikov2013; Fong et al. Reference Fong, Wang, Huang, Zhong, McKiernan, Fisher and Schneider2015) and is due to the smaller height of the roughness in our simulation. For comparison, the height of the roughness in the experiments by Fong et al. (Reference Fong, Wang, Huang, Zhong, McKiernan, Fisher and Schneider2015) and by Bountin et al. (Reference Bountin, Chimitov, Maslov, Novikov, Egorov, Fedorov and Utyuzhnikov2013) was approximately equal to 
$0.5 \delta _{X_r}$, while in our analysis, the height is 
$H_r = 0.23 \delta _{X_r}$. Mode 
$\langle 40,3\rangle$ is appreciable between the relative sonic line and the boundary-layer edge, and this region is minimally affected by the presence of the roughness which is positioned below the sonic line.
The results presented in figure 10 highlight that roughness can have a direct, local effect on the instability waves. For example, in figure 10(c,d), the growth rates of modes 
$\langle 100 , 0 \rangle$ and 
$\langle 110 , 0 \rangle$ are significantly altered in the region 
$2000 < \sqrt {Re_x} < 2100$, as these instabilities approach and travel over the roughness. As noted in § 3.1, the dominant local roughness effects in our configuration are the roughness-induced mean-flow distortion and non-parallel effects in regions where 
${\rm d} h_r / {{\rm d}\kern0.06em x}$ is large. In particular, the non-parallel effects of the optimal roughness (which is taller, more slender and more abrupt) is more pronounced compared to roughness X2p1. In addition to these local effects, the influence of roughness on transition precursors can persist downstream. For example, figure 10(e) shows that the growth rate of mode 
$\langle 40 , 3 \rangle$ in the roughed-wall case reduces dramatically for 
$\sqrt {Re_x} > 2200$. This change can be traced back to the local effects on the amplitude of modes 
$\langle 100 , 0 \rangle$ and 
$\langle 110 , 0 \rangle$: stabilization of these waves by the roughness mitigates the subsequent nonlinear interactions (figure 10i–l). Specifically, over a flat plate, mode 
$\langle 40 , 3 \rangle$ is activated near 
$\sqrt {Re_x} \approx 2200$ in a triad interaction 
$\mathcal {I}^{(4)}$ which involves 
$\langle 10 , 0 \rangle$ and 
$\langle 30 , 3 \rangle$ (figure 10l). In the rough-wall case, this triad interaction is delayed and muted, due to a delay in the generation of modes 
$\langle 10 , 0 \rangle$ and 
$\langle 30 , 3 \rangle$ from the preceding interactions. In light of these observations, we turn our attention to the local effect of the optimal roughness on the second-mode instabilities 
$\langle 100 , 0 \rangle$ and 
$\langle 110 , 0 \rangle$.
3.2.3. Local effect of optimal roughness on second-mode instabilities
Figures 11 and 12 provide flow visualizations to aid the physical interpretation of how the roughness modifies the near-wall region and alters the stability behaviour of the second modes. These figures correspond to the simulation of the boundary layer over the optimal roughness after the final iteration of the EnVar procedure. Figure 11 shows contours of
\begin{equation} \mathcal{D} = \tanh{ ( \xi_0 \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} ) } ,\end{equation}
where 
$\xi _0$ is a constant that adjusts the background colour. A few additional lines are marked on the figure: the boundary-layer edge is identified as the location where 
${\bar {u} = 0.99 u_\infty}$; the relative sonic lines are defined as the heights at which 
$c - \bar {u} = a$, where 
$c$ is the phase speed of either mode 
$\langle 100 , 0 \rangle$ or 
${\langle 110 , 0 \rangle }$. Below the relative sonic line, the alternating compressing and expanding acoustic waves travel supersonically relative to the mean flow. Inviscid, linear perturbation theory predicts that the effect of the wall roughness propagates to infinity with constant strength along the lines 
$x - [M^2_\infty - 1]^{{1}/{2}} y = \textrm {const}.$, which represent the local Mach lines (Liepmann & Roshko Reference Liepmann and Roshko2001). In figure 11, several Mach lines are plotted with slope 
$[M_{y = 35}^2 - 1]^{-{1}/{2}}$, where 
$M_{y = 35}$ is the Mach number at 
$y = 35$. These lines agree with the prediction from the theory and qualitatively capture the free stream compression zone upstream the roughness element. On top of the roughness, the slope of the compression zone progressively deviates from the theory and is followed by an expansion zone. The deviation from theory is expected in light of its restrictive assumptions, for example, it assumes that 
$\theta _{max} \sqrt {M_\infty ^2 - 1} \ll 1$, where 
$\theta _{max}$ is the maximum inclination angle of the roughness (Liepmann & Roshko Reference Liepmann and Roshko2001), while the geometrical features of the optimal roughness in our simulation yield 
$\theta _{max} \sqrt {M_\infty ^2 - 1} = 0.5 H_r L_r \sqrt {M_\infty ^2 - 1} \approx 0.6$. Figure 11 also shows that in the zone upstream of the roughness, the near-wall relative supersonic region becomes thicker whereas above the protrusion, this region becomes thinner. The trapped acoustic waves reflect back and forth at the wall and the relative sonic line. At the latter location, the waves change from compression to expansion and vice versa. As the height of the relative supersonic zone changes, the periods between consecutive reflections at the relative sonic line change, which disrupts the amplification of the instability mode.
Contours of 
$\mathcal {D}$ defined by (3.3). Dark colours are expanded regions (
$\mathcal {D} > 0$, 
$\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u} > 0$) and light colours mark compressed zones (
$\mathcal {D} < 0$, 
$\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u} < 0$). Red dash–dotted line is the boundary-layer edge, 
$\delta _{99}$; yellow lines are the local Mach lines; and white dotted and blue dashed lines mark the relative sonic line, where 
$c - \bar {u} = a$, for modes 
$\langle 100 , 0 \rangle$ and 
$\langle 110 , 0 \rangle$, respectively.
Snapshots of 
$\mathcal {D}$ and pressure fluctuations, 
$p^{\prime }$, at (a) 
$t = 0$, (b) 
$t = 418.6$ and (c) 
$t = 575.4$. Red and yellow lines are 
$p^{\prime } > 0$ and 
$p^{\prime } < 0$, respectively. Pressure fluctuations iso-lines are 
$\{0.5, 1, 2\} \times 10^{-4}$ (from outer to inner lines).
Figure 12 shows three instantaneous fields, which highlight the pressure fluctuations, 
$p^\prime = p - \bar {p}$ inside the boundary layer. While the visualized 
$p^\prime$ is the outcome of a superposition of different instability waves, rather than a single one, the figure is helpful in explaining the behaviour of the instabilities in response to the roughness element. Two separate regions in the wall-normal profile of 
$p^\prime$ can be identified in this figure: (i) near-wall (below the relative sonic line) 
$p^\prime$ contours are caused by the acoustic component of the instability waves and (ii) 
$p^\prime$ above the relative sonic line where the vortical and thermal components of the instability wave are prominent. Figure 12(a–c) track a series of compression–expansion–compression in 
$p^\prime$, and the three figures correspond to three time instances, {t1, t2, t3}. As we discuss the wall-normal profile of 
$p^\prime$, our focus will be on the phase relation below (region denoted 
$L$) and above (region denoted 
$U$) the relative sonic line. At t1 (figure 12a), the wall-normal profiles of 
$p'$ as we traverse from t1L to t1U have a fixed phase relation that is initially preserved with downstream distance. In this configuration, they are phase tuned. At time t2 (figure 12b), the 
$p^\prime$ packet has travelled to the marked location and the change in phase along the wall-normal profile of 
$p'$ across the sonic line has changed appreciably, by approximately 
$180^\circ$. We will refer to this process as phase detuning. At time t3 (figure 12c), the wall-normal profile of 
$p'$, crossing from t3L to t3U, recovers the original phase relation. Another interesting observation from figure 12(a and b) is the propagation of pressure fluctuations, 
$p^\prime$, along the Mach lines emanating above the roughness into the free stream. This process is absent at t3 (figure 12c), and has a low normalized frequency of the order of 
$F=10$. The most energetic mode at this frequency is 
$\langle 10 , 0 \rangle$, which is generated by nonlinear interaction of modes 
$\langle 110 , 0 \rangle$ and 
$\langle 100 , 0 \rangle$, and whose wavelength is 
$\lambda _{\langle 10 , 0 \rangle } / 2 W_r = 2.5$ relative to the roughness streamwise extent.
We now return to the change of phase observed in the near-wall pressure fields in figure 12. While the interpretation in terms of detuning of the individual instability waves is plausible, the same visual pattern can be caused by other effects, for example, dispersion of the waves that comprise the plotted fluctuating pressure field. To provide a more precise assessment of the influence of the roughness on the key instabilities, we perform a Fourier transform of the pressure signal in time and the span,
\begin{equation} p(\boldsymbol{x},t) = \sum_{F,k_z} \hat{p}_{\langle F , k_z \rangle }(x,y) \exp\left[{-\textrm{i} \left( \frac{\sqrt{Re_{x_0}} F }{10^6} t + \frac{2{\rm \pi} k_z}{L_z} z \right) }\right],\end{equation}
where 
$\hat {p} ( x , y )$ is the complex Fourier coefficient whose phase we will denote as 
$\phi _{\hat {p}}$. Inspired by stability theory, we introduce the ansatz
\begin{equation} \hat{p}_{\langle F , k_z \rangle } \equiv \check{p}_{\langle F , k_z \rangle } \exp \left( \int_{x_0}^x \alpha_{\langle F , k_z \rangle }\, {{\rm d}\kern0.06em x} \right),\end{equation}
where 
$\alpha (x) = \alpha _r + i \alpha _i$ is a complex streamwise wavenumber, which is evaluated from
\begin{equation} \alpha_{\langle F , k_z \rangle } = \frac{0.5}{ \mathcal{E}_{\langle F , k_z \rangle } } \int_{0}^{L_y} \left[ \bar{\rho} \left\{ \boldsymbol{\hat{u}}^* \frac{\partial \boldsymbol{\hat{u}}}{\partial x} \right\} + \frac{\bar{\rho} \bar{T}}{(\gamma-1) \gamma M_\infty^2} \left\{ \frac{ \hat{\rho}^* }{\bar{\rho}^2} \frac{\partial \hat{\rho}}{\partial x} + \frac{ \hat{T}^* }{\bar{T}^2} \frac{\partial \hat{T}}{\partial x} \right\} \right]_{\langle F , k_z \rangle }\,{{\rm d} y}. \end{equation}
The mode shape is therefore 
$\check {p}(x,y)$ and its phase will be denoted 
$\phi _{\check {p}}$.
The phases of instability modes 
$\langle 110 , 0 \rangle$ and 
$\langle 100 , 0 \rangle$ are reported in figure 13. The contour plots contrast the behaviour of 
$\phi _{\hat {p}}$ in the reference flat-plate configuration (figure 13a,b) and as the instability waves approach the optimal roughness (figure 13c,d). The results highlight the distinction between the modal pressure fluctuations that are below and above the relative sonic line, the former travelling supersonically relative to the mean flow. In figure 13(a,b), the phase change across the relative sonic line is maintained at a nearly constant value along the flat plate in the shown region in the figures. In contrast, figure 13(c,d) clearly show that the relative phase across the sonic line is disrupted as each instability wave approaches the roughness. Most importantly, the phases of the pressure modes in the trapped supersonic regions, below the relative-sonic lines, are rapidly altered. For mode 
$\langle 110 , 0 \rangle$ in figure 13(c), this change occurs at 
$\sqrt {Re_x} \approx 2000$ which, according to figure 10(c), is the location where the amplification of this instability wave is abruptly halted and it starts to decay. Similarly, for mode 
$\langle 100 , 0 \rangle$ in figure 13(d), detuning of the pressure fluctuations takes place at 
$\sqrt {Re_x} \approx 2020$, which also corresponds to the location where this instability wave reaches a local maximum amplitude after which it decays in figure 10(d). Figure 13(e, f) provide a more quantitative picture. We report the phases of the pressure in the modal representation (3.5), or 
$\phi _{\check {p}}$, along the wall. The influence of the roughness relative to the flat-plate case is evident, as well as the correlation between the changes in the phases and in the modal spectra from figure 10(c,d).
(a–d) Contours of the phase of the Fourier representation of pressure, 
$\phi _{\hat {p}}$: (a,b) reference flat-plate configuration; (c,d) optimized-roughness case. (e, f) Phase of mode-shape of pressure signal, 
$\phi _{\check {p}}$, at the wall; grey lines are reference curves for the flat-plate case, while the black lines are the results of the boundary layer over the optimal roughness. Panels (a,c,e) show mode 
$\langle 110 , 0 \rangle$ and (b,d, f) show mode 
$\langle 100 , 0 \rangle$.
The results presented in § 3.2 support our original physical argument regarding the effect of roughness on the stability behaviour of second-mode instabilities. As the thickness of the near-wall, relative supersonic region changes abruptly, the instability waves are altered: their growth rate is modulated, the phase of their acoustic component is ‘scrambled’ and their net amplification is reduced relative to the flow over a flat plate. As a result, the nonlinear interactions among key instability waves, downstream of the optimal roughness, are effectively weakened and transition to turbulence is nearly suppressed.
4. Summary
The capacity of isolated roughness to delay laminar-to-turbulent transition in a Mach 4.5 boundary layer is examined. The study is conducted using DNS of boundary layers over a flat plate, with isolated roughness elements that have various geometrical parameters mounted on the surface. The oncoming disturbance that interacts with the boundary layer in the simulations corresponds to the nonlinearly most dangerous disturbance for the prescribed energy level (Jahanbakhshi & Zaki Reference Jahanbakhshi and Zaki2019), and is comprised primarily of two planar second-mode instabilities and an oblique first-mode instability.
The roughness shape was parametrized to provide control over the location, streamwise wavenumber, height, width and abruptness of the roughness. The influence of these parameters was examined, in particular, their impact on the location of laminar-to-turbulence transition relative to the reference smooth-wall configuration. The effects of location and streamwise wavenumber were explored first, and provided the initial estimate of a roughness design that was adopted in subsequent optimization where the height, width and abruptness of the roughness were optimized to achieve maximum transition delay in our computational domain. The constrained optimization was performed using an ensemble-variational approach, in which a cost function is defined in terms of the skin-friction coefficient and is minimized to ensure the latest possible transition to turbulence. The optimal roughness that was identified was able to maintain a laminar state throughout the computational domain.
The roughness elements minimally affected the first-mode instability that was part of the inflow disturbance, while the second-mode waves and their nonlinear interactions were appreciably attenuated. Whether the second-mode waves were stabilized or destabilized depended on the location and geometrical features of the roughness. Specifically, stability is affected by (i) the relative position of the roughness and the synchronization of the slow and fast modes and (ii) the height of the near-wall region where the second-mode instability waves travel supersonically relative to the mean flow. Pre-synchronization, altering the supersonic region by a protruding roughness destabilizes the second-mode waves, while post-synchronization, the net effect is stabilizing. The change in stability is due to the shift in the phase of the trapped acoustic waves relative to the harmonic vortical and thermal waves that are beyond the relative sonic line. The net outcome is an effective delay of the instability growth, mitigation of nonlinear interaction and ultimately transition delay. The optimized roughness was more slender, slightly taller and more abrupt than the initial design. This roughness was more effective, entirely eliminating transition from the computational domain.
Depending on their shape and locations, roughness elements can have stabilizing, destabilizing or neutral net effect on the amplification of instability waves. As a result, the impact on the nonlinear stages of transition and the onset of turbulence are difficult to anticipate. Our EnVar framework provides objective guarantees that the optimized roughness delays transition to turbulence at design conditions, and we did not observe any undesirable or unexpected destabilization of originally benign disturbances. The optimized roughness does not, however, guarantee transition delay for other disturbances, away from the design conditions, e.g. since new transition mechanisms may be active. For the purpose of this work, we tested a broadband inflow spectrum of instability waves at sufficiently high amplitude to trigger transition within the computational domain over the flat plate (see Appendix B), and still observed transition delay when the optimal roughness from § 3.2 was introduced. Despite this additional test, there remains no guarantee that transition would be delayed for different inflow disturbance spectra. Essentially, the optimized roughness should not be interpreted as optimal for all inflow disturbances, but rather for the inflow condition that was adopted in the optimization procedure. By considering the nonlinearly most dangerous disturbance, our goal was to demonstrate that the optimization of the roughness can be effective, even in the most aggressive scenario.
Acknowledgements
The authors are grateful to Professor J. Larsson for sharing the Hybrid code.
Funding
This work was supported in part by the US Air Force Office of Scientific Research (grant no. FA9550-19-1-0230) and by the Office of Naval Research (grant nos. N00014-17-1-2339, N00014-21-1-2148). Computational resources were provided in part by the Maryland Advanced Research Computing Center (MARCC) and the AI.Panther HPC cluster at Florida Institute of Technology (funded by the National Science Foundation major Research Implementation grant no. MRI-2016818).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Flow over complex geometries
This appendix provides a brief description of the cut-stencil method implemented in the code Hybrid (Johnsen et al. Reference Johnsen2010) and a validation test. As shown in the schematic of figure 14, the cut-stencil is a sharp interface method where the discretization of the governing equations is modified near the solid body to explicitly enforce the no-slip boundary conditions at the fluid–solid interface. The precise location for the solid/fluid interface is adopted to ensure local, and therefore also global, conservation. The present implementation is similar to that by Greene et al. (Reference Greene, Eldredge, Zhong and Kim2016). Therefore, herein, we introduce the cut-stencil method very briefly and refer the readers to Greene et al. (Reference Greene, Eldredge, Zhong and Kim2016) for a more detailed discussion of the method. We then provide a sample validation study by comparing the results from our code with the published data for a flow simulation that is relevant to the topic of current work.
(a) Schematic of the grid near the immersed body. (b) Point classification for the cut-stencil implementation in code Hybrid.
A.1. A brief description of the cut-stencil
A 2-D schematic of the implementation is shown in figure 14. All grid points inside the solid body are discarded, and the cut-cell discretization generates four classes of fluid points that are treated separately.
(i) For interior fluid grid points that are far from the body, the Navier–Stokes equations are discretized using the interior schemes.
(ii) For the first three fluid grid points near the solid body, sixth-, fourth- and second-order schemes are used to discretize the Navier–Stokes equations on these points.
(iii) Points that are denoted as irregular fluid points in the figure are excluded from the computations along the directions where they reside too close to the solid (within 25 % of the local grid spacing). However, in the other directions, the points may be included in the computations of other neighbouring fluid points. Therefore, flow variables are required at these points and are computed via an interpolation scheme from their neighbours including the boundary points.
(iv) Boundary points are locations where the body surface intersects the grid and where the flow variables are also needed. In current implementation, the velocity of the immersed object is computed using the no-slip wall condition and the temperature is computed from the thermal boundary condition. The remaining parameters are the velocity gradients and the pressure (or density). Velocity gradients are computed using a fourth-order one-side finite-difference scheme in the direction normal the solid body. The pressure is extrapolated from the neighbouring fluid points and the equation of state for an ideal gas is used to compute the density.
The solver requires interpolation/extrapolation of flow quantities onto the irregular and boundary points, and points within the fluid domain that are not on the grid. For example, for computing 
$\partial F / \partial \eta$ at the boundary points, the value of 
$F$ is required at four points equally spaced from the surface in the direction of 
$\eta$. Similar to Greene et al. (Reference Greene, Eldredge, Zhong and Kim2016), the interpolation/extrapolation function is a second-degree polynomial in 
$x$, 
$y$ and 
$z$, while the coefficients of the polynomial are computed using least-squares fitting of the neighbouring points.
A.2. A sample validation study
The current implementation of the cut-stencil algorithm has been validated extensively against published data. The comparisons included multiple configurations for 2-D and 3-D solid bodies, and a range of Mach numbers from 0.2 to 4.8. Here, we provide the results from the validation case most relevant to the present effort, namely a high-Mach-number boundary layer over a 2-D isolated roughness.
We compare results from our cut-stencil implementation with the published data by Greene et al. (Reference Greene, Eldredge, Zhong and Kim2016) from their body-fitted curvilinear solver. The computational domain starts at 
$\sqrt {Re_{x_{0}}} = 762$ and the free stream Mach number is 
$Ma_{\infty } = 4.8$. The Blasius length scale at 
$x_0$ and free stream quantities are adopted as the reference scales. The geometry of the roughness is a smooth protrusion which is given by two hyperbolic tangent functions as
\begin{equation} h(x) = \frac{h_r}{2} [ \tanh( L_r [x - x_r + W_r] ) - \tanh ( L_r [x - x_r - W_r] )] , \end{equation}
where 
$h_r = 13.115$, 
$x_r = 1961.30$, 
$L_r = 0.152497592$ and 
$W_r = 26.23$ are respectively the non-dimensionalized height, centre location, abruptness and width. More information about the simulation parameters can be found in table 4 of Greene et al. (Reference Greene, Eldredge, Zhong and Kim2016).
Comparisons of our simulation results with those of Greene et al. (Reference Greene, Eldredge, Zhong and Kim2016) are provided in figures 15 and 16. In these figures, 
$\delta \approx 24$ is the undisturbed boundary-layer thickness at 
$x_r$. As demonstrated by figure 15, the velocity, density and pressure contours from the cut-stencil implementation and from the reference body-fitted simulation are indistinguishable. For a more detailed view, profiles of the streamwise velocity and pressure were extracted at three locations and are reported in figure 16. Agreement, in particular in the near-wall region, is evident, which verifies the accuracy of the implementation.
Contours of (a) wall-normal velocity, (b) density and (c) pressure. The contour colours and thin grey lines correspond to the current simulation and red dashed lines are the published results by Greene et al. (Reference Greene, Eldredge, Zhong and Kim2016).
(a,c,e) Streamwise velocity and (b,d, f) pressure profiles at three streamwise positions. Blue lines correspond to current simulation and red dashed lines are the published results by Greene et al. (Reference Greene, Eldredge, Zhong and Kim2016). Horizontal black lines mark the height of roughness.
Appendix B. Effect of optimized roughness on transition due to broadband spectrum
Additional simulations were carried out, where we tested the performance of the optimized roughness obtained in § 3.2 for an off-design condition. Specifically, rather than considering the nonlinearly most dangerous inflow disturbance, we adopted a broadband inflow spectrum where the total disturbance energy was distributed equally among the 400 instability waves that span the frequency and wavenumber range in figure 3(a), with randomly assigned phases.
When the total energy of the inflow instability waves was the same as in the main body of this paper, i.e. 
$\sum _{F,k_z} \mathcal {E}_{\langle F , k_z \rangle } = 2 \times 10^{-5}$, transition did not take place within the computational domain. We therefore increased the inflow total disturbance energy by two orders of magnitude, 
$\sum _{F,k_z} \mathcal {E}_{\langle F , k_z \rangle } = 10^{-3}$. Under these conditions, transition takes place within the computational domain. Results that compare flow over a flat plate and the optimized roughness from § 3.2 are shown in figure 17, where we report the skin-friction coefficient and the downstream evolution of the total disturbance energy. The results demonstrate that the roughness remains stabilizing for this configuration, shifting the onset of turbulence downstream. We stress, however, that the influence of the roughness could have been destabilizing since there is no performance guarantee when the inflow disturbance is far from the design condition. For example, the roughness may still promote transition at different (higher or lower) energy levels in the case of broadband inflow spectrum, or in the case of an entirely different inflow spectrum that interacts with the roughness in a destabilizing fashion that is not featured at design conditions.
Influence of optimal roughness on transition due to high-energy, broadband inflow disturbance spectrum: (a) skin-friction coefficient and (b) total disturbance energy versus downstream Reynolds number. Grey lines are reference curves for flat-plate case; black lines are the results of the boundary layer over the optimal roughness.
Open access
