3.1 The velocity and vorticity fields

Figure 3.

Wall-normal shear stress as function of the distance from the virtual wall (broken line). Black line: run D50; red line: D120. The symbol $\times$  indicates the value of the bottom shear stress extrapolated down to a wall-normal distance $y=y_{0}$ using the slope of the linear shear stress profile far from the wall.

Figure 4.

Profiles of (a) $\langle \overline{u}\rangle ^{+}$ and of (b) $(y^{+}-y_{0}^{+})\,\text{d}\langle \overline{u}\rangle ^{+}/\text{d}y^{+}$ as a function of the distance from the virtual wall. Lines —▫— and —○— indicate cases D50 and D120, respectively, while solid lines indicate the simulations performed over a smooth wall at $Re_{b}=2900$ and $Re_{b}=6864$ , respectively. The broken lines in (b) indicate the value of the logarithmic constant $1/\unicode[STIX]{x1D705}=2.44$ (central line) $\pm 0.12$ (upper and lower lines). Symbols ▵ (blue) and $\times$ (blue) indicate the values measured by Amir et al. (Reference Amir, Nikora and Stewart2014) in their experiments number 1 and 4, respectively.

Since in the transitionally rough regime viscous effects are still relevant over the roughness, and are presumably relevant also in the fully rough regime at least along the crevices between the roughness elements, viscous scales will be used as reference scales. Hence, the friction velocity $u_{\unicode[STIX]{x1D70F}}=\sqrt{\unicode[STIX]{x1D70F}_{w}/\unicode[STIX]{x1D71A}}$ was estimated, where the value of the wall shear stress $\unicode[STIX]{x1D70F}_{w}$ was extrapolated down to the wall-normal distance $y=y_{0}$ using the linear profile of the wall-normal total shear stress $\unicode[STIX]{x1D70F}_{tot}=\unicode[STIX]{x1D71A}\unicode[STIX]{x1D708}(\langle \overline{u}\rangle /\text{d}y)-\unicode[STIX]{x1D71A}\langle \overline{u^{\prime \prime \prime }v^{\prime \prime \prime }}\rangle$ far from the rough wall, as shown in figure 3 (Chan-Braun et al. Reference Chan-Braun, García-Villalba and Uhlmann2011). Figure 4(a) shows the velocity profiles which were computed for the runs D50 and D120 and for the respective simulations performed at the same bulk Reynolds numbers in the absence of the roughness elements. The value of $y_{0}$ for the simulation D50 was set equal to $0.8D$ according to the indications of Chan-Braun et al. (Reference Chan-Braun, García-Villalba and Uhlmann2011). Indeed, it was found that the profiles of $(y^{+}-y_{0}^{+})\text{d}\langle \overline{u}\rangle ^{+}/\text{d}y^{+}$ (which is equal to an inverse von Kármán ‘constant’) for the run D50 approached, in the logarithmic region, the curve obtained from the respective simulation over a smooth wall at $Re_{bH}\approx 2900$ (see figure 4 b). A similar agreement between the rough and smooth wall profiles was obtained for the run D120 by positioning the virtual wall $y_{0}$ at the distance $0.85D$ . This choice allows us to compare velocity profiles obtained for the smooth- and rough-wall cases in the logarithmic region and to estimate the roughness function. It is worthwhile mentioning that other methodologies adopted to estimate the position of $y_{0}$ , such as the location of the drag force centroid (Jackson Reference Jackson1981), provided approximately to the same value (see also appendix A1 of Chan-Braun et al. Reference Chan-Braun, García-Villalba and Uhlmann2011). Since, in the logarithmic region, the curves of figure 4(b) should be independent of $y^{+}$ , it was also possible to estimate the value for the von Kármán constant, $\unicode[STIX]{x1D705}$ , as the inverse of the value attained at the local minimum of the curves, as suggested (among others) by Balachandar & Ramachandran (Reference Balachandar and Ramachandran1999), Tachie et al. (Reference Tachie, Bergstrom and Balachandar2000) and George (Reference George2007). Hoyas & Jiménez (Reference Hoyas and Jiménez2006) showed that the fact of not observing a wide region where $1/\unicode[STIX]{x1D705}$ is constant is an effect of the higher-order terms usually included in the wake component of the profile. Nonetheless, these authors observed the existence of a substantial logarithmic region. Therefore, the von Kármán constant was estimated equal to 0.388 and 0.381 for the simulations D50 and D120, respectively. Slight deviations of the value of $\unicode[STIX]{x1D705}$ from the value 0.41, estimated by Coles (Reference Coles1968) for boundary layers, were often reported in the literature (e.g. Österlund et al. Reference Österlund, Johansson, Nagib and Hites2000; Marusic et al. Reference Marusic, Monty, Hultmark and Smits2013) and, in particular for the present channel-flow configuration, can be related to pressure gradient effects (George Reference George2007). However, it should be noted that the present values of $\unicode[STIX]{x1D705}$ are in the range estimated by Marusic et al. (Reference Marusic, Monty, Hultmark and Smits2013) ( $\unicode[STIX]{x1D705}=0.39\pm 0.02$ ), while a local minimum of $(y^{+}-y_{0}^{+})\text{d}\langle \overline{u}\rangle ^{+}/\text{d}y^{+}$ is attained approximately at $y^{+}-y_{0}^{+}=50$ consistent with other numerical (Hoyas & Jiménez Reference Hoyas and Jiménez2006) and experimental (e.g. Tachie et al. Reference Tachie, Bergstrom and Balachandar2000; Mckeon et al. Reference Mckeon, Li, Jiang, Morrison and Smits2004) results. Some considerations can be formulated for the simulation D50 on the basis of the velocity profiles of figure 4(a). Chan-Braun et al. (Reference Chan-Braun, García-Villalba and Uhlmann2011) found that the flow regime for the simulation D50 was transitionally rough. This also appears from the velocity profile (black square line in figure 4 a) which shows the presence of the buffer layer above the crest of the spheres. The thickness of the viscous sublayer $\unicode[STIX]{x1D6FF}_{sub}^{+}$ was estimated from the intersection of the logarithmic law and the linear law of the wall (cf. figure 4 a), which yields $\unicode[STIX]{x1D6FF}_{sub}^{+}=11.5$ in the smooth-wall case at $Re_{bH}=2900$ . In case D50 we formally apply the same procedure, although a linear law is not observed; this yields $\unicode[STIX]{x1D6FF}_{sub}^{+}=5.0$ in this case. Then, the constant $C_{I}^{+}(D^{+}\rightarrow 0)$ shown in (1.1) can be determined by extrapolating the logarithmic profile down to $y^{+}-y_{0}^{+}=1$ and it was equal to 5.5 for the smooth-wall case at $Re_{bH}\approx 2900$ . According to Ligrani & Moffat (Reference Ligrani and Moffat1986), in the transitionally rough regime, the roughness function, $\unicode[STIX]{x1D6E5}\langle \overline{u}\rangle ^{+}$ , can be estimated as follows:

which was equal to 4.4, where $C_{I}^{+}(D^{+}\rightarrow 0)$ indicates the value of $C_{I}^{+}$ for the smooth-wall case and the value of $\unicode[STIX]{x1D705}$ is that associated with the smooth-wall simulation ( $\unicode[STIX]{x1D705}\approx 0.41$ ). The value of $C_{II}^{+}$ can be also estimated on the basis of the expression (1.2) evaluated at $y=D+y_{0}$ :

If the dependence on the relative submergence is neglected for a moment, then $C_{II}^{+}$ can be interpreted as equal to $B^{+}-1/\unicode[STIX]{x1D705}\ln k_{s}/D$ . The value of $B^{+}$ can be obtained for instance from the diagram of figure 1 of Ligrani & Moffat (Reference Ligrani and Moffat1986), leaving $k_{s}^{+}$ approximately equal to 30 for the run D50. However, $C_{II}^{+}$ is affected not only by effects associated with the arrangement of the roughness elements, but also by those related to the relative submergence $H/D$ . Indeed, on the basis of the range indicated in figure 8 of Chan-Braun et al. (Reference Chan-Braun, García-Villalba and Uhlmann2011) for the case D50, the actual value of $k_{s}^{+}$ could be presumably larger than the value presently estimated. An expression similar to (3.7) can be written also for the fully rough regime:

where the values of $\unicode[STIX]{x1D705}$ and $y_{0}$ are those attained over a smooth wall, i.e. 0.41 and 0, respectively. Although the values of $\unicode[STIX]{x1D6E5}\langle \overline{u}\rangle ^{+}$ in expressions (3.7) and (3.9) are evaluated at different distances from $y_{0}$ (i.e. $\unicode[STIX]{x1D6FF}_{sub}$ and $D$ , respectively), their dependence on the distance in the range $|D-\unicode[STIX]{x1D6FF}_{sub}|$ was small (i.e. ${\sim}0.8$ for $D^{+}\sim 100$ ). Note that this dependence of $\unicode[STIX]{x1D6E5}\langle \overline{u}\rangle ^{+}$ on the distance from the wall is due to deviations of $\unicode[STIX]{x1D705}$ which are small. Similarly to boundary layers for which $k_{s}$ can be proportional to $k$ ( $k$ -roughness) or to the boundary layer thickness ( $d$ -roughness), for shallow open channels $k_{s}$ possibly depends on both $H$ and $k$ (as well as on the roughness geometrical features) and, therefore, the present estimate of $C_{II}^{+}$ should be associated only with the particular configuration of the present simulations. However, by noting that, in the fully rough regime, $\unicode[STIX]{x1D6E5}\langle \overline{u}\rangle ^{+}$ monotonically increases with $D^{+}$ for a given configuration of the roughness (e.g. see Jiménez Reference Jiménez2004), combining (3.7) and (3.9), the following inequality can be found which poses an upper limit to the value of $C_{II}$ :

