Hostname: page-component-77c89778f8-5wvtr Total loading time: 0 Render date: 2024-07-19T20:06:54.207Z Has data issue: false hasContentIssue false

Exponential roughness layer and analytical model for turbulent boundary layer flow over rectangular-prism roughness elements

Published online by Cambridge University Press:  18 January 2016

Xiang I. A. Yang
Affiliation:
Department of Mechanical Engineering, The Johns Hopkins University, Baltimore, MD 21218, USA
Jasim Sadique
Affiliation:
Department of Mechanical Engineering, The Johns Hopkins University, Baltimore, MD 21218, USA
Rajat Mittal
Affiliation:
Department of Mechanical Engineering, The Johns Hopkins University, Baltimore, MD 21218, USA
Charles Meneveau*
Affiliation:
Department of Mechanical Engineering, The Johns Hopkins University, Baltimore, MD 21218, USA
*
Email address for correspondence: meneveau@jhu.edu

Abstract

We conduct a series of large-eddy simulations (LES) to examine the mean flow behaviour within the roughness layer of turbulent boundary layer flow over rough surfaces. We consider several configurations consisting of arrays of rectangular-prism roughness elements with various spacings, aspect ratios and height distributions. The results provide clear evidence for exponential behaviour of the mean flow with respect to the wall normal distance. Such behaviour has been proposed before (see, e.g., Cionco, 1966 Tech. Rep. DTIC document), and is represented as $U(z)/U_{h}=\exp [a(z/h-1)]$, where $U(z)$ is the spatially/temporally averaged fluid velocity, $z$ is the wall normal distance, $h$ represents the height of the roughness elements and $U_{h}$ is the velocity at $z=h$. The attenuation factor $a$ depends on the density of the roughness element distribution and details of the roughness distribution on the wall. Once established, the generic velocity profile shape is used to formulate a fully analytical model for the effective drag exerted by turbulent flow on a surface covered with arrays of rectangular-prism roughness elements. The approach is based on the von Karman–Pohlhausen integral method, in which a shape function is assumed for the mean velocity profile and its parameters are determined based on momentum conservation and other fundamental constraints. In order to determine the attenuation parameter $a$, wake interactions among surface roughness elements are accounted for by using the concept of flow sheltering. The model transitions smoothly between ‘$k$’ and ‘$d$’ type roughness conditions depending on the surface coverage density and the detailed geometry of roughness elements. Comparisons between model predictions and experimental/numerical data from the existing literature as well as LES data from this study are presented. It is shown that the analytical model provides good predictions of mean velocity and drag forces for the cases considered, thus raising the hope that analytical roughness modelling based on surface geometry is possible, at least for cases when the location of flow separation over surface elements can be easily predicted, as in the case of wall-attached rectangular-prism roughness elements.

Type
Papers
Copyright
© 2016 Cambridge University Press 

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 sand-grain 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 sand-grain 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 rectangular-prism 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 Gonzalez-Flesca2003) 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 case-specific fashion.

In order to develop physics-based 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 mixing-length eddy-viscosity 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 rough-wall 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; Millward-Hopkins et al. Reference Millward-Hopkins, 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 Reynolds-averaged streamwise momentum equation in which the vertical (direction $z$ ) gradient of the momentum flux (modelled using a constant-mixing-length eddy viscosity) balances the distributed drag from roughness elements:

(1.1) $$\begin{eqnarray}\frac{\text{d}}{\text{d}z}\left[l_{m}^{2}\left|\frac{\text{d}U(z)}{\text{d}z}\right|\frac{\text{d}U(z)}{\text{d}z}\right]=C_{d}U(z)^{2}\frac{\text{d}A_{f}}{\text{d}V}.\end{eqnarray}$$

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):

