Direct numerical simulation of supersonic turbulent flows over rough surfaces

Abstract We perform direct numerical simulation of supersonic turbulent channel flow over cubical roughness elements, spanning bulk Mach numbers $M_b=0.3$–$4$, both in the transitional and fully rough regime. We propose a novel definition of roughness Reynolds number which is able to account for the viscosity variations at the roughness crest and should be used to compare rough-wall flows across different Mach numbers. As in the incompressible flow regime, the mean velocity profile shows a downward shift with respect to the baseline smooth wall cases, however, the magnitude of this velocity deficit is largely affected by the Mach number. Compressibility transformations are able to account for this effect, and data show a very good agreement with the incompressible fully rough asymptote, when the relevant roughness Reynolds number is used. Velocity statistics present outer layer similarity with the equivalent smooth wall cases, however, this does not hold for the thermal field, which is substantially affected by the roughness, even in the channel core. We show that this is a direct consequence of the quadratic temperature–velocity relation which is also valid for rough walls. Analysis of the heat transfer shows that the relative drag increase is always larger than the relative heat transfer enhancement, however, increasing the Mach number brings data closer to the Reynolds analogy line due to the rising relevance of the aerodynamic heating.


Introduction
Supersonic boundary layers are ubiquitous in aerospace systems and they feature a more complex flow physics compared with their incompressible counterpart because of the strong coupling between the velocity and thermal fields and the presence of propagating disturbances such as shock waves. The flow physics of compressible turbulent boundary layers over smooth walls is an active topic of research in the community, and important advances have been achieved on several fundamental aspects, such as the relevance of genuine compressibility effects (Bradshaw 1977;Yu, Xu & Pirozzoli 2019), compressibility transformations (Coleman, Kim & Moser 1995;Foysi, Sarkar & Friedrich 2004;Modesti & Pirozzoli 2016;Trettel & Larsson 2016;Volpiani et al. 2020;Griffin, Fu & Moin 2021), Reynolds analogy relations (Zhang et al. 2014) and internal flows (Modesti, Pirozzoli & Grasso 2019;). On the contrary, much less is known on the effect of distributed surface roughness in high-speed flows.
A typical example of supersonic flow over rough surfaces is the ablative shield of re-entry vehicles, which experience extreme thermal loads due to the intense aerodynamic heating (Candler 2019). In order to maintain their structural integrity, re-entry vehicles are equipped with tiled or ablative thermal protection systems (TPS), and in both cases the flow may experience a rough surface. Tiled TPS are constituted of carbon or ceramic tiles with a square, diamond or hexagonal shape and the gaps between the tiles form a structured roughness pattern. Ablative TPS protect the underlying structure because the material undergoes pyrolysis and the gases that are generated in this process blow the boundary layer away from the surface. The surface ablates with a non-uniform recession rate, resulting in regular (Peltier, Humble & Bowersox 2016;Wilder & Prabhu 2019) or irregular (Kocher et al. 2017) distributed roughness patterns, depending on the type of material. Another example of compressible flows over roughness is found in transonic turbines, where the blades are subjected to erosion, forming irregular surface patterns. Additionally, irregular surface patterns may form at the wing leading edge due to icing.
Despite the technological relevance of compressible flows over rough surfaces, the vast majority of studies on rough walls are limited to the incompressible flow regime, which is well documented both from the experimental (Nikuradse 1933;Perry, Schofield & Joubert 1969;Raupach, Antonia & Rajagopalan 1991;Jiménez 2004;Flack & Schultz 2014) and numerical (Leonardi et al. 2003;Cardillo et al. 2013;MacDonald et al. 2016;Busse, Thakkar & Sandham 2017) perspectives. Incompressible flows over rough walls constitute a topic of active research with several challenges and unanswered questions . However, there is at least consensus on basic aspects, whereas this is not the case for the compressible flow regime. For instance, whether the flow experiences the surface as rough or smooth depends on the roughness Reynolds number k + = k/δ v , where k is the physical roughness height, δ v = ν w /u τ the viscous length scale, u τ = √ τ w /ρ w the friction velocity, τ w the drag per plane area and ν w , ρ w the kinematic viscosity and density of the fluid at the wall, respectively. In order to compare different types of roughness the equivalent sand-grain roughness height k s is often used which is the characteristic length scale that leads to matched drag between the surface pattern of interest and the sand-grain roughness originally studied by Nikuradse (1933). For roughness Reynolds numbers k + s 5 the flow is hydraulically smooth, that is the roughness does not induce any additional drag. As the roughness Reynolds number increases (5 k + s 80) the flow becomes transitionally rough and in this regime both viscous and pressure drag are important. For k + s 80 the flow becomes fully rough, meaning the skin-friction coefficient does not depend on the Reynolds number.
Another aspect that has considerable consensus in the community is the validity of the outer layer similarity hypothesis (Townsend 1980) over roughness, namely the fact that the outer flow is not directly affected by the surface topography, but it feels the roughness through the mean wall-shear stress. A practical consequence of outer layer similarity is that the drag variation induced by the roughness can be related to the streamwise momentum deficit, namely the downward shift of the mean velocity profile, or Hama roughness function (Hama 1954;Clauser 1956), U + . The main advantage is that, unlike the drag variation, the Hama roughness function is fairly independent of the Reynolds number and therefore it allows us to use numerical simulations and experiments to estimate the drag variation at higher Reynolds numbers, typical of engineering applications. In the case of compressible flows, even these fundamental aspects are not set yet. For instance, at the present stage it is not clear if the onset of the fully rough regime occurs for the same k + s ≈ 80 also in the supersonic case, or if, for instance, the additional wave drag modifies the transition to the fully rough regime. Similarly, the effect of shock and expansion waves induced by the roughness could challenge the validity of outer layer similarity, complicating the prediction of drag in high-speed flows. At present, only few experimental studies of supersonic flows over rough walls are available. Ekoto et al. (2008) performed experiments of turbulent boundary layer at free stream Mach number M ∞ = 2.86 over distributed cubic and diamond roughness elements (k/h ≈ 11, k + s ≈ 300), with the main goal to identify similarities and differences between these two roughness shapes. They used different experimental techniques, namely particle image velocimetry, schlieren photography, Pitot measurements and pressure sensitive paint, to measure the mean velocity and Reynolds shear stress. They found that the diamond roughness elements distort the flow more than the cubes and induce shock waves which propagate to the outer wall layer. Their mean flow statistics show similarity with incompressible flows over rough walls, such as a downward shift of the mean velocity profile. However, a systematic comparison with incompressible flow data was not performed. Peltier et al. (2016) performed schlieren photography and particle image velocimetry over diamond roughness elements (k + = 160, k + s = 600) at M ∞ = 4.9, and they also observed oblique shock and expansion waves propagating from the roughness to the outer flow. The mean velocity profile in the overlap region followed a logarithmic profile with a downward shift of U + ≈ 13. This value is similar to the one of the incompressible fully rough asymptote for k + s = 600 ( U + ≈ 12), implying negligible compressibility effects and the same drag variation with the respect to a smooth wall as in the incompressible flow regime. Kocher et al. (2018) performed particle image velocimetry and schlieren photography of a supersonic boundary layer at M ∞ = 2 over diamond and realistic roughness. Similarly to other studies, the schlieren images of their rough cases reveal the presence of shock waves generated by the roughness, propagating into the boundary layer up to the free stream region. The inner-scaled mean velocity profiles show the typical downward shift, which increases in the streamwise direction from transitional to fully rough. The authors reported different values of U + than in previous studies at higher Mach number (Peltier et al. 2016), and attributed the discrepancies to the limited validity of the concept of k + s across Mach numbers. To our knowledge, only one direct numerical simulation (DNS) study of supersonic flow over roughness has been performed so far (Tyson & Sandham 2013). The authors carried out DNS of compressible turbulent channel flow at bulk Mach numbers M b = 0.3, 1.5, 3.0 over two-dimensional wavy surfaces with different wave amplitudes and wavelengths spanning both the transitional and fully rough regimes. They evaluated the van Driest transformed velocity shift and found that it decreases with the Mach number, thus contradicting previous experimental studies which reported no compressibility effects on the Hama roughness function.
From this literature survey we find that supersonic flows over rough walls have so far been mainly studied experimentally, and even the experimental studies are very limited in number. Moreover, most studies are restricted to adiabatic or nearly adiabatic wall conditions, whereas in realistic engineering applications the heat flux plays a relevant role (Bowersox 2007). As a result several fundamental aspects, which are established in the incompressible flow regime, have barely been addressed in the framework of high-speed flows. To cover this gap, in this work we present novel results from DNSs of supersonic turbulent channel flow over cubical roughness elements at various Mach and Reynolds numbers, covering both the transitionally and fully rough regime. The manuscript is organized as follows. We introduce the numerical methodology in § 2, where a description of the computational set-up and roughness configuration is provided. Results are reported in § 3 where velocity and thermal statistics are discussed, with a focus on the assessment of compressibility effects on the roughness function and the outer layer similarity. Conclusions are finally given in § 4.

