2 Computation details

2.1 Large-eddy simulations

We use LES for flow over sparsely packed cubes in a half-channel. LES resolves the relatively energetic motions and models the relatively less energetic motions and is therefore a very cost-effective computational tool (Meneveau & Katz Reference Meneveau and Katz2000). The LES code we use is LESGO. It solves the filtered incompressible Navier–Stokes equations at nominally infinite Reynolds number. The convective terms are evaluated in its rotational form. A pseudo-spectral discretization is used in both the streamwise and spanwise directions. A second-order finite difference scheme is used in the wall-normal direction. Conserved flow quantities, including the streamwise, wall-normal and spanwise velocities, are solved on a staggered grid (as opposed to a collocated one). The wall-mounted cubes are resolved using an immersed boundary method. The cube surfaces are such that they coincide with the surfaces of the pressure computational cells. The deviatoric part of the subgrid-scale stress is modelled via the scale-dependent Lagrangian dynamic model (Bou-Zeid, Meneveau & Parlange Reference Bou-Zeid, Meneveau and Parlange2005). The code has been extensively used for similar flow problems (see e.g. Cheng & Porté-Agel (Reference Cheng and Porté-Agel2015), Yang et al. (Reference Yang, Sadique, Mittal and Meneveau2016) and Giometto et al. (Reference Giometto, Lozano-Duran, Park and Moin2017)). The reader is directed to Yang & Abkar (Reference Yang and Abkar2018) for further details of the solver.

Figure 2 shows the computational domain. The height of the wall-mounted cubes is $h$ . The domain size is $L_{x}\times L_{y}\times L_{z}$ in the streamwise, spanwise and wall-normal directions, respectively. The half-channel height is $L_{z}=3.5h$ for all the cases except for A100-h, where the half-channel height is $L_{z}=7h$ . The flow is driven by a constant volumetric body force $f_{b}$ . The friction velocity is $u_{\unicode[STIX]{x1D70F}}=\sqrt{f_{b}L_{z}/\unicode[STIX]{x1D70C}}$ and is kept constant. We vary $L_{x}$ and $L_{y}$ . The streamwise and spanwise extents of the computational domain are such that they are never smaller than approximately $2\unicode[STIX]{x03C0}L_{z}$ . According to Lozano-Durán & Jiménez (Reference Lozano-Durán and Jiménez2014), our computational domain is sufficiently large. It is worth noting that domains that are much smaller were used to study rough-wall boundary layer flows at moderate surface coverage densities (Chung et al. Reference Chung, Chan, MacDonald, Hutchins and Ooi2015; MacDonald et al. Reference MacDonald, Chung, Hutchins, Chan, Ooi and García-Mayoral2016, Reference MacDonald, Chung, Hutchins, Chan, Ooi and García-Mayoral2017).

Figure 2.

A sketch of the computational domain. Periodicity is imposed in both the streamwise ( $x$ ) and spanwise ( $y$ ) directions. The cubes are of height $h$ . The size of the computational domain is $L_{x}\times L_{y}\times L_{z}$ .

Table 1 summarizes all the LES cases. By varying $L_{x}$ , $L_{z}$ and the number of wall-mounted cubes in the computational domain, we consider rough walls with surface coverage densities between $0.081\,\%$ and $4.4\,\%$ . Two extensively studied roughness arrangements are considered, namely, the aligned and the staggered arrangements. Figure 3(a,b) shows the sketches of the two arrangements. For an aligned cube array, both the streamwise and the spanwise inter-cube distances are $l=h/\sqrt{\unicode[STIX]{x1D706}_{p}}$ , which varies from approximately $4h$ to $35h$ . The streamwise distance between two cubes is $2l$ when the cubes are staggered in the spanwise direction. If we project the wall-mounted cubes onto a spanwise–wall-normal plane, for a staggered arrangement, the spanwise distance between two neighbouring cubes is $l/2$ .

Table 1.

LES cases. The names of cases are Axxx or Sxxx, where ‘A’ and ‘S’ denote ‘aligned’ and ‘staggered’ arrangements, respectively, and the three digits ‘xxx’ indicate the surface coverage density. Case A100-h is similar to A100 but in a horizontally larger and vertically higher computational domain, i.e.  $L_{x}\times L_{y}\times L_{z}=40\times 40\times 7$ . Case A100-w is similar to A100 but in a horizontally larger computational domain, i.e.  $L_{x}\times L_{y}\times L_{z}=40\times 40\times 3.5$ . Cases A100-d and A440-d are similar to A100 and A440 but with the inflows specified using a precursor channel flow simulation. The domain sizes are shown in terms of the cube height $h$ . The ‘Mesh’ column shows the grid size in million (M) cells. The column ‘N’ shows the number of wall-mounted cubes in columns (streamwise) and rows (spanwise). Note that we keep two significant digits for all the numbers.

Figure 3.

Sketches of the two roughness arrangements: (a) aligned arrangement and (b) staggered arrangement. The flow goes from left to right.

The boundary conditions are as follows. A periodic boundary condition is imposed in both the streamwise and spanwise directions for all cases except for A100-d and A440-d, where the inflow is provided by a precursor channel flow simulation. A free-slip condition is used for the top boundary. The surface shear stresses on both the ground and the cube surfaces are specified through a wall model,

(2.1) $$\begin{eqnarray}\unicode[STIX]{x1D749}/\unicode[STIX]{x1D70C}=-\left[\frac{\unicode[STIX]{x1D705}}{\log (\unicode[STIX]{x1D6E5}/z_{o})}\right]^{2}\tilde{\boldsymbol{U}}_{\Vert }|\tilde{\boldsymbol{U}}_{\Vert }|,\end{eqnarray}$$