Since the sum of the first two terms on the right-hand side of (3.10) is minimum for $\unicode[STIX]{x1D6FF}_{sub}^{+}=1/\unicode[STIX]{x1D705}$ (which is an admissible value of $\unicode[STIX]{x1D6FF}_{sub}^{+}$ in the transitionally rough regime) and is equal to 0.26 and to 0.08 for $\unicode[STIX]{x1D705}$ equal to 0.41 and to 0.38, respectively, it is possible to approximate the inequality (3.10) with the following one:

where $\unicode[STIX]{x1D705}=0.41$ . The inequality (3.11), is a bound for the value of $C_{II}^{+}$ independent of the particular flow configuration. This does not imply the independence for $C_{II}^{+}$ , which instead is a function of both the roughness geometrical features and the relative submergence. In other words, for a given value of the roughness Reynolds number $k^{+}$ , the configuration (arrangement and shape of roughness elements, relative submergence,  $\ldots$ ) which minimises the flow resistance is the configuration for which $C_{II}^{+}$ tends to $(1/\unicode[STIX]{x1D705})\ln (k^{+})$ . For the simulation D120, $\unicode[STIX]{x1D6E5}\langle \overline{u}\rangle ^{+}=6.3$ (provided that the $C_{I}^{+}(D^{+}\rightarrow 0)$ is equal to $5.2$ for the smooth-wall case at $Re_{bH}\approx 6800$ ) and the left and right sides of the inequality (3.11) are equal to 10.5 and 11.6, respectively. The fact that $C_{II}^{+}$ approaches its upper limit indicates that the present flow configuration is fairly conductive, i.e. that the quantity $\unicode[STIX]{x1D71A}{U_{bh}}^{2}/\unicode[STIX]{x1D70F}_{w}$ is relatively large (see table 1). Table 2 shows the values of $C_{II}^{+}$ and $1/\unicode[STIX]{x1D705}\ln (D^{+})$ which were also computed for five experiments of Amir et al. (Reference Amir, Nikora and Stewart2014). These authors investigated open-channel flow over a layer of spheres at rest on a smooth wall in a hexagonal pattern and assumed that $y_{0}=D$ . Although the range of values of $D^{+}$ and their choice of $y_{0}$ are somewhat different, the experimental results of Amir et al. (Reference Amir, Nikora and Stewart2014) suggest that the effect of increasing $D^{+}$ is to decrease the conductivity. In particular, at their smallest value $D^{+}=170$ , which is not too far off the value of our present case D120, the deviation of $C_{II}^{+}$ from its upper limit is similar to that computed for the present run D120. Indeed, table 2 shows that the roughness function is also very similar. The same can be said of their case number $4$ , where $D^{+}=243$ , but where the relative submergence $H/D=5$ roughly matches the present value. The values of the mean velocity measured in the experiments numbers 1 and 4 of Amir et al. (Reference Amir, Nikora and Stewart2014) are also shown in figure 4(a) as a function of $y^{+}-y_{0}^{+}$ , where $y_{0}^{+}$ is assumed equal to $0.85~D^{+}$ . A region where the mean velocity followed a logarithmic profile can be observed also in the latter two experiments. In case D50, the inequality (3.11) is no longer valid mainly because the viscous sublayer is significantly thick in the transitionally rough regime. However, it can be verified that the equality (3.10) is approximately satisfied if also the dependency of $\unicode[STIX]{x1D6E5}\langle \overline{u}\rangle ^{+}$ on the distance from the wall is taken into account. This indicates that the present square arrangement of the spheres almost maximises the flow conductivity in the transitionally rough regime.

Table 2.

Values of the parameters of the integration constant $C_{II}^{+}$ and of the shift of the velocity profile $\unicode[STIX]{x1D6E5}\langle \overline{u}\rangle ^{+}$ (computed at $y^{+}=y_{0}^{+}+D^{+}$ ) for the present simulation D120, for the simulation D50 (Chan-Braun et al. Reference Chan-Braun, García-Villalba and Uhlmann2011), and for six of the experiments carried out by Amir et al. (Reference Amir, Nikora and Stewart2014) (original numbering). Note that Amir et al. (Reference Amir, Nikora and Stewart2014) chose $y_{0}=D$ as opposed to the present choice $y_{0}=0.85D$ .

Ignoring for a moment the dependency on $H/D$ , as already proposed for the transitionally rough regime above, it would be found that $\ln k_{s}/D=(8.5-C_{II}^{+})\unicode[STIX]{x1D705}$ and the value of $k_{s}^{+}$ can be estimated to measure approximately 50. As mentioned above, this value is not a reliable indicator of the flow regime which can be significantly underestimated. In fact, experiment number 1 of Amir et al. (Reference Amir, Nikora and Stewart2014) was characterised by a value of the roughness function relatively small compared with that of Nikuradse (Reference Nikuradse1933), although the flow regime was fully rough. Also using a $k$ -criterion based on the value of the roughness parameter $k_{0}^{+}=D^{+}/30$ (also termed roughness length) to classify the flow regime, as suggested by Jayatilleke (Reference Jayatilleke1966) and Reynolds (Reference Reynolds1974), the flow of simulation D120 falls well into the fully rough range (see the diagram in figure 3.3(a) of Pimenta et al. (Reference Pimenta, Moffat and Kays1975) for a comparison). In fact, it is shown in the following that the flow structure related to the simulation D120 exhibits the typical characteristics of the fully rough regime.

Figure 5.

Panel (a) shows the profiles of the variance of the quantity $\langle \overline{\unicode[STIX]{x1D719}}\rangle _{B}$ reduced by the square product of $\langle \overline{\unicode[STIX]{x1D719}}\rangle$ , which equals the variance of $\unicode[STIX]{x1D719}^{\prime \prime }-\unicode[STIX]{x1D719}^{\prime \prime \prime }$ . Full symbols indicate $\unicode[STIX]{x1D719}\equiv u$ (normalised by $u_{\unicode[STIX]{x1D70F}}$ ) while empty symbols indicate $\unicode[STIX]{x1D719}\equiv p$ (normalised by $\unicode[STIX]{x1D71A}u_{\unicode[STIX]{x1D70F}}^{2}$ ). ▵, ▴: run D50; ○, ●: run D120. Panel (b) shows the profiles of the average velocity field computed for the run D120 at the following $(\widetilde{x},\widetilde{z})$ coordinates that are also indicated in the top view of the small inset: ○: $(0,0)$ , ▫: $(-L_{B}/4,0)$ , ▵: $(-L_{B}/2,0)$ , ▿: $(0,L_{B}/4)$ , ♢: $(0,L_{B}/2)$ .

The upper bound of the roughness sublayer was defined as the distance from the wall at which the magnitude of pressure spatial disturbance $p^{\prime \prime }-p^{\prime \prime \prime }$ , i.e. the fluctuations induced by the roughness geometry, become smaller than 1 % of its maximum value. This methodology to characterise the ‘three-dimensionality’ of the turbulence structure was also adopted by Chan-Braun et al. (Reference Chan-Braun, García-Villalba and Uhlmann2011). Figure 5(a), shows the limit of the roughness sublayer compared with the wall-normal profile of stress defined by (3.4) for the streamwise velocity component and for pressure (i.e. the variance of $p^{\prime \prime }-p^{\prime \prime \prime }$ normalised by $(\unicode[STIX]{x1D71A}u_{\unicode[STIX]{x1D70F}}^{2})^{2}$ ) for the simulations D50 and D120. It can be noted that the decay of form-induced pressure fluctuations with the distance from the crest of the spheres lies on the same line in the semi-logarithmic plot of figure 5(a). Note that this exponential decay of the three-dimensionality of pressure was also observed by Chan-Braun et al. (Reference Chan-Braun, García-Villalba and Uhlmann2011) for their case with $D^{+}=10.7$ . The curves related to the form-induced velocity fluctuations do not exhibit this exponential variation with wall distance, even though the maximum variance of both velocity and pressure fluctuations is attained in proximity of the crest of the spheres for both cases D50 and D120. The residual small variance that can be detected well above the roughness sublayer, was also observed by Florens et al. (Reference Florens, Eiff and Moulin2013) who ascribed it to time convergence error. However, this small residual could be associated with the roughness footprint that Hong et al. (Reference Hong, Katz and Schultz2011) observed on the small-scale turbulence, which persisted over the entire flow domain, although its effect on the Reynolds stress statistics was found negligible. For both the cases D50 and D120, the bound of the roughness sublayer was identified approximately at $y=1.8D$ (i.e. $y-y_{0}=0.95D$ ) which is in line with the experimental observations of Cheng & Castro (Reference Cheng and Castro2002) and of Hong et al. (Reference Hong, Katz and Schultz2011) ( $y=2k$ ). Streamwise velocity fluctuations attain 1 % of their maximum root mean square value at $y\sim 1.45D$ which is in fair agreement with the experimental results of Florens et al. (Reference Florens, Eiff and Moulin2013) who report $y=1.5k$ . The dispersion effect in the vicinity of the roughness elements can be also appreciated by considering velocity profiles obtained at different locations with respect to the centre of the spheres, as shown in figure 5(b). Although the streamwise velocity $\langle \overline{u}\rangle _{B}^{+}$ is remarkably dispersed below the crest of the spheres, it rapidly converges with growing distance above the crests until it coincides with the plane/time-averaged velocity in the upper part of the roughness sublayer.