Physical model
We solve the compressible Navier-Stokes equations for a perfect heat-conducting gas, where u i , i = 1, 2, 3 are the velocity components in the streamwise, wall-normal and spanwise directions, respectively, ρ, p, T are the the fluid density, pressure and temperature, respectively. The total energy per unit mass is denoted as E = c v T + u i u i /2, H = E + p/ρ is the total enthalpy, whereas σ ij and q j are the viscous stress tensor and heat flux vector, The dependence of the viscosity coefficient on temperature is accounted for through a power law with exponent 0.76 and k = c p μ/Pr is the thermal conductivity with Prandtl number, Pr = 0.72. We consider the channel flow configuration where the flow between two infinite isothermal walls is driven in the streamwise direction by a body force f , which is evaluated at each time step in order to discretely enforce a constant mass-flow rate with the corresponding power spent added to (2.1c). Additionally, a bulk cooling term Φ is added to the total energy equation to control the bulk flow temperature (Yu et al. 2019), that is kept constant during the simulation. In particular, Φ is evaluated at each time step such that only a fraction Θ of the bulk flow kinetic energy is converted into wall heat flux, namely T w = T b [1 + 0.5Θ(γ − 1)rM 2 b ], where γ = 1.4 is the heat capacity ratio, r = 0.89 the recovery factor, T w is the wall temperature, number of the flow, and ρ b , u b and T b are the bulk flow density, velocity and temperature (2.4a-c) 2.2. Numerical method The Navier-Stokes equations are solved using the solver STREAmS (Bernardini et al. 2021), which has been extended with immersed boundary capabilities. The nonlinear terms in the Navier-Stokes equations are discretized using a hybrid energy-conservative shock-capturing scheme in locally conservative form. In shock-free regions, we use a sixth-order energy consistent flux, which guarantees a discrete conservation of the total kinetic energy in the limit case of inviscid incompressible flow (Pirozzoli 2010). Shock-capturing is achieved through Lax-Friedrichs flux vector splitting, where the characteristic fluxes are reconstructed at the interfaces using a fifth-order, weighted essentially non-oscillatory reconstruction (Jiang & Shu 1996). To determine the local smoothness of the solution and switch between the central and shock-capturing scheme, a classical shock sensor is adopted (Ducros et al. 1999). The viscous terms are expanded into a Laplacian form and approximated with sixth-order central finite-difference formulae to avoid odd-even decoupling phenomena. Time stepping is carried out by means of Wray's three-stage third-order Runge-Kutta scheme (Spalart, Moser & Rogers 1991).
The complexity of the roughness geometry is handled using a ghost-point-forcing immersed boundary method to treat arbitrarily complex geometries (Piquet, Roussel & Hadjadj 2016;De Vanna, Picano & Benini 2020). The geometry of the solid body is provided in OFF format for three-dimensional objects, and the computational geometry library CGAL (The CGAL Project 2021) is used to perform the ray-tracing algorithm. This allows us to define the grid nodes belonging to the fluid and to the solid, and to compute the distance of each point from the interface. To retain the same computational stencil close to the boundaries, the first three layers of interface points inside the body are tagged as ghost nodes. For each ghost node, we identify a reflected point along the wall normal, lying inside the fluid domain. We interpolate the solution at the reflected point using a trilinear interpolation and use the values at the reflected points to fill the ghost nodes inside the body to apply the desired boundary condition. An extensive description of the algorithm is available in the work by De Vanna et al. (2020).