(1.2) $$\begin{eqnarray}U/U_{h}=\exp [a(z/h-1)],\end{eqnarray}$$

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 large-eddy simulations (LES) of turbulent boundary layers over arrays of wall-attached rectangular-prism 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 non-uniform 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 high-Reynolds-number 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 (non-equilibrium) 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 geometry-based 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 non-uniform 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 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 ac) $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 57 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. An analytical rough-wall 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:

(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 57 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. 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$ .

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 25. 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 25.

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

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 mixing-length 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 non-uniform 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. non-frontal) 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 rectangular-prism shape, since the sheltering region within the wake of such objects can be clearly related to the frontal cross-section 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, semi-hemispheres, 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 non-equilibrium 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 N00014-12-1-0582, 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}(y-w)]$ , where $h_{B_{1}A}=\max [0,h-l_{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}(y-w_{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)/2-l_{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}(y-w+w_{2})-\mathscr{H}(y-w)]$ . 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=(w-2w_{2})h_{B_{1}A}+2w_{2}h_{B_{2}A}$ .

Figure 32. (a) A sketch of the roughness interaction in fully staggered cube arrays. Cube $A$ can be sheltered by $B_{1}$ , $B_{2}$ and $B_{3}$ . (b) A sketch of the frontal area of $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 lower-rising and higher-rising elements separately. For the higher-rising elements $A_{1}$ , these could be sheltered by the lower-rising roughness upstream $B_{2}$ and by the higher-rising 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 high-rising roughness respectively. The height of the equivalent sheltered layer $h_{s}$ is then simply $h_{s}=(h_{s,A}+h_{B})/2$ .

Figure 33. (a) Sketch of the interaction among roughness elements for roughness of bimodal height distribution. Here, $A_{i}$ , $i=1,2,3$ , are higher-rise roughnesses and $B_{i}$ , $i=1,2,3$ , are lower-rise roughnesses. (b) Sketch of the wake interaction among roughness elements. In the case shown, the roughness of height $h_{l}$ is completely sheltered from the wake behind the roughness of height $h_{h}$ and no sheltering occurs among the higher-roughness elements. However, parts of the higher-roughness elements are sheltered by the lower-roughness elements. Because the numbers of roughnesses of height $h_{h}$ and roughnesses of height $h_{l}$ are the same, the equivalent sheltered layer height is given by $h_{s}=(h_{s,A}+h_{B})/2$ .

In general, the equivalent sheltered layer height is

(A 1) $$\begin{eqnarray}h_{s}=\mathscr{H}\left(\frac{h_{s,A}+h_{s,B}}{2}-h_{B}\right)(h_{s,A}+h_{s,B}-h_{B})+\left[1-\mathscr{H}\left(\frac{h_{s,A}+h_{s,B}}{2}-h_{B}\right)\right]\frac{h_{s,A}+h_{s,B}}{2}.\end{eqnarray}$$

To determine $h_{s,B}$ and $h_{s,A}$ , the wake interaction among the roughness elements needs to be considered. In general, the higher-roughness–higher-roughness, higher-roughness–lower-roughness and lower-roughness–lower-roughness interactions need to be considered, but since lower-roughness elements are separated by higher-roughness ones, there is no need to consider lower-roughness–lower-roughness 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).

Figure 34. A sketch of the rough wall with rectangular roughness elements. Each element is indicated by a rectangle and given a number (from 1 to 6). The ‘senders’ are highlighted with thick lines and the ‘receivers’ with thick dotted lines. The point P of the ‘receiver’ of element 5 is under consideration. The domain to search for roughness elements that could shelter P is enclosed by dashed lines. It is $3h_{max}U_{h}/u_{{\it\tau}}$ upstream and on both sides. Here, $h_{max}$ is the height of the highest roughness element. The sheltering of the ‘senders’ within the search domain is indicated by thin solid lines.

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 smooth-wall 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

(B 1) $$\begin{eqnarray}\frac{F_{s}}{F_{R}}=\frac{1}{{\it\lambda}_{f}{\it\beta}},\end{eqnarray}$$

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.