where $\unicode[STIX]{x1D749}$ is the shear stress vector, $\boldsymbol{U}_{\Vert }$ is the wall-parallel velocity vector at a distance $\unicode[STIX]{x1D6E5}$ from the wall, $\tilde{\cdot }$ is a spatial filtering operator (Bou-Zeid et al. Reference Bou-Zeid, Meneveau and Parlange2005; Yang, Park & Moin Reference Yang, Park and Moin2017b ), $z_{o}=10^{-4}h$ is a pre-specified roughness/viscous length scale, $\unicode[STIX]{x1D705}=0.4$ is the von Kármán constant (von Kármán Reference von Kármán1931; Marusic et al. Reference Marusic, Monty, Hultmark and Smits2013) and $\log$ is the natural logarithm. Equilibrium wall models including the one in (2.1) have difficulties in handling flow separations and predicting flow reattachments (Slotnick et al. Reference Slotnick, Khodadoust, Alonso, Darmofal, Gropp, Lurie and Mavriplis2014). This is less of an issue for sparsely packed roughness. For a cubical roughness element, the flow separates at the trailing edge irrespective of the wall model used. Also, because we study cube arrays that are very sparsely packed, the flow reattaches before the recirculation bubble reaches the downstream wall-mounted cube, which is clear from figure 4, where we show contours of the mean streamwise velocities on a streamwise–wall-normal plane through the centre of a column of wall-mounted cubes. In figure 4, two neighbouring cubes are $4.8h$ apart in the streamwise direction, which is the closest among the cases in table 1. Hence, it is probably safe to say that, at very low surface coverage densities, the drag forces on the wall-mounted cubes do not depend critically on near-wall modelling.

Figure 4.

Contours of the mean streamwise velocity on a streamwise–wall-normal plane through the centre of a column of wall-mounted cubes. We show results from A440. The contour levels cut off at $U=0$ to highlight the reattachment location. The flow reattached before the recirculation bubble reaches the next cube.

For all cases in table 1, we use a uniform grid, with $\unicode[STIX]{x0394}x=\unicode[STIX]{x0394}y=\unicode[STIX]{x0394}z=h/8$ . The same grid resolution has previously been used in a number of studies (Xie & Castro Reference Xie and Castro2006; Jiang et al. Reference Jiang, Jiang, Liu and Sun2008; Graham & Meneveau Reference Graham and Meneveau2012). The results are also validated by comparing with wind tunnel measurements and wall-resolved LES in appendices A and B.

2.2 Direct numerical simulation

The near-wall turbulence is modelled instead of resolved in LES. To get more reliable measurements at the wall and to study the behaviour of the skin friction, we resort to DNS. The computational set-up is the same as in the case A100. We refer to the DNS as A100-d. A structured, body-fitted mesh is used. The size of the mesh is approximately 9.6 million. We refine the mesh near the wall-mounted cubes. Figure 5(a,b) shows the mesh near a wall-mounted cube. The grid resolution is such that the local wall-normal resolution is nowhere coarser than $\unicode[STIX]{x1D6E5}_{n}^{+}=0.5$ and the local wall-parallel resolution is $\unicode[STIX]{x1D6E5}_{\Vert }^{+}\approx 5$ . Compared with resolutions used for similar flow problems (Coceal, Dobre & Thomas Reference Coceal, Dobre and Thomas2007a ; Balakumar, Park & Pierce Reference Balakumar, Park and Pierce2014), our grid resolution is quite high. Similar to A100, the flow is driven by a constant body force. The friction Reynolds number is $Re_{\unicode[STIX]{x1D70F}}=u_{\unicode[STIX]{x1D70F}}L_{z}/\unicode[STIX]{x1D708}=350$ . Case A100 is at a nominally infinite Reynolds number. Considering the difference in their Reynolds numbers, the drag and the friction in A100 and A100-d are not expected to be the same, although the LES and the DNS solutions are qualitatively alike. For the DNS, we use the in-house unstructured finite-volume flow solver CharLES (Khalighi et al. Reference Khalighi, Ham, Nichols, Lele and Moin2011; Bermejo-Moreno et al. Reference Bermejo-Moreno, Campo, Larsson, Bodart, Helmer and Eaton2014). The code solves the full compressible Navier–Stokes equation. It uses a fourth-order-accurate spatial discretization (Mahesh, Constantinescu & Moin Reference Mahesh, Constantinescu and Moin2004) and a third-order-accurate temporal discretization (Williamson Reference Williamson1980; Moin Reference Moin2010). The Mach number in the computational domain is nowhere greater than 0.15, and therefore the flow may be considered as incompressible (Lele Reference Lele1994; Park & Moin Reference Park and Moin2014). The code is well validated and has been used extensively for boundary layer flow problems (Joo et al. Reference Joo, Medic, Philips and Bose2014; Larsson et al. Reference Larsson, Laurence, Bermejo-Moreno, Bodart, Karl and Vicquelin2015; Park Reference Park2017).

Figure 5.

Grid near a wall-mounted cube: (a) on a spanwise–wall-normal plane and (b) on a streamwise–wall-normal plane.

3 Drag forces and flow physics