Flow configuration and computational parameters
In this work we consider rough walls formed by cubic elements of side k, which are representative of the structured roughness patterns forming over ablative surfaces. The cubes are placed specularly over the bottom and top channel walls with a spacing in the wall-parallel directions equal to 2k, as shown in figure 1. We develop a DNS dataset of geometrically increasing roughness, namely we keep constant the roughness size with respect to the channel half-width (k/h = 0.08), while increasing the friction Reynolds number from Re τ ≈ 500 to Re τ ≈ 1000, corresponding to roughness Reynolds numbers k + ≈ 40 and k + ≈ 80. For each Reynolds number, we consider two supersonic cases at bulk Mach number M b = 2, M b = 4 and one additional case at Re τ ≈ 500, M b = 0.3. For each rough wall case, we carry out a companion smooth wall simulation at matching bulk Mach number and approximately matching friction Reynolds number, for a total of 10 simulations, as reported in table 1.
The DNS are carried out in a rectangular box with size 6h × 2h × 3h where h is the channel half-height. The mesh spacing is constant in the wall-parallel directions, and an error-function mapping is used to cluster mesh points towards the roughness crest (i.e. y = k). A smooth transition between the two distributions is obtained using a spline function (Gregory & Delbourgo 1982). We use a standard DNS resolution for smooth wall cases, whereas a much finer mesh was adopted for the rough wall cases, with approximately 40 mesh points per roughness element in each direction (i.e. x ≈ z ≈ k/40), as shown in table 1. This resolution has been chosen on the basis of a mesh sensitivity analysis reported in the Appendix.
All computations are initiated with a parabolic velocity profile with superposed random perturbations, and with uniform values of density and temperature. As for the boundary conditions, periodicity is enforced in the homogeneous wall-parallel directions, and no-slip isothermal (Θ = 0.35) conditions are imposed at the channel walls. Here Θ is a key parameter and estimating its value in a practical flow configuration is not easy because it may depend on several aspects, such as the type of material of the solid wall, the type of cooling (ablation, tiles, transpiration) and of course Mach and Reynolds number. Yu et al. (2019) performed smooth wall simulations with Θ in the range 0.25-1 and reported increasing compressibility effects for colder walls. Our value of Θ corresponds to a relatively cold wall which is representative of the strong cooling conditions found on heat shields, and is in the range studied by Yu et al. (2019).
In the following we use both the Favre (·) and Reynolds (·) ensemble averages (i.e. averages in time and over the roughness period) defined as Ensemble averages are further decomposed into a mean and dispersive component, where angle brackets · denote intrinsic averages in the wall-parallel directions. This triple decomposition allows us to split the total Reynolds stress tensor into a turbulent and a dispersive component, which will be analysed in the next section.  (b) (a) Figure 2. Instantaneous flow field for cases S4_1000 and R4_1000. The longitudinal and cross-stream planes show the temperature field, whereas streamwise velocity fluctuations are reported in the wall-parallel planes at 10 wall units from the wall (S4_1000) and from the roughness crest (R4_1000).

