1 Introduction
Incipient particle motion is important in a wide variety of industrial and natural processes. Applications in technology include positioning (Bleil, Marr & Bechinger Reference Bleil, Marr and Bechinger2006; Sawetzki et al. Reference Sawetzki, Rahmouni, Bechinger and Marr2008) and assembling of particles (Yin et al. Reference Yin, Lu, Gates and Xia2001; Nguyen & Yoon Reference Nguyen and Yoon2009), laminar flow assays (Brooks & Tozeren Reference Brooks and Tozeren1996), cleaning of surfaces (Fan et al. Reference Fan, Soltani, Ahmadi and Hart1997; Burdick, Berman & Beaudoin Reference Burdick, Berman and Beaudoin2005), pneumatic conveying (Stevanovic et al. Reference Stevanovic, Stanojevic, Jovovic, Radic, Petrovic and Karlicic2014), removal of atmospheric pollution in cleaning of ventilating air and gravitational and centrifugal sedimentation of particle suspensions forming slurries in liquids (Derksen Reference Derksen2011) among others. Incipient particle motion is also the initial stage of several environmental phenomena such as sediment transport in rivers, granular bed erosion and dune formation (Hermann Reference Hermann2007; Groh et al. Reference Groh, Wierschem, Rehberg and Aksel2008; Wierschem et al. Reference Wierschem, Groh, Rehberg, Aksel and Kruelle2008). Given the wide range of applications, the onset of particle motion has been scope of several studies during the last century, mostly at high Reynolds numbers (Shields Reference Shields1936; Chang Reference Chang1939; Paintal Reference Paintal1971; Mantz Reference Mantz1977; Yalin & Karahan Reference Yalin and Karahan1979; Kuhnle Reference Kuhnle1993; Marsh Reference Marsh2004; Vollmer & Kleinhans Reference Vollmer and Kleinhans2007; Valyrakis et al. Reference Valyrakis, Diplas, Dancey, Greer and Celik2010). The existing data concerning the threshold for incipient motion under turbulent conditions, coming from both experimental data and theoretical or numerical models, are widely scattered (Buffington & Montgomery Reference Buffington and Montgomery1997; Marsh Reference Marsh2004). Scattering of the experimental data is mainly due to the difficulty of controlling and determining flow and particle parameters in turbulent systems. The former includes the fluctuating nature of hydrodynamic forces given by short-term nonlinear pressure gradients induced by turbulence (Vollmer & Kleinhans Reference Vollmer and Kleinhans2007) or by the duration of energetic near-bed turbulent events (Valyrakis et al. Reference Valyrakis, Diplas, Dancey, Greer and Celik2010) among others. The latter refers to size, shape, orientation, compactness and roughness of sediments, generally yielding an ambiguous definition of incipient criteria (Bravo, Ortiz & Perez-Aparicio Reference Bravo, Ortiz and Perez-Aparicio2014).
In the last decade, authors started focussing on granular motion under laminar conditions (Charru, Mouilleron & Eiffel Reference Charru, Mouilleron and Eiffel2004; Loiseleux et al. Reference Loiseleux, Gondret, Rabaud and Doppler2005; Charru et al. Reference Charru, Larrieu, Dupont and Zenit2007; Ouriemi et al. Reference Ouriemi, Aussillous, Medale, Peysoon and Guazzelli2007; Lobkovsky et al. Reference Lobkovsky, Orpe, Molloy, Kudrolli and Rothman2008; Martino, Paterson & Piva Reference Martino, Paterson and Piva2009; Ouriemi, Aussillous & Guazzelli Reference Ouriemi, Aussillous and Guazzelli2009; Derksen Reference Derksen2011; Derksen & Larsen Reference Derksen and Larsen2011; Agudo & Wierschem Reference Agudo and Wierschem2012; Seizilles et al. Reference Seizilles, Devauchelle, Lajeunesse and Metivier2013; Agudo, Dasilva & Wierschem Reference Agudo, Dasilva and Wierschem2014; Seizilles et al. Reference Seizilles, Lajeunesse, Devauchelle and Bak2014; Hong, Tao & Kudrolli Reference Hong, Tao and Kudrolli2015; Agudo et al. Reference Agudo, Luzi, Han, Hwang, Lee and Wierschem2017). Considering steady laminar flow drastically reduces difficulties in characterising the shear stress acting on the sediment and avoids the wide spectrum of length scales interacting with the bed (Derksen & Larsen Reference Derksen and Larsen2011). This regime is encountered in applications such as transportation of heavy oil (Mouilleron, Charru & Eiff Reference Mouilleron, Charru and Eiff2009) or solid gas hydrate in pipes, filtration (Olayiwola & Walzel Reference Olayiwola and Walzel2009) or microfluidics (Thompson & Bau Reference Thompson and Bau2008; Amini et al. Reference Amini, Sollier, Weaver and Dicarlo2012; Li, McKinley & Ardekani Reference Li, McKinley and Ardekani2015). Even in high-Reynolds-number flows, the particle Reynolds number may be of order one or smaller if the particles are sufficiently small (Happel & Brenner Reference Happel and Brenner1983; Guazzelli & Morris Reference Guazzelli and Morris2012).
The effect of laminar flow acting on a particle in the vicinity of a wall has been extensively studied analytically in the literature. Goldman, Cox & Brenner (Reference Goldman, Cox and Brenner1967) obtained a solution of the Stokes equation for the translational and the rotational velocity of a sphere moving parallel to a flat rigid surface in a shear flow. Later, O’Neill (Reference O’Neill1968) solved the problem for the particular case in which the sphere was in contact to the surface, providing an expression for the drag force and the corresponding torque acting on the sphere. Other authors approached the problem for semi-spheres attached to the wall (Price Reference Price1985; Sugiyama & Sbragaglia Reference Sugiyama and Sbragaglia2008) or for spherical caps exposed to an axisymmetric Stokes flow (El-Kareh & Secomb Reference El-Kareh and Secomb1996). First Saffman (Reference Saffman1965), and later Leighton & Acrivos (Reference Leighton and Acrivos1985) also considered the contribution of the lift force acting on a sphere. Taking into account lift forces, Krishnan & Leighton (Reference Krishnan and Leighton1990) obtained the translational and rotating velocity of a rough sphere rolling in a smooth horizontal plane. Yet, compared to drag force lift may be neglected at small particle Reynolds numbers (Leighton & Acrivos Reference Leighton and Acrivos1985; King & Leighton Reference King and Leighton1997). Recently, Mo et al. (Reference Mo, Zhengming, Zhipeng, Yuyun and Derksen2015) experimentally and numerically determined the critical conditions for suspension of a single sphere placed on a flat surface.
Other authors experimentally studied incipient particle motion in laminar flow over irregularly arranged sediment beds. In an annular channel, Charru et al. (Reference Charru, Mouilleron and Eiffel2004) observed an increase in compactness of a monodisperse granular bed due to the local particle arrangement. This bed armouring yielded an increase of the threshold for particle motion up to a saturated state characterised by a critical Shields number of approximately 0.12 independent of inertia. Similar values for the critical Shields number were provided in other experimental (Loiseleux et al. Reference Loiseleux, Gondret, Rabaud and Doppler2005; Ouriemi et al. Reference Ouriemi, Aussillous, Medale, Peysoon and Guazzelli2007; Seizilles et al. Reference Seizilles, Lajeunesse, Devauchelle and Bak2014) and numerical (Derksen Reference Derksen2011) studies at particle Reynolds numbers lower than one. In contrast, in a cross-sectional channel, Lobkovsky et al. (Reference Lobkovsky, Orpe, Molloy, Kudrolli and Rothman2008) reported a considerably higher value of about 0.3 at saturated conditions. Using the same experimental set-up, Hong et al. (Reference Hong, Tao and Kudrolli2015) found critical Shields numbers ranging between about 0.10 and 0.37 depending on the bed armouring level and particle size. The critical Shields number for the incipient motion of a single bead on an irregularly arranged layer of fixed spheres was found to range between 0.02 and 0.04, apparently depending on details of the bed geometry (Charru et al. Reference Charru, Larrieu, Dupont and Zenit2007).
The strong effect of the local grain arrangement on the incipient particle motion was first noticed under turbulent flow conditions. Fenton & Abbott (Reference Fenton and Abbott1977) observed that incipient conditions depend largely on the exposure of the particle to the flow. To account for the impact of geometry on the incipient motion, some authors considered geometrical factors as the angle of repose (Miller & Byrne Reference Miller and Byrne1966; Wiberg & Smith Reference Wiberg and Smith1987; Dey Reference Dey1999), exposure to the flow (Fenton & Abbott Reference Fenton and Abbott1977; Kirchner et al. Reference Kirchner, Diertrich, Iseya and Ikeda1990; Armanini & Gregoretti Reference Armanini and Gregoretti2005), relative grain protrusion (Chin & Chiew Reference Chin and Chiew1993) or streamwise bed slope (Whitehouse & Hardisty Reference Whitehouse and Hardisty1988; Dey & Debnath Reference Dey and Debnath2000; Loiseleux et al. Reference Loiseleux, Gondret, Rabaud and Doppler2005). Turning to laminar conditions, Martino et al. (Reference Martino, Paterson and Piva2009) experimentally recorded the influence of the burial degree on the onset of motion of an isolated cylinder partially exposed to the flow. Simulations of disordered assemblies of uniformly sized spheres attached to a flat bottom documented again an important impact of the local bed geometry on the hydrodynamic forces acting on the spheres (Derksen & Larsen Reference Derksen and Larsen2011). In recent experiments carried out with substrates made from regularly arranged sphere monolayers, we corroborated the strong effect of the substrate geometry on the incipient motion of a single bead (Agudo & Wierschem Reference Agudo and Wierschem2012). In laminar shear flow, we found critical Shields numbers ranging from 0.015 to 0.06 as sole function of the substrate geometry. Further experiments also revealed the role of the substrate geometry on the incipient particle motion for multi-particle systems (Agudo et al. Reference Agudo, Dasilva and Wierschem2014).
Several models including geometrical properties of the bed have been developed to predict the incipient motion, mostly under turbulent conditions. Most of them, whether analytical (Wiberg & Smith Reference Wiberg and Smith1987; Kirchner et al. Reference Kirchner, Diertrich, Iseya and Ikeda1990; Chin & Chiew Reference Chin and Chiew1993; Ling Reference Ling1995; Dey Reference Dey1999; Vollmer & Kleinhans Reference Vollmer and Kleinhans2007; Ali & Dey Reference Ali and Dey2016), numerical (Drake & Calantoni Reference Drake and Calantoni2001; McEwan & Heald Reference McEwan and Heald2001; Yergey, Beninati & Marshall Reference Yergey, Beninati and Marshall2010; Bravo et al. Reference Bravo, Ortiz and Perez-Aparicio2014), probabilistic (Cheng & Chiev Reference Cheng and Chiev1998; Papanicolaou et al. Reference Papanicolaou, Diplas, Evaggelopoulos and Fotopoulos2002; Wu & Chou Reference Wu and Chou2003) or semi-probabilistic (Soepyan et al. Reference Soepyan, Cremaschi, McLaury, Sarica, Subramani, Kouba and Gao2016), determine the threshold for incipient motion by solving the force balance on the particle or by addressing the torque balance at a contact point between particle and bed. In both cases, the hydrodynamic forces are usually modelled assuming coefficients for drag and lift forces in a certain range of Reynolds numbers. This approach has also been used for laminar flows (Loiseleux et al. Reference Loiseleux, Gondret, Rabaud and Doppler2005). Hence, a main challenge for the corresponding models is a proper formulation for the drag and lift coefficients, which can be applied in the proper range of Reynolds numbers, bed geometries, shapes or particle orientations (Mando & Rosendahl Reference Mando and Rosendahl2010). The definition of drag and lift coefficients, however, varies widely in literature. In addition, most of the theoretical and empirical approaches for the drag coefficient at defined Reynolds numbers are based on a freely falling particle (Morsi & Alexander Reference Morsi and Alexander1972; Wiberg & Smith Reference Wiberg and Smith1987; Ling Reference Ling1995; Dey Reference Dey1999; Bravo et al. Reference Bravo, Ortiz and Perez-Aparicio2014) or are only valid for a specific bed geometry (Coleman Reference Coleman and Collins1967; Vollmer & Kleinhans Reference Vollmer and Kleinhans2007; Ali & Dey Reference Ali and Dey2016). Regarding the lift coefficient, there is a much smaller theoretical or empirical basis in literature than for drag. A typical assumption is that the lift coefficient was proportional to the drag. Yet, the proportionality factor strongly depends on the considered Reynolds number range (Saffman Reference Saffman1965; James Reference James1990; Dey Reference Dey1999) and on the exposure degree (Brayshaw, Lynne & Reid Reference Brayshaw, Lynne and Reid1983; Vollmer & Kleinhans Reference Vollmer and Kleinhans2007) and, hence, on the bed geometry.
A further problem in computing the hydrodynamic forces arises from the velocity profile over the exposed area of the particle. At high Reynolds numbers, the main tendency of current models is to assume a mean logarithmic velocity distribution in the wall layer corrected by the equivalent roughness height of Nikuradse (Nikuradse Reference Nikuradse1933; van Rijn Reference van Rijn1984; Wiberg & Smith Reference Wiberg and Smith1987; Nezu & Nakagawa Reference Nezu and Nakagawa1994; Wu & Chou Reference Wu and Chou2003; Vollmer & Kleinhans Reference Vollmer and Kleinhans2007). For laminar conditions, the velocity profile close to the particle is typically assumed to vary linearly with the height (Ouriemi et al. Reference Ouriemi, Aussillous and Guazzelli2009; Derksen & Larsen Reference Derksen and Larsen2011). This linear assumption can also be applied to the viscous sublayer in turbulent flows (Ling Reference Ling1995; Dey Reference Dey1999; Bravo et al. Reference Bravo, Ortiz and Perez-Aparicio2014; Ali & Dey Reference Ali and Dey2016) at roughness Reynolds numbers below 3 (Wiberg & Smith Reference Wiberg and Smith1987). Note that the definition of the roughness Reynolds number is based on the scale length of the bed roughness and therefore is based on an average bed geometrical factor.
For turbulent conditions, different studies evidence the existence of a non-zero average velocity on the top of the sediment bed (James & Davis Reference James and Davis2001; Pokrajac, Manes & McEwan Reference Pokrajac, Manes and McEwan2007; Yergey et al. Reference Yergey, Beninati and Marshall2010). At high Reynolds numbers, the main tendency of current models is to adopt a virtual or effective zero level at a constant distance below the top of the substrate. Averaging over a region upstream and downstream of the considered particle, Kirchner et al. (Reference Kirchner, Diertrich, Iseya and Ikeda1990), for instance, defined the effective zero level as 16 % of the particle diameter below the top of the bed spheres. This definition, however, was assumed to characterise the small-scale topographic structure of a variety of mixed sphere surfaces. Other theoretical models consider the effective zero level as 25 % of the particle diameter below the top of the bed spheres (Cheng & Chiev Reference Cheng and Chiev1998; Dey Reference Dey1999; Wu & Chou Reference Wu and Chou2003; Vollmer & Kleinhans Reference Vollmer and Kleinhans2007) based on previous studies of van Rijn (Reference van Rijn1984). Later experimental observations on a mobile bed revealed a flow penetration of 21 % below the top of the bed particles for a weakly mobile bed (Dey et al. Reference Dey, Das, Gaudio and Bose2012). In laminar flow, Goharzadeh, Khalili & Jorgensen (Reference Goharzadeh, Khalili and Jorgensen2005) showed experimentally that the flow penetration is of the order of the particle diameter and depends on the particle Reynolds number. A flow penetration dependency on the Shields number was later noticed experimentally (Mouilleron et al. Reference Mouilleron, Charru and Eiff2009) and numerically (Ouriemi et al. Reference Ouriemi, Aussillous and Guazzelli2009). In a rigorous model, the level of flow penetration should be considered as a function of the bed geometry as is also evident from the numerical studies by Derksen & Larsen (Reference Derksen and Larsen2011). This premise has not been contemplated in current models so far. The level of flow penetration strongly affects the average velocity with which the hydrodynamic forces are computed and also the point of attack of the drag force (Wiberg & Smith Reference Wiberg and Smith1987; Wu & Chou Reference Wu and Chou2003; Vollmer & Kleinhans Reference Vollmer and Kleinhans2007; Ali & Dey Reference Ali and Dey2016). Consequently, the critical threshold for incipient motion is very sensitive to the flow penetration level as shown by Yergey et al. (Reference Yergey, Beninati and Marshall2010).
Another source of uncertainty in current models is related to the considered mechanism of incipient motion. At high Reynolds numbers, the mechanism of incipient motion depends on the exposure of the particle to the flow and thus on the local bed geometry. For well-exposed spheres at near-critical conditions, motion is supposed to occur by rolling (Wu & Chou Reference Wu and Chou2003; Valyrakis et al. Reference Valyrakis, Diplas, Dancey, Greer and Celik2010). For individual particles buried within the sediment bed, lifting seems to be the most appropriate mechanism (Valyrakis et al. Reference Valyrakis, Diplas, Dancey, Greer and Celik2010). Hence, most models consider rolling (White Reference White1940; James Reference James1990; Ling Reference Ling1995) or lift (Ling Reference Ling1995) as the incipient mechanism. Another mechanism that has been studied is sliding (Paintal Reference Paintal1971; Wiberg & Smith Reference Wiberg and Smith1987; Loiseleux et al. Reference Loiseleux, Gondret, Rabaud and Doppler2005) including also spherical particles rotating in place before dislodgment from their initial positions (Dey Reference Dey1999). For laminar flows, the situation apparently simplifies. Rolling and sliding are assumed to be the main incipient mechanisms (Loiseleux et al. Reference Loiseleux, Gondret, Rabaud and Doppler2005) while lift forces are usually neglected (Leighton & Acrivos Reference Leighton and Acrivos1985; Wiberg & Smith Reference Wiberg and Smith1987; Ling Reference Ling1995; Derksen Reference Derksen2011; Hong et al. Reference Hong, Tao and Kudrolli2015). For individual beads, Agudo et al. (Reference Agudo, Dasilva and Wierschem2014) showed experimentally that rolling was the only mechanism of incipient motion even for buried spheres. Recently, Kudrolli, Scheff & Allen (Reference Kudrolli, Scheff and Allen2016) have shown experimentally and theoretically that rolling stands as the mechanism of motion for well-exposed individual beads at particle Reynolds numbers up to 1500. This points to a preferential torque balance rather than a force balance to model the onset of an individual bead exposed to the flow. In the case of more than one loose bead, sliding may appear as mechanism of incipient motion depending on the bed geometry as experimentally observed by Agudo et al. (Reference Agudo, Dasilva and Wierschem2014).
Another limitation of current models arises from the evaluation of the lever arms, which strongly depend on the streamwise orientation of the particle and thus on the local bed geometry. Most models assume average or empirical geometrical factors to determine the lever arms or, more generally, the effect of the local geometry on the hydrodynamic force (Wiberg & Smith Reference Wiberg and Smith1987; Kirchner et al. Reference Kirchner, Diertrich, Iseya and Ikeda1990; Chin & Chiew Reference Chin and Chiew1993; Ling Reference Ling1995; Dey Reference Dey1999; Vollmer & Kleinhans Reference Vollmer and Kleinhans2007). It is only recently that the lever arms of hydrodynamic forces have been analysed in detail based on the three-dimensional bed geometry (Ali & Dey Reference Ali and Dey2016). The model presented by Ali & Dey (Reference Ali and Dey2016), however, considers a bead on a closely packed sediment bed only. For laminar conditions, it is clear that the resultant hydrodynamic force does not act on the particles centre of gravity since viscous forces are expected to be stronger on the surface area exposed to the flow. This has been shown experimentally by Kudrolli et al. (Reference Kudrolli, Scheff and Allen2016) in their study of an individual sphere resting on a rough surface in a cone-plate rheometer. This also holds true for high Reynolds numbers, since the line of action of the hydrodynamic force depends on the velocity distribution close to the exposed surface area of the particle (Dey Reference Dey1999; Wu & Chou Reference Wu and Chou2003; Vollmer & Kleinhans Reference Vollmer and Kleinhans2007; Kudrolli et al. Reference Kudrolli, Scheff and Allen2016).
Figure 1. Single bead resting on a quadratic substrate exposed to a flow with shear rate
$\dot{\unicode[STIX]{x1D6FE}}$
and side view with angle between the symmetry axis of the substrate and the flow direction,
$\unicode[STIX]{x1D6FD}$
, equal zero.
Given the aforementioned limitations for predicting the incipient motion as a function of the local bed geometry, the actual state of the art still provides a scattered range of thresholds depending on the premises or assumptions adopted in the models. This is shown for instance by Marsh (Reference Marsh2004) in their comparison of four different models predicting incipient motion on irregularly arranged sediment beds at high Reynolds numbers. In laminar flow, the number of models predicting incipient motion as a function of the bed geometry is scant compared to turbulent flow. This stands despite its well-documented importance on the critical Shields number.
In the present article, we address incipient particle motion of a single bead as a function of geometrical properties of the substrate at small particle Reynolds numbers. We study the steady shear flow over quadratically arranged sphere monolayers numerically and determine the effective zero level of the shear flows. To describe the impact of the substrate geometry on the critical Shields number under creeping flow conditions, we propose an analytical model. The model does not rely on any empirical value in order to evaluate the hydrodynamic forces. Hence, no calibration is needed. Nor are any empirical or average geometrical factors assumed to address the effect of the local geometry on the critical Shields number. The model is in good accordance with recent experimental results (Agudo & Wierschem Reference Agudo and Wierschem2012; Agudo et al. Reference Agudo, Dasilva and Wierschem2014) and numerics, showing a better agreement than other models in the literature.
Figure 2. Sketch of the forces acting on a single bead resting on the substrate: side view (a) and top view (b). The surface of the spherical cap exposed to the shear flow,
$S_{CAP}$
, is darkened. The red arrow in (b) represents the main flow direction and
$\unicode[STIX]{x1D6FD}$
the orientation angle of the main flow direction with respect to the symmetry axis of the substrate.
The article is organised as follows. The considered system is defined in § 2. The numerical study is described in § 3 and the analytical model is derived in § 4. The results are discussed and compared to existing models in § 5. Finally, the conclusions are summarised in § 6.
2 Formulation of the problem
We consider a single sphere of diameter
$D_{P}$
deposited on a regular substrate and exposed to a steady shear flow of a viscous incompressible Newtonian fluid of constant density,
$\unicode[STIX]{x1D70C}$
, and kinematic viscosity,
$\unicode[STIX]{x1D708}$
. The substrate consists of a layer of fixed spheres of same diameter
$D_{P}$
, as in previous experiments (Agudo & Wierschem Reference Agudo and Wierschem2012; Agudo et al.
Reference Agudo, Dasilva and Wierschem2014). Figure 1 depicts an example for a single bead resting on a quadratically arranged substrate with constant spacing,
$a$
, between beads. On the downstream side, the bead contacts the substrate spheres at point
$A$
. The equilibrium position of the single bead resting on the substrate is characterised by the angle of repose,
$\unicode[STIX]{x1D711}$
. The angle of repose, which coincides with the pivot angle for the case of non-cohesive beads, is the angle of the slope that has to be overcome by the incipient motion (Julien Reference Julien2010).
$\unicode[STIX]{x1D712}$
is the angle between the bead surface at effective zero level and the vertical axis. The orientation angle of the main flow direction with respect to the symmetry axis of the substrate is defined as
$\unicode[STIX]{x1D6FD}$
, see figure 2(b). As coordinate system we use spherical coordinates with its origin in the centre of the exposed bead as shown in figure 1. Its radial, polar and azimuthal coordinates are indicated by
$r$
,
$\unicode[STIX]{x1D717}$
and
$\unicode[STIX]{x1D719}$
, respectively.
We assume slight penetration of the shear flow into the substrate as shown for laminar flows in the previous numerical study by Derksen & Larsen (Reference Derksen and Larsen2011). They found for an arbitrarily arranged monolayer of spheres in a shear flow an effective zero level for the average flow profile that depends on the sphere occupancy, i.e. the number of spheres per area. In a Cartesian coordinate system with its origin in the centre of the mobile sphere, as shown in figure 1, the profile of the undisturbed flow is
$\dot{\unicode[STIX]{x1D6FE}}(z^{\ast }+z)\hat{e_{x}}$
, where
$\dot{\unicode[STIX]{x1D6FE}}$
is the shear rate,
$z^{\ast }$
is the vertical distance between the centre of the sphere and the effective zero level and
$\hat{e_{x}}$
is the unit vector in flow direction. Accordingly, the particle Reynolds number for the shear flow is given by
$$\begin{eqnarray}Re_{P}=\frac{\dot{\unicode[STIX]{x1D6FE}}D_{P}^{2}}{\unicode[STIX]{x1D708}}.\end{eqnarray}$$
In the following, we consider particle Reynolds numbers below one. Figure 2 illustrates the forces acting on the bead. In general, the shear flow exerts a drag force,
$F_{D}$
, and a lift force,
$F_{L}$
, on the particle. The effective weight,
$F_{G}-F_{B}$
, where
$F_{G}$
and
$F_{B}$
are gravity and buoyancy force, respectively, pushes the particle on the substrate, and interacts via the Coulomb friction force,
$F_{T}$
, and normal force,
$F_{N}$
.
$L_{D}$
,
$L_{L}$
and
$L_{G}$
, are the level arms of drag force, lift force and effective weight for the rotation axis in
$A$
;
$l_{D}$
and
$l_{L}$
are the level arms of drag force and lift force for the rotation axis in the bead centre,
$O$
, respectively. The surface of the spherical cap exposed to the shear flow is darkened in figure 2(a). The relative magnitude of the characteristic flow-induced shear force to the effective particle weight is described by the Shields number,
$\unicode[STIX]{x1D703}$
:
$$\begin{eqnarray}\unicode[STIX]{x1D703}=\frac{\unicode[STIX]{x1D708}\dot{\unicode[STIX]{x1D6FE}}}{(\unicode[STIX]{x1D70C}_{S}/\unicode[STIX]{x1D70C}-1)gD_{P}},\end{eqnarray}$$
where
$g$
is the acceleration due to gravity and
$\unicode[STIX]{x1D70C}_{S}$
is the density of the mobile bead. At
$Re_{P}$
lower than one, the critical Shields number for incipient motion,
$\unicode[STIX]{x1D703}_{C}$
, was found to be independent of relative density and viscosity, and therefore of inertia (Agudo & Wierschem Reference Agudo and Wierschem2012; Agudo et al.
Reference Agudo, Dasilva and Wierschem2014).
3 Numerical approach
Numerically we studied systems with substrates made from monolayers of quadratically arranged spheres. The monolayer is formed by 14 particles along and 18 particles perpendicular to the main flow direction. The orientation angle,
$\unicode[STIX]{x1D6FD}$
, is equal to zero. We checked that the critical Shields number is independent from side, inflow and outflow effects at these dimensions. The single bead is located on top of the substrate centre. To employ two-dimensional inflow and outflow conditions, we implement a step of height
$D_{P}$
before and after the substrate. Each step has a length of
$8.62\cdot D_{P}$
. Spacing between the substrate spheres varies from approximately zero to slightly below
$0.414\cdot D_{P}$
, which is the maximum value for placing a single sphere on top of the substrate. The shear flow is produced by a top plate moving horizontally at constant speed. The gap between the top of the substrate and the moving plate,
$h$
, is
$5\cdot D_{P}$
in accordance with experiments by Agudo & Wierschem (Reference Agudo and Wierschem2012). The geometry is shown in figure 3.
Figure 3. Geometry for the numerical calculation of the flow field around a single sphere and the determination of the critical Shields number. Flow direction is from left to right.
The steady Navier–Stokes equations were solved numerically with the finite volume method using the SIMPLE algorithm implemented in OpenFOAM-2.3.1. To solve the Poisson equation for the pressure, we used the generalised geometric algebraic multi-grid solver. The domain was modelled and discretised into a structured mesh of hexahedral elements with SALOME 7.5.1. A maximum number of 3 refinement levels were defined at the single sphere within the substrate, as shown in figure 4. Within a refinement level, the starting element edge length is halved around all the spheres. At all solid boundaries, the no-slip condition holds for the velocity and zero gradient for the pressure. The flanks are modelled by periodic boundary conditions for velocity and pressure. At the inlet, we defined a Couette flow and used the zero-gradient condition for the pressure. At the outlet, we used no backpressure and the zero-gradient condition for the velocity.
Figure 4. Structured mesh of hexahedral elements with an element edge length of
$0.0788\cdot D_{P}$
and approximately 7 million elements. The zoom depicts the meshing area of the single particle. The element edge length at the particle surface is refined to an element edge length of
$0.0197\cdot D_{P}$
.
To focus on creeping flow conditions, the particle Reynolds number was fixed at 0.01 in all numerical calculations. The particle on the substrate remained at rest. For each geometry, the critical Shields number, defined as in (2.2), is determined from the particle density for which the torque balance or the force balance is zero, respectively. The torque balance corresponds to the critical Shields number for the onset of rolling motion while the force balance to that for sliding motion. For the onset of rolling motion, the particle density reads
$$\begin{eqnarray}\unicode[STIX]{x1D70C}_{S}=\frac{F_{B}L_{G}+M_{\unicode[STIX]{x1D70F}}+M_{P}}{gV_{P}L_{G}},\end{eqnarray}$$
where
$V_{P}$
, is the particle volume.
$M_{\unicode[STIX]{x1D70F}}$
and
$M_{P}$
are the torque caused by friction and pressure exerted by the flow, respectively. Forces and torques are calculated using the OpenFOAM forces library. For the onset of sliding motion, we consider
$$\begin{eqnarray}|F_{T}|=\unicode[STIX]{x1D707}F_{N},\end{eqnarray}$$
where
$\unicode[STIX]{x1D707}$
is the friction coefficient. In this case, the particle density reads
$$\begin{eqnarray}\unicode[STIX]{x1D70C}_{S}=\frac{1}{gV_{P}}\left(F_{D}\frac{\cos \unicode[STIX]{x1D711}-\unicode[STIX]{x1D707}\sin \unicode[STIX]{x1D711}}{\sin \unicode[STIX]{x1D711}+\unicode[STIX]{x1D707}\cos \unicode[STIX]{x1D711}}+F_{L}+F_{B}\right).\end{eqnarray}$$
For sliding motion, the highest density and hence the smallest critical Shields number is obtained in the frictionless case. At the low particle Reynolds number studied, the lift force is smaller than the drag force by 3–4 decimal powers in accordance with Derksen & Larsen (Reference Derksen and Larsen2011). Due to their small contribution, we disregard the lift force in the following considerations.
We assume numerical convergence at residual values for velocity and pressure below
$10^{-6}$
. Figure 5 shows a grid convergence study for the critical Shields number for rolling motion on the substrate with a spacing of
$0.0345\cdot D_{P}$
. The critical Shields number increases with the number of volume elements and tends to an asymptotic value. Yet, the difference between 29 million and 41 million elements is approximately 0.35 %. To check grid independency, we calculated the grid convergence index (GCI) (Roache Reference Roache1994, Reference Roache1997, Reference Roache1998). The
$\text{GCI}_{12}$
of the critical Shields number for the finest grid has a value of 0.0048 %. The
$\text{GCI}_{23}$
between the grids with 29 and 19.5 million elements is 0.01 %. Both GCIs indicate that our simulations are within the asymptotic range of convergence. Based on this result we used 29 million elements for all simulations due to computational efficiency.
Figure 5. Grid convergence study for the critical Shields number for rolling motion. Substrate spacing
$a/D_{P}=0.035$
. The calculated
$\text{GCI}_{12}$
for the finest mesh refinement is
$4.8\times 10^{-5}$
. The line is to guide the eye.
Figure 6. Impact of shielding by the substrate as a function of the spacing between substrate spheres: (a) relative drag force, (b) ratio between pressure and friction forces in main flow direction and (c) distance of the lines of action of these forces above the bead centre. The shielding factor
$b$
in (a) is the ratio of the drag force to the solution for a sphere on a plane wall (Goldman et al.
Reference Goldman, Cox and Brenner1967). Solid and open symbols in (c) indicate the distance of the lines of action of friction and pressure forces from the bead centre, respectively.
To check the accuracy of our numerical implementation, we also calculated the drag force on a single sphere near a plain wall in creeping flow. Goldman et al. (Reference Goldman, Cox and Brenner1967) obtained for the factor by which the drag force deviates from that in unbounded flow a value of 1.7005, while Chaoui & Feuillebois (Reference Chaoui and Feuillebois2003) calculated the factor to be 1.7006 and O’Neill (Reference O’Neill1968) to be 1.7009. For about 29 million elements, we obtained a factor of 1.7104, which is approximately 0.6 % higher than the other values. Like O’Neill (Reference O’Neill1968) and Derksen & Larsen (Reference Derksen and Larsen2011) we also calculated the scaled drag and lift forces
$F_{D}^{\ast }$
and
$F_{L}^{\ast }$
:
$$\begin{eqnarray}\displaystyle & \displaystyle F_{D}^{\ast }=\frac{F_{D}}{\unicode[STIX]{x1D702}\dot{\unicode[STIX]{x1D6FE}}(D_{P}/2)^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle F_{L}^{\ast }=\frac{4}{Re}\frac{F_{L}}{\unicode[STIX]{x1D702}\dot{\unicode[STIX]{x1D6FE}}(D_{P}/2)^{2}}, & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D702}$
is the dynamic viscosity. We obtained for
$F_{D}^{\ast }$
a value of 32.3 and for
$F_{D}^{\ast }$
, 9.2, which corresponds to deviations of 0.3 % and 0.2 % from the respective values proposed by Derksen & Larsen (Reference Derksen and Larsen2011). We further checked the accuracy by comparing it to the study of shear flow past a bump on a plane wall by Price (Reference Price1985). While Price studied the creeping flow limit, we considered a flow at a particle Reynolds number of 0.01. Retaining from an infinite series of terms a finite number, Price obtained for the dimensionless drag force defined by (3.4) a value of 0.716 and a separation point at
$73.3^{\circ }$
. For the drag force, we obtained numerically a value of 0.695, which is 3 % below Prices result. We found the separation point at an angle of
$76.4^{\circ }$
, again quite close to Price’s results.
The monolayer partly shields the sphere against the shear flow and, hence, weakens the drag force exerted on the sphere. Figure 6(a) depicts the relative impact of shielding,
$b$
, by comparing the drag force to the solution of Goldman et al. (Reference Goldman, Cox and Brenner1967). Accordingly,
$b=F_{D}/(1.5\unicode[STIX]{x03C0}f\unicode[STIX]{x1D702}\dot{\unicode[STIX]{x1D6FE}}D_{P}^{2})$
where
$f$
is equal to 1.7005. The figure shows that shielding reduces the drag force by approximately 30 % if the substrate spheres are close to contact. With increasing spacing, the sphere is buried deeper into the substrate and the drag force is reduced to less than 20 % at the largest spacing studied. As shows figure 6(b), the ratio between contributions due to pressure,
$F_{Px}$
, and friction,
$F_{\unicode[STIX]{x1D70F}x}$
, remain close to 1:2, which is the ratio for a bead in a plug flow at creeping flow conditions (Durst Reference Durst2008). Hence, the effect of shielding is roughly the same for both contributions. Figure 6(c) depicts the distance of the lines of action of friction and pressure forces in main flow direction above the bead centre,
$l_{D\unicode[STIX]{x1D70F}}/D_{P}$
and
$l_{D\unicode[STIX]{x1D70F}}/D_{P}$
, respectively. While the line of action of the friction force moves away from the centre as the sphere is buried deeper into the substrate at higher spacing, the line of action of the pressure force remains at the sphere centre, as expected for a normal force distribution on a sphere.
The critical Shields number for rolling motion and frictionless sliding as a function of the spacing between substrate spheres is shown in figure 7. At narrow spacing, the Shields number hardly changes. As spacing widens up to its maximum, the critical Shields number ever-stronger increases. As shows figure 7, the critical Shields number for rolling motion is always smaller than for sliding. At the narrowest spacing, the critical Shields number is approximately 0.039 and it is approximately 0.44 at the largest spacing.
Figure 7. Critical Shields number as a function of the spacing between substrate spheres. Solid and open symbols indicate the critical Shields numbers for frictionless sliding and for rolling motion, respectively.
Figure 8(a) shows the velocity in main flow direction,
$U_{X}$
, averaged over four unit cells of substrate spheres at the centre of the configuration. The velocity remains about zero until a certain level below the top of the substrate. From here upwards, the velocity increases and approaches the ideal Couette profile at approximately 20 % of the particle diameter above the substrate. Hence, the mean velocity gradient at the substrate centre is slightly smaller than that at the steps. Extrapolating the linear profile over the substrate centre to zero yields an effective zero level,
$z_{0}$
, below the top of the substrate. Figure 8(b) shows that the effective zero level, obtained from a linear fit to the velocity data in the range between
$0.4\cdot D_{P}$
and
$4.4\cdot D_{P}$
above the top of the substrate spheres, increases linearly from approximately 8 % for narrow spacing to approximately 13 % at the largest spacing:
$z_{0}/D_{P}=0.077+0.13(a/D_{P})$
.
Figure 8. (a) Velocity in main flow direction averaged over four unit cells of substrate spheres at the centre of the configuration together with a linear fit indicated by solid (black) and dashed (blue) lines, respectively. The inset depicts a zoom at the top of the substrate. Extrapolation of the linear fit to zero velocity yields the effective zero level,
$z_{0}$
. Substrate spacing
$a/D_{P}:0.268$
. Effective zero level below the top of the substrate. The line is a linear fit to the data. (b) Effective zero level below the top of the substrate. The line is a linear fit to the data.
The velocity in main flow direction within the interstices of the monolayer is shown in figure 9 comparing the flow with and without a bead on top of the substrate. The presence of the bead alters the flow through the interstices noticeable in a neighbourhood of approximately 3 bead diameters. Particularly it deflects part of the flow through the void underneath the bead. Figure 9(c) depicts the volume flux underneath the bead centre plain,
$\dot{V}_{U}$
. It has been non-dimensionalised with the volume flux
$\dot{V}_{S}=\unicode[STIX]{x03C0}(D_{P}/2)^{2}\dot{\unicode[STIX]{x1D6FE}}D_{P}/2$
. The volume flux underneath the bead has a local maximum at a spacing of less than 30 % of the sphere diameter. It decreases at narrower spacing due to the proximity of the substrate spheres and at larger spacing due to the proximity of the bead to the bottom plate.
Figure 9. (a) Velocity in main flow direction and streamlines within the substrate interstices with (a) and without (b) a particle at the substrate centre. The opaque regions indicate the location of the substrate spheres. Main flow is from left to right. Substrate spacing
$a/D_{P}$
: 0.268. (c) Dimensionless volume flux underneath the sphere as a function of substrate spacing.
4 Derivation of the analytical model
As shown numerically, rolling is the incipient mechanism of motion for spherical particles. This is in agreement with previous experimental (Agudo et al. Reference Agudo, Dasilva and Wierschem2014, Reference Agudo, Luzi, Han, Hwang, Lee and Wierschem2017; Kudrolli et al. Reference Kudrolli, Scheff and Allen2016) and theoretical observations (Kudrolli et al. Reference Kudrolli, Scheff and Allen2016) for an individual bead. In what follows, we derive a model for incipient rolling motion based on the torque balance. We use the linear velocity profile with an effective zero level acting on the exposed bead as the driving mechanism. The exposure of the bead depends on the substrate geometry. To enable motion, the bead has to overcome the angle of repose, which again is a function of the substrate geometry. Lift forces are not considered in the balance of forces, yet they are negligible compared with drag forces as observed in the numerical calculations. The force and torque balances thus read
$$\begin{eqnarray}\displaystyle & \displaystyle F_{N}=F_{D}\sin \unicode[STIX]{x1D711}+(F_{G}-F_{B})\cos \unicode[STIX]{x1D711}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle F_{T}=F_{D}\cos \unicode[STIX]{x1D711}-(F_{G}-F_{B})\sin \unicode[STIX]{x1D711}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle F_{D}L_{D}=(F_{G}-F_{B})L_{G}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle |F_{T}|\leqslant \unicode[STIX]{x1D707}F_{N}. & \displaystyle\end{eqnarray}$$
At low particle Reynolds numbers, the bead moves through the troughs of the substrate and does not travel in main flow direction (Agudo & Wierschem Reference Agudo and Wierschem2012; Agudo et al.
Reference Agudo, Dasilva and Wierschem2014). This results in a reduction of the shear force exerted by the flow in the direction of particle motion if the troughs are not in line with the main flow direction. Hence, including the angle between the substrate symmetry axis and the main flow direction,
$\unicode[STIX]{x1D6FD}$
, (4.3) reads
$$\begin{eqnarray}F_{D}\cos \unicode[STIX]{x1D6FD}L_{D}=(F_{G}-F_{B})L_{G}.\end{eqnarray}$$
The effective weight acting on the centre of gravity of the sphere is given by
$$\begin{eqnarray}(F_{G}-F_{B})=(\unicode[STIX]{x1D70C}_{S}-\unicode[STIX]{x1D70C})g\frac{\unicode[STIX]{x03C0}D_{P}^{3}}{6}.\end{eqnarray}$$
The lever arm of the effective weight
$L_{G}$
can be computed analytically. Expressions for the drag force and its lever arm with respect to the axis of rotation (
$A$
) will be derived subsequently.
The drag force exerted by the linear shear flow on a sphere is given by
$$\begin{eqnarray}F_{D}={\textstyle \frac{3}{2}}\unicode[STIX]{x03C0}f\unicode[STIX]{x1D702}\dot{\unicode[STIX]{x1D6FE}}D_{P}^{2},\end{eqnarray}$$
where
$f$
is a dimensionless coefficient. A value of 1 corresponds to an infinite distance of the sphere to a wall (Goldman et al.
Reference Goldman, Cox and Brenner1967). For a sphere in contact with a plane wall, values of approximately 1.7 (Goldman et al.
Reference Goldman, Cox and Brenner1967; O’Neill Reference O’Neill1968; Chaoui & Feuillebois Reference Chaoui and Feuillebois2003) have been obtained. For a solid hemisphere (Price Reference Price1985) and for a highly viscous hemispherical droplet attached at a planar wall (Sugiyama & Sbragaglia Reference Sugiyama and Sbragaglia2008), the coefficient
$f$
is smaller than 1, which is mainly due to a lower surface area exposed to the flow.
To model the drag force and the flow-induced torque acting on the bead, we start from the solution for a single sphere on a planar wall. Shielding of the mobile bead by the substrate is taken into account by a correction factor,
$b$
, that depends on the substrate geometry. The drag force thus becomes
$$\begin{eqnarray}F_{D}={\textstyle \frac{3}{2}}\unicode[STIX]{x03C0}bf\unicode[STIX]{x1D702}\dot{\unicode[STIX]{x1D6FE}}D_{P}^{2}.\end{eqnarray}$$
In order to evaluate the shielding factor
$b$
, one may consider different geometrical quantities that characterise the fraction of the bead exposed to the flow. One of them is the exposure degree,
$E$
. It is defined as the ratio of the cross-sectional area effectively exposed to the incident flow to the total cross-sectional area of the particle. It is frequently used in pressure-driven particle motion (Wiberg & Smith Reference Wiberg and Smith1987; Wu & Chou Reference Wu and Chou2003; Vollmer & Kleinhans Reference Vollmer and Kleinhans2007). Another parameter is the protrusion, which is the distance between the top of the bead and the top of the substrate spheres (Fenton & Abbott Reference Fenton and Abbott1977; Armanini & Gregoretti Reference Armanini and Gregoretti2005).
The drag force on the bead is due to friction and pressure. For a sphere in unbounded plug flow, the Stokes law shows that viscous forces are twice as large as pressure forces (Durst Reference Durst2008). Figure 6(b) shows that this also holds true approximately for the bead on the substrate. Yet, deviations from the ratio 2 : 1 indicate that the geometry does not affect pressure and friction forces equally. Therefore, we model the effect of shielding on pressure and friction forces separately. Hence, the shielding factor
$b$
is given by
$b=(2/3)b_{\unicode[STIX]{x1D70F}}+(1/3)b_{p}$
, where
$b_{\unicode[STIX]{x1D70F}}$
and
$b_{p}$
are the shielding factors for pressure and friction forces, respectively. For the pressure force, a correction addressing the cross-sectional area is appropriate (Wiberg & Smith Reference Wiberg and Smith1987; Wu & Chou Reference Wu and Chou2003; Vollmer & Kleinhans Reference Vollmer and Kleinhans2007). For the friction force, however, considering the exposed surface area is more adequate.
The surface of the spherical cap exposed to the shear flow,
$S_{CAP}$
, is (see figure 2
a)
$$\begin{eqnarray}S_{CAP}=\int _{0}^{\unicode[STIX]{x03C0}-\unicode[STIX]{x1D712}}\int _{0}^{2\unicode[STIX]{x03C0}}\left(\frac{D_{P}^{2}}{4}\right)\sin \unicode[STIX]{x1D717}\,\text{d}\unicode[STIX]{x1D719}\,\text{d}\unicode[STIX]{x1D717}=\frac{\unicode[STIX]{x03C0}D_{P}^{2}}{2}(1-\cos (\unicode[STIX]{x03C0}-\unicode[STIX]{x1D712}))=\frac{\unicode[STIX]{x03C0}D_{P}^{2}}{2}(1+\cos \unicode[STIX]{x1D712})\end{eqnarray}$$
and the corresponding correction factor for the friction force yields
$$\begin{eqnarray}b_{\unicode[STIX]{x1D70F}}=\frac{S_{CAP}}{S}=\frac{1}{2}(1+\cos \unicode[STIX]{x1D712}),\end{eqnarray}$$
where
$S$
is the bead surface. Similarly, the cross-sectional area exposed to the flow is
$$\begin{eqnarray}b_{p}=1-\frac{\unicode[STIX]{x1D712}}{\unicode[STIX]{x03C0}}+\frac{\sin (2\unicode[STIX]{x1D712})}{2\unicode[STIX]{x03C0}}.\end{eqnarray}$$
Hence, the main Shielding factor yields
$$\begin{eqnarray}b=\frac{1}{3}\left(2-\frac{\unicode[STIX]{x1D712}}{\unicode[STIX]{x03C0}}+\cos \unicode[STIX]{x1D712}+\frac{\sin (2\unicode[STIX]{x1D712})}{2\unicode[STIX]{x03C0}}\right).\end{eqnarray}$$
Substituting (4.12) into (4.8), the drag force exerted by the shear flow results in
$$\begin{eqnarray}F_{D}=\left(\unicode[STIX]{x03C0}+\frac{\unicode[STIX]{x03C0}\cos \unicode[STIX]{x1D712}}{2}+\frac{\sin (2\unicode[STIX]{x1D712})}{3}-\frac{\unicode[STIX]{x1D712}}{2}\right)f\unicode[STIX]{x1D702}\dot{\unicode[STIX]{x1D6FE}}D_{P}^{2}.\end{eqnarray}$$
Now, we consider the level arm of the drag force with respect to the axis of rotation (
$A$
) in figure 2. Since pressure is a normal stress, the pressure force acts on the bead centre, see figure 6(c). The line of action of the friction force, however, lies above the bead centre. Due to the linear velocity profile, we obtain the viscous force exerted by the flow from integrating the linear distributed profile,
$\dot{\unicode[STIX]{x1D6FE}}z$
, over the spherical cap exposed to the flow. Hence, the line of action of the friction force is at the height
$h_{eff}$
above the effective zero level (see figure 10):
$$\begin{eqnarray}h_{eff}=\frac{\displaystyle \int _{0}^{\unicode[STIX]{x03C0}-\unicode[STIX]{x1D712}}\dot{\unicode[STIX]{x1D6FE}}z^{2}D_{P}^{2}\sin \unicode[STIX]{x1D717}\,\text{d}\unicode[STIX]{x1D717}}{\displaystyle \int _{0}^{\unicode[STIX]{x03C0}-\unicode[STIX]{x1D712}}\dot{\unicode[STIX]{x1D6FE}}D_{P}^{2}\sin \unicode[STIX]{x1D717}\,\text{d}\unicode[STIX]{x1D717}}=\frac{D_{P}}{3}(1+\cos \unicode[STIX]{x1D712}),\end{eqnarray}$$
where
$z=((D_{P}/2)\cos \unicode[STIX]{x1D717}-(D_{P}/2)\cos (\unicode[STIX]{x03C0}-\unicode[STIX]{x1D712}))$
in polar coordinates. The lever arm of the friction force can be computed from (4.14). Both, friction force and its corresponding lever arm to the axis of rotation (
$A$
), depend on the angle of repose. With the lines of action of pressure and friction forces, the lever arm of the drag force reads
$L_{D}=(L_{D\unicode[STIX]{x1D70F}}2b_{\unicode[STIX]{x1D70F}}+L_{Dp}b_{p})/(2b_{\unicode[STIX]{x1D70F}}+b_{p})$
. It can be computed analytically as a function of angle
$\unicode[STIX]{x1D712}$
and, hence, of the angle of repose. Substituting (4.6) and (4.13) into (4.5) and taking into account the definition of the Shields number, (2.2), one obtains the following general expression for the critical Shields number:
$$\begin{eqnarray}\unicode[STIX]{x1D703}_{C}=\frac{1}{3f}\frac{1}{\cos \unicode[STIX]{x1D6FD}}\frac{L_{G}}{L_{D}}\frac{1}{\displaystyle \left(2-\frac{\unicode[STIX]{x1D712}}{\unicode[STIX]{x03C0}}+\cos \unicode[STIX]{x1D712}+\frac{\sin (2\unicode[STIX]{x1D712})}{2\unicode[STIX]{x03C0}}\right)}.\end{eqnarray}$$
The ratio
$L_{G}/L_{D}$
and the angle
$\unicode[STIX]{x1D712}$
depend on the substrate geometry. For spheres of same size, the ratio
$L_{G}/L_{D}$
is a function of the angle of repose only. In what follows, we determine the critical Shields number for the quadratic arrangements of different spacing as considered in numerics. For this substrate, the orientation angle,
$\unicode[STIX]{x1D6FD}$
ranges from 0 to
$\pm \unicode[STIX]{x03C0}/4$
.
Figure 10. Schematic representation of different scenarios to define the theoretical zero bed level and the corresponding acting forces.
Since the submerged weight acts on the centre of gravity of the sphere, the lever arm
$L_{G}$
is
$$\begin{eqnarray}L_{G}^{\unicode[STIX]{x1D6F1}}=\frac{(D_{P}+a)}{4}=\frac{1}{2}\frac{\tan \unicode[STIX]{x1D711}}{\sqrt{1+2\tan ^{2}\unicode[STIX]{x1D711}}}D_{P}.\end{eqnarray}$$
The superscript indicates quadratic configuration,
$\unicode[STIX]{x1D6F1}$
. Similarly, the lever arm of the pressure force is (for derivation see appendix A)
$$\begin{eqnarray}L_{Dp}^{\unicode[STIX]{x1D6F1}}=\frac{L_{G}^{\unicode[STIX]{x1D6F1}}}{\tan \unicode[STIX]{x1D711}}=\frac{1}{2}\frac{D_{P}}{\sqrt{1+2\tan ^{2}\unicode[STIX]{x1D711}}}.\end{eqnarray}$$
Since the line of action of the friction force depends on the bead surface exposed to the flow, its lever arm,
$L_{D\unicode[STIX]{x1D70F}}$
depends on the effective zero level. For the quadratic arrangement, the effective zero level depends linearly on the substrate spacing (see figure 8
b). For the quadratic substrate with different spacing, the angle
$\unicode[STIX]{x1D712}$
is coupled to the spacing and hence to the angle of repose as follows (for derivation see appendix A):
$$\begin{eqnarray}\cos \unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6F1}}=\frac{2}{\sqrt{1+2\tan ^{2}\unicode[STIX]{x1D711}}}+2\left(\frac{z_{0}}{D_{P}}\right)-1.\end{eqnarray}$$
On the other hand, from figure 10, it can be shown that the lever arm of the friction force can be expressed by
$$\begin{eqnarray}L_{D\unicode[STIX]{x1D70F}}=h_{eff}+p_{2}-p_{1}-z_{0},\end{eqnarray}$$
where
$p_{1}$
is the protrusion of the single sphere and
$p_{2}$
is the height of the top of the mobile bead over the axis of rotation
$A$
, see figure 10. The lever arm of the drag force is given by
$$\begin{eqnarray}\displaystyle L_{D}^{\unicode[STIX]{x1D6F1}} & = & \displaystyle \frac{D_{P}}{2}\left\{\frac{(1+\cos \unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6F1}})\left[\displaystyle \frac{2}{3}(1+\cos \unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6F1}})-\frac{1}{\displaystyle \sqrt{1+2\tan ^{2}\unicode[STIX]{x1D711}}}+1-2\left(\displaystyle \frac{z_{0}}{D_{P}}\right)\right]}{2-\displaystyle \frac{\unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6F1}}}{\unicode[STIX]{x03C0}}+\cos \unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6F1}}+\displaystyle \frac{\sin (2\unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6F1}})}{2\unicode[STIX]{x03C0}}}\right\}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{D_{P}}{2}\left\{\frac{\displaystyle \frac{1}{\sqrt{1+2\tan ^{2}\unicode[STIX]{x1D711}}}\left(1-\displaystyle \frac{\unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6F1}}}{\unicode[STIX]{x03C0}}+\displaystyle \frac{\sin (2\unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6F1}})}{2\unicode[STIX]{x03C0}}\right)}{2-\displaystyle \frac{\unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6F1}}}{\unicode[STIX]{x03C0}}+\cos \unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6F1}}+\displaystyle \frac{\sin (2\unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6F1}})}{2\unicode[STIX]{x03C0}}}\right\}\end{eqnarray}$$
(for derivation see appendix A). In the quadratic configuration, the effective zero level depends on the spacing and therefore on the angle of repose. Accordingly, (4.15) remains as a sole function of the angle of repose and the orientation angle
$\unicode[STIX]{x1D6FD}$
, and can be computed analytically.
After the corresponding trigonometric derivation (see appendix A), the critical Shields number for quadratically arranged substrates reads
$$\begin{eqnarray}\displaystyle \!\left.\begin{array}{@{}l@{}}\!\!\unicode[STIX]{x1D703}_{C}^{\unicode[STIX]{x1D6F1}}=\displaystyle \frac{1}{3f}\displaystyle \frac{1}{\cos \unicode[STIX]{x1D6FD}}G^{\unicode[STIX]{x1D6F1}}(\unicode[STIX]{x1D711})\quad \text{ with}\\ \!\!\!G^{\unicode[STIX]{x1D6F1}}=\frac{\displaystyle \tan \unicode[STIX]{x1D711}}{\displaystyle \sqrt{1+2\tan ^{2}\unicode[STIX]{x1D711}}(1+\cos \unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6F1}})\left[\displaystyle \frac{2}{3}(1+\cos \unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6F1}})-\displaystyle \frac{1}{\sqrt{1+2\tan ^{2}\unicode[STIX]{x1D711}}}+1-2\left(\displaystyle \frac{z_{0}}{D_{P}}\right)\right]+\left(\displaystyle 1-\displaystyle \frac{\unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6F1}}}{\unicode[STIX]{x03C0}}+\displaystyle \frac{\sin (2\unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6F1}})}{2\unicode[STIX]{x03C0}}\right)}\\ \quad \text{where }\cos \unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6F1}}=\displaystyle \frac{2}{\sqrt{1+2\tan ^{2}\unicode[STIX]{x1D711}}}+2\left(\displaystyle \frac{z_{0}}{D_{P}}\right)-1\\ \quad \text{and }\displaystyle \frac{z_{0}}{D_{P}}=0.077+0.13\left(\displaystyle \frac{a}{D_{P}}\right)=0.077+0.13\left(\displaystyle \frac{2\tan ^{2}\unicode[STIX]{x1D711}}{\sqrt{1+2\tan ^{2}\unicode[STIX]{x1D711}}}-1\right).\end{array}\right\} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
Analogous to the quadratic geometry, for triangular substrates with spacing
$a$
between spheres, it can be shown that the critical Shields number is given by the following expression:
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}l@{}}\!\!\!\unicode[STIX]{x1D703}_{C}^{\unicode[STIX]{x1D6E5}}=\displaystyle \frac{1}{3f}\displaystyle \frac{1}{\cos \unicode[STIX]{x1D6FD}}G^{\unicode[STIX]{x1D6E5}}(\unicode[STIX]{x1D711})\quad \text{with}\\ \!\!\!G^{\unicode[STIX]{x1D6E5}}=\frac{\displaystyle \tan \unicode[STIX]{x1D711}}{\displaystyle \sqrt{1+4\tan ^{2}\unicode[STIX]{x1D711}}(1+\cos \unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6E5}})\left[\displaystyle \frac{2}{3}(1+\cos \unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6E5}})-\displaystyle \frac{\displaystyle 1}{\sqrt{1+4\tan ^{2}\unicode[STIX]{x1D711}}}+1-2\left(\displaystyle \frac{z_{0}}{D_{P}}\right)\right]+\left(\displaystyle 1-\displaystyle \frac{\unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6E5}}}{\unicode[STIX]{x03C0}}+\displaystyle \frac{\sin (2\unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6E5}})}{2\unicode[STIX]{x03C0}}\right)}\\ \quad \text{where }\cos \unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6E5}}=\displaystyle \frac{2}{\sqrt{1+4\tan ^{2}\unicode[STIX]{x1D711}}}+2\left(\displaystyle \frac{z_{0}}{D_{P}}\right)-1,\end{array}\right\} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
where the superscript
$\unicode[STIX]{x1D6E5}$
indicates the triangular configuration. Further details on the trigonometric derivation are provided in the appendix A.
5 Discussion
We studied the incipient particle motion on regular substrates in shear flows at low particle Reynolds numbers numerically and proposed a model to describe the onset. The flow-induced forces, their dependence on the substrate geometry and the effective zero level were determined numerically. To ensure grid independence we performed a grid convergence study as shown in figure 5 and calculated the grid convergence index (Roache Reference Roache1994, Reference Roache1997, Reference Roache1998). This showed that a grid resolution with approximately 29 million elements is fine enough to determine critical Shields numbers to within 1 %. Using the case of a sphere on a plane surface and of a hemisphere as benchmark confirmed this result.
In numerics, we found that the lift force is smaller by several orders of magnitude than the drag force at the particle Reynolds number studied. This is in line with previous studies by Leighton & Acrivos (Reference Leighton and Acrivos1985) and King & Leighton (Reference King and Leighton1997) for flat walls and also in agreement with experiments performed with an irregular sediment bed driven by a fluid flow in a cylindrical channel under laminar flows (Hong et al. Reference Hong, Tao and Kudrolli2015). The ratio of the two forces groups around the value for the flat plane. It is also consistent with the critical Shields number being independent of the particle Reynolds number as shown experimentally by Agudo & Wierschem (Reference Agudo and Wierschem2012). Hence, lift forces may be neglected.
While for more than one bead sliding and rolling motion has been observed (Agudo et al.
Reference Agudo, Dasilva and Wierschem2014), numerics show that rolling is always the preferred mechanism of motion for single spheres independent of substrate spacing. In accordance with previous experiments (Agudo et al.
Reference Agudo, Dasilva and Wierschem2014), this holds true even for high angles of repose and hidden spheres, since the line of action is always above the centre of gravity. The numerical calculations show a flow penetration of approximately 8–13 % of the diameter that increases linearly with the substrate spacing, see figure 8(b). The dependency of the effective zero level on the geometry has not been taken into account so far in models for incipient motion. Only a constant level of flow penetration has been considered in models at high Reynolds numbers: while Kirchner et al. (Reference Kirchner, Diertrich, Iseya and Ikeda1990) assumed a constant flow penetration of
$0.16\cdot D_{P}$
, others based on previous studies by van Rijn (Reference van Rijn1984), considered one of
$0.25\cdot D_{P}$
(Cheng & Chiev Reference Cheng and Chiev1998; Dey Reference Dey1999; Wu & Chou Reference Wu and Chou2003; Vollmer & Kleinhans Reference Vollmer and Kleinhans2007). In an irregularly arranged monolayer of spheres, Derksen & Larsen (Reference Derksen and Larsen2011) observed a decrease of flow penetration with the surface occupancy at particle Reynolds numbers similar to ours. Their values for the flow penetration are in good agreement with our numerical results: at same surface occupancy, they deviate only by approximately 5 %. The agreement between the results for the irregularly arranged monolayer and our regular substrates shows that our data can be extrapolated to irregular structures without great error. In an irregular arrangement of mobile beads submitted to a Couette flow, Ouriemi et al. (Reference Ouriemi, Aussillous and Guazzelli2009) numerically found a flow penetration of the order of the particle diameter at Shields numbers of 0.89, well above the critical Shields number which was approximately 0.06 in their study. The order of the flow penetration was in good agreement with experimental results obtained by Goharzadeh et al. (Reference Goharzadeh, Khalili and Jorgensen2005). Yet, apart from the rather large supercritical Shields numbers, we suppose that this discrepancy to Derksen & Larsen (Reference Derksen and Larsen2011) results and to ours is mainly due to the moving beads in the mobile bed, which transfer fluid motion to the downside of the beads. At onset, however, only single events have been observed, see e.g. Charru et al. (Reference Charru, Mouilleron and Eiffel2004), Wierschem et al. (Reference Wierschem, Groh, Rehberg, Aksel and Kruelle2008) or Seizilles et al. (Reference Seizilles, Lajeunesse, Devauchelle and Bak2014).
As pointed out previously, the numerical calculations may differ from the real critical Shields numbers by approximately 1 %. Figure 11 compares the critical Shields numbers from numerics to those of our model according to (4.21). For the model, we used the numerically calculated effective zero level from figure 8(b) and its average value of
$0.105\cdot D_{P}$
. The model is in very good agreement with numerics. The inset of figure 11 shows that deviations of the model to numerics are less than 10 %. The discrepancy may have several reasons: First, the characteristic shear rate in numerics is equal to the velocity of the top plate over the width between substrate crests and top plate. Hence, the Shields number does not take into account flow penetration into the substrate. For the model, this results in Shields numbers that are approximately 3 % lower than in numerics. Second, the model disregards the perturbed flow around the sphere. Furthermore, the model does not take into account the flow through the voids underneath the mobile bead (see figure 9). This flow creates drag that opposes the torque created by the main flow and hence results in a higher critical Shields number. Comparison between the volume flux underneath the mobile bead, see figure 9(c), and the trend of the deviation in the inset of figure 11 shows remarkable similarity between the two, indicating that this may be the main contribution depending on spacing.
Figure 11. Critical Shields number for the quadratic arrangement as a function of the spacing between substrate spheres. Solid and open squares represent the numerical values for frictionless sliding and for rolling motion, respectively. The (blue) solid line shows the results of the model (4.21) with linear dependency of the effective zero level on spacing. The (red) dotted line shows (4.21) assuming an average effective zero level of
$0.105\cdot D_{P}$
. Grey squares indicate the experimental results (Agudo & Wierschem Reference Agudo and Wierschem2012; Agudo et al.
Reference Agudo, Dasilva and Wierschem2014). Solid and open symbols in the inset depict the deviation of the model with linear dependency of the effective zero level on spacing and average effective zero level of
$0.105\cdot D_{P}$
to numeric. The sketches depict the dynamic position for a single mobile bead surrounded by immobile particles resting on the quadratic substrate with a spacing of
$0.035\cdot D_{P}$
.
Figure 12. Critical Shields number as a function of the angle of orientation for a single bead placed on a triangular substrate (a) and on a quadratic substrate with relative spacing of 0.035 (b). The experimental results are obtained from Agudo & Wierschem (Reference Agudo and Wierschem2012) at particle Reynolds numbers of approximately 0.005. Blue dashed and dotted lines in (a) (4.22) for triangular substrates with
$z_{0}/D_{P}=0.08$
and 0.13, respectively. Blue solid line in (b) (4.21) with linear dependency of the effective zero level on the relative spacing.
The results are also compared to previous experimental studies with uniform substrates at particle Reynolds numbers lower than one (Agudo & Wierschem Reference Agudo and Wierschem2012). These experiments were carried out placing a single sphere of
$(405.9\pm 8.7)~\unicode[STIX]{x03BC}\text{m}$
diameter on regular substrates made of spheres of same size. The substrates were arranged in triangular and quadratic order. For the latter, exposure and angle of repose was varied by using substrates with three different spacing between the beads,
$a/D_{P}$
: 0.035, 0.230 and 0.268 (Agudo & Wierschem Reference Agudo and Wierschem2012; Agudo et al.
Reference Agudo, Dasilva and Wierschem2014). To build the quadratic substrates, the glass spheres were deposited and glued onto stainless steel wire sieves of different mesh sizes and wire diameters. The laminar shear flow was set up and controlled using a parallel-disk configuration in a rotational rheometer. Error bars for spacing in the experimental data are added from the data provided by Agudo & Wierschem (Reference Agudo and Wierschem2012). Apart from deviations in the wire mesh size, other effects derived from the substrate production process, such as variation in height of the substrate, are expected to contribute considerably in the absolute uncertainty. The presence of the wire mesh is also expected to reduce the leak flow underneath the mobile bead.
Figure 13. Comparison between the critical Shields numbers, numerics and different models as a function of the angle of repose. Open squares: numerics. Blue solid, dashed, dotted lines: (4.21) with linear dependency of the effective zero level on spacing, (4.22) for triangular substrates with
$z_{0}/D_{P}=0.08$
and 0.13, respectively. Grey triangle and squares show experimental data for single exposed beads (Agudo & Wierschem Reference Agudo and Wierschem2012). Red solid symbols indicate data from different models: Wiberg & Smith (Reference Wiberg and Smith1987) (rhomboids), Ling (Reference Ling1995) (triangle), Dey (Reference Dey1999) (circle), Bravo et al. (Reference Bravo, Ortiz and Perez-Aparicio2014) with
$z_{0}/D_{P}=0.1$
(crosses) and 0.3 (pluses), Ali & Dey (Reference Ali and Dey2016) (minus). Red dash-dotted line shows the data from Vollmer & Kleinhans (Reference Vollmer and Kleinhans2007).
The data may be further compared to experiments with a single mobile bead surrounded by immobile particles resting on the quadratic substrate with a spacing of
$0.035\cdot D_{P}$
(Agudo et al.
Reference Agudo, Dasilva and Wierschem2014). For this caged bead, the initial motion ends once the bead touches its neighbour in front. This dynamic equilibrium position, as depicted in the sketch of figure 11, mimics the situation of a single sphere resting on a quadratic substrate with maximum separation between beads (i.e. about
$0.41\cdot D_{P}$
). For a proper comparison, the experimental value is corrected by
$\cos \unicode[STIX]{x03C0}/4$
to take into account the direction of the particle motion with respect to the main flow direction. Similarly, one needs to take into account that the caged bead is placed on a second layer and the gap,
$h$
, is thus slightly different from that in numerics and in the model. The top of the immobile bead forming the second layer is at a distance
$p_{1}$
over the top of the substrate spheres (see sketch of figure 11). Accordingly, the experimental Shields number for the caged bead is multiplied by the factor
$h/(h-p_{1})$
for a proper comparison.
Comparing numerics to experiment in figure 11, we find good accordance particularly at narrow spacings and for caged beads. For the substrates with a relative spacing of
$0.230\cdot D_{P}$
and
$0.268\cdot D_{P}$
, however, the experimental data result in somewhat lower values for the critical Shields number. This deviation may be partially caused for the presence of the sieve in the experimental set-up, which has a similar geometry as an underlying level of spheres. This reduces the leak flow underneath the mobile bead as shown in figure 9(c). Apart from the presence of the sieve, imperfections in substrate preparation and variations in the size of the substrate beads may result in reduced support of the exposed bead.
The model also takes into account the orientation angle between the main flow direction and the symmetry axis of the substrate,
$\unicode[STIX]{x1D6FD}$
. Figure 12 shows a direct comparison between the experimental data for triangular and quadratic substrates (Agudo & Wierschem Reference Agudo and Wierschem2012) with the model according to (4.21) with linear dependency of the effective zero level on spacing and according to (4.22). Since the penetration factor,
$z_{0}/D_{P}$
remains unknown for the triangular symmetry, we consider the two extreme cases of the quadratic configuration,
$z_{0}/D_{P}=0.08$
(dashed line), and 0.13 (dotted line). For a sphere deposited on a closely packed triangular configuration, however, the critical Shields number is hardly affected by the penetration factor (see figure 12
a). Although small deviations with experiments, the analytical model is able to track the dependence of the critical Shields number on the orientation of the substrate.
A comparison of our results to the previously proposed models by Wiberg & Smith (Reference Wiberg and Smith1987), Ling (Reference Ling1995), Dey (Reference Dey1999), Vollmer & Kleinhans (Reference Vollmer and Kleinhans2007), Bravo et al. (Reference Bravo, Ortiz and Perez-Aparicio2014), Ali & Dey (Reference Ali and Dey2016) is shown in figure 13. For ease of comparison, the critical Shields number is plotted as a function of the angle of repose at
$\unicode[STIX]{x1D6FD}$
equal to zero. Apart from the numerical results and the model for a quadratic arrangement, we also plotted the results for the critical Shields number of the triangular configuration. We considered again the two extreme cases found for the quadratic configuration. As shown in figure 13, the critical Shields number for the triangular geometry is strongly affected by the effective zero level particularly at large angles of repose. At small angles of repose, the effective zero level accounts only for a few per cent. Interestingly, the results for quadratic arrangements are within the possible range for the triangular configuration. In particular, at low angles of repose, deviations between quadratic and triangular geometries are rather small. Thus, although the orientation angle of the geometry affects the critical Shields number by up to a factor of 2 in the case of the triangular configuration, if the main flow direction coincides with the bead motion, the critical Shields number for a certain angle of repose does not depend sensitively on the considered symmetry (i.e. triangular or quadratic).
Former models in literature were derived for turbulent flows. Yet, for the case where the particle is located entirely within the viscous sublayer, they can be applied for laminar flows at small particle Reynolds numbers. Except for the models by Vollmer & Kleinhans (Reference Vollmer and Kleinhans2007) and by Ali & Dey (Reference Ali and Dey2016), the drag force is described by a drag coefficient for unbounded flow. All other models result in larger critical Shields numbers than our numerics. In what follows, we discuss possible reasons for discrepancies.
The model by Wiberg & Smith (Reference Wiberg and Smith1987) gives values that are rather close to our numerical data. They used angles of repose between
$50^{\circ }$
and
$60^{\circ }$
as average values for naturally sorted beds of uniformly sized spheres. Like in our model, they neglected lift. They, however, modelled the drag force by using only the exposure degree, i.e. by the cross-sectional area exposed to the flow and did not consider flow penetration into the substrate. Ling (Reference Ling1995) calculated the onset for a sphere deposited on a closely packed triangular configuration. Lift and flow penetration were neglected and the point of attack of the drag force was assumed to coincide with the spheres centre of gravity. As shown in figure 13, the obtained value overestimates the experimental data by almost a factor of 3. Besides neglecting flow penetration, the deviation is apparently mainly due to the smaller lever arm of the drag force.
Different from the aforementioned authors, Dey (Reference Dey1999) and Vollmer & Kleinhans (Reference Vollmer and Kleinhans2007) took into account flow penetration by considering a constant effective zero level of
$0.25\cdot D_{P}$
below the top of the substrate. Similar to our model, they approached the point of attack of the drag force by integrating the flow velocity over the exposed cross-sectional area of the sphere. The red dashed-dotted line in figure 13 is obtained by applying the model proposed by Vollmer & Kleinhans (Reference Vollmer and Kleinhans2007) to a quadratic substrate with different spacings without taking into account turbulent pressure fluctuations or shear stress excursions and assuming a linear velocity profile for particle Reynolds numbers lower than one. For sediment beds of uniformly sized spheres and using the equation for the angle of repose proposed by Ippen & Eagleson (Reference Ippen and Eagleson1955), Dey (Reference Dey1999) used an effective angle of repose of approximately
$28^{\circ }$
. For that value and using the drag coefficient for a freely falling sphere from Morsi & Alexander (Reference Morsi and Alexander1972), his model yields a critical Shields number of approximately 0.2. The lift coefficient was obtained by extensive calibration of the model with previous experimental results. On the other hand, Vollmer & Kleinhans (Reference Vollmer and Kleinhans2007) modelled the drag and lift forces by using respective coefficients experimentally obtained by Coleman (Reference Coleman and Collins1967) for a single sphere placed on a fixed stream bed. Despite the fact that the models by Dey (Reference Dey1999) and by Vollmer & Kleinhans (Reference Vollmer and Kleinhans2007) reflect very well the literature data at higher Reynolds numbers, they deviate considerably from our experimental results (see figure 13) as well as from critical Shields numbers of approximately 0.02–0.04 recently obtained in laminar flows for a single highly exposed bead on a disordered substrate (Charru et al.
Reference Charru, Larrieu, Dupont and Zenit2007). The deviation may be due to the negative lift forces of the order of drag considered by the authors at small particle Reynolds numbers.
Bravo et al. (Reference Bravo, Ortiz and Perez-Aparicio2014) also considered a flow penetration between
$0.1\cdot D_{P}$
and
$0.3\cdot D_{P}$
in their two-dimensional model. Lift forces were addressed with Saffman’s solution for a sphere in a low shear flow (Saffman Reference Saffman1965). The hydrodynamic forces, however, were assumed to act on the centre of gravity. This may explain overestimating the experimental data. Recently, Ali & Dey (Reference Ali and Dey2016) proposed a model addressing in detail the point of attack of hydrodynamic forces based on the proper geometrical analysis of a sphere placed on a closely packed triangular sediment bed. The effective zero level was set constant at
$0.21\cdot D_{P}$
below the top of the substrate according to Dey et al. (Reference Dey, Das, Gaudio and Bose2012). While the drag coefficient was obtained from Coleman’s experiments (Coleman Reference Coleman and Collins1967) inertial lift forces were evaluated with a lift coefficient of
$0.85\cdot C_{D}$
from the study by Chepil (Reference Chepil1961). Similar to us, Ali & Dey (Reference Ali and Dey2016) incorporated a correction factor lower than one into their model to account for partial particle burial. The correction factor was calibrated with available experimental data for irregularly arranged sediment beds of uniform size. At small particle Reynolds numbers, the model yields a critical Shields number of approximately 0.15. The large deviation with respect to our data for the same geometrical configuration may be due to the contribution of lift forces of the order of drag.
For granular beds, Charru et al. (Reference Charru, Mouilleron and Eiffel2004) showed that the number of moving particles decreases to a saturated value within finite time in laminar flow. The minimum Shields number to obtain a non-zero saturated value was found to be 0.12 independent of the initial bed geometry. Similar values have been reported in other studies in laminar flows (Loiseleux et al. Reference Loiseleux, Gondret, Rabaud and Doppler2005; Ouriemi et al. Reference Ouriemi, Aussillous, Medale, Peysoon and Guazzelli2007, Reference Ouriemi, Aussillous and Guazzelli2009; Derksen Reference Derksen2011; Seizilles et al. Reference Seizilles, Lajeunesse, Devauchelle and Bak2014). As stated by Charru et al. (Reference Charru, Mouilleron and Eiffel2004), exposed particles move at lower Shields numbers and get stuck in caves of the granular beds. In this armoured bed, the individual beads are then tightly packed and need to overcome high angles of repose. Other authors have reported higher critical Shields numbers for low particle Reynolds numbers at saturated conditions (Lobkovsky et al. Reference Lobkovsky, Orpe, Molloy, Kudrolli and Rothman2008; Hong et al. Reference Hong, Tao and Kudrolli2015). These authors defined the threshold condition by measuring the depth of the fluid layer after cessation of granular flow and calculated the Shields number by computing the solution of shear stress for laminar conditions in a pipe of rectangular cross-section. While Lobkovsky et al. (Reference Lobkovsky, Orpe, Molloy, Kudrolli and Rothman2008) found a non-zero saturated value of approximately 0.3, Hong et al. (Reference Hong, Tao and Kudrolli2015) observed Shields numbers ranging between 0.10 and 0.37 depending on the bed armouring level and grain size.
Since in disordered configurations, triangular structures are the most probable ones and all orientation angles are equally possible, one may consider the most unstable triangular configuration being the relevant one. Assuming for the caged beads angles of repose of above
$60^{\circ }$
and penetration factors of 0.08 to 0.13, (4.22) yields critical Shields numbers between 0.10 and 0.12 (see figure 13), very similar to those obtained for a saturated erodible bed by Charru et al. (Reference Charru, Mouilleron and Eiffel2004), Ouriemi et al. (Reference Ouriemi, Aussillous, Medale, Peysoon and Guazzelli2007). At these angles of repose, quadratic and triangular model show similar values for the minimum critical Shields number, emphasising the relevance of the studied system also for granular beds. On the other hand, the maximum critical Shields numbers found by Hong et al. (Reference Hong, Tao and Kudrolli2015) in their experimental study of granular beds is slightly smaller than ours close to maximum angle of repose, indicating that Hong et al. (Reference Hong, Tao and Kudrolli2015) indeed reached a very high armouring level.
The relation of our results to disordered configurations is also tightened by the experimental observation by Agudo et al. (Reference Agudo, Dasilva and Wierschem2014) that neighbouring beads do not have any measurable impact if they are more than 3 bead diameters apart at the low particle Reynolds numbers considered here. This shows that the regularity of the substrates considered has only a minor local impact. Variations in the neighbourhood on the level of the substrate, i.e. deviations from the regular arrangement, are supposed to have even less impact. Therefore, our model is expected to give reasonable results also for disordered beds covering the entire range from highly exposed to hidden beads.
6 Conclusions
We have studied the incipient particle motion of a single sphere as function of geometrical properties of the sediment bed at low particle Reynolds numbers. We found a preference of rolling over sliding incipient motion and showed the impact of the substrate geometry on the critical Shields number over the entire range of spacings. Depending linearly on the spacing, the effective zero level of the shear flow penetrates into the substrate by approximately 8 %–13 % of the particle diameter. We propose a model for the incipient motion as an extension of Goldman’s model for a single particle near a plane surface. Lift forces are neglected according to numerical findings. Correcting the drag force by the ratio of the relative area exposed to the shear flow and using the effective zero level from numerics, the model is able to predict the critical Shields number as a function of the sediment bed geometry. It also covers the critical Shields number dependence on the substrate orientation with respect to the main flow direction. It shows very good agreement with numerics and is in good agreement with the experimental data. Yet, although the orientation angle of the geometry affects the critical Shields number, it shows that the minimum critical Shields numbers for a certain angle of repose do not depend sensitively on the considered configuration for well-exposed particles.
Acknowledgements
The authors are thankful to N. Priebe, M. Baumann, F. Rupprecht and F. Siegel for collaborating in setting up the numerics. The support from Deutsche Forschungsgemeinschaft through WI 2672/4-1 and WI 2672/7-1 is gratefully acknowledged.
Appendix A
Figure 14. Single sphere resting on a quadratically arranged substrate with spacing
$a$
between beads (a) and on a triangularly arranged substrate with spacing
$a$
between beads (b).
From figure 14(a) it can be shown for quadratic arrangements that half of the diagonal connecting centres of opposed substrate beads is
$W=\sqrt{2}/2(D_{P}+a)$
. The distance
$T=\sqrt{D_{P}^{2}-W^{2}}=\sqrt{D_{P}^{2}-(D_{P}+a)^{2}/2}$
and the tangent of the angle of repose is given by
$$\begin{eqnarray}\tan \unicode[STIX]{x1D711}=\frac{D_{P}+a}{2T}=\frac{D_{P}+a}{2\sqrt{D_{P}^{2}-\displaystyle \frac{(D_{P}+a)^{2}}{2}}}.\end{eqnarray}$$
Rearranging (A 1) yields
$$\begin{eqnarray}\frac{D_{P}+a}{D_{P}}=\frac{2\tan \unicode[STIX]{x1D711}}{\sqrt{1+2\tan ^{2}\unicode[STIX]{x1D711}}}.\end{eqnarray}$$
For a bead resting on a triangular substrate with spacing
$a$
as depicted in figure 14(b), the radius of circumscribed and inscribed circle of the equilateral triangle with vertices
$D_{P}+a$
, is given by
$W=\sqrt{3}/3(D_{P}+a)$
and
$w=\sqrt{3}/6(D_{P}+a)$
, respectively. Accordingly, the vertical projection of the distance between the substrate sphere centre
$(O^{\prime })$
and the bead centre
$(O)$
, defined as
$T$
, results in
$T=\sqrt{D_{P}^{2}-(D_{P}+a)^{2}/2}$
and the tangent of the angle of repose yields
$$\begin{eqnarray}\tan \unicode[STIX]{x1D711}=\frac{w}{T}=\frac{\sqrt{3}}{6}\frac{D_{P}+a}{\sqrt{D_{P}^{2}-\displaystyle \frac{(D_{P}+a)^{2}}{3}}}.\end{eqnarray}$$
Rearranging (A 3), we obtain
$$\begin{eqnarray}\frac{D_{P}+a}{D_{P}}=2\sqrt{3}\frac{\tan \unicode[STIX]{x1D711}}{\sqrt{1+4\tan ^{2}\unicode[STIX]{x1D711}}}.\end{eqnarray}$$
The lever arm of gravity and buoyancy forces acts on the bead centre
$(O)$
. For the quadratic substrate, it is given by
$L_{G}^{\unicode[STIX]{x1D6F1}}=(D_{P}+a)/4$
. Substituting (A 2) yields
$$\begin{eqnarray}L_{G}^{\unicode[STIX]{x1D6F1}}=\frac{1}{2}\frac{\tan \unicode[STIX]{x1D711}}{\sqrt{1+2\tan ^{2}\unicode[STIX]{x1D711}}}D_{P}.\end{eqnarray}$$
For the triangular substrate, the lever arm is given by
$L_{G}^{\unicode[STIX]{x1D6E5}}=w/2=\sqrt{3}/12(D_{P}+a)$
. With (A 4) this yields
$$\begin{eqnarray}L_{G}^{\unicode[STIX]{x1D6E5}}=\frac{1}{2}\frac{\tan \unicode[STIX]{x1D711}}{\sqrt{1+4\tan ^{2}\unicode[STIX]{x1D711}}}D_{P}.\end{eqnarray}$$
In the case of zero spacing this reduces to
$L_{G}^{\unicode[STIX]{x1D6E5}}=(\sqrt{3}/12)D_{P}$
.
The lever arm of the drag force acts at a distance
$h_{eff}$
from the effective zero level (see figure 10). Therefore,
$L_{D}=h_{eff}+p_{2}-p_{1}-z_{0}$
with
$h_{eff}=(D_{P}/3)(1+\cos \unicode[STIX]{x1D712})$
, see (4.14). From figure 15, it can be readily seen that
$p_{1}=T$
,
$p_{2}=(T+D_{P})/2$
. and that
$T=(D_{P}+a)/(2\tan \unicode[STIX]{x1D711})=D_{P}/\sqrt{1+2\tan ^{2}\unicode[STIX]{x1D711}}$
. Similarly, for the triangular arrangement,
$T=w/(2\tan \unicode[STIX]{x1D711})=D_{P}/\sqrt{1+4\tan ^{2}\unicode[STIX]{x1D711}}$
. Using these relationships, one finally arrives at the following expressions for the lever arm of the viscous force:
$$\begin{eqnarray}L_{D\unicode[STIX]{x1D70F}}^{\unicode[STIX]{x1D6F1}}=\frac{D_{P}}{2}\left[\frac{2}{3}(1+\cos \unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6F1}})-\frac{1}{\sqrt{1+2\tan ^{2}\unicode[STIX]{x1D711}}}+1-2\left(\frac{z_{0}}{D_{P}}\right)\right]\end{eqnarray}$$
and
$$\begin{eqnarray}L_{D\unicode[STIX]{x1D70F}}^{\unicode[STIX]{x1D6E5}}=\frac{D_{P}}{2}\left[\frac{2}{3}(1+\cos \unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6E5}})-\frac{1}{\sqrt{1+4\tan ^{2}\unicode[STIX]{x1D711}}}+1-2\left(\frac{z_{0}}{D_{P}}\right)\right],\end{eqnarray}$$
for the triangular substrate.
Figure 15. Schematic representation of forces acting on a single bead resting on a quadratic (a) and on a triangular (b) substrate on varying the gap
$a$
between spheres.
The level arm of the pressure force acts on the bead centre (
$O$
). Accordingly,
$L_{Dp}=L_{D}/\tan \unicode[STIX]{x1D711}$
. Applying (A 5) and (A 6) for the lever arm of gravity and buoyancy forces, one obtains
$$\begin{eqnarray}L_{Dp}^{\unicode[STIX]{x1D6F1}}=\frac{1}{2}\frac{D_{P}}{\sqrt{1+2\tan ^{2}\unicode[STIX]{x1D711}}}\end{eqnarray}$$
and
$$\begin{eqnarray}L_{Dp}^{\unicode[STIX]{x1D6E5}}=\frac{1}{2}\frac{D_{P}}{\sqrt{1+4\tan ^{2}\unicode[STIX]{x1D711}}}.\end{eqnarray}$$
The total level arm is defined as
$L_{D}=(L_{D\unicode[STIX]{x1D70F}}2b_{\unicode[STIX]{x1D70F}}+L_{Dp}b_{p})/(2b_{\unicode[STIX]{x1D70F}}+b_{p})$
where
$b_{\unicode[STIX]{x1D70F}}$
and
$b_{p}$
are the viscous and the pressure correcting factors as defined in (4.10) and (4.11), respectively. Substituting (4.10), (4.11), (A 7) and (A 9) in the expression for the total level arm, one finally arrives to (4.20). Analogous to the quadratic case, substituting now (4.10), (4.11), (A 8) and (A 10) in the expression for the total arm yields
$$\begin{eqnarray}\displaystyle L_{D}^{\unicode[STIX]{x1D6E5}} & = & \displaystyle \frac{D_{P}}{2}\left\{\frac{(1+\cos \unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6E5}})\left[\displaystyle \frac{2}{3}(1+\cos \unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6E5}})-\frac{1}{\sqrt{1+4\tan ^{2}\unicode[STIX]{x1D711}}}+1-2\left(\frac{z_{0}}{D_{P}}\right)\right]}{2-\displaystyle \frac{\unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6E5}}}{\unicode[STIX]{x03C0}}+\cos \unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6E5}}+\frac{\sin (2\unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6E5}})}{2\unicode[STIX]{x03C0}}}\right\}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{D_{P}}{2}\left\{\frac{\displaystyle \frac{1}{\sqrt{1+4\tan ^{2}\unicode[STIX]{x1D711}}}\left(1-\frac{\unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6E5}}}{\unicode[STIX]{x03C0}}+\frac{\sin (2\unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6E5}})}{2\unicode[STIX]{x03C0}}\right)}{2-\displaystyle \frac{\unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6E5}}}{\unicode[STIX]{x03C0}}+\cos \unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6E5}}+\frac{\sin (2\unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6E5}})}{2\unicode[STIX]{x03C0}}}\right\}\end{eqnarray}$$
for the triangular case.
From figures 2 or 15, it can be shown that the level arm of drag force over the rotation axis in the bead centre, (
$O$
),
$l_{D}$
, is related to the level arms over the rotation axis in (
$A$
),
$L_{D}$
, as
$l_{D}=L_{D}-(L_{G}/\tan \unicode[STIX]{x1D711})$
.
The angle
$\unicode[STIX]{x1D712}$
, which is related to the effective zero level, is coupled to the angle of repose as follows (see figure 15):
$$\begin{eqnarray}\displaystyle \cos \unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6F1}} & = & \displaystyle \frac{2s}{D_{P}}=\frac{2}{\sqrt{1+2\tan ^{2}\unicode[STIX]{x1D711}}}+2\left(\frac{z_{0}}{D_{P}}\right)-1\nonumber\\ \displaystyle & & \displaystyle \quad \text{with }s=T-\frac{D_{P}}{2}+z_{0}=\frac{D_{P}}{2}\left[\frac{2}{\sqrt{1+2\tan ^{2}\unicode[STIX]{x1D711}}}+2\left(\frac{z_{0}}{D_{P}}\right)-1\right]\quad\end{eqnarray}$$
for the quadratic substrate and
$$\begin{eqnarray}\displaystyle \cos \unicode[STIX]{x1D712}^{\unicode[STIX]{x1D6E5}} & = & \displaystyle \frac{2s}{D_{P}}=\frac{2}{\sqrt{1+4\tan ^{2}\unicode[STIX]{x1D711}}}+2\left(\frac{z_{0}}{D_{P}}\right)-1\nonumber\\ \displaystyle & & \displaystyle \quad \text{with }s=T-\frac{D_{P}}{2}+z_{0}=\frac{D_{P}}{2}\left[\frac{2}{\sqrt{1+4\tan ^{2}\unicode[STIX]{x1D711}}}+2\left(\frac{z_{0}}{D_{P}}\right)-1\right]\quad\end{eqnarray}$$
for the triangular substrate.
