3. Flow visualisations

We begin by investigating the instantaneous and mean flow in different planes, particularly focusing on the smooth-to-rough transition region.

Figure 3. Instantaneous density field $\rho /\rho _{\infty }$ visualised in an $x$ - $y$ slice taken at $z/k = 1.5$ (centre of roughness elements) for (a) CB_A, (b) CB_S, (c) CB_R, (d) DM_S. In black: contour of the averaged velocity field where $\bar {u} =0.99 u_\infty$ .

Figure 3 shows the instantaneous density field $\rho /\rho _\infty$ in a longitudinal plane, featuring the interaction of the turbulent boundary layer and the shock wave emanated at the onset of roughness. We begin by noting that the shock wave is not directly emanated by the first row of roughness elements, but rather by the abrupt growth of the boundary layer, which has to adjust to the new surface condition, as previously noted by Cogo et al. (Reference Cogo, Modesti, Bernardini and Picano2025). After the shock wave, all cases clearly display a pattern of compression and expansion waves emanated towards the free stream, with the CB_S (figure 3 b) and CB_R (figure 3 c) cases being qualitatively more disturbed. This is clear both from the thicker dark red region, denoting an increased superposition of individual compression waves, and by the contour of $\bar {u} =0.99 u_\infty$ , shown in black, representative of the local boundary layer thickness $\delta _{99}$ . In fact, while all cases share an overshoot of the thickness $\delta _{99}$ right after the oblique shock wave, this region is more extended for the CB_S and CB_R cases. This observation is confirmed by computing the average shock angle $\bar {\beta }$ , reported in table 2, which is the highest for the CB_R case, followed by CB_S, DM_S and CB_A. Although differences in $\bar {\beta }$ between all cases are less than $1^\circ$ , this observation is in line with trends noted in various statistics discussed in the next sections. Here, we also note that $\bar {\beta }$ is consistently higher than the corresponding Mach angle at $M_\infty =2$ , which is $30^\circ$ , which confirms the suitability of the term ‘shock wave’ in lieu of ‘Mach wave’.

Table 2. Shock angle $\beta$ computed by linear fitting of the pressure-gradient contours $\partial \bar {p}/\partial x=0$

Figures 4 and 5 show density contours wall-parallel and cross-stream planes, respectively. The former is taken at a wall-normal location corresponding to half of the roughness height $k$ and clearly shows the breakdown of near-wall turbulent streaks in the smooth-wall region (before $x/\delta _{\textit{in}}$ ), due to the onset of roughness. Right before the first row of elements, density spikes towards the free-stream reference (dark red regions), before settling to a much lower value in the gaps between the elements. Regions with high density are more evident with CB_S and CB_R, which feature more intense shock waves.

Figure 4. Instantaneous density field $\rho /\rho _{\infty }$ visualised in an $x$ - $z$ slice taken at $y/k=0.5$ for (a) CB_A, (b) CB_S, (c) CB_R, (d) DM_S.

Figure 5 shows density contours in a cross-stream plane at $x/\delta _{\textit{in}}=60$ , shortly after the smooth-to-rough transition. Here, the spanwise distribution of the shock wave is visible, and its generation as the superposition of multiple spherical compression waves is apparent. We can qualitatively note that the CB_R case exhibits the highest intensity in the instantaneous value of $\rho /\rho _\infty$ due to the merging of stronger compression waves.

Figure 5. Instantaneous density field $\rho /\rho _{\infty }$ visualised in an $z$ - $y$ slice taken at $x/\delta _{\textit{in}}=60$ for (a) CB_A, (b) CB_S, (c) CB_R, (d) DM_S.

Although the presence of a first oblique shock wave is clearly observed from the instantaneous fields, there is no evidence of subsequent shock waves emanating from individual roughness elements, suggesting that the sonic line is above the roughness crest. This is in contrast to previous experimental studies, for example Kocher et al. (Reference Kocher, Kreth, Schmisseur, LaLonde and Combs2022), who noted through averaged schlieren imaging that each element generated individual shock waves and expansion fans propagating into the free stream. Figure 6 shows the numerical schlieren obtained for the present database, with an indication of the sonic line in dashed red. Here, we note that for all cases there are subsequent compressions and expansions emanating from the first few rows of roughness elements, which, however, do not appear to extend towards the free stream and more importantly severely decrease in intensity downstream. This observation is in agreement with the fact that the element crest lies on average below the sonic line for all cases (shown in red), and with the local Mach number calculated at the roughness crest, reported in table 3.

Table 3. Boundary layer properties at the selected stations. Here, $\textit{Re}_\theta = \rho _\infty u_\infty \theta / \mu _\infty$ , with $\theta = \int _0^e (\rho u)/(\rho _e u_e) (1-u/u_e)$ the momentum thickness, $U_h$ the mean streamwise velocity evaluated at the roughness crest $y=k$ and $M_h$ the Mach number evaluated at the same location using the velocity $U_h$ and the local speed of sound.