(B 2) $$\begin{eqnarray}u_{{\it\tau}}^{2}A_{T}=\frac{1+{\it\beta}{\it\lambda}_{f}}{{\it\beta}{\it\lambda}_{f}}\int _{A_{f}}C_{d}U^{2}\,\text{d}A.\end{eqnarray}$$

Equation (3.6) can also be corrected as

(B 3) $$\begin{eqnarray}d=\frac{{\it\beta}{\it\lambda}_{f}}{1+{\it\beta}{\it\lambda}_{f}}\frac{\displaystyle \int _{A_{f}}C_{d}U^{2}z\,\text{d}A}{\displaystyle \int _{A_{f}}C_{d}U^{2}\,\text{d}A}.\end{eqnarray}$$

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.

References

Akins, R. E., Peterka, J. A. & Cermak, J. E. 1977 Mean force and moment coefficients for buildings in turbulent boundary layers. J. Wind Engng Ind. Aerodyn. 2 (3), 195209.Google Scholar
Arnfield, A. J. 2003 Two decades of urban climate research: a review of turbulence, exchanges of energy and water, and the urban heat island. Intl J. Climatol. 23 (1), 126.Google Scholar
Barlow, J. F. & Coceal, O.2009 A review of urban roughness sublayer turbulence. Met Office Research and Development – Technical Report 1 (527).Google Scholar
Bhardwaj, R. & Mittal, R. 2012 Benchmarking a coupled immersed-boundary–finite-element solver for large-scale flow-induced deformation. AIAA J. 50 (7), 16381642.Google Scholar
Castro, I. P. 2007 Rough-wall boundary layers: mean flow universality. J. Fluid Mech. 585, 469485.Google Scholar
Cheng, H. & Castro, I. P. 2002 Near wall flow over urban-like roughness. Boundary-Layer Meteorol. 104 (2), 229259.Google Scholar
Cheng, H., Hayden, P., Robins, A. G. & Castro, I. P. 2007 Flow over cube arrays of different packing densities. J. Wind Engng Indust. Aerodyn. 95 (8), 715740.Google Scholar
Cionco, R. M.1966 A mathematical model for air flow in a vegetative canopy. Tech. Rep. DTIC Document.Google Scholar
Coceal, O. & Belcher, S. E. 2004 A canopy model of mean winds through urban areas. Q. J. R. Meteorol. Soc. 130 (599), 13491372.CrossRefGoogle Scholar
Colebrook, C. F. 1939 Turbulent flow in pipes, with particular reference to the transition region between the smooth and rough pipe laws. J. ICE 11 (4), 133156.Google Scholar
Coleman, S. E., Nikora, V. I., McLean, S. R. & Schlicke, E. 2007 Spatially averaged turbulent flow over square ribs. J. Engng Mech. 133 (2), 194204.Google Scholar
Coles, D. 1956 The law of the wake in the turbulent boundary layer. J. Fluid Mech. 1 (02), 191226.Google Scholar
Counehan, J 1971 Wind tunnel determination of the roughness length as a function of the fetch and the roughness density of three-dimensional roughness elements. Atmos. Environ. 5 (8), 637642.Google Scholar
Cui, J., Patel, V. C., Lin, C.-l. & Patel, V. C. 2000 Large-Eddy Simulation of Turbulent Flow over Rough Surfaces. Iowa Institute of Hydraulic Research, College of Engineering, University of Iowa.Google Scholar
Curley, A. & Uddin, M.2015 Direct numerical simulation of turbulent flow around a surface mounted cube. AIAA Paper 3431.CrossRefGoogle Scholar
Di Sabatino, S., Solazzo, E., Paradisi, P. & Britter, R. 2008 A simple model for spatially-averaged wind profiles within and above an urban canopy. Boundary-Layer Meteorol. 127 (1), 131151.Google Scholar
Flack, K. A. & Schultz, M. P. 2010 Review of hydraulic roughness scales in the fully rough regime. J. Fluids Engng 132 (4), 041203.Google Scholar
Flack, K. A., Schultz, M. P. & Shapiro, T. A. 2005 Experimental support for Townsend’s Reynolds number similarity hypothesis on rough walls. Phys. Fluids 17 (3), 035102.Google Scholar
Grimmond, C. & Oke, T. R. 1999 Aerodynamic properties of urban areas derived from analysis of surface form. J. Appl. Meteorol. 38 (9), 12621292.Google Scholar
Hagishima, A., Tanimoto, J., Nagayama, K. & Meno, S. 2009 Aerodynamic parameters of regular arrays of rectangular blocks with various geometries. Boundary-Layer Meteorol. 132 (2), 315337.Google Scholar
Hall, D., Macdonald, R. & Walker, S. 1996 Measurements of Dispersion within Simulated Urban Arrays: A Small Scale Wind Tunnel Study. Building Research Establishment.Google Scholar
Harman, I. N. & Finnigan, J. J. 2007 A simple unified theory for flow in the canopy and roughness sublayer. Boundary-Layer Meteorol. 123 (2), 339363.Google Scholar
Hussain, M. & Lee, B. E. 1980 A wind tunnel study of the mean pressure forces acting on large groups of low-rise buildings. J. Wind Engng Indust. Aerodyn. 6 (3), 207225.Google Scholar
Ikeda, T. & Durbin, P. A. 2007 Direct simulations of a rough-wall channel flow. J. Fluid Mech. 571, 235263.Google Scholar
Jackson, P. S. 1981 On the displacement height in the logarithmic velocity profile. J. Fluid Mech. 111, 1525.CrossRefGoogle Scholar
Jiang, D., Jiang, W., Liu, H. & Sun, J. 2008 Systematic influence of different building spacing, height and layout on mean wind and turbulent characteristics within and over urban building arrays. Wind Struct. 11 (4), 275289.Google Scholar
Jiménez, J. 2004 Turbulent flows over rough walls. Annu. Rev. Fluid Mech. 36, 173196.Google Scholar
Kanda, M. 2006 Progress in the scale modeling of urban climate: review. Theor. Appl. Climatol. 84 (1–3), 2333.Google Scholar
Kanda, M., Moriwaki, R. & Kasamatsu, F. 2004 Large-eddy simulation of turbulent organized structures within and above explicitly resolved cube arrays. Boundary-Layer Meteorol. 112 (2), 343368.Google Scholar
Kim, B.-G., Lee, C., Joo, S., Ryu, K.-C., Kim, S., You, D. & Shim, W.-S. 2011 Estimation of roughness parameters within sparse urban-like obstacle arrays. Boundary-Layer Meteorol. 139 (3), 457485.Google Scholar
Leonardi, S. & Castro, I. P. 2010 Channel flow over large cube roughness: a direct numerical simulation study. J. Fluid Mech. 651, 519539.CrossRefGoogle Scholar
Leonardi, S., Orlandi, P. & Antonia, R. A. 2007 Properties of $d$ - and $k$ -type roughness in a turbulent channel flow. Phys. Fluids 19 (12), 125101.Google Scholar
Leonardi, S., Orlandi, P., Smalley, R. J., Djenidi, L. & Antonia, R. A. 2003 Direct numerical simulations of turbulent channel flow with transverse square bars on one wall. J. Fluid Mech. 491, 229238.CrossRefGoogle Scholar
Lettau, H. 1969 Note on aerodynamic roughness-parameter estimation on the basis of roughness-element description. J. Appl. Meteorol. 8 (5), 828832.Google Scholar
Macdonald, R., Griffiths, R. & Hall, D. 1998 An improved method for the estimation of surface roughness of obstacle arrays. Atmos. Environ. 32 (11), 18571864.CrossRefGoogle Scholar
Macdonald, R. W. 2000 Modelling the mean velocity profile in the urban canopy layer. Boundary-Layer Meteorol. 97 (1), 2545.Google Scholar
Millward-Hopkins, J. T., Tomlin, A. S., Ma, L., Ingham, D. & Pourkashanian, M. 2011 Estimating aerodynamic parameters of urban-like surfaces with heterogeneous building heights. Boundary-Layer Meteorol. 141 (3), 443465.Google Scholar
Mittal, R., Dong, H., Bozkurttas, M., Najjar, F. M., Vargas, A. & von Loebbecke, A. 2008 A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries. J. Comput. Phys. 227 (10), 48254852.Google Scholar
Moody, L. F. 1944 Friction factors for pipe flow. Trans. ASME 66 (8), 671684.Google Scholar
Napoli, E., Armenio, V. & De Marchis, M. 2008 The effect of the slope of irregularly distributed roughness elements on turbulent wall-bounded flows. J. Fluid Mech. 613, 385394.Google Scholar
Nikuradse, J. 1950 Laws of Flow in Rough Pipes. National Advisory Committee for Aeronautics Washington.Google Scholar
Placidi, M. & Ganapathisubramani, B. 2015 Effects of frontal and plan solidities on aerodynamic parameters and the roughness sublayer in turbulent boundary layers. J. Fluid Mech. 782, 541566.Google Scholar
Pohlhausen, K. 1921 On the approximate integration of the differential equations of laminar shear layers. Z. Angew. Math. Mech. 1 (4), 252290.Google Scholar
Raupach, M. 1992 Drag and drag partition on rough surfaces. Boundary-Layer Meteorol. 60 (4), 375395.Google Scholar
Raupach, M., Finnigan, J. J. & Brunei, Y. 1996 Coherent eddies and turbulence in vegetation canopies: the mixing-layer analogy. Boundary-Layer Meteorol. 78 (3–4), 351382.Google Scholar
Raupach, M. R., Antonia, R. A. & Rajagopalan, S. 1991 Rough-wall turbulent boundary layers. Appl. Mech. Rev. 44 (1), 125.Google Scholar
Santiago, J. L., Coceal, O., Martilli, A. & Belcher, S. E. 2008 Variation of the sectional drag coefficient of a group of buildings with packing density. Boundary-Layer Meteorol. 128 (3), 445457.Google Scholar
Schultz, M., Schatzmann, M. & Leitl, B. 2005 Effect of roughness inhomogeneities on the development of the urban boundary layer. Intl J. Environ. Pollut. 25 (1–4), 105117.Google Scholar
Schultz, M. P. & Flack, K. A. 2009 Turbulent boundary layers on a systematically varied rough wall. Phys. Fluids 21 (1), 015104.Google Scholar
Shao, Y. & Yang, Y. 2005 A scheme for drag partition over rough surfaces. Atmos. Environ. 39 (38), 73517361.Google Scholar
Shao, Y. & Yang, Y. 2008 A theory for drag partition over rough surfaces. J. Geophys. Res. 113, F02S05.Google Scholar
Simpson, R. L. 1973 A generalized correlation of roughness density effects on the turbulent boundary layer. AIAA J. 11 (2), 242244.Google Scholar
Vanderwel, C. & Ganapathisubramani, B. 2015 Effects of spanwise spacing on large-scale secondary flows in rough-wall turbulent boundary layers. J. Fluid Mech. 774, R2.Google Scholar
Vardoulakis, S., Fisher, B. E., Pericleous, K. & Gonzalez-Flesca, N. 2003 Modelling air quality in street canyons: a review. Atmos. Environ. 37 (2), 155182.Google Scholar
Vedula, V., Fortini, S., Seo, J. H., Querzoli, G. & Mittal, R. 2014 Computational modeling and validation of intraventricular flow in a simple model of the left ventricle. Theor. Comput. Fluid Dyn. 28 (6), 589604.Google Scholar
Vreman, A. W. 2004 An eddy-viscosity subgrid-scale model for turbulent shear flow: algebraic theory and applications. Phys. Fluids 16 (10), 36703681.Google Scholar
Yang, X. I. A & Meneveau, C. 2016 Recycling inflow method for simulations of spatially evolving turbulent boundary layers over rough surfaces. J. Turbul. 17 (1), 7593.Google Scholar
Yang, X. I. A., Sadique, J., Mittal, R. & Meneveau, C. 2015 Integral wall model for large eddy simulations of wall-bounded turbulent flows. Phys. Fluids 27 (2), 025112.Google Scholar
You, D. & Moin, P. 2007 A dynamic global-coefficient subgrid-scale eddy-viscosity model for large-eddy simulation in complex geometries. Phys. Fluids 19 (6), 065110.Google Scholar
Zaki, S. A., Hagishima, A., Tanimoto, J. & Ikegaya, N. 2011 Aerodynamic parameters of urban building arrays with random geometries. Boundary-Layer Meteorol. 138 (1), 99120.Google Scholar
Zheng, L., Hedrick, T. L. & Mittal, R. 2013 A multi-fidelity modelling approach for evaluation and optimization of wing stroke aerodynamics in flapping flight. J. Fluid Mech. 721, 118154.CrossRefGoogle Scholar
Figure 0

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 1

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

