2.1. Case set-up

Figure 2 depicts a half-channel set-up containing a wall with cube-like roughness elements. This is a representative domain over which flow is computed for the DNS runs. The bottom surface is a rough wall comprising pits (or valleys) and protrusions (or peaks) in the form of cube-like elements. The $x$, $y$ and $z$ represent streamwise, spanwise and wall-normal directions. Periodic conditions are applied to the lateral boundaries. A stress-free condition with no penetration (${\partial u}/{\partial z} = {\partial v}/{\partial z} = w = 0$) is imposed on the upper boundary. A streamwise pressure gradient acts as the forcing that drives the flow. A friction Reynolds number $Re_\tau \approx 400$ is used for all cases in this study, where $Re_\tau$ is defined as

(2.1)\begin{equation} R e_\tau=\frac{u_\tau (L_z-h_2)}{\nu}. \end{equation}

Note that the wall-normal length $L_z$ used in (2.1) includes the depth of the negatively skewed features or valleys $h_2$. This depth $h_2$ is subtracted from $L_z$ to make a good estimate of the half-channel height. The friction velocity $u_\tau$ is obtained using

(2.2)\begin{equation} u_\tau= \sqrt{\left|\frac{\text{d} \langle \bar{p}\rangle }{\text{d}\kern0.07em x}\right|\frac{V_f}{\rho A_p}} \approx \sqrt{\left|\frac{\text{d} \langle \bar{p}\rangle }{\text{d}\kern0.07em x}\right|\frac{(L_z-h_2)}{\rho}}, \end{equation}

where ${\rm d}\langle \bar {p}\rangle /{\rm d}\kern0.07em x$ refers to the mean streamwise pressure gradient, and $\rho$, $V_f$ and $A_p$ refer to the fluid density, half-channel fluid volume and wall-parallel area, respectively. Here $\langle \cdot \rangle$ denotes the streamwise and spanwise averaging (or double-averaging) operation, and $\overline {(\cdot )}$ denotes time averaging.

Figure 2. A sketch illustrating a periodic half-channel flow set-up over cube-like peaks and valleys. Here, $h_1$ denotes peak height, $h_2$ denotes valley depth, and $b$ denotes the width of the cube-like element. The boundary layer thickness $\delta$ here is at least 6 times the maximum of $h_1$ and $h_2$.

Note that $V_f$ includes the fluid volume within valleys, and the approximation sign in (2.2) is used because the peak height is not always equal to the valley depth, causing minor differences in $V_f$ in some cases.

The size of the computational domain is determined such that $L_x>6L_z$ and $L_y>3L_z$. Lozano-Durán & Jiménez (Reference Lozano-Durán and Jiménez2014) have shown that such a domain would be sufficiently large to ensure the accuracy of first- and second-order statistics for plane channel flow. Note that $L_x$ and $L_y$ represent the streamwise and spanwise domain lengths, respectively. For rough-wall flows, this domain size should be a conservative estimate since DNS studies in Coceal et al. (Reference Coceal, Thomas, Castro and Belcher2006) and Leonardi & Castro (Reference Leonardi and Castro2010) produce good mean flow statistics with smaller domains $L_x \times L_y \times L_z = 4h \times 4h \times 4h$ and $8h \times 6h \times 8h$, respectively ($h$ being the height of the cube). Other studies, such as Chung et al. (Reference Chung, Chan, MacDonald, Hutchins and Ooi2015), have shown that even further reduction in computational domain is possible in the spanwise direction, and such minimal-span channels have been shown to capture accurate mean drag characteristics (at least for sinusoidal roughness) when $L_y > k/0.4$ and $L_y^+>100$. Furthermore, since the roughness is regular, $L_x$ and $L_y$ are also integral multiples of the length of the repeating tile to ensure a periodic domain. Here, a repeating tile represents the smallest unit that when repeated in the streamwise and spanwise directions produces the entire surface.

A uniform grid has been utilized with grid resolution such that $\varDelta _x^+ = \varDelta _y^+<5$ and $\varDelta _z^+<3$, where these represent streamwise, spanwise and wall-normal grid spacings, respectively, normalized by the viscous length scale ($\delta _v = \nu /u_\tau$). It is important to note that while $\varDelta _z^+$ may appear to be large for DNS, this resolution is based on friction velocity $u_{\tau }$ obtained from (2.2). The $\delta _v$ value computed based on this $u_{\tau }$ includes the form drag component. If $\delta _v$ is computed from the local viscous stress at the mean surface comprising the bottom wall, then the grid resolution is at most 0.66 plus units. Note that the local viscous stress $\mu \,\partial {u}/\partial {n}$ here does not pertain to the cube surfaces, where the corresponding value would spike at the leading edge due to the strong shear rate. A similar grid resolution $\varDelta _x^+ = \varDelta _y^+<4.5$ and $\varDelta _z^+\approx 2.5$ was used by Zhang et al. (Reference Zhang, Zhu, Yang and Wan2022) for deep canopy flows. The grid spacing is also comparable to that in previous studies of rough-wall DNS, where Coceal et al. (Reference Coceal, Thomas, Castro and Belcher2006) used $\varDelta ^+_{x,y,z}$ in the range 7.8–15.6, and Leonardi & Castro (Reference Leonardi and Castro2010) used $\varDelta ^+_{x,y} = 19$. For further reference, DNS studies by Forooghi et al. (Reference Forooghi, Stroh, Magagnato, Jakirlić and Frohnapfel2017) have found $\varDelta ^+_{x,y,z_{max}} = 3.5$ to be more than sufficient. Yuan & Piomelli (Reference Yuan and Piomelli2014b) utilize $\varDelta ^+_{x} = 12$ and $\varDelta ^+_{y} = 6$, whereas Jouybari et al. (Reference Jouybari, Yuan, Brereton and Murillo2021) employ $\varDelta ^+_{x} = 7.5$ and $\varDelta ^+_{y} = 6.3$ for obtaining accurate flow statistics in their rough-wall channels.