Figure 6. Numerical schlieren obtained from the averaged density field $\bar {\rho }/\rho _{\infty }$ visualised in an $x$ - $y$ slice for (a) CB_A, (b) CB_S, (c) CB_R, (d) DM_S. In dashed red: contour of the sonic line where the average local Mach number is one.

4. Mean boundary layer development

In this section, we consider the streamwise development of the boundary layer and we analyse the skin friction coefficient $C_{\!f}=\tau _w/(1/2\rho _{\infty } u_{\infty }^2)$ , the boundary layer thickness $\delta _{99}$ and the friction Reynolds number $\textit{Re}_{\tau }$ , presented in figure 7. For all cases in our database, these statistics are equal in the smooth region, thus all plots start from $x/\delta _{\textit{in}}=20$ to enhance readability of the aft portion of the domain. We also limit the upper bound of $x$ to $x/\delta _{\textit{in}}=140$ , neglecting the subsequent adjustment of the flow to a smooth region placed just before the end of the domain.

Figure 7(a) shows the streamwise evolution of $C_{\!f}$ , featuring a sharp increase at $x=55 \delta _{\textit{in}}$ for all rough cases, corresponding to the location where the roughness starts. After the initial overshoot, all cases share a relatively similar trend, although with different intensities, as each roughness shape leads to a different related added drag. In particular, rotated cubes in CB_R clearly show the highest value of $C_{\!f}$ both in the initial overshoot and in the subsequent recovery, followed by the staggered cubes, CB_S, staggered diamonds, DM_S, and aligned cubes, CB_A.

Flow case CB_R also exhibits the thicker boundary layer in terms of $\delta _{99}$ both in terms of initial overshoot and subsequent growth, followed by staggered cubes in CB_S, as shown in figure 7(b). Staggered diamonds in DM_S and aligned cubes in CB_A share a relative lower increase of $\delta _{99}$ , although it is still visible when compared with the smooth-wall reference. This aspect was previously noted by Cogo et al. (Reference Cogo, Modesti, Bernardini and Picano2025) for case CB_A, who reported an increase factor of 1.4 comparing $\delta _{99}$ after the overshoot with a reference $\delta _{99,\textit{ref}}$ computed at $x/\delta _{\textit{in}}=45$ . This fact is even more evident here for staggered and rotated cubes, which show an increase factor of approximately 1.5, which is then followed by a steeper streamwise evolution for the latter case (CB_R). Actually, for these cases the increase in $\delta _{99}$ seems to be associated with a longer streamwise extent of the overshoot region, which appears to approximately scale with the skin friction coefficient $C_{\!f}$ .

Figure 7. Mean streamwise profiles of (a) skin friction coefficient $C_{\!f}=\tau _w/(1/2\rho _{\infty } u_{\infty }^2)$ , (b) boundary layer thickness $\delta _{99}$ based on the $99\,\%$ velocity $u_{99}=0.99 u_{\infty }$ , (c) friction Reynolds number $\textit{Re}_{\tau }=\delta _{99}/\delta _{\nu }$ as a function of the streamwise coordinate $x/ \delta _{\textit{in}}$ . A smooth-wall reference case, SM_M2, from Cogo et al. (Reference Cogo, Modesti, Bernardini and Picano2025) at the same Mach and friction Reynolds numbers is also reported.

The smooth-to-rough transition and subsequent development have been a topic of major interest in many applications, not only for the prediction of the boundary layer thickness but also for the generation of an IBL. This is the region of the boundary layer that is directly influenced by the rough surface, evolving in the streamwise direction and eventually merging with the external boundary layer (Elliott Reference Elliott1958; Rouhi et al. Reference Rouhi, Chung and Hutchins2019). Multiple definitions have been proposed in recent years to detect the upper bound of the IBL, but we found that only some of them are accurate in the supersonic regime, and proposed a novel definition based on the mean wall-normal fluctuation intensity $(\bar {\rho } \widetilde {v^{\prime \prime 2}})/ \rho _{\infty } u_{\infty }^2$ (Cogo et al. Reference Cogo, Modesti, Bernardini and Picano2025). This control variable was found to be effective in both subsonic and supersonic regimes for being directly related to the influence of roughness on turbulence, which greatly enhances wall-normal velocity fluctuations in the roughness sublayer, while being weakly affected by compressibility effects, such as shocklets, which mainly influence the streamwise velocity component. In the present work, we employ this definition for all present cases without individual tuning, and we compare the evolution of the IBL thickness $\delta _{\textit{IBL}}$ with the reference value $\delta _{99}$ .

