Turbulent boundary-layer flow over regular multiscale roughness

Abstract

In this experimental study, multiscale rough surfaces with regular (cuboid) elements are used to examine the effects of roughness-scale hierarchy on turbulent boundary layers. Three iterations have been used with a first iteration of large-scale cuboids onto which subsequent smaller cuboids are uniformly added, with their size decreasing with a power-law as the number increases. The drag is directly measured through a floating-element drag balance, while particle image velocimetry allowed the assessment of the flow field. The drag measurements revealed the smallest roughness iteration can contribute to nearly 7 $\%$ of the overall drag of a full surface, while the intermediate iterations are responsible for over $12\,\%$ (at the highest Reynolds number tested). It is shown that the aerodynamic roughness length scale between subsequent iterations varies linearly, and can be described with a geometrical parameter proportional to the frontal solidity. Mean and turbulent statistics are evaluated using the drag information, and highlighted substantial changes within the canopy region as well as in the outer flow, with modifications to the inertial sublayer (ISL) and the wake region. These changes are shown to be caused by the presence of large-scale secondary motions in the cross-plane, which itself is believed to be a consequence of the largest multiscale roughness phase (spacing between largest cuboids), shown to be of the same order of magnitude as the boundary-layer thickness. Implications on the classical similarity laws are additionally discussed.

1. Introduction

Understanding turbulent flows over rough walls is paramount to a wide range of practical problems. In engineering applications, they manifest in internal flow systems such as in pipes and ducts, marine transportation or turbomachinery (Shockling, Allen & Smits 2006; Bons 2010; Monty et al. 2016). They are also ubiquitous in natural environments such as in the atmospheric boundary-layer (ABL), urban areas or in fluvial processes (Nezu, Tominaga & Nakagawa 1993; Garratt 1994; Cheng & Castro 2002; Nikora & Roy 2012). In all the aforementioned examples, roughness plays a fundamental role, influencing the overall drag, turbulent mixing and transport properties (Jiménez 2004).

One of the long-standing questions remains the characterisation of the drag by means of the surface properties; however, there is a plethora of arrangements and shapes under which roughness can be present. Following Grinvald & Nikora (1988) and Stewart et al. (2019) definitions, roughness can be classified into two categories. It can be considered discrete (regular), when presented as an amalgamation of discrete elements defined by a combination of linear length scales such as height, width and pitch (e.g. cube roughness). Alternatively, it can be interpreted as continuous (irregular), when the roughness is randomly distributed among a wider range of scales and defined by its statistical moments such as standard deviation, skewness and flatness (e.g. naturally corroded pipes). The seminal work of Nikuradse (1933) followed by the comprehensive reviews by Schlichting (1937) and Colebrook et al. (1939) served as the basis for the early drag predictions performed by Moody (1944) for pipe flows, using the equivalent sandgrain roughness height $h_{s}$. The latter has since become the common currency for predicting drag, with the aim of correlating surface properties to this aerodynamic roughness length scale $h_{s}$. Several studies have attempted to provide relations that would predict drag at high Reynolds numbers based on either laboratory measurements or numerical simulations. However, the existence of a universal correlation that would generically represent any rough surface is unlikely, due to its simultaneous dependence on various parameters such as front and plan solidities (Van Rij, Belnap & Ligrani 2002), effective slope (Napoli, Armenio & De Marchis 2008), roughness skewness (Flack & Schultz 2010), roughness directionality (Nugroho, Hutchins & Monty 2013), spanwise heterogeneity (Anderson et al. 2015) among other parameters.

There exists an extensive body of research dedicated to the study of flows over idealised two- and three-dimensional topographies such as cubes, bars, pyramids or spheres (see e.g. Cheng & Castro 2002; Schultz & Flack 2009; Volino, Schultz & Flack 2009). However, these surfaces have typically a limited range of scales through which they can interact with the turbulent flows. Fractal-like geometries on the other hand, have an inherent capacity to interact with turbulent flows through a much broader range of length scales. These flows have received special attention in the previous two decades from both the numerical and experimental communities, such as in flows encountering fractal grids or trees (see e.g. Meneveau & Katz 2000; Mazzi & Vassilicos 2004; Chester, Meneveau & Parlange 2007; Hurst & Vassilicos 2007; Gomes-Fernandes, Ganapathisubramani & Vassilicos 2012; Graham & Meneveau 2012). Particularly, realistic surfaces such as urban areas and natural landscapes are generally multiscale, exhibiting a fractal-like behaviour with a height spectral slope ranging between $-$1 and $-$3 (Huang & Turcotte 1989; Passalacqua et al. 2006; Wan & Porté-Agel 2011). These have only until the last decade started earning a wider attention from researchers, especially the modelling community. The readers are encouraged to refer to the recent study by Zhu & Anderson (2018) which offers an excellent and comprehensive review of the different contributions in the literature. One of the challenges when examining these flows is the determination of the smallest spatial resolution needed in order to accurately predict the bulk-flow quantities over surfaces featuring a multitude of scales. In fact, unless direct numerical simulations (DNS) or well-controlled laboratory measurements are accessible, computational fluid dynamics (CFD) methods such as large-eddy simulations (LES) or Reynolds-averaged Naviers–Stokes (RANS) require modelling of the unresolved roughness features.

Both ‘discrete’ and ‘continuous’ types of rough surfaces have recently been documented within the multiscale roughness framework. Using LES data obtained from an ABL flow over various irregular topographies, Wan & Porté-Agel (2011) examined the effect of the subgrid-scale roughness on the effective aerodynamic roughness length scale, and found the latter quantity to linearly increase with the root-mean-square (r.m.s.) of the unresolved roughness topography. In an LES study of a channel flow over multiscale fractal-like surfaces, Anderson & Meneveau (2011) developed a dynamic roughness model whereby the unresolved small-scale topography fluctuations are modelled using a local equilibrium wall model. The effective roughness length scale is parameterised through the r.m.s. height and a roughness parameter determined a posteriori and validated with an invariance resolution test. They showed that, by including more and more subgrid-scale roughness, the wall shear stress increased, indicating the importance of correctly modelling the unresolved roughness fluctuations.

In an effort to examine the effect of the intermediate and small topographical length scales, Mejia-Alvarez & Christensen (2010, 2013) experimentally studied a turbulent boundary-layer flow over a highly irregular surface roughness replicated from a turbine-blade damaged by deposition of foreign materials, from which they created filtered models. By comparing the original surface to its lower order representations, they observed the mean and turbulent statistics in the outer region to be self-similar as they are predominately governed by the large topographical scales. However, they reported that the filtered models failed to appropriately reproduce the contributions of the most intense Reynolds shear stress events (sweeps and ejections) within the roughness sublayer (RSL). In a similar approach to Anderson & Meneveau (2011), Barros, Schultz & Flack (2018) investigated a channel flow over systematically generated multiscale rough surfaces with varying spectral slopes, with predefined statistical moments. They showed that the surfaces with shallower spectral slopes led to higher drag, suggesting that high wavenumbers (small scales) can noticeably contribute to the overall drag. Stewart et al. (2019) have examined three different self-affine bed roughness in two different open-channel flows, and reported similar findings to Barros et al. (2018), such that decreasing the roughness spectral exponent leads to higher drag. They explained this behaviour by analytically demonstrating the existence of a relation between the effective slope with the roughness spectral exponent (also validated with experimental data), which increases as the spectral exponent decreases. They suggested that this drag increase is due to contributions from an increase in the flow separation around steeper, scaling-range roughness features.

For the discrete roughness type, Yang & Meneveau (2017) used LES to examine turbulent boundary layers over randomly distributed prism-like roughness elements with a fractal-like size distribution, and investigated the flow response to the addition of smaller roughness generations. To account for the drag that rises from the unresolved scales, they developed a dynamic-roughness model which assumes the flow to become scale-invariant when it interacts with the scale-invariant roughness, which was then validated by the independence of the mean-flow profiles on grid resolution. They also proposed an analytical model that explicitly accounts for the sheltering effect within the canopy flow; however, its predictions seem to depend on the number of generations included in the roughness. In a similar approach to Yang & Meneveau (2017), Zhu & Anderson (2018) performed a series of LES simulations of a turbulent channel flow over cube-like fractal topographies, with varying fractal dimensions (spectral slope). They assessed the drag penalty associated with changing the spectral slope and the number of smaller generations. They showed that the drag associated with the unresolved generations can be parametrised with the first few generations, and proposed a logarithmic law-based model to model their contributions. They further added that the turbulence statistics are predominantly controlled by the first generation of elements.

Despite the burgeoning interest in this area, there is still a need for experimental studies that relate to flows over such complex surfaces. To our knowledge, the systematic investigation of the multiscale roughness hierarchy effect on turbulent flows in regular roughness has only been reported in a limited number of works (e.g. Yang & Meneveau 2017; Zhu & Anderson 2018; Vanderwel & Ganapathisubramani 2019), and remain experimentally largely unexplored. It is also worth mentioning that when simulating turbulent flows over these types of surfaces, the previous LES studies usually rely on channel flow conditions, in which case the conclusions may not be fully applicable, for instance, to external flow systems such as an ABL flow. Therefore, wind-tunnel experiments in addition to the modelling techniques can prove beneficial with the different lines of enquiries mentioned above. In fact, understanding the dynamics of turbulent flows over multiscale rough surfaces will result in appropriate modelling strategies and better predictions of flows such as in urban environments.

To this end, a comprehensive experimental study is designed aimed at examining the effect of roughness-scale hierarchy on turbulent boundary layers. The design used is similar to the one employed by Yang & Meneveau (2017) and Zhu & Anderson (2018) which first builds large-scale cuboids, subsequently adding smaller iterations where their size decreases with a power-law as the number increases (see e.g. the study of Zhu & Anderson (2018) for a detailed description of an iterated function system). The main difference between the current investigation and the former studies lies in the layout of the roughness in the horizontal plane, currently chosen to be uniform instead of randomly distributed. To uncover the effects of the multiscale roughness on the turbulent boundary layer, direct wall drag as well as flow field measurements are conducted using an in-house floating-element drag balance (Ferreira, Rodriguez-Lopez & Ganapathisubramani 2018) and particle image velocimetry (PIV). The remainder of the manuscript is presented in three sections. The experimental methodology is described in § 2 depicting the multiscale roughness topography, drag and flow-field measurement techniques. The results and discussion are presented in § 3, analysing the near-wall and outer-flow regions. Finally, a summary and concluding remarks are provided in § 4.

