3.1 Flow structures

In oscillatory flows over rough surfaces, turbulence is generated by the effects of size of the roughness and also by high-amplitude oscillations. As per Jensen et al. (Reference Jensen, Sumer and Fredsøe1989), in fully developed oscillatory turbulent flows over a rough boundary, most of the flow quantities are therefore dependent upon one single parameter, namely $a/k_{s}$ . As noted in table 1, $a/k_{s}$ in the SP case is around 3.5 times larger than that in the LP case. As a result, there exist some discernible differences in the elementary processes that govern the near-bed turbulence in the large- and the small-particle cases.

For the coarse-particle LP case, $a/k_{s}$ is smaller (large roughness size) compared to that for the fine-particle SP case. The size of the individual coarse particles is large enough to cause vortex shedding that is in a way also responsible for triggering the turbulence much earlier in the flow cycle (Ghodke & Apte Reference Ghodke and Apte2016b ). Contrary to this, as a result of large $a/k_{s}$ value (smaller roughness size) in the SP case, the near-bed flow exhibits an unsteady inertial behaviour in the early acceleration cycle, until the local friction velocity is large enough to trigger the turbulence. This behaviour is evident from figure 3, where the product of instantaneous streamwise ( $u$ ) and wall-normal ( $w$ ) velocity is plotted for (a) LP and (b) SP cases. The significant vortex shedding along with sweep–burst motions in the LP case result in a strong cross-correlation signal $uw$ from the beginning of the flow cycle. This indicates that the turbulence is present during most of the oscillation cycle for the LP case. However, as shown in figure 3(b), such cross-correlation is significant for the shear-dominated SP case only from peak phase to late deceleration phases where inertial effects play an important role in triggering the turbulence. This shows that the turbulence is present only in part of the oscillation cycle for the SP case.

Figure 3.

Variation over a half cycle of the normalized product of instantaneous streamwise ( $u$ ) and wall-normal ( $w$ ) velocity plotted for (a) LP and (b) SP cases.

Figure 4.

Close-up view of the instantaneous isosurfaces of the $\unicode[STIX]{x1D706}_{2}$ parameter coloured by $u^{\prime }/U_{\infty }$ for (a) LP and (b) SP cases plotted at a peak phase, $\unicode[STIX]{x1D714}t=5\unicode[STIX]{x03C0}/10$ . The flow direction is the positive $x$ -direction.

To further elucidate the near-bed flow behaviour, instantaneous flow structures in the form of $\unicode[STIX]{x1D706}_{2}$ isosurfaces (Jeong & Hussain Reference Jeong and Hussain1995) are plotted at the peak phase in figure 4 for (a) LP and (b) SP cases. For the LP case, there exists a complex forest of highly dense near-bed structures that appear to be broken and are slightly inclined with respect to the horizontal direction. The disruption of these flow structures is a consequence of energized wall-normal velocity fluctuations in the LP case that in turn distort the streamwise structures, resulting in flow isotropization, as already discussed in Ghodke & Apte (Reference Ghodke and Apte2016b ). On the other hand, the presence of larger elongated horseshoe structures that span over more than a few particle diameters in the flow direction along with the small-scale broken structures is clearly evident in the SP case shown in figure 4(b). Contrary to the LP case, these structures appear to be less broken, suggesting a lower degree of near-bed flow isotropization in the SP case. Such near-bed structures could play a significant role in influencing spatio-temporal characteristics of hydrodynamic forces on the particle bed. Detailed quantification of the near-bed turbulence structure in the LP and SP cases is presented in Ghodke & Apte (Reference Ghodke and Apte2018).

3.2 Particle-bed force history

In the following, the nature of unsteady forces on the particle bed is analysed. The hydrodynamic force, $\boldsymbol{F}$ , acting on a particle is defined as

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}=\int _{\unicode[STIX]{x1D6E4}}\unicode[STIX]{x1D749}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}\unicode[STIX]{x1D6E4}, & \displaystyle\end{eqnarray}$$

where $\text{d}\unicode[STIX]{x1D6E4}$ is the unit surface area of the particle, $\boldsymbol{n}$ is the surface normal vector and $\unicode[STIX]{x1D749}$ is the stress, containing both viscous and pressure contributions, on the particle surface. Figure 5 shows the phase variation of instantaneous drag ( $C_{D}$ ) and lift ( $C_{L}$ ) force coefficients on each particle (thin black lines) plotted along with the respective average force coefficients on the particle bed (thick red line). The force coefficient is defined as $C_{i}=F_{i}/(0.5\unicode[STIX]{x1D70C}A_{p}U_{\infty }^{2})$ , where $i$ stands for drag or lift with $F$ as the respective force and $A_{p}$ is the projected area of the particle. For the purpose of illustration, data for four cycles are plotted. The phase variation of the streamwise velocity far from the bed is also shown for reference. Figures 5(a) and 5(b), respectively, show the drag and lift coefficients in the LP case and figures 5(c) and 5(d) show those for the SP case. Contrary to unidirectional flows, hydrodynamic forces on the particle bed in oscillatory flows vary over a cycle. Unlike drag, the sign of the lift force is the same whether the flow is in the positive or the negative $x$ -direction and has a period of half the forcing function. As seen in figure 5(a,b), the variability of the instantaneous forces on each sphere from the average value is consistent throughout the entire cycle for the larger-diameter particle case. On the other hand, for the smaller-diameter particle as shown in figure 5(c,d), such variability of the instantaneous forces is significant only from late acceleration up to early deceleration phases; notably, during the same period of the oscillation cycle where the cross-correlation signal $uw$ is significant. Such differences in the force history for the LP and SP cases are therefore mainly a result of the different nature of the near-bed turbulence in these flow cases, as discussed earlier.

Figure 5.

Variation over several cycles of the normalized instantaneous drag $C_{D}$ (a,c) and lift $C_{L}$ (b,d) for each particle in the flow set-up. Panels (a,b) are for the LP case and panels (c,d) are for the SP case. The thick red line indicates the average value of the force. The phase variation of the streamwise velocity far from the bed is shown (— ⋅ —) in the lowermost panels.

To further elucidate the differences in the nature of the forces for large- and small-particle cases, the phase variation of the average force on the particle bed in these flow cases is discussed in detail. Figure 6 shows the variation over a cycle of the phase-averaged drag ( $\overline{C}_{D}$ ) and lift ( $\overline{C}_{L}$ ) coefficients for the coarser (LP) and the finer (SP) roughness elements. The phase variation of the streamwise velocity far from the bed is also shown for reference. In the LP case shown in figure 6(a), both drag and lift coefficients peak earlier in the flow cycle, respectively exhibiting a phase lag of $\unicode[STIX]{x1D714}t=0.87$ and $0.34$ with the free-stream velocity. In comparison, as a result of increased Reynolds number $Re_{a}$ in the shear-dominated SP case, a reduced phase lag of, respectively, $\unicode[STIX]{x1D714}t=0.08$ and $0.086$ for drag and lift with the free-stream velocity is observed in figure 6(b). Note that the phase lag between the drag and lift forces also decreases as the Reynolds number is increased, as in the SP case.

Figure 6.

Variation over a cycle of the normalized phase- and spatially averaged drag $\overline{C}_{D}$ (

) and lift $\overline{C}_{L}$ (

) for (a) LP and (b) SP cases. Here, $\overline{C_{i}}=\langle \overline{F_{i}}\rangle /(0.5\unicode[STIX]{x1D70C}A_{b}U_{\infty }^{2})$ is the normalized phase- and spatially averaged force on the particle bed, $i$ stands for drag or lift, and $A_{b}$ is the total projected area of the particle bed. The phase variation of the streamwise velocity (not to scale) far from the bed is shown (— ⋅ —). Arrows in panel (b) indicate spikes in the drag and lift record for the SP case.

