PHYSICAL REVIEW E 92, 042202 (2015)

Onset and cessation of motion in hydrodynamically sheared granular beds Abram H. Clark,1 Mark D. Shattuck,2 Nicholas T. Ouellette,1,3 and Corey S. O’Hern1,4,5 1

Department of Mechanical Engineering and Materials Science, Yale University, New Haven, Connecticut 06520, USA 2 Benjamin Levich Institute and Physics Department, The City College of the City University of New York, New York, New York 10031, USA 3 Department of Civil and Environmental Engineering, Stanford University, Stanford, California 94305, USA 4 Department of Physics, Yale University, New Haven, Connecticut 06520, USA 5 Department of Applied Physics, Yale University, New Haven, Connecticut 06520, USA (Received 13 April 2015; published 16 October 2015) We performed molecular dynamics simulations of granular beds driven by a model hydrodynamic shear flow to elucidate general grain-scale mechanisms that determine the onset and cessation of sediment transport. By varying the Shields number (the nondimensional shear stress at the top of the bed) and particle Reynolds number (the ratio of particle inertia to viscous damping), we explore how variations of the fluid flow rate, particle inertia, and fluid viscosity affect the onset and cessation of bed motion. For low to moderate particle Reynolds numbers, a critical boundary separates mobile and static states. Transition times between these states diverge as this boundary is approached both from above and below. At high particle Reynolds number, inertial effects become dominant, and particle motion can be sustained well below flow rates at which mobilization of a static bed occurs. We also find that the onset of bed motion (for both low and high particle Reynolds numbers) is described by Weibullian weakest-link statistics and thus is crucially dependent on the packing structure of the granular bed, even deep beneath the surface. DOI: 10.1103/PhysRevE.92.042202

PACS number(s): 45.70.Ht, 47.57.Gc, 92.10.Wa, 92.40.Gc

I. INTRODUCTION

Fluid flowing laterally over a granular bed exerts shear stress on the grains. This occurs in many natural settings and industrial applications, such as sediment transport in riverbeds [1,2] and slurries pipes [3,4]. The ratio of the shear stress exerted by the fluid on the top of the bed to the buoyancycorrected particle weight is known as the Shields number  [5]. For small , no grain motion occurs; at sufficiently large , however, grains can be entrained by the flow [6–15]. Despite decades of research, the nature of the transition between static and mobilized granular beds is not well understood. The geometric structure of the contact network in the bed determines its mechanical strength [16–18]. Bed mobilization is also strongly affected by the complex and unsteady fluid flow above the bed, as well as how strongly the fluid flow couples to the grains, as quantified by the particle Reynolds number Rep [5,7,10,11], which measures how quickly grains equilibrate to the fluid flow. Weak stresses applied to the interior of the bed by fluid flowing through the pore spaces between grains may also play a role in bed mobilization [19–21]. Although empirical hydraulic models capture some important aspects of sediment transport problems [8], there is at present no fundamental understanding of the relative contributions of these effects on the onset and cessation of grain motion. In this paper, we study a simplified model of a fluid-driven granular bed to clarify the essential physics at the onset of bed motion. In particular, we seek to understand the nature of the mobile-to-static and static-to-mobile transitions as a function of  and Rep and to predict the parameter regime where hysteresis, defined as a finite difference between 0 , above which a static bed will begin to move, and c , below which a mobile system will come to rest, occurs. We performed molecular dynamics (MD) simulations of a two-dimensional (2D) system composed of frictionless disks subjected to a simplified fluid flow that decays from a large 1539-3755/2015/92(4)/042202(7)

value above the bed to a small value inside the bed. Although our model is highly simplified, with, for example, no explicit unsteadiness in the flow or friction between the grains, we find that c (Rep ) from the simulations is consistent with the behavior obtained from a large collection of experiments on sediment transport [7,10,11]. In particular, we find plateau values lc and hc at low and high Rep , with lc > hc , and an intermediate Rep regime that connects the two limiting values. In the low Rep limit, there is a sharp transition at c between mobile and static beds in the infinite-time and infinite-system-size limits with no hysteresis. In the large Rep limit, we find significant hysteresis, since particle inertia can sustain motion well below the 0 at which bed motion is initiated. We also find that the onset of bed motion at  > c for low Rep and  > 0 for high Rep depends strongly on system size and exhibits weakest-link statistics [22,23]. Thus, the onset of bed motion in our system depends on the bed packing structure, even deep beneath the surface. II. DETAILS OF THE MODEL