It is the conclusion of the classic drag partition theory that $C_{d}$ ( ${<}C_{R}$ ) is a non-increasing function of $\unicode[STIX]{x1D706}_{p}$ . This is often the starting point of rough-wall drag-force modelling. Here, we show results that challenge the above classic theory. Figure 6(a) shows the drag coefficient $C_{d}$ as a function of the surface coverage density $\unicode[STIX]{x1D706}_{p}$ for both the aligned and the staggered arrays. For the staggered arrays, the drag coefficient $C_{d}$ is approximately $0.49$ and stays roughly constant for $\unicode[STIX]{x1D706}_{p}<0.25\,\%$ . It increases as a function of $\unicode[STIX]{x1D706}_{p}$ from $\unicode[STIX]{x1D706}_{p}=0.25\,\%$ to $\unicode[STIX]{x1D706}_{p}\approx 5\,\%$ . For the aligned arrays, the drag coefficient increases as a function of $\unicode[STIX]{x1D706}_{p}$ from $\unicode[STIX]{x1D706}_{p}=0.08\,\%$ to approximately $\unicode[STIX]{x1D706}_{p}=1\,\%$ , and then decreases as a function of $\unicode[STIX]{x1D706}_{p}$ for $\unicode[STIX]{x1D706}_{p}>1.5\,\%$ . Figure 6(b) shows the friction coefficient $C_{f}$ as a function of the surface coverage density. For both arrangements, the skin-friction coefficient increases as a function of $\unicode[STIX]{x1D706}_{p}$ for $\unicode[STIX]{x1D706}_{p}<1\,\%$ . Note that $C_{d}$ is the drag coefficient of a single cube, not the skin-friction coefficient of the entire rough wall, the latter of which is an increasing function of the surface coverage density until approximately $\unicode[STIX]{x1D706}_{p}\approx 15\,\%{-}30\,\%$ (Jiménez Reference Jiménez2004).

Figure 6.

(a) Drag coefficient $C_{d}=f/(\unicode[STIX]{x1D70C}U_{h}^{2}h^{2})$ as a function of the surface coverage density $\unicode[STIX]{x1D706}_{p}$ for both the aligned and the staggered cube arrays. Here, $f$ is the drag force on one wall-mounted cube; and $U_{h}$ is the mean velocity at the cube height. The red circles show the results of the staggered cube arrays, and the blue squares show the results of the aligned cube arrays. (b) Friction coefficient $C_{f}=\unicode[STIX]{x1D70F}_{s}/(\unicode[STIX]{x1D70C}U_{h}^{2})$ as a function of the surface coverage density $\unicode[STIX]{x1D706}_{p}$ . Here, $\unicode[STIX]{x1D70F}_{s}$ is the friction drag per unit planar area. The symbols are the same as in panel (a).

Figure 7.

Contours of the streamwise velocity: results for aligned cases (a) A008, (b) A100 and (c) A440. The vectors show the in-plane motions (spanwise and wall-normal velocities). Both the spanwise and the wall-normal velocities are averaged in time and in the streamwise direction. The vectors are shown every four grid points in both directions in (a,b), and every other point in (c) for better presentation. The cubes are projected onto the plane. The flow is periodic in the spanwise direction. We show only one repeating flow unit.

There are a few trends that are expected from the classic drag partition theory and a number of trends that are less expected. First, we briefly discuss the trends that are expected from the theory. The decrease of the drag coefficient and the friction coefficient as a function of $\unicode[STIX]{x1D706}_{p}$ at relatively high coverage densities (although only a few per cent) is expected and is simply a result of flow sheltering (Yang et al. Reference Yang, Sadique, Mittal and Meneveau2016). For the staggered cubes, at very low surface coverage densities, i.e.  $\unicode[STIX]{x1D706}_{p}<0.25\,\%$ , the cubes act as isolated elements, and the drag coefficient stays constant, $C_{d}=C_{R}=0.49$ . The drag coefficient $C_{R}=0.49$ is close to the measurements in the literature, although the exact value may vary depending on the flow conditions, including the Reynolds number and the boundary layer height (Akins, Peterka & Cermak Reference Akins, Peterka and Cermak1977; ESDU 1986; Macdonald, Griffiths & Hall Reference Macdonald, Griffiths and Hall1998). Last (and needless to say) the flow behaves differently above aligned and staggered cubes. Next, we discuss the trends that are not expected from the drag partition theory, i.e. the increase of the drag coefficient and the skin-friction coefficient as a function of the surface coverage density. Flow sheltering leads to reduced drag forces and cannot be responsible for a drag coefficient larger than $C_{R}$ . Here, the increase of the drag coefficient as a function of $\unicode[STIX]{x1D706}_{p}$ is a result of secondary turbulent flows. Secondary flows are common in, for example, square ducts and open channels (Nikuradse Reference Nikuradse1930; Nakagawa Reference Nakagawa2017). In the context of canonical boundary layer flows, secondary flows arise often as a result of spanwise-heterogeneous ground roughness (see e.g. Wang & Cheng Reference Wang and Cheng2005; Fishpool, Lardeau & Leschziner Reference Fishpool, Lardeau and Leschziner2009; Vermaas, Uijttewaal & Hoitink Reference Vermaas, Uijttewaal and Hoitink2011), and they manifest as streamwise-elongated counter-rotating vortices that extend from the ground to the top of the boundary layer (Willingham et al. Reference Willingham, Anderson, Christensen and Barros2014; Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015; Vanderwel & Ganapathisubramani Reference Vanderwel and Ganapathisubramani2015). Although our roughness is nominally homogeneous, the secondary flows in our LES are similar to the ones above spanwise-heterogeneous roughness.