Figure 6.

Profiles of the terms of the Reynolds stress tensor (a) $\left\langle \overline{u^{\prime \prime \prime }u^{\prime \prime \prime }}\right\rangle ^{+}$ and (b) $\left\langle \overline{u^{\prime \prime }u^{\prime \prime }}\right\rangle ^{+}$ for $Re_{b}=2900$ (black lines) and $Re_{b}=6864$ (red lines) as a function of the distance from the virtual wall. Lines —▫— and —○— indicate runs D50 and D120, respectively, while solid lines indicate the corresponding smooth-wall simulations. The horizontal broken lines indicate the position of the crest of the spheres for the two runs.

The profiles of the mean square of streamwise velocity fluctuations were computed (figure 6 a,b) by considering the definitions of fluctuations given by (3.2) and (3.3), respectively, in order to evaluate the effects of the roughness geometrical features on the dimensionless streamwise normal stress (i.e. the streamwise turbulence intensity). While the profiles for the simulation D50 do not show significant discrepancies for $y^{+}>y_{0}^{+}$ by using either of the two definitions, the profiles obtained for the simulation D120 are remarkably different. This difference is related to the particular arrangement of the spheres chosen here, which, we recall, is the same for the two runs, and which allows the flow to stream along streamwise intra-roughness grooves, causing a steady and predominantly two-dimensional spanwise-periodic pattern (of periodicity $L_{B}$ , see figure 2). This interfacial flow is captured by the sphere-box average, whereas it is filtered out by the plane average. Consequently, the intensity of the groove-induced flow pattern enters the field $\boldsymbol{u}^{\prime \prime }$ , but not $\boldsymbol{u}^{\prime \prime \prime }$ . In fact, the maximum of $\langle \overline{u^{\prime \prime \prime }u^{\prime \prime \prime }}\rangle ^{+}$ is located over the crests of roughness elements for the run D50 while it falls below the crests for the run D120 (figure 6 a) and, in the latter case, the maximum disappears if the fluctuations $u^{\prime \prime }$ about the average flow field are considered (figure 6 b). Larger values of the dispersive stress (of the streamwise velocity) below the crest of the spheres for the case D120 can be also noted in figure 5 with respect to those attained for the case D50. The destruction of the maximum normal stresses is associated with that of the buffer layer and is, in our opinion, the most clear effect of the fully rough regime. In other words, it was observed that the fully rough regime was attained as $D^{+}-y_{0}^{+}$ exceeded ${\sim}$ 12, which is almost the distance attained by the maximum of $\langle \overline{u^{\prime \prime \prime }u^{\prime \prime \prime }}\rangle ^{+}$ from a smooth wall.

Figure 7.

Averaged streamwise velocity $\langle \overline{u}\rangle _{B,x}^{+}$ visualised by shadowed contour lines equispaced by 0.2 and by 2 for values smaller and larger than 1, respectively. Red and cyan lines indicate the values 0 and 1. The spheres are outlined by green lines while the horizontal line [- $+$ -] indicates the upper limit of the roughness sublayer. (a) Run D50; (b) run D120.

Figure 8.

Averaged square of the fluctuations of the averaged streamwise velocity $\langle \overline{u^{\prime \prime }u^{\prime \prime }}\rangle _{B,x}^{+}$ and $\langle \overline{u^{\prime \prime }u^{\prime \prime }}\rangle _{B,z}^{+}$ visualised in front view (a,b) and side view (c,d), respectively, by shadowed contour lines equispaced by 0.67. Cyan and red lines indicate the values 3 and 4, respectively. The spheres are outlined by green lines while the horizontal line [- $+$ -] indicates the upper limit of the roughness sublayer. White rightward arrows indicate the direction of the flow. (a,c) Run D50; (b,d) run D120.

The same picture arises also from the comparison of the average flow field computed for the simulations D50 and D120. Indeed, the streamwise component of the streamwise-averaged velocity field, $\langle \overline{u}\rangle _{B,x}^{+}$ , which is plotted in figure 7, shows that the flow of the run D120 penetrates deeper into the crevices between the spheres. However, the mean square of the corresponding fluctuations, $\langle \overline{u^{\prime \prime }u^{\prime \prime }}\rangle _{B,x}^{+}$ , is more intense in the roughness sublayer over the crest of the roughness elements in the transitionally rough case than in the fully rough case (see figure 8 a,b), although significant turbulent fluctuations are present in the roughness canopy in the latter one. It also appears that the maximum of $\langle \overline{u^{\prime \prime }u^{\prime \prime }}\rangle _{B}^{+}$ is localised over the top of the spheres, where the interaction between the shear layers produced by neighbouring spheres is smaller, and not elsewhere. This can be also deduced from figure 8(c) which shows the distribution of $\langle \overline{u^{\prime \prime }u^{\prime \prime }}\rangle _{B,z}^{+}$ , where the peaks of figure 8(a) are filtered out through the spanwise average.

Figure 9.

Separation regions where the mean flow reverses are visualised by contour surfaces of the average streamwise velocity $\langle \overline{u}\rangle _{B}^{+}$ at the values $-0.02$ (yellow) and $-0.2$ (red). Rightward arrows indicate the direction of the flow. (a) Run D50; (b) run D120.

Figure 10.

Secondary flow visualised by velocity vectors $\langle \overline{(w,v)}\rangle _{B,x}^{+}$ and shaded by their modulus. The position of the spheres is outlined by red broken lines. (a) Run D50; (b) run D120.

Flow separation was observed, both for the simulations D50 and D120, downstream of the roughness elements and over the spherical caps placed directly adjacent to the wall. Figure 9 shows that the regions where the average streamwise velocity, $\langle \overline{u}\rangle _{B}^{+}$ , is negative, increase in size and intensity for increasing values of the Reynolds number. These regions are paired with as many spanwise-oriented recirculation cells from which fluid particles can escape flowing either along the upstream side of the spheres or the downstream side of the spherical caps. A mean secondary motion arises in both the transitionally and fully rough regime simulations. Here we use the term ‘secondary motion/flow’ without implying any statement on the origin of the motion, which we believe to be due to the combined effects of the three-dimensional geometry (curvature), inertia (breaking of fore/aft symmetry) and turbulence. Note that the latter effect is expected to be relatively weak, since turbulence intensity decays rapidly when entering the interstitial below the inter-particle grooves (i.e. for $y<y_{0}/2$ ). Vectors in figure 10 highlight the direction of the secondary flow projected on the $(y,z)$ -plane: while for the run D50 (in terms of streamwise-averaged flow) the fluid is directed towards the regions characterised by negative streamwise velocity in the zone denoted by ‘A’ and it moves away from the roughness through B, the opposite picture appears for the run D120, since the secondary motion converges at C and diverges from D. This can be explained by observing that the mean flow tends to penetrate deeper into the roughness along with the secondary recirculation cells (the migration of their centre of rotation can be noted in figure 10), also encountering the spherical caps mounted on the wall. A pair of recirculation cells appears in figure 10(b) which is absent in the transitionally rough simulation, and which is possibly associated with the flow interaction with the underlying spherical caps. It is interesting to note that the direction of rotation of the recirculation cells over the spherical caps is opposite with respect to those in the vicinity of the crest of the (complete) spheres, while the average streamwise velocity has also opposite sign, i.e. being negative over the caps and positive over the spheres. This picture is confirmed by the following analysis of the average vorticity field $\langle \overline{\unicode[STIX]{x1D74E}}\rangle _{B}^{+}$ .

Figure 11.

Sketch of the top (a) and side (b) view of the region (sphere-box) over which the flow field was averaged. Broken lines indicate the position of the cross-sections at which the streamwise (red broken lines), spanwise (blue broken lines) and wall-normal (black broken lines) vorticity components were computed and shown in figures 12, 13 and 14, respectively.

Figure 12.

Streamwise component $\langle \overline{\unicode[STIX]{x1D714}_{x}}\rangle _{B}^{+}$ of the sphere-box/time-averaged vorticity field visualised by filled contours (a,b) in the vertical plane through the centre of the sphere $\tilde{x}=0$ , $\tilde{x}$ being the streamwise coordinate in the frame of reference ${\mathcal{B}}$ , (c,d) at $\tilde{x}=(D+\unicode[STIX]{x1D6E5}_{B})/4$ and (e,f) at $\tilde{x}=(D+\unicode[STIX]{x1D6E5}_{B})/2$ , i.e. at the slices $x_{A}$ , $x_{B}$ and $x_{C}$ of figure 11, respectively. The value $\langle \overline{\unicode[STIX]{x1D714}_{x}}\rangle _{B}^{+}=0$ is shown as a thick line. (a,c,e) Run D50; (b,d,f) run D120.

Figure 13.