We study a domain of width W that contains N/2 large and N/2 small disks with diameter ratio 1.4. There is no upper boundary, and the lower boundary is rigid with infinite friction so that the horizontal velocities of all particles touching it are set to zero. We use periodic boundary conditions in the horizontal direction. The total force on each particle is given by the vector sum of contact forces from other particles, a gravitational force, and a Stokes-drag-like force from a fluid that moves horizontally:  mi ai = (1) Fijc − mi g  yˆ + Bi [v0 f (r )xˆ − vi ]. j

Here, mi ∝ Di2 is the particle mass, Di is the diameter of particle i, vi and ai are the velocity and acceleration of particle

042202-1

©2015 American Physical Society

CLARK, SHATTUCK, OUELLETTE, AND O’HERN

(a)

φ i = 0.48

vf

PHYSICAL REVIEW E 92, 042202 (2015)

written as

yˆ x ˆ

i

=

Bv0 , mg 

B/m =√  . g /D i

(b) v g (×3)

vf

FIG. 1. (Color online) Layer-averaged fluid velocity vf versus depth for static (a) and mobile (b) beds. Grains are subjected to a ˆ The gray circle (a) fluid drag force and a gravitational force −mg  y. defines the area used to calculate the local packing fraction φi near the ith particle, which determines the local fluid velocity. We also show the layer-averaged grain velocity vg for a mobile bed in (b).

(2)

(3)

The Shields number  gives the dimensionless shear force at the top of a√ static bed.  is the ratio of the gravitational settling time τs = D/g  to the viscous time scale m/B. Since the particle Reynolds number Rep = v0νD , where ν is the kinematic viscosity, and the Stokes drag is proportional to ρf νD, where ρ 0 ∝ ρfg Rep (ρg is ρf is the fluid density, the ratio 2 = mv BD the mass density of the grains) compares the inertia of grains entrained in the flow to the strength of the viscous drag. To characterize flow onset and cessation in our system, we employed two protocols. In protocol A, to study the mobile-tostatic transition, we distributed particles randomly throughout the domain and set a constant value of  for a total time of roughly 105 τs . We consider the bed to be at rest when the maximum net particle acceleration amax is below a threshold athresh roughly one order of magnitude smaller than g  and roughly three orders of magnitude smaller than typical values for a moving bed. In protocol B, to understand the dynamics of the static-to-mobile transition, we begin with a static bed from protocol A and slowly increase  in increments  = 0.01. If amax < athresh after roughly one intergrain collision time, then  is increased. If amax > athresh , we keep  constant until amax < athresh . We designate a system as mobile under protocol B with slightly different criteria: amax /athresh > 10 and v¯g > 0.04Vs , where√v¯g is the average horizontal velocity of all grains and Vs = g  D is the settling velocity. These thresholds filter out small rearrangement events, keeping only states with substantial grain motion. III. RESULTS AND DISCUSSION

i, mi g  is the buoyancy-corrected particle weight, Bi ∝ Di sets the drag on disk i, v0 is a characteristic fluid velocity at the surface of a static bed, and f (r ) is the fluid velocity at r. Fijc = r r K(1 − Dijij )θ (1 − Dijij )ˆrij is the pairwise repulsive contact force on disk i from disk j , where K is the particle stiffness, rij is the separation between the centers of the particles, Dij = (Di + Dj )/2, rˆij is the unit vector connecting their centers, and θ is the Heaviside step function. f (r ), the fluid velocity profile, varies smoothly from a large value above the bed to a small value inside the bed. We choose a form that depends only on the local packing fraction φi : f (φi ) = e−b(φi −0.5) , where b controls the ratio of the magnitude of the fluid flow above and inside the bed. φi is calculated in a small circular region with diameter Di + 2Dl around each particle, as shown in Fig. 1(a). We note that f = 1 for φi = 0.5, a typical value at the bed surface. See Sec. III A for force profiles in static and mobile states. Three nondimensional numbers govern the behavior of K 3 Eq. (1). We set the nondimensional stiffness mg  > 3 × 10 to be sufficiently large that increasing it has no effect on our results. The other two nondimensional parameters can be