Flow visualizations
We begin by inspecting the instantaneous flow fields of representative flow cases S4_1000 and R4_1000 in figure 2, where the wall-parallel planes show contours of the streamwise velocity, and the wall-normal planes display the temperature field. The near-wall flow is populated by streaks both for the smooth and rough cases, but the roughness substantially changes the flow organization. The spacing between the roughness elements is large enough (2k ≈ 160δ v ) to interfere with the canonical streaks spacing (≈100δ v ), and high-speed streaks are predominantly locked between two spanwise-adjacent roughness rows, whereas low-speed streaks tend to be located on top of the roughness crests. The wall-normal planes highlight large bulges of hot fluid raising from the channel walls and protruding almost to the channel centre which are particularly evident in the cross-stream plane, whereas they get skewed in the streamwise plane under the effect of the mean shear. We note that the space between the roughness elements is primarily filled by hot fluid, whereas close to the smooth wall we can identify cold fluid patches down to the wall. Hence the roughness reduces the temperature fluctuations close to the wall. The rough wall case shows much larger temperature excursions, both at the wall and at the channel centre, which is a hint of possible effects of the roughness on the outer layer.
Another view of the instantaneous flow field is presented in figure 3, where we show the instantaneous streamwise Mach number u/c in the wall-normal planes for the smooth and rough flow cases S4_1000 and R4_1000. The smooth flow field shows a canonical organization where high-speed structures protrude down to the wall and low-speed eruptions extend towards the channel core, which is particularly evident in the cross-stream plane, figure 3(b). Over the smooth wall the flow is supersonic down to the viscous sublayer, as clear from the orange line indicating sonic conditions.
The flow organization of the rough wall case R4_1000 is substantially different from the smooth wall. The channel is essentially divided into two regions, where the high-speed core is encased by low-speed fluid close to the wall. The high-speed bundle at the channel centre appears less penetrated by the low-speed ejections from the wall as compared with the smooth wall, and the maximum Mach number is considerably higher, with localized  maxima up to u/c ≈ 8. The sonic isoline (orange) in figure 3(c,d) shows that the roughness crests are most frequently subsonic but inrushes of high-speed fluid can penetrate down to the roughness troughs and the local Mach number is supersonic.
As common for turbulent flows the instantaneous field can be very different from the mean, and in figure 4 we show the average Mach numbers associated with the streamwise and wall-normal velocities, and the mean temperature close to the roughness element. The mean flow is effectively three-dimensional close to the roughness, especially below the crest, where we note a non-uniform spatial distribution both for the Mach number and the temperature fields, whereas the mean flow becomes homogeneous in the wall-parallel directions for y 3k.
The streamwise Mach number (figure 4a-e) is affected by both Re b and M b . Increasing the Reynolds number leads to two competing effects, namely a fuller velocity profile, but also a larger k + and therefore a larger upward shift of the near wall cycle. In this case the latter effect is dominating. Increasing the bulk Mach number brings high speed fluid closer to the roughness, but for the cases under scrutiny the mean flow at the crest is never supersonic, although it reachesũ/c ≈ 0.8 for flow cases at M b = 4 (figure 4c,d).
The Mach number associated with the wall-normal velocity is largely subsonic and it does not exceedṽ/c ≈ 0.2, featuring extended separated flow regions both upstream and downstream of the cube because elements are in each other wakes.
The maximum mean temperature occurs at the crest height, but it is highly localized at the upstream corner of the roughness element, in correspondence of the stagnation point where the kinetic energy of the flow is converted into internal energy.