Spanwise component $\langle \overline{\unicode[STIX]{x1D714}_{z}}\rangle _{B}^{+}$ of the sphere-box/time-averaged vorticity field visualised by filled contours (a,b) at the vertical plane through the centre of the sphere $\tilde{z}=0$ , $\tilde{z}$ being the spanwise coordinate in the frame of reference ${\mathcal{B}}$ , (c,d) at $\tilde{z}=(D+\unicode[STIX]{x1D6E5}_{B})/2$ , i.e. at the slices $z_{A}$ and $z_{B}$ of figure 11, respectively. The value $\langle \overline{\unicode[STIX]{x1D714}_{z}}\rangle _{B}^{+}=0$ is shown as a thick line. The principal flow direction is indicated by the arrows. (a,c) Run D50; (b,d) run D120.

Figure 14.

Wall-normal component $\langle \overline{\unicode[STIX]{x1D714}_{y}}\rangle _{B}^{+}$ of the sphere-box/time-averaged vorticity field visualised by filled contours (a,b) at the wall-parallel plane $y=D+dx$ just over the crest of the spheres, (c,d) at $y=0.9D$ and (e,f) at $y=y_{0}$ , i.e. at the slices $y_{A}$ , $y_{B}$ and $y_{C}$ of figure 11, respectively. The value $\langle \overline{\unicode[STIX]{x1D714}_{y}}\rangle _{B}^{+}=0$ is shown as a thick line, while, in panels (e) and (f), red broken contour lines indicate $\langle \overline{u}\rangle _{B}^{+}=0$ . The principal flow direction is from left to right. (a,c,e) Run D50; (b,d,f) run D120.

The vorticity related to the average flow field, $\langle \overline{\unicode[STIX]{x1D74E}}\rangle _{B}^{+}$ , is shown in figures 12–14 in planes which are orthogonal to the respective vorticity components, as sketched in figure 11. For both simulations D50 and D120, the distribution of the streamwise vorticity component, $\langle \overline{\unicode[STIX]{x1D714}_{x}}\rangle _{B}^{+}$ , in figure 12 shows the presence of four streamwise-oriented vortical structures per streamwise-orientated inter-particle gap. In figure 12(a,b), which corresponds to a plane orthogonal to the $x$ -axis and crossing the centre of the spheres (referred to as plane $x=x_{A}$ in figure 11), these structures can be seen as thin sheets. Then, as we move to the next cross-section further downstream (figure 12 c,d related to the plane $x=x_{B}$ of figure 11), these structures remain practically attached to the spheres while the upper pair weakens and the lower one intensifies. Finally, at the vertical plane passing through the mid-plane between two adjacent spheres (figure 12 e,f) related to the plane $x=x_{C}$ of figure 11), the vortical structures still extend along the shear layer which forms over the recirculation region in the ‘wake’ of the reference sphere. Here we observe that while in case D50 the upper pair of vortices practically disappears, it is still present in case D120. It is also noteworthy that in the fully rough case, another quadruplet of streamwise vorticity structures similar to that previously described, but much weaker and of opposite sign, appears over the wall-mounted spherical caps (figure 12 f). The reversal of vorticity direction confirms the picture previously given for the recirculation cells related to the secondary flow (cf. figure 10 b).

The spanwise vorticity component, $\langle \overline{\unicode[STIX]{x1D714}_{z}}\rangle _{B}^{+}$ , is the most intense, since it is associated with the primary mean flow. The most intense regions of $\langle \overline{\unicode[STIX]{x1D714}_{z}}\rangle _{B}^{+}$ are located at and around the top of the spheres (figure 13 a,b corresponding to the plane $z=z_{A}$ of figure 11 a); this high vorticity zone extends clearly much further downstream along the shear layer in the fully rough case than at lower roughness Reynolds number. On the other hand, very small levels of $\langle \overline{\unicode[STIX]{x1D714}_{z}}\rangle _{B}^{+}$ are attained in the centre plane between the spheres (figure 13 c,d), plane $z=z_{B}$ of figure 11 a).

Finally, let us turn to the wall-normal component $\langle \overline{\unicode[STIX]{x1D714}_{y}}\rangle _{B}^{+}$ which is shown in figure 14 at the wall-parallel planes indicated in figure 11(b). Due to the deflection of the flow around the spheres, two regions characterised by high absolute values of the wall-normal vorticity component appear (with opposite sign) on either side of the upper part of the spheres. Again, the principal difference between the transitionally rough and the fully rough case is the extent to which the vorticity patches emanating from the spheres reach downstream. In the latter case this extent is significantly larger; e.g. at a wall-normal distance equal to the virtual origin (figure 14 d), the zones with intense values of the mean wall-normal vorticity reach all the way to the adjacent sphere.

The statistical footprint of the interaction of roughness-induced high vorticity structures with the surface of the spheres can be detected in the distribution of the stress on the roughness elements which is discussed in § 3.2.

Figure 15.

Distance between two adjacent low-(or high-)speed streaks $\unicode[STIX]{x1D706}_{uz}^{+}$ plotted as a function of the distance from the virtual wall, for the runs D50 (black line, empty circles) and D120 (red line, empty circles). The first point is located just above the crest of the spheres. $\unicode[STIX]{x1D706}_{uz}^{+}$ was estimated as twice the location of the minimum of the correlation function $R_{uu}(r_{z})$ of the fluctuations $u^{\prime }$ of the streamwise velocity component (small inset). Solid lines with filled symbols show the trend of $\unicode[STIX]{x1D706}_{uz}^{+}$ for the respective smooth-wall cases S1 and S2. Vertical broken lines indicate the distance of the crest of the spheres from $y_{0}^{+}$ , while horizontal broken lines indicate $D^{+}$ .

Figure 16.

Panels (a) and (b) show the trends of $\unicode[STIX]{x1D706}_{vz}^{+}$ and $\unicode[STIX]{x1D706}_{wz}^{+}$ obtained from $R_{vv}(z)$ and $R_{ww}(z)$ , in the same manner as $\unicode[STIX]{x1D706}_{uz}^{+}$ in figure 15. Symbols and lines are the same as those indicated in the caption of figure 15.

With the aim of determining the influence of the roughness on the turbulence structure, the distribution of the temporal fluctuations, ( $\boldsymbol{u}^{\prime }$ ), is investigated both for the transitionally and fully rough simulations. Typically, the presence of low-/high-speed streaks is associated with that of streamwise vortices. Indeed, the characteristics of velocity streaks are an effective indicator of the structure of turbulence in wall-bounded flows. It is well known that, in the smooth-wall case, the distance in the spanwise direction between two adjacent low- or high-speed streaks is approximately $\unicode[STIX]{x1D706}_{uz}^{+}=100$ immediately over the viscous sublayer ( $y^{+}\sim 5$ ), and grows at a constant rate with the distance from the wall until $y^{+}\sim 30$ (Kim et al. Reference Kim, Moin and Moser1987). The value of $\unicode[STIX]{x1D706}_{uz}^{+}$ can be estimated as twice the separation distance at which the spanwise two-point correlation of the streamwise velocity fluctuations, $R_{uu}$ , attains the maximum negative value (see the small panel of figure 15), as long as the minimum is negative. On the basis of this definition, Kim et al. (Reference Kim, Moin and Moser1987) extended the idea of a spatial scale $\unicode[STIX]{x1D706}_{uz}^{+}$ less intuitively to the other components of the velocity and found that $\unicode[STIX]{x1D706}_{wz}^{+}\simeq \unicode[STIX]{x1D706}_{uz}^{+}$ and $\unicode[STIX]{x1D706}_{vz}^{+}\simeq 2\unicode[STIX]{x1D706}_{uz}^{+}$ . Presently, a similar approach was adopted for both the smooth- and rough-wall simulations (see figure 16). As shown in figures 15, 16(a) and 16(b), respectively, the values of $\unicode[STIX]{x1D706}_{uz}^{+}$ , $\unicode[STIX]{x1D706}_{vz}^{+}$ and $\unicode[STIX]{x1D706}_{wz}^{+}$ , at $y^{+}=D^{+}$ for the simulations D50 and D120 were approximately the same as those attained for the respective smooth-wall cases at the same distance from the virtual wall, i.e. at $y_{smooth}^{+}=(D^{+}-y_{0}^{+})_{rough}$ . Figure 17 shows the distribution of instantaneous fluctuations of the streamwise velocity, $u^{\prime }$ , in the wall-parallel plane located at $y^{+}\simeq D^{+}$ for the run D120 (panel (a)) and at $y^{+}\simeq D^{+}-y_{0}^{+}$ for the smooth-wall simulation at the same bulk Reynolds number (panel (b)). Visually, the length and width of low-speed streaks in those planes are remarkably similar, as is their intensity. The value of $\unicode[STIX]{x1D706}_{wz}^{+}$ , extrapolated at $y_{smooth}^{+}=(D^{+}-y_{0}^{+})_{D120}$ for the smooth-wall simulation at $Re_{bH}=6900$ , is almost equal to that of $\unicode[STIX]{x1D706}_{wz}^{+}$ obtained for the fully rough simulation D120 (see figure 15 c). Hence, the relationships observed by Kim et al. (Reference Kim, Moin and Moser1987) between $\unicode[STIX]{x1D706}_{uz}^{+}$ and $\unicode[STIX]{x1D706}_{vz}^{+}$ and between $\unicode[STIX]{x1D706}_{uz}^{+}$ and $\unicode[STIX]{x1D706}_{wz}^{+}$ still hold both in the transitionally and fully rough regimes for the present shape and arrangement of the roughness elements, and an increase of the Reynolds number $D^{+}$ results in the increase of the slope of $\unicode[STIX]{x1D706}_{uz}^{+}$ as a function of $y^{+}$ . However, the methodology described above does not allow us to estimate the value of $\unicode[STIX]{x1D706}_{z}^{+}$ below the crest of roughness elements where velocity fluctuations are not defined continuously on wall-parallel planes.