Figure 2 shows the boundaries between mobile and static beds as a function of  and Rep from simulations with b = 2, 4, and 6. We find a curve c (Rep ) above which the particles are unable to find a stable packing under protocol A. In the low-Rep limit near c , grain motion is highly overdamped, and particles do not leave the bed. As Rep is increased, the inertial effects of mobilized particles striking the bed make finding a stable configuration more difficult, decreasing c significantly. Figure 2 shows parameter values where grain motion does (blue circles) and does not (green squares) stop under protocol A. We then apply protocol B to stopped systems and find another boundary 0 (Rep ) that specifies when grain flow can be initiated. For low Rep , 0 < c , and bed motion initiated near 0 is temporary. Thus, at low Rep , permanent grain motion is initiated at c in the large system limit, as we discuss below. However, in the high-Rep limit where particle inertia is dominant, we observe significant hysteresis: grain motion is initiated at 0 , which is well above the value of c where mobile particles colliding with the bed can sustain bed motion. We also note that this result is consistent with Ref. [24], where, in simulations of Aeolian transport at high Rep , a significant perturbation or lift force was required near

042202-2

ONSET AND CESSATION OF MOTION IN . . .

(a)

temporary grain motion

Θ

0.2

TABLE I. The flow velocity ratios at different heights and the particle Reynolds number when 0 = c . va , v0 , and vb are the fluid velocities above the bed, at the surface of the bed, and in the interior of the bed, respectively.

grains always move

0.8 0.4

PHYSICAL REVIEW E 92, 042202 (2015)

0.1

no grain motion

0.05

10 -2

10 -1

10 0

Θ0

b

va /v0

v0 /vb

va /vb

Rep when 0 = c

b=2

2 4 6

2.23 4.95 11.02

1.97 3.89 7.69

4.39 19.3 84.8

Rep ≈ 50 Rep ≈ 10 Rep ≈ 3

grains never Θc stop

2

10 1

10 2

10 3

10 4

Θ/Γ ∼ Re p (b)

grains always move

0.8 0.4

temporary grain motion

Θ0

Θ

0.2

grains never stop

no grain motion

0.1

Θc

0.05

b=4 10 -2

10 -1

10 0

2

10 1

10 2

10 3

10 4

Θ/Γ ∼ Re p (c)

0.8 0.4

temporary grain motion

grains always move

Θ0

Θ

0.2

low Rep , a decrease in the onset value at moderate Rep , and a lower plateau value at high Rep . The b parameter sets the ratios between the fluid velocity above the bed (va ), at the top of a static bed (v0 ), and in the bulk of the bed (vb ) where grains are packed densely (at φi ≈ 0.84). If a grain is well above the bed, φi ≈ 0.1 (since the grain itself contributes to the local packing fraction), so the ratio of the fluid velocity for this grain to the fluid velocity at the top of an otherwise static bed (where φi ≈ 0.5 and thus f = 1) is va /v0 ≈ e0.4b . Table I shows the ratios va /v0 , v0 /vb , and va /vb for b = 2, 4, and 6. Note that these ratios increase with b, and that varying b changes all ratios simultaneously. Obviously, a fluid flow model with an additional parameter could decouple the va /v0 and v0 /vb ratios, but the model used here was chosen as a very simple way to interpolate between a large fluid velocity above the bed and a small fluid velocity inside the bed. Figure 2 shows that the flow diagram is qualitatively the same for b = 2, 4, and 6. 0 is relatively insensitive to b, but c shifts to smaller Rep and the gap between the plateau values at large and small Rep widens with increasing b.

grains never stop

0.1 0.05

no grain motion 10

-2

10

-1

10

b=6

Θc

0

2

10

1

10

2

10

A. Protocol A: Mobile-to-static transition

3

10

4

Θ/Γ ∼ Re p FIG. 2. (Color online) The flow diagrams, Shields number  versus particle Reynolds number Rep , with (a) b = 2, (b) b = 4, and (c) b = 6. Diagonal lines of data points correspond to lines of constant . The symbols show systems that came to rest ( ) or never stopped ( ) under protocol A, and were permanently ( ) or temporarily (•) mobile as  was increased under protocol B. The dashed line shows c , above which the inertial effects from particles entrained in the flow lead to sustained motion. The solid line indicates 0 , below which the system will never be mobilized. The two large black open circles mark the parameter values we study in Fig. 6.

c to initiate grain motion. Temporary (filled black circles) and permanent (red crosses) motion under protocol B are also marked in Fig. 2. The basic nature of the flow diagram is insensitive to variations in b, although the numerical values for 0 (Rep ) and c (Rep ) change. We note the similarities of c (Rep ) in Fig. 2 to the experimental and observational data compiled in Refs. [7,10,11], even though our model is highly simplified. Both display a plateau in the onset value of  at