Compressibility effects on added drag
A crucial aspect of supersonic flows over roughness is whether the flow experiences the same added drag with the respect to the smooth wall across Mach numbers. The drag variation with respect to the smooth wall can be expressed as  where is the mean velocity at the channel centreline, R c = ρ c / ρ w , and the subscript s denotes the smooth wall. Equation (3.1) shows that the drag variation can be related to the velocity deficit at the channel centre and in the logarithmic region, where the second identity is valid if outer layer similarity holds. Hence, understanding compressibility effects on DV is intrinsically related to finding a compressibility transformation for the mean velocity. Moreover, a notable difference with respect to the incompressible case is the presence of the density ratio R c /R cs in (3.1), which requires knowledge of the density profile. After the celebrated velocity transformation proposed by van Driest (1951), several other compressibility transformations have been developed (Trettel & Larsson 2016;Volpiani et al. 2020) aiming at improving the accuracy with reference incompressible flow data, especially for isothermal walls.
A generic compressibility transformation for the mean velocity can be expressed using stretching functions for the velocity and wall-normal coordinate, where f I and g I for different transformations are reported in table 2. The transformed wall coordinate also allows us to define the equivalent incompressible Reynolds number,  Recent compressibility transformations (Trettel & Larsson 2016;Volpiani et al. 2020) have proved their ability to account for compressibility effects and the transformed mean velocity profile over smooth walls exhibits the canonical logarithmic region, where κ ≈ 0.39 is the von Kármán constant and B ≈ 5.1. In the presence of wall roughness the mean velocity profile is characterized by a downward shift, To verify the accuracy of compressibility transformations, figure 5 shows the mean velocity profile, untransformed ( The smooth wall profiles transformed according to V show a very good agreement with the nearly incompressible flow case (dashed line with plus sign). On the contrary VD and especially TL transformations are less accurate. Rough wall cases are all characterized by the typical downward shift of the mean velocity, but we find visible compressibility effects. For instance, the supersonic cases R2_500 (red long-dashed line with solid circle) and R4_500 (red long-dashed line with solid box) show a different downward velocity shift with respect to the nearly incompressible case R03_500 (solid line with plus sign), although they share the same k + . We note that none of the compressibility transformations is able to account for this discrepancy. Hence, the canonical definition of roughness Reynolds number does not fully characterize the flow regime in the case of compressible flows.
As commonly done in the low-speed regime (Ibrahim et al. 2021), rough-wall velocity profiles have been shifted to account for the effective wall origin of the flow. Over the years, several methods have been proposed to find the virtual origin, which have been recently reviewed by Chung et al. (2021). In this study we opt for a simple option and chose d = 0.9k for all cases, which allows us to substantially reduce the uncertainty in the measure of the velocity deficit. To justify this choice in figure 6, we also report the difference between smooth and rough wall for the transformed velocity profile u V as a function of the wall distance, after shifting (dash-dotted) and before shifting (solid). The figure shows that the velocity deficits for the shifted profiles are much flatter, thus increasing the confidence in evaluating U + I . Throughout this work we evaluate the velocity shifts at the nominal edge of the logarithmic region y + I = 0.3Re τ I .  To understand the role of the roughness Reynolds number we cast (3.6a,b) as where k I is the equivalent incompressible roughness height and C(k I ) is a roughness dependent function. In the incompressible regime there is no ambiguity in the definition of k I = k, whereas in the compressible case multiple definitions are possible and we consider two options. The first naturally stems from the compressibility transformation for the wall-normal coordinate (3.4), k I = y I (k).
(3.8) Equation (3.8) has the advantage of being consistent with the transformed velocity shift U + I , however, it might be difficult to estimate it from experimental data. For this reason we also consider the following length scale: whereν k is the mean kinematic viscosity evaluated at the roughness crest. Figure 7 shows the velocity shift as a function of the equivalent sand-grain roughness Reynolds number. In figure 7(a) we report the untransformed velocity shift ΔU + as a function of k + s , whereas figures 7(b-d) show the transformed shift U + I as a function of the corresponding k + sI . The relationship between the physical roughness height and k s for the present database was obtained by matching the nearly incompressible case R03_500 with the fully rough asymptote, which provides k s /k = 1.9. We also recall that the concept of sand-grain roughness Reynolds number is only valid in the fully rough regime, whereas in the transitionally rough regime different roughness geometries exhibit a different trend with k + s . Thakkar, Busse & Sandham (2018) performed DNS of incompressible channel flow over grit blasted surface which shows a U + very similar to the data of Nikuradse (1933) also in the transitionally rough regime. However, this is in general not expected especially for structured roughness patterns such as bars and cubes . The untransformed U + shows visible discrepancies with respect to the incompressible sand-grain roughness data of Nikuradse (1933, crosses), and also when compared with the more recent DNS data of Abderrahaman-Elena, Fairhall & García-Mayoral (2019) for a similar roughness geometry. This is particularly evident for the transitionally rough cases, which present substantially different values of U + , although they share the same k + s , figure 7(a). As for the effect of compressibility, we observe a slightly better agreement for ΔU + VD in figure 7(b), but transitionally rough data at M b = 4 (filled square symbol) still differ from the incompressible data of Abderrahaman-Elena et al. (2019). In figure 7(b) we also report experimental data of supersonic boundary layer over rough walls, extracted from the review of Bowersox (2007), which are in good agreement with incompressible flow data and with the fully rough asymptote. However, experiments have been carried out in nearly adiabatic wall conditions (Goddard 1959;Latin & Bowersox 2000;Ekoto et al. 2008), whereas the present DNS dataset is characterized by strong cooling at the wall, thus the thermodynamics properties variation at the crest is more significant and this is not accounted for by k + s . Figure 7(c,d) show the transformed velocity shifts U + TL and U + V as a function of the respective roughness Reynolds numbers (3.8). These transformations show similar accuracy to van Driest, and they only partially account for compressibility effects, whereas differences are still evident for the case R4_500 (filled square).
In figure 8 we report the velocity shift as a function of k + * s , as defined in (3.9). The figure shows that the scaling based on the viscosity at the crest substantially improves the agreement with incompressible flow data, when compared with the transformed roughness height k + I in figure 7. When using k + s * , we find an excellent agreement with incompressible data for all compressibility transformations ( figure 8b-d), although a slightly better result is visible Open symbols indicate cases at Re τ ≈ 1000, filled symbols at Re τ ≈ 500 for M b = 2 (circles) and M b = 4 (squares). In panel (b) experimental data of supersonic boundary layer are reported: Reda, Ketter & Fan (1975, upward triangle); Berg (1979, pentagon); Goddard (1959, downward triangle); (Latin & Bowersox 2000, star); (Ekoto et al. 2008, diamond). Incompressible data are also reported: experiments of Nikuradse (1933, crosses)  for U + V . The idea behind the length scales k * and k I is similar, namely they attempt to account for density variations between the crest and trough of the roughness, and they both seem to help the comparison with incompressible data. The transformed roughness height k I has the advantage to be consistent with its velocity transformation U + I , however, k * leads to slightly more accurate results, besides being easier to compute.
We believe that the definition of a relevant roughness Reynolds number k + I is a key aspect of compressible flows over roughness, but this topic has been often glossed over by previous studies, who might have involuntarily (but erroneously) incorporated this effect into k + s . For instance, Hill, Voisinet & Wagner (1980) reported differences between k s of their supersonic sand-grain roughness and Nikuradse's data, and attributed the discrepancies to uncertainty in the roughness manufacturing. Another example is the work by Berg (1979), who carried out experiments of bar roughness at M ∞ = 6 and computed k s by matching his U + D with Nikuradse's data. In contrast our data show that calculating k s /k by matching the transformed velocity shift U + I with the fully rough asymptote is not enough, because the roughness Reynolds number is also influenced by compressibility effects.

Turbulent fluctuations
We analyse the effect of surface roughness on turbulent fluctuations by decomposing the Reynolds stress tensor into a turbulent and a dispersive component as in (2.7). Figure 9 shows the streamwise Reynolds stress component τ 11 as a function of the viscous-scaled distance from the wall including both rough (solid lines) and smooth wall simulations (dashed lines). To facilitate comparison across Mach numbers, profiles are reported as a function of the transformed coordinate y + TL , which is usually regarded as the correct wall distance for the Reynolds stresses Huang, Coleman & Bradshaw 1995).
Smooth wall data are characterized by a lack of universality even in the near wall region where the peak of τ t 11 increases with both the Reynolds and Mach number. The former effect, widely reported also in the low-speed regime, is a consequence of the increasing relevance of outer layer structures with the Reynolds number, which provide additional energy at the low wavenumbers due to the imprinting effect on the near wall cycle (Hutchins & Marusic 2007). Instead, the latter is typically considered a genuine compressibility effect that cannot be accounted for by the density scaling (Pirozzoli & Bernardini 2011).
We note that the roughness disrupts the near wall cycle and the inner peak of τ t 11 is replaced by a much milder growth of the velocity fluctuations. The streamwise turbulent stress component has a maximum right above the crest, whose intensity is affected by both Mach and Reynolds number. At approximately matching Re τ TL , we observe that the peak of τ t 11 increases with the Mach number indicating compressibility effects in proximity of the crest. We also note that at fixed Mach number the peak of τ t 11 increases with Re τ TL , consistently with the incompressible flow regime.
The dispersive and turbulent stresses have approximately the same intensity below the roughness crest, whereas τ d 11 dominates in the region y ≈ k. The peak of τ d 11 occurs right below the roughness crest, and again a clear effect of compressibility is visible, with higher values for increasing Mach number. A rapid decay of the dispersive stress is observed moving away from the wall, implying that the mean flow becomes spatially homogeneous above the roughness crest.
Additional insight on the turbulence fluctuations can be gained from figure 10, where we report the turbulent and dispersive components of the Reynolds shear stress. Differently from the streamwise component, the density scaling is able to account for compressibility effects, and the main differences between cases are associated with the Reynolds number. As in the incompressible flow regime, the main effect of the roughness is to shift upwards the near wall cycle, thus the peak of τ t 12 occurs above the crest. As for τ t 11 , we observe a good match between rough and smooth wall in the channel core, which supports the validity of outer layer similarity for the velocity fluctuations. The dispersive component presents a maximum at the roughness crest, with only minor effect of the Mach number.
To better assess compressibility effects on turbulence in figure 11 we plot the turbulent Mach number (τ t 11 + τ t 22 + τ t 33 )/ c as a function of the wall distance for all flow cases. For the nearly incompressible case M t < 0.05, indicating genuine incompressible turbulence, as expected. For supersonic smooth flow cases at M b = 2 (figure 11a,b) the peak of M t does not exceed 0.3, which is often considered the upper edge above which compressibility effects become relevant. Smooth flow cases at M b = 4 (figure 11c,d) have a slightly higher peak M t ≈ 0.5, indicating the increasing role of compressibility in the buffer layer, which is also reflected by the higher peak of the streamwise velocity fluctuations previously discussed.
Rough wall cases show a substantially different trend with respect to the smooth walls at all Mach numbers. Besides the obvious shift of the near wall peak, rough wall cases present much higher values of M t and also a different shape of the profiles in the outer region, which cannot be accounted for by a simple virtual origin shift. This observation suggests that the roughness is able to affect the thermal field in the bulk flow region, as investigated in the next section, where we focus on the thermal flow statistics. y + (d) Figure 11. Turbulent Mach number M t = (τ t 11 + τ t 22 + τ t 33 )/ c for flow cases R2_500, S2_500 (a), R2_1000, S2_1000 (b), R4_500, S4_500(c), R4_1000, S4_1000 (d). Nearly incompressible flow data S03_500, R03_500 are also reported in panels (a,b). Symbols in table 1.