2. Experimental methodology

2.1. Facility

Measurements were carried out in an open-circuit suction-type wind tunnel at the University of Southampton. The working section follows a 7:1 contraction and extends 4.5 m in the streamwise direction, with a cross-section of $0.6\ \textrm {m} \times 0.9\ \textrm {m} \times 4.5\ \textrm {m}$ in wall-normal and spanwise directions, respectively. The test section was designed with a weak diverging cross-section to allow a constant free stream along the streamwise direction and a growth of a turbulent boundary layer with a nominally zero pressure gradient. The acceleration parameter $K = (({\nu }/{U_{\infty }}) ({\textrm {d}U_{\infty }}/{\textrm {d}x}))$, where $\nu$ and $U_{\infty }$, the air kinematic viscosity and free-stream velocity, respectively, have been observed in the range $1\text {--}5\times 10^{-8}$ from various recent studies performed in this facility (see e.g. Ferreira et al. 2018; Medjnoun, Vanderwel & Ganapathisubramani 2018). The turbulent boundary layer grows over a flat surface composed of five equally sized wooden boards onto which the roughness was mounted. In order to assess the skin friction, the boundary-layer plate at the measurement location was cut out to allow the floating-element drag balance to be inserted in. The boundary-layer plates were preceded by a ramp of 0.2 m long inclined by four degrees to the horizontal, ensuring a smooth transition of the flow from the bottom floor of the test section. The free-stream velocity can reach up to 30 ms$^{-1}$, with a turbulence level less than 0.5 $\%$, and was monitored and acquired using the Micromanometer FCO510 by Furness Controls. To account for air-density variations, temperature was also acquired, and its standard deviation for an average run was less than ${\pm }0.5\,^{\circ }\textrm {C}$. The facility schematic as well as the experimental procedures employed in this study are shown in figure 1.

Figure 1. ($a$) Schematics of the experimental arrangements including the planar- and stereo-PIV set-ups along with ($b$) a close-up on the floating-element drag balance developed by Ferreira et al. (2018). ($c$) Example of instantaneous planar and stereoscopic velocity fields performed at two streamwise and spanwise symmetry planes. Both PIV and drag-balance measurements are performed at a fetch distance $x \approx 30\delta$ from the leading edge, with $\delta$ being the boundary-layer thickness at this fetch determined from the horizontally averaged velocity profiles.

2.2. Multiscale roughness

The surface roughness considered in the current investigation is based on a multiscale distribution pattern as illustrated in figure 2. The design is inspired by a Sierpinski carpet model, which essentially builds the roughness by superimposing size-decreasing self-similar cuboid elements. The model offers an easy control of the geometrical characteristics of the roughness, allowing numerical simulations to be replicated in wind tunnels and vice versa. To generate the fractal-like roughness, a similar method employed by Zhu & Anderson (2018) is used, which considers an iterated function system to construct the successive generations with $h^{(i+1)} = rh^{(i)}$, with $h^{(i)}$ and $r$ being the cuboid height at the $i$th iteration and the scale-reduction factor, respectively. The cuboid-side width is defined as $w^{(i)} = A_{r} h^{(i)}$, with $A_{r}$ being the width-to-height aspect ratio. These roughness elements are uniformly distributed in the horizontal $(x,z)$-plane with each individual cuboid of a given iteration being equidistant from its neighbouring ones, that is $S_x^{(i)} = S_z^{(i)}=S^{(i)}$, with $S^{(i)}$ being the streamwise/spanwise spacing between two successive elements at the ith iteration. For the first iteration, $S^{(1)} = 2 w^{(1)}$, whereas for the following iterations, $S^{(i)} = w^{(i)}$. It should be noted that for scaling purposes, $S^{(1)}$ will exclusively be referred to as $S$ in the remainder of the manuscript. The subsequent self-similar cuboids following the first iteration are overlaid both at the surrounding as well as on top of the previous parent cuboid.

Figure 2. Schematics of the different multiscale roughness combinations tested in the present investigation labelled Iter$_{1}$, Iter$_{12}$, Iter$_{13}$ and Iter$_{123}$.

Due to both roughness manufacturing constraints and flow considerations, the height of the first iteration roughness element was set to $h^{(1)} = 8$ mm, so that $\delta /h^{(1)} \ge 10$. This constrained the number of iterations to three due to the lowest cuboid height that can be manufactured. The scale-reduction factor $r$ and the width-to-height aspect ratio $A_{r}$ were fixed to 1/4 and 4.2, respectively, whereas the number of elements was set to 1, 36 and 576 for the three successive iterations, namely, large, intermediate and small cuboids, respectively. The roughness elements were confined in a $100\ \textrm {mm}\times 100\ \textrm {mm}$ repeating unit tile, which was 3-D printed then replicated using a moulding-casting manufacturing process to cover the floor of the wind-tunnel test section.

In order to investigate the effect of the multiscale roughness on the turbulent boundary-layer flow and assess the contributions of different length scales, four combinations are considered. (i) Iter$_{1}$: which considers the presence of a single iteration with a cuboid height $h^{(1)}$. (ii) Iter$_{12}$: which looks at the presence of two iterations with cuboid heights $h^{(1)}$ and $h^{(2)}$. (iii) Iter$_{13}$: which examines the presence of two iterations with cuboid heights $h^{(1)}$ and $h^{(3)}$. Finally, (iv) Iter$_{123}$: which investigates the effect of the presence of all three iterations with cuboid heights $h^{(1)}$, $h^{(2)}$ and $h^{(3)}$. Schematics for the different cases are illustrated in figure 2 while a summary of their different statistical attributes are presented in table 1.

Table 1. Essential geometrical characteristics of the different multiscale rough surfaces with $h$, $\bar {h}$, $h_{rms}$, $h_{skw}$ being the maximum, mean, root-mean-square and skewness of the roughness height, while $\lambda _{f}$, $\lambda _{p}$ represent the frontal and plan solidities within one repeating unit. The parameter $\alpha = ({\bar {h}}/{h_{rms}}) \lambda _{f}$ is a geometrical parameter used to describe the aerodynamic roughness length scale of the surfaces discussed in § 3.1.3.

2.3. Floating element

The wall shear stress was directly measured by means of a floating-element drag balance, developed and validated by Ferreira et al. (2018) in both smooth- and rough-wall flow conditions. This method relies on the determination of the streamwise net force $F_{d}$ acting on a finite and structurally independent area $A$, representative of the surface being investigated. It makes use of the parallel-shift linkage principle, through pairs of bending beam transducers which allow the isolation of both the streamwise load as well as the induced pitching moment. The wall shear stress is then deduced via $\tau _{w} = F_{d}/A$. The floating element used consists of a $200 \times 200$ mm$^{2}$ surface area, on top of which four roughness tiles are placed. The balance was flush-mounted with the wind-tunnel floor and positioned around 3.1 m downstream of the leading edge. Additional care was taken when setting up the surrounding roughness owing to the tight tolerance of the air gap which is only 0.5 mm wide.

To assess whether these surfaces can reach the fully rough regime, the floating element was exposed to a series of nine free-stream speeds ranging from 9 up to 27 ms$^{-1}$. Each acquisition lasted 30 s with a sampling rate of 150 Hz (corresponding to 2500 eddy turnover times at the lowest speed), with a total of five repetitions per speed. Pre- and post-calibrations were performed for each configuration without notable discrepancies. A more detailed discussion on the complete design, acquisition procedure and uncertainty of the measurement technique can be found in Ferreira et al. (2018).

2.4. Particle image velocimetry

The turbulent boundary-layer flow was diagnosed in both the streamwise-wall-normal plane $(x,y)$ as well as in the cross-plane $(\kern-0.05pt y,z)$, using planar two-dimensional two-component (2D2C) and stereoscopic two-dimensional three-component (2D3C) PIV measurements, respectively. Both types of PIV measurements were performed at approximately 3.2 m downstream of the contraction. As illustrated in figure 1($c$), two 2D2C-$(x,y)$ measurements were conducted at the spanwise symmetry planes $z/S = 0,0.5$ (at the roughness valley and ridge), and two additional 2D3C-$(\kern-0.05pt y,z)$ measurements were also acquired at two successive streamwise symmetry planes, at ${x= 3.2}$ and 3.25 m. The flow was seeded with vaporised glycerol-water particles generated by a Magnum 1200 fog machine, then illuminated with a laser light sheet sourced by a two-pulse Litron Nd:YAG laser operating at 200 mJ. A LaVision optical system for the beam focus/expansion of the light sheet was used, which comprised convex and concave lenses in order to focus the beam, and a cylindrical lens in order to expand the sheet with relatively constant thickness in the measurement plane ($\approx$1 mm thickness).

The particle images were recorded by high-resolution LaVision Imager LX 16 MP CCD cameras fitted with 200 and 300 mm AF Micro Nikon lenses for the 2D2C- and 2D3C-PIV measurements, respectively. One camera was used for the planar-PIV set-up, and positioned at nearly 1 m away from the object plane, whereas two cameras mounted on Scheimpflug adapters to account for the oblique view angle ($\pm 42^{\circ }$) were needed for the stereo-PIV measurements, and placed at nearly 1.3 m from the object plane at either sides of the test section, also depicted in figure 1($a$). A single and a double-sided dual plane calibration target aligned with the laser light sheet were used to determine the mapping function for each set-up, using a third-order polynomial fit. This resulted in a field of view of approximately $0.9\delta \times 1.4\delta$ in the $(x,y)$-plane and $1.3\delta \times \delta$ in the $(\kern-0.05pt y,z)$-plane for the planar and stereoscopic PIV, respectively. In terms of the largest roughness length scale $S$, this was equivalent to $1.1S \times 1.4S$ in the $(x,y)$-plane, whereas in the $(\kern-0.05pt y,z)$-plane it spanned $1.5S \times 1.1S$, for the planar and stereoscopic PIV, respectively.