Figures 7 and 8 show the streamwise velocity contours on a spanwise–wall-normal plane. The streamwise velocity is averaged both in time and in the streamwise direction. For brevity, we show results for three representative surface coverage densities, i.e.  $\unicode[STIX]{x1D706}_{p}=0.08\,\%$ , $1.0\,\%$ and $4.4\,\%$ . The overlying vectors show the in-plane motions. The spanwise velocity and the wall-normal velocity are also averaged in time and in the streamwise direction. We see secondary flows. A pair of counter-rotating vortices is found on the two sides of a wall-mounted cube. For a given surface coverage density, there are twice as many windward projected cubes in staggered cases than in aligned cases; as a result, there are twice as many vortex pairs in staggered cases than in aligned cases. The counter-rotating vortices extend from the ground to the top of the boundary layer in A/S008 and A/S100, i.e. for cubes that are sparsely packed. At cube locations, the counter-rotating vortices bring high-speed fluid from the upper part of the boundary layer to the lower part, and on the two sides of the cubes, low-speed fluid in the lower part of the boundary layer is brought to the upper part. This leads to the spanwise heterogeneity in the streamwise velocity in figures 7 and 8, featuring high-speed flows above the cubes and low-speed flows on the two sides of the cubes. Because of the spanwise heterogeneity, $U_{h}$ is an underestimate of the velocity of the incoming flows to the cubes, leading to a drag coefficient $C_{d}$ that is $C_{d}>C_{R}$ .

Figure 8.

Same as figure 7, but for staggered cases (a) S008, (b) S100 and (c) S440.

The wall-normal extents of the counter-rotating vortices decrease as the inter-cube distance decreases (see figures 7 c and 8 c). For A/S440, the outer flows are nearly homogeneous. Reynolds et al. (Reference Reynolds, Hayden, Castro and Robins2007) found spanwise heterogeneity in the streamwise velocity above cube arrays of a surface coverage density $\unicode[STIX]{x1D706}_{p}=25\,\%$ . Different from the secondary flows in figures 7 and 8, whose spanwise extents are confined by two columns of cubes, the high/low-momentum pathways in Reynolds et al. (Reference Reynolds, Hayden, Castro and Robins2007) span approximately $4h$ in the spanwise direction and cover approximately two columns of cubes. As a result, the flow experiences a larger drag force at some cube locations than at others on average. This violates ergodicity. Considering that the secondary flows above spanwise-heterogeneous roughness and the secondary flows in our LES do not violate ergodicity, we think the secondary flows reported in Reynolds et al. (Reference Reynolds, Hayden, Castro and Robins2007) are probably different flow features than the ones considered here.

Next, we briefly discuss how we may go about modelling the drag coefficient at low coverage densities. Figure 9 shows a sketch of the trailing vortices and the counter-rotating vortices (secondary flows). The trailing vortices join the counter-rotating vortices and strengthen them. We define $s_{x}$ to be the streamwise inter-cube distance. (Per figure 3, for the aligned arrangement, $s_{x}=l=h/\sqrt{\unicode[STIX]{x1D706}_{p}}$ ; and for the staggered arrangement, $s_{x}=2l=2h/\sqrt{\unicode[STIX]{x1D706}_{p}}$ .) If we fix $s_{x}$ , the strength of the counter-rotating vortices is a function of $s_{y}/L_{z}$ only, where $s_{y}$ is the spanwise distance between two pairs of counter-rotating vortices, and $L_{z}$ is the half-channel height (Anderson et al. Reference Anderson, Yang, Shrestha and Awasthi2018; Yang & Anderson Reference Yang and Anderson2018). It follows that, given $s_{x}$ , then $C_{d}-C_{R}$ is a function of $s_{y}/L_{z}$ only. A possible ansatz is

(3.1) $$\begin{eqnarray}C_{d}/C_{R}-1=(s_{y}/L_{z})^{-\unicode[STIX]{x1D6FC}}g(s_{x}/h),\end{eqnarray}$$

where $g$ is a generic function and $\unicode[STIX]{x1D6FC}>0$ is a model parameter. Figure 10 shows $(C_{d}/C_{R}-1)(s_{y}/L_{z})^{\unicode[STIX]{x1D6FC}}$ for $\unicode[STIX]{x1D6FC}=0.4$ . The data collapse. For large $s_{x}$ values, the cubes are non-interacting and $C_{d}=C_{R}$ . For small $s_{x}$ values, flow sheltering dominates, and $C_{d}$ decreases as $s_{x}$ decreases (as $\unicode[STIX]{x1D706}_{p}$ increases). For $10<s_{x}/h<35$ , $C_{d}-C_{R}$ increases as $s_{x}$ decreases, and $C_{d}/C_{R}-1=(s_{y}/L_{z})^{-\unicode[STIX]{x1D6FC}}g(s_{x}/h)$ .

Figure 9.

A sketch of the relevant flow structures. Flow is from left to right. For clarity, we have sketched flow structures on only one side of the cubes.

Figure 10.

Plot of $(C_{d}/C_{R}-1)(s_{y}/L_{z})^{\unicode[STIX]{x1D6FC}}$ as a function of $s_{x}/h$ ; $\unicode[STIX]{x1D6FC}=0.4$ is a parameter.

Next, we take a brief look at the friction on the ground. Compared with the drag forces on cubes, the friction on the ground is often small at high Reynolds numbers and is not the focus of this work. The near-wall turbulence is modelled instead of resolved in our LES; hence the wall shear stresses in our LES are less reliable (Yang, Bose & Moin Reference Yang, Bose and Moin2017a ). To study the wall shear stresses, we use DNS A100-d. The DNS computes flow in a half-channel above an array of aligned cubes with a surface coverage density of $1\,\%$ . Figure 11(a) shows the mean streamwise friction on the ground. Recirculation behind each wall-mounted cube leads to a region of reduced skin friction underneath the near wake, as expected. Underneath the far wake, we see two bands of high-stress regions. This is not expected from the classic drag partition theory, but is consistent with figure 9 and is responsible for the increase of the friction coefficient as a function of the surface coverage density. To the best of our knowledge, this is the first time a high-stress region has been found downstream of a cubical roughness.