As seen in figure 6(b), there exist sharp spikes in the drag and lift records for the SP case (shown by arrows); however, such distinct spikes are absent in the LP case. As a result of significant vortex shedding along with sweep–burst motions in the LP case, there is no clear demarcation of onset of turbulence that may result in distinct spikes in the drag and lift records. Contrary to this, the occurrence of sweep–burst motions resulting in sudden transition to turbulence close to the peak velocity phase results in distinct spikes in the drag and lift records in the SP case. Note that the phase instant of these spikes in figure 6(b) clearly matches well with the occurrence of spikes in the $uw$ signal for the SP case shown in figure 3(b).

Table 2.

Force statistics for the LP and SP cases. Here, $\overline{C}_{D,max}$  ( $\overline{C}_{L,max}$ ) and $\overline{C}_{D,avg}$  ( $\overline{C}_{L,avg}$ ), respectively, show maximum and period-averaged values of the drag (lift) force coefficients over half a flow cycle; and $S_{D^{\prime }}$  ( $S_{L^{\prime }}$ ) and $K_{D^{\prime }}$  ( $K_{L^{\prime }}$ ), respectively, denote skewness and kurtosis of the drag (lift) fluctuations over half a flow cycle.

Table 2 summarizes the particle-bed force statistics for the LP and SP cases. In both flow cases, the net lift and drag coefficients remain positive and their magnitudes are decreased as the Reynolds number is increased. Higher-order moments such as skewness and kurtosis of drag/lift fluctuations over several half flow cycles are also reported in table 2. In both flow cases, the skewness values of the drag and lift forces are positive, indicating that the positive events of drag and lift are more likely to occur than the negative ones. Moreover, the high values of kurtosis of drag and lift fluctuations indicate large intermittency in the force field that results in longer tails of rare but strong events. With increase in Reynolds number (SP case) such extreme events become stronger, resulting in higher values of kurtosis. Note that kurtosis values in the SP case are almost twice as large as in the LP case. Such extreme force events might play an important role in destabilizing the particle bed. Similar observations were reported by Chan-Braun et al. (Reference Chan-Braun, Garcia-Villalba and Uhlmann2011) in their study of unidirectional flow over a rough bed.

In addition, the present results of the LP case are compared against data from Mazzuoli & Vittori (Reference Mazzuoli and Vittori2016). Figure 7 shows a comparison of normalized phase- and spatially averaged drag force for $Re_{\unicode[STIX]{x1D6FF}}=400$ (LP) case with the data from Mazzuoli & Vittori (Reference Mazzuoli and Vittori2016). As shown in figure 7, the present results compare reasonably well with Mazzuoli & Vittori (Reference Mazzuoli and Vittori2016), albeit some discrepancies are seen in the magnitude of the force in the decelerating part of the cycle. This could be due to fewer flow cycles (only three) being considered by Mazzuoli & Vittori (Reference Mazzuoli and Vittori2016) as against 15 flow cycles in the present investigation. Nevertheless, the overall phase variation of the average force, including the phase of the peak force, in the present case matches fairly well with Mazzuoli & Vittori (Reference Mazzuoli and Vittori2016).

Figure 7.

Variation over a half cycle of the normalized phase- and spatially averaged drag force for $Re_{\unicode[STIX]{x1D6FF}}=400$ (LP) case. Data from Mazzuoli & Vittori (Reference Mazzuoli and Vittori2016) are shown by symbols. Here, force is normalized by $\unicode[STIX]{x1D70C}U_{\infty }^{2}\unicode[STIX]{x1D6FF}^{2}$ .

3.3 Spatial distribution of the forces

After studying the phase variation of the drag/lift forces, their distribution on the particle surface needs to be investigated. Figure 8 shows the spatial distribution at a peak phase of drag and lift components, $\bar{\unicode[STIX]{x1D70F}}_{D}$ and $\bar{\unicode[STIX]{x1D70F}}_{L}$ , respectively, of the phase-averaged stress ( $\bar{\unicode[STIX]{x1D749}}$ ) on the particle surface. Here, $\bar{\unicode[STIX]{x1D70F}}_{D}=\bar{\unicode[STIX]{x1D749}}\boldsymbol{\cdot }\boldsymbol{e}_{1}$ and $\bar{\unicode[STIX]{x1D70F}}_{L}=\bar{\unicode[STIX]{x1D749}}\boldsymbol{\cdot }\boldsymbol{e}_{3}$ , where $\boldsymbol{e}_{1}$ and $\boldsymbol{e}_{3}$ are unit vectors in the $x$ - and $z$ -directions, respectively. Figure 8(a,b) shows isometric views and figure 8(c,d) presents top views of the particles. Figures 8(a,c) and 8(b,d) show the stress distributions on the large and small particles, respectively. Note that $\bar{\unicode[STIX]{x1D70F}}_{D}$ and $\bar{\unicode[STIX]{x1D70F}}_{L}$ plotted in figure 8 indicate phase- and spatially averaged stresses. Therefore, the spheres shown in this figure do not represent a particular particle but rather show an averaged stress field on any particle in the set-up. From the variation of forces plotted in figure 6, it can be seen that the average drag/lift forces on the particle at the peak phase in both flow cases are positive. However, as shown in figure 8, local stress components contributing to drag and lift can be either positive or negative depending on the direction of the local velocity near the particle surface.

Figure 8.

Spatial distribution of the phase- and spatially averaged local stress on the particle surface corresponding to drag, $\bar{\unicode[STIX]{x1D70F}}_{D}$ , and lift, $\bar{\unicode[STIX]{x1D70F}}_{L}$ , at a peak phase normalized by their respective standard deviations for the LP (a,c) and the SP (b,d) cases. Here, $D$ and $L$ , respectively, denote drag- and lift-related stress components. The direction of the flow is indicated by an arrow. The sizes of the particles are not to scale.

For both flow cases, there exists a region of strong positive drag component of the stress, $\bar{\unicode[STIX]{x1D70F}}_{D}$ , near the top of the sphere, a result of maximum strain in the streamwise velocity component at this location. In addition, small patches of weak positive $\bar{\unicode[STIX]{x1D70F}}_{D}$ are also found especially in the front and aft region of the large particle as shown in figure 8(a,c). These local regions of weak drag are a consequence of pore-scale flow features and are dominant mainly for the large particle due to pronounced curvature effects; however, they are clearly absent in the small-particle case shown in figure 8(b,d). Furthermore, variation of $\bar{\unicode[STIX]{x1D70F}}_{D}$ in the wall-normal direction reveals the location of the effective level at which the bulk drag force is exerted, also called the zero-displacement plane. In the present study, this zero-displacement plane is found to be located at $0.7D$  ( $0.85D$ ) for LP (SP) case. This further indicates that the total drag force is distributed on a larger area of the upper hemisphere of the large-diameter particle compared to that in the case of the small-diameter particle.

The spatial distribution of the lift component of the stress, $\bar{\unicode[STIX]{x1D70F}}_{L}$ , also shows some discernible differences between the LP and the SP case. As shown in figure 8(b,d), a region of positive magnitude of $\bar{\unicode[STIX]{x1D70F}}_{L}$ is seen mostly in the top portion of the upper hemisphere of the small-diameter particle and is a result of lower values of pressure due to flow acceleration near the particle crest. On the other hand, the $\bar{\unicode[STIX]{x1D70F}}_{L}$ distribution in the case of the large-diameter particle shows larger positive magnitudes on the side and aft region of the particle as shown in figure 8(a,c). This could again be attributed to the combined effects of enhanced three-dimensionality of the boundary layer, shielding effects in a hexagonal arrangement of the particles along with the complex pore-scale flow resulting from broken near-bed structures in the gaps between the particles. As a result, lower-pressure regions are mainly observed on the side and aft portion of the large-diameter particle that contribute positively to the lift. In comparison to the SP case, these effects are more pronounced for the large-diameter case due to strong curvature effects.