A total of 3000 statistically independent realisations of image pairs were acquired for each case at 0.6 Hz, with a time delay between pulses of 50 and 20 $\mathrm {\mu }$s for the 2D2C and 2D3C PIV, respectively. In order to allow a comparison between the different surfaces, the PIV measurements were performed at different speeds to obtain a relatively matched roughness function $\Delta U^{+}$, and ensure they reached a fully rough regime (although it will be seen later that Iter$_{1}$ never reached the fully rough condition). The free-stream speeds were adjusted from case to case, with $U_{\infty } = 18.7$, 10.3, 18.5 and 10.2 ms$^{-1}$ for the cases Iter$_{1}$, Iter$_{12}$,Iter$_{13}$ and Iter$_{123}$, respectively. The velocity vector fields were then obtained by interrogating particle images using a decreasing multipass scheme starting from $48\ \textrm {pixels} \times 48\ \textrm {pixels}$ down to a final pass of $16\ \textrm {pixels} \times 16\ \textrm {pixels}$ for the 2D2C and $24\ \textrm {pixels} \times 24\ \textrm {pixels}$ for the 2D3C. Using a 50 $\%$ interrogation window overlap, the resulting effective vector spacing ranged between 0.4 and 0.5 mm for the different cases both in the planar and stereoscopic-PIV measurements.

3. Results and discussion

This section focuses on the analysis of results reported from the floating-element drag balance with the data from both planar- and stereoscopic-PIV measurements. Section 3.1 discusses the results of the direct-drag estimates along with the flow topology in the canopy region, with the assessment of the aerodynamic parameters associated with themean-flow profiles. Section 3.2 reveals the multiscale-roughness effects on the flow topology in cross-plane, focusing on the turbulence properties in the outer region as well as the assessment of the outer-layer similarity hypothesis. The results presented in this study are openly accessible on the roughness database http://roughnessdatabase.org/ and the University of Southampton repository at: https://doi.org/10.5258/SOTON/D1765.

3.1. Inner region

3.1.1. Surface drag

Results from the floating-element drag balance are presented in figure 3. They describe the response of the wall shear stress to the four multiscale rough surfaces investigated. Figure 3($a$) examines how the skin-friction coefficient $C_{f}$ varies with respect to $Re_{x}$ (with varying free-stream speed) in a linear-logarithmic scale. At the lowest tested Reynolds numbers ($Re_{x} <3 \times 10^{6}$), the roughest surface Iter$_{123}$ is distinctively shown to experience the largest frictional drag, and is noticeably higher than Iter$_{12}$. On the other hand, the surfaces Iter$_{1}$ and Iter$_{13}$ exhibit relatively similar values and are the weakest among the four cases. Beyond a certain Reynolds number ($Re_{x} >3 \times 10^{6}$), the skin-friction coefficients are seen to weakly vary for the cases Iter$_{12}$, Iter$_{13}$ and Iter$_{123}$ with an average variation about the mean less than $\pm$2 $\%$, which is within the uncertainty of the measurements. Conversely, a substantial decay of $C_{f}$ is observed for the Iter$_{1}$ case, perhaps with a steeper slope in comparison with the smooth-wall curve represented by Schlichting's power-law $C_{f}=b\,Re_{x}^{a}$ (Schlichting 1979), highlighted with the black solid line.

Figure 3. ($a$) Skin-friction coefficient estimates for one repeating unit at various moderate $Re_{x}$, with the black solid line representing Schlichting's power-law for smooth-wall with $a \approx -1/5$ and $b \approx 0.058$, and the blue stars representing the smooth-wall data from Medjnoun et al. (2018). The red stars highlight the skin-friction estimates corresponding to the flow conditions where PIV measurements are performed. ($b$) Variation of the skin-friction coefficient and ($c$) relative drag increase for the different multiscale rough surfaces with respect to the non-dimensional roughness parameter $\alpha$. The data points in panels (b,c) represent values averaged over the different Reynolds numbers.

It is also observed that at the highest Reynolds numbers of the facility, Iter$_{12}$ experiences a reduction in $C_{f}$ of approximately 4 $\%$ from its mean value across Reynolds numbers. These two different behaviours demonstrate that Iter$_{1}$ and maybe Iter$_{12}$ still remain transitionally rough, whereas the cases Iter$_{13}$ and Iter$_{123}$ have reached the fully rough regime, i.e. Reynolds number invariance of $C_{f}$ (Jiménez 2004; Flack & Schultz 2014). As demonstrated by other studies (see e.g. Napoli et al. 2008; Yuan & Piomelli 2014), the fully rough regime stems from pressure-drag contributions prevailing compared with the viscous drag. This suggests that the viscous-drag contributions still represent a considerable fraction in the Iter$_{1}$ case (at least within the range of $Re_{x}$ tested herein), as opposed to the other cases, whose pressure-drag contributions dominate owing to the presence of additional roughness scales. Although this result seems counterintuitive since the large cuboid is essentially a sharp-edged bluff body, it should be recalled that the frontal and plan solidities of Iter$_{1}$ are noticeably smaller in comparison with the other cases. This is further investigated in § 3.1.2.

The magnitude of the skin-friction coefficient is also considerably affected by the roughness content of the surface, such that the addition of subsequent smaller scales enhances the shear stress acting at the wall. This is clear when examining the skin-friction coefficient averaged over the different Reynolds numbers $C_{F}$ as shown in figure 3($b$). It indicates that the magnitude of the skin-friction coefficient $C_{F}$ increases monotonically with increasing $\alpha$ (with $\alpha = ({\bar {h}}/{h_{rms}}) \lambda _{f}$ being a non-dimensional roughness parameter proportional to the frontal solidity). There are many alternatives and possible ways that can be used to correlate the drag with roughness statistics. In this instance, the use of this surface parameter is justified by the fact that its variation approaches a linear trend with the drag, while other parameters such as $\lambda _{p}$ would fail to capture a clear trend.

To quantify the difference between these surfaces, a relative drag-increase coefficient $\beta$ (where $\beta = 100\,\% \times C_{F}^{(i)}/C_{F}^{(4)}$ with $i=1\text {--}4$ for Iter$_{1}$, Iter$_{12}$, Iter$_{13}$ and Iter$_{123}$) is shown in figure 3($c$) for the different cases. We observe that the largest cuboid alone is responsible for over 80 $\%$ of the drag of the full surface. It is also shown that the drag increase is higher for Iter$_{12}$ than Iter$_{13}$ (93 $\%$ against 87 $\%$ of drag of the full surface), indicating that the drag increment that stems from the intermediate iteration is more prominent than the one caused by the smaller iteration despite Iter$_{12}$ and Iter$_{13}$ having a similar plan solidity $\lambda _{p}$. Therefore, it can be implied that the contributions to the full multiscale surface from the smallest cuboid roughness alone amount to $\beta ^{(3)} = \beta ^{(123)}-\beta ^{(12)} \approx 7\,\%$. Similarly, the contributions from the intermediate roughness scale alone amount to $\beta ^{(2)} = \beta ^{(123)}-\beta ^{(13)} \approx 12\,\%$. These results are in agreement with the observations reported in the previous LES studies, highlighting the importance of small roughness features and emphasising the attention needed when modelling the unresolved subgrid-scale roughness in numerical simulations (Yang & Meneveau 2017; Zhu & Anderson 2018).

3.1.2. Flow topology in the canopy

The effect of the multiscale rough surfaces on the flow topology is explored by examining the normalised mean streamwise velocity maps from the planar-PIV measurements, at the peak symmetry plane ($z/S = 0.5$). Two cases, Iter$_{1}$ and Iter$_{123}$, are illustrated for comparison in figure 4. Results show that, away from the wall, the flow is unaffected by the surface condition, while in the roughness-affected layer it undergoes strong changes in the streamwise direction (as shown by the streamline contours). Specifically within the surface canopy, the maps highlight the presence of a separation bubble past the large cuboid, whose size and intensity appear to depend on the surface condition. This is illustrated in figure 4($b$) which shows that by adding the intermediate and smaller scales, Iter$_{123}$ produces a weaker streamwise velocity deficit within the canopy.

Figure 4. Contour maps of the ($x,y$)-plane normalised mean streamwise velocity for ($a$) Iter$_{1}$ and ($b$) Iter$_{123}$ measured at the peak symmetry plane $z/S = 0.5$. The cross-sections of the roughness geometries are included at the bottom of the figures to scale for reference. The brown solid line represents the zero-velocity contour level illustrating the separation length, while the mean in-plane streamlines are superimposed to highlight the recirculation region. ($c$) Contours of the separation lengths for the different multiscale rough surfaces normalised by the first iteration cuboid height $h^{(1)}$, with $\hat {x}$ being the streamwise distance from the leeward side of the cuboid.

The impact of the multiscale roughness on the recirculation region is further assessed by examining the separation length as shown in brown in figures 4($a$) and 4($b$), and highlighted in an appropriate scaling in figure 4($c$). The results show the extent of the zero contour-level curves to be relatively conserved for the different cases, irrespective of the surface condition, and is approximately $3h^{(1)}$. It is, however, reported that the wall-normal extent subtly weakens with addition of iterations. This behaviour is believed to be caused by the increased turbulent mixing and wall drag at the top of the large cuboid, owing to the interaction of a broader range of roughness scales with the separating shear layer. To explore further the effect of the multiscale roughness on turbulence, the normalised Reynolds shear stress $-\overline {uv}^{+}$ and its associated turbulence production $P_{xy} \delta /U_{\tau }^{3}$ (with $P_{xy} = -\overline {uv} \,\textrm {d}U/{\textrm {d} y}$) fields are examined, as shown in panels (a–d) and panels (e–h) of figure 5, respectively. For Iter$_{1}$, the shear layer formed at the top surface of the large cuboid separates at its trailing edge, from which strong vortices are shed towards the canopy. By comparing the area within the encapsulated contour level in panels (a–d) of figure 5, we observe that the separated shear layer carries a weaker vortical activity as more roughness length scales are imposed. However, a plausible explanation for this behaviour can be related to the range over which turbulent and roughness length scales interact with each other. In fact, when a turbulent flow encounters a broader range of scale, such as in the case of Iter$_{123}$, the $h^{(1)}$-scale vortical structures as seen for Iter$_{1}$ break down to smaller eddies. This means that while the overall drag increases when increasing the roughness content, the shear stress activity is redistributed among scales whereby the intermediate and smaller scales are responsible for a larger portion of drag. Consequently, this results in a relatively shorter extent and a weaker separation bubble and more drag emanating from these additional scales. These observations are further substantiated by the turbulence production maps which show stronger shear layers to be associated with substantial turbulence production. In contrast, weaker shear layers are accompanied with weaker turbulence production.