Figure 8 shows the contours of $(\bar {\rho } \widetilde {v^{\prime \prime 2}})/ \rho _{\infty } u_{\infty }^2$ in a wall-normal plane for each case, where the generation and evolution of the IBL are clearly visible in yellow and initially marked by the dashed black line generated using the algorithm proposed by Cogo et al. (Reference Cogo, Modesti, Bernardini and Picano2025). In addition, the streamwise growth of the external boundary layer is reported with solid black lines, in terms of $\delta _{99}$ . Here, the cases with the most pronounced initial shock wave (see § 3), CB_S and CB_R, display a more intense wall-normal velocity fluctuations compared with the other two cases, CB_A and DM_S. Nevertheless, the intensity of wall-normal velocity fluctuations progressively dampens downstream for all cases.

Figure 8. Contours of $(\bar {\rho } \widetilde {v^{\prime \prime 2}})/ \rho _{\infty } u_{\infty }^2$ in the averaged wall-normal plane for (a) CB_A, (b) CB_S, (c) CB_R, (d) DM_S. Dashed black lines represent the detected upper bound of the IBL using the algorithm proposed by Cogo et al. (Reference Cogo, Modesti, Bernardini and Picano2025). Solid black lines indicate the evolution of $\delta _{99}$ .

Figure 9 reports the streamwise evolution of $\delta _{\textit{IBL}}$ as a function of the distance from the onset of roughness $(x-x_{\textit{ref}})$ , where $x_{\textit{ref}}=55\delta _{\textit{in}}$ . Figure 9(a) shows both the $\delta _{\textit{IBL}}$ and $\delta _{99}$ profiles as function of $(x-x_{\textit{ref}})$ , all rescaled with $\delta _{99,\textit{ref}}$ , a reference thickness of the incoming smoot- wall boundary layer selected right before the rough portion, at $x/\delta _{\textit{in}}=50$ . The streamwise growth of the IBL is calculated using the algorithm based on the wall-normal velocity fluctuations, as well as the fitting proposed in Cogo et al. (Reference Cogo, Modesti, Bernardini and Picano2025).

Figure 9. Streamwise growth of $\delta _{\textit{IBL}}$ in absolute (a) and relative (b) terms, compared with $\delta _{99}$ . Power-law extrapolations are present for each case.

The detected evolution of $\delta _{\textit{IBL}}$ , figure 9(a), approximately follows the evolution of $\delta _{99}$ discussed in § 4, with the CB_S and CB_R cases overlying CB_A and DM_S. This trend is not visible in the extrapolated power laws, figure 9(b), which, however, serve only as visual indication of the expected merging points of the two layers, and are very sensitive to small errors.

Considering the relative growth of $\delta _{\textit{IBL}}$ and $\delta _{99}$ , shown in figure 9(b), all cases are expected to approximately recover an equilibrium boundary layer in the range of $40$ – $50$ times the incoming smooth-wall boundary layer thickness $\delta _{99,\textit{ref}}$ . This is an important reference measure to inform future computational and experimental studies, although we cannot comment on Reynolds and Mach number effects.

We remark that, once an equilibrium is observed between the boundary layer and the rough region, and the separation of scales is enough to observe the fully rough regime, we do not expect this condition to hold in the long plate limit $x \rightarrow \infty$ , since the flow will progressively return to an hydraulically smooth case because the ratio $\delta _{99}/k$ increases (Pullin, Hutchins & Chung Reference Pullin, Hutchins and Chung2017). This is in contrast to internal flows, where the ratio $\delta /k$ is fixed and the fully rough regime is expected to hold from a certain location onward. In turn, the fully rough regime on a flat plate can be achieved at a given streamwise location $x$ when the unit Reynolds number $U_\infty / \nu$ is high enough (Chung et al. Reference Chung, Hutchins, Schultz and Flack2021). We note that this condition is not well defined in highly compressible flows, where the free stream can be modulated by a succession of shock/expansion waves, which are expected to be promoted by higher Mach numbers and cold walls (which lower the speed of sound in the roughness sublayer, thus increasing the local Mach number). A general schematic of the dynamics emerging from the present database is offered in figure 10, where the evolution of both $\delta _{99}$ and $\delta _{\textit{IBL}}$ is sketched, as well as the presence of the initial shock wave at the smooth-to-rough transition. In the next sections, we consider streamwise locations at around $(x-x_{\textit{ref}})/\delta _{99,\textit{ref}}=40$ , or $x/\delta _{\textit{in}}=127$ , in order to investigate turbulence statistics in a region where the flow is in equilibrium with the new surface.

Figure 10. Schematic of the main features of a supersonic boundary layer developing over a positively skewed rough surface.

5. Turbulence and thermal statistics in the rough region

