1. Introduction
The research on supersonic and hypersonic flight has gained new interest in recent years upon the development of different flight systems, such as high-speed commercial aircraft and space systems on their way to and from orbit. The development of a new class of high-speed vehicles is driven by recent advances in propulsion, manufacturing techniques and software, which enable engineers to tackle long-standing problems from a different perspective (Leyva Reference Leyva2017). However, there are still many critical aspects to consider in order to design flight systems that can withstand huge mechanical and thermal loads in a wide range of different flight regimes (Candler Reference Candler2019), along with more advanced propulsion systems that can be operated efficiently throughout the flight path (Urzay Reference Urzay2018).
One key aspect considered in this study is the detrimental effect of distributed surface imperfections and roughness, which may be a subsidiary of the presence of a thermal protection system (Uyanna & Najafi Reference Uyanna and Najafi2020), or may arise from unforeseen events such as atmospheric hazards, such as ice particles impacting the surface (Habeck et al. Reference Habeck, Kroells, Schwartzentruber and Candler2024). These factors permanently modify the wetted area of the flight system for the remainder of the trajectory, abruptly affecting the dynamics of the compressible boundary layer in a way that is not well known.
Multiple physical aspects of wall-bounded flows over rough walls have historically been studied for low-speed incompressible flows (Nikuradse Reference Nikuradse1933), with the primary goal of correlating the geometric properties of the surface with its hydrodynamic response, which has been standardised using the equivalent sand-grain roughness height
$k_s$
(Moody Reference Moody1944). However, relating the geometrical and hydrodynamic properties of a rough surface is still an open problem (Chung et al. Reference Chung, Hutchins, Schultz and Flack2021).
In this scenario, the parametrisation of different types of roughness has been of critical importance, with the aim of identifying the most influential characteristics that can be related to the equivalent sand-grain roughness height
$k_s$
, or the virtual origin
$d$
, which is a parameter that indicates a wall offset determined by the flow (Chung et al. Reference Chung, Hutchins, Schultz and Flack2021).
Among the different surface properties that have been found to have the greatest impact on the resulting drag, there is the characteristic roughness height
$k$
of individual roughness elements, the frontal solidity
$\lambda _{\!f}$
, defined as the ratio of the projected frontal area to the total plan area, and the plan solidity
$\lambda _p$
, which is the ratio of the projected horizontal area to the total plan area. Other factors, such as skewness
$Sk$
, are often considered in empirical predictive correlations (Macdonald, Griffiths & Hall Reference Macdonald, Griffiths and Hall1998; Flack, Schultz & Barros Reference Flack, Schultz and Barros2020), which, however, struggle to perform outside the specific range of considered topologies.
From a fluid dynamics perspective, the possibility of reducing the complex parameter space and estimating the added drag of a rough surface compared with a smoother one just from geometrical constraints is viable only when the roughness affects a small region in the vicinity of the wall, i.e. the separation between turbulent and roughness scales is enough to consider the flow to be in the fully rough regime, hence Reynolds number independent (Jiménez Reference Jiménez2004). This means that characteristic roughness lengths such as
$k$
have to be related to the boundary layer thickness
$\delta$
and the viscous length scale
$\delta _\nu =\nu _w/u_\tau$
, where
$u_\tau =\sqrt {\tau _w/\rho _w}$
is the friction velocity,
$\tau_w$
is the wall shear stress,
$\nu _w$
and
$\rho _w$
are the kinematic viscosity and density at the wall, respectively. From these quantities, important parameters such as the friction Reynolds number
$\textit{Re}_\tau =\delta / \delta _\nu$
and the roughness Reynolds number
$k^+=k/\delta _\nu$
can be defined, and thresholds on the expected applicability range of the fully rough regime can be formulated using these or similar control variables, such as the ratio
$k/\delta$
or
$k_s^+ = k_s /\delta _\nu$
.
Under specific circumstances, a turbulent boundary layer over a rough wall exhibits a similarity in the velocity profile in the outer layer to its smooth counterpart at the same
$\textit{Re}_\tau$
(Chung et al. Reference Chung, Hutchins, Schultz and Flack2021). However, the assessment of outer-layer similarity is an open issue even for smooth-wall flows at different Reynolds numbers, especially for boundary layers, since the correct definition of an outer length scale is still debated (Pirozzoli & Smits Reference Pirozzoli and Smits2023). One key indication for rough walls is the typical logarithmic profile in the inner layer, although with a different intercept from the smooth wall, evaluated through the roughness function
$\Delta U^+$
(Chung et al. Reference Chung, Hutchins, Schultz and Flack2021).
However, only by having access to the flow characteristics in the roughness sub-layer can we explain the specific dynamics that generated the observed added drag, which is often difficult to grasp intuitively due to mutual volumetric sheltering between different regions of the surface. This phenomenon can be more easily observed in prism-shaped roughness, where the wake of a single wall-mounted element can interact with its neighbours depending on their vicinity, so that the flow can be forced through individual gaps in some cases, and it can ’skim over’ roughness elements in others (Raupach Reference Raupach1992). The investigation of these effects, which requires the use of scale-resolving simulations, has been extremely valuable in connecting the small-scale dynamics of the roughness sub-layer with the resulting added drag, and it has been embodied in more recent reduced-order models (Yang et al. Reference Yang, Sadique, Mittal and Meneveau2016; Meneveau, Hutchins & Chung Reference Meneveau, Hutchins and Chung2024), which reported considerable improvements from previous empirical correlations and can be applied in a great variety of topologies. The growing availability of computational resources in recent years is expected to foster the systematic analysis of different roughness topologies with high-fidelity simulations, particularly in flow configurations that are not easily reproducible by experimental campaigns, such as supersonic and hypersonic flows. In this context, the study of prism-shaped roughness is not only valuable to inform models by enabling a systematic assessment of different arrangements of roughness elements of different shapes, it is also relevant for high-speed flows for being reminiscent of the behaviour of ablative materials, which develop a network of grooves and ridges on the surface that create a distinct cross-hatch pattern (Probstein & Gold Reference Probstein and Gold1970).
This fact motivated recent experimental campaigns investigating supersonic and hypersonic boundary layers over rough surfaces, consisting of three-dimensional (3-D) cubes, 2-D bars, diamond-shaped elements and realistic topologies (Ekoto et al. Reference Ekoto, Bowersox, Beutner and Goss2008; Peltier, Humble & Bowersox Reference Peltier, Humble and Bowersox2016; Williams et al. Reference Williams, Sahoo, Papageorge and Smits2021; Kocher et al. Reference Kocher, Kreth, Schmisseur, LaLonde and Combs2022), which provided the first quantitative assessments of roughness in the compressible regime. In many cases, compressibility effects were shown to be so strong they generated a pattern of shock waves emanating from individual elements, which traversed the boundary layer and emerged into the free stream. However, it is still not clear if this behaviour is due to specific geometries, such as diamond-shaped elements like the ones of Peltier et al. (Reference Peltier, Humble and Bowersox2016), or particular flow conditions such as very high friction Reynolds numbers which are still only achievable in wind tunnels. To further complicate the picture, a general consensus is lacking among these studies around the validity of compressibility transformations as a tool to recover outer-layer similarity of a rough-wall case, similarly to the recovery of the log law in smooth counterparts. Additionally, very few of these studies addressed the spatial evolution of the turbulent boundary layer transitioning from a smooth to a rough surface, thus neglecting the development of the internal boundary layer (IBL), which is a key feature of these type of flows that directly relates to the degree of equilibrium that the boundary layer has reached with the new surface (Elliott Reference Elliott1958; Rouhi, Chung & Hutchins Reference Rouhi, Chung and Hutchins2019). These points were addressed in one of the few computational studies available so far (Cogo et al. Reference Cogo, Modesti, Bernardini and Picano2025), which used Direct Numerical Simulations (DNSs) to compare subsonic and supersonic turbulent boundary layers over 3-D aligned cubes. Although limited to one specific geometry, this work provided a systematic comparison of the smooth-to-rough transition in subsonic and subsonic regimes, offering a novel definition for tracking the IBL thickness that works well for both. Additionally, turbulence and thermal statistics were evaluated in the rough region, showing that compressibility transformations help the achievement of an outer-layer similarity for the velocity-related statistics of the supersonic case, which is, however, clearly not present for the temperature field, as previously noted by Modesti et al. (Reference Modesti, Sathyanarayana, Salvadore and Bernardini2022) for compressible turbulent channels. Similar conclusions were also reached in the study of Yu et al. (Reference Yu, Liu, Tang, Yuan and Xu2023), which simulated the effects of wall disturbances on supersonic boundary layers by artificially increasing the instantaneous velocity field in a sinusoidal manner below a certain wall-normal location, thus emulating the presence of a rough wall.
These and other authors called for further assessments of compressible turbulent boundary layers over rough surfaces, noting a severe lack of studies systematically investigating the wide parameter space that includes different topologies, Mach numbers and wall temperature conditions while retaining the possibility to investigate in detail both large and small scales.
In the present work, we propose a step forward in this direction by performing DNS of zero-pressure-gradient turbulent boundary layers at free-stream Mach number
$M_\infty =2$
and friction Reynolds numbers
$\textit{Re}_\tau$
up to 2200 developing from a smooth surface to a rough one. Expanding on the work of Cogo et al. (Reference Cogo, Modesti, Bernardini and Picano2025), the current dataset includes four simulations that incorporate cubical- and diamond-shaped elements arranged in both aligned and staggered configurations. One of the cases involves cubes rotated relative to the flow direction. We conduct a general assessment of the instantaneous and mean fields, followed by a detailed analysis of the smooth-to-rough transition and internal boundary layer development. Then, turbulence and thermal statistics are investigated towards the end of the rough region, where equilibrium with the new surface is expected. Additionally, a detailed analysis of drag partition and volume sheltering is carried out in the roughness sub-layer, providing valuable insights for the development of reduced-order models.
2. Methodology and roughness properties
We solve the compressible Navier–Stokes equations for a viscous, heat-conducting gas
\begin{equation} \begin{aligned} \frac {\partial \rho }{\partial t}+\frac {\partial (\rho u_j)}{\partial x_j}=0,\\ \frac {\partial (\rho u_i)}{\partial t}+\frac {\partial (\rho u_i u_j)}{\partial x_j}+\frac {\partial p}{\partial x_i}-\frac {\partial \sigma _{\textit{ij}}}{\partial x_j}=0,\\ \frac {\partial (\rho E)}{\partial t}+\frac {\partial (\rho E u_j+pu_j)}{\partial x_j}-\frac {\partial (\sigma _{\textit{ij}}u_i-q_j)}{\partial x_j}=0, \end{aligned} \end{equation}
where
$\rho$
is the density,
$u_i$
denotes the velocity component in the ith Cartesian direction (
$i=1,2,3$
),
$p$
is the thermodynamic pressure,
$E=c_v T+u_i u_i/2$
the total energy per unit mass and
\begin{equation} \sigma _{\textit{ij}}=\mu \left ( \frac {\partial u_i}{\partial x_j}+\frac {\partial u_j}{\partial x_i}-\frac {2}{3}\frac {\partial u_k}{\partial x_k}\delta _{\textit{ij}}\right)\!, \qquad q_j=-k\frac {\partial T}{\partial x_j}, \end{equation}
are the viscous stress tensor and the heat flux vector, respectively. The molecular viscosity
$\mu$
is assumed to follow Sutherland’s law, with a reference free-stream temperature
$T_{\infty }=220.0$
K. The thermal conductivity
$k$
is related to the viscosity through the Prandtl number
$\textit{Pr}=0.72$
,
$k=c_p \mu /\textit{Pr}$
, where
$c_p$
is the specific heat at constant pressure. This model is complemented by the equation of state for a calorically perfect gas. The system of equations is solved on a Cartesian grid using the in-house code STREAmS (Bernardini et al. Reference Bernardini, Modesti, Salvadore and Pirozzoli2021, Reference Bernardini, Modesti, Salvadore, Sathyanarayana, Della Posta and Pirozzoli2023), which has been extensively validated in numerous canonical flow configurations (Bernardini, Pirozzoli & Grasso Reference Bernardini, Pirozzoli and Grasso2011; Bernardini & Pirozzoli Reference Bernardini and Pirozzoli2011; Modesti & Pirozzoli Reference Modesti and Pirozzoli2016; Cogo et al. Reference Cogo, Salvadore, Picano and Bernardini2022, Reference Cogo, Baù, Chinappi, Bernardini and Picano2023, Reference Cogo, Modesti, Bernardini and Picano2025). Convective terms are discretised using sixth-order, energy-preserving schemes applied in shock-free regions, while a high-order shock capturing scheme Weighted-Essentially Non-Oscillatory is applied when shock waves are identified by the Ducros sensor (Ducros et al. Reference Ducros, Ferrand, Nicoud, Weber, Darracq, Gacherieu and Poinsot1999). Viscous terms are discretised using a locally conservative formulation (De Vanna et al. Reference De Vanna, Benato, Picano and Benini2021) with second-order accuracy.
Schematic of the computational set-up for a turbulent boundary layer flow over different roughness topologies.
The complexity of the geometry is handled using a ghost-point-forcing immersed boundary method (Piquet, Roussel & Hadjadj Reference Piquet, Roussel and Hadjadj2016; De Vanna et al. Reference De Vanna, Picano and Benini2020), which is used to enforce no-slip adiabatic boundary conditions on the solid wall. The reader can find more details in previous applications and validation of this methodology in Modesti et al. (Reference Modesti, Sathyanarayana, Salvadore and Bernardini2022) and Cogo et al. (Reference Cogo, Modesti, Bernardini and Picano2025). In the present database, four different roughness topologies have been investigated with separate simulations sharing the same computational set-up, explained as follows. The rough-wall simulations are divided into three sections, visible in figure 1. The first is a smooth-wall region designed to allow the development of a turbulent boundary layer, generated using a recycling/rescaling method with the recycling plane located at
$ x = 40\, \delta _{\textit{in}}$
. The second region extends from
$ x = 55\, \delta _{\textit{in}}$
to
$ x = 147\, \delta _{\textit{in}}$
and features prism-shaped elements of height
$ k$
, differing in arrangement, orientation and shape. Lastly, a short smooth-wall segment is included at the downstream end of the domain, extending to
$ x = 150\, \delta _{\textit{in}}$
, serving as a buffer zone before the outflow. All cases share a computational mesh composed of
$20240\times 556\times 1408$
nodes in the
$x$
,
$y$
and
$z$
directions, respectively, and the mesh is adequately stretched to accurately resolve the roughness substrate, assuring a resolution for individual elements of at least 20 nodes in wall-parallel directions, and 40 nodes in the wall-normal direction, following the validation campaign conducted in smaller set-ups and from the work of Modesti et al. (Reference Modesti, Sathyanarayana, Salvadore and Bernardini2022).
In the rough region, we consider a combination of different prism-shaped roughness elements, focusing on highlighting differences in the arrangement, orientation and shape of roughness elements, see figure 2. We consider the same roughness geometry as that found in Modesti et al. (Reference Modesti, Sathyanarayana, Salvadore and Bernardini2022) and Cogo et al. (Reference Cogo, Modesti, Bernardini and Picano2025), namely aligned 3-D cubes of size
$k$
and separated by
$\lambda _x=\lambda _z=2k$
in both streamwise and spanwise directions, which is case CB_A shown in figure 2(a). In the second arrangement, CB_S, we consider staggered elements where every other row is shifted in the spanwise direction by
$\lambda _z/2$
, see figure 2(b). The third arrangement (CB_R) features staggered elements rotated by 45
$^\circ$
around the wall-normal axis passing through the elements’ centre, see figure 2(c). Lastly, for case DM_S we consider diamonds-shaped elements, where the streamwise axis is twice the spanwise one, see figure 2(d). We remark that the latter topology has been studied in several previous studies in the supersonic (Ekoto et al. Reference Ekoto, Bowersox, Beutner and Goss2008; Kocher et al. Reference Kocher, Kreth, Schmisseur, LaLonde and Combs2022) and hypersonic (Peltier et al. Reference Peltier, Humble and Bowersox2016) regimes, being reminiscent of the cross-hatching patterns forming on ablative surfaces.
Schematic of the different roughness patterns (flow from left to right). Cubes of size
$k$
are showed in different arrangements: aligned CB_A (a), staggered CB_S (b) and rotated of
$45^{\circ }$
CB_R (c). Panel (d) shows the diamond-shaped elements of DM_S, obtained by rescaling the CB_R elements by a factor of two in the streamwise direction (horizontal).
Even though realistic roughness patters found in high-speed vehicles are often considerably more complex than regular arrays of roughness elements, these particular shapes have been chosen to enable a systematic comparison of specific geometrical properties and their influence on the flow field, thus isolating cause–effect relationships that are fundamental to deriving engineering models. For example, the CB_A and CB_S cases share the same surface properties as those used in various predictive correlations (Chung et al. Reference Chung, Hutchins, Schultz and Flack2021), but are expected to generate different drag only due to their arrangement (aligned and staggered). At the same time, the CB_S, CB_R and DM_S cases share the same staggered arrangement, but differ in the individual element shape, thus affecting the dynamics of the flow in the roughness sublayer. These simplified patterns can then be related to more realistic configurations in future studies, incorporating additional effects, similarly to what is done in the more extensive literature concerning the incompressible flows regime (Volino, Schultz & Flack Reference Volino, Schultz and Flack2011; Yang et al. Reference Yang, Sadique, Mittal and Meneveau2016, Reference Yang, Xu, Huang and Ge2019). We also note that the diamond-shaped elements of the present database share the same half-angle of the leading edge as the study of Kocher et al. (Reference Kocher, Kreth, Schmisseur, LaLonde and Combs2022), namely
$\arctan (1/2)$
or
$26.6^\circ$
, which, however, differs in other parameters, such as the much higher friction Reynolds number
$\textit{Re}_\tau$
(approximately more than 5 times) that is not achievable with current computational resources.
Table 1 reports relevant roughness properties, which are often used by predictive correlations (Chung et al. Reference Chung, Hutchins, Schultz and Flack2021). These roughness parameters are computed on the minimum repetitive surface patch, indicated with a red rectangle in figure 2.
Roughness surface properties
Considering the height map
$h=f(x,z)$
representative of a given rough surface, table 1 reports its statistical properties
\begin{equation} k_{{a}} = \frac {1}{A_{t}} \int \left |h\right | \mathrm{d} A,\quad k_{{rms}} = \sqrt {\frac {1}{A_{{t}}} \int h^{\prime 2} \mathrm{d} A},\quad S_k =\left (\frac {1}{A_{{t}}} \int h^{\prime 3} \mathrm{\,d} A\right) /k_{{rms}}^3, \end{equation}
where
$A_t$
is the total area.
Additionally, a measure of the roughness density in a wall-parallel plane and in the cross-stream plane is reported using the plan solidity
$\lambda _p = A_p / A_t$
, where
$A_p$
is the plan area of roughness elements, and the frontal solidity
$\lambda _{{f}} = A_{{f}} / A_{t}$
. The reader can find a summary of the significance of each parameter in Chung et al. (Reference Chung, Hutchins, Schultz and Flack2021). Here, we note that for example a cube-like element oriented normal to the flow direction will have
$\lambda _{\!f}/\lambda _p=1$
, such that in our database cases CB_A and CB_S end up having the same values for both frontal and plan solidity, as well as the same statistical distribution. The real difference between the two is the arrangement, which currently lacks well-accepted metrics (Chung et al. Reference Chung, Hutchins, Schultz and Flack2021). However, this feature has been taken under consideration by recent models, noting that capturing the shielding effects between different arrangements is crucial to the prediction of drag (Yang et al. Reference Yang, Sadique, Mittal and Meneveau2016).
Throughout this study, we use the symbols
$u$
,
$v$
and
$w$
to denote the streamwise, wall-normal and spanwise velocity components and the decomposition of any variable is conducted using either the standard Reynolds decomposition (
$f=\bar {f}+f'$
) or the density-weighted (Favre) representation (
$f=\tilde {f}+f^{\prime \prime }$
), with
$\tilde {f}=\overline {\rho f}/\bar {\rho }$
. Here, mean flow statistics are averaged for at least 500
$\delta _{\textit{in}}/ u_{\infty }$
for all cases. In space, we average in the spanwise direction,
$z$
, and over a small window in the streamwise direction,
$x$
, with
$\lambda _x=1.44\delta _{\textit{in}}$
. The latter corresponds to four roughness periods for the case CB_A, or equivalently, approximately half of the mean boundary layer thickness after the onset of roughness. Below the roughness crest we perform an intrinsic average where only the fluid volume is considered.
3. Flow visualisations
We begin by investigating the instantaneous and mean flow in different planes, particularly focusing on the smooth-to-rough transition region.
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’.
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.
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.
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.
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.
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}$
.
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.
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).
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.
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
\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
\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.
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
\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.
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^+$
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
\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.
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).
(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)
\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^+$
.
(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
\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.
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.
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)
\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}
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.
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 }$
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.
7. Conclusions
In this study, we have presented a comprehensive DNS database of supersonic turbulent boundary layers consisting of four simulations spanning different types of prism-shaped roughness. The size of the computational domain offers the possibility to observe both the smooth-to-rough transition and the subsequent development of the turbulent boundary layer over a rough region, which reaches equilibrium with the new surface towards the end. Additionally, the extreme resolution allows the investigation of relevant turbulence statistics both on macroscopic (
$O(\delta _{99})$
) and microscopic (
$O(k)$
) levels, which have rarely been simultaneously considered in previous studies.
First, the analysis of instantaneous and averaged flow shows oblique shock waves at the smooth-to-rough transition generated by a sudden growth of the boundary layer. The intensity and related inclination are slightly higher for CB_R cases, followed by CB_S, DM_S and CB_A. In the rough part, the individual roughness elements do not protrude in the supersonic region, establishing a dynamics that is weakly affected by compressibility effects. Consequently, the development of the boundary layer thickness is different for each case both in terms of spatial recovery after the smooth-to-rough transition and of subsequent rate of growth. The growth of
$\delta _{99}$
in the rough region seems to be more pronounced for CB_R and CB_S, while similar trends are noted for CB_A and DM_S. This is also true for the distribution of the skin friction factor
$C_{\!f}$
, even though relative differences in
$\delta _{99}$
are more evident, it being the dominating factor that influences
$\textit{Re}_\tau$
. Similar behaviours are also noted in the evolution of the IBL, tracked with the algorithm proposed by Cogo et al. (Reference Cogo, Modesti, Bernardini and Picano2025). In general, we find that an equilibrium with the new surface in the present database is expected after roughly
$40 \delta _{99,\textit{ref}}$
from the onset of roughness, with
$\delta _{99,\textit{ref}}$
the thickness of the boundary layer incoming from the smooth wall. This is an important reference to inform future studies, but it is not yet clear if the friction Reynolds number or the Mach number can reduce or delay equilibrium.
Towards the end of the domain, the mean velocity field is preliminarily investigated by looking for outer-layer similarity across smooth and rough cases by means of the velocity defect. Here, we show that the scaling of Pirozzoli & Smits (Reference Pirozzoli and Smits2023) clearly outperforms classical definitions, and can generalise to both compressible flows and rough walls for all considered cases in combination with a suitable compressibility transformation (Van Driest Reference Van Driest1951). This results confirms that outer-layer similarity is observed for the mean velocity field in the present database across different rough-wall cases and their related smooth-wall reference. In the inner layer, mean velocity profiles of rough cases show the classical deficit compared with the smooth counterpart consistently with the observed added drag. Through the use of the compressibility transformation of Van Driest (Reference Van Driest1951) we evaluate
$\Delta U^+$
, which is used to provide an estimate of the equivalent sand-grain roughness height
$k_s$
for each topology. Outer-layer similarity is also observed for turbulent velocity fluctuations, but does not appear in the mean and fluctuating temperature statistics, in accordance to previous studies (Modesti et al. Reference Modesti, Sathyanarayana, Salvadore and Bernardini2022; Cogo et al. Reference Cogo, Modesti, Bernardini and Picano2025). Interestingly, the temperature field seems to further decrease its similarity to the smooth reference as drag increases (see for example case CB_R). On the contrary, when the mean and fluctuating temperature fields are observed in relation to their velocity counterpart through the classical Reynolds analogy theory (Gaviglio Reference Gaviglio1987; Huang et al. Reference Huang, Coleman and Bradshaw1995), a good agreement is found with reference smooth-wall data, with minor deviations in the region close to roughness sub-layer.
Finally, the overall drag generated from each geometry in the present configuration is investigated by partitioning pressure and viscous contributions, as well as considering volumetric sheltering between each element. We find that the CB_R and CB_S cases share a similar pressure contribution, which is justified by the similar wake pattern put in place by the staggered arrangement which dominates over the specific orientation of each element with the flow. However, case CB_R is able to produce more drag from viscous stress contributions on the lateral sides of roughness elements, which are inclined at an angle of
$45^\circ$
to the streamwise direction. The case CB_A is highly affected by shielding effects produced by the aligned pattern, which is able to reduce the pressure load, and thus the total drag. All cases with cubical elements show considerable peaks of very low and high friction on the bottom surface,
$\tau _S$
, denoting regions in which the flow suddenly loses and regains momentum, respectively. This scenario is not shared with the last case, DM_S, which shows very reduced wakes downstream of each element, avoiding highly recirculating flow. This feature, connected with the fact that this case has the highest planar solidity, promotes the ability of the flow to ‘skim over’ roughness elements, thus reducing the added drag.
Given these insights, we can better link the influence that the roughness sublayer dynamics has on macroscopic features of the flow, such as the increase in the boundary layer thickness and the shock-wave inclination at the smooth-to-rough transition. In this regard, we conclude that the growth rate of the boundary layer thickness is determined by the amount of drag produced by the shape and arrangement of the elements, which is strongly influenced by the mutual sheltering of neighbouring elements.
Future studies are needed to assess the relevant control variables that can promote compressibility effects in the roughness sublayer (e.g.
$\textit{Re}_\tau$
,
$T_w/T_r$
,
$M_\infty$
, positively or negatively skewed roughness). We believe that this is a key point to address in order to isolate the regime in which compressible turbulent boundary layer over rough surfaces can behave like their incompressible counterparts (if the mean property variations of the flow are taken into account with compressibility transformation), while pointing out specific factors that greatly enhance compressibility effects in the roughness sublayer (such as shock waves).
Acknowledgements
We acknowledge that the results reported in this paper have been achieved using the EuroHPC JU Extreme Scale Access Infrastructure resource Marenostrum 5 hosted at BSC-CNS, Barcelona, Spain, under project EHPC-EXT-2023E01-034. We also acknowledge the CINECA award under the ISCRA and EuroHPC initiatives (project EUHPC_E02_044), for the availability of high-performance computing resources on Leonardo booster.
Funding
We acknowledge financial support under the National Recovery and Resilience Plan (NRRP), Mission 4, Component 2, Investment 1.1, Call for tender No. 104 published on 2.2.2022 by the Italian Ministry of University and Research (MUR), funded by the European Union – NextGenerationEU – Project Title ADMIRE – CUP B53C24006770006 – Grant Assignment Decree No. 1401 adopted on 18/09/2024 by the Italian Ministry of Ministry of University and Research (MUR). This research received also financial support from ICSC – Centro Nazionale di Ricerca in ‘High Performance Computing, Big Data and Quantum Computing’, funded by European Union – NextGenerationEU. Davide Modesti is supported by the Air Force Office of Scientific Research under award number FA8655-25-1-7063.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are available upon reasonable request.
Open access