Figure 5. Contour maps of the ($x,y$)-plane normalised panels (a–d) Reynolds shear stress and panels (e–h) turbulence production for $(a{,}e)$ Iter$_{1}$, $(b{,}f)$ Iter$_{12}$, $(c{,}g)$ Iter$_{13}$ and $(d{,}h)$ Iter$_{123}$ measured at the peak symmetry plane $z/S = 0.5$. The cross-sections of the roughness geometries are included at the bottom of the figures to scale for reference.

Using the streamwise-wall-normal PIV measurements, the mean pressure distributions are additionally estimated at the peak symmetry plane ($z/S$ = 0.5), by means of the two-dimensional Reynolds-averaged Naviers–Stokes equations. For more details on the methodology as well as the numerical integration schemes employed, the reader should refer to the study by Ferreira & Ganapathisubramani (2020). The mean pressure is expressed by its non-dimensional form as $C_{p} = (P-P_{\infty })/q_{\infty }$, with $P$ and $P_{\infty }$ being the mean static and free-stream pressure, while $q_{\infty }$ is the free-stream dynamic pressure.

An example of a pressure field is shown in figure 6($a$) for the Iter$_{1}$ case. The coefficient of pressure field is shown to be dominated by an alternating high- and low-pressure region in the canopy, whereas it quickly recovers to the free-stream pressure in the outer region ($y \ge 0.2\delta$). In fact, the high pressure is seen to be associated with accelerating flow regions while low pressures are accompanied with decelerating flow regions. More specifically, the highest pressure regions are recorded at the windward side of the cuboid, while the lowest pressure magnitude is shown to occur right past the leading edge of the upper surface of the cuboid, reflecting the presence of a strong shear layer shedding off from the blunt leading edge. This negative pressure is shown to trail downstream, forming a weak patch of low pressure within the recirculation region until the reattachment point. The pressure transitions to a positive magnitude as the flow reaccelerates again past the reattachment point.

Figure 6. ($a$) Normalised mean pressure field reconstructed from the 2D2C-PIV for Iter$_{1}$ case measured at the peak symmetry plane $z/S = 0.5$. The cross-sections of the roughness geometries are included at the bottom of the figures to scale for reference. The brown solid line represents the zero-velocity contour level illustrating the separation length, while the mean in-plane streamlines are superimposed to highlight the recirculation region. ($b$) Mean streamwise pressure difference profile assessed from a single cuboid element of height $h^{(1)}$ for Iter$_{1}$.

The profile of the streamwise pressure difference across the largest cuboid is additionally examined and shown in figure 6($b$). The pressure difference is expressed as $\Delta P (\kern-0.05pt y) = P_{w}(\kern-0.05pt y)-P_{l}(\kern-0.05pt y)$ (with the $w$ and $l$ subscripts referring to the windward and leeward sides of the cuboid, respectively), normalised by the wall-normal-averaged pressure $\langle \Delta P \rangle$, and plotted against the wall-normal distance normalised by the cuboid height. Although the profile is not fully resolved down to the wall (due to the spurious region below 1 mm height), the drag profile shows similarities with the results of urban-like roughness (Ferreira & Ganapathisubramani 2020), with a maximum at nearly two-thirds of the cuboid height, and decreasing when getting close to the wall. The profile is subsequently used to get an estimate of the pressure drag produced by the large cuboid. This is done by integrating the streamwise pressure-difference profile over the cross-sectional area $h^{(1)} \times W^{(1)}$ (assuming the pressure distribution around the cuboid is uniform in the spanwise direction). The drag force and the friction velocity over the large cuboid are expressed as

(3.1a,b)\begin{equation} F^{(1)} = W^{(1)} \int_{0}^{h^{(1)}} \Delta P \,\textrm{d} y, \quad U_{\tau}^{(1)} = \left( \frac{F^{(1)}}{\rho h^{(1)} W^{(1)}} \right)^{1/2}. \end{equation}

with $\rho$ being the air density. The results reveal that the ratio between the friction velocity estimated from the largest cuboid to the total friction velocity estimated with the drag balance is approximately $60\,\%$ for Iter$_{1}$. Assuming this ratio to remain constant throughout the different cases, the pressure-drag contributions from the addition of intermediate and small roughness scales ultimately leads to the fully rough regime. In fact, for the other three surfaces, the intermediate and small roughness length scales proportionally add form-drag contributions to the overall drag, resulting in the Reynolds number invariance of $C_{f}$ observed in figure 3. In contrast, the pressure drag contribution for Iter$_{1}$ that stems from the largest cuboid alone is insufficient to result in a flow that is fully rough. Due to the viscous drag which amounts to nearly $40\,\%$ of the total drag, the skin-friction coefficient $C_{f}$ remains dependant of $Re_{x}$ as seen in figure 3. It should be noted that despite the inherent degree of uncertainty that rises from the PIV-based pressure and the assumption of the spanwise uniformity of the profiles, this method gives a reasonable indication of the expected form drag.

3.1.3. Aerodynamic parameters

In order to estimate the aerodynamic quantities that characterise the rough-wall flow, both the peak and valley symmetry planes ($z/S$ = 0, 0.5) PIV data are used. Figure 7($a$) compares the wall-normal distribution of the mean streamwise velocity at the symmetry planes for Iter$_{1}$. The angled brackets $\langle \,\rangle$ in figure 7($b$) refer to the horizontally averaged velocity between both planes, comparing profiles of the different cases with the smooth-wall. Substantial differences between the peak and valley profiles can be observed from figure 7($a$). Below the canopy layer, the velocity is shown to be higher at the valley than at the peak as expected due to the blockage of the large obstacles. At $y/\delta \approx 0.15$, it is observed that $U_{Peak}=U_{Valley}$, with the subscripts ‘peak’ and ‘valley’ referring to the velocity measured at $z/S= \pm 0.5$ and $z/S = 0$, respectively. However, beyond this point until almost two thirds of the boundary-layer thickness, the velocity above the peak becomes higher than that at the valley. This indicates that the outer flow presents a degree of spanwise heterogeneity, probably caused by the surface condition. It is also observed that at the individual planes, the velocity profiles do not exhibit a clear wall-normal logarithmic distribution. However, when both planes are averaged, a velocity range appears to vary logarithmically, and occurs approximately between $0.2$ and $0.3\delta$. It is further shown from the profiles of the different surfaces plotted in figure 7($b$), that all the cases deviate from the smooth-wall behaviour. More specifically, the richer the roughness content is, the higher the profile deficit becomes.

Figure 7. ($a$) Wall-normal distribution of the mean streamwise velocity at both symmetry planes for Iter$_{1}$. ($b$) Comparison of the horizontally averaged mean streamwise velocity profiles between the different rough surfaces.

The mean streamwise velocity profile over the rough surfaces can also be expressed using the modified law-of-the-wall,

(3.2)\begin{equation} \langle U^{+} \rangle \equiv \frac{\langle U \rangle}{U_{\tau}}=\frac{1}{\kappa} \ln \left[(\kern-0.05pt y-d)\frac{U_{\tau}}{\nu}\right]+B-\Delta U^{+}\left(\frac{h_{s}U_{\tau}}{\nu}\right), \end{equation}

where $\kappa$ and $B$ represent the slope and intercept of the logarithmic region, respectively, similar to a smooth-wall. In contrast with a flat smooth wall, rough walls will give rise to the quantities $d$ and $\Delta U^+$ in the form of wall-normal and velocity shifts in the logarithmic region, which are termed the zero-plane displacement and the roughness function, respectively. The former is interpreted as the ‘virtual’ origin representative of the height at which the mean drag acts (Jackson 1981). On the other hand, the latter provides a quantification of the momentum loss (if positive) or gain (if negative) due to surface roughness, and depends on the roughness Reynolds number $h_{s}^{+}$. The wall-normal distribution of the inner-normalised mean streamwise velocity profiles are presented in figure 8, compared with the DNS turbulent boundary-layer profile of Sillero, Jiménez & Moser (2013). As expected, the profiles are shown to be affected by the surface condition as they exhibit both wall-normal as well as velocity shifts in comparison with the smooth-wall. The rough-wall profiles appear to have relatively similar vertical shifts. However, it is worth recalling that the free-stream speeds $U_{\infty }$ were not kept constant among cases.

Figure 8. Inner scaling of the horizontally averaged streamwise velocity profile for the different multiscale rough surfaces, compared with the DNS turbulent boundary-layer profile of Sillero et al. (2013). The value of the log-law slope $\kappa$ and the smooth-wall intercept $B$ used in the current investigation are 0.39 and 4.5, respectively.

The zero-plane displacement $d$ is estimated by making use of the modified law-of-the-wall. This is achieved by taking the derivative of (3.2) with respect to $y$, to obtain what is called the indicator function, expressed as:

(3.3)\begin{equation} \varXi = (\kern-0.05pt y^{+}-d^{+}) \dfrac{\textrm{d} \langle U^{+} \rangle}{\textrm{d} y^{+}} = \dfrac{1}{\kappa}. \end{equation}

This equation is classically used to determine the extent of the inertial sublayer as well as the logarithmic slope $\kappa$ for smooth-wall flows, once $U_{\tau }$ is known a priori (Österlund et al. 2000). This consequently means that if $\kappa$ is assumed to be universal between smooth and rough walls, $d$ should be only a function of the velocity gradient and the friction velocity. Therefore, to avoid solving a two-point fit equation in the present study, the value of $\kappa$ is assumed constant.