In this section, we discuss wall-normal turbulence statistics on a fixed streamwise location selected at $x/ \delta _{\textit{in}}=127$ . This particular value is chosen to be far enough from the smooth-to-rough transition so that the flow has a good degree of equilibrium with the rough surface (see § 4), and also avoids possible interference effects from the rough-to-smooth transition near the outflow. The relevant turbulence and roughness parameters for this station are reported in table 3.

First, we are interested in characterising the mean velocity profiles of rough cases, comparing them with reference smooth profiles at the same friction Reynolds number $\textit{Re}_{\tau }$ to assess outer-layer similarity (Townsend Reference Townsend1980) and compare the results with the incompressible flow case. To this end, the mean velocity defect is one of the most important quantities used to assess the similarity of friction-scaled turbulent relative motions (Chung et al. Reference Chung, Hutchins, Schultz and Flack2021), but its evaluation on incompressible turbulent boundary layers always presented challenges in precisely determining the boundary layer edge, where the external velocity $u_e$ is computed. The general outer-layer universality assumption reads

(5.1) \begin{equation} \frac {u_{I,e}-u_I(y)}{u_0} = F\left ( \frac {y}{\delta _0}\right)\!, \end{equation}

where $u_I(y)$ is the incompressible mean velocity profile while $\delta _0$ and $u_0$ represent the outer length and velocity scales, respectively, which are classically taken as $u_0=u_\tau$ , $\delta _0=1.6 \delta _{95}$ .

However, this scaling lacks generality even across incompressible numerical and experimental results at different Reynolds numbers, and improved versions have been proposed to avoid ambiguities related to the imprecise evaluation of the boundary layer thickness. Recently, Pirozzoli & Smits (Reference Pirozzoli and Smits2023) proposed a novel scaling considering the velocity and spatial scales $u_N$ and $\delta _N$ , defined as

(5.2) \begin{equation} \begin{aligned} u_N = \frac {H-1}{H}u_e, \\ \delta _N = \frac {H}{H-1} \delta ^{*},\end{aligned} \end{equation}

where $H$ and $\delta ^*$ are the shape factor and the displacement thickness defined for the incompressible regime (Pope Reference Pope2000). Pirozzoli & Smits (Reference Pirozzoli and Smits2023) select $u_0=u_N$ , $\delta _0=2 \delta _{N}$ as the outer velocity and length scales.

To our knowledge, this improved scaling has never been evaluated and compared with the classical version in compressible flows or in the presence of roughness. Although this paradigm is strictly valid only for incompressible wall-bounded flows, the velocity profile can be scaled by means of the Van Driest (Reference Van Driest1951) transformation, which accounts for mean density variations such that $u_I =u_{\textit{VD}} = \int _0^{u} \sqrt {\bar {\rho }/\bar {\rho }_w} {\rm d}u$ .

Using the present dataset, along with DNS databases of Sillero, Jiménez & Moser (Reference Sillero, Jiménez and Moser2013) and Cogo et al. (Reference Cogo, Modesti, Bernardini and Picano2025), we can assess outer-layer similarity among incompressible and compressible smooth-wall profiles, as well as compressible rough-wall cases, using both the classical and Pirozzoli & Smits (Reference Pirozzoli and Smits2023) scalings, as reported in figure 11. Looking at figure 11(a), we can see that the classical scaling does not effectively collapse all profiles, with case CB_A showing the highest mismatch compared with the others. Figure 11(b) reports the outer-layer similarity obtained with the scaling of Pirozzoli & Smits (Reference Pirozzoli and Smits2023), which is able to collapse all cases extremely well, and agrees with the proposed compound logarithmic/parabolic fit without additional tuning of parameters. Here, we remark that the use of the Van Driest (Reference Van Driest1951) compressibility transformation is fundamental since the scaling proposed by Pirozzoli & Smits (Reference Pirozzoli and Smits2023) is based on the incompressible definition of integral length scales. These results demonstrate that the scaling of Pirozzoli & Smits (Reference Pirozzoli and Smits2023) can be generalised to both compressible flows and rough walls for all considered cases, thus confirming the presence of outer-layer similarity across different rough-wall cases and their related smooth-wall reference under the premises of this scaling. Given that this is the first time that this scaling has been used for compressible or rough flows, future works are needed to explore the range of its applicability.

Figure 11. Velocity defect profiles for the classical (a), and Pirozzoli & Smits (Reference Pirozzoli and Smits2023) (b) scalings. All velocity profiles are scaled using the transformation of Van Driest (Reference Van Driest1951). Incompressible and supersonic (SM_M2) smooth-wall data at $\textit{Re}_{\tau }=1571$ and $\textit{Re}_{\tau }=1528$ , respectively, are taken from Sillero et al. (Reference Sillero, Jiménez and Moser2013) and Cogo et al. (Reference Cogo, Modesti, Bernardini and Picano2025). The Hama fit and its parameters are given by the compound logarithmic/parabolic fit described in Pirozzoli & Smits (Reference Pirozzoli and Smits2023).