3.4 Temporal correlations of the forces

In the following, temporal characteristics of the force fluctuations in the LP and SP cases are studied. The time scales of drag and lift fluctuations in both flow cases are obtained from their respective temporal autocorrelations. In addition, cross-correlation of drag and lift fluctuations is also reported to reveal their interdependence on each other.

Following Chan-Braun et al. (Reference Chan-Braun, Garcia-Villalba and Uhlmann2013), the three-dimensional correlation function given as

(3.2) $$\begin{eqnarray}\displaystyle \displaystyle & & \displaystyle R_{\unicode[STIX]{x1D719}^{\prime }\unicode[STIX]{x1D713}^{\prime }}(\unicode[STIX]{x0394}x,\unicode[STIX]{x0394}y,z,\unicode[STIX]{x0394}t)\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{1}{N_{c}N_{p}}\mathop{\sum }_{i=1}^{N_{c}}\mathop{\sum }_{j=1}^{N_{p}}\unicode[STIX]{x1D719}^{\prime }(x_{p,j},y_{p,j},z,t_{i})\,\unicode[STIX]{x1D713}^{\prime }[(x_{p,j}+\unicode[STIX]{x0394}x),(y_{p,j}+\unicode[STIX]{x0394}y),z,(t_{i}+\unicode[STIX]{x0394}t)]\end{eqnarray}$$

is defined, where the sub-index $p$ refers to the particle position, and $N_{p}$ and $N_{c}$ are the number of particles and flow cycles, respectively. Using this function, temporal autocorrelation at zero spatial separation can be calculated by setting $\unicode[STIX]{x1D719}^{\prime }=\unicode[STIX]{x1D713}^{\prime }$ and $\unicode[STIX]{x0394}x=\unicode[STIX]{x0394}y=0$ in (3.2). Figure 9 shows the variation of the autocorrelation function of drag and lift fluctuations for (a) LP and (b) SP cases.

Figure 9.

Temporal autocorrelation of drag (——) and lift (- - - -) fluctuations for (a) LP and (b) SP plotted as a function of separation in time. Only a small fraction of the signal in which non-negligible correlation exists is plotted.

In the case of LP, autocorrelation functions for drag and lift exhibit local minima of small negative values (drag, $-0.066$ at $\unicode[STIX]{x0394}tU_{\infty }/d=3.05$ ; and lift, $-0.039$ at $\unicode[STIX]{x0394}tU_{\infty }/d=2.87$ ). For the SP case, the local minimum in the drag takes a positive value ( $0.087$ at $\unicode[STIX]{x0394}tU_{\infty }/d=2.1$ ) and that for the lift takes a negative value ( $-0.023$ at $\unicode[STIX]{x0394}tU_{\infty }/d=3.3$ ). Chan-Braun et al. (Reference Chan-Braun, Garcia-Villalba and Uhlmann2013), in their DNS study of unidirectional flow over a rough bed, also observed such local minima in drag/lift autocorrelations. They argued that the minimum peak in drag/lift autocorrelations might be linked to the effects of pressure; a result analogous to the pronounced minimum of negative value of pressure autocorrelation on a smooth wall in a turbulent channel flow. The present results also show similar trends; however, the values of the local minima are typically 3–4 times smaller than those reported by Chan-Braun et al. (Reference Chan-Braun, Garcia-Villalba and Uhlmann2013). It should be noted that the unidirectional flow study reported by Chan-Braun et al. (Reference Chan-Braun, Garcia-Villalba and Uhlmann2013) falls in the hydrodynamically smooth and transitional regime, whereas the Reynolds numbers in the present cases are dramatically higher and fall in the very rough turbulent regime.

Using the temporal autocorrelation function, two different time scales, integral ( $t_{l}$ ) and Taylor micro-scale ( $t_{\unicode[STIX]{x1D706}}$ ), for drag and lift fluctuations are calculated. The integral time scale, defined as the integral of the autocorrelation, is computed in the integration domain within which the value of the autocorrelation function monotonically decreases and goes to significantly smaller value (here, $\unicode[STIX]{x0394}tU_{\infty }/d\in [0,20]$ ). The short-time behaviour of the autocorrelation function can be quantified using the Taylor micro-scale ( $t_{\unicode[STIX]{x1D706}}$ ), defined by the zero crossing of the parabola osculating the correlation function at zero separation.

Table 3 tabulates integral time scales and temporal Taylor micro-scales of drag/lift fluctuations for the LP and SP cases. The following observations could be noted by comparing these two cases: (i) The integral time scales of force fluctuations for the LP case are typically smaller than those for the SP case. This could be attributed to the different nature of the near-bed turbulence in these flow cases. The broken streamwise structures of the time scale ${\sim}d/U_{\infty }$ in the LP case are responsible for inducing drag/lift fluctuations on the large-diameter particle, resulting in integral time scales almost close to unity. On the contrary, the presence of elongated streamwise structures of the size of several particle diameters results in larger values of integral time scales of force fluctuations in the SP case. (ii) As a result of significant flow isotropization, the Taylor micro-scales of drag/lift fluctuations for the large-diameter particle are almost comparable to their integral time scales, suggesting the absence of disparity in time scales. In the SP case, small-scale drag/lift fluctuations are dominated by small-scale turbulent motions of the time scale ${\sim}d/U_{\infty }$ , typically resulting in Taylor micro-scales of the order of unity. In general, these results indicate that the small-diameter particle typically exhibits force fluctuations of longer duration than those of the large-diameter particle case under investigation.

Table 3.

Integral time scale ( $t_{l}$ ) and temporal Taylor micro-scale ( $t_{\unicode[STIX]{x1D706}}$ ) of drag and lift fluctuations for the LP and SP cases. The time scales are non-dimensionalized using $d/U_{\infty }$ .

To further understand the behaviour of particle-bed forces, the temporal cross-correlation function of the drag and lift fluctuations is computed. This would help to investigate if drag and lift fluctuations are correlated in time and are affected by similar near-bed flow structures. Phase-averaged cross-correlation is calculated by setting $\unicode[STIX]{x1D719}^{\prime }=F_{D}^{\prime }$ , $\unicode[STIX]{x1D713}^{\prime }=F_{L}^{\prime }$ and $\unicode[STIX]{x0394}x=\unicode[STIX]{x0394}y=0$ in (3.2) and is plotted in figure 10. For no shift in time ( $\unicode[STIX]{x0394}t=0$ ), the cross-correlation is very small and the drag and lift fluctuations are weakly correlated. The cross-correlation reaches a maximum value $0.45$ ( $0.48$ ) for separation in time $\unicode[STIX]{x0394}tU_{\infty }/d=0.86$ ( $1.08$ ) for LP (SP) case, indicating that the drag and lift fluctuations, on average, are positively correlated with each other. This suggests that the positive fluctuation in the drag (result of stagnation pressure) corresponds to the simultaneous positive fluctuation in the lift. It is also interesting to note that the time lag of the maximum value of the correlation function approximately matches the Taylor micro-scale of the drag and lift fluctuations (cf. table 3). Furthermore, as shown in figure 10, the cross-correlation takes the minimum value $-0.39$ ( $-0.45$ ) for separation in time $\unicode[STIX]{x0394}tU_{\infty }/d=-0.73$ ( $-0.98$ ) for LP (SP) case. For larger negative separation times, the correlation function for the LP and SP cases increases from a minimum value to reach zero. In comparison to previous unidirectional flow studies (Hofland Reference Hofland2005; Dwivedi Reference Dwivedi2010; Chan-Braun et al. Reference Chan-Braun, Garcia-Villalba and Uhlmann2013; Amir et al. Reference Amir, Nikora and Stewart2014), there exist significant differences in the present case in terms of flow conditions, particle Reynolds number, particle shape and arrangement, methodologies for the flow data measurements, etc. However, it is promising to see that the overall shape and the behaviour of the cross-correlation function, i.e. local maximum (minimum) for positive (negative) separation, remains identical. Also, the maximum values of the cross-correlation function between $0.12$ and $0.55$ were reported in these studies and, interestingly, the present values for the LP and SP cases also fall well within this range. Such behaviour of the cross-correlation of drag and lift fluctuations seems to be very robust and is conjectured to be the effect of a strong influence of the pressure field on the particle-bed forces. This interpretation is in favour of a model approach proposed by Hofland (Reference Hofland2005) based on a convected frozen pressure field.