Figure 3 shows data from protocol A near c , which we find to be nearly independent of system size, as we discuss below. In Fig. 3(a), we show the time evolution of v¯g multiplied by the fill height N D/W and normalized by Vs . At short times (solid red curve), v¯g increases linearly with  and connects continuously to v¯g = 0. This is the expected relation for frictionless granular systems: no motion occurs below the yield stress, while above the yield stress the strain rate increases as a power law in the difference between the applied and yield stresses [25]. However, at long times we observe a discontinuity in v¯g at c , as has also been observed in experiments [12] as well as simulations of Aeolian transport [24] and sheared frictional granular media [26]. The discontinuity in v¯g moves toward c as time increases. The size of the discontinuity and the slope v¯g / both scale roughly linearly with Rep , as shown in Fig. 3(b). For Shields numbers below c , the grains settle into a stable packing in a time ts that diverges as  → c , as shown in Fig. 3(c). Figure 1 shows that mobile grains are confined to a relatively small layer at the top of the bed. So, for fill heights studied here (N D/W > 8), we find the total flow rate v¯g N D/W as a function of  and  is insensitive to the system size (recall, v¯g is the average velocity of all grains). One might think that the discontinuity in the the flow rate near c may disappear in the large-system limit: small systems of grains are able to find a stable configuration, but very large systems

042202-3

CLARK, SHATTUCK, OUELLETTE, AND O’HERN

Θc

6

10 3

m

10 2

1

4 10 1 2

δ

0 0.3

10

0.35

Θ

0.4

10 0

0.45

10 1

10 2

Θc /Γ

2

10 3

×10

-0.9 -1.0

-1.7

-2.3

-2.8

-2.9

-3.0

4 2 0 0

0.1

0.2

0.3

Θ

0.4

0.5

0.6

0.7

5

( t s − t s 0 ) /τ s

x 10

t s /τ s

3 2

Θ0

5

10

0

15 10 5 0

5

10

0

15

F¯c /mg

5

10

F¯c /mg

15

fluid force contact force

15 10 5

0

5

10

F¯c /mg

15

3

10

0.05

fluid force contact force

0

Θc

−2.5

−2

FIG. 4. (Color online) Profiles of the layer-averaged instantaneous fluid drag force (solid red lines) and the local pressure (black dashed lines) for b = 4 and  = 0.25, with (a)  = 0.2, (b) 0.35, and (c) 0.4.

−1

10

10

Θc − Θ

1 0

5

(c)

-0.2 -0.3

4

10

0

6

(d)

15

0

4

8

ts /τs

m δ

10 0

(b)

fluid force contact force

height from bottom (D)

(c)

(a)

(b)

height from bottom (D)

8

height from bottom (D)

v¯g N D/(W Vs )

(a)

PHYSICAL REVIEW E 92, 042202 (2015)

0.1

0.15

0.2

Θ

0.25

0.3

0.35

0.4

FIG. 3. (Color online) Data from protocol A near c . (a) Time evolution of the average grain velocity v¯g (t) during protocol A. The solid red curve characterizes grain motion at the time for grains to settle under no fluid flow ts0 , the black curve with open circles marks the end of the simulation, and dashed curves represent intermediate times. At long times, we observe a discontinuity in v¯g with magnitude δ and slope m at c . (b) δ and m scale roughly linearly with Rep . (c) The time ts required for the grains to come to rest as a function of  for varying Rep (left to right, blue to red, represents large to small Rep ); symbols show  = 0.01 (∗), 0.02 (), 0.05 ( ), 0.1 (), 0.15 ( ), 0.2 (), 0.25 (), 0.5 ( ), and 1 (◦). The lines show fits of (ts − ts0 ) ∝ (c − )α , where the values of α are marked next to each plot. (d) ts is plotted as a function  for b = 4,  = 0.24, and three system sizes. The symbols correspond to (W/D, N ): black circles (100, 800), blue squares (50, 400), and red triangles (50, 800). Data points show the mean of ten simulations, and error bars give the standard deviation. The inset shows a logarithmic plot of the mean of ts − ts0 . The thick dashed line corresponds to (ts − ts0 ) ∝ (c − )−2.5 .

