Direct numerical simulation of open-channel flow over a fully-rough wall at moderate relative submergence

TL;DR

This study uses direct numerical simulation to analyze open-channel flow over a fully-rough bed of spheres, revealing how flow structures and forces differ from smoother surfaces and how roughness influences turbulence and drag in the fully-rough regime.

Contribution

First direct numerical simulation of fully-rough open-channel flow over spheres, comparing flow structures and forces with experimental data and smooth-wall analogies.

Findings

Flow affected by Reynolds number and roughness geometry near the bed.

Roughness influences flow structures mainly above the virtual origin.

Hydrodynamic forces can be related to smooth-wall turbulence models.

Abstract

Direct numerical simulation of open-channel flow over a bed of spheres arranged in a regular pattern has been carried out at bulk Reynolds number and roughness Reynolds number (based on sphere diameter) of approximately 6900 and 120, respectively, for which the flow regime is fully-rough. The open-channel height was approximately 5.5 times the diameter of the spheres. Extending the results obtained by Chan-Braun et al. (J. Fluid Mech., vol. 684, 2011, 441) for an open-channel flow in the transitionally-rough regime, the present purpose is to show how the flow structure changes as the fully-rough regime is attained and, for the first time, to enable a direct comparison with experimental observations. The results indicate that, in the vicinity of the roughness elements, the average flow field is affected both by Reynolds number effects and by the geometrical features of the roughness, while at larger wall-distances this is not the case, and roughness concepts can be applied. The flow-roughness interaction occurs mostly in the region above the virtual origin of the velocity profile, and the effect of form-induced velocity fluctuations is maximum at the level of sphere crests. The spanwise length scale of turbulent velocity fluctuations in the vicinity of the sphere crests shows the same dependence on the distance from the wall as that observed over a smooth wall, and both vary with Reynolds number in a similar fashion. Moreover, the hydrodynamic force and torque experienced by the roughness elements are investigated. Finally, the possibility either to adopt an analogy between the hydrodynamic forces associated with the interaction of turbulent structures with a flat smooth wall or with the surface of the spheres is also discussed, distinguishing the skin-friction from the form-drag contributions both in the transitionally-rough and in the fully-rough regimes.

Figures