Figure 2

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$.

Figure 3

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 4

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 5

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

Figure 6

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)).

Figure 7

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$.

Figure 8

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.

Figure 9

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

Figure 10

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.

Figure 11

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.

Figure 12

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.

Figure 13

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.

Figure 14

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.

Figure 15

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. (1996), Cheng et al. (2007), Hagishima et al. (2009) and Leonardi & Castro (2010), 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 (1992), for the aligned arrangement (see appendix B).

Figure 16

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

Figure 17

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

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 (2000) and Coceal & Belcher (2004) do not depend on the relative arrangement of roughness elements, only on ${\it\lambda}_{f}$.

Figure 19

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

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

Figure 21

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.2000; Leonardi et al.2003; Coleman et al.2007; Ikeda & Durbin 2007), 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.

Figure 22

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

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$.

Figure 24

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

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.

Figure 26

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

Figure 27

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

Figure 28

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

Figure 29

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

Figure 30

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.

Figure 31

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.

Figure 32

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

Figure 33

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

Figure 34

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

Figure 35

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 36

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 37

Figure 32. (a) A sketch of the roughness interaction in fully staggered cube arrays. Cube $A$ can be sheltered by $B_{1}$, $B_{2}$ and $B_{3}$. (b) A sketch of the frontal area of $A$.