Since the friction velocity $U_{\tau }$ was directly measured through the drag balance, the zero-plane displacement is found as the value that minimises the difference between the left-hand side ($\varXi$) and the right-hand side ($1/\kappa$) of (3.3). The results are illustrated in figure 9($a$) and show the appropriate values that minimise the difference yield a good collapse of $\varXi$ for the different cases in the outer region ($y/\delta >0.2$). In the inner region, $\varXi$ is shown to depend on the surface condition with the occurrence of a peak in intensity coinciding with $y=h^{(1)}$. Unexpectedly, the current profiles are observed to result in negative values of $d$ for the cases Iter$_{1}$, Iter$_{12}$ and Iter$_{123}$, with the exception of Iter$_{13}$ which reported a positive value. These are shown to range between $-$0.2 and 0.3$h^{(1)}$. Although the values observed appear to conflict with the interpretation provided by Jackson (1981), namely the height where the mean momentum sink occurs (thus $d$ cannot be negative), these values are merely a solution of the logarithmic distribution fitting. In fact, given that both the velocity and wall shear stress were directly measured, a potential reason for this behaviour can arise from the assumed value of the logarithmic slope $1/\kappa$, hypothesised herein to be universal between smooth and rough walls. However, the discussion on the universality of $\kappa$ and $B$ between smooth and rough walls is beyond the scope of the current study. It is further reported that the indicator function exhibits a plateau $\varXi \approx 1/\kappa$ in the range $0.2\le y/\delta \le 0.3$, suggesting the inertial sublayer to have shifted farther away from the wall, consistently for all the cases. This is in contrast with the smooth-wall boundary layers which are shown to not exceed an ISL upper limit of $0.15\delta$ (equivalent to the $0.15 Re_{\tau }$ with $Re_{\tau }=\delta U_{\tau }/\nu$ reported by Marusic et al. 2013). This reveals that the roughness topography may have severely altered the mean and turbulence-flow structure in the outer region.

Figure 9. ($a$) Variation of the indicator function $\varXi$ and ($b$) the modified log-law function $\varPsi$ for different multiscale surface roughness cases, compared with the smooth-wall profile (blue curve). A second-order central-difference scheme was used in order to obtain $\varXi$.

Once the zero-plane displacement is determined, the difference between the logarithmic velocity distribution and the measured inner-normalised velocity profiles is examined to obtain

(3.4)\begin{equation} \varPsi = \langle U^{+}\rangle -\frac{1}{\kappa} \ln ( y^{+}-d^{+}) -B+\Delta U^{+}. \end{equation}

In the highlighted grey region of figure 9($b$), the profiles of $\varPsi$ reach a plateau at zero for the appropriate values of $\Delta U^{+}$. Similarly to the indicator function, the plateau of $\varPsi$ is clearly identified between $0.2$ and $0.3\delta$, indicating that the inertial sublayer of the multiscale rough surfaces occurs farther than that of a smooth wall. The values of $\varPsi$ are reported to be relatively similar between cases; however, it is recalled that the free-stream speeds were not kept constant. Their values are collated in table 2. It should be mentioned that despite the unexpected negative results of the zero-plane displacement, the values of $d$ have only a marginal influence on the roughness function. In fact, by forcing $d=0$, the values of $\Delta U^{+}$ change by less than 2 $\%$. Furthermore, the profiles exhibit a similar logarithmic slope as the smooth-wall, albeit with a noticeably less intense wake in the outer region. This observation suggests that the multiscale rough surfaces might have substantially affected the wake strength parameter $\varPi$, thereby, the outer-flow dynamics.

Table 2. Aerodynamic parameters of the turbulent boundary-layer flow over the different multiscale rough surfaces. The boundary-layer thickness $\delta$ was identified as the wall-normal distance at which the horizontally-averaged streamwise velocity reached 99 $\%$ of the free-stream speed $U_{\infty }$, while $S = S^{(1)}$. The quantities $\Delta U^{+}$, $h_{s}$ and $\varPi$ namely; roughness function, aerodynamic roughness length scale and wake strength parameter respectively are discussed in § 3.1.3.

Once $\Delta U^{+}$ is determined, the equivalent sandgrain roughness height can be deduced using the fully rough asymptote relation of Nikuradse (1933),

(3.5)\begin{equation} h_{s}^{+} = \textrm{e}^{\kappa(\Delta U^{+}+C)}. \end{equation}

The constant $C$ can be empirically determined by exposing the flow to fully rough conditions, and is defined as the difference between the fully rough intercept $B_{r}$ and the smooth-wall intercept $B$ ($C=B_{r}-B$). In the current study, velocity measurements are only available at one free-stream speed (i.e. one Reynolds number). Hence, the value of $C$ is assumed to be 3.5 which means the flow should have met the fully rough conditions (Schlichting 1979). Following the $C_{f}$ results reported in figure 3, the fully rough condition has been met for the three rough surfaces Iter$_{12}$, Iter$_{12}$ and Iter$_{123}$. At the same time, it is evident from the same plots that Iter$_{1}$ is yet to reach an asymptotic value, hence should not be ascribed a sandgrain roughness value. However, for the purpose of characterising the multiscale rough surface, an $h_{s}$ value will also be assumed for Iter$_{1}$ in order to allow a comparison across all cases.

Equation (3.5) is shown to yield $h_{s}^{+}$ values that range between 400 and 560, which are considerably beyond the $h_{s}^{+} =70$ suggested to be the lower limit for a flow to reach the fully rough regime (Jiménez 2004). Their values are tabulated in table 2 as a fraction of $\delta$, and range between $5$ and $10\,\%$ of the boundary-layer thickness. Figure 10($a$) shows the mean velocity profiles $\langle U^{+} \rangle$ plotted against the normalised wall-normal distance $(\kern-0.05pt y-d)/h_{s}$, and highlights a good degree of collapse with the blue-dashed line which follows a logarithmic distribution. More specifically, the cases Iter$_{12}$ and Iter$_{123}$ are seen to adequately collapse within $2<(\kern-0.05pt y-d)/h_{s}<3$, whereas Iter$_{1}$ and Iter$_{13}$ collapse better within $3<(\kern-0.05pt y-d)/h_{s}<6$.

Figure 10. ($a$) Variation of the inner-normalised velocity profiles with respect to the wall-normal distance normalised by the roughness length scale $h_{s}$ for the different multiscale rough surfaces. The fully rough intercept that best fits the logarithmic region is found to be $B_{r} \approx 8$. ($b$) Variation of the roughness length scale increment $\rho _{h_{s}}$ with respect to the multiscale geometrical parameter $\alpha$.

Two reasons are believed to cause this behaviour: (i) the friction Reynolds numbers for Iter$_{1}$ and Iter$_{13}$ are nearly a factor of two higher than those of Iter$_{12}$ and Iter$_{123}$, leading to a larger scale separation manifested in a slightly wider log region; (ii) the equivalent sandgrain roughness height for Iter$_{1}$ and Iter$_{13}$ is almost half that of Iter$_{12}$ and Iter$_{123}$, causing the profiles to shift away from the wall. The profiles are further shown to have a weaker departure from the logarithmic distribution in the outer region, confirming the influence of the multiscale rough surface on the outer region. Interestingly, all profiles are shown to have a fully rough intercept $B_{r} \approx 8$. It should be noted that the latter is classically reported to be 8.5 (Schlichting 1979). However, this difference is believed to stem from the fact that the 8.5 value was originally derived from pipe-flow data, for which a logarithmic slope is usually taken as $\kappa \approx 0.41$. In contrast, boundary-layer flows have a relatively smaller value of $\kappa \approx 0.39$ which subsequently results in a smaller intercept (Nagib, Chauhan & Monkewitz 2007).

The roughness function (hence the equivalent sandgrain roughness height) is known to depend on various parameters such as the surface texture, shape as well as the plan distribution. There have been numerous proposed correlations between the roughness function and the surface properties (see e.g. Flack & Schultz 2010, 2014). In the current study, we also attempt to find a relation between the aerodynamic quantity $h_{s}$ with an appropriate geometrical parameter of the roughness, such that a link between the different multiscale rough surfaces can be established. Similarly to the observations made from figure 3($c$), the variation of the relative roughness length scale increment $\rho _{h_{s}}$, with respect to the geometrical parameter $\alpha$ is presented in figure 10($b$). In line with the $C_{F}$ results, $\rho _{h_{s}}$ also shows a reasonable linear behaviour observed for the range of surfaces tested herein. This is expressed in a non-dimensional form as

(3.6)\begin{equation} \rho_{h_{s}} \equiv \frac{h_{s}^{(i)}}{h_{s}^{(1)}}-1 = a \alpha^{(i)} + b, \end{equation}

where $h_{s}^{(i)}$ and $\alpha ^{(i)}$ being the equivalent sandgrain roughness height and the geometrical parameter for the $i$th surface (with $i = 1,4$ for Iter$_{1}$ to Iter$_{123}$), respectively. This means that for these regular multiscale rough surfaces, the aerodynamic roughness length scale of the successive multiscale surfaces is linearly related to that of the previous iterations, therefore to the first iteration. This suggests that knowledge of the geometrical properties along with the aerodynamic properties of the parent rough surface are sufficient to predict the drag of a given multiscale rough surface. However, it should be noted that in the general context of multiscale roughness, this linear relation may allow the prediction of $h_{s}$ only for a limited range of scales/iterations such as from Iter$_{1}$ to Iter$_{1234}$ (with a fourth generation of smaller cuboids). In fact, if we keep adding small-scale iterations, the expected drag (hence its associated roughness length scale) that can be generated by the surface will only weakly increase (as if it reaches a certain asymptotic value of $h_{s}$), and can perhaps be more appropriately described by a power-law. In the current study, the limitation in roughness manufacturing as well as the sensitivity of the drag balance to such small magnitudes prevented the assessment of the limits of the observed linear behaviour. It is important to stress that despite the apparent suitability of this relationship for the discrete-like roughness, it currently lacks generalisation over the broader range of rough surfaces (e.g. the highly irregular continuous-type roughness). Nonetheless, while more testing is required to assess its applicability across a wider range, this relationship remains a potential candidate for a subset of topographies such as the multiscale random/regular urban-like rough surfaces.

3.2. Outer region

3.2.1. Secondary motions