Figure 10.

Temporal cross-correlation of drag and lift fluctuations plotted as a function of separation in time for the LP (- - - -) and the SP (——) cases. Only a small fraction of the signal in which non-negligible correlation exists is plotted.

3.5 Spectra

To understand the spectral content of drag and lift fluctuations, their frequency spectra $\widehat{R}_{\unicode[STIX]{x1D719}^{\prime }\unicode[STIX]{x1D719}^{\prime }}(f)$ are investigated. Owing to the finite nature of the signal, the modified periodogram method of Welch (Reference Welch1967) with $50\,\%$ overlap and applying a Hamming window to the original signal is employed. Here, spectra are averaged for all the particles and flow cycles and are normalized by the integral of spectra over all the frequencies. Figure 11 shows normalized spectra of drag (solid line) and lift (dashed line) fluctuations for (a) LP and (b) SP cases. It can be seen that, in both flow cases, a reasonably good collapse of spectra for drag and lift fluctuations is observed, with the largest spectral contributions at lower frequencies. Previously, George, Beuther & Arndt (Reference George, Beuther and Arndt1984) have performed spectral analysis of turbulent pressure fluctuations in a homogeneous constant mean-shear flow and reported the presence of $-11/3$ and $-7/3$ scaling in the pressure spectrum when turbulence–mean-shear and turbulence–turbulence interactions dominate, respectively. Interestingly, the present data for oscillatory flow over a rough bed are reminiscent of similar scalings for drag and lift fluctuations.

Figure 11.

Spectra of drag (——) and lift (- - - -) fluctuations for (a) LP and (b) SP cases.

In comparison to the LP case, the $-11/3$ scaling region for the SP case extends a little further up to higher frequencies ( ${\approx}600~\text{Hz}$ ), suggesting that a wider range of scales are dominated by turbulence–mean-shear interactions. On the other hand, as a result of significant flow isotropization in the LP case, it is possible that the turbulence–turbulence interactions play an equally important role exhibiting $-7/3$ scaling for much higher frequencies ( ${\approx}100~\text{Hz}$ ). Recently, Hofland (Reference Hofland2005) and Detert et al. (Reference Detert, Weitbrecht and Jirka2010) in their unidirectional flow studies showed that the turbulent wall pressure spectra exhibit $-7/3$ scaling, whereas Amir et al. (Reference Amir, Nikora and Stewart2014) reported the presence of $-11/3$ scaling region. Contrary to these laboratory studies, where only pressure spectra were analysed, the present results provide spectral energy content of calculated drag and lift fluctuations. Nevertheless, it is interesting to see that the present data for the LP and SP cases typically exhibit similar behaviour as that observed in unidirectional studies. This highlights the important role of pressure, indicating the universal nature of the force fluctuations on the roughness bed.

3.6 Particle forces and turbulent flow-field correlations

The aim of this investigation is to shed light on the role of near-bed flow features that are responsible for inducing unsteady fluctuating forces on the particle bed. As noted earlier, computations or measurements of drag and lift forces on the particle bed are not straightforward. Instead, they are calculated indirectly by assuming a cross-correlation with near-bed velocity fluctuations. In this context, two questions arise: What is the distribution of the near-bed velocity? And, to what spatial extent is it cross-correlated with unsteady forces on the particle bed?

To answer the former, Ghodke et al. (Reference Ghodke, Apte and Urzay2014a ,Reference Ghodke, Skitka and Apte b ) showed the probability density function (p.d.f.) of the near-bed velocity and pressure fluctuations in the region above the particle to be non-Gaussian and to follow fourth-order Gram–Charlier distribution model given by

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle f_{GC4}(\unicode[STIX]{x1D713}^{\prime })=\frac{\text{exp}(-{\unicode[STIX]{x1D713}^{\prime }}^{2}/2)}{\sqrt{2\unicode[STIX]{x03C0}}}\left[1+\frac{S_{\unicode[STIX]{x1D713}^{\prime }}}{3!}({\unicode[STIX]{x1D713}^{\prime }}^{3}-3\unicode[STIX]{x1D713}^{\prime })+\frac{K_{\unicode[STIX]{x1D713}^{\prime }}-3}{4!}({\unicode[STIX]{x1D713}^{\prime }}^{4}-6{\unicode[STIX]{x1D713}^{\prime }}^{2}+3)\right], & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D713}^{\prime }$ is the normalized fluctuating flow variable, and $S_{\unicode[STIX]{x1D713}^{\prime }}$ and $K_{\unicode[STIX]{x1D713}^{\prime }}$ are the skewness and flatness of its distribution, respectively.

To answer the latter, evaluation of the spatial extent of the cross-correlation between forces on the particle bed and the streamwise velocity and pressure fluctuations is carried out in a similar manner as done by Chan-Braun et al. (Reference Chan-Braun, Garcia-Villalba and Uhlmann2013) for unidirectional steady flow. To this end, the three-dimensional correlation function given by (3.2) is computed by taking $\unicode[STIX]{x1D719}^{\prime }$ as the local stress field corresponding to the fluctuating drag/lift force and $\unicode[STIX]{x1D713}^{\prime }$ as the flow variable fluctuations. As per (3.2), the correlation function $R_{\unicode[STIX]{x1D719}^{\prime }\unicode[STIX]{x1D713}^{\prime }}$ is averaged for all the particles and for $15$ flow cycles in each case. For the purpose of comparison between two roughness sizes, $R_{\unicode[STIX]{x1D719}^{\prime }\unicode[STIX]{x1D713}^{\prime }}$ is normalized by using maximum friction velocity $u_{\unicode[STIX]{x1D70F},max}$ and $\unicode[STIX]{x1D70C}u_{\unicode[STIX]{x1D70F},max}^{2}$ as characteristic scales for velocity and pressure fluctuations, respectively, and by using the standard deviation of the respective force fluctuations. As a result, the correlation here is not bounded by unity. The maximum values of the correlations, $|R_{\unicode[STIX]{x1D719}^{\prime }\unicode[STIX]{x1D713}^{\prime }}|_{max}$ , are given in table 4.

In the following, cross-correlations of the drag component of the stress with velocity fluctuations ( $\unicode[STIX]{x1D70F}_{D}^{\prime }-u^{\prime }$ ) and the lift component of the stress with pressure fluctuations ( $\unicode[STIX]{x1D70F}_{L}^{\prime }-p^{\prime }$ ) at the peak phase are presented. This choice is based on the following considerations. In the similar set-up of a closely packed bed studied by Mazzuoli & Vittori (Reference Mazzuoli and Vittori2016), the pressure contribution to the lift force was significantly larger than the viscous contribution. On the other hand, the viscous contribution (which is mainly related to the streamwise velocity) was significant and comparable to the pressure in the case of the drag force. In this context, drag (lift) is cross-correlated with velocity (pressure) fluctuations. Note that drag and lift fluctuations, on average, are positively correlated with a time lag (cf. figure 10) and therefore given correlation fields should suffice for complete description of the flow field.