Figure 38

Figure 33. (a) Sketch of the interaction among roughness elements for roughness of bimodal height distribution. Here, $A_{i}$, $i=1,2,3$, are higher-rise roughnesses and $B_{i}$, $i=1,2,3$, are lower-rise roughnesses. (b) Sketch of the wake interaction among roughness elements. In the case shown, the roughness of height $h_{l}$ is completely sheltered from the wake behind the roughness of height $h_{h}$ and no sheltering occurs among the higher-roughness elements. However, parts of the higher-roughness elements are sheltered by the lower-roughness elements. Because the numbers of roughnesses of height $h_{h}$ and roughnesses of height $h_{l}$ are the same, the equivalent sheltered layer height is given by $h_{s}=(h_{s,A}+h_{B})/2$.

Figure 39

Figure 34. A sketch of the rough wall with rectangular roughness elements. Each element is indicated by a rectangle and given a number (from 1 to 6). The ‘senders’ are highlighted with thick lines and the ‘receivers’ with thick dotted lines. The point P of the ‘receiver’ of element 5 is under consideration. The domain to search for roughness elements that could shelter P is enclosed by dashed lines. It is $3h_{max}U_{h}/u_{{\it\tau}}$ upstream and on both sides. Here, $h_{max}$ is the height of the highest roughness element. The sheltering of the ‘senders’ within the search domain is indicated by thin solid lines.