1 Introduction
Turbulent flows over rough surfaces have been the subject of rigorous investigation owing to their increasing technological interests in engineering and geophysical applications (Raupach & Thom Reference Raupach and Thom1981; Raupach, Antonia & Rajagopalan Reference Raupach, Antonia and Rajagopalan1991; Jiménez Reference Jiménez2004). There has been considerable work done on stationary, unidirectional turbulent flows over rough surfaces clearly identifying the effects of roughness elements on near-bed turbulence as increased bed shear stress, reduction in near-bed anisotropy, and marked changes in turbulent energy transport mechanisms (Krogstad, Antonia & Browne Reference Krogstad, Antonia and Browne1992; Krogstad & Antonia Reference Krogstad and Antonia1994; Krogstad et al. Reference Krogstad, Andersson, Bakken and Ashrafin2005; Ikeda & Durbin Reference Ikeda and Durbin2007; Chan-Braun, Garcia-Villalba & Uhlmann Reference Chan-Braun, Garcia-Villalba and Uhlmann2011; Kempe, Vowinckel & Frőhlich Reference Kempe, Vowinckel and Frőhlich2014; Yuan & Piomelli Reference Yuan and Piomelli2014b ). Presence of roughness was also shown to greatly influence near-bed sweep and ejection motions (Krogstad et al. Reference Krogstad, Antonia and Browne1992, Reference Krogstad, Andersson, Bakken and Ashrafin2005). In many practical applications, such as sediment transport in coastal flows, blood flow circulation, intake of a combustion chamber and so on, the flow conditions are mainly unsteady/non-stationary in nature. However, such detailed findings pertaining to turbulence modulation by rough-bed under oscillatory flow conditions have not been well documented. In order to improve the understanding of the mechanisms that lead to the onset of sediment erosion under coastal conditions, detailed knowledge of sediment-bed–turbulence interactions is necessary. Impact of near-bed sweep/ejection events on the onset of sediment motion can be significant and should be accounted for in the development of probabilistic models predicting sediment erosion. These models often assume Gaussian or log-normal distribution for near-bed velocity fluctuations (Einstein Reference Einstein1950; Papanicolaou et al. Reference Papanicolaou, Diplas, Evaggelopoulos and Fotopoulos2002). It is hypothesized that the near-bed flow dynamics, roughness effects, and flow through porous regions lead to non-Gaussian velocity fluctuations, and higher-order turbulence statistics, such as kurtosis and skewness, are necessary for better predictions (Wu & Kuo-Hsin Reference Wu and Kuo-Hsin2004).
Present understanding of oscillatory boundary layer over rough surfaces is mainly based on experimental investigations of Jonsson & Carlsen (Reference Jonsson and Carlsen1976), Keiller & Sleath (Reference Keiller and Sleath1976), Sleath (Reference Sleath1987) and Jensen, Sumer & Fredsøe (Reference Jensen, Sumer and Fredsøe1989), including more recent studies of Chen et al. (Reference Chen, Chen, Tang, Stansby and Li2007), Dixen et al. (Reference Dixen, Hatipoglu, Sumer and Fredsøe2008), Corvaro et al. (Reference Corvaro, Miozzi, Postacchini, Mancinelli and Brocchini2014) and Mujal-Colilles et al. (Reference Mujal-Colilles, Mier, Christensen, Bateman and Garcia2014). Keiller & Sleath (Reference Keiller and Sleath1976) first reported velocity measurements close to the rough bed sinusoidally oscillating in its own plane. Jonsson & Carlsen (Reference Jonsson and Carlsen1976) later reported velocity measurements, but not turbulence, for fully turbulent oscillatory flow over a rough bed. However, most detailed experimental studies of sinusoidally oscillating flow over rough beds are those of Sleath (Reference Sleath1987) and Jensen et al. (Reference Jensen, Sumer and Fredsøe1989). Sleath (Reference Sleath1987) presented the measurements for velocity and near-bed turbulent intensities for moderately high Reynolds numbers, whereas Jensen et al. (Reference Jensen, Sumer and Fredsøe1989) further expanded the research by carrying out experiments at higher Reynolds numbers and wider range of particle diameters. In essence, all of these studies reported a common conclusion that the presence of roughness elements significantly affects near-bed turbulence by increasing turbulent intensities and Reynolds stresses. Beyond this, however, due to experimental limitations these studies did not present detailed analysis of near-bed turbulence structure. On the other hand, owing to the enormous computational cost in resolving individual roughness elements and complexities in handling flow unsteadiness, performing detailed numerical investigations of these flows at high Reynolds numbers is still a challenge and very few numerical investigations have been reported so far (Fornarelli & Vittori Reference Fornarelli and Vittori2009; Ding & Zhang Reference Ding and Zhang2010). However, none of these studies presented detailed characterization of near-bed turbulent flow dynamics.
This paper is therefore oriented towards fundamental understanding of roughness–turbulence interactions in symmetric, sinusoidally oscillating turbulent flow over a range of Reynolds numbers. To the best of the authors’ knowledge, no such detailed characterization of near-bed turbulence in oscillatory flows over rough surfaces by means of DNS has been reported to date. The outline of the paper is as follows: the computational set-up and methodology are presented in § 2, simulation results are presented in § 3, followed by main conclusions in § 4.
2 Methodology
2.1 Computational set-up and simulation parameters
The computational domain, shown in figure 1, consists of a doubly periodic box (in
$x$
and
$y$
directions) with a smooth no-slip wall at
$z=0$
and a slip wall at
$z=30{\it\delta}$
, where
${\it\delta}=\sqrt{2{\it\nu}/{\it\omega}}$
is the Stokes-layer thickness,
${\it\nu}$
is the kinematic viscosity,
${\it\omega}=2{\rm\pi}/T$
is the oscillation frequency and
$T$
is the period of the wave. A rough bed made up of a hexagonal pack of fixed spherical roughness elements is considered corresponding to the experimental configuration of Keiller & Sleath (Reference Keiller and Sleath1976). As shown in figure 1,
$x$
,
$y$
and
$z$
are, respectively, the streamwise, spanwise and wall-normal directions and
$u$
,
$v$
and
$w$
are the velocity components in those directions. The in-plane domain length is
$24.1{\it\delta}$
and
$13.9{\it\delta}$
in the streamwise and spanwise directions, respectively. The roughness element diameter normalized with the Stokes-layer thickness is
$D=d/{\it\delta}=6.95$
, equivalent to the gravel-type roughness (Sleath Reference Sleath1987).
(a) Close-up view of the computational domain and boundary conditions. (b) A harmonic pressure forcing is imposed that results in a streamwise velocity component
$U_{\infty }\sin ({\it\omega}t)$
far from the bed.
Computational parameters and grid resolution details. Here
$k_{s}^{+}$
is computed using maximum value of the friction velocity in a flow cycle. The present study assumes low Froude numbers.
The dimensionless parameters of the problem to completely characterize the flow are outlined in table 1, and consist of the Reynolds number based on the Stokes-layer thickness,
$\mathit{Re}_{{\it\delta}}=U_{\infty }{\it\delta}/{\it\nu}$
, and the effective roughness Reynolds number,
$k_{s}^{+}=u_{{\it\tau},max}k_{s}/{\it\nu}$
. Here,
$U_{\infty }$
is the amplitude of the free-stream velocity,
$k_{s}$
is the Nikuradse roughness size and
$u_{{\it\tau},max}$
is the maximum friction velocity in a flow cycle. The corresponding particle Reynolds number,
$\mathit{Re}_{D}=U_{\infty }D/{\it\nu}$
and Reynolds number based on the wave amplitude,
$\mathit{Re}_{a}=U_{\infty }a/{\it\nu}$
are also given in table 1. Here,
$a$
is the amplitude of wave oscillation given by
$a=U_{\infty }/{\it\omega}$
. Figure 2 maps the different cases on a regime diagram (
$a/k_{s}-\mathit{Re}_{a}$
map) for oscillatory flow over roughness elements showing the laminar, transitional and rough turbulent regimes. The cases studied in this work span a broad range of Reynolds numbers, starting from transitional, where incipient vortex formation around the spheres leads to strong vertical velocities at
$\mathit{Re}_{a}=4510$
, up to a fully turbulent, very rough regime,
$\mathit{Re}_{a}=80\,000$
.
Delineation of different flow regimes for oscillatory flow illustrating previous studies (Jonsson & Carlsen Reference Jonsson and Carlsen1976; Kemp & Simons Reference Kemp and Simons1982; Sleath Reference Sleath1987; Jensen et al.
Reference Jensen, Sumer and Fredsøe1989; Krstic & Fernando Reference Krstic and Fernando2001; Chen et al.
Reference Chen, Chen, Tang, Stansby and Li2007; Dixen et al.
Reference Dixen, Hatipoglu, Sumer and Fredsøe2008; van der A et al.
Reference van der A, O’Donoghue, Davies and Ribberink2011) and present cases. Symbols represent: ●,
$\mathit{Re}_{{\it\delta}}=95$
; ♦,
$\mathit{Re}_{{\it\delta}}=150$
; ▪,
$\mathit{Re}_{{\it\delta}}=200$
; ▴,
$\mathit{Re}_{{\it\delta}}=400$
.
Table 1 tabulates the parameters and spatial resolutions used in this study. Uniform grids, although not cubic, are used in the region surrounding the roughness bed (giving roughly around 160–300 points in the wall-normal direction per particle diameter) and grids are stretched in the wall-normal direction (above
$z=10{\it\delta}$
) using a hyperbolic tangent function. In all of the cases,
${\rm\Delta}x^{+}$
and
${\rm\Delta}y^{+}$
are less than 5 and 3, respectively, and
${\rm\Delta}z^{+}<1$
in the near-bed region (up to
$z=10{\it\delta}$
). Around 15 flow cycles per case are computed to obtain statistical convergence.
A systematic grid refinement study was performed to obtain optimum grid spacing, confirming a nearly second-order accuracy of the solver as shown in Ghodke, Skitka & Apte (Reference Ghodke, Skitka and Apte2014b ). The present grid resolution was chosen as there were no discernible differences in the flow statistics of the finest and current grid (see the Appendix). In addition, energy spectra were analysed for different grid resolutions keeping the Courant number (CFL) close to 0.5 for time-accurate calculations. The streamwise and spanwise energy spectra (not shown) illustrated the adequacy of the computational grid to resolve the smallest scales of turbulence. The domain size is verified to be sufficiently large by calculating two-point spatial velocity correlation functions, in both streamwise and spanwise directions, by doubling the domain. The correlations showed no effect of domain size, and confirm the adequacy of the domain size to accommodate the largest turbulent structures. The vertical size of the domain is chosen such that all turbulent statistics decay to zero at almost half the vertical height. Increasing the vertical extent of the domain showed no effects on the near-bed flow statistics (see the Appendix). Similar computational domain sizes were used in the previous work by Fornarelli & Vittori (Reference Fornarelli and Vittori2009) and Ding & Zhang (Reference Ding and Zhang2010).
2.2 Validation study
The numerical solver used in this work is based on a co-located grid, fractional time-stepping, finite volume algorithm based on energy conserving principles and has been developed and validated thoroughly for a range of canonical and complex turbulent flows (Apte et al. Reference Apte, Mahesh, Moin and Oefelein2003; Apte, Mahesh & Moin Reference Apte, Mahesh and Moin2009a ; Apte et al. Reference Apte, Mahesh, Moin and Gorokhovski2009b ; Mahesh et al. Reference Mahesh, Constantinescu, Apte, Iaccarino, Ham and Moin2006; Moin & Apte Reference Moin and Apte2006). This solver was further extended to incorporate the fictitious domain method to handle arbitrary shaped immersed objects without requiring the need for body-fitted grids (Apte, Martin & Patankar Reference Apte, Martin and Patankar2008; Apte & Finn Reference Apte and Finn2013). The fictitious domain approach facilitates the solution of freely moving particles for a wide range of fluid–particle density ratios, although the particles are held fixed in the present study. The solver is fully validated for a range of test cases including flow over a cylinder and sphere for different Reynolds numbers, flow over touching spheres at different orientations, flow developed by an oscillating cylinder, laminar and turbulent flow through porous media, among others. The details of the algorithm as well as very detailed verification and validation studies have been published elsewhere (Apte et al. Reference Apte, Martin and Patankar2008; Apte & Finn Reference Apte and Finn2013). In addition, the solver was also used to perform direct one-to-one comparison with a body-fitted solver with known second-order accuracy for steady inertial, unsteady inertial and turbulent flow through porous media (Finn & Apte Reference Finn and Apte2013). The details of this comparison focused on issues such as grid resolution needed near the wall, issues related to touching spheres and contact points, quality of solution compared with the body-fitted solver, comparison with experimental work of Suekane, Yokouchi & Hirai (Reference Suekane, Yokouchi and Hirai2003) on flow through a packed array of half spheres and establishing the grid resolution requirement per particle diameter for predicting the unsteady flow-field in packed porous beds within 5 % uncertainty based on the grid convergence index (GCI). Turbulent flow at a pore Reynolds number of 600 was also computed in the same configuration and compared with the body-fitted approach to obtain good predictive capability of the present fictitious domain solver. This same solver was also used to study turbulent flow in a randomly packed bed of 51 spheres capturing complex porescale flow features in agreement with published data.
In addition, in order to verify and validate the solver for the specific case under investigation, oscillatory flow over a particle bed corresponding to
$\mathit{Re}_{{\it\delta}}=95$
case is validated against the experimental data of Keiller & Sleath (Reference Keiller and Sleath1976). Figure 3 shows variation of (a) normalized peak fluid-frame resultant velocity magnitude
$U_{R}$
, and (b) phase at which this peak velocity is recorded, plotted at various heights above the roughness crest for
$\mathit{Re}_{{\it\delta}}=95$
. Data from the experimental work of Keiller & Sleath (Reference Keiller and Sleath1976) along with recent simulations of Fornarelli & Vittori (Reference Fornarelli and Vittori2009) and Ding & Zhang (Reference Ding and Zhang2010) are also plotted. Existing literature pertaining to this topic frequently uses an alternative frame of reference, consistent with that of the classic experiments conducted by Keiller & Sleath (Reference Keiller and Sleath1976), which features an oscillating particle bed in a fixed fluid frame. Figure 3 is presented in the fluid frame to be consistent with this convention. As seen in figure 3, current simulations data are generally in a better agreement with the experimental work of Keiller & Sleath (Reference Keiller and Sleath1976), suggesting the present choice of parameters is able to provide accurate descriptions of the oscillatory flow phenomena. Note that, the vertical trends in both the plots starting at around
$(z-z_{c})/{\it\delta}=0.5$
are artifacts of the near-bed incipient ejection motions resulting in strong vertical velocity and are well captured in the current simulations.
For
$\mathit{Re}_{{\it\delta}}=95$
, variation of (a) normalized peak fluid frame resultant velocity magnitude
$U_{R}$
, and (b) phase at which the peak velocity is recorded, plotted at various heights above the roughness crest location
$z_{c}$
normalized by the Stokes-layer thickness
${\it\delta}$
. Here
$U_{R}$
is the resultant velocity magnitude of streamwise and vertical velocity components. Symbols represent: ●, (Keiller & Sleath Reference Keiller and Sleath1976); ♦, (Ding & Zhang Reference Ding and Zhang2010); ▪, (Fornarelli & Vittori Reference Fornarelli and Vittori2009); ▴, present data.
3 Results and discussion
The majority of the results presented here are for
$\mathit{Re}_{{\it\delta}}=400$
, which is in the fully turbulent, very rough regime. Similar trends are observed for other turbulent flow conditions
$\mathit{Re}_{{\it\delta}}=150$
and 200. Following Mignot, Bartheleemy & Hurther (Reference Mignot, Bartheleemy and Hurther2009), a double-averaging procedure is employed for the analysis described in the following that decomposes a flow quantity,
${\it\phi}$
, into double average
$\langle \overline{{\it\phi}}\rangle$
(where overbar and brackets respectively denote phase and homogeneous spatial averages), the spatial disturbance of the phase average
$\widetilde{{\it\phi}}$
and the turbulent fluctuation
${\it\phi}^{\prime }$
,
$$\begin{eqnarray}{\it\phi}(x,y,z,{\it\omega}t)=\langle \overline{{\it\phi}}\rangle (z,{\it\omega}t)+\widetilde{{\it\phi}}(x,y,z,{\it\omega}t)+{\it\phi}^{\prime }(x,y,z,{\it\omega}t).\end{eqnarray}$$
3.1 The mean flow field and Reynolds stresses
Figure 4 shows profiles of double-averaged streamwise velocity in wall units,
$u^{+}=\langle \overline{u}\rangle /u_{{\it\tau}}$
, plotted against
$A(z-z_{b})/k_{s}$
in a semi-log plot; where
$z_{b}$
is the zero-displacement plane (Dixen et al.
Reference Dixen, Hatipoglu, Sumer and Fredsøe2008). The logarithmic law for flow over a rough bed can be written as
$$\begin{eqnarray}\langle \overline{u}(z,{\it\omega}t)\rangle =\frac{u_{{\it\tau}}}{{\it\kappa}}\ln \left(\frac{A(z-z_{b})}{k_{s}}\right),\end{eqnarray}$$
with the von Karman constant
${\it\kappa}=0.41$
and
$A=30$
as given by (Raupach & Thom Reference Raupach and Thom1981). As discussed in Dixen et al. (Reference Dixen, Hatipoglu, Sumer and Fredsøe2008), the zero-displacement plane
$z_{b}$
and Nikuradse’s equivalent roughness
$k_{s}$
were determined by fitting the double-averaged velocity profile to the log-law given by (3.2). The values of
$z_{b}=0.7D$
and
$k_{s}=2d$
, are in a good agreement with the literature (Sleath Reference Sleath1987; Jensen et al.
Reference Jensen, Sumer and Fredsøe1989; Dixen et al.
Reference Dixen, Hatipoglu, Sumer and Fredsøe2008). As seen in figure 4, a significant portion of the velocity curve that follows the log-law given by (3.2) is present from early-acceleration until mid-deceleration phases, i.e. from
$2{\rm\pi}/10$
until
$7{\rm\pi}/10$
. This is due to the fact that late deceleration phases are characterized by the absence of near-bed turbulence production as discussed later; during this period, a new laminar boundary layer starts to develop in the reverse direction and continues to grow until early acceleration phases after which local friction velocity is large enough to trigger the turbulence. For
${\it\omega}t=8{\rm\pi}/10$
, velocity close to the bed is therefore negative due to the near-bed flow reversal. Also, it is interesting to note that some portion of the velocity curve at this phase seem to follow log-law, although with the larger value of the slope, resulting in the departure from universal value of the von Karman constant
${\it\kappa}=0.41$
.
Phase evolution of double-averaged streamwise velocity profile normalized by
$u_{{\it\tau}}$
for
$\mathit{Re}_{{\it\delta}}=400$
at: (a)
$3{\rm\pi}/10$
, (b)
$4{\rm\pi}/10$
, (c)
$5{\rm\pi}/10$
, (d)
$6{\rm\pi}/10$
, (e)
$7{\rm\pi}/10$
, (f)
$8{\rm\pi}/10$
. The dashed line represents the log-law.
Double-averaged profiles of the Reynolds stresses and TKE normalized by
$U_{\infty }^{2}$
are plotted in figure 5. From early-acceleration until late-deceleration phases (
${\it\omega}t=2{\rm\pi}/10$
up to
$8{\rm\pi}/10$
), presence of near-bed turbulence is evident as seen in figure 5. Reynolds stresses and TKE increase away from the effective bed location
$(z-z_{b})/D=0$
and peak close to the roughness crest level
$(z-z_{b})/D=0.3$
, and decay to zero further in the outer region above
$(z-z_{b})/D>1$
. The streamwise component of the Reynolds stress
$\langle \overline{u^{\prime 2}}\rangle$
contributes most to the TKE and peaks at the roughness crest location. On the other hand, as a consequence of the near-bed turbulent motions, the Reynolds shear stress
$\langle \overline{u^{\prime }w^{\prime }}\rangle$
peaks slightly above the crest level.
Phase evolution of double-averaged Reynolds stress and TKE for
$\mathit{Re}_{{\it\delta}}=400$
at: (a)
$3{\rm\pi}/10$
, (b)
$4{\rm\pi}/10$
, (c)
$5{\rm\pi}/10$
, (d)
$6{\rm\pi}/10$
, (e)
$7{\rm\pi}/10$
, (f)
$9{\rm\pi}/10$
. Symbols represent: ●, TKE; ▪, streamwise Reynolds stress
$\langle \overline{u^{\prime 2}}\rangle$
; ♦, spanwise Reynolds stress
$\langle \overline{v^{\prime 2}}\rangle$
; ▸, wall-normal Reynolds stress
$\langle \overline{w^{\prime 2}}\rangle$
; ▴, Reynolds shear stress
$\langle \overline{u^{\prime }w^{\prime }}\rangle$
. The dash-dotted line at
$(z-z_{b})/D=0.3$
shows the roughness crest level.
Phase variation of the friction factor
$f_{w}$
. Symbols represent: ●,
$\mathit{Re}_{{\it\delta}}=95$
; ♦,
$\mathit{Re}_{{\it\delta}}=150$
; ▪,
$\mathit{Re}_{{\it\delta}}=200$
; ▴,
$\mathit{Re}_{{\it\delta}}=400$
.
Variation with
$a/k_{s}$
of (a) normalized boundary layer thickness
${\it\delta}^{\prime }/k_{s}$
; (b) maximum wave friction factor, given as,
$f_{w,max}=2(u_{{\it\tau},max}/U_{\infty })^{2}$
(Bagnold Reference Bagnold1946; Kamphuis Reference Kamphuis1975); (c) peak period-averaged streamwise and (d) wall-normal turbulence intensity. Symbols represent: ●,
$\mathit{Re}_{{\it\delta}}=95$
; ♦,
$\mathit{Re}_{{\it\delta}}=150$
; ▪,
$\mathit{Re}_{{\it\delta}}=200$
; ▴,
$\mathit{Re}_{{\it\delta}}=400$
. The proposed relation (dotted line) between
${\it\delta}^{\prime }/k_{s}$
and
$a/k_{s}$
in (a) is based on Sleath (Reference Sleath1987) and is slightly modified to accommodate lower
$a/k_{s}$
values from Dixen et al. (Reference Dixen, Hatipoglu, Sumer and Fredsøe2008).
For
$\mathit{Re}_{{\it\delta}}=400$
at peak phase
$5{\rm\pi}/10$
, variation of (a) production and dissipation terms where symbols represent ●,
$P_{s}$
; ▪,
$P_{m}$
; ♦,
$P_{w}$
; ▸, viscous dissipation; and (b) transport terms where symbols represent ▸, turbulent transport; ▴, viscous transport; ▪, pressure transport. Bed-induced turbulent transport is negligible and is not plotted. The dash-dotted line at
$(z-z_{b})/D=0.3$
shows the roughness crest level. All of the terms are normalized by
$u_{{\it\tau},max}^{4}/{\it\nu}$
. (c) Schematic representation of the energy transfer mechanisms. The dashed arrow path is not studied in this work.
Phase evolution of the structure parameter,
$a_{1}=\langle \overline{u^{\prime }w^{\prime }}\rangle /2k$
for
$\mathit{Re}_{{\it\delta}}=400$
at: (a)
$2{\rm\pi}/10$
, (b)
$3{\rm\pi}/10$
, (c)
$4{\rm\pi}/10$
, (d)
$5{\rm\pi}/10$
, (e)
$6{\rm\pi}/10$
, (f)
$7{\rm\pi}/10$
, (g)
$8{\rm\pi}/10$
, (h)
$9{\rm\pi}/10$
. Here,
$k$
is the TKE. The dash-dot and dashed lines show the roughness crest level at
$(z-z_{b})/D=0.3$
and the value
$a_{1}=0.15$
, respectively.
To validate the numerical predictions and to further verify the adequacy of the computational domain, grid resolution and rigid body treatment, the boundary layer thickness,
${\it\delta}^{\prime }$
, and maximum friction factor,
$f_{w,max}=2(u_{{\it\tau},max}/U_{\infty })^{2}$
, were calculated. Here, boundary layer thickness
${\it\delta}^{\prime }$
is defined as the elevation of maximum overshoot in the streamwise velocity at a peak phase (Jensen et al.
Reference Jensen, Sumer and Fredsøe1989) measured from
$z_{b}$
. Following Yuan & Piomelli (Reference Yuan and Piomelli2014a
), the friction velocity
$u_{{\it\tau}}$
is given as
$$\begin{eqnarray}u_{{\it\tau}}^{2}=-\overline{F_{t}}/{\it\rho}A,\end{eqnarray}$$
where
$A$
is the horizontal area and
$F_{t}$
is the total streamwise drag acting on the roughness elements along with that on the smooth wall and is computed explicitly in the solver by integrating viscous and pressure forces in the flow direction. Alternatively,
$u_{{\it\tau}}$
can also be computed by fitting the mean velocity profile to the log-law as briefly described by Dixen et al. (Reference Dixen, Hatipoglu, Sumer and Fredsøe2008); the latter yielding approximately the same values as those obtained by the former method for the flow under consideration. Figure 6 shows phase variation of the friction velocity for range of Reynolds numbers in terms of dimensionless friction factor defined as
$f_{w}=2(u_{{\it\tau}}/U_{\infty })^{2}$
. Consistent with Jensen et al. (Reference Jensen, Sumer and Fredsøe1989), the evolution of friction factor and free-stream velocity are not in phase. However, as a result of increased momentum transfer due to turbulence, this phase delay decreases with the increase in Reynolds number.
Figure 7(a) and (b), respectively, show variation of the normalized boundary layer thickness,
${\it\delta}^{\prime }/k_{s}$
, and maximum friction factor,
$f_{w,max}$
with the roughness parameter
$a/k_{s}$
for various Reynolds numbers. The present predictions show very good match against the experimental data (Jonsson & Carlsen Reference Jonsson and Carlsen1976; Sleath Reference Sleath1987; Jensen et al.
Reference Jensen, Sumer and Fredsøe1989; Dixen et al.
Reference Dixen, Hatipoglu, Sumer and Fredsøe2008). With the increase in Reynolds number, the roughness parameter
$a/k_{s}$
also increases, resulting in turbulent and rough bed conditions. Period-averaged turbulence intensities for various values of
$a/k_{s}$
corresponding to different values of Reynolds numbers are also plotted in figure 7(c,d) and compared with Sleath (Reference Sleath1987). The numerical data tends to collapse onto a unique curve with increase in Reynolds number, confirming Sleath’s observation. It can also be seen that, increase in
$a/k_{s}$
results in decrease in streamwise intensity,
$\sqrt{\overline{u^{\prime 2}}}$
, and increase in the wall-normal intensity,
$\sqrt{\overline{w^{\prime 2}}}$
. If compared with the fully turbulent, smooth-wall case from Jensen et al. (Reference Jensen, Sumer and Fredsøe1989) (test 10 for
$\mathit{Re}_{{\it\delta}}\approx 3300$
, not shown here), it is observed that for the present rough-bed cases, the ratio of wall-normal and spanwise turbulent intensities to streamwise component is higher than that in the smooth-wall case of Jensen et al. (Reference Jensen, Sumer and Fredsøe1989). These results indicate the tendency of the roughness elements to redistribute the energy from streamwise to spanwise and wall-normal components, more so at the peak phase and, therefore, decrease the overall large-scale near-bed anisotropy as also observed by Sleath (Reference Sleath1987) in oscillatory flows.
3.2 TKE budget
The TKE budget is studied to analyse the role and relative importance of different terms involved in energy transport mechanisms. Following Mignot et al. (Reference Mignot, Bartheleemy and Hurther2009), double-averaged TKE budget equation for flow over roughness with streamwise–spanwise homogeneity is given as
$$\begin{eqnarray}\displaystyle \frac{\partial \langle \overline{u_{i}^{\prime }u_{i}^{\prime }}\rangle /2}{\partial t} & = & \displaystyle \underbrace{-\langle \overline{u_{i}^{\prime }w^{\prime }}\rangle \frac{\partial \langle \overline{u_{i}}\rangle }{\partial z}}_{P_{s}}+\left[\underbrace{-\!\left\langle \widetilde{u_{i}^{\prime }u_{j}^{\prime }}\frac{\partial \widetilde{u}_{i}}{\partial x_{j}}\right\rangle }_{P_{w}}-\underbrace{\langle \overline{u_{i}^{\prime }u_{j}^{\prime }}\rangle \left\langle \frac{\partial \widetilde{u_{i}}}{\partial x_{j}}\right\rangle }_{P_{m}}\right]\nonumber\\ \displaystyle & & \displaystyle -\,\frac{\partial }{\partial z}\left[\underbrace{\langle \overline{u_{i}^{\prime }u_{i}^{\prime }w^{\prime }}\rangle /2}_{F_{k}}+\underbrace{\langle \widetilde{u_{i}^{\prime }u_{i}^{\prime }}\widetilde{w}\rangle /2}_{F_{w}}\right]-\frac{1}{{\it\rho}}\frac{\partial }{\partial z}\langle \overline{p^{\prime }w^{\prime }}\rangle +{\it\nu}\frac{\partial ^{2}}{\partial z^{2}}\langle \overline{u_{i}^{\prime }u_{i}^{\prime }}\rangle /2-\langle \overline{{\it\epsilon}}\rangle \qquad \quad\end{eqnarray}$$
where the eight terms on the right-hand side are: the shear production term,
$P_{s}$
, that represents the work of the double-averaged velocity against the double-averaged shear; the wake production term,
$P_{w}$
, that is the work of wake-induced velocity fluctuations against the bed-induced shear;
$P_{m}$
, that is the work of the bed-induced velocity fluctuations against double-averaged shear;
$\partial F_{k}/\partial z$
, that is turbulent transport, whereas
$\partial F_{w}/\partial z$
is the bed-induced turbulent transport followed by the sixth term of pressure transport and the seventh term of viscous transport of TKE. The last term on the right represents viscous dissipation. Note that the terms
$P_{m}$
,
$P_{w}$
and
$F_{w}$
arise as a result of spatial heterogeneity at the roughness element length scale. Figure 8(a) and (b) respectively show variation of the production/dissipation terms and TKE transport terms.
As seen in figure 8(a), a peak in the shear production,
$P_{s}$
, is seen just below the crest of the particle at
$0.98D$
, a result of the formation of shear layers near the crest. The dissipation peak is observed slightly below the production peak at
$0.96D$
. The ratio of peak production to the dissipation rate in the cycle is around 1.9. Variation of transport terms plotted in figure 8(b) shows that the excess energy generated in the vicinity of the crest of the roughness elements is transported away mostly by turbulent diffusion and pressure transport, along with a small amount of viscous diffusion. Also, pressure transport is significant mostly below the crest of the particle. These findings are consistent with the existing literature on unidirectional flows (Raupach et al.
Reference Raupach, Antonia and Rajagopalan1991; Finnigan Reference Finnigan2000; Mignot et al.
Reference Mignot, Bartheleemy and Hurther2009; Yuan & Piomelli Reference Yuan and Piomelli2014b
).
Contrary to the earlier studies of unidirectional flows; however, bed-induced production terms,
$P_{m}$
and
$P_{w}$
, are far from negligible and peak at around same location as
$P_{s}$
. Figure 8(c) shows important pathways and energy transfer mechanisms observed in the present work. It is interesting to note that the wake production term,
$P_{w}$
, is negative close to the roughness crest (see II in figure 8
a), indicating the conversion of TKE to the wake kinetic energy (WKE) which is given as
$\langle \widetilde{u_{i}}\widetilde{u_{i}}\rangle /2$
. At this location, streamwise component of the Reynolds stress,
$\langle \overline{u^{\prime 2}}\rangle$
, contributes most to the TKE and is typically associated with near-bed turbulent structures of length scales larger than the roughness scale. The work of these large-scale structures against the pressure drag results in effective transfer of energy from TKE to WKE, therefore resulting in a negative peak in the production term,
$P_{w}$
. On the other hand, a secondary positive peak in
$P_{w}$
is seen just above the effective bed location at
$(z-z_{b})/D\approx 0.1$
(see I in figure 8
a), indicating conversion of WKE to TKE. At this location, all three normal components of the Reynolds stress are comparable and are associated with turbulent scales smaller than the roughness size; therefore, resulting in conversion of WKE to TKE giving a positive peak in
$P_{w}$
.
Further away from the particle crest, for
$(z-z_{b})/D>0.5$
, all of the terms except shear production
$P_{s}$
and dissipation
$\langle \overline{{\it\epsilon}}\rangle$
decay to zero, establishing equilibrium outer layer where the rate of production balances the rate of dissipation. The presence of local equilibrium is also evident from figure 4 as the flow in this region obeys the logarithmic law. Another measure to effectively quantify the presence of local equilibrium is the Townsend structure parameter,
$a_{1}=\langle \overline{u^{\prime }w^{\prime }}\rangle /2k$
(Townsend Reference Townsend1961). As shown in figure 9, it was found that
$a_{1}\sim 0.15$
, similar to the logarithmic layer for steady, unidirectional flows, in most of the fluid column over an entire oscillation cycle, confirming the presence of equilibrium turbulence and validating the approximate self-similarity in the log-law region for oscillatory flows over a rough bed.
For
$\mathit{Re}_{{\it\delta}}=400$
at peak phase: (a) instantaneous isosurfaces of the
${\it\lambda}_{2}$
parameter coloured by
$u^{\prime }/U_{\infty }$
, where the flow direction is the positive
$x$
-direction; (b) the anisotropy invariant map. Here,
${\it\xi}=(\text{III}_{b}/2)^{1/3}$
and
${\it\eta}=(-\text{II}_{b}/3)^{1/2}$
, where
$\text{II}_{b}$
and
$\text{III}_{b}$
are, respectively, second and third principal invariants of the Reynolds stress anisotropy tensor. Symbols represent: ○,
$(z-z_{b})/D\lesssim 0.1$
; ▫,
$0.1<(z-z_{b})/D\lesssim 0.8$
; ▵,
$(z-z_{b})/D>0.8$
.
3.3 Near-bed turbulence structure
Instantaneous near-bed flow structures are plotted in figure 10(a) in the form of
${\it\lambda}_{2}$
isosurfaces (Jeong & Hussain Reference Jeong and Hussain2006) for the peak phase. There exists a complex forest of highly dense near-bed structures, slightly inclined with respect to the rough bed. Contrary to the presence of near-wall horseshoe vortex structures in canonical smooth-wall flows, turbulent structures in this case appear to be ‘broken’. This could be due to the fact that, presence of roughness elements typically energizes wall-normal fluctuations as discussed earlier, that in turn distort these near-bed streamwise structures and cause flow isotropization near the bed; a phenomenon evident from these broken structures. The presence of these structures was speculated by Keiller & Sleath (Reference Keiller and Sleath1976) and Sleath (Reference Sleath1987); however, were not visualized based on the experimental data. To further quantify this tendency of flow isotropization and study the ‘shape’ of near-bed turbulence, local anisotropy in the Reynolds stresses by constructing the anisotropy invariant map (also known as Lumley’s triangle) (Choi & Lumley Reference Choi and Lumley2001) is quantified, which contains all of the possible turbulent states. As seen in figure 10(b), the data reveals following observations: close to the effective bed location, for
$(z-z_{b})/D\lesssim 0.1$
, a two-component type of turbulence is present as wall-normal fluctuations are much weaker compared with streamwise and spanwise fluctuations. Further above this region,
$0.1<(z-z_{b})/D\lesssim 0.8$
, wall-normal fluctuations gain in strength as discussed earlier, decreasing near-wall anisotropy and the turbulent state crosses over from upper to the lower branch of the triangle (note that the data in this region is almost parallel to
${\it\eta}$
axis indicating reduction in the anisotropy). Further away from the particle crest
$(z-z_{b})/D>0.8$
, the roughness effects are not pronounced and the cigar-shaped structures dominate, as streamwise fluctuations stay significantly larger than spanwise and wall-normal fluctuations.
For
$\mathit{Re}_{{\it\delta}}=400$
(a) phase variation of the quadrant contribution at a location close to the roughness crest; (b) spatial variation of the quadrant contribution plotted at peak phase. Symbols represent: ●, outward; ▸, ejection; ♦, inward; ▪, sweep. Here,
$Q$
indicates quadrant number. (c) Instantaneous isosurfaces of the
${\it\lambda}_{2}$
parameter at the peak phase coloured by ejection (red) and sweep (blue) motions.
Quadrant analysis (Wallace, Eckelmann & Brodkey Reference Wallace, Eckelmann and Brodkey1972) is also performed to understand the near-bed turbulence structure. Figure 11(a) shows phase variation of the quadrant contributions towards the Reynolds stress from ejection, sweep, inward and outward motions collected at a spatial location very close to the roughness crest. In general, significant contributions come from sweep and ejection motions, whereas contributions from inward and outward motions are minimal. The wall-normal distributions of the quadrant contributions at peak phase are plotted in figure 11(b), which shows that the sweep motions are more significant below the roughness crest and ejection motions play important role above the crest level. On the other hand, overall contributions from inward and outward motions throughout the fluid column remain insignificant. Figure 11(c) shows instantaneous isosurfaces of the
${\it\lambda}_{2}$
parameter (Jeong & Hussain Reference Jeong and Hussain2006) coloured by the indicator function for ejections (red) and sweep (blue) plotted for the peak phase, confirming that the upper portion of the near-bed structures relate much closely to the ejection, whereas the lower portion of these structures are due to sweep motions. Based on these results, question arises as to what is the distribution of the near-bed flow variables and if an assumption of Gaussian distribution is valid. The probability-density function (PDF) of the near-bed (a) velocity and (b) pressure fluctuations in the region above the rough bed at a peak phase for
$\mathit{Re}_{{\it\delta}}=400$
, are plotted in figure 12. As hypothesized, these distributions are highly non-Gaussian, and the PDFs fit well with the fourth-order Gram–Charlier model, given by
$$\begin{eqnarray}f_{GC4}({\it\psi}^{\prime })=\frac{\exp (-{\it\psi}^{\prime 2}/2)}{\sqrt{2{\rm\pi}}}\left[1+\frac{S_{{\it\psi}^{\prime }}}{3!}({\it\psi}^{\prime 3}-3{\it\psi}^{\prime })+\frac{F_{{\it\psi}^{\prime }}-3}{4!}({\it\psi}^{\prime 4}-6{\it\psi}^{\prime 2}+3)\right],\end{eqnarray}$$
where
${\it\psi}^{\prime }$
is the normalized fluctuating flow variable,
$S_{{\it\psi}^{\prime }}$
and
$F_{{\it\psi}^{\prime }}$
are skewness and flatness of its distribution, respectively. The data are recorded at around
$0.5D$
above the effective bed location
$z_{b}$
. The choice of this location is based on the high values of correlation between flow parameters and destabilizing hydrodynamic forces on the rough bed as reported by Ghodke, Apte & Urzay (Reference Ghodke, Apte and Urzay2014a
), a subject presently under investigation. Similar distributions are observed for all other Reynolds numbers under consideration. The peakedness as well as the long positive tails of such distributions can play critical role in destabilizing the particle bed and therefore should be included in probabilistic models for onset of erosion. Equation (3.5) facilitates incorporation of higher-order flow statistics describing the effect of near-bed bursting phenomena into such models.
For
$\mathit{Re}_{{\it\delta}}=400$
at a peak phase, PDF normalized by standard deviations: (a) streamwise velocity and (b) pressure fluctuations. DNS data, recorded at around
$0.5D$
above the effective bed location
$z_{b}$
, is represented by histograms and solid line represents fourth-order Gram–Charlier distribution (3.5). Here
$\text{Skewness}=-0.93$
,
$\text{flatness}=3.5$
for velocity, whereas
$\text{skewness}=-0.54$
and
$\text{flatness}=4.5$
for pressure fluctuations.
3.4 Summarizing near-bed oscillatory turbulence
A schematic view summarizing important near-bed flow features of oscillatory turbulent flow over a rough bed is shown in figure 13.
Schematic representation (not to scale) to show important near-bed processes.
Although the statistics of turbulent oscillatory flow differ significantly from that of unidirectional flow, the nature of near-bed processes in oscillatory flow is almost the same as that in unidirectional flow; especially with regards to significant contribution from ejection and sweep motions and reduction in large-scale anisotropy. As indicated in figure 13, quadrant analysis revealed the dominance of sweep motions below the roughness crest and ejection motions above the crest level. Ejections are therefore responsible for flux of TKE away from the wall, whereas sweeps will cause energy diffusion towards the wall. These sweep–burst cycles, therefore, may directly influence the sediment particle entrainment and deposition in coastal flows. Similar observations were also reported by Krogstad et al. (Reference Krogstad, Antonia and Browne1992) and Yuan & Piomelli (Reference Yuan and Piomelli2014b ) in their study of unidirectional flow.
Also, as in the case of steady flows, the presence of roughness is seen to greatly influence these near-bed turbulent structures by redistribution of the energy from streamwise fluctuations to spanwise and wall-normal fluctuations, causing near-bed flow isotropization and therefore resulting in an overall reduction of large-scale anisotropy. Again, these observations are consistent with the unidirectional steady flow literature (Krogstad & Antonia Reference Krogstad and Antonia1994; Antonia & Krogstad Reference Antonia and Krogstad2001; Ikeda & Durbin Reference Ikeda and Durbin2007; Yuan & Piomelli Reference Yuan and Piomelli2014b ).
Despite the similar nature of these near-bed flow motions, there exist some differences in the elementary processes that maintain the turbulence production in oscillatory flows. Contrary to the unidirectional flows, bed-induced production terms due to spatial inhomogeneity are comparable with the shear production term. Present data for the gravel-type roughness also indicates energy transfer from WKE to TKE and also from TKE to WKE; latter typically found absent in unidirectional flows (Raupach & Thom Reference Raupach and Thom1981; Raupach et al. Reference Raupach, Antonia and Rajagopalan1991; Mignot et al. Reference Mignot, Bartheleemy and Hurther2009).
Furthermore, unlike unidirectional flows, all of these important near-bed processes are pronounced for a range of phases in a flow cycle, more so, close to the peak phase for high Reynolds number oscillatory flows. In other words, depending on the friction velocity, turbulence is present only in a part of the oscillation cycle followed by flow re-laminarization. Such time varying nature of near-bed turbulence in oscillatory flows greatly affects destabilizing forces on particle-bed, a subject of future investigation.
4 Conclusions
Particle-resolved DNS of transitional and turbulent oscillatory flows over a rough bed were performed. A double-averaging technique is employed to study spatial heterogeneity at the roughness length scale.
Characterization of oscillatory flow in terms of mean velocity distribution, Reynolds stress variation, TKE budget, near-bed flow structures and PDF distributions of velocity and pressure perturbations is reported.
-
(i) Double-averaged streamwise velocity followed logarithmic law from early-acceleration until mid-deceleration phases (
${\it\omega}t=2{\rm\pi}/10$
to
$7{\rm\pi}/10$
), indicating the presence of fully-developed equilibrium turbulence; also confirmed by the uniform value of structure parameter (
$a_{1}\approx 0.15$
) and the ratio of shear production to dissipation rate close to unity in this region. -
(ii) The spatial inhomogeneity at the roughness length scale introduced additional production and transport terms in double-averaged TKE budget. Contrary to the unidirectional, steady flows over rough beds, bed-induced production terms, both
$P_{m}$
and
$P_{w}$
, for the present gravel-type roughness, were comparable to the shear production terms,
$P_{s}$
. Negative peak in
$P_{w}$
at the roughness crest level was observed and could be attributed to the conversion of turbulent kinetic energy to wake kinetic energy as a result of work of large-scale structures associated with
$\langle \overline{u^{\prime 2}}\rangle$
at this location against pressure drag of roughness elements. Secondary positive peak in
$P_{w}$
is observed close to the effective bed location, indicating the conversion of wake kinetic energy to turbulent kinetic energy of scales smaller than roughness length scale at this location. Bed-induced transport term was found to be negligible. -
(iii) Increasing the effective roughness energized wall-normal fluctuations while dampened the streamwise ones. This led to break-up of streamwise near-bed structures and caused reduction in large-scale anisotropy. As a result, broadly two separate regions, a near-bed region of two-component turbulence and an outer region of cigar-shaped turbulence were observed.
-
(iv) Quadrant analysis showed dominance of ejection and sweep motions over inward/outward interactions. The wall-normal distributions of the quadrant contributions showed that the sweep motions are more significant below the roughness crest and ejection motions play important role above the crest level.
-
(v) The PDF distributions of the velocity and pressure fluctuations showed a non-Gaussian behaviour that followed a fourth-order Gram–Charlier distribution, indicating importance of higher-order turbulence statistics. This finding is especially critical for the development of probabilistic models of sediment erosion.
Acknowledgements
This work was supported by the National Science Foundation (NSF-CBET-1133363). The computations were performed on the Texas Advanced Computing Center’s Stampede machine. Authors thank Dr J. Urzay of Stanford University for fruitful discussions during the CTR Summer Program 2014. Authors also thank anonymous reviewers for useful suggestions leading to the improvement of the manuscript.
Appendix
Table 2 presents the parameters used to perform grid resolution study for
$Re_{{\it\delta}}=400$
. Uniform grids are used in the region surrounding the roughness bed (up to
$z=10{\it\delta}$
) and grids are stretched in the wall-normal direction (above
$z=10{\it\delta}$
) using a hyperbolic tangent function. Case C1 represents current simulation parameters, while case C2 is finer grid resolution. To limit the computational requirements in case C2, grid points are clustered up to
$z=8{\it\delta}$
region with around 375 points per particle diameter. Furthermore, appropriateness of vertical extent of the domain was also tested by running case C3 with the increased domain height of
$45{\it\delta}$
(grid resolution in the near-bed region is kept same as that of case C1).
For
$\mathit{Re}_{{\it\delta}}=400$
, phase evolution of normalized double-averaged streamwise velocity profile at: (a)
$3{\rm\pi}/10$
, (b)
$5{\rm\pi}/10$
, (c)
$7{\rm\pi}/10$
. Symbols represent: ——, C1 (present grid resolution and domain); ●, C2 (refined grid resolution); ♦, C3 (present grid resolution and larger domain).
Grid resolution study for
$\mathit{Re}_{{\it\delta}}=400$
.
As seen in figure 14, the velocity profile is unchanged even after increasing the number of grid points (case C2), therefore confirming the adequacy of the present grid resolution (case C1). It can also be seen that, increasing the vertical extent of the domain (case C3) has no impact on the near-bed velocity distribution. In addition, the normalized boundary layer thickness, maximum friction factor and also the maximum values of period-averaged turbulence intensities showed no discernible impact of increasing the number of grid points.