Figures 11($a$) and 11($b$) show the cross-stream fields of the normalised mean streamwise velocity for the cases Iter$_{1}$ and Iter$_{123}$, respectively. These maps have been obtained by averaging two cross-planes at the streamwise locations $x= 3.2$ and $x=3.25$ m downstream of the contraction. Unlike the streamwise-wall-normal maps shown in figure 4, the mean cross-plane flow exhibits significant spanwise heterogeneity above the canopy forming alternating high- and low-momentum pathways (HMP and LMP) between the peaks and valleys, respectively. This spanwise undulation in the mean flow is shown to be accompanied by two streamwise rotating cells flanking the sides of the large cuboids emphasised with the in-plane velocity vector plot superimposed on top of the maps. These time-averaged structures were identified using the vorticity-signed swirling strength $\lambda _{ci}$ as shown in figure 12 and indicate the presence of considerable streamwise circulation in nearly half the boundary-layer thickness irrespective of the multiscale rough surface.

Figure 11. Contour maps of the ($y,z$)-plane of the normalised mean streamwise velocity averaged over the two planes at $x=3.2$ and 3.25 m for ($a$) Iter$_{1}$ and ($b$) Iter$_{123}$, respectively. The cross-sections of the roughness geometries are included at the bottom of the figures to scale for reference. The black-dashed line represents the velocity contour level $U=0.8U_{\infty }$. The mean in-plane velocity vector plot $V$ and $W$ are superimposed to highlight the vortical structures.

Figure 12. Effect of the multiscale rough surfaces on the normalised vorticity-signed swirling strength distributions for ($a$) Iter$_{1}$ and ($b$) Iter$_{123}$. The black-dashed line represents the velocity contour level $U=0.8U_{\infty }$. The mean in-plane velocity vector plot $V$ and $W$ are superimposed to highlight the vortical structures on top of the swirling motions.

Figure 12 highlights the rotational sense of the vortical structures, with upwash motions occurring at the valley, while a downwash motion is observed above the peaks. These cross-plane motions are known to be a manifestation of Prandtl's secondary flows of the second kind since they arise from turbulence anisotropy, as opposed to the first kind which stems from mean flow curvature (Prandtl 1952). These features have recently come to scrutiny and are now well-documented for instance in boundary layers (Mejia-Alvarez & Christensen 2013; Nugroho et al. 2013; Willingham et al. 2014; Anderson et al. 2015), channels and open-channel flows (Wang & Cheng 2006; Yang & Anderson 2017; Chung, Monty & Hutchins 2018; Zampiron, Cameron & Nikora 2020) as well as pipe flows (Chan et al. 2018).

A common feature in all these studies is that when the surface spanwise characteristic length scale is comparable to the dominant length scale of the flow (boundary-layer thickness, channel half-height or pipe radius), large-scale secondary motions manifest with generally upwash and downwash motions occurring above ‘low-’ and ‘high-roughness’ respectively. In the present investigation, the largest surface length scale corresponds to the largest phase in the multiscale roughness which is $S$ in both the streamwise and spanwise direction. This incidentally leads to $\delta /S \approx 1.1$ for all different multiscale rough surfaces, which appears to line up with the previous studies, hence, the observed large-scale secondary motions.

The location of the HMP and LMP has recently been questioned since there has been apparent disagreement between some studies that reported opposite trends, where upwash and downwash motions can occur above low- and high-roughness, respectively. In a recent study by Medjnoun, Vanderwel & Ganapathisubramani (2020), it has been highlighted that these conflicting trends are likely to arise from two different driving mechanisms representing non-equivalent types of surface heterogeneity conditions, which can result in opposite observations. It is argued that turbulent secondary flows of the second kind can be generated over two main types of heterogeneous surfaces, namely, topography-driven and skin-friction-driven heterogeneity.

For surfaces purely driven by skin friction heterogeneity (spanwise-uniform height distribution), Willingham et al. (2014) and Chung et al. (2018) have demonstrated that HMP and LMP (hence downwash and upwash) systematically occur above regions of high and low skin friction, respectively. For topographically driven heterogeneity (alternating peaks and valleys), both scenarios have been observed. Several numerical and experimental studies showed upwash and downwash to occur above low- and high-roughness (Mejia-Alvarez & Christensen 2013; Barros & Christensen 2014; Yang & Anderson 2017; Awasthi & Anderson 2018), which is in line with the current observations. Other studies on the other hand highlighted the opposite behaviour (Nezu & Nakagawa 1984; Wang & Cheng 2006; Vanderwel & Ganapathisubramani 2015; Hwang & Lee 2018). However, one of the key differences between these two sets of studies is the presence/absence of streamwise heterogeneity in the topography. For surfaces where HMP and LMP have been reported above peaks and valleys, respectively, the roughness clearly exhibited wake producing protrusions which lead to form-drag contributions, similar to the present study. In contrast, studies that showed upwash and downwash occurring above the elevated and recessed regions, mainly used streamwise homogeneous surfaces. It is worth recalling that the location of the HMP (downwash) and LMP (upwash) is a consequence of the averaging procedure. The instantaneous structures on the other hand have been shown to be able to reverse this behaviour as both scenarios can occur in a non-periodic (chaotic) fashion, as recently highlighted by Anderson (2019).

Hence, despite the current surfaces being topographically spanwise-heterogeneous, the presence of the streamwise heterogeneity leads the surfaces to act more like the skin-friction-type heterogeneous surfaces due to the additional pressure drag, unlike the other set of studies which used surfaces that are predominantly viscous drag producing surfaces. Therefore, care should be taken when designing numerical simulations or laboratory experiments, by accounting for both the roughness height spectral content as well as the phase information, which under certain conditions can dynamically interact with the $\delta$-scale flow structures, leading large modifications in the boundary-layer characteristics.

3.2.2. Turbulent and dispersive stresses

The effect of the multiscale rough surfaces on the turbulence organisation is examined using the triple decomposition performed on the flow field in the cross-plane. This method is usually employed in flows where spatial heterogeneity is prevalent, in order to quantify the magnitude of the stresses originating from the mean flow inhomogeneity. In this scenario, the dispersive stresses can have substantial contributions to the total stresses and play a major role in the transport of momentum fluxes (Meyers, Ganapathisubramani & Cal 2019; Nikora et al. 2019). In this scheme, the streamwise velocity field can be written as

(3.7)\begin{equation} \underbrace{u_{i}(\kern-0.05pt y,z,t)}_{\substack{\textit{Instantaneous}\\\textit{velocity}}} = \overbrace{\underbrace{\langle U(\kern-0.05pt y) \rangle_{S} }_{\substack{\textit{Spatial-time-}\\\textit{averaged}\\\textit{velocity}}} + \underbrace{\tilde{u}(\kern-0.05pt y,z)}_{\substack{\textit{Dispersive}\\\textit{velocity}\\\textit{component}}}}^{\substack{\textit{Time-averaged velocity}}} +\underbrace{u(\kern-0.05pt y,z,t)}_{\substack{\textit{Random}\\\textit{velocity}\\\textit{component}}}, \end{equation}

where $u_{i}$ is the instantaneous velocity field measured at a fixed streamwise location. Here, $\langle U(\kern-0.05pt y) \rangle _{S}$ is the time and horizontally averaged velocity profile over the spanwise wavelength $S$, $\tilde {u}(\kern-0.05pt y,z)$ is the time-invariant spatial deviation field and $u(\kern-0.05pt y,z,t)$ is the time- and space-dependant fluctuating part from the Reynolds double decomposition. Using this method, the different turbulent and dispersive stress tensor terms can be evaluated and compared. More specifically, the total shear stress term can then be computed as

(3.8)\begin{equation} \underbrace{\overline{\tau_{xy}}}_{\substack{\textit{total}\\\textit{shear}\\\textit{stress}}} = \underbrace{\nu\frac{\partial U }{\partial y}}_{\substack{\textit{viscous}\\\textit{shear}\\\textit{stress}}} - \underbrace{\overline{uv} }_{\substack{\textit{turbulent}\\\textit{shear}\\\textit{stress}}} - \underbrace{\overline{\tilde{u} \tilde{v}}}_{\substack{\textit{dispersive}\\\textit{shear}\\\textit{stress}}}, \end{equation}

which typically accounts for viscous, turbulent and dispersive shear stress contributions. The viscous contribution is only present very near the wall, and very negligible away from the wall. On the other hand, Nikora et al. (2019) have shown that the dispersive stresses can be decomposed into roughness-induced and large-scale secondary motion-induced contributions. They have highlighted that the former can have a notable part in the total dispersive stress within the canopy region, whereas the secondary motion-induced part is predominant above the canopy layer. In the present study, we have noticed that these roughness-induced stresses are very negligible throughout the entirety of the resolved flow, as the heterogeneity is mainly driven by the secondary motion-induced stresses.

Contour maps of the turbulent, dispersive and total normal stress terms of Iter$_{1}$ are presented in columns of figure 13, while the shear stress terms are shown in figure 14, respectively. The different normal turbulent stress maps illustrated in figure 13 clearly exhibit spanwise heterogeneity due to the surface condition. In fact, two behaviours can be distinguished, in the near canopy region (at $y/\delta \approx 0.15$), the turbulence intensities are shown to be much higher above the roughness elements while weaker magnitudes are seen between the ridges. In the outer region (at $y/\delta \approx 0.4$), the opposite behaviour occurs with higher intensities occurring between the ridges while weaker magnitudes are located above the roughness elements. These heterogeneous cross-plane turbulence intensity distributions are seen to be related to the rotational sense of the secondary motion, which itself is imposed by the surface condition. The high-turbulence fluctuations observed near the ridges (at $y/\delta \approx 0.15$) are advected laterally towards the mid canopy, and upwards resulting in the LMP. This is further accompanied by the low-turbulence fluctuations being advected from the free stream towards the ridges, inducing the observed HMP above the ridges (clearly emphasised when examining the activity of $\overline {vv}$ distribution). This spanwise turbulence inhomogeneity has previously been shown to be associated with an imbalance between turbulence production and viscous dissipation, causing these amplified lateral motions and sustaining the streamwise vorticity observed in the form of secondary flows (Hinze 1973; Anderson et al. 2015; Hwang & Lee 2018).