Figure 17.

Instantaneous visualisation of velocity fluctuations $u^{\prime }$ for (a) the simulation D120 at $y^{+}-y_{0}^{+}=20$ and (b) the smooth-wall simulation with $Re_{bH}=6900$ at $y^{+}=20$ . White broken lines delimit the original computational domain, extended periodically in the streamwise and spanwise direction up to the dimensions of the domain of the simulation D120 in order to facilitate the visual comparison. The principal flow direction is from left to right.

It emerges from this picture that the effect of roughness is not to change sharply the structure of turbulence above the roughness elements, but to select a particular scale related to the roughness which then appears more pronounced than in an equivalent turbulent open-channel flow developing over a smooth wall at the same bulk Reynolds number. In the fully rough simulation, the distance between low- and high-speed streaks, $\unicode[STIX]{x1D706}_{uz}^{+}/2$ , at $y^{+}=D^{+}$ , was equal to 85.

3.2 Force and torque acting on the roughness elements

The stress acting on each roughness element was calculated on the basis of the volume forces associated with the immersed boundary method (Chan-Braun et al. Reference Chan-Braun, García-Villalba and Uhlmann2011). The square arrangement of the spheres allows us to refer to them with the indices $(i,j)$ as in a two-dimensional array, where $i=1,\ldots ,n_{c}$ and $j=1,\ldots ,n_{r}$ enumerate the streamwise indices (columns) and the spanwise indices (rows) of the array, respectively. For the present simulations $n_{c}=64$ and $n_{r}=16$ . We apply the sphere-box/time-average operator described in § 3.1 to the total stress tensor $\unicode[STIX]{x1D749}$ (containing both the viscous and pressure contributions, as defined in (3.5)). At the sphere surface, the projection of the average stress tensor upon the outward-facing normal vector $\widetilde{\boldsymbol{n}}$ yields the average stress vector, viz.

Then the sphere-box/time-averaged hydrodynamic force is defined as:

while the time average of the force acting on the $(i,j)$ th sphere is:

where $\widetilde{{\mathcal{S}}}$ denotes the surface of the sphere contained in the sphere-box ${\mathcal{B}}$ (see figure 2 and the definition (A 8) in the appendix A), ${\mathcal{S}}^{(i,j)}$ indicates the surface of the $(i,j)$ th sphere centred in $\boldsymbol{x}_{c}^{(i,j)}$ and $\boldsymbol{n}$ the surface-normal unit vector.

The fluctuations $\boldsymbol{F}^{\prime \prime }$ around $\boldsymbol{F}$ can be defined in a similar way as those of the velocity given in (3.2). The components $(F_{x},F_{y},F_{z})$ of the force $\boldsymbol{F}$ are the mean drag, vertical lift and lateral lift forces acting on roughness elements, respectively, which can be normalised by the reference force $F_{R}=\unicode[STIX]{x1D71A}u_{\unicode[STIX]{x1D70F}}^{2}L_{B}^{2}$ . The values of the mean dimensionless force components for the simulations D50 and D120 are provided in table 3 along with their second, third and fourth moment statistics, the angle $\unicode[STIX]{x1D6FC}_{F}$ between the average drag and lift forces, and the values of the ratios of the viscous components of drag and lift forces, $F_{x}^{(\unicode[STIX]{x1D708})}$ and $F_{y}^{(\unicode[STIX]{x1D708})}$ , to the respective total forces. The mean dimensionless drag force remains essentially unchanged between the simulations D50 and D120 as well as the value of $U_{bh}/u_{\unicode[STIX]{x1D70F}}$ (see table 1). This can be explained in terms of the Darcy–Weisbach friction factor, since $H/D$ (which is the inverse of the relative roughness) had nearly the same value in the two simulations, while the values of bulk Reynolds number were sufficiently close and high to keep the friction factor almost constant. Contrarily, the mean vertical lift experiences a significant increase, resulting in the increase of the angle $\unicode[STIX]{x1D6FC}_{F}$ . Chan-Braun et al. (Reference Chan-Braun, García-Villalba and Uhlmann2011) observed that $F_{y}/F_{R}$ increased as an effect of increasing $D^{+}$ , but since the relative submergence $H/D$ was simultaneously changed, they could not conclude on a pure scaling with particle Reynolds number $D^{+}$ . Despite the fact that the dimensionless drag remains essentially unchanged from D50 to D120, the pressure contribution to the mean drag force significantly increases (table 3). To a lesser extent this is also true for mean wall-normal lift. The kurtosis of the three components of $\boldsymbol{F}^{\prime \prime }$ decreases for increasing values of $D^{+}$ , tending towards the value 3 expected for a Gaussian distribution. Indeed, fluctuations of small length scale, at large values of $D^{+}$ , are filtered out through the surface integral, thereby reducing the probability of values far from the mean and consequently the value of  ${\mathcal{K}}$ .

Table 3.

Statistics of the hydrodynamic force acting on the top layer particles in runs D50 and D120. $F_{i}/F_{R}$ denotes the normalised mean force component in the $x_{i}$ -direction,  $\unicode[STIX]{x1D6FC}_{F}=\arctan ({\mathcal{C}}_{F_{y}}/{\mathcal{C}}_{F_{x}})$ denotes the angle of the resulting force with respect to the $x$ -axis,  $\unicode[STIX]{x1D70E}_{F_{i}}$ , ${\mathcal{S}}_{F_{i}}$ and ${\mathcal{K}}_{F_{i}}$ are the standard deviation, the skewness and the kurtosis of the $i$ th component of the force, respectively. The reference force is defined as $F_{R}=\unicode[STIX]{x1D71A}u_{\unicode[STIX]{x1D70F}}^{2}L_{B}^{2}$ .

The average torque acting on the roughness elements is defined as follows:

where the distance vector is denoted by $\widetilde{\boldsymbol{r}_{c}}=(\widetilde{\boldsymbol{x}}-\widetilde{\boldsymbol{x}}_{c})$ , with $\widetilde{\boldsymbol{x}}_{c}=(0,D/2,0)$ . For the sake of completeness, the values of $\boldsymbol{T}$ normalised by the reference torque $T_{R}=F_{R}(y_{0}-D/2)$ , as well as the statistics of torque fluctuations $\boldsymbol{T}^{\prime \prime }$ , were computed and reported in table 4 for the simulation D120. The values of the standard deviation of torque fluctuations in the three directions, $(\unicode[STIX]{x1D70E}_{T_{x}},\unicode[STIX]{x1D70E}_{T_{y}},\unicode[STIX]{x1D70E}_{T_{z}})/T_{R}$ , are nearly the same whereas the counterparts for force fluctuations are significantly different. This suggests that the fluctuations of the shear stress (associated with viscous effects) are less anisotropic than pressure fluctuations, which dominate the hydrodynamic force and which do not contribute to the torque. This is true despite the fact that only the most exposed part of the sphere surface contributes significantly to the surface shear stress, i.e. this is essentially the upstream-facing part of the upper hemisphere (figure not shown).

Table 4.

Statistics of torque acting on the roughness elements for cases D50 and D120. $T_{i}/F_{R}$ denotes the normalised mean torque component in the $x_{i}$ -direction,  $\unicode[STIX]{x1D70E}_{T_{i}}$ , ${\mathcal{S}}_{F_{i}}$ and ${\mathcal{K}}_{T_{i}}$ are the standard deviation, the skewness and the kurtosis of the $i$ th component of the torque, respectively.

Figure 18.

Distribution of the mean pressure, $\langle \overline{p}_{tot}\rangle _{B}$ , evaluated at the sphere surface and normalised by $F_{R}/A_{sph}$ . (a,c) Show the top view while (b,d) show the side view of the sphere. (a,b) Run D50; (c,d) run D120. Thin contour lines at values $\pm [3,6,9]$ . The thick contour line indicates the value 0.

Figure 19.

Root mean square of the pressure fluctuations $p_{tot}^{\prime \prime }$ at the sphere surface, normalised by $F_{R}/A_{sph}$ . (a,b) Show the top and side views of the sphere for the run D120. Contour lines are equispaced by 0.5.

The distribution of the average surface pressure $\langle \overline{p}_{tot}\rangle _{B}$ (i.e. the surface-normal component of the stress vector $\widetilde{\unicode[STIX]{x1D749}}_{n}$ with the sign reversed), normalised by $F_{R}/A_{sph}$ , where $A_{sph}=\unicode[STIX]{x03C0}D^{2}$ , is shown in figure 18. Two regions of strongly negative and positive mean pressure appear on the upper hemisphere (please recall the definition of pressure in (3.6)). While the negative-valued region is slightly shifted downstream of the sphere, the positive region (i.e. with the normal stress directed towards the centre of the sphere) is located somewhat upstream. In case D50 this latter (high pressure) region is almost confined above $y=y_{0}$ , while it is elongated in the spanwise direction in case D120, reaching to smaller wall-normal distances on the sides of the sphere. As it was previously noted in § 3.1, the mean flow penetrates deeper through the grooves between the roughness elements for increasing values of $Re_{bH}$ , causing the steady spanwise vorticity structures to squeeze on the top of the spheres. This explains the lateral spreading of the two regions in figure 18(c,d) characterised by peaks of the average pressure. By virtue of the distribution shown in figure 18(c,d), it is possible to deduce that, for the square arrangement of spherical roughness elements, an effective position for two pressure probes in an analogous experiment should be in the centre of the two regions presently highlighted.