We now turn our attention to the inner layer, evaluating mean velocity profiles for smooth and rough walls in absolute and relative terms. Figure 12 shows mean velocity profiles for smooth and rough configurations in the untransformed and transformed cases, visible in figures 12(a) and 12(b), respectively. We note that we still employ the Van Driest (Reference Van Driest1951) scaling in order to be consistent with the velocity defect profiles computed in figure 11, for which more recent transformations have not been validated. For each case, we compute the resulting velocity deficit function, shown in figures 12(c) and 12(d), as

(5.3) \begin{equation} \Delta u^+ = \tilde {u}^+_S (y^+)-\tilde {u}^+_R (y^+), \end{equation}

where $\tilde {u}^+_S$ and $\tilde {u}^+_R$ are the smooth and rough profiles, respectively. The smooth-wall profiles are generated at exactly matched $\textit{Re}_\tau$ using the model of Hasan et al. (Reference Hasan, Larsson, Pirozzoli and Pecnik2024) in order to reduce errors due to wake effects, and have been validated against the smooth-wall reference of Cogo et al. (Reference Cogo, Modesti, Bernardini and Picano2025). The virtual origin $d$ is calculated iteratively to minimise the standard deviation of $\Delta u^+_{\textit{VD}}(y^+)$ in the range $100\lt y^*\lt 0.2Re_\tau$ . The nominal velocity deficit $\Delta U^+$ is calculated in the same range as the averaged value and reported in table 4 alongside the virtual origin shift.

Table 4. Boundary layer properties considering the reference smooth-wall cases. Here, $d/k$ is the virtual origin wall-normal location relative to $k$ and $\Delta U^+$ is the velocity deficit, which is used to compute the equivalent sand-grain roughness height by inverting the relation $\Delta U^{+}=1 / \kappa \ln k_s^{+}+A-B_s$ with $A=5.2$ , $\kappa =0.41$ and $B_s=8.5$ . The ratio $k_s/k$ is obtained as $k_s^+ /k^+$

Figure 12. Mean velocity profiles for smooth and rough-wall cases obtained at stations listed in table 3. Panels (b–d) show the profiles scaled with the velocity transformation of Van Driest (Reference Van Driest1951). All rough cases are shifted by $d$ according to table 4.

Flow case CB_R shows the most pronounced velocity shift, followed by CB_S, and with CB_A and DM_S attaining more similar values, confirming the trend of the skin friction coefficient.

The value of $\Delta U^+$ is then related to the equivalent sand-grain roughness height $k_s$ by means of the equation

(5.4) \begin{equation} \Delta U^+ = 1/ \kappa \ln {k_s^+} + A - B_s, \end{equation}

where $\kappa =0.41$ is the von Kármán constant, $A=5.2$ and $B_s=8.5$ (Chung et al. Reference Chung, Hutchins, Schultz and Flack2021). The computed values of $k_s^+$ are reported in table 4, together with the ratio $k_s/k$ , which gives an idea of the equivalent sand-grain roughness that can generate a similar incompressible velocity deficit.

Another indicator of the presence of outer-layer similarity for turbulent statistics can be found by looking at turbulent Reynolds stresses, reported in figure 13 as a function of the wall-normal coordinate, both in inner, figure 13(a), and outer, figure 13(b), units. In the incompressible theory, Reynolds stresses of rough-wall cases are expected to recover a smooth-wall-like behaviour in the outer layer. In the present database, the comparison is made with a single reference smooth profile obtained from Cogo et al. (Reference Cogo, Modesti, Bernardini and Picano2025) at $M_\infty =2$ and $\textit{Re}_\tau =1528$ , hence slight discrepancies are expected for rough-wall cases that differ from this value, see table 3. Nevertheless, we observe a general good agreement between rough-wall cases and the smooth reference, an additional indication that outer-layer similarity is present for the velocity field. We only note a minor mismatch for the streamwise component $\tau _{11}^+$ , as observed by Cogo et al. (Reference Cogo, Modesti, Bernardini and Picano2025), who attributed this effect to the modulation of turbulent structures on the outer layer induced by the initial shock wave. The general agreement between smooth and rough cases of the turbulent velocity fluctuation profiles is in agreement with the computational study of Yu et al. (Reference Yu, Liu, Tang, Yuan and Xu2023), who emulated the effect of roughness for supersonic boundary layers through a synthetic velocity boundary condition. However, different experiments (Latin & Bowersox Reference Latin and Bowersox2000; Ekoto et al. Reference Ekoto, Bowersox, Beutner and Goss2008; Kocher et al. Reference Kocher, Kreth, Schmisseur, LaLonde and Combs2022) reported clear deviations from the smooth-wall Reynolds stresses attributed to the presence of shock waves propagating from the roughness crests to the outer layer. In the present dataset the roughness crest is subsonic, thus the validity of outer-layer similarity at higher Mach number remains an open question.