Figure 11.

(a) Contours of the temporally averaged streamwise skin friction in case A100-d. (b) A zoomed view of the region underneath the far wake.

To briefly summarize: we find secondary flows on homogeneously placed cubical roughness elements; the secondary flows manifest as counter-rotating vortices; as a result of the secondary flows, the spatially and temporally averaged wind speed at the cube height underestimates the velocity of the incoming flow to the cubes, and is more so when the secondary flows become stronger; as a result, the drag coefficient $C_{d}$ takes values $C_{d}>C_{R}$ . Detailed analysis shows that, for aligned and staggered cube arrays, the increase in the drag coefficient scales as $C_{d}/C_{R}-1=g(s_{x}/h)(s_{y}/L_{z})^{-\unicode[STIX]{x1D6FC}}$ .

4 Further discussion on secondary flows

Turbulent secondary flows in rough-wall boundary layers have received a lot of attention in the recent literature. They are found above spanwise-heterogeneous roughness (see e.g. Mejia-Alvarez & Christensen Reference Mejia-Alvarez and Christensen2013; Barros & Christensen Reference Barros and Christensen2014) and spanwise-heterogeneous superhydrophobic surfaces (Jelly, Jung & Zaki Reference Jelly, Jung and Zaki2014; Lee, Jelly & Zaki Reference Lee, Jelly and Zaki2015). Spanwise heterogeneity is often imposed by placing streamwise strips of closely packed roughness elements (see e.g. Nugroho, Hutchins & Monty Reference Nugroho, Hutchins and Monty2013; Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015; Vanderwel & Ganapathisubramani Reference Vanderwel and Ganapathisubramani2015; Yang & Anderson Reference Yang and Anderson2018), and the spanwise-heterogeneous surface roughness changes the flow dynamics in the outer layer. High-momentum pathways (HMPs) form above roughness strips, and low-momentum pathways (LMPs) form between two strips (Mejia-Alvarez & Christensen Reference Mejia-Alvarez and Christensen2010). The nomenclature, i.e. HMPs and LMPs, is meant to distinguish the secondary flows from the high-speed and low-speed regions, which are transient in nature (see the detailed discussions in Meinhart & Adrian (Reference Meinhart and Adrian1995), Hutchins & Marusic (Reference Hutchins and Marusic2007) and de Silva, Hutchins & Marusic (Reference de Silva, Hutchins and Marusic2016)). Although the surface roughness in this work is nominally homogeneous, the secondary flows in our LES are very similar to the ones above spanwise-heterogeneous roughness. Figures 12 and 13 show the time-averaged streamwise velocity on a streamwise–spanwise plane at a distance $3h$ from the wall. Except for S440 (in figure 13 c), LMPs and HMPs are found in all the cases. HMPs are found at cube locations, and LMPs are found on the sides of the cube locations. The HMPs and the LMPs are slightly compressed in the spanwise direction at locations downstream of the cubes, leading to weak streamwise heterogeneity. It is worth noting that the HMPs and LMPs could not be identified in instantaneous flow fields (see figures 14 and 15, where we show the contours of the instantaneous streamwise velocities on a spanwise–wall-normal plane at an equal distance from two rows of roughness elements).

Figure 12.

Contours of the time-averaged streamwise velocity on a streamwise–spanwise plane located at a distance $3h$ from the wall for aligned cases (a) A008, (b) A100 and (c) A440. The contour levels are kept the same among the plots.

Figure 13.

Same as figure 12, but for staggered cases (a) S008, (b) S100 and (c) S440.

Figure 14.

Contours of the instantaneous streamwise velocity on a spanwise–wall-normal plane located at an equal distance from two rows of cubes for aligned cases (a) A008, (b) A100 and (c) A440. The spanwise extent of the computational domain is 70. For a better presentation, we show only half of the computational domain.

Figure 15.

Same as figure 14, but for staggered cases (a) S008, (b) S100 and (c) S440.

The sizes of the secondary flows are confined by both the wall-mounted cubes and the half-channel height. At a relatively high surface coverage density, e.g. A440, the cubes dictate the sizes of the secondary flows; and at a low surface coverage density, e.g. A008, the half-channel height dictates the sizes of the secondary flows. Comparing A100, where the computational domain size is $L_{x}\times L_{y}\times L_{z}=20h\times 20h\times 3.5h$ , A100-w, where the computational domain size is $L_{x}\times L_{y}\times L_{z}=40h\times 40h\times 3.5h$ , and A100-h, where the computational domain size is $L_{x}\times L_{y}\times L_{z}=40h\times 40h\times 7h$ , we study the effect of the domain size at a somewhat intermediate surface coverage density. Figure 16(a) shows the streamwise- and time-averaged longitudinal velocity of A100-h. Comparing with figures 7(a) and 16(a), the centres of the secondary vortices are at the same wall-normal locations. In the absence of secondary flows, a roughness sublayer, which is usually $3h$ to $5h$ (Lee et al. Reference Lee, Sung and Krogstad2011), can be defined, and the flow is statistically homogeneous above the roughness sublayer. For A100-h, because the secondary vortices span the whole channel, the flow is nowhere statistically homogeneous. The measured drag coefficient is $C_{d}=0.52$ from A100-h, which agrees reasonably well with that measured from A100, i.e.  $C_{d}=0.53$ . Figure 16(b) shows the streamwise- and time-averaged longitudinal velocity of A100-w. Comparing A100 and A100-w, extending the computational domain in the streamwise and the spanwise directions does not seem to make any difference. The drag coefficient is $C_{d}=0.53$ from A100-w, and is exactly the same as that measured from A100. Last, figure 17 shows the mean velocity profiles of A100, A100-h and A100-w. The three profiles agree reasonably well. This is particularly the case for the flow in the roughness layer.