For the simulation D120, a region of intense pressure fluctuations was detected on the upper upstream quarter of the sphere surface, in correspondence to the region where the average pressure $\langle \overline{p}_{tot}\rangle _{B}$ is positive. More specifically, figure 19 shows the root mean square value of the pressure fluctuations around the time-and-sphere-box average, $p_{tot}^{\prime \prime }$ , at the sphere surface. The region of high pressure fluctuation intensity in case D120 is found to coincide largely with the high mean pressure region observed in figure 18(c,d), while the fluctuations are of relatively small intensity over a considerably large region around the top of the sphere. Although this evidence was not exhaustively investigated, it could be inferred that in the fully rough regime and for the present arrangement of the spheres, most of the sphere–turbulence interaction occurs in the upstream region at the top of the sphere, which might be a signature of the vortices shed from the sphere located just upstream.

Figure 20.

Distribution of the streamwise component of $\widetilde{\unicode[STIX]{x1D749}}_{n}$ normalised by $F_{R}/A_{sph}$ . (a,c) Show the top view while (b,d) show the side view of the sphere. Black crosses and blue broken lines indicate the distance from the wall at which $\unicode[STIX]{x1D70F}_{tot}=0$ and the values $|\unicode[STIX]{x1D749}_{n}-(\unicode[STIX]{x1D749}_{n}\boldsymbol{\cdot }\boldsymbol{n})\boldsymbol{n}|A_{sph}/F_{R}=10^{-2}$ at the sphere surface, respectively. (a,b) Run D50; (c,d) run D120. Thin contour lines at values $[-0.5,0.5,3.5,5.0,7.5]$ . The thick contour line indicates the value 0.

Figure 20 shows the projection of the surface stress vector $\widetilde{\unicode[STIX]{x1D749}}_{n}$ upon the streamwise direction, i.e. the map of the local surface stress contribution to drag, $\widetilde{\unicode[STIX]{x1D749}}_{n}\boldsymbol{\cdot }\boldsymbol{e}_{x}$ . It can be seen that in case D120 (figure 20 c,d) this quantity is largest in roughly the same region on the upper, upstream part of the sphere, where the mean surface pressure is large and positive (figure 18 c,d), and where the surface pressure fluctuations were found to be large (figure 19). In the transitionally rough case D50, on the contrary, the maximum contribution to drag is provided by the shear stress near the crest of the roughness elements (figure 20 a,b), in a zone less affected by high pressure (figure 18 a,b). Therefore, although the values of the dimensionless drag force and the intensity of its fluctuations are similar in the transitionally and in the fully rough regimes, they clearly originate from different processes: the first is associated with the skin friction, while the latter is due to pressure. This picture can be interpreted as a manifestation of a ‘form-drag mechanism’ in case D120. It is also supported by the considerations previously formulated in § 3.1 about the shift of $\langle \overline{u^{\prime \prime \prime }u^{\prime \prime \prime }}\rangle ^{+}$ beneath the crest of the spheres in the fully rough regime.

Figure 21.

Cumulative streamwise (——) and wall-normal (– – –) components of $f_{\unicode[STIX]{x1D6FC}}(y)$ normalised by $F_{\unicode[STIX]{x1D6FC}}$ , where $\unicode[STIX]{x1D6FC}$ is replaced by $x$ and $y$ , respectively. The horizontal broken line is located at $y_{0}/D$ . The lines $(\cdot \cdot$ ○ $\cdot \cdot )$ and $(\cdot \cdot +\cdot \cdot )$ indicate the viscous contributions $f_{\unicode[STIX]{x1D6FC}}^{(\unicode[STIX]{x1D708})}(y)$ to the cumulative streamwise and wall-normal force, respectively. (a) Run D50; (b) run D120.

In order to further investigate the contributions of pressure and viscous stresses to the net force on the spheres, let us define a cumulative force $\boldsymbol{f}(y)$ as the integral of the stress vector over the sphere surface up to height $y$ , viz.

where $\unicode[STIX]{x1D702}$ denotes the wall-normal coordinate defined in the range $[0,D]$ and $\unicode[STIX]{x1D703}$ the azimuthal angle. Obviously $\boldsymbol{f}(D)=\boldsymbol{F}$ , as defined in (3.13). The streamwise and wall-normal components of $\boldsymbol{f}(y)$ , normalised by the respective components of the total force $\boldsymbol{F}$ , are shown in figure 21. For the purpose of comparison, the same figure also shows the viscous contributions $\boldsymbol{f}^{(\unicode[STIX]{x1D708})}(y)$ to the cumulative force in the same scaling. It can be observed that the cumulative viscous contributions to both drag and lift are negligible for small wall distances. At larger wall distances (for $y\gtrsim 0.6D$ in case D50 and for $y\gtrsim 0.7D$ in case D120) the viscous contribution to drag increases monotonically up to the final value listed in table 3. Contrarily, in both cases the cumulative viscous contribution to lift is negative for small wall distances, and changes its sign only very close to the top of the spheres, most of the positive lift contribution being generated for $y\gtrsim 0.85D$ . Note that the different distribution of the pressure with the distance from the wall in the present cases does not manifest itself in any significant difference in the normalised cumulative lift force profile, even though the vertical dimensionless lift $F_{y}/F_{R}$ is much larger in case D120.

Figure 22.

Standard deviation of the torque acting on a smooth square tile of side $s^{+}$ , with an ideal arm of length $s^{+}/2$ , originated by the shear stress, $\unicode[STIX]{x1D748}_{{\mathcal{T}}}$ (solid lines), or by the shear stress and pressure, $\unicode[STIX]{x1D70E}_{{\mathcal{T}}_{x}^{(p)}}$ (broken lines), normalised by $\unicode[STIX]{x1D71A}u_{\unicode[STIX]{x1D70F}}^{2}s^{3}/2$ . The standard deviation of the torque acting on the sphere $\widetilde{{\mathcal{S}}}$ , normalised by $T_{R}$ , is also visualised with full symbols for the simulations (a) D50 and (b) D120. Symbols indicate the streamwise (—▵—, ▴), the wall-normal (—▿—, ▾) and spanwise (—▫—, ▪) components of torque. (a): $Re_{bH}\sim 2900$ , run D50; (b): $Re_{bH}\sim 6900$ , run D120.

Chan-Braun et al. (Reference Chan-Braun, García-Villalba and Uhlmann2011) have introduced an analogy between the force and torque acting upon a square region of the wall in smooth-wall channel flow and the corresponding force/torque components acting upon a wall-mounted sphere (cf. the sketch in their figure 9). This ‘smooth-wall analogy’ was formulated on the basis of the observation that in the transitionally rough regime essentially only the crests of the spheres were exposed to significant flow, and that viscous effects were dominating the surface stress. Hence, in their case the action of turbulent fluctuations on the upper part of the sphere surface could be modelled by analogy with that occurring on a smooth-wall tile $\unicode[STIX]{x1D6E4}_{s}$ of side length $s$ comparable to the sphere diameter $D$ . In order to test to which extent this model still holds in the fully rough regime, we use the data from the simulation of open-channel flow over a smooth wall (at the same bulk Reynolds number as D120) in order to compute the statistics of the wall force and torque fluctuations as functions of $s$ . In particular, the standard deviations of the torque components $\unicode[STIX]{x1D70E}_{T_{x}}/T_{R}$ , $\unicode[STIX]{x1D70E}_{T_{y}}/T_{R}$ and $\unicode[STIX]{x1D70E}_{T_{z}}/T_{R}$ for the simulations D50 and D120 are compared with those of the torque:

and