Table 4.

Maximum amplitude of the correlation between particle force and the flow field, $|R_{\unicode[STIX]{x1D719}^{\prime }\unicode[STIX]{x1D713}^{\prime }}|_{max}$ , in the LP case and the SP case. Here, $|R_{\unicode[STIX]{x1D719}^{\prime }\unicode[STIX]{x1D713}^{\prime }}|_{max}$ is normalized by using maximum friction velocity $u_{\unicode[STIX]{x1D70F},max}$ and $\unicode[STIX]{x1D70C}u_{\unicode[STIX]{x1D70F},max}^{2}$ as characteristic scales for velocity and pressure fluctuations, respectively, and by using standard deviation of the respective force fluctuations.

Figure 12.

Close-up view of isocontours of correlation function $R_{\unicode[STIX]{x1D719}^{\prime }\unicode[STIX]{x1D713}^{\prime }}$ of drag ( $\unicode[STIX]{x1D70F}_{D}^{\prime }$ ) and streamwise velocity ( $u^{\prime }$ ) fluctuations for the LP (a,b) and the SP (c,d) cases plotted at a peak phase. Panels (a,c) show planes at zero spanwise shift ( $\unicode[STIX]{x0394}y=0$ ) and (b,d) show planes at zero streamwise shift ( $\unicode[STIX]{x0394}x=0$ ). Here, correlation is normalized by $u_{\unicode[STIX]{x1D70F},max}$ and standard deviation of drag fluctuations $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D70F}_{D}^{\prime }}$ at a peak phase. Solid (dashed) lines indicate positive (negative) correlation.

Figure 12 shows isocontours of the correlation function of the drag component of the stress and velocity fluctuations for the LP (a,b) and the SP (c,d) cases plotted at a peak phase. Planes at zero spanwise shift ( $\unicode[STIX]{x0394}y=0$ ) and zero streamwise shift ( $\unicode[STIX]{x0394}x=0$ ) are shown in figures 12(a,c) and 12(b,d), respectively. A close-up view of the domain where significant correlations exists is presented. As shown in figure 12(a), drag is positively correlated with velocity fluctuations at around one diameter distance upstream and downstream from the large-diameter particle. The spanwise view shown in figure 12(b) reveals the presence of ‘wing-shaped’ positive correlation between drag and velocity fluctuations. The presence of such a wing shape could be attributed to the complex shielding effects due to the hexagonal packing of the bed that results in disruption of the typical pattern of high-speed/low-speed streaks found in the square arrangement of particles. As shown in figure 12(c,d), a similar shape of the correlation is observed in the case of the small-diameter particle. However, in comparison to the LP case, the spatial extent of the correlation is typically 2–3 times larger in both streamwise and spanwise directions for the SP case. As discussed earlier, this could be attributed to the presence of longer near-bed flow structures found in the SP case. On the contrary, as a result of significant flow isotropization in the LP case, only local flow structures around the particle are responsible for inducing hydrodynamic forces on the particle.

As seen in figure 12(a,b), there exists a region of negative correlation between drag and velocity fluctuations just above the crest of the large-diameter particle; however, in the case of smaller particles, the correlation is positive at this location, as shown in figure 12(c). It is conjectured that such a region of negative correlation in the LP case is a result of local near-bed flow deceleration before the end of the peak phase and is typically less predominant in the high-Reynolds-number SP case, a behaviour also evident from the reduced phase lag between particle forces and outer flow as discussed in § 3.2. These effects are attributed to the oscillatory nature of the flow and were not seen in a unidirectional flow study (Chan-Braun et al. Reference Chan-Braun, Garcia-Villalba and Uhlmann2013).

Figure 13 presents isocontours of correlation function of lift component of the stress and pressure fluctuations for the LP (a,b) and the SP (c,d) cases plotted at a peak phase. The overall shape of the correlation is similar for large- and small-diameter particles. In both the cases, a negative correlation between lift and pressure fluctuations is observed in the region above the particle crest, indicating that negative fluctuation in the pressure results in positive fluctuation in the lift. This region extends up to $1D$ distance upstream and downstream in the streamwise direction from the large and the small particle centre as shown in figure 13(a,c). Similar to figure 12(b,d), the presence of wing-shaped positive correlation between lift and pressure fluctuations for both LP and SP cases is clearly evident in figure 13(b,d). Again, the overall shape of this correlation is comparable for the large- and small-diameter particles, although a little wider for the latter, as shown in figure 13(d). Also, as reported in table 4, it is interesting to see that the maximum values of correlation between particle forces and flow variables are almost of the same order of magnitude for both flow cases.

In summary, the presence of such regions of significant correlations between the forces on particles and the turbulent flow field around them further provides evidence of the direct link between near-bed flow structures and unsteady hydrodynamic forces on the roughness bed. The knowledge of such regions of significant correlations is beneficial and could be used to indirectly calculate the particle forces by measuring the flow variables at such locations.

Figure 13.

Close-up view of isocontours of correlation function $R_{\unicode[STIX]{x1D719}^{\prime }\unicode[STIX]{x1D713}^{\prime }}$ of lift ( $\unicode[STIX]{x1D70F}_{L}^{\prime }$ ) and pressure ( $p^{\prime }$ ) fluctuations for the LP (a,b) and the SP (c,d) cases plotted at a peak phase. Panels (a,c) show planes at zero spanwise shift ( $\unicode[STIX]{x0394}y=0$ ) and (b,d) show planes at zero streamwise shift ( $\unicode[STIX]{x0394}x=0$ ). Here, correlation is normalized by $\unicode[STIX]{x1D70C}u_{\unicode[STIX]{x1D70F},max}^{2}$ and standard deviation of lift fluctuations $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D70F}_{L}^{\prime }}$ at a peak phase. Solid (dashed) lines indicate positive (negative) correlation.

3.7 Particle force distributions

In order to further explore the possibility of influencing probabilistic models for sediment incipient motion, it would be of interest to analyse the p.d.f. of near-bed flow parameters and relate it with the hydrodynamic forces on the particles. Given the high correlation of particle forces with the near-bed flow parameters, it is expected that they follow a similar distribution as that of the velocity and pressure fluctuations. In that case, statistics of forces on the particle bed could be obtained by using the information of measurable near-bed flow variables.

To this end, p.d.f.s of drag and lift fluctuations are constructed for the LP and SP cases. Also, pressure and velocity data are collected in the regions of high correlations as shown in figures 12 and 13, and their higher-order moments such as skewness and flatness are used to build a fourth-order Gram–Charlier distribution given by (3.3). Figure 14 shows p.d.f.s of (a,c) drag and (b,d) lift fluctuations represented by histograms for the (a,b) LP and (c,d) SP cases plotted at a peak phase. The p.d.f.s of (a,c) velocity and (b,d) pressure fluctuations constructed by using a fourth-order Gram–Charlier distribution are also plotted (solid line) in figure 14. In both cases, p.d.f.s of drag (lift) fluctuations, in particular, peakedness and long tails, match remarkably well with fourth-order Gram–Charlier distributions of velocity (pressure) fluctuations. These findings imply that the distribution of hydrodynamic forces on the particle bed is influenced by near-bed flow structures and can be constructed through higher-order flow statistics, further re-emphasizing the need to incorporate the effects of near-bed bursting and sweep phenomena into probabilistic models for onset of erosion.