Figure 13. Contour maps of the normalised normal components of the streamwise ($a$,d,g), wall-normal ($b$,e,h) and spanwise ($c$,f,i) stresses respectively for the Iter$_{1}$ case. Panels (a–c) represent the turbulent, (d–f) the dispersive and (g–i) the total stresses respectively. The colour scales have been adjusted to highlight the features for the dispersive stresses.

Figure 14. Contour maps of the normalised shear components of the streamwise-wall-normal ($a$,d,g), streamwise-spanwise (b,e,h) and spanwise-wall-normal (c,f,i) stresses, respectively, for the Iter$_{1}$ case. Panels (a–c) represent the turbulent, panels (d–f) the dispersive and panels (g–i) the total stresses, respectively. The colour scales have been adjusted to highlight the features for the dispersive stresses.

Furthermore, the normal dispersive stresses shown in panels (d–f) of figure 13 are shown to be relatively small with respect to the turbulent ones, however, not negligible. They are shown to extend for nearly half of the boundary-layer thickness although with a varying intensity along the spanwise direction. Their magnitude is shown to be relatively weak in comparison with the turbulent fluctuations (shown with a smaller colourmap scale), while their spatial distribution also exhibits spanwise variation different to that observed for the turbulent components. They are in fact shown to be more localised with upwash regions accompanied with intense activity whereas HMP are associated with smaller magnitudes.

The streamwise and wall-normal dispersive stresses are shown to vanish at spanwise locations between HMP and LMP which coincide with weak upwash/downwash motions. On the other hand, the spanwise dispersive stress is shown to be negligible at the LMP and HMP except in between. This region is seen to coincide with enhanced spanwise fluctuations, albeit with a relatively weaker magnitude in comparison with the normal and streamwise dispersive stresses. The total normal stresses expressed as the sum of the turbulent and dispersive stress (with the viscous stress contributions above the canopy being negligible) is illustrated in panels (g–i) of figure 13. The results show that globally, the total stresses remain relatively similar to the turbulent ones except with an increased activity between the large roughness elements corresponding to the upwash regions, which are amplified because of the dispersive stress contributions.

For the shear stress terms presented in figure 14, their spatial distributions are also seen to undergo a spanwise modulation due to the wall condition. The spatial distribution of $\overline {uv}$ presents enhanced shear stress activity above the ridges (in the near canopy) while a weaker activity is shown in the outer region. Between the ridges, the opposite trend happens, with weak activity near the canopy and increased shear stress in the outer layer. At the same time, the dispersive shear stress distribution shows a similar behaviour to that reported for the normal dispersive stress term $\overline {\tilde {u}\tilde {u}}$, with localised patches above the ridges and in valleys, while its magnitude weakens in between. Interestingly, a positive dispersive shear stress activity is reported near the canopy, more specifically above the ridges and slightly less at the valley.

The reason for this behaviour comes from the sign of the product of both quantities $\tilde {u}$ and $\tilde {v}$. These quantities ($\tilde {u}$ and $\tilde {v}$) are shown to have opposite signs at the HMP and LMP, with $\tilde {u} > 0$ and $\tilde {v} < 0$ at the HMP while with $\tilde {u} < 0$ and $\tilde {v} > 0$ at the LMP. However, the latter is true except near the canopy where $\tilde {u}$ exhibits a sign reversal such that $\tilde {u} < 0$ right above the ridge (flow experiences increased local friction) and, at a similar height at the valley $\tilde {u} > 0$ (flow experiences weaker local friction), resulting in the observed $\overline {\tilde {u}\tilde {v}} >0$. This leads to a relative increase of the total shear stress $\overline {\tau _{xy}}$ in the outer region and a reduction near the canopy. The turbulent, dispersive and total shear stress distributions of the streamwise-spanwise and the wall-normal-spanwise shear components are shown in figures 14($b$,e,h) and 14($c$,f,i), respectively. Similarly to the previous stress terms, their spatial distributions are shown to experience spanwise heterogeneity with varying intensity along the span. The sign of the stress maps is shown to be switching at the symmetry planes ($z/S = 0$ and ${\pm }0.5$) which coincide with locations of zero spanwise fluctuations. The extent of the dispersive shear stress terms is seen to reach at least half of the boundary-layer thickness but with a weaker contribution to the total stresses.

3.2.3. Outer-layer similarity

Following the observations made above, outer-layer similarity is also examined. The similarity between the rough- and smooth-wall structural paradigms stems from an established concept of flow universality known as Townsend's wall-similarity, which states that the influence of viscosity is limited to the near-wall region (Townsend 1976). In this analogy, outer-layer similarity extends this concept to include surface roughness with its validity relying on two main conditions: (i) a large scale separation ($Re_{\tau } \gg 1$) and (ii) a small relative roughness height to the depth of the flow ($h/\delta \ll 1$), as emphasised by Jiménez (2004). This is generally expressed in the form of the streamwise velocity deficit and the turbulence quantities which are hypothesised to have a universal form of

(3.9a,b)\begin{equation} \frac{U_{\infty} - \langle U \rangle}{U_{\tau}}= f\left(\frac{y}{\delta}\right), \quad \frac{\langle\overline{u_iu_j} \rangle}{U^2_{\tau}}= g_{ij}\left(\frac{y}{\delta}\right). \end{equation}

This implies that apart from affecting the near-wall region, the main role of roughness is to set the wall drag ($U_{\tau }$) and adjust the boundary-layer thickness ($\delta$), while the mean and turbulence structure remain unaffected by the surface condition. To examine the impact of the multiscale rough surfaces on the outer-layer similarity hypothesis, the horizontally averaged mean streamwise velocity profiles are plotted in defect form and shown in figure 15($a$), while the streamwise turbulence intensity profiles are depicted in figure 15($b$). The mean velocity profiles indicate a clear lack of collapse in the outer region, showing a weaker deficit than that of the smooth-wall. However, it is interesting to point out that between the cases, a reasonable degree of self-similarity is maintained. The absence of similarity is also shown in the inner-normalised turbulence intensity profiles, both with respect to the smooth-wall profile as well as among cases. The profiles are shown to recover similarity only near the edge of the boundary-layer thickness.

Figure 15. Wall-normal distribution of ($a$) the streamwise velocity deficit and ($b$) the streamwise variance profiles for the different multiscale rough surfaces, compared with the DNS turbulent boundary-layer profile.

In order to quantify the changes in the outer flow due to the multiscale rough surfaces, the wake strength parameter $\varPi$ is examined. To this end, it is common to assume the flow in the outer region to be fully described by the composite log-wake profile. This is expressed as

(3.10)\begin{equation} \langle U_{composite}^{+} \rangle= \langle U_{log\text{-}law}^{+} \rangle +\frac{\varPi}{\kappa} \left[2- w \left(\frac{y}{\delta}\right)\right], \end{equation}

where $U_{log\text {-}law}^{+}$ represents (3.2) and $w$ an analytical expression known as Cole's wake function (Coles 1956). There are numerous expressions that have been proposed in the literature, and each one yields a value specific to that function (see e.g. Jones, Marusic & Perry 2001; Junke, Julien & Meroney 2005; Nagib et al. 2007). The wake strength parameter is then determined through fitting the measured velocity profiles to the assumed ‘universal’ expression. In order to circumvent the use of a particular expression for the outer region, the wake strength parameter is instead determined using (3.4), and varies proportionally with the maximum of $\varPsi$. This is expressed as

(3.11)\begin{equation} \varPi = \frac{\kappa}{2} \max (\varPsi) \equiv \frac{\kappa}{2} \max \left[\langle U^{+} \rangle - \left( \frac{1}{\kappa} \ln (\kern-0.05pt y^{+}-d^{+}) + B - \Delta U^{+} \right) \right]. \end{equation}

Their values are summarised in table 2 and are shown to be considerably weaker than that of a smooth-wall (0.57), with $\varPi$ = 0.18, 0.27, 0.20 and 0.23 for the cases Iter$_{1}$, Iter$_{12}$, Iter$_{13}$ and Iter$_{123}$, respectively. These values of $\varPi$ seem to agree with the defect velocity profiles which were shown to have a relatively similar deficit across cases, while being smaller than that of a smooth-wall. The results suggest that universal outer-layer similarity cannot hold for flows developing over this type of surface. In other words, flows that harbour large-scale secondary motions, yielding a highly three-dimensional flow are unlikely to satisfy similarity, as they could have substantially altered the wake intermittency and the uniform momentum zones which together result in the wake strength parameter (Krug, Philip & Marusic 2017).

Using the dispersive stress components, the inner-normalised profiles of the total stress terms $\langle \tau _{xx} \rangle ^{+}$ and $\langle \tau _{zz} \rangle ^{+}$ are examined in figure 16($a$), while $\langle \tau _{yy} \rangle ^{+}$ and $\langle \tau _{xy} \rangle ^{+}$ are shown in figure 16($b$). The spanwise-averaged shear stresses $\langle \tau _{xz} \rangle ^{+}$ and $\langle \tau _{yz} \rangle ^{+}$ are null because of the spanwise antisymmetry at the symmetry plane $z/S =0$, as illustrated in figures 14($b$) and 14($c$). Despite the contributions of the dispersive stresses, overall, the normal and the shear stress profiles show a smaller magnitude when compared with their equivalent smooth-wall. This indicates that the increase in the stresses due to the turbulent and dispersive fluctuations caused by the roughness is not necessarily accompanied with a proportional increase in skin friction, since the profiles do not collapse with the smooth-wall (i.e. lack of outer-layer similarity). A degree of collapse seems to be recovered only beyond $y/\delta >0.6$, which coincides with the wall-normal distance at which the dispersive stresses are mitigated as reported in figures 13 and 14.

Figure 16. Wall-normal distribution of $(a{,}b)$ the inner-normalised total normal and shear stresses for the different multiscale rough surfaces compared with the smooth-wall in blue-solid line. Wall-normal distribution of the ($c$) inner- and ($d$) outer-normalised dispersive shear stress, highlighting the extent of the roughness sublayer.