will always have a weak spot that does not allow the system to stop. However, this is not the case. Figure 3(d) shows that the behavior near  = c is insensitive to the system size. The stopping time as a function of  is shown for three different system sizes, and the three curves are virtually identical. The mean and fluctuations of the stopping time both diverge as a power law as  → c . The data shown is for  = 0.24 and b = 4, where c > 0 , and the power law exponent is roughly 2.5. As shown in Fig. 3(c), this exponent is roughly 3 in the low Rep limit, and it decreases with increasing Rep

to less than 1 in the high Rep limit. 0 is also shown, and it corresponds reasonably well to where the divergence begins, which suggests that c − 0 is related to the power law exponent of the divergence. In Fig. 4, we show typical force profiles during protocol A for grains under three different conditions. All panels have b = 4 and  = 0.25, but with varying . Solid red curves show the average force exerted by the fluid on the grains at a particular height. Black dashed curves show the average pressure due to grain-grain contacts as a function of height. To calculate i the pressure, we first calculate the force moment tensor Mαβ at particle i by summing over the particles j that contact particle i  ij ij i to obtain Mαβ = R1i j Fα rβ , where Ri is the particle radius, ij

α and β represent Cartesian components, and rβ represents the β component of the branch vector connecting the center of particle i with the point of contact with particle j . The mean of the eigenvalues of this tensor provides a grain-scale estimate of the local pressure, which we then average over the horizontal direction to obtain the average contact force F¯c . We plot F¯c in Fig. 4 as a function of height from the lower boundary as a black dashed line, which has roughly slope 1 in all plots. Thus, the contact forces below the surface are dominated by the weight of grains above a particular layer. Figure 4(a) shows  = 0.2, which is below both 0 and c , so grains are not moving. Figure 4(b) shows  = 0.35, so the flow is metastable and grains will eventually come to a stop. Figure 4(c) shows  = 0.4, so grains continue to move indefinitely. B. Protocol B: Static-to-mobile transition

For  > c , a sufficiently large perturbation will lead to sustained bed motion. At high Rep , this motion never

042202-4

ONSET AND CESSATION OF MOTION IN . . .

pf (yf /D)

(a)

PHYSICAL REVIEW E 92, 042202 (2015)

fluid force at the top of the bed to that inside the bed v0 /vb , as given in Table I. Thus, while large-scale particle motion always begins at the surface of the bed, this motion is often correlated to small particle rearrangements that occur deep in the bed. Figures 6(a) and 6(d) show that the probability distribution P (f ) approaches δ(f − c ) for low Rep and δ(f − 0 ) for high Rep in the large system limit, where we varied both N and W to change the system size. The insets show that the distributions P (f ) collapse when rescaled by their mean values. The system-size dependence suggests a “weakest link” picture: at a given value of excess stress above c or 0 (for low and high Rep , respectively), there is a better chance of finding a sufficiently weak local grain arrangement somewhere in a larger system. If we consider the bed to be a composite system of M uncorrelated subsystems that fails if any of the subsystems fail, the cumulative distribution CM () for failure is related to that of a single subsystem C() by

probability of failure fluid force profile

0.06 0.04 0.02 0 -35

-30

-25

-20

-15

-10

-5

0

-5

0

yf /D (b)

0.1 probability of failure fluid force profile

pf (yf /D)

0.08 0.06 0.04 0.02 0 -35

-30

-25

-20

-15

-10

1 − CM () = [1 − C()]M .

yf /D

If we assume a Weibull distribution [22,23,27]     − c α , C() = 1 − exp β

(c) 0.25 probability of failure fluid force profile

pf (yf /D)

0.2 0.15 0.1 0.05 0 -35

-30

-25

-20

-15

yf /D

-10

-5

0

FIG. 5. (Color online) The probability pf that bed failure occurs at depth yf (where yf = 0 is the top of the bed) is proportional to the local fluid force for (a) b = 2 and  = 0.1, (b) b = 4 and  = 0.25, and (c) b = 6 and  = 0.5. The bed failure depth is determined as the depth of the particle whose acceleration first exceeds athresh (or a weighted average of depths, if multiple grains meet this condition simultaneously). The agreement between pf (yf ) and the local fluid force is good for b = 4 and 6 cases, but there is a deviation for the b = 2 case, such that pf is smaller than expected near the lower boundary and bed surface. One explanation for this is that failure events are less localized, causing the weighted average of ai > athresh to be more likely in the middle of the system.