Figure 14.

P.d.f.s, normalized by standard deviations at a peak phase, of (a,c) drag (velocity) fluctuations and (b,d) lift (pressure) fluctuations represented by histograms (solid line). Panels (a,b) and (c,d) are for LP and SP cases, respectively. The p.d.f. of velocity and pressure fluctuations (solid line) is constructed by using a fourth-order Gram–Charlier distribution (3.3).

3.8 Concept of impulse

Diplas and co-workers (Diplas et al. Reference Diplas, Dancey, Celik, Valyrakis, Greer and Akar2008; Celik et al. Reference Celik, Diplas, Dancey and Valyrakis2010; Valyrakis et al. Reference Valyrakis, Diplas, Dancey, Greer and Celik2010; Celik et al. Reference Celik, Diplas and Dancey2013; Celik, Diplas & Dancey Reference Celik, Diplas and Dancey2014) introduced the concept of impulse indicating the importance of both magnitude and duration of fluctuating hydrodynamic forces on sediment erosion. Although such an approach holds promise in accounting for the intermittent dynamics of turbulent flows, there are certain difficulties in describing such a criterion. For example, detecting the peak events to calculate impulse from the time signal of a force history requires a priori specification of a threshold or a critical value of the force. It is difficult, however, to determine the correct choice of such a threshold that is best suited to predict the incipient motion. The present investigation deals with a fixed particle bed and it is not possible to detect the peak events resulting in particle erosion. Nevertheless, to shed some light on the topic, statistical analyses to provide a brief description of intermittency in fluctuating forces by means of counts and duration of peak events, impulse magnitude and their p.d.f.s for a range of threshold conditions are provided. Such characterization of force intermittency might serve as a useful tool in making a reasonable choice of a threshold criterion in future studies pertaining to this topic.

Consider a sample time history of the drag force ( $F_{D}$ ) on a large-diameter particle as shown in figure 15. In this case, although the period-averaged drag force (dashed line) on the particle is less than the specified threshold (dash-dotted line), an instantaneous force exceeding the threshold for a long duration could still possibly result in particle dislodgement (Diplas et al. Reference Diplas, Dancey, Celik, Valyrakis, Greer and Akar2008). As shown in figure 15, such extreme force events are identified when the instantaneous force value exceeds the threshold and then the durations of such events ( $t_{I_{1}}$ , $t_{I_{2}}$ , $t_{I_{3}}$ , etc.) are computed. Based on this, the impulse on a particle for a given duration is calculated using

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle \text{impulse}=I(\unicode[STIX]{x1D719}>\unicode[STIX]{x1D719}_{thr})=\int _{t_{I}}\unicode[STIX]{x1D719}\,\text{d}t, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D719}_{thr}$ is the specified threshold value of the variable $\unicode[STIX]{x1D719}$ . Implementing this procedure for the LP and SP cases, the number of such peak force events, their duration and impulse magnitude are calculated for various threshold conditions ( $0\leqslant \unicode[STIX]{x1D719}_{thr}/\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D719}}\leqslant 5$ ). Subsequently, averaging is carried out for number of particles and flow cycles. Note that change in the sign of the drag force after flow reversal is taken into consideration while identifying peak force events.

Figure 15.

Representative time history (——) of the drag force ( $F_{D}$ ) in the LP case. Period-averaged drag force (- - - -) along with the specified threshold, $\unicode[STIX]{x1D719}_{thr}/\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D719}}=1$ (— ⋅ —) are also indicated. Local peaks in $F_{D}$ exceeding the threshold are detected as shown. The widths of the local peaks indicate their duration $t_{I}$ .

Figure 16.

Variation with threshold values for (a) LP and (b) SP, of average number of peak force events, $N_{I}$ , experienced by a particle per flow cycle. The lines indicate drag (——) and lift (- - - -).

Figures 16(a) and 16(b) show average number of intervals, $N_{I}$ , of peak force events per flow cycle for the LP and SP cases, respectively. In both flow cases, the number of intervals for the lift force is typically larger than that for the drag force, indicating higher intermittency in the lift. In the case of LP, the number of intervals $N_{I}$ initially increases (more rapidly for the lift) with increase in the threshold, indicating the presence of relatively low-magnitude and high-frequency local peaks. Note that these local peaks are otherwise counted as a single peak for small threshold values under the present peak detection scheme as represented in figure 15. With further increase in the threshold value, $N_{I}$ for drag and lift begins to decline gradually isolating extreme force events. As shown in figure 16(a), peak value $N_{I}$ for the lift is almost twice as high as that for the drag. Therefore, in comparison to the drag, the number of high-frequency peak events is typically more in the case of lift. Moreover, the peak in the magnitude of $N_{I}$ for drag (lift) is seen at $\unicode[STIX]{x1D719}_{thr}/\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D719}}\approx 1~(2)$ . This further shows that, on average, extreme events in the lift force for the large-diameter particle also have higher magnitudes compared to those in the drag.

In the case of SP, this behaviour is more intricate. Consistent with the LP case, $N_{I}$ for the lift in the SP case initially increases rapidly with increase in the threshold, indicating the presence of low-magnitude local peaks; however, $N_{I}$ for the drag decreases slowly, exhibiting a local maximum close to $\unicode[STIX]{x1D719}_{thr}/\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D719}}\approx 0.6$ . Based on figure 16(b), two observations for the SP case can be drawn, as follows. First, similar to the LP case, the peak in the $N_{I}$ for the lift is almost 1.5–2 times higher than that for the drag. And second, both the peaks are observed at around the same threshold value ( $\unicode[STIX]{x1D719}_{thr}/\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D719}}\approx 0.6$ ). The former suggests that the lift data consist of more high-frequency peak events than those found in the drag, whereas the latter indicates that the peaks in the lift and drag force have comparable magnitudes.

Furthermore, all the profiles for $N_{I}$ in both LP and SP cases exhibit a recognizable change in the slope after which the drop in the value of $N_{I}$ is more gradual (e.g. after $\unicode[STIX]{x1D719}_{thr}/\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D719}}\approx 1$ in the SP case). Such flattening of profiles can be interpreted as the presence of well-pronounced isolated peaks in the drag and lift forces that keep getting detected even with the increase in the threshold value.

Figures 17(a) and 17(b) respectively show variation of average duration, $t_{I}$ , of peak force events per flow cycle for the LP and SP cases. In the case of LP, the time duration of lift is typically 1.5–2 times longer than that of the drag force. This corroborates figure 16(a,b), suggesting the presence of more high-frequency peak events in the lift that results in larger time duration compared to drag, albeit, for the small-diameter particle, this difference in drag and lift duration is indeed smaller. In both LP and SP cases, average time duration decreases with increase in the threshold value. Also, as shown in figure 17(a,b), there exists a sudden change in the slope of the profiles. As the threshold value increases, it is evident that the change in $t_{I}$ is more gradual. This is in accordance with figure 16(a,b) and further provides the evidence for the presence of pronounced isolated peaks in drag and lift forces.

Based on this, an average impulse on a particle per flow cycle is calculated using (3.4) and is plotted in figures 18(a) and 18(b) for the LP and SP cases, respectively. Consistent with the earlier observations, the impulse based on the lift force is higher than that based on the drag. With increase in the threshold, the impulse related to the lift for LP (SP) particle is seen to decrease; however, that of the drag almost stays constant up to the threshold value of 1 (0.7). It is possible that the impulse based on the drag is dominated by a few local events that get detected until this threshold value. Again, the profiles of average impulse show very gradual decrease for high threshold values. For the purpose of illustration, the distribution of impulse events at a threshold value $\unicode[STIX]{x1D719}_{thr}/\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D719}}=1$ is plotted in figure 19 for (a) LP and (b) SP case. All distributions are heavily and positively skewed and are well described by a generalized extreme value distribution.