Figure 16($c$) illustrates the inner-normalised dispersive shear stress profiles with respect to the wall-normal distance normalised by the roughness length scale $h_{s}$. It shows that the maximum magnitude of the dispersive shear stresses rises up to nearly 0.16 $\%$ of the friction velocity, whereas its location changes across cases. It is observed that the wall-normal location of this maximum reduces as the roughness content ($h_{s}$) increases, with $\langle \tilde {u} \tilde {v} \rangle ^{+}_{max}$ occurring at $y \approx$ 2$h_{s}$ for Iter$_{123}$ while its occurrence is seen at $y \approx$ 5$h_{s}$ for Iter$_{1}$. The reason for the latter behaviour rises from the increase in $h_{s}$ when increasing the roughness content. This leads to a lack of similarity between cases suggesting that the increment in the aerodynamic roughness length scale does not necessarily induce a proportional change in the secondary motions. However, it is worth noting that $\langle \tilde {u} \tilde {v} \rangle ^{+}_{max}$ remains nearly constant across cases. Beyond their maximum value, the dispersive shear stresses appear to have reasonably similar decay rates in the outer region, despite the offset in the wall-normal direction.

The quantification of flow heterogeneity for the different surfaces is also examined in outer-scaling as presented in figure 16($d$), illustrating the relative contribution of the dispersive to the total shear stress. The profiles are shown to have a weak magnitude near the canopy, exhibiting even negative values right above the large cuboid for Iter$_{1}$, in line with the observations reported in figure 14. Farther away, the dispersive shear stresses reach their maximum between $0.2$ and $0.25\delta$, corresponding to the wall-normal location at which the secondary motions are most significant. In this scaling, the dispersive stresses are shown to amount to a substantial fraction of the total shear stresses with roughly 22 $\%$ at their maximum, while they seem to have a relatively better collapse in the outer region as opposed to the scaling presented in figure 16($c$). Their contributions rise up to approximately $0.6\delta$, corresponding to the turbulent stresses collapsing with the total shear stress profiles. This means that the inertial sublayer is completely submerged within the roughness sublayer (indicated in figure 16($d$) as the RSL). This behaviour is in contrast with homogeneous rough-wall flows, which generally exhibit a RSL as a small portion of the flow depth, and does not exceed the upper edge of the ISL. The current results also indicate that the dispersive stresses are mostly a function of the largest length scale features of the surface ($S$), while they only weakly vary with addition of small-scale roughness. The multiscale rough surfaces have therefore impacted not only the wall-normal location of the inertial sublayer as seen in figure 9, but also the overall structure organisation in the outer region.

Castro, Segalini & Alfredsson (2013) surveyed a wide range of rough-wall data published in the literature to examine the outer-flow behaviour using the diagnostic-plot method. This method has first been introduced by Alfredsson, Segalini & Örlü (2011) for smooth-wall flows, and consists of examining the dependence of the streamwise turbulence intensity profile with respect to the mean flow. Their assessment revealed the existence of a region spanning from the inertial sublayer up to almost the entire wake which varies linearly, whose extent increases with increasing scale separation. Accordingly, Castro et al. (2013) examined the applicability of this method over a comprehensive range of rough-wall data to reveal whether smooth and rough walls also exhibit outer-layer similarity in this form.

Their results showed that just like the smooth-wall flow, rough-wall flows also presented a linear variation except with a higher slope, shown to be independent of the roughness morphology. This rise in the turbulence intensity level was shown to be primarily due to the roughness function $\Delta U^{+}$, with marginal effect from the wake strength parameter. Using a similar approach, we try to assess whether a self-similar region still exists for the current multiscale rough surfaces using this scaling.

The total normalised streamwise r.m.s. profiles $\langle \sqrt {\tau _{xx}}/ U \rangle = \sqrt { \langle \overline {uu} + \tilde u \tilde u \rangle }/\langle U \rangle$ are plotted against the normalised streamwise velocity $\langle U \rangle /U_{\infty }$, as shown in figure 17($a$). The dispersive stress component has been included to account for the presence of spanwise mean flow heterogeneity, which would have been non-existent in the case of a two-dimensional boundary-layer flow. The profiles of the different cases have been plotted along with a smooth-wall profile and two asymptotes: the smooth-wall asymptote of Alfredsson et al. (2011) and the fully rough asymptote of Castro et al. (2013). The results show a good collapse both between the different cases as well as with the fully rough line, indicating that all the cases exhibit similar attributes as a fully rough flow, including Iter$_{1}$. This suggests that despite the lack of collapse observed in the profiles in figure 15, the surface condition influences both the turbulence intensity and the mean-flow proportionally, regardless of the presence of smaller or intermediate roughness length scales. As demonstrated by Castro et al. (2013), the increase in the streamwise stress level (slope) for rough-wall flows stems from the roughness function $\Delta U^{+}$. As such, they showed that by rescaling the turbulence-intensity profiles with the modified quantities $U' = U + \Delta U^{+}$ and $U'_{\infty } = U_{\infty } + \Delta U^{+}$, all profiles collapsed on to the smooth-wall asymptote.

Figure 17. Effect of the different multiscale rough surfaces on ($a$) the total normalised streamwise rms profiles in the diagnostic-plot form as proposed by Alfredsson et al. (2011). ($b$) The same data with the modified scaling $U' = U + \Delta U^{+}$ and $U'_{\infty } = U_{\infty } + \Delta U^{+}$ proposed by Castro et al. (2013).

After applying the modified scaling approach as seen in figure 16($b$), we notice in fact that the slopes of the different profiles line up relatively well with the smooth-wall asymptote, as reported by Castro et al. (2013). It should be noted that if the dispersive stress component $(\langle \overline {\tilde {u}\tilde {u}}\rangle )$ is omitted, the profiles will result in shallower slopes. Figure 17($b$) also suggests that the roughness function embodies both stress contributions, namely, the turbulence which stems from the temporal fluctuations and dispersive stresses which arise from spatial heterogeneities. This means that despite the large ramifications they cause to the base flow, the secondary motions are capable of reorganising the flow such that the mean and turbulence structures maintain a degree of local similarity between smooth and rough walls. However, it is worth recalling that despite this apparent collapse in the diagnostic plot, Townsend's similarity hypothesis remains clearly violated as a consequence of the large-scale secondary motions. Hence, the validity and application of outer-layer similarity hypothesis should be carefully considered when examining flows developing over such surfaces.

4. Conclusions

In this experimental study, the characteristics of a turbulent boundary-layer flow developing over regular multiscale rough surfaces have been examined. Four surfaces representing a multiscale distribution of roughness elements, based on a Sierpinski fractal-like pattern, made of superimposing size-decreasing self-similar cuboids are considered. A floating-element drag balance along with PIV measurements allowed us to assess the impact of these multiscale rough surfaces on the skin friction as well as the mean and turbulent flow properties.

Drag-balance measurements revealed that the skin-friction coefficients reached their asymptote with Reynolds number (i.e. fully rough regime) for surfaces containing at least the intermediate or smaller scales (Iter$_{12}$ or Iter$_{13}$), in contrast with the surfaces containing only the largest roughness element (Iter$_{1}$) that was shown to be Reynolds number dependent (despite the high value of $\Delta U^+$). The results showed that the magnitude of the skin friction is dependant of the scale hierarchy of the surface roughness, and increases with frontal solidity of the added roughness ($\alpha$).

The impact of the multiscale roughness on the turbulent boundary layer is also examined using velocity field measurements. In the canopy region, a strong shear layer formed above the largest cuboid separates at its trailing edge leading to the formation of a separation bubble past the obstacle. The streamwise extent of this bubble (reattachment length) was seen to be insensitive to the roughness content while its intensity decreases with the presence of smaller roughness length scales. This behaviour is caused by the increased turbulent mixing and wall drag at the top of the largest cuboid, owing to the increased interactions of the shear layer with a broader range of roughness scales. Additionally, the PIV-based pressure fields showed the canopy to be dominated by alternating high- and low-pressure zones associated with the flow accelerating and decelerating upstream and downstream of the first roughness iteration, respectively. This allowed the assessment of the cross-sectional drag, and revealed that these large-scale roughness elements are responsible for nearly 60 $\%$ of the surface drag.

The aerodynamic parameters are examined using the horizontally averaged mean velocity profiles. The virtual origin as well as the roughness function have been determined by means of the indicator function. The latter highlighted an inertial sublayer occurring between $0.2$ and $0.3 \delta$, which is farther in comparison with the classical smooth and rough walls. The roughness function estimates were then used to determine the equivalent sandgrain roughness height $h_{s}$ whose values indicated that all the cases can be considered in the fully rough regime, including Iter$_{1}$ ($h_{s}^{+} > 400$). A ratio of the aerodynamic roughness length scale between successive iterations ($\rho _{h_{s}}$) was defined and shown to vary linearly with a suitable geometrical parameter of the rough surface ($\alpha$). This relationship together with the aerodynamic properties of the first iteration can be sufficient to describe the roughness length scale of the subsequent multiscale surfaces (and hence the surface drag).

The cross-plane velocity fields revealed a significant spanwise heterogeneity in the outer layer for all surfaces examined. This spanwise flow modulation is shown to be caused by the presence of Prandtl's secondary flows of the second kind, which created HMP and LMP alternating above the ridges and the valleys, respectively. As shown from previous studies, this feature occurs as a consequence of the dominant length scale of the flow interacting with $\delta$-scale roughness features. In this study, the largest characteristic length scale of the roughness is the spacing $S$ between the largest cuboid which results in $\delta /S \approx O(1)$.

Triple decomposition of the velocity fields allowed the qualitative and quantitative assessment of the spanwise heterogeneity in the form of dispersive stresses. The results showed that, at their maximum, these form-induced stresses amount to more than $20\,\%$ of the total stresses and can extend up to 0.6$\delta$. Additionally, the validity of Townsend's similarity hypothesis was assessed by examining the velocity deficit as well as the turbulence intensity profiles. Lack of outer-layer similarity in both the mean and turbulence statistics is reported despite the required conditions being satisfied ($4600 \leq Re_{\tau } \leq 9700$ and $0.05 \leq h_{s}/\delta \leq 0.1$), while a good degree of collapse between cases is observed in the mean flow.