Figure 13. Turbulent velocity fluctuations $\tau _{\textit{ij}}=\widetilde {u_i^{\prime }u_j^{\prime }}$ scaled with the wall shear stress $\tau _w$ as a function of the wall-normal distance in wall units (a) $y^+, y^+-d^+$ and outer units (b) $y/\delta _{99},(y-d)/\delta _{99}$ . Rough-wall cases are adjusted using a virtual origin shift according to table 4. The smooth-wall reference SM_M2 is from Cogo et al. (Reference Cogo, Modesti, Bernardini and Picano2025).

Figure 14. (a) Mean temperature profiles as a function of the wall-normal distance $y^+$ . (b) Temperature fluctuations scaled with the wall temperature as a function of $y^+$ . The smooth-wall reference SM_M2 is from Cogo et al. (Reference Cogo, Modesti, Bernardini and Picano2025).

Figure 14 shows the mean and root-mean-square temperature profiles, respectively in figures 14(a) and 14(b). Even for these statistics, the rough-wall cases are compared with the smooth reference from Cogo et al. (Reference Cogo, Modesti, Bernardini and Picano2025). Differences between smooth and rough flow cases are apparent both for averaged and Root mean square temperature profiles, and they seem enhanced for the cases with highest added drag. It follows that outer-layer similarity between rough and smooth cases is not observed for the mean temperature nor for the fluctuations. This fact was observed by previous works, see Modesti et al. (Reference Modesti, Sathyanarayana, Salvadore and Bernardini2022) and Cogo et al. (Reference Cogo, Modesti, Bernardini and Picano2025), however, it does not imply that the relative coupling between the kinetic and thermal fields is inherently different across smooth and rough cases, especially in the outer layer (Modesti et al. Reference Modesti, Sathyanarayana, Salvadore and Bernardini2022). For the temperature fluctuations, the absence of outer-layer similarity is attributed to the roughness disrupting the near-wall peak of the temperature fluctuations, leading to a nearly constant temperature layer below the roughness crest. The mean temperature gradient is a driving factor in the generation of the temperature fluctuations, and this is enhanced in the outer layer, leading to the formation of temperature fluctuations away from the wall (Cogo et al. Reference Cogo, Modesti, Bernardini and Picano2025).

Following the above discussion, we now turn our attention to the relation between the mean and fluctuating temperature and velocity fields, which is classically characterised under the framework of the Reynolds analogy. The fundamental relationship between the mean total enthalpy and velocity fields, derived for compressible boundary layers over smooth walls, can be written as (Walz Reference Walz1969)

(5.5) \begin{equation} \bar {H}_r-\bar {H}_w = U_w \bar {u}, \end{equation}

where $\bar {H}_r=c_pT+r u^2/(2)$ is the recovery total enthalpy, $r\approx Pr^{1/3}$ is the recovery factor and $U_w$ is a constant velocity scale, which under adiabatic conditions is equal to zero. We note that more recent formulations of this relation, such as the one of Zhang et al. (Reference Zhang, Bi, Hussain and She2014), fall back on 5.5 for adiabatic walls, which characterises the present database. Figure 15(a) reports the term $(\bar {H}_r-\bar {H}_w)$ for the present database and compares it with the smooth-wall reference of Cogo et al. (Reference Cogo, Modesti, Bernardini and Picano2025), which are observed to agree reasonably well for $y/\delta _{99}\gt 0.5$ . Even though below this point the rough-wall profiles deviate from the smooth-wall reference, they are still close to zero, which is the value expected for smooth adiabatic walls. This is reflected in the resulting mean temperature–velocity relation, reported in figure 15(b), which follows the same trend. From similar observations, Modesti et al. (Reference Modesti, Sathyanarayana, Salvadore and Bernardini2022) linked the absence of outer-layer similarity for the mean temperature field to the existence of an approximate quadratic relationship between the temperature and velocity even in the presence of roughness, leading to a temperature difference between rough and smooth cases which is a function of $y^+$ .

Figure 15. (a) Left-hand side of (5.5) as a function of the wall-normal coordinate $y/\delta _{99}$ . (b) Mean temperature profiles as a function of the mean velocity.

For the fluctuating temperature–velocity relationship, we consider first the strong Reynolds analogy (SRA) relation formulated for adiabatic flows (Gaviglio Reference Gaviglio1987; Huang, Coleman & Bradshaw Reference Huang, Coleman and Bradshaw1995), which reads

(5.6) \begin{equation} \frac {\big(\widetilde {T^{\prime \prime 2}}\big)^{1 / 2} / \tilde {T}}{(\gamma -1) \tilde {M}^2\big(\widetilde {u^{\prime \prime 2}}\big)^{1 / 2} / \tilde {u}} \approx 1. \end{equation}

