2.1. Simulation set-up

We use the in-house finite-difference code Vicar3D to solve the incompressible filtered Navier–Stokes equations, including the Vreman subgrid-scale model (Vreman Reference Vreman2004; You & Moin Reference You and Moin2007) and the iWMLES model for the near-wall closure (Yang et al. Reference Yang, Sadique, Mittal and Meneveau2015). Details of the code and the LES set-up can be found in Mittal et al. (Reference Mittal, Dong, Bozkurttas, Najjar, Vargas and von Loebbecke2008) and Yang et al. (Reference Yang, Sadique, Mittal and Meneveau2015). The code has been extensively validated (see, e.g., Mittal et al. Reference Mittal, Dong, Bozkurttas, Najjar, Vargas and von Loebbecke2008; Bhardwaj & Mittal Reference Bhardwaj and Mittal2012; Zheng, Hedrick & Mittal Reference Zheng, Hedrick and Mittal2013; Vedula et al. Reference Vedula, Fortini, Seo, Querzoli and Mittal2014). A validation of particular interest can be found in Yang et al. (Reference Yang, Sadique, Mittal and Meneveau2015), where simulation results for channel flow with the bottom wall mounted with cubes have been compared with experimental measurements, showing good agreement in the predicted mean velocity distribution between the cubes. Considering the fact that in LES the subgrid-scale turbulence is not explicitly computed, we restrict the current analysis to the spatially/temporally averaged velocity profiles, which are expected to be less sensitive to subgrid-scale modelling than, say, second-order statistics, and have been shown to be reproduced accurately by the code.

The roughness elements considered have the shape of rectangular prisms. The simulated flow is a spatially growing boundary layer, for which the boundary layer height grows downstream. This enables us to sample results at various downstream distances, and thus test that the roughness parameters to be obtained (i.e. $z_{o}$ , $d$ ) are independent of the outer boundary layer conditions and the Reynolds number which vary slowly as a function of  $x$ . The inflow is generated via an extended rough-wall rescaling–recycling technique (Yang & Meneveau Reference Yang and Meneveau2016). The boundary layer thickness at the inflow, ${\it\delta}_{0}$ , is $4h_{m}$ (for single-height roughness cases) to $6h_{m}$ (for varying-height roughness cases), where $h_{m}$ is the averaged roughness height. Length scales are expressed in terms of the incoming boundary layer thickness (i.e. ${\it\delta}_{0}=1$ ) used for the reference case (simulations with single-height roughness). In these units, the value of the mean height of roughness elements is always $h_{m}=0.25$ , while the computational domain size is $16\times 6\times 6$ in the streamwise, spanwise and wall normal directions respectively. The Cartesian mesh size for LES is $256\times 96\times 96$ . The average roughness height $h_{m}$ is resolved by eight grid points in the vertical direction. For the cases with varying height, the inflow boundary layer is increased to ${\it\delta}_{0}^{\prime }=1.5$ (to keep the largest height from approaching the height of the incoming boundary layer thickness). The Reynolds number based on the inlet boundary layer thickness and free-stream velocity is $Re_{{\it\delta}_{0}}={\it\delta}_{0}U_{0}/{\it\nu}=10^{6}$ , i.e. sufficiently high for the flow to be in the fully rough regime ( $U_{0}$ is the free-stream velocity and ${\it\nu}$ is the kinematic viscosity). The surfaces of the roughness elements themselves, as well as the bottom wall, are assumed to be smooth, with no subgrid-scale roughness. A total simulation time of 100 flow-through times is computed, using a time step of $\text{d}t=0.33\,\text{d}x/U_{0}$ , where $\text{d}x$ is the grid spacing in the streamwise direction. The flow-through time is defined as $T_{0}=L_{x}/U_{0}$ , where $L_{x}=16$ is the computational domain size. The code is massively parallelizable and MPI is used for interprocessor communication.

2.2. The roughness geometries considered

Large-eddy simulations of developing turbulent boundary layer flow over three sets of rough walls are conducted, for various roughness configurations. For set I, the solidity is systematically varied for both aligned and staggered arrangements of cubic roughness elements. Figure 1 shows the repeating tiles for this set of rough walls. Set II consists of aligned rectangular roughness with 11 % surface coverage and non-uniform roughness heights. The roughness height is selected from a Gaussian distribution with several values of standard deviation: ${\it\sigma}_{h}=0.24h_{m}$ , $0.35h_{m}$ , $0.5h_{m}$ . The probability-density function (PDF) of the roughness heights used in the LES is compared with the Gaussian distribution in figure 2. For each ${\it\sigma}_{h}$ , four realizations of rough walls are randomly generated, thus resulting in 12 rough walls in total in this set. Each case in this set is denoted with LXXStdXXA/S-X, where the two digits following L denote the value equal to $100{\it\lambda}_{f}$ , the digits following Std represent $100{\it\sigma}_{h}/h_{m}$ , ‘A’ and ‘S’ stand for aligned and staggered arrangement respectively and the last digit is the simulation index associated with the same ${\it\sigma}_{h}$ . Figure 3 shows a visualization of the surface for case L11Std50A-1 (figure 3 a), and instantaneous and averaged streamwise velocities on a plane at height $z=h_{m}$ (figure 3 b,c), for flow case L11Std50-1 in rough-wall set II.

Figure 1. Repeating tiles for the first set of simulations denoted by ‘I’. The roughness elements are cubic and are coloured black. Each case in this set is denoted with LXXS/A, where L stands for ${\it\lambda}_{f}$ , the number following L is $100{\it\lambda}_{f}$ , S stands for ‘staggered’ and ‘A’ stands for ‘aligned’: (a) L03S, (b) L03A, (c) L06S, (d) L06A, (e) L11S, (f) L11A, (g) L25S, (h) L25A. The surface coverage density is systematically varied within this set of rough walls for both perfectly aligned and staggered roughness arrangements.

Figure 2. The discretized roughness height PDF compared with the Gaussian PDF. The standard deviation in roughness height is (from a–c) $0.24h_{m}$ , $0.35h_{m}$ , $0.5h_{m}$ ; $h_{m}\pm \text{std}(h)$ are marked with dashed lines.

Figure 3. Flow visualization from LES of flow over an array of wall-attached rectangular prisms, for case L11Std50A-1. (a) Bottom surface geometry and transparent isovelocity surface, (b) instantaneous and (c) averaged streamwise velocities on a plane at height $z=h_{m}$ , for cubic roughness L11Std50A-1. Here, $h_{m}$ is the mean height of the roughness element. The velocity is normalized with the free-stream velocity $U_{0}$ and the length is normalized with ${\it\delta}_{0}=4h_{m}=4h$ .

Using set III, the hydrodynamic response to the lateral displacement of roughness elements is examined. In set III, we systematically vary the staggering at two roughness solidities, ${\it\lambda}_{f}=0.06,0.11$ . The repeating tiles for this set of rough walls are shown in figure 4.

Figure 4. Repeating tiles for set III. Roughness elements are cubic and are coloured black. Each case in this set is denoted as LXXSXX, where L stands for ${\it\lambda}_{f}$ and S stands for percentage of staggering: (a) L06SXX, (b) L11SXX. The percentage of staggering is defined as $200l_{s}\sqrt{{\it\lambda}_{f}}/h$ , where $h$ is the cube height and $l_{s}$ is the roughness displacement in the spanwise direction. A perfectly aligned arrangement is 0 % staggering and a perfectly staggered arrangement is 100 % staggering. The number following L is $100{\it\lambda}_{f}$ and the number following S is the percentage of staggering. For ${\it\lambda}_{f}=0.06$ , six different percentages of staggering, spaced equally between 0 % to 100 %, are considered, and for ${\it\lambda}_{f}=0.11$ , seven different cases of staggering are considered. Hence, this set includes 13 separate cases.

Figure 5. Comparison of simulated velocity profiles within the roughness layer and the exponential fit (thin red lines), for set I. The hollow symbols are for the staggered arrangement and the filled symbols are for the aligned arrangement.