b), the fact that the Reynolds number in the present DNS is substantially lower should not matter. Figure [25](#S3.F25) shows the temporal cross-correlation of drag and lift force fluctuations, denoted by $R_{F^{\prime}_{x}F^{\prime}_{y}}(\Delta r_{x}=i\,L_{B},\Delta r_{z}=j\,L_{B},\Delta_{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](#S3.SS2), the lift force is associated with pressure fluctuations

b), the fact that the Reynolds number in the present DNS is substantially lower should not matter. Figure [25](#S3.F25) shows the temporal cross-correlation of drag and lift force fluctuations, denoted by $R_{F^{\prime}_{x}F^{\prime}_{y}}(\Delta r_{x}=i\,L_{B},\Delta r_{z}=j\,L_{B},\Delta_{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](#S3.SS2), the lift force is associated with pressure fluctuations

b), the fact that the Reynolds number in the present DNS is substantially lower should not matter. Figure [25](#S3.F25) shows the temporal cross-correlation of drag and lift force fluctuations, denoted by $R_{F^{\prime}_{x}F^{\prime}_{y}}(\Delta r_{x}=i\,L_{B},\Delta r_{z}=j\,L_{B},\Delta_{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](#S3.SS2), the lift force is associated with pressure fluctuations

Tables

Table 1. Table 1: Parameters of the simulations. N x subscript 𝑁 𝑥 N_{x} , N y subscript 𝑁 𝑦 N_{y} and N z subscript 𝑁 𝑧 N_{z} denote the number of grid points in the stramwise, wall-normal and spanwise directions, respectively. T o ​ b ​ s subscript 𝑇 𝑜 𝑏 𝑠 T_{obs} denotes the simulation time (excluding the initial transient), normalised by the bulk time-unit T b ​ H = H / U b ​ H subscript 𝑇 𝑏 𝐻 𝐻 subscript 𝑈 𝑏 𝐻 T_{bH}=H/U_{bH} . The data of case D50 is from Chan-Braun et al., ( 2011 ) .

$\mathrm{run}$	$\frac{U_{bh}}{u_{\tau }}$	$\frac{H}{D}$	$R{e}_{bH}$	$R{e}_{\tau }$	$D^{+}$	${\Delta }_x^{+}$	$\frac{L_x}{H}$	$\frac{L_y}{H}$	$\frac{L_z}{H}$	$\frac{D}{{\Delta }_x}$	$N_x$	$N_y$	$N_z$	$\frac{T_{obs}}{T_{bH}}$	source
D50	$12.2$	$5.6$	$2872$	$234$	$49$	$1.1$	$12$	$1$	$3$	$56$	$3072$	$256$	$768$	120	CB2011
D120	$12.4$	$5.4$	$6865$	$544$	$119$	$1.1$	$12$	$1$	$3$	$106$	$6912$	$576$	$1728$	60	present

Table 2. Table 2: Values of the parameters, of the integration constant C I ​ I + superscript subscript 𝐶 𝐼 𝐼 C_{II}^{+} and of the shift of the velocity profile Δ ​ ⟨ u ¯ ⟩ + Δ superscript delimited-⟨⟩ ¯ 𝑢 \Delta\left\langle\overline{u}\right\rangle^{+} (computed at y + = y 0 + + D + superscript 𝑦 superscript subscript 𝑦 0 superscript 𝐷 y^{+}=y_{0}^{+}+D^{+} ) for the present simulation D120, for the simulation D50 (Chan-Braun et al.,, 2011 ) , and for six of the experiments carried out by Amir et al., ( 2014 ) (original numbering). Note that Amir et al., ( 2014 ) chose y 0 = D subscript 𝑦 0 𝐷 y_{0}=D as opposed to the present choice y 0 = 0.85 ​ D subscript 𝑦 0 0.85 𝐷 y_{0}=0.85D .

$\mathrm{run}$	$\frac{H}{D}$	$R{e}_{bH}$	$D^{+}$	$\frac{1}{\kappa }\ln \left(D^{+}\right)$	$C_{II}^{+}$	$\Delta {⟨\overline{u}⟩}^{+}$	arrangement	source
D50	$5.6$	$2872$	$49$	$9.5$	$11.0$	$4.0$	square	CB2011
D120	$5.4$	$6865$	$119$	$11.6$	$10.5$	$6.3$	square	present
no. 1	$2.5$	$4915$	$170$	$12.5$	$11.1$	$6.5$	hexagonal	AM2014
no. 4	$5.0$	$15508$	$243$	$13.4$	$11.1$	$7.4$	hexagonal	AM2014
no. 2	$2.5$	$10678$	$376$	$14.5$	$10.9$	$8.7$	hexagonal	AM2014
no. 5	$5.0$	$32373$	$526$	$15.3$	$10.6$	$9.8$	hexagonal	AM2014
no. 8	$7.5$	$63305$	$658$	$15.8$	$10.2$	$10.7$	hexagonal	AM2014
no. 6	$5.0$	$45763$	$724$	$16.1$	$10.4$	$10.8$	hexagonal	AM2014

Table 3. Table 3: Statistics of the hydrodynamic force acting on the top-layer particles in runs D50 and D120. F i / F R subscript 𝐹 𝑖 subscript 𝐹 𝑅 F_{i}/F_{R} denotes the normalised mean force component in the x i subscript 𝑥 𝑖 x_{i} -direction, α F = arctan ⁡ ( 𝒞 F y / 𝒞 F x ) subscript 𝛼 𝐹 subscript 𝒞 subscript 𝐹 𝑦 subscript 𝒞 subscript 𝐹 𝑥 \alpha_{F}=\arctan{(\mathcal{C}_{F_{y}}/\mathcal{C}_{F_{x}})} denotes the angle of the resulting force with respect to the x 𝑥 x -axis, σ F i subscript 𝜎 subscript 𝐹 𝑖 \sigma_{F_{i}} , 𝒮 F i subscript 𝒮 subscript 𝐹 𝑖 \mathcal{S}_{F_{i}} and 𝒦 F i subscript 𝒦 subscript 𝐹 𝑖 \mathcal{K}_{F_{i}} are the standard deviation, the skewness and the kurtosis of the i 𝑖 i -th component of the force, respectively. The reference force is defined as F R = ϱ ​ u τ 2 ​ L B 2 subscript 𝐹 𝑅 italic-ϱ subscript superscript 𝑢 2 𝜏 superscript subscript 𝐿 𝐵 2 F_{R}=\varrho u^{2}_{\tau}L_{B}^{2} .

Run	$F_x/F_R$	$F_y/F_R$	$F_z/F_R$	${\sigma }_{F_x}/F_R$	${\sigma }_{F_y}/F_R$	${\sigma }_{F_z}/F_R$	${\alpha }_F$	$F_x^{\left(\nu \right)}/F_x$	$F_y^{\left(\nu \right)}/F_y$
D50	$1.14$	$0.37$	$0.0002$	$1.31$	$0.65$	$1.26$	${18.1}^{\circ }$	$0.55$	$0.06$
D120	$1.17$	$0.48$	$0.002$	$1.57$	$0.76$	$1.40$	${22.2}^{\circ }$	$0.42$	$0.02$
	${𝒮}_{F_x}$	${𝒮}_{F_y}$	${𝒮}_{F_z}$	${𝒦}_{F_x}$	${𝒦}_{F_y}$	${𝒦}_{F_z}$
D50	$0.056$	$0.27$	-$0.014$	$5.07$	$5.76$	$4.32$
D120	$0.09$	$0.13$	-$0.012$	$3.64$	$4.17$	$3.91$

Table 4. Table 4: Statistics of torque acting on the roughness elements for cases D50 and D120. T i / F R subscript 𝑇 𝑖 subscript 𝐹 𝑅 T_{i}/F_{R} denotes the normalised mean torque component in the x i subscript 𝑥 𝑖 x_{i} -direction, σ T i subscript 𝜎 subscript 𝑇 𝑖 \sigma_{T_{i}} , 𝒮 F i subscript 𝒮 subscript 𝐹 𝑖 \mathcal{S}_{F_{i}} and 𝒦 T i subscript 𝒦 subscript 𝑇 𝑖 \mathcal{K}_{T_{i}} are the standard deviation, the skewness and the kurtosis of the i 𝑖 i -th component of the torque, respectively.

Run	$T_x/T_R$	$T_y/T_R$	$T_z/T_R$	${\sigma }_{T_x}/T_R$	${\sigma }_{T_y}/T_R$	${\sigma }_{T_z}/T_R$
F50	-$0.0006$	-$0.00043$	-$0.72$	$0.17$	$0.11$	$0.27$
D120	-$0.0011$	$0.0005$	-$0.44$	$0.10$	$0.10$	$0.11$
	${𝒮}_{T_x}$	${𝒮}_{T_y}$	${𝒮}_{T_z}$	${𝒦}_{T_x}$	${𝒦}_{T_y}$	${𝒦}_{T_z}$
F50	-$0.010$	-$0.0082$	-$0.75$	$3.78$	$4.92$	$3.35$
D120	$0.009$	-$0.001$	-$0.52$	$3.48$	$3.51$	$3.38$

Equations

$${⟨\overline{u}⟩}^{+}(y^{+})=\frac{1}{\kappa }\ln (y^{+}-y_0^{+})+C_I^{+}(k^{+})$$

$${⟨\overline{u}⟩}^{+}(y)=\frac{1}{\kappa }\ln \frac{y-y_0}{k}+C_{II}^{+}.$$

$${\mathbf{u}}^{\prime }$$

$${\mathbf{u}}^{\prime \prime }$$

$${\mathbf{u}}^{\prime \prime \prime }$$

$$⟨\overline{\left({\varphi }^{\prime \prime }-{\varphi }^{\prime \prime \prime }\right)\left({\psi }^{\prime \prime }-{\psi }^{\prime \prime \prime }\right)}⟩=⟨\overline{{⟨\overline{\varphi }⟩}_B{⟨\overline{\psi }⟩}_B}⟩-⟨\overline{\varphi }⟩⟨\overline{\psi }⟩,$$

$$\mathbf{\tau }={\mathbf{\tau }}_{\nu }-p_{tot}\mathbf{I},$$

$$p_{tot}(\mathbf{x},t)=p_l(x,t)+p(\mathbf{x},t),$$

$$\Delta {⟨\overline{u}⟩}^{+}=C_I^{+}(D^{+}\to 0)-{\delta }_{sub}^{+}+\frac{1}{\kappa }\ln {\delta }_{sub}^{+}$$

$$C_{II}^{+}={⟨\overline{u}⟩}^{+}\left(D+y_0\right).$$

$$\Delta {⟨\overline{u}⟩}^{+}=C_I^{+}(D^{+}\to 0)-C_{II}^{+}+\frac{1}{\kappa }\ln (D^{+})$$

$$C_{II}^{+}\le {\delta }_{sub}^{+}-\frac{1}{\kappa }\ln {\delta }_{sub}^{+}+\frac{1}{\kappa }\ln (D^{+}).$$

$$C_{II}^{+}\lesssim \frac{1}{\kappa }\ln (D^{+})$$

$${\tilde{\mathbf{\tau }}}_n={\widetilde{⟨\overline{\mathbf{\tau }}⟩}}_B\cdot \tilde{\mathbf{n}}.$$

$$\mathbf{F}={\int }_{\tilde{\mathcal{S}}}{\tilde{\mathbf{\tau }}}_{n}dS$$

$${\mathbf{F}}^{(i,j)}={\int }_{{\mathcal{S}}^{(i,j)}}\overline{\mathbf{\tau }}\cdot {\mathbf{n}}^{(i,j)}dS,$$

$$\mathbf{T}={\int }_{\tilde{\mathcal{S}}}\widetilde{{\mathbf{r}}_c}\times {\tilde{\mathbf{\tau }}}_{n}dS$$

$$\mathbf{f}(y)=\frac{D}{2}{\int }_0^y{\int }_0^{2\pi }{\tilde{\mathbf{\tau }}}_n(\eta ,\theta )d\theta d\eta ,$$

$$\mathcal{T}={\int }_{-s/2}^{s/2}{\int }_{-s/2}^{s/2}{\tilde{\mathbf{r}}}_s\times \left({\tilde{⟨\overline{{\mathbf{\tau }}_{\nu }}⟩}}_B{|}_{y=0}\cdot {\mathbf{e}}_y\right)dxdz$$

$${\mathcal{T}}^{(p)}={\int }_{-s/2}^{s/2}{\int }_{-s/2}^{s/2}{\tilde{\mathbf{r}}}_s\times \left({\widetilde{⟨\overline{\mathbf{\tau }}⟩}}_B{|}_{y=0}\cdot {\mathbf{e}}_y\right)dxdz$$

$$E_1$$

$$E_2$$

$$E_3$$

$$\overline{q}(\mathbf{X})=\frac{1}{\mathcal{T}}{\int }_{\mathcal{T}}qdt$$

$${\widetilde{⟨q⟩}}_B(\tilde{\mathbf{x}})$$

$$\phi\left(\boldsymbol{X}\right)=\left\{\begin{array}[]{l l}1&\mathrm{if}\quad\left|\boldsymbol{X}-\boldsymbol{x}^{(j)}_{c}\right|>D/2\quad\mathrm{for}\>\mathrm{every}\quad j=1,\ldots,N_{s}\\ 0&\mathrm{otherwise}\>\>.\end{array}\right.$$

$${⟨q⟩}_B(\mathbf{x})$$

$$\mathrm{or}$$

$${\widetilde{⟨q⟩}}_B(\tilde{\mathbf{x}})$$

Full text

Different statistical tools were used to investigate the flow field in the roughness sublayer and in the logarithmic region.

The results indicate that, in the vicinity of the roughness elements, the average flow field is affected both by Reynolds number effects and by the geometrical features of the roughness,

while at larger wall-distances this is not the case, and roughness concepts can be applied. Thus, the roughness function is computed which in the present set-up can be expected to depend on the relative submergence.

The flow-roughness interaction occurs mostly in the region above the virtual origin of the velocity profile, and the effect of form-induced velocity fluctuations is maximum at the level of sphere crests. In particular,

the root mean square of fluctuations about the streamwise component of the average velocity field

reflects the geometry of the spheres in the roughness sublayer and attains a maximum value just above the roughness elements. The latter is significantly weakened and shifted towards larger wall-distances as compared to the transitionally-rough regime or the case of a smooth wall. The spanwise length scale of turbulent velocity fluctuations in the vicinity of the sphere crests shows the same dependence on the distance from the wall as that observed over a smooth wall, and both vary with Reynolds number in a similar fashion. Moreover, the hydrodynamic force and torque experienced by the roughness elements are investigated and the footprint left by vortex structures on the stress acting on the sphere surface is observed. Finally, the possibility either to adopt an analogy between the hydrodynamic forces associated with the interaction of turbulent structures with a flat smooth wall or with the surface of the spheres is also discussed, distinguishing the skin-friction from the form-drag contributions both in the transitionally-rough and in the fully-rough regimes.

1   Introduction

The present investigation is motivated by the problem of erosion and deposition of sediment at the bottom of fluvial and estuarine environments as well as along hill slopes, which result from the action of the surface water flow. The ultimate goal is to understand the morphological evolution of the sediment bed. Other applications of channel flows can be found in chemistry, biological fluid dynamics and industrial engineering.

Providing reliable predictions of the river-bed evolution requires a clear picture of the interaction between the flow and the bottom roughness due to sediments and bedforms. Typically, in the regimes of practical interest, the flow is turbulent and the bottom is not smooth, but it is characterised by the presence of natural or artificial protrusions that affect the structure of the flow over a region few times thicker than the size of the protrusions, namely the roughness sublayer. Hence, a detailed description of turbulence in the vicinity of individual roughness elements is necessary to comprehend the dynamics of solid-fluid interaction and possibly formulate consistent relationships with quantities observable at larger scales. This knowledge becomes crucial in the case of shallow open-channels which are characterised by small values of the relative submergence (also termed inundation ratio) $H/k$, where $H$ denotes the open-channel height and $k$ the roughness size, since velocity fluctuations originating at the bottom may significantly affect the entire flow field. However, even in the simple case in which monosized regular roughness elements are considered, obtaining accurate measurements of the velocity and pressure fluctuations in the crevices between the roughness elements is extremely difficult (e.g. Hong et al.,, 2011; Amir et al.,, 2014), and discrepant pictures of the origin of the turbulent vortices and of their influence on the flow structure have been provided (e.g. Marusic et al.,, 2010, for a review).

It is well established that an open-channel flow (or similarly a boundary layer) over a smooth wall at sufficiently high Reynolds numbers develops a viscous sublayer, dominated by viscous effects, and a logarithmic region, where the fluid viscosity plays a negligible role, matched together through a buffer layer characterised by strong normal Reynolds stresses. Let us consider roughness elements located on a flat smooth wall with a certain arrangement,

defining the geometrical roughness size $k$ as the average distance from the wall to the crest of the roughness elements.

Let, for the moment, large values of $H/k$ be considered,

for which effects associated with the free-slip boundary condition at the free surface of the open channel can be neglected in the vicinity of the bottom wall. Therefore, for increasing values of the roughness Reynolds number $k^{+}=ku_{\tau}/\nu$, where $u_{\tau}$ and $\nu$ denote the friction velocity and the kinematic viscosity of the fluid, three flow regimes can be identified: the hydraulically-smooth regime, the transitionally-rough regime and the fully-rough regime (Jiménez,, 2004).

Let henceforth $y$ and $y_{0}$ denote the wall-normal coordinate and the position of the virtual origin; the latter one is defined as the plane where a smooth wall should be placed (in absence of the roughness elements) to observe the logarithmic region originating at the same distance from $y_{0}$ as in the rough-wall case.

In the hydraulically-smooth regime,

the roughness elements are entirely contained within the viscous sublayer and the velocity profile, as a function of the distance from the virtual wall,

practically collapses upon the profile that could be obtained over a smooth wall at the same bulk Reynolds number

$Re_{bh}=hU_{bh}/\nu$, where $h$ equals $H-y_{0}$ and

$U_{bh}$

denotes the bulk velocity defined as $U_{bh}=1/h\int^{H}_{y_{0}}\left\langle\overline{u}\right\rangle\,dy$, the operator $\left\langle\overline{\cdot}\right\rangle$ indicating the statistical average defined more precisely below. Then, by gradually increasing the Reynolds number $k^{+}$ until the transitionally-rough regime is attained, the viscous sublayer is significantly thinned with respect to the hydraulically-smooth regime, while the mean velocity, normalised by $u_{\tau}$, is reduced (shifted) in the logarithmic region as an effect of the increasing momentum transfer to the roughness elements. Finally, for $k^{+}\gtrsim 55-90$, the buffer layer disappears and the fully-rough flow regime is reached (Ligrani and Moffat,, 1986).

The shift of the velocity profile in the logarithmic region with respect to that observed in absence of the roughness elements (i.e. over a smooth-wall), namely the roughness function $\Delta\left\langle\overline{u}\right\rangle^{+}$, can be used as a universal parameter to classify the flow regime in wall bounded turbulent flows characterised by large values of the relative submergence so that the flow structure becomes independent of $H/k$ ($H/k>40$, Jiménez,, 2004). Instead, for small values of the relative submergence, some mechanisms of the wall-turbulence are possibly affected while the effects of the geometrical features of the roughness elements can be recognised in the entire domain. These cases can be suitably studied as flows over obstacles, since roughness concept cannot be applied to any region of the flow field. However, for moderate relative submergence, the effects associated with individual roughness elements tend to vanish far from the bottom and a logarithmic region can be clearly distinguished over the roughness sublayer (e.g. Bayazit,, 1976). In the latter case, which is the object of the present investigation, the hydraulically-smooth, transitionally- and fully-rough regimes can be identified on the basis of the roughness function and of the relative submergence $H/k$ which becomes a parameter of the problem.

The roughness characteristics of a natural bed (for instance a river channel) are not homogeneous (due to the non-regular shape and arrangement of roughness elements and to the presence of multiple scales defining the roughness geometry) and can vary with time (e.g. due to sediment transport and bedform evolution). We presently consider the particular case of fixed and identical roughness elements arranged with a regular pattern, which limits the geometrical characterization of the roughness to that of a minimal roughness unit (i.e. a single roughness elements and its closest neighbours) and allows us to consider the roughness characteristics as spatially homogeneous and constant with time.

Although this approach limits the scope of application of the present results to specific problems, it allows us to identify clearly the mechanism of flow-roughness interaction.

Let us focus our attention on the flow structure in the

logarithmic region, where

the velocity profile behaves like a logarithmic similarity law (i.e. Prandtl’s celebrated “law of the wall”) which can be expressed in terms of wall units as follows:

[TABLE]

or equivalently in terms of $k$:

[TABLE]

It turns out that $\frac{d\left\langle\overline{u}\right\rangle^{+}}{dy}$ only depends on $y$, $k$ and $u_{\tau}$, while $(y^{+}-y^{+}_{0})\frac{\partial\left\langle\overline{u}\right\rangle^{+}}{\partial y^{+}}$ is constant, namely the inverse of the Von Kármán constant $\kappa$, $y_{0}$ denoting the aforementioned virtual origin of the wall-normal coordinate. The integration constant $C^{+}_{I}$ was experimentally found to tend to $5.1$ for $k^{+}\rightarrow 0$ (hydraulically-smooth regime) and to $C^{+}_{II}-\frac{1}{\kappa}\ln k^{+}$, with $C^{+}_{II}$ constant, for $k^{+}\rightarrow\infty$ (fully-rough regime) (Nikuradse,, 1933; Schlichting,, 1968; Pimenta et al.,, 1975; Ligrani and Moffat,, 1986). The value of the constant $C^{+}_{II}$, in the fully-rough regime

at large relative submergence,

depends on the shape and the arrangement of roughness elements and tends to $8.5$ for the case of a boundary layer over a plane, well-packed, sandy bottom (Ligrani and Moffat,, 1986). The fact that $\left\langle\overline{u}\right\rangle^{+}$ is independent of $k^{+}$ in the fully-rough regime follows the disappearance of the buffer layer and the attenuation of viscous effects above the roughness (Tani,, 1987). This feature is in line with Townsend’s hypothesis, according to which in a boundary layer at high values of $k^{+}$ the turbulence structure above the roughness sublayer is practically unaffected by the roughness shape and arrangement. From equation (2) it can be observed that, once the fully-rough regime is attained,

if the relative submergence is sufficiently large that a logarithmic region can be identified,

the velocity profile experiences a shift which increases proportionally to the logarithm of $k$ (e.g. Nikuradse,, 1933; Ligrani and Moffat,, 1986). This suggests that the flow structure may remain basically unaffected by further increases of $k^{+}$, and that increases of $k^{+}$ result in the progressive truncation of the velocity profile at distance approximately $k^{+}$ from the wall.

Deviations of such a distance from $k^{+}$ are associated with the position of $y_{0}$ that in turn reflects the geometrical features of the roughness.

After the pioneering systematic work of Nikuradse, (1933) on rough pipe flows, Schlichting, (1936) was the first who experimentally investigated the effects on open-channel flow of combining different (regular, homogeneous) arrangements of spheres mounted on a smooth wall and varying their size. Since then, large efforts have been devoted to investigate the relevance of the roughness geometrical features on the turbulence structure.

Providing an exhaustive roughness characterization is still one of the major challenges in the field.

A number of recent experimental studies showed that the size, shape and arrangement of roughness elements can significantly affect the bottom drag and the flow structure in the near bottom region

(Amir and Castro,, 2011; Cooper et al.,, 2013; Florens et al.,, 2013; Willingham et al.,, 2014; Placidi and Ganapathisubramani,, 2015; Bossuyt et al.,, 2017, to cite a few examples)

In the attempt of synthesizing the complexity of the roughness geometry within a single parameter, the sand-grain roughness (or effective roughness), $k_{s}$, is commonly used which is either found to be proportional to the roughness size $k$ or proportional to the channel height (boundary layer thickness), depending on the precise geometrical features of the rough surface.

In the fully rough regime, for large values of $H/k$, $k_{s}$ is defined as the roughness height which produces the same roughness function as that measured by Nikuradse, (1933), cf. Flack and Schultz, (2010). Concerning the value of $k_{s}$, Schlichting, (1936) observed that, for mono-sized, spherical roughness elements in a hexagonal arrangement, $k/k_{s}$ ranged from $0.26$ to $4.41$ only by varying the distance between the grains. Using a single parameter to characterise rough surfaces is undoubtedly convenient as long as the relative submergence and the roughness Reynolds number are sufficiently high ($H/k>40$ and $k_{s}^{+}>50$, Jiménez, (2004)). However, a single parameter does not suffice, at moderate relative submergence, to characterize the roughness function, which ultimately depends also on $H/k$ and on the other length scales characterising the roughness geometry.

At moderate relative submergence, the value of the roughness function at which the fully-rough regime is attained can be different from that measured by Nikuradse, (1933) and, consequently, the parameter $k_{s}^{+}$ can be no longer used unambiguously to determine the flow regime. For instance, Amir et al., (2014) have recently carried out experiments of a moderately shallow open-channel flow in the fully-rough regime. They could observe the mean velocity profile to follow a logarithmic curve in the core of the flow field and, in one of their tests, they measured the roughness function equal to $6.5$, which corresponds to a value of $k_{s}^{+}$ barely larger than $50$, although the flow regime was fully-rough.

Indeed, the quantity $k_{s}$ is defined heuristically, leaving us free to interpret it as the hydrodynamic response of the flow to the disturbance induced by the roughness. For example, in this line, Orlandi et al., (2003) and Flores and Jimenez, (2006) investigated the effect of superimposing a disturbance of the velocity field in the vicinity of a smooth-wall over an otherwise undisturbed channel-flow and observed the development of turbulent fluctuations associated with the disturbance that were similar to those induced by a physical roughness.

In fact, the velocity close to the wall could be locally and instantaneously nullified by the disturbance.

Although the bottom was not rough, Orlandi et al., (2003) and Flores and Jimenez, (2006) simulated the fully-rough regime and, in principle, they could have estimated the value of

the roughness function and of $k_{s}$.

In the present

direct numerical

simulations,

the effect of the relative submergence cannot be neglected, and

the grain size $k$ will be used as the length scale instead of $k_{s}$,

since the results cannot be generalised to other rough-wall flows characterised by the same roughness function.

Nonetheless, Chan-Braun et al., (2011), who performed two direct numerical simulations of moderately shallow open-channel flow in the hydraulically-smooth and the transitionally-rough regimes, noted that the values of the roughness function were in the range of those obtained for rough boundary layers at the same roughness Reynolds numbers $k_{s}^{+}$. Chan-Braun et al., (2011, 2013) could observe the presence of a buffer layer just above the crest of the roughness elements in which the velocity field is affected by viscosity, the friction velocity, the channel height and bulk velocity as well as the roughness size. It is widely recognised that the features of the turbulence structure related to the geometrical characteristics of the roughness are lost at a certain distance from the wall in well developed boundary layers at sufficiently high Reynolds numbers. Experimental results confirm a fair agreement with the “law of the wall” (e.g Ligrani and Moffat,, 1986; Bandyopadhyay,, 1987; Schultz and Flack,, 2007). However, some of the effects associated with the roughness or the presence of a pressure gradient were not found completely to disappear in experiments carried out in plane channel flow (e.g. Grass et al.,, 1991; Hong et al.,, 2011, 2012) and open-channel flow (e.g. Balachandar and Ramachandran,, 1999; Tachie et al.,, 2000; Amir et al.,, 2014).

In particular, Hong et al., (2011) found that Townsend’s hypothesis for the Reynolds stress statistics was supported above the roughness sublayer ($y\gtrsim 2k$) but the presence of roughness-related small-scale turbulence affected the dissipation rate in the entire flow field, while George, (2007) showed that the effect of the mean pressure gradient on turbulence statistics in pressure-gradient-driven channel (or pipe) flows can be assumed negligible only in a certain region in the vicinity of the bottom and at moderate Reynolds numbers. Such evidences has questioned the universality of Townsend’s similarity hypothesis, thereby challenging researchers to define its limitations more clearly. Even though this lies outside the purpose of the present contribution, it is worthwhile to mention that the influence of the roughness on the turbulence structure could be presumably amplified if the size of the roughness elements is of the same order of magnitude as the open-channel height, without necessarily degenerating into the flow around a sequence of obstacles. For instance, Amir and Castro, (2011) observed that, in a boundary layer over genuinely three-dimensional roughness, inner and outer scales were distinguishable and separated as long as the size of the roughness did not exceed $15\%$ of the boundary layer thickness, which is well above the value ($2.5\%$) indicated by Jiménez, (2004).

From these considerations it follows that an investigation of the flow-roughness interaction at the scale of the roughness elements is needed in order to push significantly further our understanding of the dynamics of moderately shallow open-channel flow. Indeed, an exhaustive description of the flow-roughness interactions, in particular between the transitionally-rough and the fully-rough regimes, is still missing and should in principle require the systematic exploration of different wall configurations and relative submergence, as well as the possibility to make accurate measurements of velocity and pressure fields in the vicinity of the roughness elements.

An effect of roughness, in the fully-rough regime, is to introduce and sustain turbulent fluctuations characterised by length and time scales of the same order of magnitude as $k$ and $k/U_{bh}$, respectively. In fact, in a series of experiments of open-channel flow in the fully-rough regime in which a “random” distribution of closely-packed spheres resting on the wall was used, Grass et al., (1991) found that the fluctuations of the streamwise velocity component in the vicinity of sphere crests were spanwise-correlated over a specific wavelength, similarly to what happens over a smooth wall.

However, in that case, the wavelength was proportional to the size of the roughness elements.

In the experiment of Grass et al., (1991) the arrangement of the spheres on the wall was not regular, thus the measurements of velocity fluctuations in the vicinity of sphere crests could be affected also by the local geometrical configuration of the rough wall. A more general result could be obtained by systematically computing the fluctuations about the time-averaged flow field

at many points distributed in the vicinity of the roughness elements.

This task could be more easily addressed with a numerical approach, as Grass et al., (1991) themselves suggested.

Amir et al., (2014) investigated the hydrodynamic force exerted on spheres mounted in hexagonal pattern on the wall, by instrumenting several spheres with four-point pressure probes.

They computed the statistics of pressure fluctuations for several cases differing in the bed slope and flow depth and were able to estimate the integral time scale of drag and lift fluctuations on the basis of measurements made from spheres equispaced in a row. However, by adopting the four-point pressure technique, Amir et al., (2014) did not consider the distribution of the stress on the surface of individual spheres, and they could not account for the distribution

of the shear stress and how it is affected by Reynolds number effects.

Also this task can be more easily achieved numerically.

In the present work the study of Chan-Braun et al., (2011, 2013) has been extended to higher Reynolds numbers such that the fully-rough regime is attained. We have kept the geometrical particle arrangement identical as in the earlier study, i.e. a square pattern with a relative submergence of $H/k\simeq 5.5$. By increasing the bulk Reynolds number, the diameter of the spherical roughness elements in the present simulation well exceeds the value of $100$ wall units.

The organization of the paper is as follows. After describing the chosen setup and the numerical approach, an analysis of the flow structure below and over the crest of the spheres is provided in section 3.1. A comparison between length- and time-scales of the turbulent flow over the roughness elements and that obtained over a flat smooth wall for the same value of the bulk Reynolds number was also made. Moreover, the distribution of the stress on the surface of the roughness elements is analysed in section 3.2.

Finally, although the Reynolds numbers investigated by Amir et al., (2014) are larger than those reproducible nowadays by DNS, a comparison between their experimental results and present numerical results in the fully-rough regime was possible and is shown in section 3.3. The paper closes with the description of a conceptual model of the interaction between turbulent structures and the roughness elements, and, in section 4, with the conclusions of the present results.

2   Flow configuration and numerical approach

Direct numerical simulations of the incompressible Navier-Stokes equation were performed over a rectangular domain of dimensions $L_{x}$, $L_{y}$ and $L_{z}$ in the streamwise, wall-normal and spanwise directions, respectively. Let us indicate hereafter the simulations of the transitionally- and fully-rough open-channel flow with D50 and D120, respectively. The flow configuration is sketched in figure 1, details of the roughness geometry are shown in figure 2. The arrangement of roughness elements is the same as that adopted by Chan-Braun et al., (2011) for their simulations, hence it consists in one layer of $1024$ monosized rigid spheres crystallised on the wall at the vertices of a square grid of side $L_{B}=D+\Delta_{B}$ aligned with the $x$- and $z$-directions (yellow spheres in figure 2), and a second layer of spherical caps (highlighted in red in figure 2) with the same arrangement as the first one shifted by $(L_{B}/2,y_{2},L_{B}/2)$, where $y_{2}=D/2-\sqrt{2}(D/2+\Delta_{x})$

and $\Delta_{x}$ denotes the computational grid spacing.

Figure 2 also shows that, by exploiting the symmetry properties of the roughness arrangement, a cuboidal sub-domain $\mathcal{B}=[-L_{B}/2,L_{B}/2[\times[0,H]\times[-L_{B}/2,L_{B}/2[$ can be defined in the local coordinates system $(\widetilde{x},y,\widetilde{z})$ with the same orientation as $(x,y,z)$ and origin in the projection on the wall of the top-layer sphere center. The sub-domain $\mathcal{B}$ is geometrically periodic in the streamwise and spanwise directions and invariant to the exchange between $\widetilde{x}$- and $\widetilde{z}$-axis, i.e. a square arrangement. There are no gaps between the spheres and the wall, while the minimal distance $\Delta_{B}$ is required between the spheres by the immersed boundary method proposed by Uhlmann, (2005). Since roughness elements are spherical, hereinafter the roughness size and the relative Reynolds number will be referred to as $D$ and $D^{+}$, instead of $k$ and $k^{+}$, respectively.

The numerical approach used by Chan-Braun et al., (2011) was also used for the present simulations. It consists in a second-order accurate fractional-step method. In particular, a semi-implicit scheme is employed for the viscous terms along with a three-step (low storage) Runge-Kutta method for the nonlinear terms. Standard centred second-order finite-difference approximations of the spatial derivatives are used over a staggered uniform Cartesian grid of spacing $\Delta_{x}=\Delta_{y}=\Delta_{z}$. The spherical roughness elements are represented by means of an immersed boundary technique, and further details on the numerical method can be found in Uhlmann, (2005) and Uhlmann, (2006). This numerical procedure underwent an extensive validation and has been recently used to perform the direct numerical simulation of flows in which particles were either fixed or free to move (Uhlmann,, 2008; Chan-Braun et al.,, 2011; Garcia-Villalba et al.,, 2012; Uhlmann and Doychev,, 2014; Kidanemariam and Uhlmann, 2014a, ; Kidanemariam and Uhlmann, 2014b, ; Chouippe and Uhlmann,, 2015; Mazzuoli et al.,, 2016; Uhlmann and Chouippe,, 2017). The spacing $\Delta_{B}$ was set equal to $2\Delta_{x}$. We do not expect the results to be significantly affected by the value of $\Delta_{B}$, as along as it remains a small fraction of the particle diameter $D$ and, consequently, it does not entail any significant change in the width (and spanwise separation) of the inter-particle grooves.

The size of the computational domain, the small-scale resolution and further details are shown in table 1. The present parameter point corresponds to $D^{+}=119$ and $Re_{bH}=6865$, and it will henceforth be denoted as case D120. The table also shows the parameters of case D50 of Chan-Braun et al., (2011) which will be used frequently in the discussion below. The continuity and momentum equations were solved numerically, obtaining the velocity components and pressure, denoted by $(u,v,w)$ and $p$ respectively, throughout the whole computational domain including the space occupied by the spheres. Periodicity conditions were applied at the boundaries of the domain in the streamwise and spanwise directions, and the no-slip condition was forced at the fluid-spheres interface by means of the immersed boundary method proposed by Uhlmann, (2005), i.e. through force terms directly added to the momentum equations (direct forcing method). In the immersed boundary method, the velocity field at the end of a full time step slightly deviates from the desired value (zero) at the surface of the spheres due to the effect of the projection step. The magnitude of this error amounts to $6.6\cdot 10^{-3}\,u_{\tau}$ on average, with the local instantaneous maximum measuring $5.8\cdot 10^{-1}\,u_{\tau}$. We believe that this level of error is not significant with respect to the statistical analysis performed in the present manuscript. Finally, no-slip and free-slip boundary conditions were imposed to the flow field at the wall $(y=0)$ and at the open-surface $(y=H)$ of the computational domain, respectively.

The mass flow is maintained steady throughout the simulations by a uniform pressure gradient which drives the flow and is updated at each time step. Thus, while the bulk velocity defined as $U_{bH}=1/H\int^{H}_{0}\left\langle\overline{u}\right\rangle\,dy$ and consequently the bulk Reynolds number $Re_{bH}=U_{bH}H/\nu$ were constant, the Reynolds numbers $Re_{\tau}$ and $D^{+}$ fluctuated about the average value indicated in table 1.

The definition of the Reynolds number $R_{bh}=U_{bh}h/\nu$ based on the effective flow depth $h=H-y_{0}$ and on the bulk velocity $U_{bh}=1/h\int_{0}^{h}\left\langle\overline{u}\right\rangle\,dy$ is recalled here for the sake of clarity. The ratio between $U_{bh}$ and $U_{bH}$ was found equal to $1.16$ and $1.17$ for the simulations D50 and D120, respectively.

The grid spacing was sufficiently small to resolve the vortices associated with the dissipative turbulent scales ($\Delta^{+}_{x}=1.1$) while the size of the computational box was found large enough to include the large vortex structures of the flow. In fact, the pre-multiplied spectra of turbulent velocity fluctuations (figure omitted) show that only weak energy is associated with structures larger than the domain presently considered which, anyway, should not influence the spectra in the vicinity of the bottom (Del Alamo et al.,, 2004; Chan-Braun,, 2012).

Moreover, additional DNS of open-channel flow over a smooth wall have been performed for the same box-size and bulk Reynolds number as case D120 using a pseudo-spectral method (Kim et al.,, 1987). These smooth-wall data, as well as those reported for $Re_{bH}=2870$ by Chan-Braun et al., (2011), will be used for the purpose of comparison below.

3   Results and discussion

The main concern of the present investigation is how an increase of the bulk Reynolds number changes the flow structure such that the flow regime becomes fully-rough. Prior to defining the roughness sublayer for the present simulation and exploring the numerical results, some considerations are formulated on the basis of the geometrical configuration of the solid boundary, in order to enable appropriate statistical tools to investigate the fluid-roughness interaction.

Turbulence in the vicinity of a spherical roughness element, namely in the roughness sublayer, is clearly neither homogeneous nor isotropic and in principle turbulence statistics should be described as functions of three space dimensions. For the present square particle arrangement in a doubly-periodic computational domain statistics are invariant with respect to spanwise or streamwise shifts by integer multiples of the distance between two spheres. Furthermore, statistics are symmetric with respect to the $(x,y)$-plane through the center of any sphere. These features are exploited by the averaging operators used here, which are formally defined in appendix A.

What will be hereafter referred to as average flow field and indicated with $\left\langle\overline{\boldsymbol{u}}\right\rangle_{B}$, can be defined in the sphere-boxes, namely the cuboidal periodic subdomain $\mathcal{B}$ introduced in section 2 (figure 2), and estimated by combining the time-average operator and the sphere-box-average operator. Note that $\left\langle\overline{\boldsymbol{u}}\right\rangle_{B}$ is a three-dimensional quantity. Indeed, $\mathcal{O}(100)$ flow fields were obtained systematically collecting one snapshot every $0.6$ bulk time-units ($T_{b}=H/U_{bH}$). Thus, the average flow field $\left\langle\overline{\boldsymbol{u}}\right\rangle_{B}$ was obtained in post-processing phase over $\mathcal{O}(10^{5})$ samples. Additionally, the plane-averaged velocity and pressure fields,

as well as the variance of their fluctuations,

were computed and collected during run time with a frequency $100$ times larger than the sampling of the snapshots. Then, plane-averaged samples were also averaged over time in post-processing phase. The simulations were preliminarily run until turbulence was fully developed before starting the sampling procedure. Turbulent fluctuations around the time-averaged ($\overline{\boldsymbol{u}}$), sphere-box/time-averaged ($\left\langle\overline{\boldsymbol{u}}\right\rangle_{B}$) and plane/time-averaged ($\left\langle\overline{\boldsymbol{u}}\right\rangle$) velocity fields can be defined as follows

[TABLE]

respectively, each helping to interpret physical phenomena at different scales. In a similar way, the fluctuations of pressure and other scalar or vectorial quantities can be also defined. The definition (3) of time-fluctuations can be useful to investigate the evolution of turbulence structures of size much larger than $D$ (i.e. comparable with that of the computational domain). The definition (4) seems more appropriate to study the interaction of individual roughness elements with the flow, selecting the turbulence scales larger than the size of $\mathcal{B}$. In the general case of an arbitrary arrangement of monosized spheres, fluctuations (4) would not contain the information associated with the geometrical properties of individual roughness elements (individual form contribution) because the spheres are identical, but that related to their arrangement (collective form contribution). Since the present arrangement is regular, the collective form contribution is also not present in fluctuations (4), while it definitely affects the average flow field $\left\langle\overline{\boldsymbol{u}}\right\rangle_{B}$.

In fact, the quantity $\left\langle\overline{\boldsymbol{u}}\right\rangle_{B}-\left\langle\overline{\boldsymbol{u}}\right\rangle$, which is equal to $\boldsymbol{u}^{\prime\prime}-\boldsymbol{u}^{\prime\prime\prime}$, is also called spatial disturbance (or form disturbance) and depends on the size of the sub-domain $\mathcal{B}$ (Nikora et al.,, 2001). It can be verified that the stress related to the spatial disturbance of the generic velocity components $\phi$ and $\psi$ is equal to the dispersive stress (or form-induced stress) reduced by the product of the respective plane/time-averaged velocity, which reads

[TABLE]

where the dispersive stress is the first term on the right side of (6). The same definition of the spatial disturbance can be extended also to the fluctuations of pressure or any scalar quantity.

Hence,

turbulent fluctuations defined by (5) contain both the time-fluctuations and those associated with the flow pattern around the roughness elements.

At this point let us also define the stress tensor for future reference, viz.

[TABLE]

where $\boldsymbol{\tau}_{\nu}=\varrho\nu\left(\nabla\mathbf{u}+\left(\nabla\mathbf{u}\right)^{T}\right)$ is the viscous stress (with $\varrho$ the fluid density), and $\mathbf{I}$ is the identity tensor. The pressure field in fully-developed plane-channel flow can be written as follows

[TABLE]

where $p_{l}$ corresponds to the linear variation in the streamwise direction due to the imposed driving pressure gradient (with $\langle p_{l}\rangle_{x}=0$) and $p$ is the instantaneous fluctuation whose box-average is zero.

3.1 The velocity and vorticity fields

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_{\tau}=\sqrt{\tau_{w}/\varrho}$ was estimated, where the value of the wall shear stress $\tau_{w}$ was extrapolated down to the wall-normal distance $y=y_{0}$ using the linear profile of the wall-normal total shear stress

$\tau_{tot}=\rho\nu\frac{\langle\overline{u}\rangle}{dy}-\rho\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.,, 2011) . Figure 4a 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 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., (2011). Indeed, it was found that the profiles of $(y^{+}-y^{+}_{0})\frac{d\left\langle\overline{u}\right\rangle^{+}}{dy^{+}}$ (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 4b). 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 to mention that other methodologies adopted to estimate the position of $y_{0}$, like the location of the drag force centroid (Jackson,, 1981), provided approximately to the same value (see also Appendix A1 of Chan-Braun et al.,, 2011).

Since, in the logarithmic region, the curves of figure 4b should be independent of $y^{+}$, it was also possible to estimate the value for the von Kármán constant, $\kappa$, as the inverse of the value attained at the local minimum of the curves, as suggested (among others) by Balachandar and Ramachandran, (1999), Tachie et al., (2000) and George, (2007). Hoyas and Jiménez, (2006) showed that the fact not to observe a wide region where $1/\kappa$ 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 $\kappa$ from the value $0.41$, estimated by Coles, (1968) for boundary layers, were often reported in the literature (e.g. Österlund et al.,, 2000; Marusic et al.,, 2013) and, in particular for the present channel-flow configuration, can be related to pressure gradient effects (George,, 2007).

However, it should be noted that the present values of $\kappa$ are in the range estimated by Marusic et al., (2013) ($\kappa=0.39\pm 0.02$), while a local minimum of $(y^{+}-y^{+}_{0})\frac{d\left\langle\overline{u}\right\rangle^{+}}{dy^{+}}$ is attained approximately at $y^{+}-y^{+}_{0}=50$ consistently with other numerical (Hoyas and Jiménez,, 2006) and experimental (e.g. Mckeon et al.,, 2004; Tachie et al.,, 2000) results.

Some considerations can be formulated for the simulation D50 on the basis of the velocity profiles of figure 4a. Chan-Braun et al., (2011) found that the flow regime for the simulation D50 was transitionally-rough. This also appears from the velocity profile (black-square line in figure 4a) which shows the presence of the buffer layer above the crest of the spheres. The thickness of the viscous sub-layer $\delta^{+}_{sub}$ was estimated from the intersection of the logarithmic law and the linear law of the wall (cf. figure 4a), which yields $\delta^{+}_{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 $\delta^{+}_{sub}=5.0$ in this case. Then, the constant $C^{+}_{I}(D^{+}\rightarrow 0)$ shown in equation (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 and Moffat, (1986), in the transitionally-rough regime, the roughness function, $\Delta\left\langle\overline{u}\right\rangle^{+}$, can be estimated as follows:

[TABLE]

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 $\kappa$ is that associated with the smooth-wall simulation ($\kappa\approx 0.41$). The value of $C^{+}_{II}$ can be also estimated on the basis of the expression (2) evaluated at $y=D+y_{0}$:

[TABLE]

If the dependence on the relative submergence is neglected for a moment, then $C^{+}_{II}$ can be interpreted as equal to $B^{+}-\frac{1}{\kappa}\ln\frac{k_{s}}{D}$. The value of $B^{+}$ can be obtained for instance from the diagram of figure 1 of Ligrani and Moffat, (1986), 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., (2011) for the case D50,

the actual value of $k^{+}_{s}$ could be presumably larger than the value presently estimated. An expression similar to (9) can be written also for the fully-rough regime:

[TABLE]

where the values of $\kappa$ and $y_{0}$ are those attained over a smooth-wall, i.e. $0.41$ and [math], respectively. Although the values of $\Delta\left\langle\overline{u}\right\rangle^{+}$ in expressions (9) and (11) are evaluated at different distances from $y_{0}$ (i.e. $\delta_{sub}$ and $D$, respectively), their dependence on the distance in the range $|D-\delta_{sub}|$ was small (i.e. $\sim 0.8$ for $D^{+}\sim 100$). Note that this dependence of $\Delta\left\langle\overline{u}\right\rangle^{+}$ on the distance from the wall is due to deviations of $\kappa$ 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, $\Delta\left\langle\overline{u}\right\rangle^{+}$ monotonically increases with $D^{+}$ for a given configuration of the roughness (e.g. see Jiménez,, 2004), combining equations (9) and (11), the following inequality can be found which poses an upper limit to the value of $C_{II}$:

[TABLE]

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

[TABLE]

where $\kappa=0.41$. The inequality (13), 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 $\frac{1}{\kappa}\ln(k^{+})$.

For the simulation D120, $\Delta\left\langle\overline{u}\right\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 (13) 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 $\varrho U_{bh}^{2}/\tau_{w}$ is relatively large (see table 1).

Table 2 shows the values of $C_{II}^{+}$ and $\frac{1}{\kappa}\ln(D^{+})$ which were also computed for five experiments of Amir et al., (2014). 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., (2014) 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., (2014) are also shown in figure 4a as a function of $y^{+}-y^{+}_{0}$, where $y^{+}_{0}$ is assumed equal to $0.85\leavevmode\nobreak\ 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 (13) 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 (12) is approximately satisfied if also the dependency of $\Delta\left\langle\overline{u}\right\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.

Ignoring for a moment the dependency on $H/D$, as already proposed for the transitionally-rough regime above, it would be found that $\ln\frac{k_{s}}{D}=(8.5-C_{II}^{+})\kappa$ 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, the experiment number 1 of Amir et al., (2014) was characterized by a value of the roughness function relatively small compared with that of Nikuradse, (1933), 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, (1966) and Reynolds, (1974), the flow of simulation D120 falls well into the fully-rough range (see the diagram in figure 3.3a of Pimenta et al.,, 1975, 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.

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 characterize the “three-dimensionality” of the turbulence structure was also adopted by Chan-Braun et al., (2011). Figure 5a, shows the limit of the roughness sublayer compared with the wall-normal profile of stress defined by equation (6) for the streamwise velocity component and for pressure (i.e. the variance of $p^{\prime\prime}-p^{\prime\prime\prime}$ normalized by $(\varrho u_{\tau}^{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 5a. Note that this exponential decay of the three-dimensionality of pressure was also observed by Chan-Braun et al., (2011) 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., (2013) who ascribed it to time convergence error. However, this small residual could be associated with the roughness footprint that Hong et al., (2011) 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 and Castro, (2002) and of Hong et al., (2011) ($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., (2013) 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 center of the spheres, as shown in figure 5b. Although the streamwise velocity $\left\langle\overline{u}\right\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.

The profiles of the mean square of streamwise velocity fluctuations were computed (figures 6a and 6b) by considering the definitions of fluctuations given by (4) and (5), 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 $\left\langle\overline{u^{\prime\prime\prime}u^{\prime\prime\prime}}\right\rangle^{+}$ is located over the crests of roughness elements for the run D50 while it falls below the crests for the run D120 (figure 6a) and, in the latter case, the maximum disappears if the fluctuations $u^{\prime\prime}$ about the average flow field are considered (figure 6b).

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 $\left\langle\overline{u^{\prime\prime\prime}u^{\prime\prime\prime}}\right\rangle^{+}$ from a smooth wall.

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, $\left\langle\overline{u}\right\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, $\left\langle\overline{u^{\prime\prime}u^{\prime\prime}}\right\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 figures 8a,b), although significant turbulent fluctuations are present

in the roughness canopy

in the latter one. It also appears that the maximum of $\left\langle\overline{u^{\prime\prime}u^{\prime\prime}}\right\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 8c which shows the distribution of $\left\langle\overline{u^{\prime\prime}u^{\prime\prime}}\right\rangle^{+}_{B,z}$, where the peaks of figure 8a are filtered out through the spanwise-average.

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, $\left\langle\overline{u}\right\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 center 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 10b 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 vinicity 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 $\left\langle\overline{\boldsymbol{\omega}}\right\rangle^{+}_{B}$.

The vorticity related to the average flow field, $\left\langle\overline{\boldsymbol{\omega}}\right\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, $\left\langle\overline{\omega_{x}}\right\rangle^{+}_{B}$, in figure 12 shows the presence of four streamwise-oriented vortical structures per streamwise-orientated inter-particle gap. In figure 12a,b, which corresponds to a plane orthogonal to the $x$-axis and crossing the center 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 12c,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 12e,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 12f). 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, $\left\langle\overline{\omega_{z}}\right\rangle^{+}_{B}$, is the most intense, since it is associated with the primary mean flow. The most intense regions of $\left\langle\overline{\omega_{z}}\right\rangle^{+}_{B}$ are located at and around the top of the spheres (figure 13a,b

corresponding to the plane $z=z_{A}$ of figure 11a);

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 $\left\langle\overline{\omega_{z}}\right\rangle^{+}_{B}$ are attained in the center plane between the spheres (figure 13c,d,

plane $z=z_{B}$ of figure 11a).

Finally, let us turn to the wall-normal component $\left\langle\overline{\omega_{y}}\right\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 14d), 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 section 3.2.

With the aim of determining the influence of the roughness on the turbulence structure, the distribution of the temporal fluctuations, ($\mathbf{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 $\lambda_{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.,, 1987). The value of $\lambda_{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., (1987) extended the idea of a spatial scale $\lambda^{+}_{uz}$ less intuitively to the other components of the velocity and found that $\lambda^{+}_{wz}\simeq\lambda^{+}_{uz}$ and $\lambda^{+}_{vz}\simeq 2\lambda^{+}_{uz}$. Presently, a similar approach was adopted for both the smooth- and rough-wall simulations (see figure 16). As shown in figure 15, figure 16a and figure 16b, respectively, the values of $\lambda^{+}_{uz}$, $\lambda^{+}_{vz}$ and $\lambda^{+}_{wz}$, at $y^{+}=D^{+}$ for the simulations D50 and D120 were roughly 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 is remarkably similar, as is their intensity.

The value of $\lambda^{+}_{wz}$, extrapolated at $y^{+}_{smooth}=(D^{+}-y^{+}_{0})_{D120}$ for the smooth-wall simulation at $Re_{bH}=6900$, is almost equal to that of $\lambda^{+}_{wz}$ obtained for the fully-rough simulation D120 (see figure 15c). Hence, the relationships observed by Kim et al., (1987) between $\lambda^{+}_{uz}$ and $\lambda^{+}_{vz}$ and between $\lambda^{+}_{uz}$ and $\lambda^{+}_{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 $\lambda^{+}_{uz}$ as a function of $y^{+}$. However, the methodology described above does not allow us to estimate the value of $\lambda^{+}_{z}$ below the crest of roughness elements where velocity fluctuations are not defined continuously on wall-parallel planes.

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, $\lambda^{+}_{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.,, 2011). 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 section 3.1 to the total stress tensor $\boldsymbol{\tau}$ (containing both the viscous and pressure contributions, as defined in equation 7). 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.

[TABLE]

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

[TABLE]

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

[TABLE]

where $\widetilde{\mathcal{S}}$ denotes the surface of the sphere contained in the sphere-box $\mathcal{B}$ (see figure 2 and the definition (28) in the Appendix), $\mathcal{S}^{(i,j)}$ indicates the surface of the $(i,j)$-th sphere centered in $\boldsymbol{x}^{(i,j)}_{c}$ 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 (4). 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}=\varrho u_{\tau}^{2}L^{2}_{B}$.

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 forth moment statistics,

the angle $\alpha_{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}^{(\nu)}$ and $F_{y}^{(\nu)}$, 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_{\tau}$ (see table 1). This can be explained in terms of Darcy-Weisbach’s 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 $\alpha_{F}$. Chan-Braun et al., (2011) 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}$.

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

[TABLE]

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, $(\sigma_{T_{x}},\sigma_{T_{y}},\sigma_{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).

The distribution of the average surface pressure $\langle\overline{p}_{tot}\rangle_{B}$ (i.e. the surface-normal component of the stress vector $\widetilde{\boldsymbol{\tau}}_{n}$ with the sign reversed), normalised by $F_{R}/A_{sph}$, where $A_{sph}=\pi 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 equation 8). While the negative-valued region is slightly shifted downstream of the sphere, the positive region (i.e. with the normal stress directed towards the center 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 section 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 figures 18c,d characterised by peaks of the average pressure. By virtue of the distribution shown in figures 18c,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 center 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 r.m.s. 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 shows the projection of the surface stress-vector $\widetilde{\boldsymbol{\tau}}_{n}$ upon the streamwise direction, i.e. the map of the local surface-stress contribution to drag, $\widetilde{\boldsymbol{\tau}}_{n}\cdot\mathbf{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 section 3.1 about the shift of $\left\langle\overline{u^{\prime\prime\prime}u^{\prime\prime\prime}}\right\rangle^{+}$ beneath the crest of the spheres in the fully-rough regime.

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.

[TABLE]

where $\eta$ denotes the wall-normal coordinate defined in the range $[0,D]$ and $\theta$ the azimuthal angle. Obviously $\boldsymbol{f}(D)=\boldsymbol{F}$, as defined in (15). 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}^{(\nu)}(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.

Chan-Braun et al., (2011) 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 $\Gamma_{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 $\sigma_{T_{x}}/T_{R}$, $\sigma_{T_{y}}/T_{R}$ and $\sigma_{T_{z}}/T_{R}$ for the simulations D50 and D120 are compared with those of the torque:

[TABLE]

and

[TABLE]

obtained for the smooth-wall square tile $\widetilde{\Gamma}_{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{\Gamma}_{L_{B}}$. The symbol $\boldsymbol{e}_{y}$ in the equations (19-20) indicates the wall-normal unit vector of the canonical base while $\boldsymbol{\tau}$ and $\boldsymbol{\tau}_{\nu}$ denote the total stress tensor and the viscous stress tensor (i.e. without the pressure contribution), respectively, cf. the definition in (7). 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 $\sigma_{\mathcal{T}^{(p)}_{x}}(L^{+}_{B})$ and $\sigma_{\mathcal{T}^{(p)}_{z}}(L^{+}_{B})$, obtained considering also the contribution of pressure, normalised by $\frac{1}{2}\varrho u_{\tau}s^{3}$, 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 22b shows graphically that the standard deviations of the three components of torque, $\boldsymbol{\sigma}_{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 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$-compontent of $\widetilde{\boldsymbol{\tau}}_{n}$ (this is the global wall-normal direction), which is visualised in figure 23, exhibits essentially the same distribution as the pressure.

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 section 3.1 and shown in figure 14. Here we relate these features to the distribution of the spanwise component of $\widetilde{\boldsymbol{\tau}}_{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 $\sigma_{F_{z}}/F_{R}$ and of $\sigma_{T_{y}}/T_{R}$ (cf. tables 3-4). In particular, the increase of the value of $\sigma_{T_{y}}/T_{R}$ in the fully-rough regime along with the decrease of $\sigma_{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., (2014), 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_{\tau}^{2}$. By normalizing the root mean square value of pressure fluctuations with $\varrho u_{\tau}^{2}$ they defined the coefficient $k_{L}$, which was independent of $U_{bh}/u_{\tau}$ 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., (2014), we have computed the value of the coefficient $k_{L}$ in an analogous manner (i.e. as the normalized r.m.s. 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., (2014), 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., (2014) 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., (2014) (which is different from that used by Chan-Braun et al.,, 2013), 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.72\,U_{bH}$ and $0.64\,U_{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., (2014) for their experiments characterised by the relative submergence $H/D=5$. Since Amir et al., (2014) showed that the speed $U_{c}$ was essentially independent of $U_{bh}/u_{\tau}$ (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^{\prime}_{x}F^{\prime}_{y}}(\Delta r_{x}=i\,L_{B},\Delta r_{z}=j\,L_{B},\Delta_{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., (2014) could equip each sphere with a one-direction pressure probe only, they could estimate $R_{F^{\prime}_{x}F^{\prime}_{y}}$ 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^{\prime}_{x}F^{\prime}_{y}}(0,L_{B},\Delta_{T})$ was computed for the simulations D50 and D120 (see broken lines in figure 25). Although the present value of $R_{F^{\prime}_{x}F^{\prime}_{y}}$ is approximately two times larger than that observed by Amir et al., (2014) (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.14\,U_{bH}/H$ and to $0.16\,U_{bH}/H$, respectively. Moreover, the single-sphere cross-correlation function $R_{F^{\prime}_{x}F^{\prime}_{y}}(0,0,\Delta_{T})$ was computed (see solid thick lines in figure 25) which revealed the same time lag of the positive (negative) maximum ($\Delta_{T}=0.16\,H/U_{bH}$ and $\Delta_{T}=-0.14\,H/U_{bH}$ for cases D50 and D120), and the presence of a secondary negative (positive) peak at $\Delta_{T}=0.5\,H/U_{bH}$ ($\Delta_{T}=-0.5\,H/U_{bH}$). The trend of $R_{F^{\prime}_{x}F^{\prime}_{y}}$ 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, (2005), Dwivedi, (2010) as well as by Amir et al., (2014), while the cross-correlation function for case D50 was also shown in figure 8 of Chan-Braun et al., (2013).

In case D120 the average time interval between positive and negative fluctuations of lift, equal to $\Delta_{T}=0.33\,H/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,\Delta_{T})$ (not shown here). In Chan-Braun et al., (2013) the value of $\Delta_{T}=0.42H/U_{bH}$ is measured for their case D50. Consistently, nearly after the same time interval ($\Delta_{T}=0.37\,H/U_{bH}$), the cross-correlation $R_{F^{\prime}_{x}F^{\prime}_{y}}$ is found to vanish in both cases D50 and D120. Furthermore, the symmetry of $R_{F^{\prime}_{x}F^{\prime}_{y}}$ with respect to the point $\Delta_{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.,, 2013). Thereby, it is possible to estimate the length scale, (henceforth denoted as $\ell$) as the product of the time interval $\Delta_{T}=0.3\,H/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.22\,H$ ($1.24\,D$) and $0.19\,H$ ($1.03\,D$), respectively. Furthermore, in the spanwise direction the two-point correlation of lift and drag was found negligible at the separation distance $2\,D$ 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.

With the purpose to help 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 center. In this model, turbulent coherent structures induce pulsations of the intensity of the vortices characterizing 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., (2013) 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 $\Delta_{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^{\prime}_{x}F^{\prime}_{y}}$, 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 $\Delta_{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 $\Delta^{(F)}_{T}-\Delta^{(E)}_{T}$ equal to $0.30$ and $0.36$ for the simulations at $Re_{bH}=2900$ and $6900$, respectively, where $\Delta^{(E)}_{T}$ and $\Delta^{(F)}_{T}$ 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.71\,U_{bH}$ and $0.69\,U_{bH}$, respectively. Hence, the length scale of the “smaller” turbulent structures was found equal to $0.20\,H$ and to $0.25\,H$ for the simulations at $Re_{bH}=2900$ and $6900$, respectively, which are nearly equal to $L_{B}$. Indeed, Chan-Braun et al., (2011) 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.

4   Conclusions

Direct numerical simulation of open-channel flow has been performed over a bed of spheres in square arrangement for values of the bulk and roughness Reynolds numbers of approximately $6900$ and $120$, respectively, for which the fully-rough flow regime was attained. The objective of the present work is to quantify the differences in the flow structure which arise when passing from the transitionally-rough flow regime to the fully-rough one, providing us with an enhanced picture of the flow-roughness interaction. Therefore, one of the DNSs performed by Chan-Braun et al., (2011) in the transitionally-rough regime ($Re_{bH}=2900$ and $D^{+}=50$), characterised by the same shallowness and arrangement of the spheres as the present DNS, is extensively referred to. Furthermore, smooth-wall data from simulations at the same bulk Reynolds numbers and computational domain size are used for the purpose of comparison. The two DNSs, in the transitionally- and fully-rough regimes, were compared by using different statistical tools. In particular, three average operators were employed, namely the simple time-average, the periodic-sphere-box-average and the wall-normal plane-average (the latter two combined with the time-average), each extracting different statistical features of the flow.

When increasing the Reynolds number it is observed that the flow penetrates deeper into the crevices between the roughness elements, causing a secondary flow in the vicinity of the spheres which consists of two pairs of recirculation cells. A downward shift of the mean flow profile is observed, which results in the disappearance of the viscous sublayer and of the buffer layer, while viscous effects were confined to the region below the crests of the spheres where velocity fluctuations were found to be strongly correlated. In fact, velocity fluctuations around the periodic-sphere-box-averaged flow field (i.e. those which exclude the average flow field at the roughness element scale) were found to be of nearly negligible intensity below the crests, whereas those around the plane-average were large. This fact indicates that vortex structures in the interstitial flow are weakly turbulent, but mostly associated with the average shear-layer originating from the sphere surface. Although, the damping of turbulent fluctuations in the interstitial fluid region does not significantly affect the flow structure in the transitionally-rough regime, because normal stresses attain the maximum intensity over the crest of the roughness elements, it definitely does in the fully-rough regime, since the peak of normal stresses would have been placed below the crest of the spheres. This is a clear indicator of the destruction of the buffer layer and, consequently, that the flow regime is fully-rough (only a small relative maximum of the root mean square of the streamwise velocity is present above the crest of the spheres).

Moreover, the flow in the log-law region was investigated. It was found that, for the present configuration, the von Kármán constant deviates somewhat from the value that is attained in absence of the roughness, attaining the value $0.381$ in the fully-rough regime.

Note that this value is still within the bounds ($\kappa=0.39\pm 0.02$) determined by Marusic et al., (2013) from a data-set including various boundary layer and pipe flow experiments as well as measurements in the atmospheric surface layer, covering a considerable range of Reynolds numbers under nominally smooth and, in the latter case, transitionally-rough conditions ($k_{s}^{+}\approx 21$).

By comparing the velocity profile obtained from simulations over a smooth and a fully-rough wall it was possible to quantify the conductivity of the present configuration with respect to the ideal case of maximum conductivity. The integration constant $C_{I}^{+}$ appearing in the law of the wall (1) tends to zero with increasing particle Reynolds number while the conductivity $C_{II}^{+}$ (which includes the effects associated with the geometrical configuration of the domain) approaches $11.6$. Presently, the value of $C_{II}^{+}$ was found to be equal to $10.7$, indicating that, for the present values of $H/D$ and for a well-packed square arrangement of spheres, the roughness geometry is highly conductive.

The vorticity field, averaged in time and over periodic boxes around the spheres, has been discussed in detail. It was found that on average the vortical structures in the spheres’ wakes reach significantly further downstream at the larger particle Reynolds number, clearly reconnecting with the surface of the downstream neighbor sphere. Furthermore, in case D120 a second set of four counter-rotating vortices appears in the inter-particle grooves due to the fact that the primary mean flow velocities are enhanced in that region.

Low- and high-speed streaks were observed over the crest of the spheres similarly to those forming over a smooth wall. These structures are characterised by a spanwise spacing increasing linearly with the distance from the spheres. In particular, the spacing between the streaks at the crest level ($y^{+}=D^{+}$) is almost the same as that observed over a smooth wall at the same Reynolds number when the same distance from the virtual wall is chosen (i.e. $y^{+}_{smooth}=(D^{+}-y_{0}^{+})_{rough}$). In general, the relationships between the minima locations of the spanwise two-point correlations of the different velocity components $\lambda_{ux}^{+}$, $\lambda_{uy}^{+}$ and $\lambda_{uz}^{+}$ shown by Kim et al., (1987) for channel-flow over a smooth wall were found to hold also over a fully-rough wall, but their individual magnitudes are modulated by the presence of the roughness.

The action of the flow on individual roughness elements has been investigated in detail, including an investigation of the spatial distribution of stress statistics on the surface of the spheres. The coefficient of the average drag force acting on the spheres remains almost constant when the bulk Reynolds number is increased from $2900$ to $6900$ (for a fixed value of the relative submergence), which is in line with the Darcy-Weisbach equation. At the same time the average lift coefficient increases by approximately $30\%$ in this interval, mostly due to the contribution of pressure acting on the upper upstream part of the sphere surface. Concerning the distribution of the streamwise component of the stress vector at the sphere surface, we have observed that in the fully-rough regime the most intense contribution shifts towards smaller wall-distances (downwards) and towards the upstream-facing part of the upper side of the spheres, as compared to the transitionally-rough case. This same region approximately coincides with the region where the average pressure on the sphere’s surface attains relatively large positive values (note the definition of pressure in 8). Compared to the transitionally-rough reference case, the contribution of the skin-friction to the net drag force acting on the roughness elements is reduced by nearly $15\%$ in the fully-rough regime, while the pressure contribution becomes dominant.

An analogy with the force/torque acting upon an element of a smooth-wall channel flow, which was successfully used by Chan-Braun et al., (2011) to explain e.g. the statistical moments of the torque acting on the roughness elements, is found to break down in the fully-rough regime. The breakdown of the analogy with a smooth wall, which is not unexpected, can be attributed to the roughness-induced destruction of the buffer layer and to the reduction of the importance of viscous effects in the vicinity of the sphere crests.

The results obtained through DNS in the fully-rough regime were compared with those of experiments by Amir et al., (2014) in an open-channel flow over a bed of immobile spheres in a hexagonal arrangement at larger particle Reynolds numbers. The convective velocity of turbulent structures, $U_{c}$, was estimated from the time correlation of the hydrodynamic forces acting on neighbouring spheres and showed values in fair agreement with those obtained by Amir et al., (2014). This supports the observation that the convective velocity depends only on the relative submergence while it is independent of the Reynolds number and, on the basis of the present results, it also appears to be largely independent of the precise arrangement of the roughness elements. Indeed, this was, to the knowledge of the authors, the first time that results obtained by DNS were used to supplement experimental observations in a fully-rough open-channel flow.

Finally, we have investigated the cross-correlation between fluctuations of drag and wall-normal lift forces, which are known to be highly correlated for non-zero separation times. Our data reveals a periodic sequence of events with an average period of approximately $0.64H/U_{bH}$. Such a cyclic process was explained as a manifestation of turbulence below the crest of the spheres, consisting in pulsations of the intensity of the spanwise-oriented vortices adhering to the sphere surface which were identified in the sphere-box/time-averaged flow field.

Although, DNS in the fully-rough regime at even larger Reynolds numbers than presently considered are already possible, we belive that it would be more relevant for the purpose of studying the flow resistance in natural rivers to try to remove the condition of regularity of the arrangement of the roughness elements. Roughness elements of similar size and shape in a “random” arrangement will presumably have a different interaction with the flow with respect to the present cases. Moreover, the shelter effect of neighbouring roughness elements on each other could be investigated and quantified when considering a random bed.

This study has been funded by the Deutsche Forschungsgemeinschaft (project UH 242/4-2). We acknowledge the generous support from the Leibniz-Rechenzentrum (LRZ, Munich) for the extensive computational resources on SuperMUC (project pr87yo) and from the Steinbuch Centre for Computing (SCC, Karlsruhe), in particular for generous allotment of storage space. The authors wish to thank Clemens Chan-Braun, Manuel Gracía-Villalba and Aman Kidanemariam for the helpful discussions throughout the elaboration of the present work. Thanks is also due to Clemens Chan-Braun for initially setting up the simulation D120.

Appendix A Definition of average operators

Let $q(\boldsymbol{X},T)$ be a quantity defined in the present computational domain as a function of the random variables $\boldsymbol{X}$ and $T$, and let us consider the following events

[TABLE]

where $\boldsymbol{x}^{(i)}_{c}$ are the coordinates of the center of the $i-$th sphere, $\mathcal{B}$ is the rectangular domain defined in section 2, and $\phi$ denotes an indicator function which measures unity at points occupied by fluid and vanishes otherwise. Note that in the statement (21) the position of the (fixed) sphere $\boldsymbol{x}^{(i)}_{c}$ is not a random variable, and, due to the perfectly square arrangement, turbulence statistics are independent of the position of the $i-$th sphere, and the probabilty $\mathbb{P}(E_{1})$ is constant and equal to $1/N_{s}$. Moreover, $q$ is assumed to be an ergodic process such that its statistics converge for a “sufficiently long” sampling time and the dependence on time can be removed by applying the temporal average operator

[TABLE]

where $\mathcal{T}$ denotes the time interval during which $q$ is sampled. For the sake of clarity, hereafter the dependence on time is implied if an overline is not present.

Then, two space average operators are defined on the basis of the events (21-23):

• •

The sphere-box-average $\widetilde{\left\langle q\right\rangle}_{B}(\widetilde{\boldsymbol{x}})$ is the expected value of $\phi q$ when the conditioned event $E_{1}|E_{2}$ occurs

[TABLE]

where $\widetilde{\boldsymbol{x}}\in\mathcal{B}$ and $\phi$ is

[TABLE]

Because of the streamwise and spanwise periodicity of the sphere arrangement and of the boundary conditions, we note that the function $\phi\left(\boldsymbol{X}\right)$ is unaffected by the translation of coordinates $\phi(\boldsymbol{X}-\boldsymbol{x}^{(i)}_{c})$ because $\phi$ is the periodic extension of the function $\widetilde{\phi}$ defined in the random variable $\widetilde{\boldsymbol{X}}\in\mathcal{B}$. Similarly, $\widetilde{\left\langle q\right\rangle}_{B}$ can be extended periodically in the streamwise and spanwise directions such that the resulting sphere-box-averaged flow field $\left\langle q\right\rangle_{B}(\boldsymbol{x})$ equals $\widetilde{\left\langle q\right\rangle}_{B}(\boldsymbol{x}-\boldsymbol{x}^{(i)}_{c})$ for every $i=1,\ldots,N_{s}$. Thus, since $\mathbb{P}(\phi=1)=\mathbb{P}(\widetilde{\phi}=1)$ and $\mathbb{P}(E_{2}|E_{1})=\mathbb{P}(E_{2})$, the sphere-box-average of $q$ can be equally expressed in the two following forms

[TABLE]

• •

The plane-average $\left\langle q\right\rangle(y)$ is defined as the expected value of $\phi q$ when the conditioned event $E_{3}|E_{1}$ occurs, i.e. the mean value of $\phi q$ over wall-parallel planes excluding the points where $\phi$ vanishes. Similarly the streamwise- and spanwise-averages can be defined under the further conditions $Z=z$ and $X=x$ and indicated by $\left\langle q\right\rangle_{x}(y,z)$ and $\left\langle q\right\rangle_{z}(x,y)$, respectively.