The left-hand side of (5.6) is reported in figure 16 for all cases in the present database, along with the smooth-wall reference. Here, we observe an excellent agreement between all rough cases and the smooth reference, with deviations only visible very close to the roughness elements, $y/\delta _{99}\lt 0.1$ . Again, we note that, similarly to the discussion on the mean field, the absence of outer-layer similarity for temperature fluctuations does not directly imply that even their relation to the velocity fluctuating field does not agree with the reference smooth-wall case.

Figure 16. Classical formulation of the SRA, left-hand side of (5.6), as a function of the wall-normal coordinate $y/\delta _{99}$ .

6. Drag partition and shielding effects

In this section, we investigate the drag partitioning resulting from pressure and viscous contributions applied to roughness elements and to the underlying flat surface. From this investigation, we aim to provide a bigger picture of how the shape and arrangement of the roughness elements are able to modulate the dynamics of compressible boundary layers.

Although the analysis is facilitated by the fact that there is a clear separation between the bottom surface and the individual wall-mounted elements, there are mainly two aspects that can mislead our intuition for this analysis. The former is related to the fact that the wake of an individual roughness element on a rough surface is more complex than that of an obstacle in a non-turbulent free stream Raupach (Reference Raupach1992). The latter comes from the fact that the wake interaction between subsequent elements (in both the streamwise and spanwise directions) can create patterns that are inherently different from what we can predict just by knowing the dynamics of a single isolated wall-mounted element, due to the mutual sheltering between roughness elements.

This is visible in figure 17, where the mean streamwise velocity in different planes is used to highlight the 3-D wake. Volumetric sheltering is most evident for the CB_A case, given the aligned pattern, which gives rise to high-speed regions of relatively undisturbed flow. A much more complex interaction pattern is observed for the CB_S and CB_R cases, showing that subsequent wakes are oriented in both the streamwise and spanwise directions, following the staggered arrangement. However, this is not visible for the DM_S case, for which wakes are present in much more confined regions and volumetric sheltering does not occur.

Figure 17. Averaged velocity field $\bar {u}/u_{\infty }$ over $x$ - $y$ and $x$ - $z$ slices for (a) CB_A, (b) CB_S, (c) CB_R, (d) DM_S. The $x$ - $z$ slices are selected at $y/k=0.5$ and $y/k=0.99$ , respectively above and below the horizontal solid red line. The horizontal dashed blue line indicates the $z$ location where $x$ - $y$ slices are selected. Black contours indicate locations indicates locations in which $\bar {u} =0$ .

In order to investigate how the element sheltering influences the total drag, we focus on a unit roughness patch of area $A_t$ (see figure 2), on which the flow exerts a force in the streamwise direction $F_x$ . The total stress at the wall $\tau _{\textit{tot}}= F_x /A_t$ is split into the one acting on the flat bottom wall $\tau _S = F_{x_S}/A_t$ , and the one acting on each roughness element $\tau _R = F_{x_R}/A_t$ (Raupach Reference Raupach1992; Shao & Yang Reference Shao and Yang2008). Additionally, $\tau _R$ is further split into the pressure $\tau _{R,p}$ and the viscous contribution $\tau _{R,v}$ (see the sketch of figure 18 b)

(6.1) \begin{equation} \tau _{\textit{tot}} = \frac {F_x}{A_t}= \underbrace {\frac {F_{x_S}}{A_t}}_{\tau _S}+\underbrace {\frac {F_{x_{R,p}}}{A_t}}_{\tau _{R,p}}+\underbrace {\frac {F_{x_{R,v}}}{A_t}}_{\tau _{R,v}}. \end{equation}

Figure 18. Total stress partition on a unit roughness patch (located in the black box displayed in figure 19). Here, $\tau _S$ is the friction at the bottom surface, $\tau _R$ is the drag of roughness elements, split into pressure and shear stress contributions $\tau _R= \tau _{R,p}+ \tau _{R,v}$ . All stresses have been rescaled with the dynamic pressure $q_{\infty }=1/2 \rho _{\infty } u_{\infty }^2$ in order to be comparable to $C_{\!f}=\tau /q_{\infty }$ . Green bars corresponding to the $\tau _{R,v}$ contribution display a darker region, representative of the viscous drag on the element due only to the top surface (parallel to the bottom surface).