Figure 16.

Contours of the streamwise- and time-averaged longitudinal velocity of (a) A100-h and (b) A100-w. The arrows show the in-plane motion.

Figure 17.

Mean velocity profiles of A100, A100-h and A100-w. Here $U_{xy}$ is the streamwise-, spanwise- and time-averaged velocity.

Figure 18.

(a) Contours of time-averaged streamwise velocity at a wall-normal height $z=1.2h$ in case A100-d. The wall area near the outlet is not mounted with cubes and is not shown. The cube locations are indicated with thin solid lines. (b) Normalized drag forces at each row for case A100-d. (c) Same as panel (a) but for A440-d. (d) Same as panel (b) but for A440-d.

So far, we have considered flows that are fully developed, for which a periodic boundary condition is imposed in the streamwise direction. Next, we consider flows that are developing, where the inflow is provided by a precursor channel flow simulation. We examine if secondary flows form shortly downstream and if the drag forces increase as a result. The two cases we will consider are A100-d and A440-d. Their fully developed counterparts, i.e. A100 and A440, show strong and weak secondary flow effects, respectively. Figure 18(a) shows the time-averaged streamwise velocity at a wall-normal distance $z=1.2h$ . A few rows downstream from the inlet, HMPs form above the wall-mounted cubes and LMPs form in between. This is also clear from figure 19(a,c,e), where we show contours of the time-averaged streamwise velocity at various streamwise locations. Figure 18(b) shows the normalized drag force at each row. After a decrease in the normalized drag from the first row to the second due to flow sheltering (where the secondary vortices have not emerged yet), the drag force gradually increases as a result of secondary vortices, and at row 8, the drag is already above that of the first row. The results of A440-d are shown in figures 18(c,d), and 19(b,d,f). The flow field is much less heterogeneous in the spanwise direction, and the normalized drag increases only slightly from row 2 to row 17. Hence, it is probably safe to say that the secondary vortices we see in figures 7 and 8 are not artifacts of the periodic boundary condition.

Figure 19.

(a,c,e) Time-averaged streamwise velocity at $x=20h$ , $x=40h$ and $x=60h$ , respectively, from A100-d. (b,d,f) Time-averaged streamwise velocity at $x=21h$ , $x=40h$ and $x=59h$ , respectively, from A440-d. The streamwise locations are such that they are at the midpoint of two rows of roughness elements.

Having verified that the secondary vortices in our LES are not artifacts of periodic boundary conditions nor a confined computational domain, we take a closer look at these flow structures. A defining feature of secondary flows is that they carry a momentum flux. This could potentially pose a challenge to drag-force modelling. Although Raupach (Reference Raupach1992) made no assumption about the mean flow, neglecting dispersive stresses in the Reynolds-averaged streamwise momentum equation is a common practice for rough-wall drag-force modelling (Macdonald Reference Macdonald2000; Coceal & Belcher Reference Coceal and Belcher2004; Harman & Finnigan Reference Harman and Finnigan2007).

Here we derive the equation used in rough-wall modelling. The streamwise momentum equation reads

(4.1) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}t}+u\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}x}+v\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}y}+w\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}z}=-\frac{1}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x}+\unicode[STIX]{x1D708}\left(\frac{\unicode[STIX]{x2202}^{2}u}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}^{2}u}{\unicode[STIX]{x2202}y^{2}}+\frac{\unicode[STIX]{x2202}^{2}u}{\unicode[STIX]{x2202}z^{2}}\right)-f,\end{eqnarray}$$

where $f$ represents the drag force. Averaging in time, equations (4.1) leads to

(4.2) $$\begin{eqnarray}\displaystyle & & \displaystyle U\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}x}+V\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}y}+W\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}z}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{\unicode[STIX]{x2202}\langle (u-U)(u-U)\rangle }{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}\langle (u-U)(v-V)\rangle }{\unicode[STIX]{x2202}y}+\frac{\unicode[STIX]{x2202}\langle (u-U)(w-W)\rangle }{\unicode[STIX]{x2202}z}\nonumber\\ \displaystyle & & \displaystyle \quad =-\frac{1}{\unicode[STIX]{x1D70C}}\frac{\text{d}P}{\text{d}x}+\unicode[STIX]{x1D708}\left(\frac{\unicode[STIX]{x2202}^{2}U}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}^{2}U}{\unicode[STIX]{x2202}y^{2}}+\frac{\unicode[STIX]{x2202}^{2}U}{\unicode[STIX]{x2202}z^{2}}\right)-F,\end{eqnarray}$$

where $F$ is the time average of $f$ . For simplicity, let us consider a periodic computational domain with constant half-channel height. Averaging in both the streamwise and the spanwise directions, equations (4.2) leads to

(4.3) $$\begin{eqnarray}\displaystyle & & \displaystyle -\frac{1}{\unicode[STIX]{x1D70C}}\frac{\text{d}P}{\text{d}x}+\unicode[STIX]{x1D708}\frac{\text{d}^{2}U_{xy}}{\text{d}z^{2}}-F_{xy}\nonumber\\ \displaystyle & & \displaystyle \quad -\,\frac{\text{d}\langle (U-U_{xy})(W-W_{xy})\rangle _{xy}}{\text{d}z}-\frac{\text{d}\langle (u-U)(w-W)\rangle _{xy}}{\text{d}z}=0,\end{eqnarray}$$