Figure 17.

Variation with threshold values of average time duration of peak force events, $t_{I}U_{\infty }/d$ , for a particle per flow cycle. Line styles and panels are the same as in figure 16.

Figure 18.

Variation with threshold values of average impulse, $I$ , for a particle per flow cycle. Here, $F_{R}=\unicode[STIX]{x1D70C}u_{\unicode[STIX]{x1D70F},max}^{2}A_{s}U_{\infty }/d$ , where $A_{s}$ is the projected area of a single spherical particle given by $\unicode[STIX]{x03C0}d^{2}/4$ . Line styles and panels are the same as in figure 16.

Figure 19.

Distribution of impulse normalized by standard deviations at a threshold $\unicode[STIX]{x1D719}_{thr}/\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D719}}=1$ for (a,c) drag force and (b,d) lift force represented by histograms. Panels (a,b) and (c,d) are for the LP and SP cases, respectively.

Based on these results, it is reasonable to suggest that using a threshold value $\unicode[STIX]{x1D719}_{thr}/\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D719}}>1.5$  (2) for drag and lift in the case of LP (SP) might be sufficient to identify well-pronounced peak events. As noted earlier, the present study considers a roughness bed made up of fixed particles and therefore it is only conjectured that these events could possibly result in particle dislodgement, a subject of future investigation. Nevertheless, these results are typically in agreement with recently reported unidirectional studies (Chan-Braun Reference Chan-Braun2012; Celik et al. Reference Celik, Diplas and Dancey2014). For example, Celik et al. (Reference Celik, Diplas and Dancey2014) used a threshold value based on the pressure, $\unicode[STIX]{x1D719}_{thr}/\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D719}}=2$ for $D^{+}$ values between 330 and 440. The distribution of impulse events reported by Celik et al. (Reference Celik, Diplas, Dancey and Valyrakis2010, Reference Celik, Diplas and Dancey2013, Reference Celik, Diplas and Dancey2014) also showed highly skewed nature similar to that observed in figure 19.

3.9 Summarizing implications of flow oscillations

The oscillatory nature of the near-bed turbulence governs the overall structure of hydrodynamic forces on the particle bed. As expected in the case of wall-bounded oscillatory flows, the near-bed flow depending on the Reynolds number is typically out of sync with the external pressure gradient, resulting in early deceleration of local flow in the vicinity of the particle bed. As a consequence of this, the drag and lift record shows a phase delay with the free-stream velocity. Furthermore, contrary to unidirectional flows, turbulence is present only in some part of the oscillation cycle. In the present investigation, especially for the small-diameter particle case, the sudden occurrence of turbulence followed by flow relaminarization is evident; the former provides favourable conditions for dislodgement of the particles, whereas the latter typically offers stability to the particle bed. Such behaviour of near-bed flow also greatly influences the spatial structure of the forces on the particles, resulting in local patches of positive/negative forces. As a result, strong non-Gaussian behaviour of both fluctuating near-bed flow and force statistics is observed in oscillatory flows.

On the other hand, although the overall structure of turbulent oscillatory flow differs significantly from that of unidirectional flow, there are some similarities in the characteristics of drag/lift forces that appear to be robust. For example, cross-correlation of drag and lift depicts qualitatively similar variation as that observed in other unidirectional studies, albeit with differences in the amplitudes of maxima and minima. In addition, scaling regions observed in force spectra as well as the distribution of impulse events suggest the universal nature of force fluctuations on the particle bed.

3.10 Implications for the sediment incipient motion models

In the present investigation, a rough bed made up of a single layer of hexagonally packed fixed spherical particles is considered. Although using such an idealized flow set-up may limit the scope of application of the present results to specific problems, it still provides us with a detailed understanding of particle–turbulence interactions, which are in turn responsible for sediment incipient motion. The novelty of the present investigation is to provide data, without requiring any ad hoc models, on the temporal variations in the magnitude of drag and lift forces on the sediment bed, the time scales associated with these variations, and their correlations with near-bed turbulence in an oscillatory flow environment. The three-dimensional correlation function between particle forces and near-bed turbulent flow reveals the presence of a significant correlation region, the shape of which remained unchanged between large- and small-diameter particle cases. Velocity and pressure fluctuations measured in such high-correlation regions follow a fourth-order Gram–Charlier distribution, and it is demonstrated that the distribution of hydrodynamic forces on the particle bed could be constructed using such higher-order flow statistics. In other words, it is possible to obtain the particle force distribution by simply measuring the velocity or pressure fluctuations in such high-correlation regions. In this section, a brief summary giving a few examples of how such detailed flow information could be used in developing sediment erosion models is presented.

Computations of sediment-laden flows typically employ simplistic models to study sediment transport, in particular, onset of sediment erosion (Bagnold Reference Bagnold1966; Fredsøe & Deigaard Reference Fredsøe and Deigaard1992; Nelson et al. Reference Nelson, Schmeeckle, Shreve and McLean2000; Nino et al. Reference Nino, Lopez and Garcia2003). Historically, these models are based on the Shields parameter, $\unicode[STIX]{x1D703}=\bar{\unicode[STIX]{x1D70F}}_{w}/(\unicode[STIX]{x1D70C}_{p}-\unicode[STIX]{x1D70C})gd$ (Shields Reference Shields1936), which is the ratio of average bed shear stress ( $\bar{\unicode[STIX]{x1D70F}}_{w}$ ) to the stabilizing settling force due to gravitational acceleration ( $g$ ). Here, $\unicode[STIX]{x1D70C}_{p}$ and $\unicode[STIX]{x1D70C}$ are sediment and fluid densities, respectively, and $d$ is the average particle size. Based on this, rates of sediment pickup ( $P_{k}$ ) could be expressed with a relationship, $P_{k}\propto (\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}_{c})$ . For example, the sediment pickup rate reported by Van Rijn (Reference Van Rijn1984) has been widely used in many sediment erosion studies and is given by

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{P_{k}}{\sqrt{\displaystyle \frac{(\unicode[STIX]{x1D70C}_{p}-\unicode[STIX]{x1D70C})}{\unicode[STIX]{x1D70C}}gd}}=\unicode[STIX]{x1D6FC}\left(\frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}_{c}}{\unicode[STIX]{x1D703}_{c}}\right)^{\unicode[STIX]{x1D6FE}}\left(\frac{(\unicode[STIX]{x1D70C}_{p}-\unicode[STIX]{x1D70C})gd^{3}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D708}^{2}}\right)^{\unicode[STIX]{x1D6FD}}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D703}_{c}$ is the critical Shields parameter, and constants $\unicode[STIX]{x1D6FC}$ , $\unicode[STIX]{x1D6FD}$ and $\unicode[STIX]{x1D6FE}$ are obtained using a fit to experimental data. Assuming the erosion is initiated by lift forces alone, a simple force balance between hydrodynamic lift force and the buoyant weight of the particle yields an expression for the critical Shields parameter and is given as (Chan-Braun et al. Reference Chan-Braun, Garcia-Villalba and Uhlmann2011)

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}_{c}=\frac{2}{3\widehat{C}_{L}}, & \displaystyle\end{eqnarray}$$