Figure 6. Comparison of simulated velocity profiles within the roughness layer and the exponential fit (thin red lines), for set II. Cases corresponding to different ${\it\sigma}_{h}$ are displaced horizontally, while four random generated rough walls with the same ${\it\sigma}_{h}$ (from bottom to top: cases L11StdXXA-1/2/3/4) are displaced vertically for clarity of presentation: (a)  ${\it\sigma}_{h}=0.24h_{m}$ , (b)  ${\it\sigma}_{h}=0.35h_{m}$ , (c)  ${\it\sigma}_{h}=0.50h_{m}$ .

2.3. Mean velocity profile within the roughness layer

The velocity profile is examined in detail within the roughness layer, i.e. below $z=H$ , where $H=h_{m}+{\it\sigma}_{h}$ and ${\it\sigma}_{h}$ is the standard deviation in roughness height. Here, $H$ is referred to as the generalized roughness layer height, which for the case of uniform roughness height is simply $H=h=h_{m}$ . For the general cases, the choice of $H=h_{m}+{\it\sigma}_{h}$ is guided by examination of the data and by practicality. First, the results for the case of a bimodal distribution of heights will clearly show that the break between the logarithmic and exponential behaviours occurs at the height of the larger roughness elements. This has also been observed by Jiang et al. (Reference Jiang, Jiang, Liu and Sun2008) and Hagishima et al. (Reference Hagishima, Tanimoto, Nagayama and Meno2009). Since for a bimodal distribution with equal probabilities the standard deviation is equal to the difference between the mean and the maximum (or minimum) height, the height of the larger element equals $h_{m}+{\it\sigma}_{h}$ . Mean velocity profiles for the case with Gaussian height distributions further confirm that this is a good approximation of where the mean velocity profile transition takes place.

The velocity profile within the roughness layer is obtained via spatially averaging the temporally averaged velocity field for two repeating tiles (for perfectly aligned arrangements) or one repeating tile (for all other arrangements) in the streamwise direction at the middle of the computational domain. Figures 5–7 show the LES computed profiles within the roughness layer for all three sets of rough walls. In each of the individual plots, an exponential has been fitted in a range of $z$ between $H/2$ and $H$ . Fits are constrained to pass through $(U/U_{H}=1,z/H=1)$ , thus, the parameter fitted is the attenuation coefficient $a$ for each case.

Figure 7. Comparison of the simulated velocity profiles within the roughness layer and the exponential fit (thin red line), for set III. The cases with two different solidities are displaced horizontally. Within each panel, the percentage of staggering increases from bottom to top (from bottom to top: cases L06S00/25/50/75/100 ( ${\it\lambda}_{f}=0.06$ ) (a), L11S00/17/33/50/67/83/100 ( ${\it\lambda}_{f}=0.11$ ) (b)).

In figure 8, we summarize all cases simulated by plotting $\ln (U/U_{H})/a$ against the wall normal distance, using the fitted value of $a$ . The linear behaviour indicates that for all rough walls considered in our LES, except for the near-wall region, an exponential profile is a very good representation of the velocity profile within the roughness layer (the collapse among the cases is due to the fitting and is thus not surprising in this plot).

Figure 8. Velocity profiles for all LES cases plotted together within the roughness layer, in linear–log scale. The collapse onto a line confirms exponential behaviour, although each case is characterized by a different (fitted) parameter  $a$ .

These observations provide good support for the generality of the exponential mean velocity profile within a roughness layer, at least for the specific class of roughness elements that have a well-defined length scale, each with reasonably uniform cross-sectional area, and well-defined flow separation points.

Figure 9. Sketch of layers of the assumed mean velocity profile representing the horizontally averaged velocity of the fluid in the flow domain (excluding the roughness elements). From top to bottom a standard logarithmic layer characterized by roughness and displacement lengths, and friction velocity (and further above possibly a wake), and an exponential layer characterized by the roughness element height and an attenuation coefficient.

3.1. Assumed shape function for mean velocity profile

The proposed approach is inspired by the von Karman–Pohlhausen method (Pohlhausen Reference Pohlhausen1921) and begins by postulating a shape function for the vertical distribution of the horizontally averaged velocity profile $U(z)$ . We divide the height into two layers. (1) The roughness layer for heights below the roughness elements (for now we use the symbol $h$ , but for cases with varying heights this can be replaced by $H$ which includes the mean and standard deviation of the element heights). Based on the results presented in § 2, for the first layer we use an exponential profile. (2) Above $z=h$ (or $H$ ), we assume a standard logarithmic law (Raupach, Antonia & Rajagopalan Reference Raupach, Antonia and Rajagopalan1991; Jiménez Reference Jiménez2004), arising from a constant stress layer in which turbulence alone conveys momentum in the vertical direction. In sum, the following shape function is assumed in replacement of the integration of the RANS equations:

(3.1) $$\begin{eqnarray}U(z)=\left\{\begin{array}{@{}ll@{}}U_{h}\exp [a(z-h)/h], & 0<z<h,\\ u_{{\it\tau}}/{\it\kappa}[\log [(z-d)/z_{o}]+{\it\Pi}W(z/{\it\delta})], & h\leqslant z<{\it\delta},\end{array}\right.\end{eqnarray}$$

where ${\it\kappa}$ is the von Karman constant, ${\it\delta}$ is the boundary layer thickness and ${\it\Pi}W(z/{\it\delta})$ is a wake function (Coles Reference Coles1956), with the parameter ${\it\Pi}$ typically of order unity (we take ${\it\Pi}=0.2$ in this study) and $W(1)=2$ (Jiménez Reference Jiménez2004). The parameter $U_{h}$ is the mean velocity at the roughness element height. Above $z={\it\delta}$ , we take $U=U_{0}$ , the free-stream velocity. Figure 9 illustrates the various layers. It should be noted that the profile does not vanish at $z=0$ , but for now we are not interested in an additional, possibly very thin, further layer in which the velocity decreases rapidly from $\exp (-a)U_{h}$ down to zero. While the numerical results displayed in figures 5–7 show that the exponential fit deviates from the data in the bottom 20–25 % of the roughness layer, experimental and DNS data typically show a more sudden drop closer to the surface (a thinner near-wall layer, see, e.g., figure 4 in Macdonald (Reference Macdonald2000) and figure 4 in Leonardi & Castro (Reference Leonardi and Castro2010)). In our simulations, the deviations from the exponential layer only occur for the 1–2 grid points nearest to the bottom surface, where the LES results can be affected by the wall model and numerical errors. Hence, we take the view that the validity of the exponential profile extends even closer to the surface than our LES results show.

It should be noted that at $z=h$ , (3.1) implies a sign change in the curvature of the assumed velocity profile. An inflection point near the roughness layer top is typically observed in rough-wall boundary layers and canopy flows (see, e.g., Raupach, Finnigan & Brunei Reference Raupach, Finnigan and Brunei1996), and is the cause of the instability that in turn leads to mixing among the roughness layer and the inertial layer. Often, hyperbolic tangent-like profiles are used to model this mixing-layer-type profile, and one could imagine introducing a third intermediate layer in the present model. Here, for simplicity, connecting an exponential and a logarithmic profile on both sides serves a similar purpose and provides a good approximation to the mixing-layer-type profile for the present purposes.

It should be noted that there are five unknown parameters in (3.1), i.e. $U_{h}$ , $u_{{\it\tau}}$ , $d$ , $z_{o}$ and $a$ . It will be assumed that we know the velocity outside the boundary layer $U_{0}$ , the boundary layer height ${\it\delta}$ , the element height $h$ (or $H$ ) and the geometrical distribution of the element on the surface. Using this information as input, we require five constraints to express the five unknown parameters $U_{h}$ , $u_{{\it\tau}}$ , $d$ , $z_{o}$ and $a$ as functions of $U_{0}$ , ${\it\delta}$ and $h$ , and knowledge about the spatial distribution of the roughness elements, including the parameter ${\it\lambda}_{f}$ .

Three fundamental constraints can be found, one based on the basic principle of momentum balance and two from continuity of the velocity profile. A fourth constraint is based on relating the displacement height $d$ to the centre of force, the height at which one may consider the effective wall stress to be applied by the roughness elements onto the flow (Jackson Reference Jackson1981). The fifth constraint determines the exponential attenuation coefficient $a$ and is discussed separately in § 3.2, based on considerations of flow sheltering.

First, we take the control volume that contains vertically the whole roughness layer and part of the inertial layer. The downward momentum flux within the inertial layer is balanced by the form drag due to the wall roughness. This leads to the vertically integrated momentum equation:

(3.2) $$\begin{eqnarray}A_{T}u_{{\it\tau}}^{2}=\int _{A_{f}}C_{d}U(z)^{2}\,\text{d}A_{f},\end{eqnarray}$$

where the integration is performed on the projected frontal area $A_{f}$ within the lot area $A_{T}$ . For the right-hand side we have employed the quadratic law for the form drag, i.e. $\text{d}F=C_{d}U(z)^{2}\,\text{d}A_{f}$ . Here, $C_{d}$ is the drag coefficient, for which a typical value $C_{d}=1$ is used (Coceal & Belcher Reference Coceal and Belcher2004). The drag coefficient can depend on the vertical distance, the detailed distribution of the ground roughness, the details of the local flow, etc. (Santiago et al. Reference Santiago, Coceal, Martilli and Belcher2008). The simplification of employing a constant drag coefficient is made here, following common practice in rough-wall modelling (Macdonald Reference Macdonald2000; Coceal & Belcher Reference Coceal and Belcher2004; Harman & Finnigan Reference Harman and Finnigan2007; Di Sabatino et al. Reference Di Sabatino, Solazzo, Paradisi and Britter2008). The viscous skin friction is not considered on the right-hand side of (3.2), based on the assumption of a ‘fully rough’ regime (transitional roughness will be considered in appendix B based on partitioning the drag to include viscous skin friction). For rectangular-prism roughness elements, integrating the exponential velocity profile between $z=0$ and $z=h$ , (3.2) leads to

(3.3) $$\begin{eqnarray}u_{{\it\tau}}^{2}=\frac{1-\exp (-2a)}{2a}{\it\lambda}_{f}C_{d}U_{h}^{2}.\end{eqnarray}$$

Commenting on the accuracy of the exponential fit near the bottom surface, if the integration were performed between $z=0.2h$ and $h$ (the typically observed range of validity of the exponential fit in our LES) instead of between 0 and $h$ , for a representative value of $a=1$ , the resulting drag force from the integration would lead to a difference in $u_{{\it\tau}}$ of only ${\sim}4\,\%$ because the integral is dominated by the profile near $z\sim h$ .

The second condition related to velocity profile continuity imposed at $z=h$ leads to

(3.4) $$\begin{eqnarray}\frac{U_{h}}{u_{{\it\tau}}}=\frac{1}{{\it\kappa}}\log \left(\frac{h-d}{z_{o}}\right),\end{eqnarray}$$

where we have assumed $W(h/{\it\delta})\approx 0$ (an assumption that requires $h/{\it\delta}\ll 1$ but that will be checked based on data in § 4 to hold even in cases where $h/{\it\delta}\sim 0.2$ ). Continuity of total stress at $z=h$ is usually used in rough-wall models to constrain the mixing length (see, e.g., Coceal & Belcher Reference Coceal and Belcher2004). Because we do not invoke the mixing length in the formulation, continuity in total stress is not explicitly used (although it could be used to derive the implied mixing length from the roughness layer at $z=h$ ).

Third, we require the velocity at $z={\it\delta}$ to be the free-stream velocity $U_{0}$ :

(3.5) $$\begin{eqnarray}\frac{U_{0}}{u_{{\it\tau}}}=\frac{1}{{\it\kappa}}\left[\log \left(\frac{{\it\delta}-d}{z_{0}}\right)+2{\it\Pi}\right].\end{eqnarray}$$

Fourth, we use the viewpoint proposed in Jackson (Reference Jackson1981) that the displacement height $d$ can be set equal to the centroid height of the distributed drag force, namely

(3.6) $$\begin{eqnarray}d=\frac{\displaystyle \int _{A_{f}}C_{d}U(z)^{2}z\,\text{d}A_{f}}{\displaystyle \int _{A_{f}}C_{d}U(z)^{2}\,\text{d}A_{f}}.\end{eqnarray}$$

For single-height rectangular-prism roughness after integration of the exponential profile between $z=0$ and $z=h$ , (3.6) simply leads to

(3.7) $$\begin{eqnarray}\frac{d}{h}=\frac{1}{1-\exp (-2a)}-\frac{1}{2a}.\end{eqnarray}$$

Combining (3.3), (3.4) and (3.7), $z_{o}$ for single-height rectangular-prism roughness can be expressed as

(3.8) $$\begin{eqnarray}\frac{z_{o}}{h}=\left(1-\frac{d}{h}\right)\exp \left[-{\it\kappa}\bigg/\sqrt{\frac{1}{2a}C_{d}{\it\lambda}_{f}(1-\text{e}^{-2a})}\right].\end{eqnarray}$$

Replacing (3.8) into (3.4) and (3.5), we obtain the expression for $u_{{\it\tau}}$ :

(3.9) $$\begin{eqnarray}\frac{u_{{\it\tau}}}{U_{0}}=\left[\frac{1}{{\it\kappa}}\log \frac{{\it\delta}-d}{h-d}+\sqrt{\frac{2a}{1-\exp (-2a)}}\frac{1}{\sqrt{C_{d}{\it\lambda}_{f}}}+\frac{2{\it\Pi}}{{\it\kappa}}\right]^{-1},\end{eqnarray}$$

and

(3.10) $$\begin{eqnarray}\frac{U_{h}}{U_{0}}=\left[1+\frac{1}{{\it\kappa}}\left(\log \frac{{\it\delta}-d}{h-d}+2{\it\Pi}\right)\sqrt{C_{d}{\it\lambda}_{f}\frac{1-\exp (-2a)}{2a}}~\right]^{-1},\end{eqnarray}$$

where $d$ is obtained from (3.7). Again, it should be noted that (3.7)–(3.10) are valid only for single-height rectangular-prism roughness elements. Equations (3.2) and (3.4)–(3.6) can be used for general rectangular roughness. The fifth condition to determine $a$ incorporates the effects from mutual sheltering among the roughness elements and is discussed in the next section.

Figure 10. (a) Sketch of the sheltering effect and simplified model in which the complicated flow between roughness elements is assumed to consist of two regions with characteristic velocities $U_{h}$ (unsheltered region) and small velocity in the sheltered region. The case shown is for when the downstream element overlaps with the sheltered region, i.e. when $h_{s}>0$ . (b) A sketch of the volume within the wake of a rectangular-prism roughness that has reduced momentum. The roughness is $h$ in height, $w$ in width and $b$ in streamwise length. The streamwise length of the sheltered region is $L_{s}$ and based on the linear expansion model is given by $L_{s}=h/\tan {\it\theta}$ .

3.2. Volumetric sheltering

In this section, we model the reduction in the momentum in the wakes of rectangular-prism roughness elements and its effects upon the drag on neighbouring roughness elements. This reduction is denoted as volumetric sheltering and has been considered in various prior studies (see, e.g., Raupach Reference Raupach1992; Shao & Yang Reference Shao and Yang2005, Reference Shao and Yang2008). Qualitatively, the extremes have been denoted as ‘ $d$ -type’ and ‘ $k$ -type’ roughness, where in the first case much sheltering occurs and the flow can ‘skim over’ the elements, whereas in ‘ $k$ -type’ roughness each element produces considerable drag (Leonardi, Orlandi & Antonia Reference Leonardi, Orlandi and Antonia2007).

Consider the situation as sketched in figure 10(a). In the wake of rectangular-prism roughness elements there is a sheltered region in which the velocity is lower than in the unsheltered region. Depending on the spacing between the roughness elements, one may be in a sheltered, unsheltered or ‘just sheltered’ condition (the latter is defined as when the sheltering region ‘just’ begins to intersect the downstream element, i.e. when $h_{s}$ changes from $h_{s}=0$ to $h_{s}>0$ ). We consider the case of ‘just sheltered’ and determine the attenuation parameter in the mean velocity profile based on momentum balance, namely that the integrated distributed drag equals the drag on the fully exposed frontal area of an individual roughness element (here we take $A_{f}=wh$ , where $w$ is the element width). We assume the latter to be $C_{DH}U_{h}^{2}wh/2$ , where $C_{DH}$ is the drag coefficient (assumed to be known) that expresses the total drag on the element when using the tip velocity $U_{h}$ as the velocity scale. In other words, for a rectangular prism we can write

(3.11) $$\begin{eqnarray}\int _{A_{f}}C_{d}U(z)^{2}\,\text{d}A_{f}=C_{d}U_{h}^{2}wh\frac{1}{2a}[1-\exp (-2a)]=\frac{1}{2}C_{DH}U_{h}^{2}wh,\end{eqnarray}$$

where for $C_{DH}$ we use the drag coefficient for a surface-mounted cube for which data exist, and (again) $C_{d}$ is assumed to be constant. It is known that approximately $C_{DH}\approx 1.4$ (Akins, Peterka & Cermak Reference Akins, Peterka and Cermak1977; Hussain & Lee Reference Hussain and Lee1980; Curley & Uddin Reference Curley and Uddin2015). Essentially we are assuming that the drag coefficient $C_{DH}$ appropriate for an unsheltered or ‘just sheltered’ cuboid in the array is equal to that of an isolated element. Since we are using $C_{d}=1$ , from this equality we can find the attenuation of the velocity profile. The condition for $a$ for the ‘just sheltered’ condition is given by

(3.12) $$\begin{eqnarray}\frac{1}{2a}[1-\exp (-2a)]=\frac{C_{DH}}{2C_{d}}\approx 0.7.\end{eqnarray}$$

The solution is $a\approx 0.4$ . As roughness elements are placed further distances apart (even lower ${\it\lambda}_{f}$ ), momentum balance (for fully rough conditions) implies that $a$ is unchanged and remains at 0.4. Thus, we refer to this value as $a_{min}=0.4$ . Direct measurements of the attenuation coefficient $a$ characterizing the nearly exponential velocity profile are presented in § 4, and the measurements confirm the numerical value $a_{min}\approx 0.4$ characterizing the limiting attenuation in the limit of small ${\it\lambda}_{f}$ . As elements are placed closer together leading to a sheltering effect, we expect $a$ to increase above 0.4, since the velocity profile will be increasingly attenuated to values smaller than $U_{h}$ as $z$ decreases from $z=h$ down towards the wall.

In order to determine the value of $a$ that depends on the degree of sheltering, we employ again the momentum balance as before. However, because $C_{DH}$ is measured for unsheltered cubes, the momentum balance is only applied for the top portion above the sheltered region for $z>h_{s}$ as if the element only protrudes a height $h_{s}$ above the surface, thus neglecting the contribution of the sheltered region to drag. The sheltering height $h_{s}$ depends on the expansion rate of the wake, as shown in figure 10(b), and will be modelled later. For the momentum balance with sheltering we obtain the condition

(3.13) $$\begin{eqnarray}\int _{h_{s}}^{h}C_{d}U(z)^{2}w\,\text{d}z=\frac{1}{2a}(1-\exp [-2a(1-h_{s}/h)])whC_{d}U_{h}^{2}=\frac{1}{2}C_{DH}U_{h}^{2}w(h-h_{s}),\end{eqnarray}$$

leading to

(3.14) $$\begin{eqnarray}\frac{1}{2a(1-h_{s}/h)}(1-\exp [-2a(1-h_{s}/h)])=\frac{1}{2a_{min}}[1-\exp (-2a_{min})].\end{eqnarray}$$

The solution is simply

(3.15) $$\begin{eqnarray}a=\frac{a_{min}}{1-h_{s}/h}.\end{eqnarray}$$

Hence, knowing $h_{s}/h$ allows us to determine $a$ . For unsheltered cases ( $h_{s}=0$ ), we use $a=a_{min}=0.4$ .

Next, we discuss the determination of $h_{s}$ . Figure 10(b) sketches the volume within the wake of a rectangular-prism element that has reduced-momentum fluid. Wake expansion in the vertical direction can bring in fluid with high momentum; a reduction in the volume of low-momentum fluid would thus be expected. Wake expansion in the spanwise direction, on the other hand, expands the volume of low-momentum fluid. This configuration is physical for near-wall roughness with aspect ratio $h/w$ not large. Naturally, we expect the wake be eaten up from the side if the roughness geometry is ‘stick-like’, which is not captured by the configuration in figure 10(b). The horizontal convective velocity in the roughness sublayer is of the order of $U_{h}$ , while the turbulent transport velocity scale in the vertical and spanwise directions is of the order of the friction velocity $u_{{\it\tau}}$ . As a result, we estimate the wake expansion rate as

(3.16) $$\begin{eqnarray}\tan {\it\theta}=C_{{\it\theta}}\frac{u_{{\it\tau}}}{U_{h}},\end{eqnarray}$$

where $C_{{\it\theta}}$ is a coefficient of order unity that may depend upon element geometry (and in the case of rectangular prisms, the aspect ratio). Any portion of a downstream roughness element that falls in this region of reduced momentum is considered to be sheltered. Determination of $h_{s}$ thus proceeds by calculating the expansion rate using (3.16), and using it to evaluate the area of the sheltered region $A_{s}$ as the frontal area of the closest downstream rectangular-prism roughness elements. For cases in which sideways growth of the wake causes partial sheltering of the width, the momentum argument presented above is equivalent to setting $h_{s}=A_{s}/w$ in (3.13) and thus in (3.15) (a more precise definition is given in (3.21)).

Depending on roughness element placement on the surface, the calculation of the sheltered area can be more or less involved. The most general procedure is to, first, find all upstream roughness elements that could shelter the roughness under consideration and, second, calculate the sheltered frontal area due to each roughness element found in the first step. To ensure that no roughness interactions are missed, a large enough upstream search domain of size $\ell \times (2\ell +w)$ is used, with $\ell =3h_{max}U_{h}/u_{{\it\tau}}$ , as shown in figure 11(a). Multiple roughness elements could shelter the roughness element under consideration, e.g. in figure 11(b), the roughness element under consideration is sheltered by three upstream roughness elements. In general, if $n$ roughness elements upstream leave sheltering imprints on the roughness under consideration, we denote each of the sheltering areas by $S_{i}(y)$ ( $i=1,2,3\ldots n$ ). For example, in figure 11(b), $S_{1}(y)=h_{1}[\mathscr{H}(y-y_{1a})-\mathscr{H}(y-y_{1b})]$ , where $\mathscr{H}(\cdot )$ is the Heaviside function and $y_{1a}$ and $y_{1b}$ are the beginning and end locations for the sheltering area of the $i=1$ element. The combined sheltered area of the roughness element under consideration is ${\rm\Delta}A_{s}=\int _{0}^{w}\max [S_{1}(y),S_{2}(y),\ldots ]\,\text{d}y$ . It can be noticed that in this way we do not double count the sheltering areas.

Figure 11. (a) Sketch of the size of the upstream search domain for roughness elements that could shelter the roughness under consideration (black cube). To determine $\ell$ , $h_{max}$ is the maximum height of all roughness on the wall, $U_{h}/u_{{\it\tau}}$ is iterated. (b) Sketch of the frontal area of the roughness under consideration. The rectangles $i=1{-}3$ are sheltering areas caused by three upstream roughness elements.

For simple cases, this procedure can be performed analytically. For example, in aligned cube arrays, in which the interaction is limited to the roughness element under consideration and its nearest roughness element upstream (no need to search for a large area), the entire width of the element under consideration is in the sheltered region if sheltering occurs. Then, $h_{s}$ is simply given by

(3.17) $$\begin{eqnarray}h_{s}=\max [h-C_{{\it\theta}}(u_{{\it\tau}}/U_{h})L_{x},0],\end{eqnarray}$$

where $L_{x}$ is the horizontal distance between roughness elements. Substitution of (3.17) into (3.15) leads to

(3.18) $$\begin{eqnarray}a=0.4\max \left[1,C_{{\it\theta}}^{-1}\frac{h}{L_{x}}\frac{U_{h}}{u_{{\it\tau}}}\right],\end{eqnarray}$$

which, together with (3.7)–(3.10), determines the model. More generally, this procedure of calculating the sheltering can be performed numerically (see appendix A).

Once $h_{s}$ is determined from such geometrical arguments, (3.15) provides the final condition to determine the unknowns in (3.1). In the appendix, some other evaluations of $h_{s}$ and $a$ appropriate for some sample geometries considered later in this paper are also provided (e.g. staggered cubes, non-uniform heights, etc.). These evaluations of $h_{s}$ are only relevant for cases in which the spacing between elements is smaller than $L_{s}$ , and thus $a>a_{min}$ . We still have an as of yet undetermined parameter, $C_{{\it\theta}}$ .

3.3. The wake expansion rate

The wake expansion rate is expressed as $\tan {\it\theta}=C_{{\it\theta}}u_{{\it\tau}}/U_{h}$ . We consider here rectangular-prism roughness elements, allowing for different height/side/width ratios. We recognize that the expansion rate can be different for different ratios. Consider cubes and 2D transverse ribs as examples to be contrasted. Most data available on rectangular roughness suggest that the region affected by sheltering is shorter for 3D cubes than for 2D ribs. A similar trend is known to exist for recirculation regions which are typically shorter for 3D objects compared with 2D obstructions (although we clarify that the sheltering region is different from the recirculation region). Thus, we expect stronger volumetric sheltering effects for 2D roughness compared with 3D roughness elements, which is related to the relatively weaker spreading rate of the sheltering region for 2D roughness elements.

To estimate the differing expansion rates $\tan {\it\theta}$ for objects of different aspect ratios, we consider the momentum balance for the case of ‘just sheltered’ but for various aspect ratios $w/h$ , as sketched in figure 12. The total drag on the element, $C_{DH}U_{h}^{2}hw/2$ , is equated to the vertical momentum flux in the rectangular area of length $h/\tan {\it\theta}$ and width $w+2h$ that is physically associated with wake expansion before the sheltering region has been reduced vertically such that the wake ‘touches’ the ground (dashed line in figure 12). The width $2h+w$ comes from the fact that the rate of the wake being ‘eaten up’ from the top equals the rate at which the wake expands sideways, and by the time the entire wake is ‘eaten away’ from the top (i.e. by a distance $h$ ), the side expansion is then also $h$ , on both sides. The vertical flux of momentum per unit area is $u_{{\it\tau}}^{2}$ . Thus, the balance can be written as

(3.19) $$\begin{eqnarray}\frac{h}{\tan {\it\theta}}(2h+w)u_{{\it\tau}}^{2}=\frac{1}{2}C_{DH}hwU_{h}^{2}.\end{eqnarray}$$

We replace $\tan {\it\theta}=C_{{\it\theta}}(u_{{\it\tau}}/U_{h})$ and solve for $C_{{\it\theta}}$ to obtain

(3.20) $$\begin{eqnarray}C_{{\it\theta}}=\left(\frac{6}{C_{DH}}\frac{u_{{\it\tau}}}{U_{h}}\right)\left(\frac{1}{3}+\frac{2h}{3w}\right)\approx O(1)\times \left(\frac{1}{3}+\frac{2h}{3w}\right)=1-\frac{2}{3}\left(1-\frac{h}{w}\right),\end{eqnarray}$$

where the prefactor $O(1)$ arises because $C_{DH}=0.7$ , and while the ratio $(u_{{\it\tau}}/U_{h})$ depends on details of the flow and roughness configurations, in most cases it is of the order of ${\sim}O(10^{-1})$ . Therefore, as a model for $C_{{\it\theta}}$ , we simply take $C_{{\it\theta}}=1/3+2h/3w$ . As a result, for a given ratio $u_{{\it\tau}}/U_{h}$ the expansion rate for 2D bar roughness elements ( $h\ll w$ ) is predicted to be one-third of that for cubic roughness elements ( $h=w$ ). The validity of (3.20) is associated with the validity of the assumed sheltering region shape in figure 10(b). Because such a shape is not physically reasonable for slender ‘stick-like’ roughness elements (e.g. a canopy of long wall-attached vertical cylinders or tall buildings for which the sideways expansion should transition into an inward reduction at some distance downstream), the use of (3.20) should be limited to roughness elements of an aspect ratio $h/w$ that does not exceed an upper threshold (here we typically have $h/w\lesssim 2$ ).

This completes the description of the analytical roughness model for rectangular prisms, since based solely on geometric considerations one may find $C_{{\it\theta}}$ , $\tan {\it\theta}$ , $h_{s}$ and thus $a$ , followed by the roughness parameters $z_{0}$ , $d$ , $U_{h}$ and $u_{{\it\tau}}$ .

Figure 12. Top view of the rough wall. Consider the case of ‘just touching’. The roughness is shown by the thick solid line. The roughness width is $w$ , its length is $b$ and its height is $h$ . The rectangular region enclosed by the dashed line is argued to carry the vertical flux of momentum associated with the wake growth and vertical reduction of the sheltered region.

3.4. Summary of the wall and sheltering model

We briefly summarize the rough-wall model and outline how it can be applied to rough walls mounted with rectangular-prism roughness elements. First, an initial guess for the attenuation coefficient $a$ is made, e.g. $a=a_{min}=0.4$ . Then, the four unknowns $U_{h}$ , friction velocity $u_{{\it\tau}}$ , effective roughness height $z_{o}$ and zero-plane displacement $d$ are obtained by solving the constraint equations (3.2) and (3.4)–(3.6). With $u_{{\it\tau}}$ and $U_{h}$ known, their ratio is used to determine the angle $\tan ({\it\theta})$ by means of (3.16) and (3.20).

The next iteration of the attenuation coefficient $a$ is determined using (3.15), but this expression requires the equivalent sheltered layer height, $h_{s}$ . A condition to find $h_{s}$ that can be applied in general configurations of rectangular-prism roughness elements can be written as

(3.21) $$\begin{eqnarray}\int _{0}^{h_{s}}w_{t}(z)\,\text{d}z=A_{s},\end{eqnarray}$$

where $A_{s}=\sum {\rm\Delta}A_{s}$ is the total sheltered frontal area in a given region of interest and $w_{t}(z)$ is the total flow direction projected width of roughness elements at height $z$ , also over the same area of interest. Equation (3.21) reduces to $h_{s}={\rm\Delta}A_{s}/w$ for single-height single-aspect-ratio roughness. We need (3.21) mainly to account for roughness elements with non-uniform height distribution. In such a case, in fact $h_{s}$ could be larger than the height of some roughness elements.

The sheltered frontal area $A_{s}$ needs to be calculated geometrically using the sheltering model. The volumetric sheltering behind a cube was sketched in figure 10, and the wake expansion rate is given by (3.16) and (3.20), i.e. $\tan ({\it\theta})=[1/3+2h/3w](u_{{\it\tau}}/U_{h})$ , where $h$ is the height of the rectangular roughness and $w$ is its width. When the surface contains elements with different aspect ratios, $\tan ({\it\theta})$ is calculated for each roughness element aspect ratio. That is to say, while the ratio $(u_{{\it\tau}}/U_{h})$ is an averaged quantity, $\tan ({\it\theta})$ is dependent on each individual element.

The steps to implement the sheltering model are as follows.

1. (i) Begin with an initial guess for $u_{{\it\tau}}/U_{h}$ (e.g. 0.1) and use it as an initial guess of $\tan ({\it\theta})$ .

2. (ii) For every roughness element under consideration, identify all upstream roughness elements that could shelter it.

3. (iii) Calculate the area sheltered by each element identified in step (ii) and for all elements under consideration to obtain $A_{s}$ .

4. (iv) Calculate the sheltering height $h_{s}$ with (3.21).

5. (v) Use (3.15) and (3.21) to obtain $a$ . (Note that $a>a_{min}=0.4$ .)

6. (vi) Solve for $U_{h}$ , $u_{{\it\tau}}$ , $d$ , $z_{o}$ using (3.2) and (3.4)–(3.6).

7. (vii) Obtain corrected $u_{{\it\tau}}/U_{h}$ and repeat previous steps until convergence.

Typically this procedure leads to convergence in just a few iterations.

The model parameters are the von Karman constant ${\it\kappa}=0.4$ , the sectional drag coefficient $C_{d}=1$ (for rectangular prisms), the outer-wake correction strength ${\it\Pi}=0.2$ and the minimum attenuation coefficient $a_{min}=0.4$ . The model inputs include the boundary layer height ${\it\delta}$ , the free-stream velocity $U_{0}$ and the geometric information about the surface roughness elements and their positioning. As it turns out, the predictions on $z_{o}$ and $d$ are independent of ${\it\delta}$ and $U_{0}$ , while $u_{{\it\tau}}$ and $U_{h}$ do depend on these outer-flow conditions.

4.1. Aligned and staggered cubic arrays

In this subsection the model predictions are compared with relatively recent experimental and numerical datasets from Hall, Macdonald & Walker (Reference Hall, Macdonald and Walker1996), Cheng et al. (Reference Cheng, Hayden, Robins and Castro2007), Hagishima et al. (Reference Hagishima, Tanimoto, Nagayama and Meno2009) and Leonardi & Castro (Reference Leonardi and Castro2010) as well as the results from the current LES. The LES code has already been described in § 2. The types of roughness considered in this section are the aligned/staggered arranged cubic roughness elements with various surface coverage densities (see § 2, set I).

The mean velocity profile within and above the roughness canopy is obtained by averaging over one repeating tile (for staggered arrays) or two repeating tiles (for aligned arrays) centred around $x=9{\it\delta}_{0}$ for all cases ( $x=0$ corresponds to the inlet location). The form drag leading to the wall stress from the roughness elements (denoted as ${\it\tau}_{w}$ ) is measured within the simulation at each time step by integrating the pressure over the immersed boundary. The friction velocity is $u_{{\it\tau}}=\sqrt{{\it\tau}_{w}/{\it\rho}}$ . We have neglected the contribution from viscous skin friction, which in all simulated cases is very small. We have checked from the integral wall model (which includes viscous wall stress at the wall as part of the model) that the contributions are in all cases less than 2 % of the form drag on the roughness elements.

Having available the friction velocity $u_{{\it\tau}}$ , we take the mean velocity profiles above the roughness elements and fit a logarithmic law determining the hydrodynamic roughness height $z_{o}$ and displacement $d$ from the regression, as done also in Cheng & Castro (Reference Cheng and Castro2002) and Kanda, Moriwaki & Kasamatsu (Reference Kanda, Moriwaki and Kasamatsu2004). It should be noted that in figure 13 the logarithmic scaling extends down to almost $z=h$ (the first points shown in the plots correspond to the heights of the first LES grid points above $z=h$ ). In order to provide fits that are least affected by possible deviations from logarithmic scaling near the roughness elements, the log laws are fitted between $z=1.5h$ to $z=2.5h$ . A value of ${\it\kappa}=0.4$ is used in this procedure. Comparisons between the simulated mean profile above the cubic roughness ( $z/h>1$ ) and the resulting logarithmic fits (dashed lines) are shown in figure 13. As can be seen, the profiles obtained from the LES follow the log law quite well. We also observe that close to the roughness height any wake correction is negligible, and it is therefore valid to neglect the wake term in (3.4). Finally, we recall that the attenuation coefficient $a$ has been obtained by performing an exponential profile regression on the measured mean velocity profile in the range of heights $0.5<z/h<1$ .

Figure 13. Mean velocity profile above the roughness elements ( $z/h>1$ ) from LES, and the log-law fits. The measured friction velocity $u_{{\it\tau}}$ , displacement height $d$ and hydrodynamic roughness height $z_{o}$ are listed in table 1.

Table 1 lists the relevant quantities determined for each case. The attenuation parameter $a$ for all cases is plotted against ${\it\lambda}_{f}$ in figure 14. It is observed that $a\rightarrow 0.4$ for ${\it\lambda}_{f}\rightarrow 0$ , providing empirical evidence for the result presented in the previous section that for unsheltered conditions $a_{min}\approx 0.4$ .

Figure 14. Measured attenuation coefficient from fits to the LES results plotted against the packing density for all cases. Solid symbols are for aligned cases and hollow symbols are for staggered cases.

Table 1. List of the relevant quantities for each case for cubic roughness of staggered or aligned arrangement. Normalization is via the boundary layer thickness at the inlet and the free-stream velocity. The friction velocity $u_{{\it\tau}}$ is obtained directly in the simulations. The roughness height is $h=0.25{\it\delta}_{0}$ . The boundary layer thickness ${\it\delta}$ is measured at $x=9{\it\delta}_{0}$ downstream of the simulation inlet.

Next, the analytical model is applied to these cases. We mainly study the dependence of various quantities on the solidity ${\it\lambda}_{f}$ . The boundary layer thickness is kept constant in the model, ${\it\delta}/{\it\delta}_{0}=1.3$ (a typical value for the boundary thickness listed in table 1), and $h=0.25$ . A moderate wake correction of ${\it\Pi}=0.2$ is added (Castro Reference Castro2007). The von Karman constant ${\it\kappa}$ is set to ${\it\kappa}=0.4$ , and $C_{d}=1$ throughout. A comparison of the model-predicted $z_{o}$ and $d$ with the experimental and numerical data is shown in figures 15 and 16. It is noted that the data in the literature exhibit high scatter. Compared with relatively early experiments, for example Hall et al. (Reference Hall, Macdonald and Walker1996) (which Macdonald (Reference Macdonald2000) and Coceal & Belcher (Reference Coceal and Belcher2004) have used for model calibration), the hydrodynamic roughness lengths reported by more recent experiments and simulations (Cheng et al. Reference Cheng, Hayden, Robins and Castro2007; Jiang et al. Reference Jiang, Jiang, Liu and Sun2008) tend to be smaller. We follow the usual convention and examine the behaviour of $z_{o}$ and $d$ as functions of the packing density and compare with the data available from the prior studies. Moreover, model predictions for $U_{h}$ and $u_{{\it\tau}}$ are compared with LES data in figure 17. Since they are typically not reported in the literature, only data from the present LES study are included. The model predictions and overall trends agree well with the experimental and simulation data, considering the scatter in the data. The model-predicted $z_{o}$ and $d$ are independent of the boundary layer thickness, following the ‘fully rough’ assumption, yet $u_{{\it\tau}}$ and $U_{h}$ depend on ${\it\delta}$ . The dominant dependence is on the solidity, and dependence on ${\it\delta}$ is very weak.

Figure 15. (a) Comparison of the hydrodynamic roughness length $z_{o}$ predicted by the new model with experimental and simulation data from Hall et al. (Reference Hall, Macdonald and Walker1996), Cheng et al. (Reference Cheng, Hayden, Robins and Castro2007), Hagishima et al. (Reference Hagishima, Tanimoto, Nagayama and Meno2009) and Leonardi & Castro (Reference Leonardi and Castro2010), and with LES data from this study (upright triangles), for cubic roughness. Solid symbols, aligned arrangement; empty symbols, staggered arrangement. The thick and thin lines are the model predictions for the aligned and staggered arrangements respectively. (b) Comparison of the model prediction for $z_{o}$ with (solid line) and without (dashed line) the correction based on the drag partition by Raupach (Reference Raupach1992), for the aligned arrangement (see appendix B).

Figure 16. The same as figure 15 but for the zero-plane displacement $d$ .

Figure 17. Comparison of model-predicted (a) velocity $U_{h}$ at the top of the cubic roughness and (b) friction velocity $u_{{\it\tau}}$ with LES data. Solid symbols, aligned arrangement; empty symbols, staggered arrangement. The thick and thin lines are the model predictions for the aligned and staggered arrangements respectively.

Figure 18 compares the present model with the roughness wall models of Macdonald (Reference Macdonald2000) and Coceal & Belcher (Reference Coceal and Belcher2004). Because these prior models are insensitive to the relative arrangement of roughness elements, and the parameters were chosen to fit the data for staggered cube arrays, we compare only the predictions for staggered arrays. The comparison shows that all models give quite similar predictions, but comparison with figure 15 shows that the present model yields significantly smaller values of $z_{0}$ for the aligned cases, consistent with the data (solid symbols in figure 15).

Figure 18. A comparison of various model predictions for $z_{o}$ , for the case of staggered arrays at different ${\it\lambda}_{f}$ . It should be noted that the models of Macdonald (Reference Macdonald2000) and Coceal & Belcher (Reference Coceal and Belcher2004) do not depend on the relative arrangement of roughness elements, only on ${\it\lambda}_{f}$ .

Finding experimental and simulation data for ${\it\lambda}_{f}>0.4$ is challenging, but the model prediction in densely packed regions can be assessed via asymptotic behaviours. It is expected that when ${\it\lambda}_{f}\rightarrow 1$ , the flow becomes entirely skimming over the roughness elements, approaching an elevated flat-plate turbulent boundary layer. Hence, $z_{o}$ must reduce to 0, while the displacement height $d\rightarrow h$ . Those behaviours can indeed be observed in figure 16. The expected independence of $z_{o}$ and $d$ with respect to ${\it\delta}$ and $U_{o}$ is also satisfied.

4.2. Transverse 2D ribs

In this subsection, the model is applied to transverse (2D) ribs, which, as is well known, can exhibit very different behaviour from 3D roughness elements. The height of the ribs is $0.25{\it\delta}_{0}$ . The roughness solidity studied includes ${\it\lambda}_{f}=0.125$ and $0.25$ . The computational domain is of size $16{\it\delta}_{0}\times 4{\it\delta}_{0}\times 4{\it\delta}_{0}$ and the mesh size is $256\times 64\times 64$ . The velocity profile is averaged over one repeating tile centred around $8{\it\delta}_{0}$ downstream of the inlet. Drag is measured within the simulation, and $z_{o}$ and $d$ are fitted from the measured profile.

Sample instantaneous and mean streamwise velocity fields are shown in figure 19. A comparison between the simulated mean profile above the cubic roughness ( $z/h>1$ ) and the log fit is shown in figure 20. Qualitatively, it can be seen that the vertical rate at which the sheltered region is decreasing downstream between the 2D bars appears to be slower than downstream of 3D cubic objects (comparing, e.g., with figure 3). This trend is consistent with the discussion of expansion rates as a function of roughness element aspect ratio in § 3.3.

Figure 19. Sample instantaneous (a) and mean (b) streamwise velocity field at $z/h=0.5$ for boundary layer flow over transverse square ribs. The velocity is normalized with the free-stream velocity. This case corresponds to ${\it\lambda}_{f}=0.25$ .

Figure 20. Symbols: mean velocity profile above ( $z/h>1$ ) the transverse rib roughness elements from LES. Lines: log-law fit plotted with fitted displacement heights and roughness lengths. The displacement heights are $0.13{\it\delta}_{0}$ and $0.17{\it\delta}_{0}$ , the friction velocities are $0.12$ and $0.11$ (normalized by the free-stream velocity), and the hydrodynamic roughness lengths are $0.053{\it\delta}_{0}$ and $0.030{\it\delta}_{0}$ , for the ${\it\lambda}=0.125$ and $0.25$ cases respectively. The rib height is $0.25{\it\delta}_{0}$ .

Following prior practice, results are first compared as a function of the distance between the centres of two transverse square ribs, ${\it\lambda}=L_{x}/h+1$ . Results are shown as a function of the additive constant $B$ in the expression $U(z)/u_{{\it\tau}}=1/{\it\kappa}\log [(z-d)/h]+B$ . Figure 21 shows the comparison between the model and extensive datasets.

Figure 21. Comparison of the additive constant $B$ from the model prediction ( $B={\it\kappa}^{-1}\ln (h/z_{o})$ ) and from experiments and numerical simulations (Cui et al. Reference Cui, Patel, Lin and Patel2000; Leonardi et al. Reference Leonardi, Orlandi, Smalley, Djenidi and Antonia2003; Coleman et al. Reference Coleman, Nikora, McLean and Schlicke2007; Ikeda & Durbin Reference Ikeda and Durbin2007), for transverse ribs. Two points are from the LES of this study. The height of the transverse square rib is ${\approx}\!1/5$ of the boundary layer height.

In order to provide a comparison with the results for cubic roughness elements, in figure 22 the model prediction for 2D ribs is compared with the model predictions for aligned and staggered cubes, shown before in figure 15. Overall, the shape and order of magnitude of the results is comparable, although the peak occurs at lower ${\it\lambda}_{f}$ for the 2D ribs (sheltering effects still occur at larger distances between the elements, i.e. smaller ${\it\lambda}_{f}$ , due to the decreased spreading rate in the 2D case).

Figure 22. Comparison of the roughness length $z_{o}$ from the model prediction for aligned and staggered cubes, as well as 2D ribs.

4.3. Rectangular roughness with non-uniform bimodal height distributions

The cases presented in the previous sections considered roughness elements of uniform height. The next level of complexity is to consider roughness elements with non-uniform height distributions. Varying-height distributions were studied in experiments as described in Cheng & Castro (Reference Cheng and Castro2002), Hagishima et al. (Reference Hagishima, Tanimoto, Nagayama and Meno2009) and Zaki et al. (Reference Zaki, Hagishima, Tanimoto and Ikegaya2011). Jiang et al. (Reference Jiang, Jiang, Liu and Sun2008) studied the behaviour of the hydrodynamic roughness length $z_{o}$ by systematically varying the degree of roughness height non-uniformity while keeping the solidity ${\it\lambda}_{f}=0.11$ . The results showed that $z_{o}$ increases monotonically with the standard deviation of the roughness height distribution, although it is not known whether this trend would hold for roughness at a different solidity, or with different roughness arrangements.

The goals of this subsection are twofold. First, to provide LES data on rectangular-prism roughness with non-uniform heights for various solidities and different arrangements (staggered and aligned). For each case the hydrodynamic roughness length $z_{o}$ and friction velocity $u_{{\it\tau}}$ are measured from the LES. Second, we aim to compare the measured values with the predictions from the analytical model described in § 3. As a first step in considering non-uniform heights, we consider bimodal height distributions in which the mean height $h_{m}$ of the elements is kept constant but the spread around it, quantified using the standard deviation of the element height, is varied. It should be noted that the roughness height is either $h_{h}=h_{m}+{\it\sigma}_{h}$ or $h_{l}=h_{m}-{\it\sigma}_{h}$ . The results of roughness with Gaussian height distribution are presented in § 4.4. The set-up for these cases (denoted as LXXS/A-StdXX) was already described in § 2. The repeating tiles for the rough walls are presented in figure 23.

Figure 23. Repeating tiles for the cases considered. Here, L stands for ${\it\lambda}_{f}$ , the two digits following represent $100{\it\lambda}_{f}$ , the last letter S is for ‘staggered’ and A is for ‘aligned’: (a) L06S, (b) L11S, (c) L11A, (d) L25S. The roughness is coloured with black/grey. Black is for elements of height $h_{h}=h_{m}+{\it\sigma}_{h}$ and grey is for height $h_{l}=h_{m}-{\it\sigma}_{h}$ . The subscript ‘ $h$ ’ stands for ‘high’, ‘ $l$ ’ stands for ‘low’ and the mean height is $h_{m}=(h_{h}+h_{l})/2$ .

The mean velocity profile above the roughness is obtained via averaging the temporally averaged flow field in the spanwise $y$ direction, and in the streamwise direction $x$ between two repetitive tiles starting from $x=4{\it\delta}_{0}$ . Figure 24 shows the resulting set of velocity profiles for all cases. The log law is fitted between $z=h_{h}+0.5h_{m}$ and $z=h_{h}+2h_{m}$ . To reduce the uncertainty in the log-law fitting, we use the displacement height determined from the rough-wall model ((3.6), and note that in this case $A_{f}$ is not a constant but a function of the wall normal distance). The values obtained from the model are listed in tables 2–5. Figure 25 shows the velocity profiles in a log-linear scale. The measured (fitted) values for the hydrodynamic roughness length $z_{o}$ , friction velocity $u_{{\it\tau}}$ and velocity at the top of the roughness canopy $U_{h}=U(z=h_{h})$ are listed in tables 2–5.

Figure 24. The mean velocity profiles for all of the cases: (a) Lf06S, (b) Lf11S, (c) Lf11A, (d) Lf25S. The indicated standard deviation of the roughness elements is normalized with the mean roughness height.

Figure 25. The log-law-fitted velocity profile for all of the cases: (a) Lf06S, (b) Lf11S, (c) Lf11A, (d) Lf25S. Only data above $z>h_{h}+0.5h_{m}$ are shown. The standard deviation of the roughness elements is normalized with the mean roughness height.

Table 2. Displacement height $d$ (from model), hydrodynamic roughness length $z_{o}$ , friction velocity $u_{{\it\tau}}$ and velocity at the top of the roughness canopy $U(z=h_{h})$ for case Lf06S. Here, Std stands for the standard deviation, length is normalized with $h_{m}$ and velocity is normalized with the free-stream velocity $U_{0}$ .

Table 3. The same as table 2 for case Lf11S.

Table 4. The same as table 2 for case Lf11A.

Table 5. The same as table 2 for case Lf25S.

We then compare the LES results with the model predictions. The analytical model requires careful evaluation of $h_{s}$ (to obtain $a$ , etc.), and some intermediate results of the geometric calculations are shown in appendix A. Figure 26 shows a comparison of the values of $z_{o}$ , $U_{h}$ and $u_{{\it\tau}}$ predicted by the analytical model with those measured from LES. We observe that the analytical model captures the measured results quite well over a significant range of parameters.

Figure 26. Comparison of the analytical model predictions of $z_{o}$ , $d$ , $U_{h}$ and $u_{{\it\tau}}$ as a function of the height standard deviation (line) with the LES measurements (symbols). The LES measurements are denoted with symbols: L06S, solid squares; L11S, empty squares; L11A, empty circles; L25S, solid circles. The lines are model predictions: L06S, solid line; L11S/A, dashed line (no difference is observed in the model predictions for L11S and L11A for this roughness configuration); L25S, dot-dashed line.

Several observations can be made about the model. First, both $z_{o}$ and $u_{{\it\tau}}$ increase with the standard deviation of the roughness height. An up to threefold increase in $z_{o}$ can be observed for Lf25S from $\text{std}(h/h_{m})=0$ to $\text{std}(h/h_{m})=1$ . Second, while $z_{o}$ and $u_{{\it\tau}}$ increase steadily for Lf06S and Lf11S as $\text{std}(h/h_{m})$ increases, for Lf25S a faster increase in $z_{o}$ and $u_{{\it\tau}}$ is observed at low $\text{std}(h/h_{m})$ . A similar rate of increase with the standard deviation is observed at high standard roughness deviation for all three sets of cases considered. Third, an increase in $U_{h}$ is observed at low standard deviation of roughness height, and then the value stays almost constant, although the exact value varies from case to case. The increase in $U_{h}$ at low $\text{std}(h/h_{m})$ could due to the fact that we have defined $U_{h}$ at $U(z=h)$ . Because of height non-uniformity, at $h_{h}$ the planar area density ${\it\lambda}_{p}$ is only half of that for the case when roughness has uniform height.

Next, a qualitative explanation is attempted to elucidate the trends in the changing slopes of $z_{0}$ and $d$ as functions of $\text{std}(h/h_{m})$ . It is postulated that the change in behaviour indicates a change in roughness regime, changing from ‘ $d$ ’-type roughness to ‘ $k$ ’-type roughness at larger $\text{std}(h/h_{m})$ . Consider figure 27, which shows several vertical plane cuts from cases L06S-Std00, L06S-Std50, L25S-Std00 and L25S-Std50. A $d$ -type roughness behaviour is observed in L25S-Std00. For ‘ $d$ ’-type roughness, the flow above the roughness skims over the roughness elements and as a result less drag is generated. As the standard deviation in the roughness height increases, the roughness for L25S clearly changes from ‘ $d$ ’-type roughness to ‘ $k$ ’-type roughness (figure 27 d), allowing the flow above to impinge some of the roughness elements. This transition in flow regime is consistent with a more rapid increase in $u_{{\it\tau}}$ and $z_{o}$ at low $\text{std}(h/h_{m})$ . Once $\text{std}(h/h_{m})$ is above some value (for the cases studied here it seems to be approximately $\text{std}(h/h_{m})\sim 0.3$ ), the increase is less rapid since the entire behaviour follows ‘ $k$ ’-type behaviour without the rapid increase in drag associated with the transition from ‘ $d$ ’-type to ‘ $k$ ’-type roughness. For low initial solidities, even the cases with uniform roughness height are already in the ‘ $k$ ’-type regime and thus display a more uniform increase with standard deviation.

Figure 27. Vertical plane cuts of the averaged streamwise velocity from (a) L06S-Std00, (b) L06S-Std50, (c) L25S-Std00 and (d) L25S-Std50 through the middle plane of the roughness.

4.4. Rectangular roughness with Gaussian distribution

In this subsection, the model predictions are compared with the LES measurements of rectangular roughness with Gaussian distribution (set II in § 2.2). Both the mean roughness height and the surface coverage (11 %) are kept constant. The standard deviation in roughness height is varied according to values ${\it\sigma}_{h}/h=0.24$ , $0.35$ and $0.50$ . For each ${\it\sigma}_{h}$ , four realizations of rough walls are randomly generated for the LES runs. The LES set-up has been described in § 2.1.

The temporally averaged velocity is further averaged within a streamwise section of length $6h_{m}$ centred at $x=24h_{m}$ to obtain the mean velocity profile. Figure 28 shows the mean profiles for all of the cases considered in this subsection. Very little difference is observed between cases with the same roughness height distribution. The log law is fitted between $z=1.5h_{m}+{\it\sigma}_{h}$ and $z=2.5h_{m}+{\it\sigma}_{h}$ using the displacement $d$ from the rough-wall model (see below) and the von Karman constant ${\it\kappa}=0.4$ . The measured values of $z_{o}$ , $u_{{\it\tau}}$ and the boundary layer thickness ${\it\delta}$ are listed in table 6. Figure 29 compares the log law and the fitted profiles.

Figure 28. The mean streamwise velocity profile for $\text{std}(h/h_{m})=0.24,0.35,0.50$ (a–c). For each ${\it\sigma}_{h}$ , four LES consisting of four realizations of randomly generated surfaces are conducted.

Figure 29. A comparison of the log law and the fitted mean streamwise velocity profile for all cases. From (a–c) $\text{std}(h/h_{m})=0.24,0.35,0.50$ .

Table 6. Hydrodynamic roughness height $z_{o}$ , displacement height $d$ , friction velocity $u_{{\it\tau}}$ and 99 % boundary layer thickness ${\it\delta}$ for roughness with different height variations. Here, $z_{o}$ , $d$ and ${\it\delta}$ are normalized with the mean roughness height $h_{m}$ and the friction velocity is normalized with the free-stream velocity $U_{0}$ .

Next, to obtain the model predictions for $z_{o}$ , $d$ , $u_{{\it\tau}}$ and $U_{h}$ , we generate 512 realizations (sufficiently many to obtain converged statistics) of rough surfaces for a large range of roughness height variances $\text{std}(h/h_{m})$ between 0 and 0.5. We apply the sheltering model to each surface. Geometrically computing the sheltered frontal area is non-trivial for such complex surfaces, and therefore in this case we must use a numerical code to obtain $A_{s}$ . The method is as described in more detail in appendix A. Figure 30 is a graphical visualization of the modelled sheltering regions for one particular realization with ${\it\sigma}_{h}=0.5h$ . We then average the model predictions for $z_{o}$ , $d$ , $u_{{\it\tau}}$ and $U_{h}$ and plot the results as solid lines in figure 31.

Figure 30. Visualization of the sheltering regions among the roughness elements with Gaussian height distribution in a realization, with ${\it\sigma}_{h}=0.5h_{m}$ . Periodicity is assumed in the spanwise ( $y$ ) and streamwise ( $x$ ) directions.

Figure 31. A comparison of the model predictions and the LES measurements. Here, Mdl stands for model and the symbols stand for LES results. For each symbol, the largest deviation from the mean value observed in the four LES with the same $\text{std}(h/h_{m})$ is shown as the extension on either side of the error bar.

Figure 31 compares the model predictions and the mean values from the four LES realizations for each value of $\text{std}(h)$ . The boundary layer thickness is set to be ${\it\delta}=6.4h_{m}$ in the model, and a wake correction of ${\it\Pi}=0.5$ is used. As can be seen, the model predictions agree well with the LES data. The uncertainty in $z_{o}$ and $u_{s}$ due to the randomness in the roughness height distribution is not strong. We observe that with an increase in the roughness height variation, $z_{o}$ increases considerably, while $u_{{\it\tau}}$ and $U_{h}$ increase but only by a small amount.