1. Introduction
The problem of predicting the friction drag exerted by turbulent flow on a surface based on knowledge about the surface geometry has received extensive and enduring attention. Efforts on this topic date back to the early experimental work of Colebrook (Reference Colebrook1939) and Nikuradse (Reference Nikuradse1950). Since 1944 the Moody diagram that relates the friction drag with the equivalent sandgrain roughness height $k_{s}$ (Moody Reference Moody1944) has been the most commonly employed engineering tool. There have been many further efforts since then to correlate the surface topology with the hydrodynamic response through equivalent sandgrain roughness height, friction factor (Simpson Reference Simpson1973), effective slope (Napoli, Armenio & De Marchis Reference Napoli, Armenio and De Marchis2008; Schultz & Flack Reference Schultz and Flack2009), surface skewness factor (Flack & Schultz Reference Flack and Schultz2010) and examinations of Reynolds number similarity for rough walls (see, e.g., Flack, Schultz & Shapiro Reference Flack, Schultz and Shapiro2005). For reviews, see Grimmond & Oke (Reference Grimmond and Oke1999) and Flack & Schultz (Reference Flack and Schultz2010).
Researchers working on modelling urban canopy flows have been interested in parameterizations of the hydrodynamic roughness length ( $z_{o}$ ), displacement height ( $d$ ) and drag coefficient for the specific case of rectangularprism roughness elements, due to the typically cubiform shapes of buildings. Some of the extensive efforts have been reviewed in Grimmond & Oke (Reference Grimmond and Oke1999) and Barlow & Coceal (Reference Barlow and Coceal2009). For particular applications of various models for $z_{o}$ and $d$ , see Arnfield (Reference Arnfield2003), Vardoulakis et al. (Reference Vardoulakis, Fisher, Pericleous and GonzalezFlesca2003) and Kanda (Reference Kanda2006). Among the parameters expected to be the most important for $z_{o}$ and $d$ are the solidity ${\it\lambda}_{f}$ (defined as the projected frontal area per unit lot area), the planar density ${\it\lambda}_{p}$ (defined as the projected horizontal area per unit lot area) and the characteristic height $h$ of individual roughness elements (Grimmond & Oke Reference Grimmond and Oke1999). Morphometric models usually include explicit dependences of, e.g., $z_{o}$ on these parameters determined through empirical approaches such as fitting with experimental or numerical data. Models that fall into this category include those introduced in Lettau (Reference Lettau1969), Counehan (Reference Counehan1971), Macdonald, Griffiths & Hall (Reference Macdonald, Griffiths and Hall1998) and Kim et al. (Reference Kim, Lee, Joo, Ryu, Kim, You and Shim2011). Calibrated for certain surface roughness and over a given range of ${\it\lambda}_{f}$ , ${\it\lambda}_{p}$ , such models can provide reasonably accurate predictions of $z_{o}$ and $d$ for practical applications. Nevertheless, to a great extent they remain dependent upon much empirical input. Thus, the problem of relating hydrodynamic and geometrical roughness remains open, since local flow conditions may depend on details of the roughness elements in a highly casespecific fashion.
In order to develop physicsbased models for flows over surfaces with attached roughness elements, some detailed understanding of the averaged velocity profile within the roughness layer (defined here as the region between the surface and the top of the roughness elements) must be developed. This is similar to the situation where knowledge about the logarithmic law has led to physically based drag laws for smooth boundary layers. One option is to examine the differential momentum equation, as is done in differential urban canopy models (e.g. 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) which focus on predicting the horizontally ( $x{}y$ ) and temporally averaged velocity profile inside and above the canopy using integration of the momentum differential equation. To obtain the mean velocity profile as a function of height, the modelling task focuses on the Reynolds stress (which arises due to temporal averaging over turbulence), the dispersive stress (which arises due to spatial averaging of the mean velocity across spatial heterogeneities in the temporal mean velocity distribution) and the form drag (which arises due to the direct momentum extraction by the roughness elements interacting with the flow). The Reynolds stress is usually modelled with a Prandtl mixinglength eddyviscosity model, the dispersive stress is commonly neglected and the drag is typically modelled using a quadratic law for a body force associated with the form drag $F=C_{d}{\it\rho}U^{2}{\it\delta}A_{f}/{\it\delta}V$ , where ${\it\rho}$ is the fluid density, $C_{d}$ is the drag coefficient, $U$ is a mean streamwise reference velocity and ${\it\delta}A_{f}$ is the projected frontal area within a volume ${\it\delta}V$ . The mixing length (denoted as $l_{m}$ ) and the drag coefficient $C_{d}$ need to be specified. With these models for stresses and the force, the spatially and temporally averaged streamwise momentum equation can be integrated in the wall normal direction to obtain the mean velocity profile.
To increase the accuracy of the predictions from such models, various empirical inputs have been employed by different authors. Examples include the approach of Macdonald (Reference Macdonald2000), requiring an attenuation factor to be specified a priori as a function of the solidity ${\it\lambda}_{f}=A_{f}/A_{T}$ (where $A_{f}$ and $A_{T}$ are the projected frontal and horizontal lot areas respectively), while Coceal & Belcher (Reference Coceal and Belcher2004) employ an empirical function to model the displacement height $d$ . As mentioned before, it is quite well established that the solidity ${\it\lambda}_{f}$ is the most important parameter characterizing the surface morphology, and most currently available roughwall models are insensitive to more detailed characterizations of the roughness distribution (see the reviews by Grimmond & Oke Reference Grimmond and Oke1999; Barlow & Coceal Reference Barlow and Coceal2009). While typically the wake interactions are not explicitly modelled, there have been some attempts (e.g. Jiménez Reference Jiménez2004; MillwardHopkins et al. Reference MillwardHopkins, Tomlin, Ma, Ingham and Pourkashanian2011) to model the mutual sheltering between roughness elements.
Thus, while significant progress has been achieved in roughness modelling over the past few years, shortcomings in the models reviewed above can be identified. Morphometric and urban canopy models depend significantly on empirical input such as ad hoc specifications of the mixing length $l_{m}$ on ${\it\lambda}_{f}$ , the height from the bottom surface, etc. Moreover, although differential urban canopy models can make more detailed predictions than morphometric models, being differential instead of algebraic (e.g. Coceal & Belcher Reference Coceal and Belcher2004; Harman & Finnigan Reference Harman and Finnigan2007), they are more costly to evaluate. This can be an obstacle when attempting to combine these models with numerical weather prediction codes. Lastly, mutual sheltering among roughness elements, while being a commonly accepted concept, lacks a clear operational definition, and there is still no simple yet accurate model that can include mutual sheltering in the context of a practical roughness model.
In this study we attempt to resolve some of these difficulties by adopting the von Karman–Pohlhausen integral approach (Pohlhausen Reference Pohlhausen1921), in which a functional form for the mean velocity is assumed, including parameters that must be obtained from physical constraints. At the simplest level, there have been expectations that an exponential profile occurs within the roughness layer. This expectation can be motivated (Cionco Reference Cionco1966; Macdonald Reference Macdonald2000; Coceal & Belcher Reference Coceal and Belcher2004) by writing the Reynoldsaveraged streamwise momentum equation in which the vertical (direction $z$ ) gradient of the momentum flux (modelled using a constantmixinglength eddy viscosity) balances the distributed drag from roughness elements:
Here, $\text{d}A_{f}/\text{d}V$ is the projected frontal area per unit volume, and $l_{m}$ and $C_{d}$ are the mixing length and drag coefficient respectively. These quantities can in principle be $z$ dependent, but at this stage they are assumed to be constant in order to obtain a simple solution, whose validity must then be tested empirically. It can be readily seen that an exponential function of $z$ , $U(z)\sim \exp (z)$ , solves (1.1):
where $h$ represents the height of the roughness elements, $U_{h}$ is the velocity at $z=h$ and $a$ is an attenuation factor (which depends upon the parameters in (1.1)).
Our first goal is to establish additional empirical evidence to confirm or refute the expectation of an exponential velocity profile under more realistic conditions in which some of the underlying assumptions used in (1.1) may not hold exactly. Therefore, as a first step we determine the mean flow behaviour within the roughness layer via a series of largeeddy simulations (LES) of turbulent boundary layers over arrays of wallattached rectangularprism roughness elements. Various height distributions, arrangements and surface coverage densities are considered. The results will be shown to support a generic form for the velocity profile within the roughness layer with exponential behaviour, (1.2). This observation motivates us to employ an exponential shape function for the velocity profile with free parameters ( $a,U_{h},\ldots$ ) that depend on the local surface morphometric properties and flow configurations in the roughness sublayer. Since we have used the shape function in place of the momentum equation within the roughness sublayer, the model is algebraic rather than differential. A similar approach has been used recently to develop a new wall model for LES (Yang et al. Reference Yang, Sadique, Mittal and Meneveau2015). In order to determine the most important profile parameter, namely the attenuation coefficient $a$ in the exponential function, an additional model for mutual sheltering effects among roughness elements must be included. While the model is developed for possibly more general rough surfaces, it is motivated, validated and applied here only for the specific, but important, case of 3D and 2D rectangular roughness elements with uniform and nonuniform height distributions. As mentioned before, such surfaces are commonly found in urban canopies but also in flow over engineered surfaces such as electronic circuit boards, etc. In developing an analytically tractable model for hydrodynamic roughness, we are motivated by the continued relevance of such models when computational flow predictions spanning multiple scales and levels are required. For example, in the newly introduced integral wall model for LES (Yang et al. Reference Yang, Sadique, Mittal and Meneveau2015), for highReynoldsnumber flow over rough surfaces in which the roughness falls much below the first grid point affordable via LES near the surface, an analytical model for the roughness length $z_{0}$ must still be provided to the (nonequilibrium) wall model. Modelling of geophysical flows also requires such parameterizations, while engineering design tools require the ability to cover many cases rapidly, often precluding numerical simulations.
The remainder of the paper is organized as follows. Numerical simulations (LES) that resolve the roughness elements within the roughness layer are used to determine the mean velocity profile in the roughness layer in § 2. Motivated by the simulation results, the assumed shape function for the velocity profile is briefly described in § 3 with various constraints that must be satisfied in order to build an analytical model. The shape function provides a horizontally averaged description of the flow field, yet some additional details about the flow within the roughness sublayer in between roughness elements must also be taken into account to build the model. The flow within the roughness sublayer is examined and modelled in § 3.2, where a geometrybased sheltering model is proposed. Detailed comparisons of the model predictions with experimental and numerical data for 3D cubic roughness elements and 2D elements (bars) are presented in § 4. The model predictions for roughness with a nonuniform height distribution are compared with LES results in §§ 4.3 and 4.4. In order for the model to display the correct asymptotic behaviour at low surface coverage, a correction for drag partitioning is included in appendix B. Conclusions are given in § 5.
2. Mean flow profiles within the roughness layer from LES
2.1. Simulation setup
We use the inhouse finitedifference code Vicar3D to solve the incompressible filtered Navier–Stokes equations, including the Vreman subgridscale model (Vreman Reference Vreman2004; You & Moin Reference You and Moin2007) and the iWMLES model for the nearwall closure (Yang et al. Reference Yang, Sadique, Mittal and Meneveau2015). Details of the code and the LES setup 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 subgridscale 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 subgridscale modelling than, say, secondorder 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 roughwall rescaling–recycling technique (Yang & Meneveau Reference Yang and Meneveau2016). The boundary layer thickness at the inflow, ${\it\delta}_{0}$ , is $4h_{m}$ (for singleheight roughness cases) to $6h_{m}$ (for varyingheight 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 singleheight 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 freestream 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 freestream 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 subgridscale roughness. A total simulation time of 100 flowthrough 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 flowthrough 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
Largeeddy 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 nonuniform 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 probabilitydensity 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/SX, 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 L11Std50A1 (figure 3 a), and instantaneous and averaged streamwise velocities on a plane at height $z=h_{m}$ (figure 3 b,c), for flow case L11Std501 in roughwall set II.
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.
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.
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 nearwall 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).
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 welldefined length scale, each with reasonably uniform crosssectional area, and welldefined flow separation points.
3. An analytical roughwall flow model
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:
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 freestream 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 nearwall 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 roughwall 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 tangentlike profiles are used to model this mixinglayertype 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 mixinglayertype 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:
where the integration is performed on the projected frontal area $A_{f}$ within the lot area $A_{T}$ . For the righthand 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 roughwall 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 righthand 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 rectangularprism roughness elements, integrating the exponential velocity profile between $z=0$ and $z=h$ , (3.2) leads to
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
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 roughwall 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 freestream velocity $U_{0}$ :
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
For singleheight rectangularprism roughness after integration of the exponential profile between $z=0$ and $z=h$ , (3.6) simply leads to
Combining (3.3), (3.4) and (3.7), $z_{o}$ for singleheight rectangularprism roughness can be expressed as
Replacing (3.8) into (3.4) and (3.5), we obtain the expression for $u_{{\it\tau}}$ :
and
where $d$ is obtained from (3.7). Again, it should be noted that (3.7)–(3.10) are valid only for singleheight rectangularprism 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.
3.2. Volumetric sheltering
In this section, we model the reduction in the momentum in the wakes of rectangularprism 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 rectangularprism 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
where for $C_{DH}$ we use the drag coefficient for a surfacemounted 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
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
leading to
The solution is simply
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 rectangularprism element that has reducedmomentum fluid. Wake expansion in the vertical direction can bring in fluid with high momentum; a reduction in the volume of lowmomentum fluid would thus be expected. Wake expansion in the spanwise direction, on the other hand, expands the volume of lowmomentum fluid. This configuration is physical for nearwall roughness with aspect ratio $h/w$ not large. Naturally, we expect the wake be eaten up from the side if the roughness geometry is ‘sticklike’, 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
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 rectangularprism 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}(yy_{1a})\mathscr{H}(yy_{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.
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
where $L_{x}$ is the horizontal distance between roughness elements. Substitution of (3.17) into (3.15) leads to
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, nonuniform 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 rectangularprism 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
We replace $\tan {\it\theta}=C_{{\it\theta}}(u_{{\it\tau}}/U_{h})$ and solve for $C_{{\it\theta}}$ to obtain
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 onethird 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 ‘sticklike’ roughness elements (e.g. a canopy of long wallattached 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}}$ .
3.4. Summary of the wall and sheltering model
We briefly summarize the roughwall model and outline how it can be applied to rough walls mounted with rectangularprism 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 zeroplane 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 rectangularprism roughness elements can be written as
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 singleheight singleaspectratio roughness. We need (3.21) mainly to account for roughness elements with nonuniform 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.

(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})$ .

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

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

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

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

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