Thermal statistics
One of the peculiar features of compressible supersonic flows is the active coupling between momentum and heat transport. The study of heat transfer over rough walls has been limited to the incompressible flow regime where temperature can be considered a passive scalar (MacDonald, Hutchins & Chung 2019; Peeters & Sandham 2019), which allows direct comparison with the velocity. In the compressible case, feedback coupling between the temperature and momentum equations leads to a substantially different effect of the roughness on the thermal field.  Figure 12 shows the mean temperature profile normalized with its value at the centreline T c , as a function of the inner scaled wall distance. As typical of supersonic boundary layers with isothermal cold walls, the smooth-wall profiles exhibit a peak of the static temperature within the wall layer, at y + ≈ 5, with a positive temperature gradient at the wall, meaning that the fluid releases heat to the wall.
The rough wall temperature profiles retain the same qualitative trend of the smooth wall cases, but they also present relevant differences. The maximum of the temperature is now located approximately at the roughness crest, and it is higher than for the smooth wall. The rough wall profiles show higher values of T /T c over the entire wall layer, suggesting that outer layer similarity does not hold for the temperature field.
In order to understand why outer layer similarity holds for the mean velocity but not for the mean temperature we investigate the temperature-velocity relation. The coupling between the momentum and energy equations gives rise to a well known quadratic relationship between temperature and velocity (Smits & Dussauge 1996), and several analytical expressions have been proposed which accurately approximate this functional form (Walz 1959;Modesti & Pirozzoli 2016). A temperature-velocity relation for isothermal walls has been proposed by Zhang et al. (2014), where T rg = T c + r g ũ c 2 /(2C p ), r g = 2C p (T w − T c )/ ũ c 2 − 2Prq w /( ũ c τ w ), and the subscript c indicates quantities at the centreline. Figure 13 shows the mean temperature as a function of the mean velocity for all flow cases, compared with the analytical relation by Zhang et al. (2014). For flow cases at M b = 2 (figure 13a,b) the agreement between DNS data and (3.10) is excellent, for both rough and smooth wall flow cases, whereas some discrepancies are visible at higher Mach number (figure 13c,d). Despite these differences, figure 13 shows that the quadratic relation between velocity and temperature is a robust flow feature and it is also valid for rough walls. Therefore, both for rough and smooth walls it is possible to write where a, b are parameters which depend on Mach, Reynolds number, Prandtl number and the boundary conditions. We further assume that the compressible mean velocity profiles of smooth and rough walls exhibit a logarithmic region, where κ , B are not universal as in the incompressible case but depend on the Mach number and boundary conditions. This is essentially equivalent to assuming that outer layer similarity holds for the untransformed mean velocity profile, which is a rather accurate approximation as shown in figure 5. Substituting (3.12a,b) into the respective temperature-velocity relation (3.11), and taking the difference we obtain namely, unlike for the mean velocity, the difference between the smooth wall temperature and the rough wall temperature in the logarithmic region is a function of the wall distance. This clearly breaks the outer layer similarity for the mean temperature. Hence, for compressible flows over rough walls the lack of outer layer similarity can be traced back to the nonlinear relationship between the mean temperature and the mean velocity, which is a direct consequence of the aerodynamic heating. Additionally, we report the temperature fluctuations normalized by the friction temperature T τ , where q w is the heat flux per plane area. Figure 14 shows the temperature fluctuations in viscous units for smooth-and rough-wall cases. For the smooth wall we note small values of the temperature fluctuations up to y + ≈ 10, indicating a thicker viscous sublayer for the temperature than for velocity, and as a result the near wall peak of the temperature is shifted to y + ≈ 25-30.
The temperature fluctuations of the rough wall cases show a substantially different behaviour from the smooth wall. The roughness disrupts the near wall cycle and it shifts turbulence upwards, thus temperature fluctuations start increasing above the roughness crest.
The rough wall cases show the emergence of a peak of the temperature fluctuations in the outer layer, which is particularly evident for the high Reynolds number cases, figure 14(b). This peak cannot be associated with an upward shift of the near wall cycle because it occurs well into the outer layer and it increases with the Reynolds number. Therefore, we conclude that, unlike for the velocity field, the roughness is able to modify the temperature fluctuations in the outer layer.  Figure 13. Temperature-velocity relation for flow cases R2_500, S2_500 (a), R2_1000, S2_1000 (b), R4_500, S4_500 (b), R4_1000, S4_1000 (d). Dashed lines denote smooth wall cases and solid lines rough wall cases. Symbols in table 1. The dotted (smooth) and dash-dotted (rough) lines refer to (3.10). The analysis of the one-point statistics of the temperature field reveals that T is substantially affected by the presence of the roughness, therefore it is worth comparing the wall heat flux of smooth and rough walls.
We report the wall heat flux in terms of Stanton number, where T r is the recovery temperature based on the bulk Mach number, and we compare the relative heat transfer increase with the relative drag increase using the skin-friction coefficient . 942 A44-20 Half-filled symbols refer to the values obtained from the empirical correlation (3.16) by Hill et al. (1980) for the corresponding DNS data. Symbols style in table 1.
In figure 15(a) we report the ratio St/St s as a function of C f /C f s and compare the present data with the case of forced convection in incompressible flows (MacDonald et al. 2019) and to the engineering correlation of Hill et al. (1980), developed for hypersonic flows, where β = 0.4. Note that (3.16) requires the explicit knowledge of the skin-friction coefficient of the rough case and therefore in figure 15(a) each DNS flow case (symbols in table 1) has a corresponding point evaluated with this empirical correlation (half-filled symbols). We observe that both St/St s and C f /C fs > 1, namely the rough surface experiences a higher drag and heat transfer than the baseline smooth wall. The dashed line represents the Reynolds analogy line, that is the case in which the drag increase is equal to the heat transfer increase. All data lie below the axis bisector, implying that drag increases more than heat transfer, due to the additional pressure drag induced by the roughness, which does not have an equivalent for the temperature.
Our nearly incompressible flow case follows quite closely the data of egg-carton roughness by MacDonald et al. (2019), although these cases refer to forced thermal convection. A very similar trend holds for supersonic flow cases at M b = 2 (circles), which are very close to incompressible data because aerodynamic heating plays a minor role at this Mach number. Instead, flow cases at M b = 4 (squares) becomes closer to the Reynolds analogy line, especially for the lower k + , which indicates that aerodynamic heating tends to compensate for the additional pressure drag caused by the roughness. In figure 15(a) we also plot the values obtained from Hill's correlation (3.16) (half-filled symbols), which show that the empirical formula is very accurate for cases at M b = 4 as symbols essentially lay on top of DNS data, whereas we observe large discrepancies at M b = 2, especially at higher k + .
In figure 15(b) we report the Reynolds analogy factor (2St/C f ) as a function of the roughness Reynolds number, both for rough and smooth wall (k = 0) cases, and we compare DNS data with the experiments of Hill et al. (1980) (black circles) and to the empirical formula (3.16). Smooth wall data from DNS have a Reynolds analogy factor larger than one, whereas for experimental data are closer to unity, or lower, which can be probably associated with the much higher Reynolds number of the experiments (Re τ ≈ 5000). Despite the very different flow conditions of the experiments (higher Reynolds, M ∞ = 9, different roughness geometry) we note a similar trend with DNS for increasing k + , indicating that the Reynolds analogy factor mainly depends on the roughness Reynolds number.

