1. Introduction
The capability to predict the onset of particle motion is fundamental to numerous natural and industrial processes, ranging from sediment transport in rivers, estuaries and coastal waters to slurry flows in dredging and mining operations (Van Rijn Reference Van Rijn1993; Wierschem et al. Reference Wierschem, Groh, Rehberg, Aksel and Krülle2008). Identifying the threshold conditions for the transport of pollutants, such as microplastics (1 μm−5 mm), in environmental settings also requires a deep understanding of particle dynamics, which, for slightly buoyant particles, are particularly sensitive to hydrodynamic forces (Kane et al. Reference Kane, Clare, Miramontes, Wogelius, Rothwell, Garreau and Pohl2020). Therefore, a detailed examination of the early stages of particle motion is necessary (Nielsen Reference Nielsen1992; Buffington & Montgomery Reference Buffington and Montgomery1997; Vowinckel et al. Reference Vowinckel, Jain, Kempe and Fröhlich2016; Eyal et al. Reference Eyal, Enzel, Meiburg, Vowinckel and Lensky2021). The motion threshold has been extensively studied, driven by applications in civil and coastal engineering (Wilcock Reference Wilcock1993; Agudo & Wierschem Reference Agudo and Wierschem2012; Petit et al. Reference Petit, Houbrechts, Peeters, Hallot, Van Campenhout and Denis2015).
The criterion introduced by Shields (Reference Shields1936) is often used to determine the conditions of incipient sediment transport. According to Shields’ criterion, the motion of a sediment layer exposed to steady flow begins once the hydrodynamic force acting tangentially to the bottom surface exceeds the frictional resistance between the mobile layer and the underlying substrate, which is proportional to the layer’s submerged weight (see figure 1 a). This balance leads to the definition of the Shields parameter, given by the ratio between tangential and vertical forces. At the threshold of incipient motion, this parameter equals the sediment friction coefficient. The critical value of the Shields parameter depends solely on the granular Reynolds number, which quantifies the relative contributions of inertial and viscous forces at the sediment grain scale. By definition, the Shields parameter is inherently a statistical measure due to the randomness of the bottom geometry, as well as flow–particle and particle–particle interactions. In practice, the threshold is typically determined empirically for an ensemble of grains on a rough bed, with critical values accompanied by subjective definitions of different stages of motion (Breusers & Schukking Reference Breusers and Schukking1971).
(a) Side view of a hypothetical uniform layer of mobile spheres lying on a regular substrate, i.e. the Shields approach. (b) Top view and (c) side view of the present particle arrangement, where grey spherical particles form a fixed substrate, and the red sphere represents the mobile particle. The blue curve illustrates the Stokes velocity profile. Red arrows indicate the hydrodynamic drag, lift and effective weight forces (added mass and imposed pressure gradient forces are not shown in this sketch). The corresponding lever arms are shown as black double-headed arrows, referenced to the downstream-side contact points (which coincide in this side view). The symbols are categorised and described in table 1.
Recently, there has been a growing emphasis on developing more precise criteria for the motion threshold of a single particle and characterising particle dynamics immediately following the motion onset. These efforts have primarily involved combinations of detailed experiments (Charru et al. Reference Charru, Larrieu, Dupont and Zenit2007; Agudo & Wierschem Reference Agudo and Wierschem2012) and analytical methods (Agudo et al. Reference Agudo, Illigmann, Luzi, Laukart, Delgado and Wierschem2017; Topic et al. Reference Topic, Retzepoglu, Wensing, Illigmann, Luzi, Agudo and Wierschem2019). Unlike the statistical nature of the Shields criterion, these approaches provide well-defined thresholds for the motion of individual spherical particles, directly linking them to specific flow conditions, material properties, local bed structure and the influence of other surrounding particles (Agudo et al. Reference Agudo, Dasilva and Wierschem2014, Reference Agudo, Illigmann, Luzi, Laukart, Delgado and Wierschem2017).
It should be stressed that the previously mentioned studies focus exclusively on unidirectional steady flows. However, the insights gained from these studies have only limited applicability in settings with transient ambient flows. Oscillatory flows, in particular, are critical in many engineering and environmental settings, such as in the periodic agitation of submerged particles in multiphase mixtures or the back-and-forth motion near sediment bedforms induced by coastal waves (Sleath Reference Sleath1984; Nielsen Reference Nielsen1992). Oscillatory flows are well known to be more complex than unidirectional flows, even under laminar conditions. Nonlinear residual flows, such as steady streaming flows, can significantly affect particle dynamics on small and (very) large time scales compared with the period of oscillation (van Overveld et al. Reference van Overveld, Shajahan, Breugem, Clercx and Duran-Matute2022b ). The complex particle–fluid interactions can lead to self-organisation into a wide range of patterns, including chains and bands in diluted systems (van Overveld, Clercx & Duran-Matute Reference van Overveld, Clercx and Duran-Matute2023) or vortex ripples in dense systems (Sleath Reference Sleath1984). Despite their importance, oscillatory flows have received far less attention than unidirectional flows. To our knowledge, a similarly rigorous analytical approach for oscillatory flows has yet to be developed to determine the conditions at the onset of particle motion.
Overall, the unique description of the physics in an oscillatory flow needs one additional dimensionless parameter compared with the unidirectional case, since the system has an additional degree of freedom. This is analogous to general particle motion in oscillatory flows (Klotsa et al. Reference Klotsa, Swift, Bowley and King2007; van Overveld et al. Reference van Overveld, Breugem, Clercx and Duran-Matute2022a ,Reference van Overveld, Shajahan, Breugem, Clercx and Duran-Matute b ). Consequently, reported values for the Shields criterion (or its single-particle equivalent described by Agudo et al. Reference Agudo, Illigmann, Luzi, Laukart, Delgado and Wierschem2017) will depend not only on the granular Reynolds number, but also on an additional dimensionless parameter related to the oscillatory flow component, such as the frequency-dependent viscous length scale relative to the exposed particle’s diameter (Klotsa et al. Reference Klotsa, Swift, Bowley and King2007; Mazzuoli et al. Reference Mazzuoli, Kidanemariam, Blondeaux, Vittori and Uhlmann2016; van Overveld et al. Reference van Overveld, Shajahan, Breugem, Clercx and Duran-Matute2022b , Reference van Overveld, Clercx and Duran-Matute2023).
Due to the larger parameter space, the dynamics of particle motion in oscillatory flows may differ significantly from those in unidirectional flows. In the latter, once a particle starts moving, it is slightly lifted from its pocket in the bed, exposing a larger surface area to the ambient flow and increasing the drag force induced by the ambient flow, accelerating the particle even further. This feedback mechanism produces a sharp transition between stationary and dynamic states, assuming laminar flow conditions without fluctuations. In contrast, oscillatory flows introduce several additional complexities, including unsteady forces such as added mass and Basset history forces, which induce complex dynamics (van Overveld et al. Reference van Overveld, Clercx and Duran-Matute2023). Notably, analytical solutions to the problem have been derived for specific cases, such as the periodic motion of a small particle in an unbounded viscous fluid (Coimbra & Rangel Reference Coimbra and Rangel2001), which was later experimentally validated (Coimbra et al. Reference Coimbra, L’esperance, Lambert, Trolinger and Rangel2004). More generally, analytical approaches for particle motion in time-dependent flows have been developed, provided the velocity field is sufficiently smooth (Van Hinsberg, ten Thije Boonkkamp & Clercx Reference Van Hinsberg, ten Thije Boonkkamp and Clercx2011). Moreover, the velocity profile in an oscillatory flow is often non-monotonic, like for a Stokes boundary layer over a flat bottom. As a result, both the drag force and associated lever arm are transient quantities, yielding a complex time-varying torque. The particle’s dynamics can become even more intricate after the onset of motion. Depending on the flow’s oscillation period relative to the particle’s settling time, the particle may not settle in or even reach the next pocket in the bed before the flow reverses.
From a conceptual perspective, we classify particle behaviour into several progressive modes of motion. Initially, the particle may remain static when the hydrodynamic forces are too weak to induce any movement. As the forcing increases, the particle may exhibit a wiggling motion, characterised by brief periodic movements, typically near moments of maximum drag or shear stress. By definition, the amplitude of this motion is small, as the particle remains within a single pocket of the bed, falling back into its original position and coming to rest before the ambient flow reverses direction. This motion is expected to be periodic under laminar flow conditions where fluctuations are absent. As the driving forces continue to increase, the particle may start to roll between pockets in the bed, always maintaining contact with at least one other particle. Rolling is expected to precede any slipping or sliding, consistent with predictions for unidirectional flow (Agudo et al. Reference Agudo, Illigmann, Luzi, Laukart, Delgado and Wierschem2017). Even under laminar flow conditions, this rolling motion can be aperiodic, with the particle approaching a neighbouring pocket but not coming to rest before the flow direction reverses. When the driving force is sufficiently large to lift the particle temporarily from the substrate, it may hop between pockets (Topic et al. Reference Topic, Agudo, Luzi, Czech and Wierschem2022). As the driving force increases further, the particle may eventually become suspended for extended periods, moving freely without touching the substrate.
In this work, we comprehensively explore the threshold of particle motion under oscillatory flow conditions and the subsequent dynamics post-initiation. We first analyse the governing equations, provide theoretical predictions for the relevant dimensionless parameters and make quantitative predictions for the conditions at the onset of motion. We then extend our study using laboratory experiments, exploring the parameter space to identify the threshold for different flow conditions using a combination of particle tracking and particle image velocimetry (PIV). We complement our experimental findings with direct numerical simulation (DNS) to give insight into the transient forces acting on the particle and to explore the three-dimensional flow field near the onset of motion. The flow field around the particles in the simulations is fully resolved, with the particles accounted for using an immersed boundary method. To faithfully reproduce the motion onset, the numerical framework is supplemented by a particle–particle contact model that accounts for the build-up and release of elastic energy around contact points when the particle is still at rest. Throughout our broad approach, we particularly pay attention to the implications of the additional degree of freedom, distinguishing it from the more commonly studied unidirectional flow cases.
The work is structured as follows. We give an overview of the system in § 2, followed by the experimental approach in § 3 and the numerical approach in § 4. The results with accompanying discussion and comparison to other criteria are given in § 5, and concluding remarks in § 6.
2. Formulation of the problem
We first present a systematic overview of the physical system consisting of a single mobile, spherical particle resting on top of one fixed layer of monosized spherical particles, hereafter referred to as substrate. The oscillatory flow of an incompressible Newtonian fluid is driven by an externally imposed harmonic pressure gradient over a horizontal flat wall, resulting in an oscillatory boundary layer (OBL). Figures 1(b) and 1(c) provide accompanying sketches of the particle and flow configuration. For a regular substrate, the relevant dimensional parameters are the particle diameter
$D$
, the amplitude of velocity oscillations far from the bottom
$U_0$
, the angular oscillation frequency
$\omega$
, the kinematic viscosity of the fluid
$\nu$
, the fluid density
$\rho _{\!f}$
, the particle density
$\rho _s$
, the gravitational acceleration
$g$
, the exposure height of the particle to the ambient flow
$h$
and the elevation
$l$
of the particle centre relative to its contact points with the bed. These quantities and other main parameters are listed in table 1.
Overview of the main parameters describing material properties, oscillation characteristics, geometric configuration (cf. figure 1), and relevant forces and associated lever arms.
2.1. Governing equations
2.1.1. Fluid motion
The flow is governed by the continuity equation
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0 \end{align}
and by the incompressible Navier–Stokes equation
\begin{align} \frac {\partial \boldsymbol{u}}{\partial t} + \left (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\right )\boldsymbol{u} = -\frac {1}{\rho _{\!f}}\boldsymbol{\nabla } p + \nu {\nabla} ^2\boldsymbol{u} + U_0\omega \cos (\omega t)\widehat {\boldsymbol{x}}, \end{align}
where
$p$
is the pressure,
$t$
is the time and
$\boldsymbol{u}=(u,v,w)$
is the fluid velocity with the components referred to the Cartesian coordinates
${\boldsymbol{x}}=(x,y,z)$
. The unit vector
$\widehat {\boldsymbol{x}}$
points in the
$x$
-direction, aligned with the oscillatory flow, and
$z$
denotes the wall-normal coordinate pointing upwards. The last term in (2.2) corresponds to the imposed oscillatory pressure gradient.
In the laminar regime, the velocity in the oscillatory boundary layer over a flat, smooth wall located at
$z=z_0$
, is given by the Stokes boundary layer solution
\begin{align} \boldsymbol{u} = U_0\left [\sin (\omega t)-\exp \left (-\frac {z-z_0}{\delta _S}\right )\sin {\left (\omega t - \frac {z-z_0}{\delta _S}\right )}\right ]\widehat {\boldsymbol{x}}, \end{align}
where
$\delta _S=\sqrt {2\nu /\omega }$
denotes the boundary layer thickness, i.e. the characteristic viscous length scale (Landau & Lifshitz Reference Landau and Lifshitz1987; Acheson Reference Acheson1990).
The presence of bottom roughness adds a remarkable complexity to the problem, even when the roughness elements are monosized spherical particles arranged regularly on a horizontal wall (Mazzuoli et al. Reference Mazzuoli, Kidanemariam, Blondeaux, Vittori and Uhlmann2016; Agudo et al. Reference Agudo, Illigmann, Luzi, Laukart, Delgado and Wierschem2017). Notably, roughness can amplify nonlinear effects, such as steady streaming (Lyne Reference Lyne1971; Riley Reference Riley2001), and thereby influence the stability of the oscillatory boundary layer, eventually causing a transition to turbulence (Mazzuoli & Vittori Reference Mazzuoli and Vittori2016; Kaptein et al. Reference Kaptein, Duran-Matute, Roman, Armenio and Clercx2020). Presently, we consider a compact square arrangement of the bottom substrate.
The influence of substrate roughness on the undisturbed flow is reflected in the velocity field overlapping with the top layer of the substrate. In other words, the flow does not vanish at the top of the substrate, but instead penetrates into it to some extent. When modelling the flow using a Stokes profile, this effective penetration is accounted for by adjusting
$z_0$
, the virtual wall location, as illustrated in figure 1(c). In practice,
$z_0$
is chosen such that far above the substrate, the Stokes velocity profile matches the velocity field obtained from the DNS or experiments. Note that the mean velocity at
$z=z_0$
in the DNS or experiments does not necessarily vanish. The influence of substrate roughness on the flow experienced by the mobile particle is quantified by the particle exposure height
\begin{align} h=z_b+D-z_0, \end{align}
where
$z_b$
is the bottom elevation of the mobile particle.
2.1.2. Particle motion
The translational velocity
$\boldsymbol{u}_s$
of a spherical particle is described by Newton’s second law
\begin{align} \rho _sV_s \frac {{\rm d}\boldsymbol{u}_s}{{\rm d}t} = \boldsymbol{F}_{\kern-1.5pt B} + \boldsymbol{F}_{\kern-1.5pt C} + \boldsymbol{F}_{\kern-1.5pt S}, \end{align}
where
$V_s=\pi D^3/6$
is the particle volume,
$\boldsymbol{F}_{\kern-1.5pt B}$
is the resultant of body forces,
$\boldsymbol{F}_{\kern-1.5pt C}$
is the resultant of inter-particle contacts due to interactions with the substrate and
$\boldsymbol{F}_{\kern-1.5pt S}$
is the resultant of surface forces
\begin{align} \boldsymbol{F}_{\kern-1.5pt S}=\displaystyle \int _{\mathcal S}\left (-(p+P)\boldsymbol{n} + \boldsymbol{\tau }_\nu \right ) {\rm d}S, \end{align}
where
$P$
denotes the imposed pressure,
$\boldsymbol{\tau }_\nu$
the viscous stress tangential to the sphere surface
$\mathcal S$
, and
$\boldsymbol{n}$
the surface-normal unit vector. In the
$x$
-direction, following Maxey–Riley–Gatignol’s approach, the hydrodynamic force
$\boldsymbol{F}_{\kern-1.5pt S}$
can be modelled as the sum of the drag force and contributions from added mass, Basset force and imposed pressure gradient (Gatignol Reference Gatignol1983; Maxey & Riley Reference Maxey and Riley1983). While not strictly exact, this decomposition provides a reasonable approximation for small particles, where drag is predominantly viscous and the effects of flow unsteadiness, captured by the added mass, Basset and pressure gradient terms, can be effectively treated as additional contributions to the steady-state drag model. In the
$z$
-direction, the hydrodynamic force corresponds to the lift force.
The net body force
$\boldsymbol{F}_{\kern-1.5pt B}$
, proportional to the volume occupied by the particle, results from the sum of contributions due to gravity and buoyancy
\begin{align} \boldsymbol{F}_{\kern-1.5pt B} = - (\rho _s-\rho _{\!f})V_sg\hat {\boldsymbol{z}}, \end{align}
where
$g$
is the gravitational acceleration pointing in the negative
$z$
-direction,
$-\hat {\boldsymbol{z}}$
.
The rotation of a spherical particle is described by Euler’s second law
\begin{align} \frac {\rho _s\pi D^5}{60}\frac {{\rm d}\boldsymbol{\varOmega }_s}{{\rm d}t} = \boldsymbol{T}_{\kern-1.5pt C} + \boldsymbol{T}_{\kern-1.5pt S}, \end{align}
where
$\boldsymbol{\varOmega }_s$
denotes the particle angular velocity. The torques acting on the particle arise from interactions with the substrate
$\boldsymbol{T}_{\kern-1.5pt C}$
and hydrodynamic viscous stresses
\begin{align} \boldsymbol{T}_{\kern-1.5pt S}=\displaystyle \int _{\mathcal S} \left (\boldsymbol{x}-\boldsymbol{x}_c\right )\times \boldsymbol{\tau }_\nu \; {\rm d}S, \end{align}
with
$\boldsymbol{x}_c$
denoting the particle centre coordinates. For a rolling particle in the compact square arrangement (cf. figure 1
b), rotation occurs around the contact point at elevation
$z_D$
. The distance between
$z_D$
and the particle centre elevation is given by
\begin{align} l = z_b+\frac {D}{2}-z_D, \end{align}
as shown in figure 1(c).
2.2. Dimensional considerations
Let us now consider the following dimensionless variables of the hydrodynamic problem, indicated with asterisks,
\begin{align} \boldsymbol{u^*}=\dfrac {\boldsymbol{u}}{U_0}, \quad t^*=\omega t, \quad \boldsymbol{x}^* = \dfrac {\boldsymbol{x}}{\delta _S} \quad \boldsymbol{\nabla }^* = \delta _S\boldsymbol{\nabla }, \quad p^*=\dfrac {p}{\rho _{\!f} \delta _S U_0\omega }, \end{align}
which were also considered by Mazzuoli & Vittori (Reference Mazzuoli and Vittori2016). Moreover, the following dimensionless quantities characterising the particle dynamics are introduced
\begin{align} \boldsymbol{x}_s^* &= \dfrac {\boldsymbol{x}_s}{D}, \quad \boldsymbol{u}_s^* = \frac {\boldsymbol{u}_s}{U_0}, \quad \boldsymbol{\varOmega }_s^* = \frac {D\boldsymbol{\varOmega }_s}{U_0}, \quad \boldsymbol{F}_{\kern-1.5pt C}^* = \frac {\boldsymbol{F}_{\kern-1.5pt C}}{{\mathcal W}_s}, \quad \boldsymbol{T}_{\kern-1.5pt C}^* = \frac {2\boldsymbol{T}_{\kern-1.5pt C}}{{\mathcal W}_sD}, \nonumber \\ \quad \boldsymbol{F}_{\kern-1.5pt S}^* &= \frac {\boldsymbol{F}_{\kern-1.5pt S}}{\tau _{\nu 0}\pi D^2}, \quad \boldsymbol{T}_{\kern-1.5pt S}^* = \frac {2\boldsymbol{T}_{\kern-1.5pt S}}{\tau _{\nu 0}\pi D^3}, \quad \boldsymbol{F}_{\kern-1.5pt B}^* = \frac {\boldsymbol{F}_{\kern-1.5pt B}}{{\mathcal W}_s}, \quad \boldsymbol{T}_B^* = \frac {2\boldsymbol{T}_B}{{\mathcal W}_sD}, \end{align}
where the contact forces are normalised by the particle submerged weight
${\mathcal W}_s=(\rho _s-\rho _{\!f})gV_s$
, and the hydrodynamic forces are normalised by the characteristic viscous shear stress
$\tau _{\nu 0}=\rho _{\!f}U_0\omega \delta _S/2$
multiplied by the area of the particle surface
$\pi D^2$
.
We now use (2.11) and (2.12) to non-dimensionalise the governing equations for both the fluid phase and the particle motion. Hence, the dimensionless form of (2.1) and (2.2) reads
\begin{align} \boldsymbol{\nabla }^*\boldsymbol{\cdot }\boldsymbol{u}^*= 0 \end{align}
and
\begin{align} \frac {\partial \boldsymbol{u}^*}{\partial t^*} + \frac {1}{2}{\textit{Re}}_\delta (\boldsymbol{u}^*\boldsymbol{\cdot }{\nabla }^* )\boldsymbol{u}^* = -{\nabla }^*p^* + \frac {1}{2}{\nabla} ^{*2}\boldsymbol{u}^* + \cos (t^*)\hat {\boldsymbol{x}}, \end{align}
respectively, where
\begin{align} {\textit{Re}}_\delta = \frac {U_0 \delta _S}{\nu } \end{align}
is the Reynolds number of the oscillatory boundary layer. Substitution into (2.5) and (2.8) yields
\begin{align} s\frac {{\rm d}\boldsymbol{u}_s^*}{{\rm d}t^*} = \frac {1}{\varGamma }\left (\boldsymbol{F}_{\kern-1.5pt C}^* - \hat {\boldsymbol{z}}\right ) + 3\delta \boldsymbol{F}_{\kern-1.5pt S}^* \end{align}
for the particle translation and
\begin{align} \frac {s}{5}\frac {{\rm d}\boldsymbol{\varOmega }^*_s}{{\rm d}t^*} = \frac {1}{\varGamma }\boldsymbol{T}_{\kern-1.5pt C}^*+3\delta \boldsymbol{T}_{\kern-1.5pt S}^* \end{align}
for the particle rotation, where
\begin{align} s = \frac {\rho _s}{\rho _{\!f}} \end{align}
is the particle–fluid density ratio,
\begin{align} \delta = \frac {\delta _S}{D} \end{align}
is the normalised viscous length scale, and
\begin{align} \varGamma = \frac {\rho _{\!f}U_0\omega }{(\rho _s-\rho _{\!f})g} \end{align}
is the ratio between the amplitude of the oscillatory acceleration far from the bottom and the specific gravitational acceleration. From this point onward, all variables are non-dimensional unless otherwise indicated and the asterisks are omitted for notational convenience.
According to (2.16), when lift forces (i.e. the vertical component of
$\boldsymbol{F}_{\kern-1.5pt S}$
) are relatively small, vertical motion is primarily determined by the balance between contact forces and effective weight, with the dimensionless acceleration being inversely proportional to
$s\varGamma$
. Horizontal motion is governed by the horizontal components of the contact force
$\boldsymbol{F}_{\kern-1.5pt C}$
and the hydrodynamic force
$\boldsymbol{F}_{\kern-1.5pt S}$
, the latter comprising contributions from the imposed pressure gradient and drag, which are proportional to
$\delta /s$
.
The dynamic problem involves nine distinct dimensional parameters (
$U_0$
,
$D$
,
$\omega$
,
$\rho _{\!f}$
,
$\rho _s$
,
$\nu$
,
$g$
,
$h$
,
$l$
) and is, therefore, uniquely characterised by six independent dimensionless parameters. Four of these already appeared in the governing equations:
$\delta$
,
${\textit{Re}}_\delta$
,
$s$
and
$\varGamma$
. In addition, the dynamics of the exposed particle depend on two geometric parameters,
$H=h/D$
and
$L = l/D$
, which relate to the substrate arrangement and are defined by (2.4) and (2.10), respectively.
Combinations of these numbers result in other dimensionless quantities that are frequently considered in sediment transport problems, such as the particle mobility number
\begin{align} \psi = \frac {\rho _{\!f}U_0^2}{(\rho _s-\rho _{\!f})gD} = \frac {1}{2}\delta {\textit{Re}}_\delta \varGamma , \end{align}
which is the ratio between the convective force and the submerged weight of the particle (Mazzuoli et al. Reference Mazzuoli, Kidanemariam, Blondeaux, Vittori and Uhlmann2016, Reference Mazzuoli, Kidanemariam and Uhlmann2019; Vittori et al. Reference Vittori, Blondeaux, Mazzuoli, Simeonov and Calantoni2020).
Finally, the bulk flow velocity
$U_0$
does not necessarily represent the local flow around the mobile particle, especially when
$\delta$
is large. The flow conditions near the substrate are better captured by the Shields parameter, representing the dimensionless shear stress at the bed as
\begin{align} \theta = \frac {\tau _0}{(\rho _s-\rho _{\!f})gD}, \end{align}
where
$\tau _0$
denotes the instantaneous bottom shear stress. According to the Shields criterion, this stress represents the horizontal hydrodynamic force per unit area acting on the top layer of particles resting on the substrate. This force is approximately equal to the shear stress evaluated at an elevation
$z_0$
, which typically lies within one particle diameter
$D$
above the base of the top layer
$z_b$
(see figure 1
a). In our configuration, where the mobile layer consists of a single particle and the particle-induced disturbance to the flow is negligible,
$z_0$
is close to
$z_b$
(see figure 1
c). Consequently, the shear stress at
$z_0$
closely approximates the actual shear stress acting at the base of the mobile particle. In the absence of a substrate (i.e. a single sphere lying on a smooth wall),
$\tau _{0}$
can be computed by evaluating the derivative of the Stokes profile (2.3) at
$z=z_0=0$
, yielding
$\theta = (\sqrt {2}/2)\delta \varGamma \sin (t+\pi /4 )$
. Finally, we stress that, in contrast to steady flow conditions, the Shields parameter
$\theta$
in (2.22) oscillates over time, implying that both the magnitude and the duration of the hydrodynamic forcing determine the onset of motion.
2.3. Model for the motion threshold
The detailed derivation underlying the model for the motion threshold is given in Appendix A, while the main concepts are summarised in this section.
Under the assumption that an exposed particle is more likely to roll rather than slip or slide out of a pocket in the bed (Agudo et al. Reference Agudo, Illigmann, Luzi, Laukart, Delgado and Wierschem2017), the motion threshold is governed by a torque balance. This balance includes a stabilising contribution from the submerged weight and typically destabilising contributions from drag, lift, added mass and the imposed pressure gradient. However, these torques may be out of phase, meaning some may have stabilising effects during parts of the oscillation cycle. The Basset history term is omitted from the present balance, since its contribution is expected to be small for the range of
${\textit{Re}}_\delta$
considered. In the absence of a substrate, its magnitude is approximately
$10\,\%{-}20\,\%$
of drag, while precise evaluation in the present geometry remains challenging. Moreover, since the history force is related to flow acceleration, its phase does not coincide with that of the drag force, further limiting its influence on the onset of particle motion.
Lubrication forces are not included, as prior studies have shown good agreement without them (Agudo et al. Reference Agudo, Illigmann, Luzi, Laukart, Delgado and Wierschem2017), and they are not expected to be significant compared with dominant contributions such as drag and pressure gradients in the present parameter range. The dimensionless spanwise component of the particle torque
$\boldsymbol{T}$
referred to an arbitrary point is defined by
\begin{align} T_{y} = T_{\textit{Sy}}^{\textit{drag}}+T_{\textit{Sy}}^{\textit{lift}}+T_{\textit{Sy}}^{\textit{AM}}+T^{\boldsymbol{\nabla}\kern-1pt p}_{\textit{Sy}}+ T_{\textit{By}} + T_{Cy}, \end{align}
where the terms on the right-hand side denote contributions of drag, lift, added mass, pressure gradient, submerged weight and inter-particle contacts, respectively. In this section, let us evaluate the torque contributions referring to the contact point elevation at the incipient rolling conditions, i.e. when
$T_{y}=0$
and the mobile particle has only two contact points aligned in the spanwise direction, such that also
$T_{Cy}=0$
. Therefore, in the definition of the spanwise component of the hydrodynamic torque
$T_{\textit{Sy}}$
provided by (2.9), the coordinates
$(x_c,z_c)$
are replaced by
$(x_D,z_D)$
.
The primary challenge lies in relating the hydrodynamic torque contributions to the ambient flow, for which we adopt an approach similar to that of Agudo et al. (Reference Agudo, Illigmann, Luzi, Laukart, Delgado and Wierschem2017) for steady shear flows. In general, there is no analytical formulation for the velocity field
$u$
that captures all substrate-induced variations. Therefore, we introduce a horizontal plane-average to obtain an effective velocity profile that represents the mean flow. This averaged profile is defined as
\begin{align} \langle u\rangle (z) \equiv \frac {1}{\mathcal A}\iint _{\mathcal A} u\, {\rm d}x \,{\rm d}y, \end{align}
where
$\mathcal A$
denotes the substrate area. The averaging procedure accounts for spatial velocity variations induced by the substrate, which are present even in non-turbulent flow conditions. Notably, while the ambient flow is free of turbulent fluctuations, random velocity fluctuations may still appear in the wake downstream of the mobile particle.
Although
$\langle u\rangle$
cannot be computed analytically, our numerical simulations provide full access to the three-dimensional velocity field, allowing us to evaluate the plane-average directly. The resulting mean velocity profiles are well approximated by a vertically shifted Stokes solution (2.3) (see figure 16 in Appendix A) for the parameter values presently considered. The torque due to the hydrodynamic drag is defined as the moment of the drag force with respect to the contact point between the particle and the substrate. For the spherical particle, this torque is modelled as
\begin{align} T_D\propto \pi D^2\,u_1, \end{align}
according to (A5), where
\begin{align} u_1 = \dfrac {1}{D}\int _{z_b}^{z_b+D} (z-z_D) \langle u\rangle \,{\rm d}z \end{align}
is the first moment of the streamwise ambient velocity (see (A3)) with
$z_b$
and
$z_D$
denoting the
$z$
-coordinates of the lowest point of the exposed particle and of the rotation axis, respectively (see figure 1
c). Similarly, the contribution of lift to the hydrodynamic torque is based on Saffman’s (Reference Saffman1965) expression, see (A15). Notably, both the expressions of drag and lift in dimensionless form contain correction coefficients that depend on the instantaneous particle Reynolds number
\begin{align} {\textit{Re}}_s(t) = \dfrac {u_cD}{\nu } = \frac {{\textit{Re}}_\delta }{\delta }\left [\sin (t)-e^{-\zeta }\sin \left (t-\zeta \right )\right ]\!, \end{align}
where
$u_c$
denotes the ambient velocity at the same elevation as the particle centre and
$\zeta =(H-1/2)/\delta$
is the normalised elevation of the particle centre above the zero level of the velocity profile
$z_0$
.
The ratio
$\varUpsilon$
between the aforementioned destabilising contributions and the stabilising contribution due to submerged weight is given by
\begin{align} \varUpsilon &\equiv \dfrac {T_{\textit{Sy}}^{\textit{drag}}+T^{\boldsymbol{\nabla}\kern-1pt p}_{\textit{Sy}}+T_{\textit{Sy}}^{\textit{AM}}+T_{\textit{Sy}}^{\textit{lift}}}{T_{\textit{By}}} \nonumber\\ &=4\varGamma \left [ 9C'\delta ^2H g\left (t,\frac {H}{\delta }\right )\,L_D + h(t,\zeta )\,L + \frac {9.66}{\pi \sqrt {2}}\,C^{\prime \prime }\dfrac {\delta ^3}{{\textit{Re}}_\delta }\sqrt {{\textit{Re}}_s^3}L_L \right ]\!, \end{align}
where
$C'$
and
$C''$
are correction coefficients for the drag and lift forces, respectively, and
$L_D$
and
$L_L$
are their corresponding lever arms, as detailed in Appendix A. The function
\begin{align} g\left (t,\frac {H}{\delta }\right ) = \sin (t)-\frac {\sqrt {2}}{2}\frac {\delta }{H}\left [\sin \left (t-\frac {\pi }{4}\right )- e^{-H/\delta }\sin \left (t - \frac {H}{\delta } - \frac {\pi }{4}\right )\right ] \end{align}
captures the effect of shielding by the substrate, and
\begin{align} h(t,\zeta ) = \frac {3}{2}\cos (t)-\frac {1}{2}e^{-\zeta }\cos \left (t-\zeta \right ) \end{align}
groups the contributions of the pressure gradient and added mass. Equation (2.28) allows us to define the criterion for the initiation of rolling motion, depending on (dimensionless) time
$t$
and five dimensionless numbers:
\begin{align} \varUpsilon \left (t;\,{\textit{Re}}_\delta ,\,\delta ,\,\varGamma ,\,H,\,L\right ) \geqslant 1. \end{align}
Note that the contribution of
$s$
to the threshold is contained within the parameter
$\varGamma$
. Once the threshold is exceeded,
$s$
becomes an independent control parameter that governs the relative importance of particle inertia in the dynamics.
Both
$H$
and
$L$
are primarily determined by the geometry of the bottom substrate. For a substrate composed of spheres in a compact square arrangement,
$L = \sqrt {2}/4 \approx 0.35$
. The value of
$H$
is related to the effective position of the virtual origin
$z_0$
of the velocity profile within the substrate crevices. In a similar substrate configuration but under steady flow conditions, Agudo et al. (Reference Agudo, Illigmann, Luzi, Laukart, Delgado and Wierschem2017) and Topic et al. (Reference Topic, Agudo, Luzi, Czech and Wierschem2022) positioned
$z_0$
at a distance
$0.077D$
below the top of the substrate, corresponding to
$H=2L+0.077\approx 0.78$
. In our DNS results,
$H\approx 0.84$
(cf. table 2).
To get familiar with the
$\varUpsilon$
-criterion (2.28)–(2.31) and distinguish the different contributions to the particle torque, figure 2(a) shows the time variation of
$\varUpsilon$
for constant parameters
$H=0.8$
,
$L=0.35$
and
$\varGamma =0.45$
, which are representative of the present experimental and numerical conditions.
(a) Temporal evolution of the motion threshold parameter
$\varUpsilon$
, as defined by (2.28), for typical parameter values in the present experiments and simulations (specifically matching DNS run 3, see table 2):
$\delta =0.12$
,
${\textit{Re}}_\delta =41$
,
$\varGamma =0.45$
,
$s=1.81$
,
$H=0.8$
and
$L=0.36$
(black line). Dotted lines represent relative contributions of hydrodynamic drag (blue), lift (green) and added mass plus imposed pressure gradient (red). (b) Normalised time interval
$\Delta t_{\varUpsilon \gt 1}$
during which the motion threshold condition (2.31) is exceeded, as a function of
$\varGamma$
(proportional to the forcing strength) for various values of
$\delta$
, with other parameters as in panel (a). Crosses mark the characteristic time scale
$\tau$
(2.33) for the particle to reach the crest of the substrate.
2.4. Particle motion just above the motion threshold
Assuming that the particle rolls during the early stages of its motion, as verified in § 5.1, we define the spanwise rotation angle
$\phi _s$
as the angle swept by the particle while rolling from its initial position to the crest of the substrate particle, before settling into the adjacent pocket. For the present substrate configuration, a compact square arrangement, we find
$\phi _s=\pi /2-\arctan (\sqrt {2})\approx 0.62$
. The angle
$\phi _s$
may alternatively be obtained by integrating Euler’s second law (2.17) twice with respect to time, which yields
\begin{align} \phi _s = \iint _0^\tau \frac {{\rm d}\varOmega _{sy}}{{\rm d}t} {\,{\rm d}t}' \,{\rm d}t \approx \frac {5\delta {\textit{Re}}_\delta }{4s}\tau ^2\left (3 \delta T_{\textit{Sy}}-\frac {1}{\varGamma }T_{Cy} \right )\!, \end{align}
where
$\tau$
denotes the characteristic (dimensionless) time for the particle to reach the substrate crest, hereafter referred to as the ‘time to crest’, and
$T_{\textit{Sy}}$
and
$T_{Cy}$
represent the characteristic values of the respective torque components during this time. Close to the motion threshold, the difference between stabilising and destabilising torques is small. Based on the scaling analysis in § 2.2, the resultant dimensionless torque (i.e. the left-hand side of (2.17)) is expected to be of order unity during the early stages of motion. Furthermore, the same analysis suggests that
$3\delta T_{\textit{Sy}}$
and
$T_{Cy}/\varGamma$
, i.e. the terms between brackets in (2.32), are typically of the same order of magnitude. Then, the characteristic time to crest scales as
\begin{align} \tau \sim \sqrt {\frac {s}{\delta {\textit{Re}}_\delta }} = \sqrt {\frac {s}{2K_C}}, \end{align}
where
$K_C=U_0/(\omega D)$
is the Keulegan–Carpenter number, describing the relative importance of hydrodynamic drag to inertial forces in oscillatory flows. The range of
$K_C$
values covered in this study spans from 2.5 to 80 in the simulations and from 0.58 to 10 in the experiments.
Figure 2(b) shows the normalised time interval
$\Delta t_{|\varUpsilon |\gt 1}$
, representing the fraction of the oscillation period during which the motion threshold (2.31) is met. The prediction of
$\Delta t_{|\varUpsilon |\gt 1}$
is crucial for determining whether a particle can accelerate sufficiently to overcome substrate barriers and settle into neighbouring pockets before flow reversal. At low forcing (small
$\varGamma$
), the threshold is typically not exceeded (
$\Delta t_{|\varUpsilon |\gt 1}=0$
). As
$\varGamma$
increases,
$\Delta t_{|\varUpsilon |\gt 1}$
grows rapidly and asymptotically approaches
$2\pi$
, corresponding to continuous particle motion throughout the oscillation cycle. The dependence on
$\delta$
reflects the effect of flow uniformity at the particle scale, where larger values of
$\delta$
(more shear-like flow) promote particle motion at lower
$\varGamma$
, for otherwise identical parameters. The figure also shows the predicted time to crest (2.33), marked by crosses. Below each cross, the duration during which
$|\varUpsilon |\geq 1$
is too short for the particle to roll over the crest of the substrate particle, resulting in wiggling motion. Above the cross, the particle can roll out of its pocket and reach a neighbouring one. Notably, for small
$\delta$
, the time to crest diverges. In this regime, the oscillation frequency is sufficiently high that flow reversal occurs before the particle can accelerate sufficiently, preventing it from rolling over the substrate.
3. Experimental approach
The experiments are carried out in an oscillatory flow tunnel (OFT) and serve a dual purpose. First, they enable the characterisation of flow conditions and particle behaviour as a function of the problem’s parameters close to the motion threshold. Second, they validate and inform the numerical simulations.
3.1. Experimental set-up
Figure 3 shows a schematic overview of the set-up. We use a custom-made OFT consisting of an acrylic tank with inner dimensions
$(100\times 50\times 30)\,\textrm {cm}$
in the streamwise, spanwise and vertical directions, respectively. A U-shaped lid is suspended within this tank, creating a central section with two parallel flat plates spaced
${30}\,\textrm {mm}$
apart. Each end of the tank contains a vertical section with a free surface, measuring
${50}\,\textrm {mm}$
in width. The tank is filled with tap water with mass density
$\rho _{\!f}=(0.999\pm 0.001)\times {10^{3}}{\,\textrm {kg}\,\textrm {m}^{- 3}}$
and kinematic viscosity
$\nu =(1.05\pm 0.01)\times {10}^{-6}{\,\textrm {m}^2\,\textrm {s}^{- 1}}$
. Special care has been taken to ensure a watertight seal between the lid and the side walls of the outer tank at the front and back (out of plane). A partially submerged piston with a width of
${48}\,\textrm {mm}$
is placed in the left column, covering approximately
$77\,\%$
of the free surface area (as seen from above). The piston is driven by a PID-controlled linear motor (LinMot P01-37×120F/100×180-HP) in a vertical oscillatory motion. This motion modulates the water level in the left column, generating a hydrostatic pressure difference between the left and right columns, which in turn drives the flow through the central horizontal section. The piston does not fully cover the free surface in the lateral direction, but it is sufficiently large to generate the necessary flow through vertical oscillations of the water level. A tight seal or full interfacial coverage is not required, as this would involve extreme forces and special care to prevent friction.
(a) Schematic of the oscillatory flow tunnel (OFT), where a harmonically oscillating piston (grey) drives flow between two parallel flat plates. On the bottom, a single mobile particle (white) lies on top of a fixed monolayer (black), with a vertical laser sheet (green) illuminating the cross-section parallel to the oscillation direction. (b) Experimental snapshot (
$\delta \approx 0.12$
and
$s=1.81$
) showing the substrate, exposed particle and laser-illuminated tracer particles. (c) A contrast-enhanced snapshot overlaid with velocity vectors obtained using PIV, scaled and colour-coded by magnitude (green to red for increasing velocity magnitude).
The substrate consists of a monolayer of spheres with a diameter of
$(5.950\pm 0.005)\,\textrm {mm}$
, arranged in a square lattice with
${6.0}\,\textrm {mm}$
centre-to-centre spacing, forming a
$25\times 26$
particle grid aligned with the oscillation direction. The spheres are securely positioned in circular holes (
${3.0}\,\textrm {mm}$
in diameter) on a bottom plate, with a 5 mm-high surrounding frame preventing the layer from shifting. The bottom plate and frame are CNC milled to ensure a highly regular spacing and alignment of the substrate spheres. A single spherical particle with diameter
$D={(5.950\pm 0.005)}\,\textrm {mm}$
is placed on top of the substrate approximately in the centre of the grid. Its density is either
$\rho _s={(1.09\pm 0.05) \times {10}^{3}}{\,\textrm {kg}\,\textrm {m}^{-3}}$
or
$\rho _s={(1.81\pm 0.05)\times {10}^{3}}{ \,\textrm {kg}\,\textrm {m}^{-3}}$
for light and heavy particles, respectively. All spheres used in the experiments are plastic BB pellets made of polylactic acid (PLA) and are specially designed to be smooth and highly spherical.
We use PIV to characterise the oscillating flow field, specifically to measure the mean flow velocity
$U_0$
far from the bed (e.g. as in figure 4) and to obtain the velocity profile near the exposed particle (cf. figure 16). A vertical two-dimensional slice of the velocity field is captured above the bed and around the exposed particle using a laser sheet approximately
${5}\,\textrm {mm}$
wide, similar to the particle’s diameter. The sheet is centred on the exposed particle and aligned with both the substrate grid and the oscillation direction, which limits out-of-plane reflections. The laser illuminates non-fluorescent polycrystalline tracer particles (Optimage PIV Seeding Powder with nominal diameter of
${100}\, {\unicode{x03BC} } \rm{m}$
) and operates in a double-pulsed mode, emitting two pulses separated by
${5.0}\,\textrm {ms}$
, with pulse pairs generated at a frequency of 15 Hz. A RedLake MegaPlus II camera equipped with a Nikon 28 mm
$f/2.8$
lens records the experiment from the side (along the
$y$
-axis), perpendicular to the laser sheet plane (
$xz$
) to further minimise reflections from the substrate and the mobile sphere. Nevertheless, some reflections near the particle are unavoidable due to the finite thickness of the sheet. Figure 3 shows a typical snapshot, where the vertical laser sheet illuminates the substrate at the bottom of the image, the mobile particle on top and the tracer particles dispersed in the fluid for PIV analysis.
(a) Temporal evolution of the space-averaged flow velocity above the substrate,
$u$
(blue), and the (normalised) particle position,
$x_s$
(red), over approximately three oscillation periods, obtained from the OFT experiments (
$s=1.81$
,
$\delta \approx 0.11$
,
${\textit{Re}}_\delta \approx 150$
,
$\varGamma \approx 0.13$
,
$\psi \approx 1.1$
). The velocity amplitude is gradually ramped up to identify the conditions for incipient motion. The complementary DNS data (black) show the trajectory of a rolling particle just above the motion threshold for run 4 (
$\psi =1.1$
). Roman numerals mark every eight frames, corresponding to the snapshots in panel (b). The inset shows the particle motion and velocity over a longer duration, with the highlighted region corresponding to the main plot.
The PIV recordings are analysed with the image processing software PIVview2C (v3.9.3). For each image pair, instantaneous velocity vectors are computed within interrogation windows of
$32 \times 32$
pixels (roughly corresponding to
$3\times {3}\,\textrm {mm}$
, approximately the particle radius), with 50 % overlap between adjacent interrogation windows (Prasad Reference Prasad2000). Each window typically contains 10–20 tracer particles. Figure 3(c) shows a typical frame after PIV analysis. A high-pass filter and pixel value threshold are applied to the original image to enhance the contrast of the tracer particles. The coloured arrows represent velocity vectors within each interrogation window.
Finally, a preferential path for the mobile particle may arise due to secondary flows or minor substrate irregularities. While the latter contributions are minimised by using CNC-milled components and uniform particles, factors such as imperfect levelling of the tank and inertial effects from the piston’s motion (e.g. vortex shedding) may still contribute to minor asymmetries in the flow. While these asymmetries are typically irrelevant, they can influence particle motion near the threshold, where small torque differences may result in a preferred direction of particle movement.
3.2. Measurement approach
For a given particle–fluid combination and substrate geometry, the density ratio
$s$
and the geometrical parameter
$L$
remain constant across all experiments. The parameter
$H$
, primarily determined by the substrate geometry (as discussed in § 2.3), is also assumed to remain approximately constant. Three independent degrees of freedom remain (
$\delta$
,
${\textit{Re}}_\delta$
,
$\varGamma$
), all depending on the flow conditions. For comparison with DNS results, we also consider the particle mobility number
$\psi =\delta {\textit{Re}}_\delta \varGamma /2$
(see (2.21)). In each experiment, the oscillation frequency
$\omega$
is fixed while gradually increasing the piston’s stroke, thereby increasing the fluid velocity amplitude
$U_0$
. The constant frequency implies that the value of
$\delta$
remains constant within each experiment, but varies across different trials as the frequency is adjusted. At each stroke value, the piston oscillates for ten periods before its peak-to-peak amplitude is slightly increased: by
${5}\,\textrm {mm}$
at lower amplitudes and by
${1}\,\textrm {mm}$
near the motion threshold. The oscillation period ranges between
$1.0$
s and
${1.8}\,\textrm {s}$
, corresponding to the Stokes boundary layer thickness
$\delta _S\approx 0.58{-}{0.78}\,\textrm {mm}$
, which is significantly smaller than the
${30}\,\textrm {mm}$
spacing between the parallel plates. This ensures that the velocity profile above the bottom remains unaffected by the OFT’s upper boundary (van Overveld et al. Reference van Overveld, Breugem, Clercx and Duran-Matute2022a
). The ranges of dimensionless parameters explored near the motion threshold are listed in table 2.
Values of the dimensionless parameters that characterise the DNS cases, alongside experimental ranges near the motion threshold (corresponding to the symbols in figure 6). The Reynolds number based on the particle diameter is denoted by
${\textit{Re}}_D$
. The final column reports the threshold value of the mobility number
$\psi _{{\textit{tr}},w}$
for the onset of wiggling motion. In the DNS, the density ratio is fixed at
$s=1.813$
and the value of
$\psi$
is explicitly prescribed. In the experiments,
$\psi$
is not an independent control parameter but follows from the other dimensionless parameters.
4. Numerical approach
Direct numerical simulations (DNSs) are conducted to extend the experimental results and to determine the dynamics of the mobile particle during its early stages of motion. The simulations are designed to closely replicate the experiments at
${\textit{Re}}_\delta \approx 164$
, using a matching geometry with the exception that the substrate in the DNS extends across the entire horizontal bottom plane. The incompressible Navier–Stokes equations (2.2) are solved numerically using a second-order central-difference scheme to discretise the spatial derivatives and a Runge–Kutta-based fractional-step method to advance in time. At the top boundary (far above the substrate), a stress-free condition is imposed, while periodic boundary conditions are applied in both horizontal directions. At the bottom of the numerical domain, i.e. on the plane supporting the fixed bed, a no-slip/no-penetration condition is enforced. The no-slip/no-penetration boundary conditions at the fluid–solid interfaces of both the fixed and mobile particles are enforced using an immersed boundary method (IBM) developed for particulate flows by Uhlmann (Reference Uhlmann2005). This method employs a direct forcing approach, adding a localised volume force term to the momentum equations. The additional forcing term is explicitly computed at each time step as a function of the particle position and velocity, providing the stress distribution at each particle’s surface. The dynamics of the mobile particle are determined by Newton’s second law for its translation and by Euler’s second law for its rigid-body rotation, according to (2.5) and (2.8), respectively. The hydrodynamic forces and torques result from the zeroth and first angular moments of the stresses, respectively, integrated over the particle surface. The contribution of inter-particle contacts is given by the sum of the contact forces between a considered particle and its neighbours, computed using a discrete element method (DEM) similar to that proposed by Kidanemariam & Uhlmann (Reference Kidanemariam and Uhlmann2014). The present approach, which is described in Appendix B, accounts for inter-particle static friction, simultaneously incorporating elastic, inelastic and frictional contributions to the tangential force at each contact point. Correct evaluation of the inter-particle contact force is necessary for describing the early stages of particle motion because the mobile particle may either roll or slide. Tables 2 and 3 show the values of the parameters characterising the flow conditions and the computational domain, respectively. The density ratio is fixed at
$s=1.813$
in all simulations, matching experiments with heavier particles. The lighter particles, with density ratios closer to unity, could not be accurately simulated due to stability limitations of the numerical approach.
Simulation domain parameters, where
$\varLambda$
represents the size of the computational domain, while
$n$
indicates the number of grid points used to discretise the domain along the corresponding axis.
$\Delta x$
and
$\Delta t$
denote the grid spacing and time step used in the simulations, respectively.
5. Results and discussion
5.1. Description of the early stages of particle motion
The flow velocity within the OFT is averaged over the PIV measurement area (e.g. as shown in figure 3), excluding regions near the top plate, the substrate and the exposed particle to minimise boundary effects. This yields a representative, spatially averaged velocity for the undisturbed flow field. Figure 4(a) shows the space-averaged flow velocity (blue curve), oscillating harmonically with an amplitude
$U_0\approx {(0.23\pm 0.01)}{\,\textrm {m}\,\textrm {s}^{- 1}}$
, where the error follows from PIV correlations and subsequent averaging procedure, for
$\delta =0.11$
and
$s=1.81$
. Figure 4(b) shows some frames of the particle motion.
The particle’s horizontal centre position is determined by filtering out tracers and identifying the median position of pixels with intensities above half the maximum value, restricted to the visible portion of the particle above the substrate. This method allows us to accurately and consistently resolve the particle’s centre position to within a few percent of its diameter. The red curve in figure 4(a) shows the particle position over time under conditions near the motion threshold. The DNS results (black curve) closely match the experimental results, exhibiting similar motion characteristics, including the phase of motion initiation, typical excursion lengths and movement directions. This agreement suggests that, in the absence of significant velocity fluctuations in the ambient flow, deterministic predictions of the incipient particle motion are possible. Quantitative differences arise due to the system’s sensitivity to small perturbations near the motion threshold, where minor positional variations alter the particle’s exposure to the flow, leading to significant long-term deviations in its trajectory, such as settling into a different pocket in the substrate.
For the parameter values considered in the experiments, the particle motion was observed to begin around the maximum velocity phases (frames I, IV, IX). In some cases, the particle rolled without escaping its initial pocket (frame II), a behaviour we refer to as ‘wiggling motion’. In others, sustained motion continued until the flow reversal (frame V), after which the particle rolled across multiple pockets in the opposite direction (frames VI, VII), and occasionally overshot (frames V, VII, X) before settling. The range of flow velocities (captured in the forcing parameters
${\textit{Re}}_\delta$
,
$\varGamma$
and
$\psi$
) for which wiggling motion is observed is relatively narrow, as even a slight increase of the hydrodynamic torque intensity or duration above the critical threshold can cause the particle to roll over the substrate. As the particle leaves its initial position and falls into one of the neighbouring pockets, the effects of the nonlinearities associated with the inter-particle contacts reflect on the hydrodynamic force, which becomes impossible to predict deterministically. Therefore, from now on, only the early stages of particle motion are considered.
While figure 4(b) provides a detailed view of the particle translation, it does not resolve its rotation due to the laser sheet illumination. Therefore, supplementary experiments are conducted with the laser turned off, as shown in figure 5. These snapshots support that the particle rotates throughout its motion, also when it moves between the pockets of the substrate, which is consistent with Agudo et al.’s (Reference Agudo, Illigmann, Luzi, Laukart, Delgado and Wierschem2017) experiments in steady flow conditions, where rolling always characterised the initiation of motion.
Consecutive snapshots from an experiment with
$s=1.81$
,
$\delta \approx 0.12$
,
${\textit{Re}}_\delta \approx 190$
and
$\varGamma \approx 0.13$
, showing that the particle (with black dots drawn on its surface) rolls when it moves over the substrate (dark grey). The brightness has been enhanced to improve substrate visibility, and the time interval between snapshots is 0.05 s.
To confirm the deterministic nature of the motion threshold and wiggling, figure 6 shows the streamwise centre position of a wiggling particle over ten oscillation cycles, projected onto two cycles for clarity. The particle’s centre position is resolved to within a few percent of its diameter, which is sufficient to reliably detect the wiggling motion and demonstrate its repeatability. When the particle comes to rest before the flow reverses, its motion remains highly repeatable. This is reflected in the tight clustering of red dots and minimal fluctuations in the mean position, which also serve as the basis for estimating the phases of motion onset and cessation in § 5.3. In contrast, if the particle has not settled in time, it remains exposed after flow reversal and is likely to remain displaced from its equilibrium position in the subsequent half-cycle, resulting in a trajectory that is not exactly predictable. Additionally, a slight asymmetry in the set-up, as discussed in § 3, leads to asymmetric wiggling, with motion predominantly directed towards one side. While this asymmetry does not notably affect the overall repeatability of the motion, it introduces a small directional preference in the particle’s trajectory.
Particle position during ten consecutive oscillation cycles with constant flow amplitude, projected onto two oscillation cycles, from an experiment with
$s=1.09$
and
$\delta =0.127$
. The spatially averaged flow velocity (blue dots) and its sinusoidal fit (blue curve) are shown alongside the exposed particle position (red dots) and its mean trajectory over the ten cycles (black curve).
To further quantitatively investigate the nature of the particle motion, DNS run
$4$
replicates the experimental conditions, with three main objectives:
$(i)$
examine contact-point dynamics to verify that the mobile particle does not slide;
$(ii)$
characterise the local flow field around the mobile particle and
$(iii)$
quantify the hydrodynamic force acting on the particle.
Topic et al. (Reference Topic, Agudo, Luzi, Czech and Wierschem2022) showed that a spherical particle on a rough substrate may either roll, maintaining continuous contact with the substrate, or slide, if the tangential force at the contact point exceeds the Coulomb static friction threshold. In the former case, the minimum torque required to initiate rotation can be predicted based on the hydrodynamic torque, as described in § 2.3. In the latter case, including sliding, the particle rotation angle no longer matches the angular displacement of its centre, invalidating the model described in § 2.3.
Figure 7(a–c) presents the contact forces (split into normal and tangential components) at the four contact points of the mobile particle during the early motion stages. The three simulations have constant parameters
${\textit{Re}}_\delta =164$
and
$\delta =0.12$
, while the mobility number
$\psi$
increases, corresponding to a stationary case (
$\psi =0.99$
, panels
$a$
and
$d$
), a case with wiggling motion (
$\psi =1.1$
, panels
$b$
and
$e$
) and a case where the particle rolls over the substrate (
$\psi =1.46$
, panels
$c$
and
$f$
).
Force and motion characteristics for three representative cases from run 4 (with
${\textit{Re}}_\delta =164$
and
$\delta =0.12$
): a stationary particle at
$\psi =0.99$
(
$\varGamma =0.10$
, panels
$a$
and
$d$
), a wiggling particle at
$\psi =1.1$
(
$\varGamma =0.11$
, panels
$b$
and
$e$
) and a rolling particle at
$\psi =1.46$
(
$\varGamma =0.15$
, panels
$c$
and
$f$
).
$(a{,}b{,}c)$
Normal (
$F_{C\perp }$
, black) and tangential (
$F_{C\parallel }$
, red) components of the contact force between the mobile (red) particle and particles
$j=$
#1 (
$\square$
), #2 (
$\triangledown$
), #3 (
$\bigcirc$
), #4 (
$\vartriangle$
). Forces are normalised by the particle’s submerged weight (cf. (2.12)). The inset in panel
$(d)$
defines the contact point indices.
$(d{,}e{,}f)$
Time evolution of the mobile particle’s wall-normal centre position
$z_s$
(black) and spanwise angular velocity
$\varOmega _y$
(red). In panel
$(e)$
, the dashed black curve indicates the velocity far above the substrate, while the dashed red curve represents the particle’s angular acceleration until it comes to rest in a pocket. Panel
$(f)$
includes an inset showing a longer time interval.
Initially, the contact forces at the four contact points can be distinguished. As the flow accelerates from left to right, the magnitude of the contact force at points
$j=3,\,4$
increase, while those at
$j=1,\,2$
rapidly diminish. For sufficiently large hydrodynamic torque, these forces even vanish and the particle starts to roll (cf. figure 7
b,c). Otherwise, when the hydrodynamic force is not strong enough to initiate rolling, the particle remains within its original pocket and the contact forces at
$j=1,\,2$
recover their initial values (cf. figure 7
a,d). Notably, the difference between the normal and tangential components of the contact force is set by the static friction coefficient.
Importantly, no sliding is observed at the contact points, consistent with the experimental observations. This is most clearly reflected in the smooth temporal evolution of the angular velocity in figure 7(f), which shows the particle rolling continuously from one pocket to the next. These results support the conclusion that the motion threshold is governed by a torque balance, as outlined in § 2.3.
The hydrodynamic force acting on the mobile particle is governed by a combination of viscous and advective effects, depending on the surrounding flow regime and the development of a wake downstream. As the maximum particle Reynolds number
${\textit{Re}}_D={\textit{Re}}_\delta /\delta$
(see table 2) increases, different flow behaviours emerge that may significantly influence the particle dynamics. In the viscous-dominated regime, the drag force scales linearly with the ambient flow velocity, as is the case in run
$3$
(
${\textit{Re}}_\delta =41$
,
$\delta =0.12$
,
$\varGamma =0.45$
), where the flow remains laminar and the boundary layers developing on both the substrate and the mobile particle stay attached throughout the oscillation cycle (cf. figure 8
a,b). As inertial effects become more pronounced, such as in run
$6$
(
${\textit{Re}}_\delta =164$
,
$\delta =0.96$
,
$\varGamma =0.019$
), the flow separates downstream of the mobile particle, resulting in the formation of a wake that, while not strictly steady due to the transient nature of the oscillating flow, remains coherent without shedding vortices (cf. figures 8
c,d and 9
a). In this case, the particle size is similar to the Stokes boundary layer thickness (
$\delta \approx 1$
), and the vicinity of the bottom significantly affects the resulting hydrodynamic force acting on the particle. At higher particle Reynolds numbers
${\textit{Re}}_D$
, advective effects dominate. In run
$4$
(
${\textit{Re}}_\delta =164$
,
$\delta =0.12$
,
$\varGamma =0.10$
), the wake becomes unsteady, leading to a transitional regime where vortices are randomly shed downstream of the mobile particle (cf. figures 8
e,f and 9
b). To account for these advective effects in the torque balance model described in § 2 and Appendix A, the Schiller–Naumann correction (A13) is applied to the hydrodynamic drag.
Contour lines of the spanwise vorticity component in the vertical symmetry plane crossing the centre of the exposed particle (yellow) for:
$(a{,}b)$
run 3 (
$\delta =0.12$
,
${\textit{Re}}_\delta =41$
,
${\textit{Re}}_D=340$
,
$\varGamma =0.45$
,
$\psi =1.1$
);
$(c{,}d)$
run 6 (
$\delta =0.96$
,
${\textit{Re}}_\delta =164$
,
${\textit{Re}}_D=170$
,
$\varGamma =0.019$
,
$\psi =1.5$
) and
$(e{,}f)$
run 4 (
$\delta =0.12$
,
${\textit{Re}}_\delta =164$
,
${\textit{Re}}_D=1400$
,
$\varGamma =0.10$
,
$\psi =0.98$
). The snapshots correspond to the phases of
$(a{,}c{,}e)$
flow reversal and
$(b{,}d{,}f)$
maximum velocity far above the substrate. Contours represent positive (black) and negative (red) values of vorticity normalised by the oscillation frequency, with contour levels spaced by
$(a{,}b)$
$0.05$
,
$(c{,}d)$
$0.4$
and
$(e{,}f)$
$0.3$
. The zero contours are omitted for clarity.
Vortex structures visualised by isosurfaces of
$\lambda _2=-0.3\,U_0^{2}/\delta ^{2}$
, where
$\lambda _2$
is the second-largest eigenvalue of a tensor derived from the velocity gradient, used to distinguish vortical regions from purely straining ones (Jeong & Hussain Reference Jeong and Hussain1995). The visualisations show the flow shortly before the black particle starts to roll, for (a) run
$6$
(
$\delta = 0.96$
,
${\textit{Re}}_\delta =164$
,
$\psi =1.53$
,
$\varGamma =0.019$
) and (b) run
$4$
(
$\delta = 0.12$
,
${\textit{Re}}_\delta =164$
,
$\psi =1.1$
,
$\varGamma =0.11$
). These runs illustrate two distinct wake regimes despite having the same value of
${\textit{Re}}_\delta$
: (a) a coherent wake without vortex shedding (mean flow directed leftward) and (b) a transitional regime with unsteady wake dynamics and vortex shedding (mean flow directed rightward).
5.2. Characterisation of the parameter space
The onset and modality of the early particle motion (whether the particle remains stationary, wiggles within its pocket or rolls over the substrate) are predicted for each experimental and numerical case using the torque-balance criterion (2.31) presented in § 2 and fully derived in Appendix A. First, figure 10 maps these cases in the parameter space
$(\psi ,\delta ,{\textit{Re}}_\delta )$
. In the experiments, the oscillation frequency is held constant while the velocity amplitude is gradually increased, such that both
${\textit{Re}}_\delta$
and
$\psi$
increase simultaneously (figure 10
b). Circles indicate the cases where the particle is never set into motion throughout the oscillation cycle, while crosses and squares mark wiggling and rolling motion, respectively. In the experiments, the values of
$\psi$
separating the threshold conditions between stationary, wiggling and rolling are approximately equispaced by
$\Delta \psi =0.1$
.
Phase space overview of particle behaviour as a function of
$\delta$
,
${\textit{Re}}_\delta$
and
$\psi$
. Crosses and squares mark wiggling and rolling cases, respectively, while circles indicate cases just below the motion threshold. Experimental results are shown in blue (
$s=1.09$
) and red (
$s=1.81$
), whereas black symbols represent simulation results for small (
$\delta \approx 0.96$
, open symbols) and large (
$\delta \approx 0.12$
, filled symbols) particles. Error bars indicate measurement uncertainties. Horizontal error bars for
$\delta$
are smaller than the symbol sizes and thus not visible. The uncertainty in
${\textit{Re}}_\delta$
is mainly due to the velocity measurement uncertainty, while the error in
$\psi$
is additionally due to the uncertainty in the density ratio
$s$
, which particularly amplifies the upper bound of
$\psi$
for the light particles (blue).
For the numerical simulations, a series of independent runs is performed at increasing values of
$\psi$
, while keeping the same values of the other parameters indicated in table 2. We find good agreement between the numerical and experimental results at
$s=1.81$
, as shown in figure 10. For these conditions, the motion threshold is estimated at
$\psi \approx 1.0$
for wiggling and
$\psi \approx 1.1$
for rolling. Wiggling thus occurs within a narrow range of values of
$\psi \approx 1.0{-}1.1$
, indicating that the transition from stationary to rolling is relatively sensitive to small changes in flow conditions. This behaviour was already observed in figure 4. It is worth stressing that the wiggling motion is peculiar to oscillatory flow and has no equivalent in steady flow conditions.
For lighter particles (
$s=1.09$
), the motion threshold is reached at smaller values of
${\textit{Re}}_\delta$
at similar values of
$\delta$
, which corresponds to weaker hydrodynamic torque. The corresponding
$\psi$
-values for wiggling and rolling motion are approximately four times lower than those of the denser particles. The data further show a significant scatter, particularly for the light particles. This is partly due to their increased sensitivity to small velocity fluctuations near the motion threshold, which can arise from minor substrate irregularities, flow asymmetries or transient effects during velocity ramp-up. As a result, the measured velocity
$U_0$
at the onset of motion has a larger uncertainty. Moreover, the upper bound of the uncertainty in
$\psi$
is amplified for light particles because
$\psi$
is inversely proportional to
$(s-1)$
. As
$s$
approaches unity (i.e. as the movable particle becomes nearly neutrally buoyant), the stabilising torque due to weight vanishes and
$\psi$
diverges. In our experiments, this means that even small variations in
$s$
of the order of the measurement precision (
$\pm 0.05$
) can lead to significant discrepancies in
$\psi$
. While
$\psi$
reflects the characteristic pressure amplitude, it does not solely determine the motion threshold. As shown in figure 10(b), the threshold occurs at different
$\psi$
-values depending on the flow regime, e.g.
$\psi \sim 0.25$
for
${\textit{Re}}_\delta \sim 10{-}40$
and
$\psi \sim 1.0$
for
${\textit{Re}}_\delta \approx 150$
. This demonstrates that the critical value of
$\psi$
is not universal for the cases considered here, but instead varies with other dimensionless parameters.
Figure 11 presents the experimental and numerical maximum values of
$\vert \varUpsilon \vert$
attained during the oscillation cycle and calculated using (2.28), and compares them to the theoretical motion threshold of
$|\varUpsilon |=1$
. The figure shows good agreement between experiments, simulations and the theoretical threshold. For dense particles (red), the obtained values of
$\vert \varUpsilon \vert$
for both wiggling and rolling regimes deviate by only
$2\,\%$
and
$4\,\%$
, respectively, from the predicted threshold at unity. These deviations fall well within the overall uncertainty in predicted
$\psi$
(approximately
$10\,\%$
). Notably, the numerically determined threshold for wiggling motion is
$\vert \varUpsilon \vert =1.02$
, which closely matches both the experimental data and the model prediction. Figure 11 further highlights that
$\varUpsilon$
is highly sensitive to small variations in the density ratio
$s$
, particularly as it approaches unity. Even a slight correction of the experimental value from
$s=1.09$
to
$s=1.04$
, within the bounds of measurement uncertainty, yields a significant shift in
$\varUpsilon$
in the diagram and substantially improves the prediction of incipient motion for these lighter particles. As with
$\psi$
, the increased scatter for the light particles is further amplified by their sensitivity to velocity fluctuations near the motion threshold.
Threshold parameter
$\varUpsilon$
as a function of (a)
$\psi$
, (b)
$\delta$
and (c)
${\textit{Re}}_\delta$
, with circles, crosses and squares as in figure 10. Experimental results are shown in red for
$s=1.81$
, and in blue for
$s=1.09$
(dark symbols) and
$s=1.04$
(light symbols). The latter corresponds to a slight modification of
$s$
to the lower bound of measurement uncertainty, highlighting the model’s sensitivity to lighter particles, due to the scaling with
$s-1$
. Simulation results are shown in black, with open symbols for small particles (
$\delta \approx 0.96$
) and filled symbols for large particles (
$\delta \approx 0.12$
). The vertical error bars represent uncertainty in
$\varUpsilon$
, primarily due to velocity measurement errors (
${\sim} 5\,\%$
) for dense particles (red), and additionally due to uncertainty in the density ratio
$s$
for light particles (blue), which particularly amplifies the upper bounds. Error bars for
$\psi$
and
${\textit{Re}}_\delta$
are omitted for clarity but are provided in figure 10.
The accuracy of our model is further confirmed by comparing the transient values of
$\varUpsilon$
over multiple oscillation cycles and identifying the instants when the particle motion starts, as shown in figure 12. Notably, in simulations where the maximum value of
$\varUpsilon$
exceeds unity, the mobile particle is consistently observed to roll. A key advantage of the
$\varUpsilon$
-criterion is that it predicts the onset of wiggling or rolling motion using a threshold value based on mechanistic understanding, while also determining the phase within each half-cycle at which the motion starts.
Computed values of
$\varUpsilon$
for:
$(a)$
run
$6$
at
$\psi =1.53,\,1.58,\,1.67,\,2.23$
;
$(b)$
run
$3$
at
$\psi =1.11$
and
$(c)$
runs
$4{-}5$
at
$\psi =0.99,\,1.11$
. Black and red curves correspond to simulations in which the mobile particle is static or starts to roll, respectively. The black triangles in panel (a) mark the phases of incipient particle motion observed in the DNS, which intercept the corresponding model prediction (red curve). Note that the red curve for
$\psi =1.58$
nearly overlaps with the black curve for
$\psi =1.53$
.
5.3. Phase information
The phases of motion initiation and cessation are determined from mean particle trajectories near the motion threshold, e.g. as shown in figure 6. The resulting phase difference
$\Delta t$
in figure 13 exhibits a linear correlation with
$\sqrt {s/K_C}$
, similar to that of the time to crest
$\tau$
predicted by (2.33), which is particularly evident for light particles. Unlike in the calculation of
$\varUpsilon$
, the results in figure 13 are much less sensitive to small variations in
$s$
, as the scaling here involves
$s$
rather than
$s-1$
. For all the present simulations, the dimensionless particle angular acceleration remained of order one (e.g. see figure 7
e for run
$4$
), leading to an equality for the time to crest
$\tau$
as given by (2.33).
Duration of particle motion within each half-cycle for conditions just above the motion threshold. Red and blue symbols correspond to experiments with relatively heavy (
$s=1.81$
) and light (
$s=1.09$
) particles, respectively. Crosses and squares indicate wiggling and rolling motion, respectively. The dashed line shows a linear fit through the data. Horizontal error bars reflect measurement uncertainty in the velocity amplitude and density ratio, while vertical error bars indicate the precision in phase measurements.
To further investigate the phase of motion initiation, we analyse cases well above the motion threshold by increasing the value of
$\psi$
in the DNS, as shown in figure 14. The experimental data for dense particles and the DNS results from run 6 (triangles) complement each other, with our model accurately predicting the initiation time observed in the DNS. As
$\psi$
increases and the condition
$|\varUpsilon |\gt 1$
is satisfied for a longer duration within the oscillation cycle, the onset of motion shifts to earlier phases. Remarkably, for sufficiently strong forcing, motion begins even before the reversal of the flow far above the substrate. This is attributed to viscous contributions within the oscillatory boundary layer, which induce a phase lead of up to
$\pi /2$
in the near-bed velocity relative to the bulk flow. Since the boundary layer is relatively thick in this case (
$\delta =0.96$
), this phase lead has a pronounced effect on the phase of the hydrodynamic force. For larger particles, but still in laminar flow conditions (run 3, blue circles), motion initiation in the DNS occurs precisely at the flow reversal, aligning with our model predictions.
(a) Initiation time
$t_{{in}}$
at which
$\vert \varUpsilon \vert$
exceeds unity and the particle begins to move. DNS results (filled symbols), experimental results (filled grey symbols) and model predictions (open symbols) are compared for parameter values corresponding to run
$4$
(red squares), run
$3$
(blue circles) and run
$6$
(black triangles). The inset shows the initiation time projected onto part of the oscillation cycle centred around the flow reversal, highlighting that for some cases, initiation occurs at flow reversal (blue circles) or even before (leftmost black triangles). The experimental measurement uncertainties for
$\psi$
and
$t_{{in}}$
are previously shown in figures 10 and 13, respectively. (b) Modelled time interval during which
$\vert \varUpsilon \vert \gt 1$
. The threshold values separating wiggling (filled symbols) from rolling particles (open symbols) are indicated by the horizontal red and black dashed lines, as predicted by the time to crest scaling (2.33) for runs
$4$
and
$6$
, respectively.
Comparing the type of particle motion with the duration
$\Delta t_{|\varUpsilon |\,\gt\, 1}$
for which
$\vert \varUpsilon \vert \gt 1$
(figure 14
b), we find that this interval generally exceeds the characteristic time scale
$\tau$
from (2.33) which separates the wiggling and rolling regimes. The only case where wiggling is observed corresponds to
$\Delta t_{| \Upsilon |\gt1}\lt \tau$
, supporting the relevance of
$\tau$
as the critical time scale for a particle to roll out of a pocket in the bed. Additionally, the sharp increase in
$\Delta t_{| \Upsilon |\gt1}$
with particle mobility (i.e. increasing
$\psi$
in figure 14
b) explains the relatively abrupt transition from static to rolling particles, with only a narrow regime of wiggling in between, as previously observed in figure 10. Finally, the scaling
$\tau \sim \delta ^{-0.5}$
predicts that
$\tau$
asymptotically vanishes for small particles (
$\delta \gg 1$
), which is also the reason why it is challenging to observe the wiggling motion in the case of small particles (run 6, black triangles in figure 14
b). The wiggling motion state is thus typically characteristic of large particles (e.g. run 3, red squares in figure 14
b).
For these cases with relatively large particles or fast flow reversal (small
$\delta$
), distinguishing between wiggling and rolling requires accounting for the duration over which the force exceeds the threshold (i.e.
$\Delta t_{| \Upsilon |\gt1}$
). This suggests that combining
$\varUpsilon$
with a characteristic time scale into an angular impulse-based criterion could yield predictive value even during early stages of particle motion. This approach aligns with the impulse-based framework of Diplas et al. (Reference Diplas, Dancey, Celik, Valyrakis, Greer and Akar2008) and Valyrakis et al. (Reference Valyrakis, Diplas, Dancey, Greer and Celik2010), where the product of force and duration is used as a predictor for particle dislodgement in turbulent conditions.
5.4. Comparison with other criteria
The method presented in this study enables the deterministic prediction of the initiation of rolling motion for a single sphere exposed to an oscillatory flow. This is conceptually different from predicting the incipient motion of the surface layer of a horizontal sediment layer. In the latter case, the arrangement of the particles in the substrate, as well as the disturbance in the flow due to neighbouring mobile particles and their contacts, introduce randomness and additional resistance. Hence, criteria based on the balance of horizontal forces acting on the mobile layer, such as the Shields criterion (balancing drag and bottom friction), are expected to fail to predict the incipient motion of individual sediment particles. Moreover, for steady flows, the torque balance criterion (2.28) can be reformulated as a critical Shields parameter by expressing it in terms of a critical particle density and substituting this into the general Shields expression (2.22), as shown by Agudo et al. (Reference Agudo, Illigmann, Luzi, Laukart, Delgado and Wierschem2017). However, this approach is not practical for unsteady flows, where it leads to both the criterion and the Shields parameter becoming time-dependent, complicating their interpretation and application. Only in the limit of
$\delta \gg 1$
and
${\textit{Re}}_D\lt 1$
, our model reduces to a formulation closely aligned with the Shields criterion, as derived for uniform shear flows by Agudo et al. (Reference Agudo, Illigmann, Luzi, Laukart, Delgado and Wierschem2017). In an OBL, the imposed oscillatory pressure gradient by itself can lead to the initiation of particle motion, making the Shields parameter
$\theta$
, as defined in (2.22), incomplete since it only accounts for the contribution of drag. Figure 15 compares incipient motion predictions based on
$\theta$
and
$\varUpsilon$
(as computed from (2.28)) for runs
$6$
and
$4$
. The critical Shields parameter
$\theta _{cr}$
was estimated using the ‘lower-limit’ formula obtained by Paphitis (Reference Paphitis2001), based on experimental unidirectional flow data. Paphitis (Reference Paphitis2001) adopted a statistical approach to identify the critical conditions for the transport of quartz particles and proposed lower- and upper-limit curves that enveloped a large heterogeneous dataset. In the present work, we compare with the ‘lower-limit’ value because our substrate geometry, a closely packed square arrangement, supports early onset of particle motion with respect to other geometries. While both quantities relate to the onset of motion, their time evolution differs significantly, especially for run 4, making it difficult to infer the exact phase of incipient rolling motion directly from
$\theta$
. More importantly, the threshold
$\theta _{cr}$
substantially overestimates the actual onset of motion by up to a factor of five, as indicated by the difference between the black curve and the black dashed threshold line in figure 15. Thus, it mispredicts the actual conditions for particle transport in unsteady flows. Beyond the Shields parameter, Frank et al. (Reference Frank, Foster, Sou, Calantoni and Chou2015) proposed to consider also the Sleath parameter, which normalises the force due to the imposed pressure gradient by the particle weight, to predict incipient sediment transport within an OBL with relatively small
$\delta$
. However, none of the models proposed so far accounted for the variation in the velocity distribution during the oscillation cycle and the consequent effects on the torque. The present mechanistic approach additionally eliminates the need to determine empirical threshold values for different criteria, which typically depend on flow conditions.
Time evolution of
$\varUpsilon$
(red curve) as computed from (2.28) and the Shields parameter
$\theta$
(black curve) for cases just above the motion threshold in: (a) run
$6$
and (b) run
$4$
. The horizontal dashed lines indicate the respective critical values, where
$\theta _{cr}$
is estimated using the ‘lower-limit’ formula proposed by Paphitis (Reference Paphitis2001).
6. Concluding remarks
The threshold and characteristics of particle motion on top of a structured substrate driven by hydrodynamics have been studied, focusing on a uniform monolayer of spheres arranged in a square grid with a single mobile particle on top. While prior results have primarily focused on steady unidirectional flows, this work concentrates on the added complexity introduced by oscillatory flows. The problem is subsequently governed by an additional degree of freedom: the oscillation frequency introduces a time scale that enters the problem through the frequency-dependent viscous length scale
$\delta$
, corresponding to the normalised thickness of the oscillatory boundary layer.
Our model for the motion threshold, based on a balance of torques, provides excellent quantitative agreement with both our experimental and numerical data, with time-averaged relative errors consistently of the order of a few per cent (see figure 18). The model effectively captures the threshold differences due to variations in particle size relative to the viscous length scale, which sets the relative importance of the pressure gradient and hydrodynamic drag. Notably, the model performs well despite relying solely on theoretical velocity profiles and derived forces in the streamwise direction, with only two tuneable constants. Our mechanistic approach provides an accurate estimate of the motion threshold without requiring consideration of the shear stress at the substrate, particularly for relatively dense particles. In these cases, the observed values of
$\varUpsilon$
deviate by only
$2\,\%{-}4\,\%$
(for wiggling and rolling, respectively) from the theoretically predicted threshold, which is well within the experimental uncertainty (
${\sim} 10\,\%$
). In other words, it replaces the need for empirically determined threshold values at different flow conditions with a criterion based on the mechanistic balance of stabilising and destabilising torques.
Comparison of streamwise velocity profiles at different phases of the oscillation cycle from DNS results for: (a) run 6 (
$\delta =0.96$
,
${\textit{Re}}_\delta =164$
); (b) run
$3$
(
$\delta =0.12$
,
${\textit{Re}}_\delta =41$
) and (c) run
$4$
(
$\delta =0.12$
,
${\textit{Re}}_\delta =164$
). All velocity profiles are plotted relative to the virtual wall elevation
$z_0$
(black curves) and compared with the analytical Stokes solution shifted vertically by
$z_0$
(red curves). Black dashed lines mark the crest elevation of the fixed substrate particles, located at
$z=z_0+0.06\,D$
. The lower red dashed lines mark the base elevation of the mobile particle, expressed as
$(z_b-z_0)/D$
, which equals
$-0.14$
(run
$6$
),
$-0.18$
(run
$3$
) and
$-0.16$
(run
$4$
). The upper red dashed lines indicate the exposure height
$h$
of the mobile particle above the substrate. In panel (c), symbols represent experimental PIV data at five oscillation phases, shown for comparison with the numerical results. Interrogation windows of
$32\times 8$
pixels with 50 % overlap were used to enhance vertical resolution. Note that the velocity amplitudes far above the substrate and the vertical shift in the decay towards zero velocity are reasonably captured in the experimental data. However, the detailed boundary layer profile near the bed is not resolved, likely due to smoothing and averaging errors caused by scattering and reflections.
Temporal evolution of the dimensionless functions governing the hydrodynamic torque: the lever arm
$f(\omega t,H/\delta )$
, the zeroth-order moment of the velocity profile
$g(\omega t,H/\delta )$
, and the contributions from added mass and imposed pressure gradient
$h(\omega t,\zeta )$
. Note that
$f(\omega t,H/\delta )$
diverges at the zero crossings of
$g(\omega t,H/\delta )$
, corresponding to the flow reversal, where the instantaneous drag force vanishes.
Comparison between model predictions (red) and DNS results (black) for: (a) run
$6$
(
${\textit{Re}}_\delta =164$
,
$\delta =0.96$
); (b) run
$3$
(
${\textit{Re}}_\delta =41$
,
$\delta =0.12$
) and (c) run
$4$
(
${\textit{Re}}_\delta =164$
,
$\delta =0.12$
). Solid and dashed lines represent the zeroth velocity moment
$u_0$
and the drag lever arm
$L_D$
, respectively.
Above the threshold, the particle motion in oscillatory flows is richer than in steady flows, as explored by Agudo & Wierschem (Reference Agudo and Wierschem2012) and Topic et al. (Reference Topic, Agudo, Luzi, Czech and Wierschem2022). In the steady case, only the motion threshold itself is relevant: once exceeded, a positive feedback loop between hydrodynamic drag and exposure to the flow starts (Agudo et al. Reference Agudo, Illigmann, Luzi, Laukart, Delgado and Wierschem2017). In contrast, the additional degree of freedom in the oscillatory case adds a dynamic complexity to the particle behaviour, emerging as a competition between the particle inertia as it rolls towards a neighbouring pocket and the flow reversal time. When the reversal occurs too quickly, the particle may move but fail to reach the next pocket, resulting in a wiggling state that is not observed in the steady case.
Our findings suggest that the present approach can be extended to various flow and substrate conditions, e.g. higher Reynolds numbers, and even to flows that are not strictly laminar but in the transitional regime, provided the relevant length scale characterising the boundary layer can be accurately predicted. In such cases, the pressure gradient has a dominant contribution to the destabilising torque, for which the model accurately accounts. However, the predictive power reaches its limits when hydrodynamic fluctuations become significant, such as those induced by chaotic vortex shedding in turbulent conditions. These fluctuations introduce variability in the force and torque acting on the particle, making its motion inherently unpredictable.
Finally, our work provides extensive opportunities to further explore incipient particle motion in unsteady flows and extend our modelling approach. The effect of bottom geometry warrants further study: while its role has been examined in unidirectional flows (Agudo & Wierschem Reference Agudo and Wierschem2012; Agudo et al. Reference Agudo, Illigmann, Luzi, Laukart, Delgado and Wierschem2017; Topic et al. Reference Topic, Retzepoglu, Wensing, Illigmann, Luzi, Agudo and Wierschem2019), its effect on particle dynamics in oscillatory flows remains largely unexplored. An especially intriguing aspect is the potential for cross-stream particle transport relevant to the mixing and dispersion processes of both sediment grains and microplastics. The combination of symmetrical forcing, which on average eliminates net streamwise transport, and anisotropy in the substrate geometry could lead to lateral transport and particle–particle interactions. This mechanism may contribute to the formation of sediment bedforms, such as ripples, alongside other self-organisation processes like steady streaming flows (Klotsa et al. Reference Klotsa, Swift, Bowley and King2009; Mazzuoli et al. Reference Mazzuoli, Kidanemariam, Blondeaux, Vittori and Uhlmann2016; van Overveld et al. Reference van Overveld, Clercx and Duran-Matute2023).
Funding
Marco Mazzuoli was supported under the subaward no. SUB00004180 of the University of Florida, Gainesville, FL, USA (PTE Federal award no. N00173-21-2-C900). The simulations were partly performed using the CINECA HPC facilities under the ISCRA-b project SEAWAVES no. HP10BPMJZR.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are (soon) openly available in 4TU.ResearchData.
Appendix A. Analytical model for the torque balance
This appendix presents the full derivation of the analytical model describing the torque balance on the mobile spherical particle, evaluated about its contact points. To avoid ambiguity, all quantities are expressed in dimensional form. The model distinguishes between stabilising and destabilising contributions. The onset of particle motion is assessed by evaluating the ratio of these contributions, yielding a criterion for incipient rolling. The destabilising torques are associated hydrodynamic forces such as drag, lift, added mass and the imposed pressure gradient, the latter two being related to the flow unsteadiness. The stabilising torque is primarily due to the submerged weight of the particle, which helps preserve its stable position.
A.1. Drag and lift force models
First, we analyse the hydrodynamic drag and lift forces acting on the mobile particle, which are expressed as
\begin{align} F_D &= {\displaystyle \int }_{{\mathcal S}} (\mathbb{T}\boldsymbol{\cdot }{\boldsymbol{n}})\boldsymbol{\cdot }{\boldsymbol{e}}_x \,{\rm d}S = \pi D{\displaystyle \int }_{z_b}^{z_b+D} f_x \,{\rm d}z, \end{align}
\begin{align} F_L &= {\displaystyle \int }_{{\mathcal S}} (\mathbb{T}\boldsymbol{\cdot }{\boldsymbol{n}})\boldsymbol{\cdot }{\boldsymbol{e}}_z \,{\rm d}S = \pi D{\displaystyle \int }_{x_c-D/2}^{x_c+D/2} f_z \,{\rm d}x, \end{align}
where
$\mathbb{T}$
is the fluid stress tensor,
$x_c$
the streamwise coordinate of the particle centre and
$f_x$
(
$f_z$
) the streamwise (wall-normal) component of the hydrodynamic stress acting on vertical (horizontal) strips of infinitesimal thickness
${\rm d}z$
(
${\rm d}x$
), respectively. We further define the
$n$
th moment
$u_n$
of the ambient streamwise velocity
$\langle u\rangle (z)$
relative to the elevation
$z=z_D$
of the contact point, as
\begin{align} u_n = \dfrac {1}{D}{\displaystyle \int }_{z_b}^{z_b+D} (z-z_D)^n \langle u\rangle \,{\rm d}z. \end{align}
Assuming
$f_x$
is proportional to the local undisturbed velocity, the drag force and its torque about
$z=z_D$
follow as
\begin{align} F_D &\propto \pi D{\displaystyle \int }_{z_b}^{z_b+D} \langle u\rangle\, {\rm d}z = \pi D^2u_0, \end{align}
\begin{align} T_D &\propto \pi D{\displaystyle \int }_{z_b}^{z_b+D} (z-z_D)\langle u\rangle\, {\rm d}z = \pi D^2 u_1, \end{align}
where define the drag lever arm
$L_D$
as
\begin{align} L_D=\dfrac {T_D}{F_D} = \dfrac {u_1}{u_0}. \end{align}
In the case of a uniform ambient flow, the non-uniformity of
$f_x$
over the particle height is solely due to the varying orientation (altitude angle) of the local surface element, since the strip surface remains constant at
$\pi D {\rm d}z$
. However, in more complex flows, particularly if convective transport dominates (e.g. at moderate-to-high Reynolds numbers or turbulent wake conditions), the actual value of
$L_D$
may deviate from (A6).
Similarly, the lift lever arm
$L_L$
is defined by
\begin{align} L_L=\dfrac {T_L}{F_L}, \end{align}
where the lift torque about the particle centre is given by
\begin{align} T_L &= \pi D{\displaystyle \int }_{x_c-D/2}^{x_c+D/2} f_z(x-x_D) \,{\rm d}x. \end{align}
For a spherical particle resting on a flat wall (
$z_D=z_b=0$
) exposed to a linear shear flow, the zeroth velocity moment reads
\begin{align} u_0 &= \dfrac {1}{D}\int _0^{D} 2\dfrac {u_c}{D}z\, {\rm d}z = u_c, \end{align}
where
$u_c$
denotes the ambient flow velocity at the particle centre elevation. This quantity plays a central role in the model, as it allows generalisation to unsteady flow conditions beyond the canonical oscillatory boundary layer (OBL). In the low-Reynolds-number limit, the drag force reduces to Stokes’ drag,
\begin{align} F_D = 3\pi \rho _{\!f}\nu u_0D = 3\pi \rho _{\!f}\nu u_cD, \end{align}
which is consistent with the assumption (A4). The corresponding lever arm equals
$L_D=(2/3)D$
.
For more complex flow conditions, where
$u_0$
deviates from
$u_c$
, empirical corrections are introduced. Following Agudo et al. (Reference Agudo, Illigmann, Luzi, Laukart, Delgado and Wierschem2017) for a particle resting on a rough, structured substrate, corrections account for: (i) the effective elevation of the particle (
$z_D-z_b$
) and (ii) the shielding effect of the substrate (e.g. related to
$z_0-z_b$
). The normalised geometrical parameters
$H=h/D$
and
$L=l/D$
, from (2.4) and (2.10), capture these effects. Additionally, a correction coefficient
$C^\prime$
incorporates effects due to finite-size and nonlinear shear flows, such that the drag force is modified to
\begin{align} F_D = 3\pi C^\prime \rho _{\!f}\nu u_0D, \end{align}
where
$C^\prime$
depends on the instantaneous particle Reynolds number
\begin{align} {\textit{Re}}_s(t) = \dfrac {u_cD}{\nu } \end{align}
and the relative boundary layer thickness
$\delta$
. For
${\textit{Re}}_s\lesssim 1000$
, the Schiller–Naumann correction (Schiller Reference Schiller1933) is used:
\begin{align} C'= 1 + 0.15{\textit{Re}}_s^{0.687}. \end{align}
For
$\delta \gt 1$
, the correction proposed by Zeng et al. (Reference Zeng, Najjar, Balachandar and Fischer2009) applies for
${\textit{Re}}_s\lt 200$
:
\begin{align} C'= 1.7005\big (1 + 0.104{\textit{Re}}_s^{0.753}\big ). \end{align}
In the present simulations, (A13) is applied for cases with
$\delta =0.12$
(runs 1–5) and (A14) for run
$6$
(see table 2). The Schiller–Naumann correction generally holds for uniform flows and is not strictly applicable in shear-dominated flows. However, for the cases where (A13) is used, the small value of
$\delta =0.12$
implies that the velocity profile experienced by the mobile particle is nearly uniform, making the correction a reasonable approximation. Near-wall effects under such conditions (small
$\delta$
) can be compared with the case of a particle moving parallel to a wall in a quiescent fluid (Zeng et al. Reference Zeng, Balachandar and Fischer2005, Reference Zeng, Najjar, Balachandar and Fischer2009), where the drag coefficient diverges as the particle nears the wall, while the Reynolds dependence diminishes. Since it is not directly applicable to our system due to the more complex bottom geometry and flow conditions, we will explicitly validate the modelled drag force against the present DNS results.
We now turn to the lift force. For small particles, we adopt the Saffman lift force expression (Saffman Reference Saffman1965)
\begin{align} F_L = 3.22 C^{\prime \prime } D \rho _{\!f}\sqrt {\frac {\nu D}{2\vert u_c\vert }} u_c^2, \end{align}
where
$C^{\prime \prime }=1$
for
$\delta \gtrsim 1$
. For
${\textit{Re}}_s\gt 40$
, we use the correction proposed by McLaughlin (Reference McLaughlin1991) and Mei (Reference Mei1992):
\begin{align} C^{\prime \prime } = 0.0741\sqrt {{\textit{Re}}_s}. \end{align}
To prevent overprediction at high values of
${\textit{Re}}_s$
, we cap
$C^{\prime \prime }$
at the Saffman value (
$C^{\prime \prime }\!=1$
) when
${\textit{Re}}_s\gtrsim 1100$
, namely when advective effects on the particle dynamics became dominant (see Appendix A.2). A similar cut-off was also predicted by Schiller (Reference Schiller1933) for drag at high
${\textit{Re}}_s$
in steady flow conditions, though it is not observed here for the OBL for the parameter values presently considered.
Since the lift force distribution
$f_z(x)$
along the streamwise coordinate is not directly known, the lift lever arm
$L_L$
is estimated from our DNS results reported in the following section. Its validity is thus specific to the OBL under conditions close to ours.
A.2. Torque balance in an oscillatory boundary layer
We now examine the torque acting on a particle exposed to OBL flow, in the configuration shown in figure 1.
Assuming that the mean velocity
$\langle u\rangle$
is reasonably approximated by the Stokes solution (2.3), referenced from the effective zero level
$z_0$
, each torque contribution can be expressed in closed form using the previously introduced force models. Figure 16 supports this assumption, showing good agreement between DNS velocity profiles and the vertically shifted Stokes solution, particularly for low to moderate
${\textit{Re}}_\delta$
. At higher values
${\textit{Re}}_D={\textit{Re}}_\delta /D$
(e.g. run 4), deviations remain moderate and can be compensated by adjusting the effective viscous length scale. Corrections for such advective effects associated with bottom roughness at moderate Reynolds numbers are further discussed in Appendix A.3.
The instantaneous particle Reynolds number is now given by
\begin{align} {\textit{Re}}_s(t) = \frac {{\textit{Re}}_\delta }{\delta }\frac {u_c}{U_0} = \frac {{\textit{Re}}_\delta }{\delta }\left [\sin (t)-e^{-\zeta }\sin \left (t-\zeta \right )\right ]\!, \end{align}
where
$\zeta =(H-1/2)/\delta$
denotes the dimensionless elevation of the particle centre above the zero level of the velocity profile
$z_0$
.
The zeroth moment of the velocity across the particle height is
\begin{align} \dfrac {u_0(t,H,\delta )}{U_0} = \dfrac {1}{U_0D}\int _{z_b}^{z_b+D} u\, {\rm d}z = \dfrac {1}{U_0D}\int _{z_0}^{z_b+D} u\, {\rm d}z = H\,g\left (\omega t,\frac {H}{\delta }\right )\!, \end{align}
with
$H=(z_b+D-z_0)/D$
the dimensionless exposure height of the particle. The function
$g(\omega t,H/\delta )$
, plotted in figure 17(b), quantifies the effect of shielding by the substrate and is defined as
\begin{align} g\left (\omega t,\frac {H}{\delta }\right ) = \sin (\omega t)-\frac {\sqrt {2}}{2}\frac {\delta }{H}\left [\sin \left (\omega t-\frac {\pi }{4}\right )- e^{-H/\delta }\sin \left (\omega t - \frac {H}{\delta } - \frac {\pi }{4}\right )\right ]\!. \end{align}
In the limit
$H/\delta \gg 1$
(large particle),
$u_0/U_0\approx H\sin (\omega t)$
, reflecting a quasi-uniform profile. In contrast, for
$H/\delta \ll 1$
(small particle), the expression reduces to
$u_0/U_0\approx (\sqrt {2}/2 ) (H^2/\delta )\sin (\omega t+\pi /4 )$
, which is maximum during the phases of maximum bottom shear stress. The latter case further shows the torque being dominated by bottom shear stress rather than imposed pressure gradient.
With the drag force (A11) now reading
\begin{align} F_D = 3\pi C'\rho _{\!f}\nu U_0DH\,g\left (\omega t,\frac {H}{\delta }\right )\!, \end{align}
the hydrodynamic torque associated with drag is given by the product of the drag force
$F_D$
and its effective lever arm
$L_D$
, illustrated in figure 1 and approximated by (A6). Using the Stokes boundary layer profile (2.3), the normalised lever arm can be written as
\begin{align} \frac {L_D}{D}=\frac {u_1}{D\,u_0} = H\left [f\left (\omega t,\frac {H}{\delta }\right )-1\right ] + L + \frac {1}{2}, \end{align}
where
$L=l/D = 1/2 + (z_b-z_D )/D$
accounts for the vertical offset between the particle centre and the contact points, effectively incorporating the influence of the substrate geometry (cf. figure 1
b).
The time-dependent function
$f(\omega t,H/\delta )$
, shown in figure 17(a), describes the first moment of the velocity distribution and is given by
\begin{align}& f\left (\omega t,\frac {H}{\delta }\right ) \nonumber\\& =\dfrac { \sin (\omega t) + \left (\!\dfrac {\delta }{H}\!\right )^2 \cos\, (\omega t ) +\dfrac {\delta }{H}{\textrm{e}}^{-H/\delta }\left [\!\sqrt {2}\sin \left (\!\omega t-\dfrac {H}{\delta }-\dfrac {\pi }{4}\!\right )-\dfrac {\delta }{H}\cos \!\left (\!\omega t-\dfrac {H}{\delta }\!\right )\!\right ] }{ 2 g\left (\!\omega t,\dfrac {H}{\delta }\!\right ) } . \end{align}
In the limiting case of a particle resting directly on a flat surface (
$L=1/2$
and
$H=1$
), the lever arm simplifies to
$L_D/D=f(\omega t, H/\delta )$
. For a substrate of spherical particles in a square arrangement, we obtain
$L=\sqrt {2}/4$
. In the large particle regime (
$H/\delta \gg 1$
), the lever arm reduces to
$L_D=1/2$
, while in the small particle regime (
$H/\delta \ll 1$
), it tends to
$L_D\approx 2/3$
, consistent with predictions for steady linear shear flow (Agudo et al. Reference Agudo, Illigmann, Luzi, Laukart, Delgado and Wierschem2017). Thus, the torque associated with the drag force, normalised according to (2.12), reads
\begin{align} \frac {T_{\textit{Sy}}^{\textit{drag}}}{\tau _{\nu 0}\pi D^3/2} = 6 C'H\delta \left [\left (f\left (\omega t,\frac {H}{\delta }\right )-1\right )H+L+\frac {1}{2}\right ]g\left (\omega t,\frac {H}{\delta }\right )\!, \end{align}
where the terms inside the brackets correspond to the time-dependent lever arm
$L_D$
defined in (A21).
The analytical predictions for the zeroth moment of the velocity
$u_0$
and the drag lever arm
$L_D$
show excellent agreement with DNS results from runs
$3$
,
$4$
and
$6$
(figure 18), with time-averaged relative errors of 3.8 % and 1.2 % (figure 18
a), 2.5 % and 3.9 % (figure 18
b), and 8.0 % and 5.7 % (figure 18
c) for
$u_0$
and
$L_D$
, respectively. This confirms the accuracy of the expressions derived from the Stokes profile referenced to an effective zero level. Furthermore, figure 19 shows that the same model accurately reproduces the time-resolved drag and lift forces for both small and large particle cases (
$\delta =0.96$
and
$\delta =0.12$
, respectively). This close agreement between model and DNS results across a range of flow conditions supports the validity of the force expressions and associated torque calculations within the present parameter space. In particular, the good agreement between the modelled force contributions and the DNS results in figure 19 lends further confidence to our chosen expressions for the correction factors
$C'$
and
$C''$
. Furthermore, the overall expression for
$\varUpsilon$
(A28) is not particularly sensitive to them, with only the drag and lift terms scaling linearly with these coefficients. The added mass and pressure gradient contributions (especially relevant for larger particles, as shown in figure 2) are independent of
$C'$
and
$C''$
, thereby reducing their influence on the final result. Although the lift force model does not explicitly account for substrate effects, these are partially incorporated via the torque lever arm derived from DNS data. The agreement between model and DNS is generally strong (e.g. run 6 in figure 19
a), but minor deviations for cases with thinner boundary layers (e.g. run 4 in figure 19
b) suggest that substrate effects not captured by the model may play a role. A more refined treatment, e.g. based on particles in wall-bounded shear flow (Zeng et al. Reference Zeng, Najjar, Balachandar and Fischer2009), could enhance predictive accuracy in future extensions.
Comparison between model predictions (red) and DNS results (black) of the hydrodynamic forces acting on the particle for: (a) run 6 and (b) run 4. Solid and dashed lines indicate drag and lift force components, respectively. In run
$6$
, the Stokes profile is referenced from a virtual origin located at
$z_0=z_b+0.15\,D$
.
In addition to viscous drag, the unsteady nature of the oscillatory flow gives rise to additional hydrodynamic forces, namely those due to the imposed external pressure gradients and the added mass effect. These also contribute to the horizontal torque on the particle. The torques associated with these forces are expressed by
\begin{align} T_{\boldsymbol{\nabla}\kern-1pt p} &= F_{\boldsymbol{\nabla}\kern-1pt p} L_{\boldsymbol{\nabla}\kern-1pt p} = \rho _{\!f} V_s U_0 \omega \cos (\omega t)\,L, \end{align}
\begin{align} T_{\textit{AM}} &\approx \frac {1}{2}\rho _{\!f}V_s\frac {\partial u_c}{\partial t} L, \end{align}
respectively.
The torques associated with forces in the vertical direction include the gravitational torque, depending on the submerged weight of the particle, and the lift force. The corresponding lever arm associated with the particle’s weight is
$L_{{w}}=D/4$
, while the lever arm of the lift force,
$L_L$
, was obtained directly from the DNS. In all cases, the lift-force lever arm was found to be smaller than
$L_w$
. During phases of the oscillation cycle when the lift force is significant,
$L_L\simeq 0.62\,L_w$
corresponding to a dimensionless value
$L_L/D=0.155$
(see figure 20). This reduced value indicates that the lift force acts further from the particle centroid than the gravitational force, which is reflected in the net torque balance. In fact, during these phases, the stress distribution over the sphere surface is asymmetric between the upstream and downstream sides. As a result, the net force is shifted downstream, as illustrated in figures 20(a) and 20(c), which correspond to relatively small and large particles, respectively.
Spatial distribution of the vertical force density
$f_z$
(defined in (A2)) along the normalised horizontal coordinate for runs: (a) 6 and (c) 4, at the phase of maximum velocity far above the substrate. The vertical red dashed lines mark the location of the resultant lift force, as defined by (A7). (b,d) Temporal evolution of the lift lever arm
$L_L$
normalised by the weight lever arm
$L_W$
(red) and the lift force
$F_{\textit{lift}}$
normalised by the submerged weight
${\mathcal{W}}_{\mathcal{\int }}$
(black), for the same runs. The horizontal red dashed lines indicate the phase-averaged value of the normalised lever arm during the phases of maximum lift force.
Consequently, the torques due to the submerged weight and the lift force are given by
\begin{align} T_{\textit{By}} &= - \frac {{\mathcal W}_sD}{2} \:, \end{align}
\begin{align} T_{\textit{Sy}}^{\textit{lift}} &= 3.22 C'' \rho _{\!f}\nu ^2\sqrt {{\textit{Re}}_s^3} L_L, \end{align}
where
${\mathcal W}_s$
is the submerged weight of the particle.
Finally, the ratio
$\varUpsilon$
between destabilising torques (due to hydrodynamic drag, pressure gradient, added mass and lift force) and the stabilising torque (due to the particle’s submerged weight) is defined as
\begin{align} \varUpsilon &\equiv \dfrac {T_{\textit{Sy}}^{\textit{drag}}+T^{\boldsymbol{\nabla}\kern-1pt p}_{\textit{Sy}}+T_{\textit{Sy}}^{\textit{AM}}+T_{\textit{Sy}}^{\textit{lift}}}{T_{\textit{By}}} \nonumber \\ &= 4\varGamma \left [ 9C'\delta ^2H g\left (\omega t,\frac {H}{\delta }\right )\,L_D + h(\omega t,\zeta )\,L + \frac {9.66}{\pi \sqrt {2}}\,C^{\prime \prime }\dfrac {\delta ^3}{{\textit{Re}}_\delta }\sqrt {{\textit{Re}}_s^3}L_L \right ]\!, \end{align}
where the function
$h(\omega t,\zeta )$
groups the contributions of the pressure gradient and added mass:
\begin{align} h(\omega t,\zeta ) = \frac {3}{2}\cos (\omega t)-\frac {1}{2}e^{-\zeta }\cos \left (\omega t-\zeta \right )\!. \end{align}
The onset of particle motion is predicted when
$\varUpsilon$
reaches or exceeds unity:
\begin{align} \varUpsilon \left (\omega t;\,{\textit{Re}}_\delta ,\,\delta ,\,\varGamma ,\,H,\,L\right ) \geqslant 1. \end{align}
A.3. Substrate effects at moderate Reynolds numbers
The model described so far relies on the assumption that the ambient flow velocity is well approximated by the analytical solution for Stokes’ oscillatory boundary layer, which is valid under laminar, smooth-wall flow conditions. This assumption holds when the viscous boundary layer thickness
$\delta$
is sufficiently large compared with the substrate roughness.
However, for
$\delta \lt 1$
and moderate values of
${\textit{Re}}_\delta$
, the flow can become increasingly influenced by the bottom roughness. In this regime, the effective length scale of the flow grows beyond
$\delta$
and can approach the particle diameter
$D$
as
${\textit{Re}}_\delta$
increases (Jensen, Sumer & Fredsøe Reference Jensen, Sumer and Fredsøe1989). Consequently, the actual velocity profile deviates from the Stokes solution, as shown in figure 21(a).
Nonetheless, as long as the flow does not become turbulent, a reasonable approximation of the velocity profile can be recovered by introducing an effective boundary layer thickness, defined as
$\delta _S=\beta \sqrt {2\nu /\omega }$
, where
$\beta \gt 1$
is an empirical correction factor. For example, in run
$4$
, characterised by
${\textit{Re}}_D\gt 1500$
, setting
$\beta = 1.22$
yields a much improved agreement between the model predictions and DNS results (see figure 21
b). This empirical correction extends the applicability of the torque-based motion criterion discussed in Appendix A.2 to moderate
${\textit{Re}}_D$
values, while retaining the same theoretical framework.
Velocity profiles from run
$4$
(black curves) compared with analytical Stokes profiles (red curves) shifted upward by
$z_0$
. The Stokes profiles are scaled using either (a) the viscous length scale
$\delta$
or (b) the amplified length scale
$\beta \delta$
, where the amplification coefficient is set to
$\beta =1.22$
.
Finally, the present approach can be generalised to different substrate configurations. For instance, in a compact square arrangement, the lever arm of the hydrodynamic torque is given by
$L=(\sqrt {2}/4)D$
, while the lever arm of the submerged weight
$L_{{w}}$
varies between
$D/4$
and
$(\sqrt {2}/4)D$
, depending on the flow orientation relative to the substrate geometry. In a closely packed hexagonal arrangement, these values become
$L = (\sqrt {6}/6)D$
and
$L_{{w}} = [1,\,2]\,(\sqrt {3}/12)D$
. Further investigations will be required to incorporate the effects of substrate geometry and packing density (i.e. inter-particle spacing), which introduce additional degrees of freedom to the problem.
Appendix B. Contact model
A soft-sphere collision model is employed to describe particle contact dynamics. When two spherical particles, with diameter
$D$
and centre positions at
$\boldsymbol{X}^{(i)}$
and
$\boldsymbol{X}^{(j)}$
, approach each other, contact occurs when the distance between their centres,
$\vert \boldsymbol{X}^{(i)}-\boldsymbol{X}^{(j)}\vert$
, becomes smaller than
$D+\varDelta$
. Here,
$\varDelta$
denotes a safety distance, which is taken to be equal to the grid spacing
$\Delta x$
used in the IBM, as illustrated in figure 22.
The force developing at the contact point,
$\boldsymbol{F}^{(\textit{ij})}$
, is decomposed into surface normal and tangential components,
$\boldsymbol{F}_n^{(\textit{ij})}$
and
$\boldsymbol{F}_t^{(\textit{ij})}$
, respectively. Each of these components is modelled as a spring-dashpot system, where the spring term accounts for the elastic contributions, while the dashpot term represents the dissipative, inelastic contributions. Thus, the normal and tangential components of the force that the
$i$
th particle exerts on the
$j$
th particle are given by the sum of their respective elastic and dissipative contributions:
\begin{align} \boldsymbol{F}_n^{(\textit{ij})} &= \boldsymbol{F}_n^{{el}} + \boldsymbol{F}_n^{\textit{dis}}, \end{align}
\begin{align} \boldsymbol{F}_t^{(\textit{ij})} &= \boldsymbol{F}_t^{{el}} + \boldsymbol{F}_t^{\textit{dis}}, \\[9pt] \nonumber \end{align}
respectively. The dissipative contributions
$\boldsymbol{F}_n^{\textit{dis}}$
and
$\boldsymbol{F}_t^{\textit{dis}}$
are modelled as the product of the respective components of the relative velocity
$\boldsymbol{V}^{(\textit{ij})}$
at the contact point, the effective mass (which is
$( {1}/{2})m$
for particles of equal mass
$m$
), and the respective damping coefficients,
$\mu _{n}$
and
$\mu _{t}=0.5\,\mu _{n}$
. The relative velocity
$\boldsymbol{V}^{(\textit{ij})}$
is equal to
\begin{align} \boldsymbol{V}^{(\textit{ij})} = \boldsymbol{V}^{(i)}-\boldsymbol{V}^{(j)} + L^{(\textit{ij})}(\boldsymbol{\varOmega }^{(i)} + \boldsymbol{\varOmega }^{(j)}) \times \boldsymbol{n}^{(\textit{ij})} , \end{align}
where
$\boldsymbol{V}^{(i)}$
and
$\boldsymbol{V}^{(j)}$
are the linear velocities of the particle centres,
$L^{(\textit{ij})}$
is the distance from the contact point to the centre of the
$i$
th particle,
$\boldsymbol{\varOmega }^{(i)}$
and
$\boldsymbol{\varOmega }^{(j)}$
are the angular velocities of the particles, and
$\boldsymbol{n}^{(\textit{ij})}$
is the surface-normal unit vector pointing towards the
$j$
th particle. The damping coefficients depend on the normal restitution coefficient
$e_n$
, the effective mass and the normal stiffness coefficient
$k_{n}$
(Silbert et al. Reference Silbert, Ertaş, Grest, Halsey, Levine and Plimpton2001). In particular, the normal damping coefficient is assumed to be twice the value of the tangential damping coefficient. The parameters for the contact model are summarised in table 4.
Contact model parameters. The parameters denote the Coulomb coefficient of dynamic friction
$\mu _c$
, the dimensionless normal stiffness coefficient
$\tilde k_n=6\,k_n^*\Delta x/(\pi \,\varrho _sgD^{3})$
, the dry restitution coefficient
$e_n$
, the dimensionless penalty force range
$\varDelta /\Delta x$
and number of grid points per particle diameter
$D/\Delta x$
. The reported values are commonly used in sediment transport studies, such as those reported by Mazzuoli et al. (Reference Mazzuoli, Kidanemariam, Blondeaux, Vittori and Uhlmann2016).
Schematic of the inter-particle contact kinematics. Particles
$i$
and
$j$
are assumed to have equal radii.
The normal and tangential elastic contributions to the contact force are proportional to the respective normal and tangential displacements,
$\boldsymbol{\delta }_n^{(\textit{ij})}$
and
$\boldsymbol{\delta }_t^{(\textit{ij})}$
, at the contact point, and are given by
\begin{align} \boldsymbol{F}_n^{{el}} &= -\boldsymbol{\delta }_n^{(\textit{ij})} k_n, \end{align}
\begin{align} \boldsymbol{F}_t^{{el}} &= -\boldsymbol{\delta }_t^{(\textit{ij})} k_t , \\[9pt] \nonumber \end{align}
where
$\boldsymbol{\delta }_n^{(\textit{ij})}=\delta _n\,\boldsymbol{n}^{ij}$
, and the normal deformation is
$\delta _n=d+\Delta -2\,L^{(\textit{ij})}$
(see figure 22). The tangential stiffness coefficient
$k_t$
is assumed to be
$k_t=(2/5)k_n$
, following the approach by Silbert et al. (Reference Silbert, Ertaş, Grest, Halsey, Levine and Plimpton2001). The tangential displacement is computed by accounting for the possibility that elastic energy accumulated during the contact can be released when the corresponding tangential elastic force exceeds the Coulomb friction threshold. The particles do not slide, as long as the condition
\begin{align} \left\vert \boldsymbol{F}_t^{(\textit{ij})}\right\vert \leqslant \mu _c^{(\textit{ij})}\,\left\vert \boldsymbol{F}_n^{(\textit{ij})}\right\vert \end{align}
is satisfied, where
$\mu _c^{(\textit{ij})}$
denotes the Coulomb coefficient of dynamic friction and
$\boldsymbol{F}_t^{(\textit{ij})}$
is the static friction force determined by (B2), following a widely adopted approach in granular flow models (e.g. Syamlal, Rogers & OBrien Reference Syamlal, Rogers and OBrien1993).
To compute the elastic contribution to the static friction force,
$\boldsymbol{F}_t^{{el}}$
, the tangential displacement
$\boldsymbol{\delta }^{(\textit{ij})}_t$
is updated at each time step
$\Delta t$
according to the following expression:
\begin{align} \boldsymbol{\delta }^{(\textit{ij})}_t(t+\Delta t) = \hat {\boldsymbol{\delta }}^{(\textit{ij})}_t(t) + \delta _t\,\boldsymbol{t}^{(\textit{ij})}, \end{align}
where
$\delta _t=V_{t}^{(\textit{ij})}\Delta t$
and
$\boldsymbol{t}^{(\textit{ij})}=\boldsymbol{V}_{t}^{(\textit{ij})}/\vert \boldsymbol{V}_{t}^{(\textit{ij})}\vert$
denote the increment of the tangential relative displacement and the displacement unit vector at time
$t$
, respectively. The initial value of
$\delta _t$
is given by the smaller of two values:
$V_{t}^{(\textit{ij})}\Delta t$
or
$\delta _nV_{t}^{(\textit{ij})}/V_{n}^{(\textit{ij})}$
, where
$V_{n}^{(\textit{ij})}=\boldsymbol{V}^{(\textit{ij})}\boldsymbol{\cdot }\boldsymbol{n}^{(\textit{ij})}$
represents the normal component of the relative velocity. In (B7), the accumulated tangential displacement
$\boldsymbol{\delta }^{(\textit{ij})}_t(t)$
is adjusted for rotation of the contact plane during the time interval
$\Delta t$
by the expression
\begin{align} \hat {\boldsymbol{\delta }}^{(\textit{ij})}_t(t) &= \big(\boldsymbol{\delta }^{(\textit{ij})}_t(t)\boldsymbol{\cdot }\boldsymbol{\zeta } \big)\boldsymbol{\zeta } + \big(\boldsymbol{\delta }^{(\textit{ij})}_t(t)\boldsymbol{\cdot }\boldsymbol{t}_{{old}}\big)\boldsymbol{t}_{{new}} , \end{align}
where the unit vectors
$\boldsymbol{\zeta }$
,
$\boldsymbol{t}_{{old}}$
and
$\boldsymbol{t}_{{new}}$
are defined by
\begin{align} \boldsymbol{\zeta } &= \dfrac {\boldsymbol{n}^{(\textit{ij})}(t)\times \boldsymbol{n}^{(\textit{ij})}(t+\Delta t)}{\vert \boldsymbol{n}^{(\textit{ij})}(t)\times \boldsymbol{n}^{(\textit{ij})}(t+\Delta t)\vert },\quad \boldsymbol{t}_{{old}} = \dfrac {\boldsymbol{\zeta }\times \boldsymbol{n}^{(\textit{ij})}(t)}{\vert \boldsymbol{\zeta }\times \boldsymbol{n}^{(\textit{ij})}(t)\vert } ,\quad \boldsymbol{t}_{{new}} = \dfrac {\boldsymbol{\zeta }\times \boldsymbol{n}^{(\textit{ij})}(t+\Delta t)}{\vert \boldsymbol{\zeta }\times \boldsymbol{n}^{(\textit{ij})}(t+\Delta t)\vert }. \end{align}
When the condition (B6) no longer holds, the two particles slide relative to each other at the contact point, and the static friction force, expressed by (B2), is replaced by the dynamic friction force
\begin{align} \boldsymbol{F}_t^{(\textit{ij})} =\left \lbrace \begin{array}{ll} -\mu _c^{(\textit{ij})}\, \big\vert \boldsymbol{F}_{n}^{(\textit{ij})}\big\vert \,\boldsymbol{t}^{(\textit{ij})} & \text{if}\:\boldsymbol{t}^{(\textit{ij})} \neq 0 ,\\[4pt]-\mu _c^{(\textit{ij})}\, \big\vert \boldsymbol{F}_{n}^{(\textit{ij})}\big\vert \,\frac {\boldsymbol{\delta }_t}{\vert \boldsymbol{\delta }_t\vert } & \text{if}\:\boldsymbol{t}^{(\textit{ij})} = 0,\:\boldsymbol{\delta }_t\neq 0,\\ 0 & \text{otherwise}. \end{array} \right . \end{align}
Finally, the total contact force and torque acting on the
$i$
th particle are obtained by summing the contributions from all surrounding particles in contact, as given by
\begin{align} \boldsymbol{F}_c^{(i)} &= \displaystyle \sum _{j\neq i} \boldsymbol{F}^{(\textit{ij})}, \end{align}
\begin{align} \boldsymbol{T}_c^{(i)} &= \displaystyle \sum _{j\neq i} \big (L^{(\textit{ij})}\boldsymbol{n}^{(\textit{ij})}\times \boldsymbol{F}^{(\textit{ij})}\big ) , \end{align}
respectively. To solve the contact dynamics effectively, the fixed time substep
$\Delta t_{ctc}$
used to solve (2.16) and (2.17) is chosen to be
${\mathcal O}(10^2)$
smaller than the fixed time step
$\Delta t$
adopted to advance the flow field.
Open access