All simulations are performed using LESGO: a parallel pseudo-spectral large-eddy simulation code (see https://lesgo-jhu.github.io/lesgo), a solver of the filtered Navier–Stokes equations with pseudo-spectral discretization in the streamwise and spanwise directions, and a second-order finite difference in the wall-normal direction. Code LESGO and modified versions of the code have been used extensively for studies involving channel flows, including Bou-Zeid, Meneveau & Parlange (Reference Bou-Zeid, Meneveau and Parlange2005), Anderson & Meneveau (Reference Anderson and Meneveau2011), Abkar & Porté-Agel (Reference Abkar and Porté-Agel2012), Zhu et al. (Reference Zhu, Iungo, Leonardi and Anderson2017), Yang et al. (Reference Yang, Xu, Huang and Ge2019), and several others. Immersed boundary conditions (Anderson Reference Anderson2013) are used at the solid boundaries comprising the rough wall. A Courant–Friedrichs–Lewy number (${\rm CFL} = u\,\Delta t/\Delta x$) of 0.2 is employed to automatically adjust the time step, which is advanced via the second-order Adams–Bashforth scheme.

2.2. Roughness generation

Roughness elements are distributed in various arrangements to generate different rough surfaces. We will vary one statistic at a time, but repeat the process and vary many roughness statistics. Figures 3(a), 3(c) and 3(d) depict these elements aligned in the spanwise (AY), streamwise (AX) and $45^\circ$ (AXY) directions, respectively. Their staggered equivalents (S and SXY) are shown in figures 3(b) and 3(e). Figure 3( f) shows another surface (NV) where the negatively skewed features (valleys) are removed. The rough surfaces also include variations in $\lambda _p$, as indicated by figures 3(g) and 3(h), peak height to width ratio ($h_1/b$), and valley depth to width ratio ($h_2/b$). Table 1 provides the naming convention used, and table 2 lists all 36 different surfaces considered for the DNS study. For each case, the nomenclature is as follows: [Arrangement][Packing density][Roughness height], where arrangement could be AX, AY, AXY, S, SXY, NV, packing density could be L1, L2, L3, and roughness height could be H1, H2, H3. For example, the L1, L2 and L3 in AYL1H1, AYL2H1 and AYL3H1 correspond to $\lambda _p$ magnitudes 6.25 %, 11.1 % and 25 %, respectively. Similarly, the H1, H2 and H3 in AYL1H1, AYL1H2 and AYL1H3 denote variations in the heights of the peaks as $h_1 = b, 1.2b, 0.8b$. Note that in this context, the conditions $h_1/h_2 = 1$, $h_1/h_2 > 1$ and $h_1/h_2 < 1$ correspond to zero skewed, positively skewed and negatively skewed surfaces. The separation lengths between the roughness elements are listed as $4b$, $3b$ and $2b$ for L1, L2 and L3 configurations, with $b/L_z$ assuming values 0.1375, 0.1111 and 0.1429 for H1, H2 and H3 cases, respectively.

Figure 3. Top view of various rough-surface types considered for the DNS study: (a) AYL1H1, (b) SL1H1, (c) AXL1H1, (d) AXYL1H1, (e) SXYL1H1, ( f) NVL1H1, (g) AYL2H1 and (h) AYL3H1. Here, $k$ is the elevation of the surface. The peaks have positive $k$ (depicted as yellow), and the valleys have negative $k$ (depicted as blue).

Table 1. Nomenclature for the 36 configurations of the parameter space studied.

Table 2. Details of the DNS case-set-up, where $\lambda _p$, $h_1/b$, $h_2/b$ denote planar packing density, peak height to width, and valley depth to width ratios, respectively; $N_{x,y,z}$ denotes the number of grid cells in the streamwise ($x$), spanwise ($y$) and wall-normal ($z$) directions; $N$ denotes the number of peaks plus valleys in the case. Note that the surface coverage densities of the NV cases are half the value of their AX (and AY) counterparts.

The definitions for a few roughness statistics $k_{rms}$, $k_a$, $k_t$, $Sk$, $ES$, $Ku$ and $\lambda _p$ (which are discussed further in table 3) for the rough surfaces thus generated are listed as follows:

(2.3)$$\begin{gather} \bar{k}=\frac{1}{L_x L_y}\int_{L_x}\int_{L_y} k(x,y)\,{\rm d}\kern0.07em{x}\,{\rm d}{y}, \end{gather}$$

(2.4)$$\begin{gather}k_a=\frac{1}{L_x L_y}\int_{L_x}\int_{L_y} |k(x,y)-\bar{k}|\,{\rm d}\kern0.07em{x}\,{\rm d}{y}, \end{gather}$$

(2.5)$$\begin{gather}k_{rms}=\sqrt{\frac{1}{L_x L_y}\int_{L_x}\int_{L_y}|k(x,y) -\bar{k}|^2\,{\rm d}\kern0.07em{x}\,{\rm d}{y}}, \end{gather}$$

(2.6)$$\begin{gather}k_{t}=\text{max}(k(x, y)) - \text{min}(k(x,y)), \end{gather}$$

(2.7)$$\begin{gather}Sk=\frac{1}{L_x L_y}\int_{L_x}\int_{L_y}(k(x,y)-\bar{k})^3\,{\rm d}\kern0.07em{x}\,{\rm d}{y}/k_{rms}^3, \end{gather}$$

(2.8)$$\begin{gather}Ku=\frac{1}{L_x L_y}\int_{L_x}\int_{L_y}|k(x,y)-\bar{k}|^4\,{\rm d}\kern0.07em{x}\,{\rm d}{y}/k_{rms}^4, \end{gather}$$

(2.9)$$\begin{gather}ES=\frac{1}{L_x}\int_{L_x} \left|\frac{\partial k}{\partial x}\right|{\rm d}\kern0.07em x, \end{gather}$$

(2.10)$$\begin{gather}\lambda_p=\frac{A_p}{L_x L_y}. \end{gather}$$

Table 3. Hydrodynamic properties of roughness. The geometric parameters listed include average and root mean square roughness height $k_a$ and $k_{rms}$, skewness $Sk$, kurtosis $Ku$, maximum peak to trough height $k_t$, effective slope (along the $x$ direction) $ES$, streamwise and spanwise correlation lengths $Rl_{x}$ and $Rl_y$ (defined in (4.6), (4.7)). Here, $d$ denotes zero-plane displacement height, and $z_0$ is the effective roughness height.

The packing density $\lambda _p$ in (2.10) is defined in the same way as plan solidity in Chung et al. (Reference Chung, Hutchins, Schultz and Flack2021). Here, $A_p$ comprises the plan area of roughness elements, which in this case includes the base area of peaks and valleys.

Additional constraints have been placed to ensure that the grid surface matches the cube surface for accurate roughness representation in the simulations. This is done by enforcing $b/\varDelta _x$, $b/\varDelta _y$, $h_1/\varDelta _z$ and $h_2/\varDelta _z$ to be integers, where $\varDelta _x$, $\varDelta _y$ and $\varDelta _z$ stand for the grid spacing in the streamwise, spanwise and wall-normal directions. This has been implemented for all cases except the $45^\circ$ alignment rough surfaces.

As can be seen from figure 4, the grid is in perfect alignment with the roughness element in AYL1H1 as opposed to AXYL1H1. As a result, there are minor variations in the average roughness height $k_a/b$, root mean square roughness height $k_{rms}/k_a$, and kurtosis $Ku$ for AXYL$i$H1 and SXYL$i$H1 surfaces when compared with AYL$i$H1 and SL$i$H1, where $i=1, 2, 3$. This can be noticed in table 3. A finer resolution has been opted for in all such cases to minimize these variations. The streamwise and spanwise grid resolutions, in terms of number of cells per roughness element, are given by $b/\varDelta _x = b/\varDelta _y = 12\unicode{x2013}14$ for AX, AY, S and NV surfaces. These grid resolutions are approximately 16–18 for AXY and SXY surfaces. The wall-normal grid resolution is in the range $b/\varDelta _z = 20\unicode{x2013}25$ for all cases.

Figure 4. Top and side views of DNS grids used for two cases: (a,c) AYL1H1 and (b,d) AXYL1H1. The side views are streamwise wall-normal ($x$–$z$) planes cut at the middle of the cubes. The top views are wall-parallel ($x$–$y$) planes cut near the roughness peak.

3.1. Instantaneous flow field

We begin by presenting the instantaneous flow field as an introduction to the dataset. From the contours in figure 5, an overall decrease in instantaneous bulk velocity can be observed with increasing $\lambda _p$. A region of reduced streamwise velocity can be noted in the immediate wake of the protrusions. The flow in the pits can be observed to be less energetic, and their interactions with the outer layer are minimal.

Figure 5. Streamwise instantaneous velocity contours at $y/L_y \approx 0.43$ for (a) AYL1H1, (b) AYL2H1 and (c) AYL3H1.

Figure 6 shows the velocity contours on a wall-parallel plane near the peak height for different roughness element arrangements at the same $\lambda _p$. By visual inspection, contours in figures 6(a) and 6(b) are observed to have similar and intermediate average streamwise velocities. Figure 6(c) contains the highest and figure 6(d) the lowest average streamwise velocities amongst these four cases. In all cases, streamwise streaks of high and low velocities can be observed near the roughness peak. These streaks are noticeably longer in certain cases, such as in figure 6(c), as the wakes overlap each other. These wake interactions become more pronounced at higher $\lambda _p$ and result in a flow-sheltering effect as observed in densely packed surfaces such as in Xu et al. (Reference Xu, Altland, Yang and Kunz2021).

Figure 6. Streamwise instantaneous velocity contours at $z/b \approx 1.0$ for (a) AYL1H1, (b) SL1H1, (c) AXL1H1 and (d) AXYL1H1.

Although we will pursue a quantitative analysis in the following subsections, here we may already conclude first that the valleys do not contribute significantly to the drag, and second that two-point rough-wall statistics are important to the characterization of the equivalent sand-grain roughness height $k_s$.

3.2. Statistical convergence

Before presenting the flow statistics, we first check the statistical convergence of our data. The statistical convergence can be evaluated with the help of the mean momentum budget following Coceal et al. (Reference Coceal, Thomas, Castro and Belcher2006), Leonardi & Castro (Reference Leonardi and Castro2010), Xu et al. (Reference Xu, Altland, Yang and Kunz2021) and Zhang et al. (Reference Zhang, Zhu, Yang and Wan2022). The equation is comprised of viscous diffusion, turbulent transport, dispersive stress and pressure-gradient terms:

(3.1)\begin{equation} \nu\,\frac{\text{d}\langle \bar{u} \rangle}{\text{d}z} - \langle \overline{u^\prime w^\prime} \rangle - \langle u^{\prime \prime}w^{\prime \prime} \rangle = \frac{1}{\rho}\,\frac{\text{d} \langle\bar{p}\rangle}{\text{d}\kern0.07em x}\,(z-\delta). \end{equation}

It is considered to be statistically converged when the sum of these stresses is a linear function of the wall-normal distance ($z$), i.e. when the budget balances without the unsteady term. Note that (3.1) is obtained by integrating the $x$ component of the double-averaged momentum equation

(3.2)\begin{equation} \frac{\text{d}}{\text{d}z}\left(\nu\,\frac{\text{d}\langle \bar{u} \rangle}{\text{d}z} - \langle \overline{u^\prime w^\prime} \rangle - \langle u^{\prime \prime}w^{\prime \prime} \rangle \right) - \frac{1}{\rho}\,\frac{\text{d} \langle\bar{p}\rangle}{\text{d}\kern0.07em x} - f = 0 \end{equation}

found in (15) of Nikora et al. (Reference Nikora, Ballio, Coleman and Pokrajac2013), with static roughness in place of the mobile roughness assumed by the authors. Here, $f$ represents the drag force. As $f$ is non-zero within the roughness-occupied region, (3.1) applies to the region outside the roughness only. The averaging procedure corresponds to the double-averaging method, common in the roughness literature, introduced by Raupach & Shaw (Reference Raupach and Shaw1982) for flows within vegetation canopies. The spatial averaging used in the context of (3.1) and (3.2) is superficial, as per the definition for superficial averaging in Schmid et al. (Reference Schmid, Lawrence, Parlange and Giometto2019).

Figure 7 shows these terms for a few cases. These cases have been selected so that they include different $\lambda _p$, $Sk$ and roughness arrangements while maintaining brevity. We see that the total stresses in all cases follow a linear function of $z/b$, therefore our DNS are statistically converged. In addition, we make the following observations. The Reynolds stress $R_{13}$ can be observed to be the largest in magnitude as compared to the viscous stress $\tau _v$ and dispersive stress $D_{13}$ outside the roughness-occupied layer. The roughness-occupied layer is marked by the vertical line in figure 7. In the roughness-occupied layer, i.e. below this marking, $\tau _v$ is observed to contain two maxima, one caused by the surface at $z=0$, and the other near the cube height $z=h_1$. The latter is also observed in other DNS studies by Xu et al. (Reference Xu, Altland, Yang and Kunz2021) and Zhang et al. (Reference Zhang, Zhu, Yang and Wan2022). Meanwhile, $R_{13}$ attains a maximum just above the peak height, and declines to zero as $z$ approaches the outer layer; $D_{13}$ is small and contained within the roughness region for most cases. For $z/b < 0$, i.e. below the surface in the region of valleys, all stresses are observed to be negligible. Note that figure 7(e) starts from $z/b = 0$ as the case NVL1H1 contains no valleys.

Figure 7. The stress budget plot comprising Reynolds stress $R_{13}$, dispersive stress $D_{13}$, viscous stress $\tau _v$ and total stress $\tau _t$ for cases: (a) AYL3H1, (b) SL1H2, (c) SL2H3, (d) AXL2H3, (e) NVL1H1 and ( f) AXYL3H1. The black solid line corresponds to $1-z/\delta$.The vertical dashed line denotes the peak roughness height.

It might be worth noting certain differences in dispersive stress $D_{13}$ in figure 7. For instance, most noticeable would be the differences in cases AYL3H1 and NVL1H1: $D_{13}$ is observed to be non-negligible within the roughness-occupied region for AYL3H1, whereas $D_{13}$ for NVL1H1 inside the sublayer is almost zero. Figure 8 aims to provide an intuitive explanation for this. Here, the dispersive streamwise and wall-normal fluctuating components, $u^{\prime \prime }$ and $w^{\prime \prime }$, are shown for averaged repeating tiles for cases NVL1H1 and AYL3H1. Three planes, one within and the others at and above the roughness crest, are chosen to display their scatter. When $u^{\prime \prime }$ and $w^{\prime \prime }$ predominantly lie in the second or fourth quadrant, the dispersive stress $-\langle u^{\prime \prime }w^{\prime \prime }\rangle$ becomes positive. The negative slope of the linear fit is also indicative of this. This can be seen in figure 8(c), which corresponds to NVL1H1 just above the roughness, and figure 8(d), which corresponds to AYL3H1 within the roughness. These may indirectly imply mean flow inhomogeneity, which is due to the presence of surface roughness and/or secondary flows. On the other hand, when the points are more evenly spread across all quadrants, as is the case in figures 8(a), 8(b) and 8(e), or too close to the origin as in figure 8( f), this dispersive stress becomes negligible.

Figure 8. Quadrant analysis for the dispersive stress components of $D_{13}$ for cases NVL1H1 at (a) $z/b=0.5$, (b) $z/b=1$, (c) $z/b=2$ and AYL3H1 at (d) $z/b=0.5$, (e) $z/b=1$ and ( f) $z/b=2$. The solid black line depicts a linear fit on the data points.

3.3. Mean velocity profiles

The mean flow statistics are presented in this subsection. Additionally, the relevant rough surface parameters and the hydrodynamic properties determined from mean velocity profiles are also listed in table 3.

Figure 9(a) shows the inner scaled mean velocity profiles, comparing the spanwise-aligned and staggered cases, i.e. the AY cases and the S cases. The following observations can be duly noted. First, the higher $\lambda _p$ cases of surface coverage density 25 %, comprising AYL3H1, AYL3H2, AYL3H3 and their staggered counterparts SL3H1, SL3H2 and SL3H3, produce lower magnitudes of mean velocity. This is expected for roughness in the k-type regime. Second, the staggered arrangement produces similar velocity profiles to the spanwise-aligned arrangement. This is not unexpected, and similar observations were made in Yang et al. (Reference Yang, Sadique, Mittal and Meneveau2016). Due to a large spanwise distance between neighbouring roughness elements in the spanwise direction, staggering the roughness element in the streamwise direction has little effect on the mean flow. Finally, at constant $\lambda _p$, cases with higher skewness ($Sk$) produce lower mean velocity profiles. This role of $Sk$ is apparent as we notice increased drag with increased $Sk$ (consistent with systematic studies performed by Flack et al. Reference Flack, Schultz and Barros2020), and we observe this effect in all our cases, irrespective of arrangement. For example, comparing AYL1H1, AYL1H2 and AYL1H3 from figure 9(a) or AXL1H1, AXL1H2 and AXL1H3 from figure 9(b) at the same $z^+$, one can notice that the mean velocities in the log layer are slightly lower for the surface with higher $Sk$ (H3 cases being lower than H2 and H1 cases). In the present dataset, variations in $Sk$ are introduced by varying the height of the peaks, while maintaining the depth of the dips constant, therefore any changes due to $Sk$ could also be viewed as due to the roughness height. Similarly, as variations in $k_{rms}/k_a$ are introduced by varying the surface coverage density, any changes due to $k_{rms}/k_a$ could also be viewed as due to surface coverage density.

Figure 9. Mean streamwise velocity profile comparing (a) spanwise-aligned and staggered cases, and (b) spanwise-aligned and streamwise-aligned cases.

Similarly, figure 9(b) compares mean velocity profiles for roughness with streamwise and spanwise alignment, i.e. AX cases and AY cases. When the roughness features are aligned in the streamwise direction, higher mean velocities can be observed as compared to those aligned in the spanwise direction. An intuitive explanation for this is the unobstructed passage of fluid between two rows of roughness elements giving the fluid less drag when they are aligned in the streamwise direction. This is also evident from the instantaneous velocity contours in figure 6(c).

Figure 10 includes the mean velocity profiles from the $45^\circ$ arrangements (XY), with figure 10(a) comparing AY, AX, AXY, and figure 10(b) comparing S and SXY. For the same $\lambda _p$, the surfaces follow ${\rm AX} > {\rm AY} >{\rm AXY}$ in figure (a), and ${\rm S} > {\rm SXY}$ in figure (b) for the order of mean velocity magnitudes in the log-law region. The orientation of the roughness element is found to be an important factor here. The AXY and SXY surfaces possess a higher frontal area (i.e. higher solidity) as compared to their counterparts AX, AY or S, which is the reason for the higher drag observed for these surfaces.

Figure 10. Mean streamwise velocity profile comparing (a) spanwise, streamwise and $45^\circ$ aligned cases, and (b) staggered and staggered $45^\circ$ cases.

Figure 11(a) shows the mean velocity profiles comparing three pairs of surfaces with and without valleys. It can be observed that the valleys do contribute to mean velocity reduction, but by only a small amount. This reduced mean velocity can be inferred to be caused by the additional pressure drag contributions from the valleys. In support of this claim, figures 11(b) and 11(c) are shown, which contain mean pressure contours that are averaged over repeating tiles in both $(x, y)$ directions in AYL1H1 and NVL1H1. In these pressure contours, which are normalized by $\rho u_{\tau }^2$, the volume-averaged mean pressure is taken as the reference pressure. The additional pressure drop can be noticed clearly in valleys. Note that the observation here cannot be generalized to all dips. In particular, dips or valleys will play a significant role if they occupy a considerable part of the surface area, in which scenario a new bottom surface forms and the protrusions between the dips can be viewed as surface roughness.

Figure 11. Details of the effects of valleys. (a) Mean streamwise velocity profile comparing cases with and without valleys. The solid and dotted lines depict spanwise-aligned and no-valleys cases, respectively. (b) Mean pressure contours for a repeating tile in NVL1H1. (c) Mean pressure contours for a repeating tile in AYL1H1.

Next, we measure the hydrodynamic properties of the roughness. The effective roughness height $z_0$ and zero-plane displacement height $d$ for the rough surfaces are determined from the mean streamwise velocity profiles as

(3.3)\begin{equation} \frac{\langle \bar{u} \rangle}{u_\tau}=\frac{1}{\kappa} \log \frac{z-d}{z_0}, \end{equation}

where $\langle \bar {u} \rangle$ and $z$ are taken from the log-law region, which has been taken to be in the range $(z-h_1)^+ \approx 30\unicode{x2013}100$. Further evidence of the existence of a log law in this range will be discussed later with the help of figure 13. Linear regression of (3.3) yields the resulting $z_0$ and $d$ reported in table 3. Here, the von Kármán constant is set to $\kappa = 0.4$. There can be alternate approaches to calculate $d$ from the centroid of the total force distribution following Jackson (Reference Jackson1981) and Coceal et al. (Reference Coceal, Thomas, Castro and Belcher2006). However, it is also observed that the centre of the force could sometimes underestimate $d$ when the flow primarily interacts with the top portion of the roughness, as seen in figure 7 of Xu et al. (Reference Xu, Altland, Yang and Kunz2021). Here, we prefer regression fitting over this approach for the aforementioned reason.

Since we assume the value of the von Kármán constant $\kappa = 0.4$ for our regression procedure, checking the changes in measured $z_0/b$ due to variations in $\kappa$ is important. Figure 12 shows the sensitivity of calculated $z_0/b$ and $d/b$ with various values of $\kappa$. It can be seen that $z_0/b$ decreases and $d/b$ increases with $\kappa$. While we do see noticeable variations in these values, it is worth noting that these changes would not be significant on the log scale, which is the relevant scale for the drag produced by these surfaces. Also, for a given value of $\kappa$, the differences in $z_0/b$ and $d/b$ are almost the same, which means that it is reasonable to compare these values for the different surfaces.

Figure 12. Plot showing the sensitivity to variations in von Kármán constant $\kappa$ for spanwise-aligned (AY) and staggered (S) cases for the calculated (a) effective roughness height $z_0/b$, and (b) zero-plane displacement height $d/b$.

Further, we note that determining the roughness height $z_0$ or equivalent sand-grain roughness height $k_s$ requires the flow to be in the fully rough regime. Based on Jiménez (Reference Jiménez2004), flow is considered fully rough when $k_s^+ > 80$. In Jouybari et al. (Reference Jouybari, Yuan, Brereton and Murillo2021), the authors have chosen $k_s^+ > 50$ to be their definition for the fully rough regime. We have verified that all surfaces except AXL1H3 and AXL3H3 satisfy this latter condition, the $k_s^+$ values for which stand at 42.6 and 36 respectively. However, it is important to note that the exact value of $k_s^+$ for which this transition to fully rough flow happens is unknown and also depends on the type of roughness being considered Flack & Schultz (Reference Flack and Schultz2010). Even so, we could still measure $\Delta U^+$ and $z_0$, and discuss their variations with respect to the roughness as we are making these comparisons at a fixed Reynolds number ($Re_{\tau }=400$).

Figure 13 illustrates the quality of this log-law fit. The mean velocity plots can be observed to show a good collapse in the log layer region for all cases. The collapse is also a sign of mean flow universality, wherein such two-parameter forms of mean velocity profiles have been know to be adequate irrespective of the rough surface for $h/\delta \approx 0.03\unicode{x2013}0.5$. Figure 14 shows mean velocity profiles in defect form. This is being presented here as evidence for outer-layer similarity. Reasonably well collapsed regions can be observed farther away from the wall for the smooth- and rough-wall cases, underscoring a universality in mean flow behaviour in the outer layer irrespective of surface conditions for fully rough flow.

Figure 13. Plot showing the log-law collapse for all cases: (a) AY, S cases; (b) AX, NV, AXY and SXY cases.

Figure 14. Velocity defect profiles for all surfaces: (a) AY, S cases; (b) AX, NV, AXY and SXY cases.

The mean velocity profile for the smooth wall in figures 14(a) and 14(b) was obtained from DNS of a half-channel with $Re_{\tau } = 400$. Similarly, figure 15 depicts outer-layer similarity in the streamwise Reynolds stress profiles.

Figure 15. Normalized streamwise Reynolds stress profiles for all surfaces: (a) AY, S cases; (b) AX, NV, AXY and SXY cases.

Figure 16 summarizes the $z_0$ magnitudes for all cases, and the data are consistent with the various mean velocity profiles observed so far. Figure 16(a) compares the spanwise-aligned (AY) and staggered (S) cases. The $z_0/b$ magnitudes are observed to be similar between the AY and S cases when comparing cases with similar $\lambda _p$ and $Sk$ magnitudes. When cases with similar packing densities $\lambda _p$ and orientation but different $Sk$ are compared – for instance, AYL1H1, AYL1H2 and AYL1H3 – $z_0/b$ increases with $Sk$. This holds true for all cases listed in table 3. When cases with similar $Sk$ and orientation but different $\lambda _p$ are compared – for instance, AYL1H1, AYL2H1 and AYL3H1 – $z_0/b$ decreases with decreasing $k_{rms}/k_a$. This is observed true for all cases except the AX surfaces, which will be discussed subsequently. It should be of interest to note that low $k_{rms}/k_a$ magnitudes actually correspond to high $k_{rms}/b$ and $\lambda _p$. The normalization with $k_a$ causes this as an increment in $\lambda _p$ increases both $k_a$ and $k_{rms}$, but $k_a$ increases at a rate faster than $k_{rms}$. Figure 16(b) compares the spanwise-aligned (AY) and streamwise-aligned (AX) cases. The effects due to changes in roughness element alignment are evident as the AX surfaces in figure 16(b) show notably lower $z_0/b$ magnitudes when compared to AY surfaces. For instance, even the densely packed surfaces AXL3H1, AXL3H2 and AXL3H3, which have higher solidities, are lower in $z_0/b$ than the coarsely packed surfaces AYL1H1, AYL1H2 and AYL1H3. Moreover, it can be noticed that $z_0/b$ first increases then decreases with decreasing $k_{rms}/k_{a}$ (or increasing $\lambda _p$) for H1 and H3 cases in AX surfaces. For higher $\lambda _p$ surfaces AXL2H2 and AXL3H2, the increase in $z_0/b$ is not as significant as it is for AYL2H2 and AYL3H3. These observations showcase a flow-sheltering effect in AX surfaces. Trends with changes in $Sk$ for AX surfaces are similar to those observed in AY ones, with $z_0/b$ increasing with increasing $Sk$. Figure 16(c) reports the effect of valleys on $z_0/b$. The presence of valleys in the AY cases leads to higher $k_{rms}/k_a$ compared to the NV cases, but the $z_0/b$ magnitudes are very close between the AY and NV cases, with some noticeable difference at low $k_{rms}/k_a$ magnitudes. Figure 16(d) presents spanwise-aligned (AY), $45^\circ$ alignment (AXY) and their respective staggered equivalents (S and SXY). A significant increase in $z_0/b$ can be observed as alignment is changed to $45^\circ$ when comparing AY and AXY. Also, SXY falls in the intermediate, with $z_0/b$ increasing due to a change in alignment when compared with S. However, this effect is less pronounced as $\lambda _p$ increases or $k_{rms}/k_a$ decreases.

Figure 16. The effective roughness length to width ratio ($z_0/b$) for the various cases: (a) spanwise-aligned versus staggered; (b) spanwise-aligned versus streamwise-aligned; (c) spanwise-aligned with and without valleys; (d) spanwise-aligned/staggered and their $45^\circ$ rotated equivalents.

3.4. Reynolds and dispersive stresses

We present the Reynolds stress and dispersive stress results. These results are usually not relevant to roughness modelling but are fundamental aspects of the flow, and we present them here for completeness.

Figure 17 shows the normal components of the horizontally averaged Reynolds and dispersive stresses for the five spanwise-aligned (AYL1H1, AYL2H1, AYL3H1, AYL1H2, AYL1H3) and staggered (SL1H1, SL2H1, SL3H1, SL1H2, SL1H3) cases. At first glance, all stresses are negligible in the valleys below the surface at $z=0$, and the spanwise and wall-normal components of the dispersive stresses are small throughout the channel. Second, it can be observed that the corresponding streamwise components of both stresses ($uu$) are higher than the respective spanwise ($vv$) and wall-normal ($ww$) components for all cases. Third, the maximum value for dispersive stress lies within the roughness-occupied layer, whereas for the Reynolds stress, the maxima lie in the neighbourhood of the roughness peak. Moreover, the dispersive stress spike for staggered cases is slightly higher than for the aligned ones. This spike also decreases with packing density for this set of cases. This can be attributed to the reduced inhomogeneity in the mean velocity field as the spacing between cubes reduces. Comparing the zero skewed surface (figure 17a) with the corresponding positively skewed (figure 17d) and negatively skewed (figure 17e) ones, it can be seen that the dispersive stress spike gets narrower and increases with decreasing $Sk$. The narrowing is an indirect consequence of the decrease in roughness height ($h_1/b$) as $Sk$ decreases. The Reynolds stresses, on the other hand, are similar irrespective of packing density or skewness.

Figure 17. Reynolds $\langle \overline {u_i'u_j'}\rangle ^+$ and dispersive $\langle \overline {u_i''u_j''}\rangle ^+$ stress comparisons for spanwise-aligned and staggered cases. The solid and dashed lines depict aligned and staggered arrangements, respectively. Cases: (a) AYL1H1 and SL1H1; (b) AYL2H1 and SL2H1; (c) AYL3H1 and SL3H1; (d) AYL1H2 and SL1H2; and (e) AYL1H3 and SL1H3. The vertical line indicates peak roughness height.

Similarly, plots comparing the stresses for spanwise-aligned (AY) and streamwise-aligned (AX) surfaces are shown in figure 18. The higher magnitudes of dispersive stress for AX surfaces are due to the presence of high momentum pathways between two rows of cubes and low momentum pathways over cube locations. Unlike in figure 17, the dispersive stress spike here increases with packing density for the AX surfaces. To further probe these observations, the streamwise mean velocity contours for AYL1H1, AYL2H1, AXL1H1 and AXL2H1 near the roughness peaks are shown in figures 19(a), 19(b), 19(c) and 19(d) respectively. Indeed, we see more inhomogeneity in the mean flow in figures 19(c) and 19(d) due to the overlapping wake regions and the high momentum pathways formed in between the cubes. With respect to skewness, the increasing trend of dispersive stress with decreasing skewness for AX surfaces is similar to that observed in AY surfaces.

Figure 18. Reynolds $\langle \overline {u_i'u_j'}\rangle ^+$ and dispersive $\langle \overline {u_i''u_j''}\rangle ^+$ stress comparisons for spanwise-aligned and streamwise-aligned cases. The solid and dashed lines depict spanwise-aligned and streamwise-aligned arrangements, respectively, for (a) AYL1H1 and AXL1H1, (b) AYL2H1 and AXL2H1, (c) AYL3H1 and AXL3H1, (d) AYL1H2 and AXL1H2, and (e) AYL1H3 and AXL1H3. The vertical line indicates peak roughness height.

Figure 19. Streamwise mean velocity contours at $z/b \approx 1.0$ for cases (a) AYL1H1, (b) AYL2H1, (c) AXL1H1 and (d) AXL2H1.

The normal stress comparisons in figure 20 show slight variations when valleys are removed from AYL1H1, AYL2H1 and AYL3H1. Note that the NV surfaces have very different roughness statistics in $k_{rms}/k_a$, $Sk$, $Ku$ and $ES$, but very similar Reynolds and dispersive stresses (figure 20), mean velocity profiles (figure 11) and $z_0/b$ (figure 16c) with respect to the AY and S counterparts.

Figure 20. Reynolds $\langle \overline {u_i'u_j'}\rangle ^+$ and dispersive $\langle \overline {u_i''u_j''}\rangle ^+$ stress comparisons for spanwise-aligned cases with and without valleys: (a) AYL1H1/NVL1H1, (b) AYL2H1/NVL2H1 and (c) AYL3H1/NVL3H1. The solid and dashed lines depict spanwise-aligned and no-valleys cases, respectively. The vertical line indicates peak roughness height.

In figures 17, 18 and 20, it can be observed that the streamwise dispersive component is still measurable near the half-channel height. These effects hint at the presence of secondary flows within these domains. Figure 21 shows the mean streamwise averaged contours of $\bar {u}^+$ with secondary flow $\bar {v}^+$ and $\bar {w}^+$ for a few cases. One can observe small but non-zero secondary flow near the half-height, and noticeable counter-rotating vortical structures near the roughness peaks. In most cases, the secondary structures can be seen to penetrate beyond the roughness sublayer, which is in alignment with other studies such as Vanderwel & Ganapathisubramani (Reference Vanderwel and Ganapathisubramani2015).

Figure 21. Mean streamwise-averaged constant $x$ plane showcasing $\bar {u}^+$ contours with secondary flows $\bar {v}^+$ and $\bar {w}^+$ (denoted by vectors) within the domain in (a) AYL1H1, (b) SL1H1, (c) AXYL1H1 and (d) SXYL1H1.

4.1. The task of roughness modelling

In this subsection, we try to answer the following question: what is the nature of roughness modelling? We begin by comparing the $z_0$ magnitudes in table 3 to the estimates provided by the existing rough-wall correlations. Here, the correlations in Forooghi et al. (Reference Forooghi, Stroh, Magagnato, Jakirlić and Frohnapfel2017) and Flack et al. (Reference Flack, Schultz and Barros2020) are taken as examples. The two correlations are given in (4.1) and (4.2) respectively:

(4.1)$$\begin{gather} k_s/k_t = 1.07\times(1-{\rm e}^{{-}3.5 \, ES})(0.67Sk^2 + 0.93Sk +1.3), \end{gather}$$

(4.2)$$\begin{gather}k_s = \left\{\begin{array}{@{}ll@{}} 2.11 k_{rms}, & Sk = 0, \\[2pt] 2.48 k_{rms}(1+Sk)^{2.24}, & Sk > 0, \\[2pt] 2.73 k_{rms}(2+Sk)^{{-}0.45}, & Sk < 0, \end{array} \right. \end{gather}$$

where roughness statistics $k_{rms}$, $Sk$, $k_t$ and $ES$ are invoked as inputs to model $k_s$. We note that the effective roughness height ($z_0$) in table 3 is linked directly to $k_s$ by the expression

(4.3)\begin{equation} k_s = z_0\,{\rm e}^{(\kappa A)}, \end{equation}

where the von Kármán constant is $\kappa = 0.4$, and $A = 8.5$.

Figure 22 shows the predicted $k_s$ for the rough surfaces in this study using the aforementioned correlations. We see that the surfaces in the present study do not fit into these existing correlations. Further, we see from figure 22(a) that surfaces with different arrangements (AX, AY, S), although having very different measured $k_s$, produce the same estimates according to (4.1) as they have the same $k_t$, $Sk$ and $ES$ (as indicated by the dashed lines). Similarly, in figure 22(b), the surfaces AXY and SXY, although having very different measured $k_s$, join the other arrangements with identical $k_{rms}$ and $Sk$, resulting in identical $k_s$ predictions according to (4.2). These observations indicate the inherent lack of roughness statistics that distinguish roughness surfaces in these correlations. Although the discussion here is limited to the correlations in Forooghi et al. (Reference Forooghi, Stroh, Magagnato, Jakirlić and Frohnapfel2017) and Flack et al. (Reference Flack, Schultz and Barros2020), the same inadequacies are prevalent in other empirical correlations as well, as observed by Yang et al. (Reference Yang, Zhang, Yuan and Kunz2023).

Figure 22. Existing correlations by (a) Forooghi et al. (Reference Forooghi, Stroh, Magagnato, Jakirlić and Frohnapfel2017) and (b) Flack et al. (Reference Flack, Schultz and Barros2020), applied to surface parameters obtained from the current study. The plots show the predicted equivalent sand-grain roughness height $k_s$ normalized by cube width $b$.

A remedy, as advocated and explained in Chung et al. (Reference Chung, Hutchins, Schultz and Flack2021) and other references cited therein, is to expand the list of inputs and include other roughness statistics. Here, statistics that help to distinguish different arrangements of elements are needed. Among other statistics, two-point correlations are the simplest ones for this purpose. The conventional definition of correlation length measures the distance at which the autocorrelation, i.e.

(4.4)$$\begin{gather} C_x(\Delta x) = \frac{1}{L_x L_y k_{rms}^2}\int_{L_y} \int_{L_x} (k(x, y) - \bar{k})(k(x + \Delta x, y) - \bar{k})\, {{\rm d}\kern0.07em x}\,{{\rm d}y}, \end{gather}$$

(4.5)$$\begin{gather}C_y(\Delta y) = \frac{1}{L_x L_y k_{rms}^2}\int_{L_y} \int_{L_x} (k(x, y) - \bar{k})(k(x, y + \Delta y) - \bar{k}) \,{{\rm d}\kern0.07em x}\,{{\rm d}y}, \end{gather}$$

drops to a certain value, typically $1/{\rm e}$. Here, $k(x,y)$ contains the height information for the rough surface at location $(x, y)$. In other words, the correlation lengths $Rl_x$ and $Rl_y$ are such that

(4.6)$$\begin{gather} C_x(Rl_x) = 1/{\rm e}, \end{gather}$$

(4.7)$$\begin{gather}C_y(Rl_y) = 1/{\rm e}. \end{gather}$$

The above definitions are suitable for irregular surfaces, but when applied to regular (e.g. cuboidal) roughness, thus-defined correlations simply give the roughness element width. Intuitively, we need a distance measure that is representative of the inter-repeating-tile length. This can be made possible by choosing the peak to peak distance in the correlation plot. Figure 23 demonstrates these definitions, with an irregular surface A0020 from Forooghi et al. (Reference Forooghi, Stroh, Magagnato, Jakirlić and Frohnapfel2017), and a regular surface R25 from Xu et al. (Reference Xu, Altland, Yang and Kunz2021) as examples. Extending them to both streamwise and spanwise directions, one can distinguish between staggered and aligned, as well as surfaces with different arrangements.

Figure 23. An illustration of the definition of correlation length for (a) an irregular surface and (b) a regular surface.

Having defined the correlation lengths, it is now possible to distinguish between the rough surfaces involved in this study. We demonstrate this in figure 24, where we look at 27 surfaces from our DNS study, namely AYL$i$H$j$, SL$i$H$j$ and AXL$i$H$j$, with $i, j = 1,2,3$. The number of input roughness statistics is progressively increased so that figure 24(a) contains $k_{rms}/k_a$ and $Sk$, and figure 24(b) contains $k_{rms}/k_a$, $Sk$ and $Rl_x/Rl_y$. We observe that the data points overlap in figure 24(a) as AY/AX/SL$i$H$j$ occupy the same location in the space, but in figure 24(b), the rough surfaces do not overlap. Now, if one were to build a rough-wall model using only $Sk$ and $k_{rms}$ as its inputs, then one would be seeking different outputs from the same input, which would inevitably result in errors, as observed in figure 22. However, if $Sk$, $k_{rms}$ and $Rl_x/Rl_y$ are invoked as inputs, then one would be seeking a single-valued function that maps from the input roughness statistics space to the output $k_s$ space. Following this line of thinking, we argue that the task of roughness modelling can be reduced to identifying a set of roughness statistics that distinguishes the rough surfaces in question.

Figure 24. Rough surfaces in two different parameter spaces: (a) $k_{rms}/k_a$ and $Sk$; (b) $k_{rms}/k_a$, $Sk$ and $Rl_x/Rl_y$.

This conclusion has many implications for rough-wall modelling. First, considering that there is no finite set of statistics that allows unique identification of all rough surfaces, i.e. two-dimensional functions $k(x,y)$, constructing a universal roughness-statistics-based rough-wall model could be a very difficult problem to solve. Consequently, all roughness-statistics-based rough-wall models are bound to have predictive power for surfaces that are closely registered around its calibration dataset only. Second, roughness modelling must consider the calibration dataset, and a universal list of roughness statistics for all rough-wall models is highly unlikely. One can argue that certain roughness types are more important than others, therefore certain roughness statistics are more important than others, but such judgements have to be subjective.

4.2. Selecting roughness statistics

Considering that all rough-wall models will have their applicable range, selecting inputs to the rough-wall models is a highly relevant issue. In this subsection, we discuss how one might go about selecting the roughness statistics that distinguish the rough surfaces in a given group of rough surfaces.

Without loss of generality, we limit the discussion to the seven roughness statistics in table 3. Figure 25 is a representation of the 36 rough surfaces in the seven-dimensional roughness statistics space represented by $k_{rms}/k_a$, $Sk$, $ES$, $Rl_x/k_a$, $Rl_y/k_a$, $Ku$ and $k_t/k_a$. In the figure, every roughness statistic corresponds to a vertical axis, and every rough surface corresponds to a line in the figure. We will make use of this plot to select the roughness statistics that are most relevant to the modelling of the rough surfaces in this study. We can, of course, use all roughness statistics. The objective, however, is to identify as few roughness statistics as possible such that the identified roughness statistics distinguish the rough surfaces in question. It is also worth noting that for random roughness, $k_t$ might not be a suitable parameter to represent the surface as it tends to increase as the sampling size gets larger, instead of converging.

Figure 25. High-dimensional parallel plots of the surface parameter space for (a) 36 surfaces in the current study, (b) $ES$ parameter space, and (c) $Ku$ parameter space. The values shown are scaled by the standard deviation.

First, we make the following observation. For a roughness statistic $s$ that represents a vertical axis in figure 25, rough surfaces that do not overlap in the axis can readily be distinguished by that statistic. For example, consider two surfaces AYL1H1 and SL1H1. All features excluding $Rl_y/k_a$ would remain the same for both these surfaces, hence this parameter would serve as the distinguishing feature. However, depending on how accurately one could measure the features, the features need not necessarily be identical for them to be overlapping. We may mathematically define ‘overlap’ as $|s_i-s_j|< c\,{\rm STD}_s$. That is, if $s_i$ and $s_j$ are close, then they are considered indistinguishable or overlap. Here, $s_i$ is the roughness statistic of the $i$th rough surface, STD$_s$ is the standard deviation of statistic $s$ in the rough surfaces, and $c$ is a constant. Again, the choice of $c$ will depend on how well one can measure the roughness statistic in question: if one can measure the statistic with precision $0.01$STD$_s$, then $c$ can take value 0.02, i.e. $+0.01-(-0.01)$. Here, we set $c=0.1$, which should be a very conservative choice. Following this line of thinking, a roughness statistic that distinguishes more rough surfaces is more relevant to the modelling of the rough surfaces in question. Take the statistics and the rough surfaces in figure 25(a) as our illustrative example. Kurtosis is not an extremely relevant roughness statistic for the modelling of these rough surfaces as it distinguishes only five groups of rough surfaces (see figure 25c). On the other hand, from figure 25(b) it can be observed that ES is a more relevant roughness statistic as it distinguishes around 12 groups of rough surfaces. In other words, the more the spread in the feature space, the higher its differentiating ability. Note that surfaces indistinguishable by a parameter may be indistinguishable in terms of drag. In such cases, the model would remain accurate even if it does not distinguish the two surfaces.

Following the discussion above, the process of selecting roughness statistics that distinguish rough surfaces in a given group of rough surfaces can follow a recursive greedy algorithm. To illustrate the algorithm, we take our current study as an example. For this group of rough surfaces, we have 36 surfaces. We will represent each rough surface by a number, and the 36 rough surfaces are numbered from 1 to 36. First, we will split each feature into bins of size 0.1 STD, and the number of surfaces present in each bin is counted. The feature with the greatest number of non-zero bins is the most distinguishing feature. For the 36 rough surfaces, $Rl_x$ is found to be that feature. If a surface is the only surface in one of the $Rl_x$ bins, then that surface is readily distinguishable by $Rl_x$. However, if there are a couple of rough surfaces in one $Rl_x$ bin, then these surfaces, which we call a cluster, cannot be distinguished by $Rl_x$, and we must invoke another roughness statistic to distinguish them. From figure 26, it can be seen that there are 10 such clusters, with one cluster containing 6 surfaces, three clusters containing 4 surfaces each, one cluster containing 3 surfaces, and five clusters containing 2 surfaces each. The purpose of the next stage is to find the next feature that best distinguishes these surfaces. A similar procedure is repeated for this stage for each cluster, and $ES$ is identified as the next most distinguishing feature. The process continues until either no clusters emerge or all features have been selected. In this case, this stopping criterion is attained after the third roughness statistics, $Rl_y$, is selected. Figure 26 illustrates the procedure.

Figure 26. Feature selection procedure illustration using rough surface indices. Each number here is representative of a rough surface: I denotes the initial stage involving all surfaces in the group; II denotes the next stage after one feature $Rl_x/k_a$ was selected; and III denotes the third stage when $ES$ is selected.

We can apply this recursive greedy algorithm to other groups of surfaces. Table 4 shows some of these examples, where this feature selection procedure is applied to the rough surfaces from Jouybari et al. (Reference Jouybari, Yuan, Brereton and Murillo2021), Womack et al. (Reference Womack, Volino, Meneveau and Schultz2022), Xu et al. (Reference Xu, Altland, Yang and Kunz2021) and the current study. Although we do not consider all surface geometries and other roughness parameters such as spanwise effective slope $ES_z$ and surface porosity $P_0$ mentioned in Jouybari et al. (Reference Jouybari, Yuan, Brereton and Murillo2021), it is interesting to note that we obtain one variable as the distinguishing feature within each of the three pairs $(E_x, E_z)$, $(P_0, Sk)$ and $(k_{rms}, Ku)$ that contain strong correlations in their work, and the algorithm does not pick up two strongly correlated features. In the six rough surfaces in Xu et al. (Reference Xu, Altland, Yang and Kunz2021), where the primary varying parameter is the planar packing density of the cubes, just one parameter, $k_{rms}/k_a$ in this case, seems to be sufficient. An interesting observation is that $k_{rms}/k_a$ is not always a distinguishing roughness statistic. For the surfaces with multiple roughness arrangements, such as in our study, the dimensionless two-point correlation lengths $Rl_x/k_a$ and $Rl_y/k_a$ emerge as having better distinguishing capability.

Table 4. Features selected for different groups.

A scenario might occur where this feature selection methodology, when extended to a larger set of surfaces, would exhaust all features. In this case, certain surfaces are indistinguishable from the available roughness statistics, yet they have different $k_s$. For example, figure 27 depicts two clusters of such surfaces from Forooghi et al. (Reference Forooghi, Stroh, Magagnato, Jakirlić and Frohnapfel2017) and Yang et al. (Reference Yang, Stroh, Chung and Forooghi2022), taken from a larger set of about 143 rough surfaces in the aforementioned roughness database. The three rough surfaces in figures 27(a)–27(c) from Forooghi et al. (Reference Forooghi, Stroh, Magagnato, Jakirlić and Frohnapfel2017) are referred to as A3588, A7088 and A1588. These are indistinguishable from the seven roughness statistics in table 3. The $k_s/k_a$ magnitudes in surfaces A1588, A3588 and A7088 from Forooghi et al. (Reference Forooghi, Stroh, Magagnato, Jakirlić and Frohnapfel2017) stand at approximately 11.1, 9.98 and 8.86, respectively, which are different, whereas the deviations in the aforementioned seven statistics are less than 0.1 standard deviation (STD) for these surfaces. The same applies to surfaces N14 and N18 from Yang et al. (Reference Yang, Stroh, Chung and Forooghi2022) whose $k_s/k_a$ magnitudes are close to 5.15 and 4.57, respectively.

Figure 27. Examples of surfaces with very similar statistics when considering the seven roughness statistics in table 3. The surfaces (a) A1588, (b) A3588 and (c) A7088 from Forooghi et al. (Reference Forooghi, Stroh, Magagnato, Jakirlić and Frohnapfel2017) belong to one group. The surfaces (d) N14 and (e) N18 from Yang et al. (Reference Yang, Stroh, Chung and Forooghi2022) belong to another group.

4.3. Building a rough-wall model

Having identified the important roughness statistics required to build a rough-wall model, this subsection discusses the next step involved, which is developing rough-wall models themselves. This is quite straightforward, and there are a number of approaches. Here, we present two such approaches: a multi-variate polynomial regression (MPR), which gives interpretable explicit expressions but has less descriptive power, and a feedforward neural network (FNN) based approach, which has great descriptive power but is a black box.

For the MPR model, the following steps were used to generate the polynomial terms: (i) minmax scaling and polynomial feature generation; (ii) determination of the degree and the polynomial features involved; and (iii) linear regression of the polynomial features. The regression follows an iterative bagging-based approach. At every iteration, the training dataset is randomly selected maintaining a test–train split of 0.7, and $m$ random polynomial terms are selected to develop a regression fit. Ridge regression is employed here, and the mean-squared test and training errors are computed. This step is repeated iteratively for different polynomials with the selected degree $d$ and number $m$ of polynomial features until we find a fit where the minimum in both training and test datasets occurs. These fits for the groups Jouybari et al. (Reference Jouybari, Yuan, Brereton and Murillo2021), Womack et al. (Reference Womack, Volino, Meneveau and Schultz2022), Xu et al. (Reference Xu, Altland, Yang and Kunz2021) and the current study are given by expressions (4.8), (4.9), (4.10) and (4.11), respectively (the reader is directed to Kleinbaum et al. (Reference Kleinbaum, Kupper, Nizam and Rosenberg2013) for further details of MPR):

(4.8)\begin{align} k_s/k_a &={-}0.53(k_{rms}/k_a)^3 -0.26\,Sk^3 +0.69\,ES^3 +0.78(k_{rms}/k_a)^2\,Sk + 2.22\,Sk^2\,ES \nonumber\\ &\quad + 2.45(k_{rms}/k_a)^2\,ES +1.44(k_{rms}/k_a)\,Sk^2 -0.85\,Sk\,ES^2 \nonumber\\ &\quad +2.86(k_{rms}/k_a)\,Sk \,ES -0.90(k_{rms}/k_a)^2 -1.51\,ES^2 + 0.5\,Sk\,ES \nonumber\\ &\quad -0.44(k_{rms}/k_a)\,Sk +0.65\,ES -0.12k_{rms}/k_a, \end{align}

(4.9)\begin{align} k_s/k_a &= 0.95\,Sk^2 + 0.77\,ES\,Rl_x/k_a -2.63\,ES +1.79\,Sk -2.58k_{rms}/k_a, \end{align}

(4.10)\begin{align} k_s/k_a &= 0.88(k_{rms}/k_a)^2 + 0.13k_{rms}/k_a, \end{align}

(4.11)\begin{align} k_s/k_a &= 0.03\,{ES}^2 -1.19({Rl}_x/k_a)^2 + 1.00({Rl}_y/k_a)^2 -0.08({Rl}_x/k_a)({Rl}_y/k_a)\nonumber\\ &\quad + 0.06\,ES + 1.71\,Rl_x/k_a -1.41\,Rl_y/k_a. \end{align}

Figure 28 shows the true versus predicted plots for the various groups, with $R^2$ standing for the coefficient of determination, which is a measure of accuracy. The $R^2$ values shown in figure 28 have been calculated from the training data, holding values 0.98, 0.98, 0.99 and 0.85 for the respective groups. The corresponding $R^2$ values based on the test data stand as 0.77, 0.99, 0.99 and 0.64. It is possible to choose polynomial forms of higher $R^2$ on the training data, but we may risk overfitting. Hence these models are adopted, maintaining a balance between generalization and accuracy. It is further noted that the errors shown in $k_s$ predictions here would generate much smaller errors in drag prediction, as the skin-friction coefficient $C_f$ is proportional to $[\log (k_s^+)]^2$.

Figure 28. The MPR-based regression for groups: (a) Jouybari et al. (Reference Jouybari, Yuan, Brereton and Murillo2021), (b) Womack et al. (Reference Womack, Volino, Meneveau and Schultz2022), (c) Xu et al. (Reference Xu, Altland, Yang and Kunz2021), and (d) current study. Here, $R^2$ stands for the coefficient of determination and is calculated on the training dataset. The red colour stands for test data, and the blue stands for training data. The $k_s/k_a$ values are normalised using a minmax scaler.

For the FNN model, a standard three-layer network has been trained with 10 neurons in each layer. The size of the neural network is small, but as we will see, it is sufficient. Figure 29 shows the parity plots for these networks. Here, each FNN utilizes 100 % of its group for its training data. This is done to show the perfect fitting of the data in each case, which demonstrates the effectiveness of our feature selection strategy and the strong descriptive power of FNNs.

Figure 29. The FNN-based regression for groups: (a) Jouybari et al. (Reference Jouybari, Yuan, Brereton and Murillo2021), (b) Womack et al. (Reference Womack, Volino, Meneveau and Schultz2022), (c) Xu et al. (Reference Xu, Altland, Yang and Kunz2021), and (d) current study.