where the first term represents a mean pressure force, the second term represents the viscous stress, the third term represents the drag force, the fourth term represents the dispersive stress and the fifth term represents the Reynolds stress. If one neglects dispersive stresses, the Reynolds-averaged streamwise momentum equation reads

(4.4) $$\begin{eqnarray}-\frac{1}{\unicode[STIX]{x1D70C}}\frac{\text{d}P}{\text{d}x}+\unicode[STIX]{x1D708}\frac{\text{d}^{2}U_{xy}}{\text{d}z^{2}}-\frac{\text{d}\langle (u-U)(w-W)\rangle _{xy}}{\text{d}z}-C_{ds}U_{xy}^{2}\frac{\text{d}A_{f}(z)}{\text{d}z}=0,\end{eqnarray}$$

where $C_{ds}$ is the sectional drag coefficient and is often assumed to be constant, $A_{f}(z)$ is the frontal area under height $z$ per unit planar area, and the subscript $xy$ indicates streamwise and spanwise averaging. At high Reynolds numbers and in the absence of a streamwise pressure gradient, equation (4.4) reduces to

(4.5) $$\begin{eqnarray}-\frac{\text{d}\langle (u-U)(w-W)\rangle _{xy}}{\text{d}z}-C_{ds}U_{xy}^{2}\frac{\text{d}A_{f}(z)}{\text{d}z}=0,\end{eqnarray}$$

which is often the starting point of rough-wall drag-force modelling.

Dispersive stresses are neglected in (4.5). Although dispersive stresses are not negligible in the roughness layer (Lien & Yee Reference Lien and Yee2004; Martilli & Santiago Reference Martilli and Santiago2007; Barlow & Coceal Reference Barlow and Coceal2008), their effects can be modelled and absorbed in the drag coefficient $C_{ds}$ ; and therefore neglecting dispersive stresses is usually not problematic for drag-force modelling.

Figure 20.

The Reynolds stresses, the dispersive stresses and the sum of the two for cases (a) A008, (b) A100, (c) A440, (d) S008, (e) S100 and (f) S440. The dot-dashed lines are the Reynolds stresses, the dashed lines are the dispersive stresses and the thick solid lines are the total stresses. The thin solid line corresponds to $\unicode[STIX]{x1D70F}^{+}=-1+z/L_{z}$ . The red line shows the dispersive stress in a rough-wall channel, where the roughness is staggered cubes of $25\,\%$ surface coverage density (Leonardi & Castro Reference Leonardi and Castro2010).

This is not the case for sparsely mounted cubes. Figure 20 shows the Reynolds stress $\unicode[STIX]{x1D70F}_{xz,R}=\langle (u-U)(w-W)\rangle _{xy}$ , the dispersive stress $\unicode[STIX]{x1D70F}_{xz,R}=\langle UW\rangle _{xy}-U_{xy}W_{xy}$ and the ‘total’ stress $\unicode[STIX]{x1D70F}_{xz,t}=\unicode[STIX]{x1D70F}_{xz,R}+\unicode[STIX]{x1D70F}_{xz,D}$ for cases A/S008, A/S100 and A/S440. Here, the subscript $xz$ in $\unicode[STIX]{x1D70F}_{xz}$ denotes the tensor direction, and the subscript $xy$ in $U_{xy}$ and $W_{xy}$ denotes streamwise and spanwise averaging. The ‘total’ stress does not contain the subgrid stress nor the viscous stress. Subgrid stresses should be small if the flow is well resolved. Viscous stresses should be small if the flow is at high Reynolds numbers. Figure 20 shows the stresses above $z/h=0.25$ . The numerical solutions are less susceptible to numerical errors away from the wall (Kawai & Larsson Reference Kawai and Larsson2012; Larsson et al. Reference Larsson, Kawai, Bodart and Bermejo-Moreno2016). We make a few observations. First, the dispersive stress is practically zero at a coverage density $\unicode[STIX]{x1D706}_{p}=25\,\%$ , and therefore it is not responsible for momentum transfer. Second, the total stresses in our LES follow $\unicode[STIX]{x1D70F}^{+}=-1+z/L_{z}$ closely above $z/h=1$ in all cases, and therefore the energy-containing scales in flows are ‘well resolved’ above the cube crest and the simulations are well converged. Third, the dispersive stresses in our LES are comparable with the Reynolds stresses. For cases A/S008, the dispersive stresses are responsible for more than half of the momentum fluxes at a distance approximately $h$ above the cubes. For cases A/S100 and A/S440, the Reynolds stresses are the dominant component but the dispersive stresses are not negligible. If one includes the dispersive stress, the streamwise momentum equation is

(4.6) $$\begin{eqnarray}-\frac{\text{d}\langle (u-U)(w-W)\rangle _{xy}}{\text{d}z}-\frac{\text{d}(\langle UW\rangle _{xy}-U_{xy}W_{xy})}{\text{d}z}-C_{ds}U_{xy}^{2}\frac{\text{d}A_{f}(z)}{\text{d}z}=0,\end{eqnarray}$$

where we have neglected the viscous and the pressure forces.

Figure 21.

Same as figure 20, but for the streamwise component.