Conclusions
We have presented DNSs of supersonic flows over cubical roughness elements, which are representative of the structured patterns found on the ablative shields of high-speed vehicles. We used these data to assess basic aspects of compressible flows over rough walls, which are well established in the incompressible regime, whereas they have been less explored for supersonic flows.
First, we analysed compressibility transformations and their ability to collapse the mean velocity shift onto the incompressible rough asymptote. The DNS data show that using the transformed U + I alone does not allow us to fully account for compressibility effects because of the ambiguity in defining a relevant roughness Reynolds number for the supersonic flow cases. We find important compressibility effects on the standard roughness Reynolds number k + , which therefore is not a good candidate for comparing flows across Mach numbers. These effects have been always ignored in previous studies of supersonic flows over rough walls or, most probably, involuntarily incorporated in the equivalent roughness Reynolds number k + s . We propose a transformed roughness Reynolds number 942 A44-22 k + * which accounts for viscosity variations between the crest and the trough, yielding excellent agreement with incompressible data, especially when used in conjunction with the recent velocity transformations by Volpiani et al. (2020).
Second, velocity statistics of rough and smooth walls show a very good match away from the wall, which supports the validity of outer layer similarity also for supersonic flows over rough walls. However, this does not seem to hold for the thermodynamic variables, and the temperature of rough and smooth wall cases is substantially different in the whole wall layer. We show that this can be traced back to the quadratic temperature-velocity relation which characterizes compressible flows both on smooth and rough walls.
Finally, we assessed the accuracy of the engineering correlation developed by Hill et al. (1980) for the relative Stanton number increase with respect to a smooth wall. The accuracy of the correlation increases with the Mach number and it yields nearly perfect results already at M b = 4, even if it was developed using sand-grain roughness and much higher Mach and Reynolds numbers.
This study is a first step towards the understanding of compressibility effects in flows over rough walls but it only explored a small part of the large parameter space characterizing this problem, in particular a single roughness geometry has been investigated. Future work will be dedicated to the study of different roughness geometries, as diamonds elements and sand-grain roughness, and also different wall-thermal conditions in both channel and boundary layer configurations.
DNS data are available at http://doi.org/10.4121/19403864. Immersed boundary methods are known to be sensitive to the mesh resolution and often show non-monotonic grid convergence. In this study we carried out simulations using streamwise and spanwise mesh spacings which guarantee approximately 40 points inside the roughness element (i.e. x ≈ z ≈ k/40), and used y + ≈ 0.5 below the roughness crest. In order to prove that this mesh resolution is adequate, we carried out two simulations of turbulent channel flow at bulk Mach number M b = 0.2 and bulk Reynolds number Re b = 10293 using the same computational domain of the main dataset. The first simulation M02A has N x × N y × N z = 3088 × 384 × 1536 mesh points, and the second M02B N x × N y × N z = 4096 × 384 × 2048, corresponding to 40 and 60 mesh points per roughness element, respectively. Figure 16 shows a comparison between M02A (solid) and M02B (circles) which shows that the two solutions are in nearly perfect agreement. We observe only minor deviations for the wall-normal and spanwise Reynolds stress components, but these are limited to less than 1 %, and therefore, the mesh refinement study confirms that using 40 points per roughness element allows us to achieve grid independence of the flow statistics.