(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 outerwake 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 freestream 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 outerflow conditions.
4. Applications
In this section, we report results of LES of flow over (1) aligned and staggered cube arrays, (2) 2D transverse ribs, (3) rectangular roughness with bimodal height distribution and (4) rectangular roughness with Gaussian height distribution. The purposes of this section are twofold: first, we provide data on roughness distributions that have not been fully studied; second, we compare measurements of $z_{o}$ , $d$ , $u_{{\it\tau}}$ and $U_{h}$ from LES with the predictions from the sheltering model.
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$ .
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$ .
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 modelpredicted $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 modelpredicted $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 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).
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 flatplate 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.
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 [(zd)/h]+B$ . Figure 21 shows the comparison between the model and extensive datasets.
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).
4.3. Rectangular roughness with nonuniform 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 nonuniform height distributions. Varyingheight 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 nonuniformity 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 rectangularprism roughness with nonuniform 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 nonuniform 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 setup for these cases (denoted as LXXS/AStdXX) was already described in § 2. The repeating tiles for the rough walls are presented in figure 23.
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 loglaw fitting, we use the displacement height determined from the roughwall 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 loglinear 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.
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.
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 nonuniformity, 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 L06SStd00, L06SStd50, L25SStd00 and L25SStd50. A $d$ type roughness behaviour is observed in L25SStd00. 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.
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 setup 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 roughwall 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.
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 nontrivial 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 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.
5. Conclusions
In this study, the mean velocity distribution within the roughness layer of turbulent boundary layer flow over arrays of rectangular prisms is examined. It is found that for a wide range of element placement morphologies, a generic exponential velocity profile with respect to the wall normal distance is a very good description of the mean velocity within the upper 70–80 % portion of the roughness layer. To model the hydrodynamic effects of rough walls, a shape function for the velocity profile is proposed as a replacement for a mixinglength closure and integration of the temporally/spatially averaged momentum equation. There are five unknowns in the shape function: the hydrodynamic roughness length $z_{o}$ , the displacement height $d$ , the friction velocity $u_{{\it\tau}}$ , the velocity at the top of the canopy $U_{h}$ and an attenuation coefficient $a$ . Four constraints, including the momentum balance, are available from first principles. In addition, a geometric sheltering model is developed to provide the fifth condition. Moreover, the drag coefficient $C_{d}$ is set to unity and the sheltering expansion rate is set to $C_{{\it\theta}}u_{{\it\tau}}/U_{h}$ , with a coefficient $C_{{\it\theta}}$ that is dependent on the roughness element aspect ratio. Differently from earlier analytical roughness models, the model developed in this study is responsive to the roughness distribution because the model is coupled with a geometric sheltering model. The model is applied to cubic roughness distributions of various solidities, to the case of arrays of transverse 2D square ribs and to roughnesses with nonuniform height distributions. The model predictions compare well with experimental and numerical datasets from other authors, and with the LES results from this study. Correct asymptotic behaviours are obtained at both ${\it\lambda}_{f}\rightarrow 1$ and ${\it\lambda}_{f}\rightarrow 0$ (in the latter case including a correction via drag partition described in appendix B). The sheltering model developed here can also be responsive to additional variations in the spatial roughness distribution. Comparisons of the model predictions with rectangular roughnesses oriented at angles (i.e. nonfrontal) with the flow, or different incident flow directions and arrays with elements displaced in the spanwise direction for arbitrary distances will be presented in a separate report.
The analytical model is, for now, restricted to roughness with rectangularprism shape, since the sheltering region within the wake of such objects can be clearly related to the frontal crosssection of the object. For more general shapes, such as surfaces covered with LEGO blocks (Placidi & Ganapathisubramani Reference Placidi and Ganapathisubramani2015; Vanderwel & Ganapathisubramani Reference Vanderwel and Ganapathisubramani2015), as well as with as cones, semihemispheres, frustums, etc., where the flow separation point may not be easily identified, further generalizations and detailed comparisons with data are required to establish wider generality and applicability of the basic ingredients of the model. Moreover, inclusion of spatially changing roughness configurations and effects of nonequilibrium conditions (Schultz, Schatzmann & Leitl Reference Schultz, Schatzmann and Leitl2005) would be of considerable interest.
Acknowledgements
The authors wish to thank Professor M. Schultz for useful discussions on the topic of turbulent flow over rough surfaces, and the Office of Naval Research (grant number N000141210582, Dr R. Joslin, programme director) for financial support. Simulations were performed using the DoD system.
Appendix A. Further examples of sheltered area evaluations
We have already considered aligned cube arrays in § 3.2. Here, we provide an additional two examples of calculating the sheltered area analytically for regular simple roughness arrangements. Then, we briefly discuss how to implement a simple numerical code to perform this geometrical calculation for more complicated cases.
First, we consider fully staggered cube arrays. Because all roughness elements will be equally sheltered, we only need to consider one element. We begin by identifying the upstream elements that could shelter the particular element under consideration. Figure 32 sketches the interactions that need to be considered. The sheltering due to $B_{1}$ is $S_{B_{1}A}=h_{B_{1}A}[\mathscr{H}(y)\mathscr{H}(yw)]$ , where $h_{B_{1}A}=\max [0,hl_{x}\tan ({\it\theta})]$ , $l_{x}=2h/\sqrt{{\it\lambda}_{f}}h$ , $w=h$ . The sheltered area due to $B_{2}$ is $S_{B_{2}A}=h_{B_{2}A}[\mathscr{H}(y)\mathscr{H}(yw_{2})]$ , where $h_{B_{2}A}=\max [0,h\tan ({\it\theta})(l_{x}h)/2]$ , $w_{2}=\max [\tan ({\it\theta})(l_{x}h)/2l_{y},0]$ , $l_{y}=h/\sqrt{{\it\lambda}_{f}}h$ . Because of symmetry, the sheltering due to $B_{3}$ is $S_{B_{3}A}=h_{B_{2}A}[\mathscr{H}(yw+w_{2})\mathscr{H}(yw)]$ . The sheltered area of $A$ is then determined by $A_{s,A}=\int _{0}^{w}\max [S_{B1A},S_{B2A},S_{B3A}]\,\text{d}y=(w2w_{2})h_{B_{1}A}+2w_{2}h_{B_{2}A}$ .
Second, we consider roughness with staggered bimodal height distribution. The interactions that need to be considered are sketched in figure 33(a). This time, we need to calculate the sheltered area for both the lowerrising and higherrising elements separately. For the higherrising elements $A_{1}$ , these could be sheltered by the lowerrising roughness upstream $B_{2}$ and by the higherrising element further upstream $A_{3}$ , as well as by the one that is diagonally upstream $A_{2}$ . One can notice that the interactions among $A_{1}$ , $A_{2}$ and $A_{3}$ are quite like the interactions among the staggered arrays. The interactions between $B_{2}$ , $A_{1}$ and $A_{2}$ , $B_{1}$ , on the other hand, are just like aligned cubes.
The flow sheltering in the canopy layer is sketched in figure 33(b). The sheltered projected frontal area in the example sketched in figure 33 is $A_{f,s}/A_{T}={\it\lambda}_{f}(h_{B}w+h_{s,A}w)/(h_{B}w+h_{A}w)$ , where $h_{A}=h_{m}+\text{std}(h)$ , $h_{B}=h_{m}\text{std}(h)$ , and $h_{s,B}$ and $h_{s,A}$ are the heights of the sheltered areas for the low and highrising roughness respectively. The height of the equivalent sheltered layer $h_{s}$ is then simply $h_{s}=(h_{s,A}+h_{B})/2$ .
In general, the equivalent sheltered layer height is
To determine $h_{s,B}$ and $h_{s,A}$ , the wake interaction among the roughness elements needs to be considered. In general, the higherroughness–higherroughness, higherroughness–lowerroughness and lowerroughness–lowerroughness interactions need to be considered, but since lowerroughness elements are separated by higherroughness ones, there is no need to consider lowerroughness–lowerroughness interactions. With $u_{{\it\tau}}/U_{h}$ known as a function of $a$ , $h_{s}$ is expressed with $a$ via (A 1). Equation (3.15) can be then used to solve $a$ . With $a$ known, $U_{h}$ , $u_{{\it\tau}}$ , $z_{o}$ and $d$ can again be solved via (3.2) and (3.4)–(3.6).
For more complex cases, such as the example with a Gaussian distribution of heights, the geometric calculations must be performed by a code. The method is based on discretizing the upstream and downstream edges of the base of roughness elements. A sketch is provided in figure 34. Because of the rectangular shape of the roughness, only the base plane needs to be considered. First, for all rectangular roughness elements, the windward faces are identified, i.e. their projection on the ground, a segment of line. These lines will be called ‘receivers’ (because they receive sheltering). The height $h$ of the receiver element is associated with it in order to make sure that sheltering does not go beyond the height of the sheltered roughness. Second, we identify all leeward faces, and group their projections on the ground into an ‘emitters’ set. We also keep track of their heights to calculate the sheltering height to any downstream roughness element (if the downstream (receiver) element is ${\it\delta}x$ downstream, the sheltering height on the receiver segment is $h_{s}=h{\it\delta}x\tan ({\it\theta})$ ). Both receiver and emitter lines are discretized (here we use 100 points per line).
The task now is to loop through all points in each member of the ‘receivers’ set and calculate how it is sheltered by each member in the set of ‘emitters’. To calculate how a particular ‘emitter’ member S is sheltering a particular point P in a line that belongs to the ‘receiver’ set, we take the minimum between the height of the volumetric sheltering of S at P and the height of element P. This height is then compared among sheltering by all ‘emitter’ members and the maximum value of the sheltering height at P is the sheltering height at P. By doing this, the sheltering height for each of the ‘receiver points’ can be calculated and an integral of sheltering height at each point across the receiver line leads to $A_{s}$ .
For example, consider the point P in figure 34. Within the searching domain (the domain enclosed by the dashed box which is $3h_{max}U_{h}/u_{{\it\tau}}$ upstream and on both sides), there are four ‘sender’ members, 1–4. Element 2 cannot shelter P because $\tan ({\it\theta}_{2,P})>\tan ({\it\theta}_{2})$ , where $\tan ({\it\theta}_{2})$ is the wake expansion rate of roughness 2. Element 1 cannot shelter P because $x_{P}x_{S1}>h_{1}/\text{tan}({\it\theta}_{1})$ , where $x_{P}$ and $x_{S1}$ are the streamwise coordinates of point P and the ‘sender’ of element 1. Both elements 4 and 3 can shelter P. Their sheltering heights at P can be calculated: $S_{4,P}=\min [h_{4}(x_{P}x_{S4})\tan ({\it\theta}_{4}),h_{5}]$ , $S_{3,P}=\min [h_{3}(x_{P}x_{S3})\tan ({\it\theta}_{3}),h_{5}]$ , where $S_{i,P}$ is the sheltering from $i$ to P, $h_{i}$ is the height of element $i$ and $i$ is 3, 4 or 5. The sheltered height of P is $h_{P}=\max [h_{4,P},h_{3,P}]$ .
Appendix B. Including surface friction drag for ${\it\lambda}_{f}\rightarrow 0$
A conceptual difficulty occurs in the model for $d$ as ${\it\lambda}_{f}\rightarrow 0$ , namely that $d$ does not tend to 0. This difficulty arises because we assumed the bottom surface to be smooth, and we have made the assumption of a fully rough flow regime in which friction drag is entirely neglected. As ${\it\lambda}_{f}\rightarrow 0$ viscous friction drag on the smooth portions of the surface must become comparatively relevant. Therefore, we must consider how the drag is partitioned between roughness elements and the underlying surface (Raupach Reference Raupach1992).
We describe the contribution of the bottom surface and the top of the roughness elements (i.e. the planform surfaces) to the overall momentum loss as coming from either ‘unresolved roughness’ form drag or viscous drag on the bottom surface. We assume that the hydrodynamic roughness height of that unresolved roughness, $z_{o}^{\prime }$ , is known a priori, or in the case of viscous drag it can be determined iteratively as $z_{o}^{\prime }={\it\nu}/u_{{\it\tau}}\exp ({\it\kappa}B_{0})$ (where $B_{0}=5$ is the usual offset of the smoothwall log law and $u_{{\it\tau}}$ is part of the solution). The force on the overall planform surface is modelled as $F_{s}=C_{s}U_{h}^{2}A_{T}$ , with $C_{s}=[{\it\kappa}/\ln (h/z_{0}^{\prime })]^{2}$ . Moreover, the force on the frontal surface is modelled according to $F_{R}=C_{HD}U_{h}^{2}A_{f}/2=C_{R}U_{h}^{2}{\it\lambda}_{f}A_{T}$ , with $C_{R}=C_{DH}/2$ , valid for low coverage fraction according to the derivation presented in § 3.2. Thus, the ratio of forces is
where ${\it\beta}=C_{R}/C_{s}=2[{\it\kappa}/\!\ln (h/z_{0}^{\prime })]^{2}/C_{DH}$ , with $C_{DH}=1.4$ as before. This simple estimate of the force or stress partition derivation is valid for low packing densities. In fact, a similar estimate to (B 1) can be justified also for higher packing ratios since then sheltering reduces both the numerator and the denominator in a similar fashion (Raupach Reference Raupach1992). The momentum balance (3.2) is then augmented by considering $u_{{\it\tau}}^{2}A_{T}=F_{R}(F_{s}/F_{R}+1)=F_{R}(1+{\it\beta}{\it\lambda}_{f})/({\it\beta}{\it\lambda}_{f})$ , i.e.
Equation (3.6) can also be corrected as
Equations (3.4) and (3.5) remain unaffected provided that the corrected values of $z_{0}$ and $d$ are used.
For illustration purposes, consider for example that the underlying surface includes an unresolved roughness with element height $h^{\prime }=0.02h$ . The typical rule of thumb is $z_{o}^{\prime }\sim 0.06h^{\prime }$ (Grimmond & Oke Reference Grimmond and Oke1999), leading to $C_{s}=0.00357$ . With $C_{R}=1.4$ the ratio is, in this case, ${\it\beta}=196$ . The predictions of the model with drag partition correction in this case are plotted in figures 15(b) and 16(b). It is observed that with this correction, $\lim _{{\it\lambda}_{f}\rightarrow 0}d=0$ , $\lim _{{\it\lambda}_{f}\rightarrow 0}z_{o}=z_{o}^{\prime }$ , while this correction is negligible for packing densities ${\it\lambda}_{f}>0.1$ where the fully rough condition is more closely reproduced.