where $\widehat{C}_{L}$ is the lift coefficient with a slightly different normalization and is defined as $\widehat{C}_{L}=F_{L}/(\unicode[STIX]{x1D70C}u_{\unicode[STIX]{x1D70F}}^{2}d^{2}\unicode[STIX]{x03C0}/4)$ . Calculation of $\unicode[STIX]{x1D703}_{c}$ therefore requires information of force statistics. Most of the erosion models typically assume a constant value for the lift coefficient. However, the present DNS investigation reveals a large intermittency in the lift forces, resulting in longer tails of force distributions. Such DNS investigation is therefore useful in providing accurate descriptions of force statistics to calculate instantaneous critical values of the Shields parameter. In addition, most of these models are deterministically based on mean values of the bed shear stress, and do not consider particle size, shape, bed arrangements, or multi-scale interactions between turbulent flow and sediment particles. Using DNS results, such as these, however, it is possible to incorporate the effects of particle–turbulence interactions and provide time-dependent wall shear stress information to calculate the instantaneous pickup rate, $P_{k}(t)$ , of sediment grains. For example, it is possible to include higher-order flow statistics to calculate $\unicode[STIX]{x1D703}_{c}$ and lift coefficients. By doing so, it would then be possible to incorporate effects of fluctuations, in addition to mean stress, into an erosion criterion.

On the other hand, stochastic models of incipient motion, based on modelling the probability of particle entrainment, hold promise in accounting for the intermittent dynamics of turbulent flows (Einstein Reference Einstein1950; Mingmin & Qiwei Reference Mingmin and Qiwei1982; Cheng & Chiew Reference Cheng and Chiew1998; Papanicolaou et al. Reference Papanicolaou, Diplas, Evaggelopoulos and Fotopoulos2002; Wu & Lin Reference Wu and Lin2002; Cheng & Law Reference Cheng and Law2003; Wu & Chou Reference Wu and Chou2003; Wu & Kuo-Hsin Reference Wu and Kuo-Hsin2004; Hofland & Battjes Reference Hofland and Battjes2006). Such models often depend on many factors, including particle force statistics, instantaneous bed shear stress, near-bed flow statistics, critical impulse, among others. For example, the probability of particle entrainment ( $P_{E}$ ) can be determined by integrating p.d.f.s of lift, drag or shear stress acting on the particle bed over the range of forces that would result in particle dislodgement and is given as

(3.7) $$\begin{eqnarray}\displaystyle & \displaystyle P_{E}=1-\int _{-F_{cr}}^{F_{cr}}f(F)\,\text{d}F, & \displaystyle\end{eqnarray}$$

where $F$ is the force in question (e.g. lift) and $F_{cr}$ is the critical value of the force at which sediment motion occurs (Papanicolaou et al. Reference Papanicolaou, Diplas, Evaggelopoulos and Fotopoulos2002; Wu & Kuo-Hsin Reference Wu and Kuo-Hsin2004). Calculation of $F$ and $F_{cr}$ requires information about near-bed velocity ( $u_{b}$ ) and force coefficients such as $C_{D}$ or $C_{L}$ . For a well-packed particle bed, such as the present one, movement of the particle typically occurs if the lift force exceeds the submerged weight of the particle (Papanicolaou et al. Reference Papanicolaou, Diplas, Evaggelopoulos and Fotopoulos2002). In such cases, the pickup probability of lifting can be defined as (Wu & Lin Reference Wu and Lin2002)

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle P_{L}=P(F_{L}>W)=P(u_{b}^{2}>B_{L}^{2})=P(u_{b}>B_{L})+P(u_{b}<-B_{L}), & \displaystyle\end{eqnarray}$$

where $F_{L}=C_{L}(\unicode[STIX]{x03C0}d^{2}/4)(\unicode[STIX]{x1D70C}u_{b}^{2}/2)$ is the instantaneous lift force, $W=(\unicode[STIX]{x1D70C}_{p}-\unicode[STIX]{x1D70C})g(\unicode[STIX]{x03C0}d^{3}/6)$ is the submerged weight of the particle and $B_{L}=(4(\unicode[STIX]{x1D70C}_{p}-\unicode[STIX]{x1D70C})gd/3\unicode[STIX]{x1D70C}C_{L})^{0.5}$ denotes the lifting threshold. Based on this, the probability of entrainment in the lifting mode can be calculated as

(3.9) $$\begin{eqnarray}\displaystyle & \displaystyle P_{L}=1-\int _{-B_{L}}^{B_{L}}f(u_{b})\,\text{d}u_{b}. & \displaystyle\end{eqnarray}$$

Calculation of such a probability therefore often requires knowledge of the near-bed velocity distribution along with force statistics. Most of the present models typically assume Gaussian or log-normal distribution for near-bed velocity fluctuations (Papanicolaou et al. Reference Papanicolaou, Diplas, Evaggelopoulos and Fotopoulos2002; Wu & Lin Reference Wu and Lin2002; Cheng & Law Reference Cheng and Law2003; Wu & Chou Reference Wu and Chou2003). Also, as mentioned earlier, computations or measurements of time-varying drag and lift force coefficients are not straightforward. Instead, in many instances, they are assumed to be constant. The present DNS study provides many details, including an accurate description of near-bed flow that follows a fourth-order Gram–Charlier distribution along with time-accurate force statistics. Using such information in stochastic models would therefore facilitate incorporation of higher-order turbulent flow statistics to account for the effects of particle–turbulence interactions, flow oscillations along with near-bed turbulent bursting phenomenon. These results also verify the validity of the fourth-order Gram–Charlier based theoretical model reported by Wu & Kuo-Hsin (Reference Wu and Kuo-Hsin2004). In addition, Diplas et al. (Reference Diplas, Dancey, Celik, Valyrakis, Greer and Akar2008), Celik et al. (Reference Celik, Diplas, Dancey and Valyrakis2010) and Valyrakis et al. (Reference Valyrakis, Diplas, Dancey, Greer and Celik2010) demonstrated the suitability of the impulse concept, providing an additional degree of freedom, time, to predict the particle dislodgement under turbulent flow conditions. In accordance with earlier theoretical studies (Celik et al. Reference Celik, Diplas, Dancey and Valyrakis2010; Valyrakis et al. Reference Valyrakis, Diplas, Dancey, Greer and Celik2010), the present DNS results also show impulse events to be positively skewed and to follow a generalized extreme value distribution; results that could be used to find the probability of particle entrainment as reported by Celik et al. (Reference Celik, Diplas, Dancey and Valyrakis2010). Also, characterization of impulse-generating force events suggests that a threshold value of $\unicode[STIX]{x1D719}_{thr}/\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D719}}\approx 2$ might be sufficient to identify well-pronounced peak events that could play an important role in dislodging the particle. These results are in agreement with recently reported unidirectional studies (Chan-Braun Reference Chan-Braun2012; Celik et al. Reference Celik, Diplas and Dancey2014).

Studying high-Reynolds-number flows, the majority of sediment transport calculations often employ a large-eddy simulation (LES) approach. Such LES calculations of sediment erosion with an Eulerian two-fluid type model typically require a boundary condition for the flux of the sediment and is given by

(3.10) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x1D708}_{T}}{Pr_{T}}\frac{\unicode[STIX]{x2202}\overline{C}}{\unicode[STIX]{x2202}z}=-P_{k}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D708}_{T}$ is the eddy viscosity, $Pr_{T}$ is the turbulent Prandtl number, $\overline{C}$ is the grid-filtered sediment concentration and $z$ represents the wall-normal direction. Using time-accurate flow information from DNS investigation, it is possible to specify turbulent flow properties along with an instantaneous sediment pickup rate that is also non-uniform and varies spatially, $P_{k}(t,x,y)$ .

The present DNS data could therefore be useful in many such ways to provide inputs to existing sediment erosion models. Such detailed information would possibly allow us to devise a refined erosion criterion for oscillatory flows. In order to study actual erosion events, however, a DNS investigation with moving particles is needed and is a subject of future investigation.