In (4.6), both the dispersive stress and the Reynolds stress need be modelled. To model the Reynolds stress, one could resort to an eddy viscosity model, and $-\unicode[STIX]{x1D70F}_{R,xz}=\unicode[STIX]{x1D708}_{T}\,\text{d}U/\text{d}z$ . As $-\unicode[STIX]{x1D70F}_{R,xz}$ is positive definite, the eddy viscosity is also positive definite and therefore well defined. An eddy viscosity model is not possible for the dispersive stress, because it is not a positive definite quantity (see figure 20 a), and one has to take a different approach. A similarity solution was recently proposed in Anderson et al. (Reference Anderson, Yang, Shrestha and Awasthi2018). The similarity solution models the flow above the roughness and may be useful for modelling the dispersive stress, although, in its present form, the similarity solution assumes streamwise homogeneity and will have to be adapted for the flows considered in this work. In addition to the $xz$ component of the dispersive stress, the $xx$ component quantifies spanwise heterogeneity as well. Figure 21 shows the root-mean-square (r.m.s.) of the streamwise component of the dispersive stress tensor $\unicode[STIX]{x1D70F}_{xx,D}$ as a function of the wall-normal distance; $\unicode[STIX]{x1D70F}_{xx,D}$ is not negligible above sparsely packed cubes but is small if the cubes are packed closely ( $\unicode[STIX]{x1D706}_{p}=25\,\%$ ).

The dispersive stress (the $xz$ component) shows up in the $x$ -momentum equation and therefore affects the mean flow scaling. Figure 22 shows the mean velocity profiles in cases A/S008, A/S100 and A/S440. For cases A/S008, the cubes are very sparsely packed and the mean flow is not very different from that of a flat plate. The mean flows in cases A/S440 show typical behaviours as in rough-wall boundary layers (see e.g. Raupach et al. Reference Raupach, Antonia and Rajagopalan1991; Squire et al. Reference Squire, Morrill-Winter, Hutchins, Schultz, Klewicki and Marusic2016). While the secondary flows do not affect much the mean flows in A/S440 and A/S008, the mean flows in A/S100 are noticeably modified by the secondary flows. Owing to the secondary-flow-induced mixing, the mean velocity profile is nearly flat above $z/h=1$ in A100.

Figure 22.

Mean velocity profiles. Squares are used for aligned cubes. Crosses are used for staggered cubes. Cases A/S008 are shown in blue. Cases A/S100 are shown in red. Cases A/S440 are shown in yellow. The thin solid line corresponds to the expected mean flow when there are no wall-mounted cubes, i.e.  $U^{+}=1/\unicode[STIX]{x1D705}\log (z/z_{o})$ .

5 Concluding remarks

We conduct LES and DNS of flow in a half-channel with aligned and staggered cubes at the bottom wall. We study the properties of cube-roughened surfaces with a coverage density between $\unicode[STIX]{x1D706}_{p}=0.08\,\%$ and $\unicode[STIX]{x1D706}_{p}=4.4\,\%$ . Although flow over cube arrays is a well-studied model problem for rough-wall turbulent boundary layers, there have not been many studies on rough walls with such a low coverage density ( $\unicode[STIX]{x1D706}_{p}<1\,\%$ ).

Our findings are summarized below.

1. (i) For both aligned and staggered arrangements, we see an increase of the drag coefficient $C_{d}$ and the friction coefficient $C_{f}$ as a function of $\unicode[STIX]{x1D706}_{p}$ at low surface coverage densities. Since the seminal work by Raupach (Reference Raupach1992), rough-wall drag-force modelling has been predominantly focusing on modelling the effects of flow sheltering, which only leads to a reduced drag coefficient and friction coefficient. Hence, the findings here pose a challenge to rough-wall drag-force modelling.

2. (ii) The observed increases in both $C_{d}$ and $C_{f}$ are the result of secondary turbulent flows. The secondary flows lead to HMPs at cube locations, which increases the incoming wind speed to wall-mounted cubes, and in turn leads to an increase in the drag coefficient. Visually, the secondary flows manifest as streamwise-elongated counter-rotating vortices in the mean flow fields. These secondary vortices are different from the commonly observed turbulent streaks, which are transient in nature (in time or in space).

3. (iii) The amount of increase in the drag coefficient is determined by the strength of the counter-rotating vortices and the distance between two neighbouring pairs of counter-rotating vortices, which may be parametrized with $s_{x}/h$ and $s_{y}/L_{z}$ , respectively. Detailed analysis shows that $(C_{d}/C_{R}-1)(s_{y}/L_{z})^{\unicode[STIX]{x1D6FC}}$ collapses as a function of $s_{x}/h$ for both aligned and staggered cube arrays, thereby providing a possible approach to drag-force modelling at low surface coverage densities.

4. (iv) The counter-rotating vortices enhance mixing and carry a significant part of the wall-normal momentum flux, leading to a non-negligible dispersive stress $\langle UW\rangle _{xy}-U_{xy}W_{xy}$ . This poses yet another challenge to rough-wall drag-force modelling, as the common practice is to neglect the dispersive stress in the momentum balance.

5. (v) We also found two streamwise bands of high-stress regions beneath the far wake of a wall-mounted cube. As wakes have velocity deficits, the above finding is rather counter-intuitive.

The problem of flow over cube-roughened surfaces at low coverage densities can be approached from many perspectives. We have approached this problem from the perspective of drag-force modelling, and we have focused predominantly on first-order flow quantities. Rough surfaces with low coverage densities is still an under-explored area; we hope that this work will provide some preliminary yet useful information about the flow in that regime.

Figure 23.

Contours of the instantaneous streamwise velocity on an $x$ – $y$ plane at a distance $0.5h$ from the wall, a $y$ – $z$ plane through the centre of a column of cubes, and an $x$ – $z$ plane through the centre of a row of cubes.