obtained for the smooth-wall square tile $\widetilde{\unicode[STIX]{x1D6E4}}_{L_{B}}$ (i.e. of side $s=L_{B}$ ) about the ideal point $\widetilde{\boldsymbol{x}}_{s}=(0,-s/2,0)$ such that $\widetilde{\boldsymbol{r}}_{s}=\widetilde{\boldsymbol{x}}-\widetilde{\boldsymbol{x}}_{s}$ with $\widetilde{\boldsymbol{x}}\in \widetilde{\unicode[STIX]{x1D6E4}}_{L_{B}}$ . The symbol $\boldsymbol{e}_{y}$ in (3.17)–(3.18) indicates the wall-normal unit vector of the canonical base while $\unicode[STIX]{x1D749}$ and $\unicode[STIX]{x1D749}_{\unicode[STIX]{x1D708}}$ denote the total stress tensor and the viscous stress tensor (i.e. without the pressure contribution), respectively, cf. the definition in (3.5). The comparison is visualised in figure 22, where the full symbols indicate the values related to the simulations D50 and D120, respectively. It is to be pointed out that the contribution of pressure fluctuations to $\boldsymbol{{\mathcal{T}}}$ was not accounted for (solid lines in figure 22) in order to preserve the analogy, since pressure fluctuations do not contribute to torque on a sphere. In fact, the trend for the standard deviations $\unicode[STIX]{x1D70E}_{{\mathcal{T}}_{x}^{(p)}}(L_{B}^{+})$ and $\unicode[STIX]{x1D70E}_{{\mathcal{T}}_{z}^{(p)}}(L_{B}^{+})$ , obtained considering also the contribution of pressure, normalised by $\unicode[STIX]{x1D71A}u_{\unicode[STIX]{x1D70F}}s^{3}/2$ , are definitely far from $T_{x}$ and $T_{z}$ obtained for the spheres (broken lines in figure 22). Indeed, this simple approach was found surprisingly satisfactory for flows in the transitionally rough regime, but it can be seen that it fails for the simulation D120, because the mechanism of sphere–flow interaction in the fully rough regime is not correctly interpreted by the model. In particular, figure 22(b) shows graphically that the standard deviations of the three components of torque, $\unicode[STIX]{x1D748}_{T}/T_{R}$ , tend to collapse on the same value, as discussed above. In other terms, the failure of the ‘smooth-wall analogy’ in the fully rough regime can be interpreted as an effect of the destruction of the buffer layer and the disappearance of significant viscous effects in the vicinity of the sphere crests.

Figure 23.

Distribution of the $y$ -component of $\widetilde{\unicode[STIX]{x1D749}}_{n}$ normalised by $F_{R}/A_{sph}$ . (a,c) Show the top view while (b,d) show the side view of the sphere. (a,b) Run D50; (c,d) run D120. Thin contour lines at values $\pm [0.5,1.5,4.5,7.5,10.5]$ . The thick contour line indicates the value 0.

Figure 18 showed that most of the lift force can be attributed to the contributions of pressure in the region near the top of the sphere for both the simulations D50 and D120. In fact, the $y$ -component of $\widetilde{\unicode[STIX]{x1D749}}_{n}$ (this is the global wall-normal direction), which is visualised in figure 23, exhibits essentially the same distribution as the pressure.

Figure 24.

Distribution of the spanwise component of $\widetilde{\unicode[STIX]{x1D749}}_{n}$ normalised by $F_{R}/A_{sph}$ . (a,c) Show the top view while (b,d) show the side view of the sphere. (a,b) Run D50; (c,d) run D120. Thin contour lines at values $\pm [1,2,3,4]$ . The thick contour line indicates the value 0.

Finally, let us return to the discussion of the effect which the structure of the average flow field in the vicinity of the spheres has upon the stress distribution on the spheres’ surface. Please recall the strong wall-normal vorticity structures described in § 3.1 and shown in figure 14. Here we relate these features to the distribution of the spanwise component of $\widetilde{\unicode[STIX]{x1D749}}_{n}$ on the upper sides of the sphere, which is shown in figure 24. The latter figure confirms that these average vortical structures, which are more intense at higher Reynolds numbers, and which tend to adhere closer to the sphere surface, cause stronger lateral force contributions locally. From these considerations it can also be understood that the intensity of turbulent fluctuations of the lateral force increases from D50 to D120, resulting in the observed large values of $\unicode[STIX]{x1D70E}_{F_{z}}/F_{R}$ and of $\unicode[STIX]{x1D70E}_{T_{y}}/T_{R}$ (cf. tables 3 and 4). In particular, the increase of the value of $\unicode[STIX]{x1D70E}_{T_{y}}/T_{R}$ in the fully rough regime along with the decrease of $\unicode[STIX]{x1D70E}_{T_{z}}/T_{R}$ provide further pieces of evidence of the breakdown of the ‘smooth-wall analogy’.

3.3 What can be inferred from the force acting on the roughness elements?

In the following, a comparison of the results obtained in the present DNS with those from laboratory experiments is performed. As already pointed out, it is still not entirely clear how the details of the flow structure scale with the bulk Reynolds number, the particle Reynolds number and with the relative submergence. However, Amir et al. (Reference Amir, Nikora and Stewart2014), who experimentally investigated the forces acting on wall-mounted spheres in hexagonal arrangement, identified quantities that are almost independent of some of these parameters in the fully rough regime. These authors estimated the temporal cross-correlation between neighbouring spheres on the basis of differential pressure measurements. In particular, they estimated the root mean square of the lift force fluctuations on the basis of the values of the difference between pressure measured at the top and at the bottom of the spheres (as indicated with $A$ and $B$ in figure 19), and showed that it was almost directly proportional to $u_{\unicode[STIX]{x1D70F}}^{2}$ . By normalising the root mean square value of pressure fluctuations with $\unicode[STIX]{x1D71A}u_{\unicode[STIX]{x1D70F}}^{2}$ they defined the coefficient $k_{L}$ , which was independent of $U_{bh}/u_{\unicode[STIX]{x1D70F}}$ and also of $H/D$ for $H/D>5$ , and found that it was equal to 3.4. In order to facilitate the comparison with the data of Amir et al. (Reference Amir, Nikora and Stewart2014), we have computed the value of the coefficient $k_{L}$ in an analogous manner (i.e. as the normalised root mean square pressure difference between the two poles of the sphere at points $A$ , $B$ in figure 19), and found that $k_{L}=1.1$ in case D120. It is reasonable that the present value is smaller than the value obtained by Amir et al. (Reference Amir, Nikora and Stewart2014), because the present spheres in a square arrangement shelter each other more strongly than in the hexagonal arrangement investigated by those authors.

A quantity which can be expected to be more robust with respect to variations of the arrangement of the roughness elements is the convection velocity $U_{c}$ of larger vortex structures. Amir et al. (Reference Amir, Nikora and Stewart2014) found a fair correlation of drag and lift between spheres separated by a distance of either 1 or 2 diameters according to the availability of instrumented neighbouring spheres (only a few spheres were instrumented with one-dimensional pressure probes each). Mimicking the technique described by Amir et al. (Reference Amir, Nikora and Stewart2014) (which is different from that used by Chan-Braun et al. Reference Chan-Braun, Garcia-Villalba and Uhlmann2013), the convection velocity $U_{c}$ of force fluctuations which are induced by the interaction of roughness elements with turbulent vortex structures, was estimated on the basis of the lag of the peak of the temporal cross-correlation of drag force between neighbouring spheres aligned along the streamwise direction. In particular, $U_{c}$ is found to measure $0.72U_{bH}$ and $0.64U_{bH}$ for the runs D50 and D120, respectively. These values are in fair agreement with the value $0.66U_{bH}$ obtained by Amir et al. (Reference Amir, Nikora and Stewart2014) for their experiments characterised by the relative submergence $H/D=5$ . Since Amir et al. (Reference Amir, Nikora and Stewart2014) showed that the speed $U_{c}$ was essentially independent of $U_{bh}/u_{\unicode[STIX]{x1D70F}}$ (see their figure 19b), the fact that the Reynolds number in the present DNS is substantially lower should not matter. Figure 25 shows the temporal cross-correlation of drag and lift force fluctuations, denoted by $R_{F_{x}^{\prime }F_{y}^{\prime }}(\unicode[STIX]{x0394}r_{x}=iL_{B},\unicode[STIX]{x0394}r_{z}=jL_{B},\unicode[STIX]{x1D6E5}_{T})$ with $i=0,\ldots ,n_{c}-1$ and $j=0,\ldots ,n_{r}-1$ , for the runs D50 and D120 (thick lines) as well as for the respective simulations performed over a smooth wall (thin lines). In the spirit of the smooth-wall analogy discussed in § 3.2, the lift force is associated with pressure fluctuations while drag force fluctuations are associated with those of the viscous shear stress. Since Amir et al. (Reference Amir, Nikora and Stewart2014) could equip each sphere with a one-direction pressure probe only, they could estimate $R_{F_{x}^{\prime }F_{y}^{\prime }}$ only from measurements made on distinct spheres, separated by a non-negative distance either in the streamwise or spanwise directions. According to their hexagonal arrangement, the spanwise direction was more convenient. Presently, $R_{F_{x}^{\prime }F_{y}^{\prime }}(0,L_{B},\unicode[STIX]{x1D6E5}_{T})$ was computed for the simulations D50 and D120 (see broken lines in figure 25). Although the present value of $R_{F_{x}^{\prime }F_{y}^{\prime }}$ is approximately two times larger than that observed by Amir et al. (Reference Amir, Nikora and Stewart2014) (possibly due to the different geometrical configuration of the roughness), the time lags of the negative (or positive) peaks are found to be essentially the same, equal to $0.14U_{bH}/H$ and to $0.16U_{bH}/H$ , respectively. Moreover, the single-sphere cross-correlation function $R_{F_{x}^{\prime }F_{y}^{\prime }}(0,0,\unicode[STIX]{x1D6E5}_{T})$ was computed (see solid thick lines in figure 25) which revealed the same time lag of the positive (negative) maximum ( $\unicode[STIX]{x1D6E5}_{T}=0.16H/U_{bH}$ and $\unicode[STIX]{x1D6E5}_{T}=-0.14H/U_{bH}$ for cases D50 and D120), and the presence of a secondary negative (positive) peak at $\unicode[STIX]{x1D6E5}_{T}=0.5H/U_{bH}$ ( $\unicode[STIX]{x1D6E5}_{T}=-0.5H/U_{bH}$ ). The trend of $R_{F_{x}^{\prime }F_{y}^{\prime }}$ indicates that a positive (negative) fluctuation of lift is most probably preceded by a positive (negative) fluctuation of drag and followed by a negative (positive) fluctuation of drag, but also that lift fluctuations are practically uncorrelated with any drag fluctuation at the same instant. Similar statistics were observed also by Hofland (Reference Hofland2005), Dwivedi (Reference Dwivedi2010) as well as by Amir et al. (Reference Amir, Nikora and Stewart2014), while the cross-correlation function for case D50 was also shown in figure 8 of Chan-Braun et al. (Reference Chan-Braun, Garcia-Villalba and Uhlmann2013).

