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

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}$ :

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

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

where $\unicode[STIX]{x1D707}$ is the friction coefficient. In this case, the particle density reads

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 }$ :

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

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

The effective weight acting on the centre of gravity of the sphere is given by

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

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

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)

and the corresponding correction factor for the friction force yields

where $S$ is the bead surface. Similarly, the cross-sectional area exposed to the flow is

Hence, the main Shielding factor yields

Substituting (4.12) into (4.8), the drag force exerted by the shear flow results in

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):

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:

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

The superscript indicates quadratic configuration, $\unicode[STIX]{x1D6F1}$ . Similarly, the lever arm of the pressure force is (for derivation see appendix A)

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):

On the other hand, from figure 10, it can be shown that the lever arm of the friction force can be expressed by

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

(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

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:

where the superscript $\unicode[STIX]{x1D6E5}$ indicates the triangular configuration. Further details on the trigonometric derivation are provided in the appendix A.