begins for  < 0 . But as the red crosses in Fig. 2 show, stable motion is not always initiated at these minimum values but depends on the particular arrangement of grains in the bed. We show that sustained grain motion is consistent with Weibullian weakest link statistics [22,23] and that failure events are always initiated at c for low Rep and 0 for high Rep for sufficiently large systems. First, we note that although the first large-scale motion of the grains always occurs at the top layer of grains, failure events do not always originate there. If we measure the depth of the particle whose acceleration first exceeds athresh at f , we find that it can occur at any depth yf below the bed surface, with a probability distribution pf that is proportional to the local applied fluid force, as shown in Fig. 5. Essentially, the ratio of the probability of failure at the surface to the probability of failure in the bed is roughly given by the ratio between the

(4)

(5)

then CM () will have the same form with αM = α and βM = βM −1/α . Figures 6(a) and 6(d) show that P (f ) = dC/df does indeed obey a Weibull distribution with shape parameter α ≈ 2.6. Thus, if Eq. (5) fits the data, and if α is constant but β ∝ M −1/α as M is varied, then Eq. (4) applies, and global failure is caused by the failure of a single member of a collection of uncorrelated subsystems. Figures 6(b) and 6(e) confirm this, showing that the mean Shields number at flow −1/α ¯ f − c ) ∝ Meff ¯ f scales as ( for small Rep and onset  −1/α ¯ f − 0 ) ∝ Meff for large Rep , where Meff = Weff Heff / ( D 2 is the effective system size. This scaling means that larger systems are more likely to fail near c or 0 for small and large Rep , respectively. However, systems that fail near these minimum values are very slow to become fully mobilized. To demonstrate this, we consider the mobilization time tm , defined as the time between the initial force imbalance that leads to failure and the time when v¯g reaches its asymptotic value. As shown in Figs. 6(c) and 6(f), this time scale also diverges as  → c for low Rep and  → 0 for high Rep , and is independent of system size. To calculate Meff , we determine Weff and Heff as follows. Since the vertical symmetry is broken by the fluid forcing profile, we calculate Heff by integrating the probability of failure over the depth of the system, which is equal to the fluid force profile. That is, Heff is the integral of the profiles shown in Fig. 5 over the depth of the bed, since the relevant system size for the Weibullian weakest link scaling has to do with the probability of local failure, and, as previously mentioned, particles near the surface cause failure with a likelihood that is greater than particles beneath the surface by a factor v0 /vb . We find that this form for Heff collapses the data for systems that are sufficiently large in the horizontal (periodic) direction. We also check the form of Heff for a case where b = 2 and  = 0.1, such that v0 /vb ≈ 1.97, and we find it to collapse the data well according to the Weibullian scaling. The form

042202-5

CLARK, SHATTUCK, OUELLETTE, AND O’HERN

(a)

PHYSICAL REVIEW E 92, 042202 (2015)

(b)

P (Θf )

10

5

10

−1/α

0.5

1

1.5

2

2.5

¯ f − Θc ) (Θf − Θc )/(Θ

1

2

(d)

10 2

10

5

0.4

0.5

(f) 0.4

1

10 0

−1/α

7

× 102

6

Θ0

5

0.2

0.1

4

1

1.5

Θ

400 200 100 50 25

-0.9

0.05 0.1 0.2 0.4

3

0 0

1

2

¯ f − Θ0 ) (Θf − Θ0 )/(Θ

0 0.2

10 -1

10 3

tm /τs

15

2

Meff

¯ f − Θ0 Θ

¯ f − Θ0 )P (Θf − Θ0 ) (Θ

Θ0 P (Θf )

10 1

(e) 20

4 6

2

1.5

Θf

-0.43

8

Θc

4

10 0

0 0.5

8

10 -1

0

× 103

10

1

0

× 103

tm /τs

Θc

(c) 0

¯ f − Θc Θ

¯ f − Θc )P (Θf − Θc ) (Θ

15

0.6

Θf

0.8

1

2

0.05 1 10

10

2

10

1

3

0.2

Meff

0.4

Θ

0.6

0.8

FIG. 6. (Color online) Onset of bed motion is governed by Weibullian weakest-link statistics. Panels (a)–(c) correspond to low Rep with  = 0.25, and (d)–(f) correspond to high Rep with  = 0.1. Panels (a) and (d) show the probability distribution of the Shields number at bed failure f for many different system sizes (W/D from 3.125 to 200, N from 25 to 1600, and W N/D from 8 to 80). The dashed vertical ¯ f − c , where  ¯ f is the mean of each distribution. lines define (a) c and (d) 0 . The insets show that P (f ) collapses when rescaled by  ¯ f − 0 for ¯ f − c for low Rep and  The solid line is a Weibull distribution with shape parameter α ≈ 2.6. Panels (b) and (e) indicate that  −1/α high Rep scale with the effective system size Meff . Panels (c) and (f) show that tm , the mobilization time after bed failure, diverges near c and 0 , respectively, independent of system size. The insets show a logarithmic plot of tm − tm,0 versus (c)  − c and (f)  − 0 . The dashed lines show tm − tm,0 ∝ ( − c )−0.43 (c) and tm − tm,0 ∝ ( − c )−0.9 (f), and the thin vertical dashed line indicates c . Symbols ) correspond to different values of N (25, 50, 100, 200, 400, 800, and 1600, respectively) with varying W/D. Each data ( point represents an average of 20 simulations.

¯ f − c )−α ζ (W/D) = 0.91(

but with a smaller value of ξ ≈ 3. The form of Weff suggests a horizontal correlation length of roughly ξ , which is larger for low Rep (ξ ≈ 17) than for high Rep (ξ ≈ 3). 1.4 1.2 1

Weff /W

of Heff confirms our observation that surface grain motion can be initiated from deep beneath the surface, as it is also proportional to the local applied shear force. The horizontal dimension is symmetric, but we find finite system-size effects when W/D < ξ , where ξ is a horizontal correlation length that varies with Rep . As shown in Fig. 6, we ¯ f − c )−α for large systems, where Heff find Meff = 0.91( is calculated using the method described previously. However, when the horizontal dimension becomes small, we observe ¯ f − c ) is larger than expected, which corresponds to that ( an effective system width Weff that is smaller than the real system width W . ¯ f − c )−α = Heff Weff /D 2 , and we asIf Meff = 0.91( sume that Weff /W = ζ (W/D), then we can write (6)

A plot of this quantity is shown in Fig. 7, and the fit line corresponds to   W . (7) ζ (W/D) = 1 − exp − ξD As Fig. 7 shows, the form Weff /W = 1 − exp [− ξWD ] captures the finite-size effects. We find a similar result for the high Rep case in Figs. 6(d)–6(f), where  = 0.1 and b = 4,

0.6 0.4 0.2 0 0

2

D . Heff W

0.8

50

100

W/D

150

200

FIG. 7. (Color online) This plot shows the finite system size effects in the horizontal periodic direction (with b = 4 and  = 0.25, which are the same values for the low Rep data shown in Figs. 6(a)– 6(c)). The horizontal axis is the system width in particle diameters, ¯f − W/D. As we show in Figs. 6(a)–6(c), we find Meff = 0.91( c )−α for large systems. The vertical axis shows this quantity, ¯ f − c )−α , divided by the product of the effective fill height 0.91( Heff /D and the width in particle diameters W/D. This shows a good collapse, and the fit line corresponds to Weff /W = 1 − exp [− ξWD ], with ξ = 16.7, and 95% confidence interval of roughly 15 < ξ < 20.

042202-6

ONSET AND CESSATION OF MOTION IN . . .

PHYSICAL REVIEW E 92, 042202 (2015)

IV. SUMMARY

In summary, we performed numerical simulations of a granular bed subjected to a simple fluid shear flow to understand general features of the initiation and cessation of grain motion. The critical Shields number for the onset of grain motion c (Rep ) from the simulations is consistent with the behavior from a large body of experimental results [7,10,11]. At low Rep , c (Rep ) separates mobile and static beds, but at high Rep , we observe significant hysteresis as a consequence of particle inertia. We find that the onset of grain motion is directly connected to the packing structure,

[1] H. A. Einstein, The Bed-load Function for Sediment Transportation in Open Channel Flows, Technical bulletin (U.S. Department of Agriculture, 1950). [2] F. Charru, B. Andreotti, and P. Claudin, Sand ripples and dunes, Annu. Rev. Fluid Mech. 45, 469 (2013). [3] P. Doron, D. Granica, and D. Barnea, Slurry flow in horizontal pipesexperimental and modeling, Int. J. Multiphase Flow 13, 535 (1987). [4] J. Capecelatro and O. Desjardins, Eulerian-Lagrangian modeling of turbulent liquid-solid slurries in horizontal pipes, Int. J. Multiphase Flow 55, 64 (2013). [5] A. Shields, Application of similarity principles and turbulence research to bed-load movement, Technical Report No. 167, California Institute of Technology, Pasadena, CA (1936). [6] P. L. Wiberg and J. D. Smith, Calculations of the critical shear stress for motion of uniform and heterogeneous sediments, Water Resour. Res. 23, 1471 (1987). [7] J. M. Buffington and D. R. Montgomery, A systematic analysis of eight decades of incipient motion studies, with special reference to gravel-bedded rivers, Water Resour. Res. 33, 1993 (1997). [8] F. Charru, H. Mouilleron, and O. Eiff, Erosion and deposition of particles on a bed sheared by a viscous flow, J. Fluid Mech. 519, 55 (2004). [9] M. Ouriemi, P. Aussillous, M. Medale, Y. Peysson, and L. Guazzelli, Determination of the critical shields number for particle erosion in laminar flow, Phys. Fluids 19, 061706 (2007). [10] A. A. Beheshti and B. Ataie-Ashtiani, Analysis of threshold and incipient conditions for sediment movement, Coastal Eng. 55, 423 (2008). [11] S. Dey and A. Papanicolaou, Sediment threshold under stream flow: A state-of-the-art review, KSCE J. Civil Eng. 12, 45 (2008). [12] E. Lajeunesse, L. Malverti, and F. Charru, Bed load transport in turbulent flow at the grain scale: Experiments and modeling, J. Geophys. Res. 115, F04001 (2010). [13] J. J. Derksen, Simulations of granular bed erosion due to laminar shear flow near the critical shields number, Phys. Fluids 23, 113303 (2011).

even deep in the bed where there is a weak but nonzero fluid stress [19–21]. Our results from this simple model clarify the essential physics governing the transition between mobile and static beds. In future work, additional effects such as turbulent flow, intergrain friction, and nontrivial particle shape can be added one-by-one to determine their distinct effects on the onset and cessation of grain motion. ACKNOWLEDGMENT

This work was supported by the U.S. Army Research Office under Grant No. W911NF-14-1-0005.

[14] A. Hong, M. Tao, and A. Kudrolli, Onset of erosion of a granular bed in a channel driven by fluid flow, Phys. Fluids 27, 013301 (2015). [15] M. Houssais, C. P. Ortiz, D. J. Durian, and D. J. Jerolmack, Onset of sediment transport is a continuous transition driven by fluid shear and granular creep, Nat. Commun. 6, 6527 (2015). [16] L. E. Silbert, D. Ertas¸, G. S. Grest, T. C. Halsey, and D. Levine, Geometry of frictionless and frictional sphere packings, Phys. Rev. E 65, 031304 (2002). [17] P.-E. Peyneau and J.-N. Roux, Frictionless bead packs have macroscopic friction, but no dilatancy, Phys. Rev. E 78, 011307 (2008). [18] C. F. Schreck, N. Xu, and C. S. O’Hern, A comparison of jamming behavior in systems composed of dimer- and ellipseshaped particles, Soft Matter 6, 2960 (2010). [19] H. E. Rose, On the resistance coefficientreynolds number relationship for fluid flow through a bed of granular material, Proc. Inst. Mech. Eng. 153, 154 (1945). [20] G. S. Beavers and D. D. Joseph, Boundary conditions at a naturally permeable wall, J. Fluid Mech. 30, 197 (1967). [21] A. E. Scheidegger, The Physics of Flow Through Porous Media (University of Toronto Press, Toronto, 1974). [22] W. Weibull, A Statistical Theory of the Strength of Materials, Ingeni¨orsvetenskapsakademiens handlingar (Generalstabens litografiska anstalts f¨orlag, 1939). [23] W. Weibull, A statistical distribution function of wide applicability, J. Appl. Mech. 18, 293 (1951). [24] M. V. Carneiro, T. P¨ahtz, and H. J. Herrmann, Jump at the Onset of Saltation, Phys. Rev. Lett. 107, 098001 (2011). [25] N. Xu and C. S. O’Hern, Measurements of the yield stress in frictionless granular systems, Phys. Rev. E 73, 061303 (2006). [26] M. Otsuki and H. Hayakawa, Critical scaling near jamming transition for frictional granular particles, Phys. Rev. E 83, 051301 (2011). [27] S. V. Franklin, Extensional rheology of entangled granular materials, Europhys. Lett. 106, 58004 (2014).

042202-7

Onset and cessation of motion in hydrodynamically sheared granular beds.

We performed molecular dynamics simulations of granular beds driven by a model hydrodynamic shear flow to elucidate general grain-scale mechanisms tha...
566B Sizes 1 Downloads 6 Views