Figure 18 shows the individual stress contributions along with the total. We remark that these values are calculated on a single roughness repeating unit, although they are in close agreement with the spanwise-averaged values reported in figure 7. In fact, we observe that $\tau _{\textit{tot}}$ is similar for aligned cubes and diamonds (slightly higher for the latter case), while staggered and rotated cubes have progressively more drag. Staggered, CB_S, and rotated, CB_R, cubes have very similar contributions from pressure drag, while this is lower for aligned cubes, CB_A, a sign that mutual sheltering effects play a more important role for the latter case. It is also interesting to note that $\tau _S$ is also similar between the CB_S and CB_R cases, while the viscous drag applied to each element $\tau _{R,v}$ is higher for the rotated case in both absolute terms and percentage (relative to the corresponding $\tau _{\textit{tot}}$ ).

From this fact, we can preliminarily state that the main factor contributing to the added drag of rotated cubes (compared with staggered ones) is the viscous stress applied to the roughness elements $\tau _{R,v}$ . On this point, we can further specify that the shear stress is higher on the lateral sides of the cube, whereas it has comparable intensity on the top side (dark green region).

We note that the fraction of $\tau _{R,v}$ applied on the top side of each roughness element (dark green region) is comparable between the three cases featuring cubical elements, which can be explained by the fact that they all have the same area. Following the same reasoning, diamond-shaped elements are also expected to have a double amount of shear stress on their top side, given that their area is greater by a factor of two ( $2 k^2$ ).

It is interesting to note how diamond-shaped elements, DM_S, which are characterised by the lowest pressure contribution of all cases in absolute and relative terms, are able to produce a total stress comparable $\tau _{\textit{tot}}$ to the one of aligned cubes, CB_A. This is mainly due to a considerable increase in the wall shear stress applied to the sides of the element $\tau _{R,v}$ , which takes up to $34.2\,\%$ of the total stress, which can be expected from a larger total surface exposed to the flow, but also a considerable amount of shear stress applied to the bottom surface $\tau _S$ ( $14.3\,\%$ ). The latter figure also has to be supported by the fact that the diamond-shaped elements are characterised by a plan solidity that is twice that of all other cases, leaving less area available to the bottom surface.

Although it is important to note that the computed values of $C_{\!f}$ , and their relative contributions, strictly refer to the Reynolds number under scrutiny, we remind the reader that all cases feature an identical average state of the incoming smooth-wall turbulent boundary layer, thus relative contributions are directly related to differences in roughness geometry.

In figure 19 we visualise the contours of $\tau _S$ , where negative values highlight the wake region behind each element. We can infer the relative importance of the separated regions by comparing the areas subjected to negative wall shear stress. It is interesting to note how each roughness geometry establishes a unique interaction pattern between adjacent elements, which is known to highly impact the resulting drag prediction (Yang et al. Reference Yang, Sadique, Mittal and Meneveau2016; Meneveau et al. Reference Meneveau, Hutchins and Chung2024).

Aligned cubes show preferential gaps in which the flow has a positive wall shear stress but displays a very intense wake denoted by strongly negative values of $\tau _S$ . Staggered cubes show a more complex pattern, since each cube’s wake interacts with two successive cubes. This creates a large region of $\tau _S$ close to zero, with local spots of negative values (in front of each cube), and positive values (spanwise gaps between each cube). Rotated cubes show a similar wake pattern as in the staggered case, this time with higher positive values (spanwise gaps) and fewer negative ones (front side of each cube). Lastly, diamond-shaped elements exhibit very small recirculation areas which, compared with all other cases, do not seem to interact directly with subsequent elements.

Table 5. Averaged shear stress on the bottom surface $\tau _S/q_\infty$ for the entire roughness patch and conditioned to positive ( $\tau _S\gt 0$ ) and negative ( $\tau _S\lt 0$ ) values. Data collected on a unit surface patch located at $x/\delta _{\textit{in}}=127$ . All stresses have been rescaled with the dynamic pressure $q_{\infty }=1/2 \rho _{\infty } u_{\infty }^2$ in order to be comparable to $C_{\!f}=\tau /q_{\infty }$

Figure 19. Averaged shear stress ( $x$ component) $\tau _S/q_{\infty }$ on the bottom surface for (a) CB_A, (b) CB_S, (c) CB_R, (d) DM_S. Black rectangles indicate a sample unit repeating patch.

The positive and negative contributions of $\tau _S$ represent a key indicator of the wake interaction pattern, and are reported with conditional averages in table 5. Here, we can see that all roughness geometries that feature cubical elements have a strong contribution of negative values of $\tau _S$ , which amounts to almost $68\, \%$ of the mean value of $\tau _S$ for the staggered case (the highest). On the other hand, diamond-shaped elements do not display such negative values and retain a nearly zero value of $\tau _S$ downstream of each element, thus high recirculation regions are mostly avoided. We argue that this feature is promoted by both the slender shapes of diamond elements and by the fact that their topology exhibits the highest planar solidity of all cases, thus reducing the physical space between wall-mounted elements. For these reasons, in the DM_S case the flow is more prone to ’skim over’ the elements and is less forced around gaps and troughs, thus reducing the added drag.