Figure 25.

Cross-correlation $R_{F_{x}^{\prime }F_{y}^{\prime }}(0,\unicode[STIX]{x0394}r_{z},\unicode[STIX]{x1D6E5}_{T})$ between drag and lift acting on the spheres as a function of time and for spanwise separation $\unicode[STIX]{x0394}r_{z}=0$ (thick solid lines) and $\unicode[STIX]{x0394}r_{z}=L_{B}$ (thick broken lines) for the runs D50 (black) and D120 (red). Thin lines indicate the function $R_{F_{x}^{\prime }F_{y}^{\prime }}$ referred to a square element of smooth wall (in absence of the spheres) of side length $s^{+}=51$ (red) and $s^{+}=129$ (black) for the smooth-wall simulations at $Re_{bH}\sim 2900$ and ${\sim}6900$ , respectively. For the sake of reference in the text, symbols $E$ and $F$ indicate the location of the zero crossings for positive delay $\unicode[STIX]{x1D6E5}_{T}$ in these two latter curves.

In case D120 the average time interval between positive and negative fluctuations of lift, equal to $\unicode[STIX]{x1D6E5}_{T}=0.33H/U_{bH}$ , was obtained from the occurrence of the global minimum of the auto-correlation function of lift, $R_{F_{y}F_{y}}(0,0,\unicode[STIX]{x1D6E5}_{T})$ (not shown here). In Chan-Braun et al. (Reference Chan-Braun, Garcia-Villalba and Uhlmann2013) the value of $\unicode[STIX]{x1D6E5}_{T}=0.42H/U_{bH}$ is measured for their case D50. Consistently, nearly after the same time interval ( $\unicode[STIX]{x1D6E5}_{T}=0.37H/U_{bH}$ ), the cross-correlation $R_{F_{x}^{\prime }F_{y}^{\prime }}$ is found to vanish in both cases D50 and D120. Furthermore, the symmetry of $R_{F_{x}^{\prime }F_{y}^{\prime }}$ with respect to the point $\unicode[STIX]{x1D6E5}_{T}=0$ suggests that forthcoming and past fluctuations have similar spatial and temporal scales as well as the same intensity. The fact that highly auto-correlated lift fluctuations and cross-correlated drag–lift fluctuations are synchronised suggests that they are associated with the same turbulent event. Thus, a linear relationship between length and time scales of turbulent fluctuations through the convection velocity $U_{c}$ , i.e. the validity of Taylor’s frozen turbulence hypothesis, can be assumed (e.g. Chan-Braun et al. Reference Chan-Braun, Garcia-Villalba and Uhlmann2013). Thereby, it is possible to estimate the length scale, (henceforth denoted as $\ell$ ) as the product of the time interval $\unicode[STIX]{x1D6E5}_{T}=0.3H/U_{bH}$ between drag (or lift) fluctuations of opposite sign and the convection velocity $U_{c}$ previously estimated. For the runs D50 and D120, $\ell$ was equal to $0.22H$ ( $1.24D$ ) and $0.19H$ ( $1.03D$ ), respectively. Furthermore, in the spanwise direction the two-point correlation of lift and drag was found negligible at the separation distance $2D$ in both the cases D50 and D120 (not shown here). Finally, since $\ell$ depends on $U_{c}$ , which is in turn affected by the value of $H/D$ , the drag-lift cross-correlation function is possibly affected by  $H/D$ .

Therefore, although the evolution of the coherent turbulent structures propagating over the roughness elements has not been presently studied, the curves of figure 25 suggest that the fluctuations of drag and wall-normal lift forces are predominantly caused by the interaction between the roughness elements and vortices of specific size. In particular, the size of these vortices is comparable to the length scale $\ell$ which we have estimated in the aforementioned way, and, therefore, to the size of the spheres.

Figure 26.

Conceptual scheme of the sequence of positive-to-negative drag and lift fluctuations ideally induced by disturbances propagating with the speed $U_{c}$ . The values of time intervals were obtained for the simulation D120.

With the purpose of helping the interpretation of figure 25, a conceptual model of the response of the spheres to the action of these vortices is sketched in figure 26 where, in a schematic side view, spheres are represented by circles, vortices by circular arrows and hydrodynamic forces by vectors applied to the sphere centre. In this model, turbulent coherent structures induce pulsations of the intensity of the vortices characterising the average flow field. These pulsations propagate downstream with speed $U_{c}$ and produce on the spheres the sequence of fluctuations of drag and wall-normal lift selected by the cross-correlation curves of figure 25. Vortices involved in the intensity pulsation at different propagation phases are highlighted in black in the sequences of figure 26. Hence, the period of the pulsation can be defined as the time interval between two consecutive positive (or negative) peaks of the cross-correlation function in figure 25, which measures $0.64H/U_{bH}$ . Let us note that disturbances which do not result in highly cross-correlated fluctuations of the net drag and lift forces acting on the spheres are not considered by the model. The conceptual model of interaction between turbulence structures and roughness elements sketched in figure 26 is consistent with the statistical evidences related to the hydrodynamic forces and with the topology of the average vorticity field described above.

However, so far our argument only relates to pulsations of the mean flow in the vicinity of the spheres, not to the actual coherent structures which cause the pulsations. In order to go further, a detailed analysis of the correlation between the particle forces and the flow field should be carried out, as has been done by Chan-Braun et al. (Reference Chan-Braun, Garcia-Villalba and Uhlmann2013) for the transitionally rough case. In the present case, one should specifically identify those turbulent events which have a maximum correlation with the occurrence of a strong positive lift force fluctuation, and which are the same events that have a maximum correlation with a strongly positive (negative) particle drag fluctuation at a fixed negative (positive) delay time of $\unicode[STIX]{x1D6E5}_{T}\sim 0.15H/U_{bH}$ . This additional analysis, which is beyond the scope of the present work, would yield information on the flow structure responsible (in an average sense) for the said cross-correlation.

Finally, the cross-correlation function of drag and lift, $R_{F_{x}^{\prime }F_{y}^{\prime }}$ , was also calculated for the smooth-wall simulations at $Re_{bH}=2900$ and 6900 in order to estimate length and time scales of force fluctuations acting on a square tile (with side length $L_{B}$ ) of a smooth channel wall at a corresponding Reynolds number. The thin solid lines of figure 25 show the cross-correlation for the smooth-wall cases. Both the curves for the smooth-wall simulations at $Re_{bH}=2900$ and 6900 show a negative correlation at the time interval $\unicode[STIX]{x1D6E5}_{T}=0$ , which means that a lift fluctuation is most probably associated with a drag fluctuation of opposite sign. Indeed, these combinations bring to mind the action of sweep-like and ejection-like events which dominate the Reynolds stress statistics in a turbulent wall-bounded flow over a smooth wall. As opposed to the rough-wall cases, the cross-correlation of lift and drag for the smooth-wall simulations do not exhibit the symmetry with respect to the origin of the axes, suggesting that forthcoming and past turbulent events are in general characterised by different length and time scales. However, from the diagrams of figure 25, it is possible to distinguish large turbulent structures, characterised by length scale much larger than $L_{B}$ , from ‘smaller’ turbulent structures characterised by the time scale $\unicode[STIX]{x1D6E5}_{T}^{(F)}-\unicode[STIX]{x1D6E5}_{T}^{(E)}$ equal to 0.30 and 0.36 for the simulations at $Re_{bH}=2900$ and 6900, respectively, where $\unicode[STIX]{x1D6E5}_{T}^{(E)}$ and $\unicode[STIX]{x1D6E5}_{T}^{(F)}$ denote the time intervals at the points $E$ and $F$ , marked by a star in figure 25, when the cross-correlation function vanishes. Then, by using the same procedure described for the rough bottom simulations, the convective velocity of turbulent structures $U_{c}$ was estimated for the smooth-wall cases equal to $0.71U_{bH}$ and $0.69U_{bH}$ , respectively. Hence, the length scale of the ‘smaller’ turbulent structures was found equal to $0.20H$ and to $0.25H$ for the simulations at $Re_{bH}=2900$ and 6900, respectively, which are nearly equal to $L_{B}$ . Indeed, Chan-Braun et al. (Reference Chan-Braun, García-Villalba and Uhlmann2011) noted that the average operator over a tile of size $L_{B}$ filtered out turbulent structures much smaller or much larger than $L_{B}$ . Such a filter effect was also actually produced by the sphere surface in the rough bottom cases. However, in the limits of the cross-correlation analysis, the present results show that the mechanism of flow–roughness interaction both in the transitionally and fully rough regimes appears significantly different from the interaction between the turbulent flow and a smooth wall.