Patterns in Flowing Sand: Understanding the Physics of Granular Flow ... tions and experiments has led to a solid understanding of steady ﬂow states i...

0 downloads 22 Views 2MB Size

PHYSICAL REVIEW LETTERS

PRL 103, 178302 (2009)

Patterns in Flowing Sand: Understanding the Physics of Granular Flow Tama´s Bo¨rzso¨nyi,1,2,* Robert E. Ecke,2 and Jim N. McElwaine3 1

Research Institute for Solid State Physics and Optics, Post Office Box 49, H-1525 Budapest, Hungary Condensed Matter and Thermal Physics Center for Nonlinear Studies, Los Alamos National Lab, New Mexico 87545, USA 3 Department of Applied Mathematics and Theoretical Physics, University of Cambridge, United Kingdom (Received 11 June 2009; published 23 October 2009)

2

Dense granular flows are often unstable and form inhomogeneous structures. Although significant advances have been recently made in understanding simple flows, instabilities of such flows are often not understood. We present experimental and numerical results that show the formation of longitudinal stripes that arise from instability of the uniform flowing state of granular media on a rough inclined plane. The form of the stripes depends critically on the mean density of the flow with a robust form of stripes at high density that consists of fast sliding pluglike regions (stripes) on top of highly agitated boiling material—a configuration reminiscent of the Leidenfrost effect when a droplet of liquid lifted by its vapor is hovering above a hot surface. DOI: 10.1103/PhysRevLett.103.178302

PACS numbers: 47.57.Gc, 45.70.n

Granular flows are ubiquitous in the environment and in industry, but there are still no known equations for general granular systems. For flow on an inclined plane, however, progress has been made [1–3], and current theories give a reasonable explanation of uniform flowing states. The simplicity of the inclined-plane geometry plays a crucial role in such descriptions because boundary effects, which can be exceedingly complicated in granular flows, are well defined. This system is also well suited to discrete element method (DEM) simulations [4–6] because it can be treated periodically, and steady flows often result, so good statistics can be obtained. The combination of detailed simulations and experiments has led to a solid understanding of steady flow states in this system where experiments and simulations can be well summarized [2,5,7,8] by a recently proposed model [7]. Little is known, however, about the instabilities of this model or indeed of most granular flows, and the model has only been well tested in simple shear states. In contrast, the Navier-Stokes equation of fluid flow has been known for over a century, and its accuracy has been repeatedly tested by comparing the results of linear and weakly nonlinear stability analysis to experimental systems displaying an instability from a simple state to one with a distinct pattern [9]. Such an approach, that is, investigating the growth of instabilities from their respective steady states [10,11], will certainly be very useful in testing granular constitutive models and will provide critical tests for emerging theories of granular flow. The steady and fully developed state of a rapid, dilute granular flow on a rough inclined plane was shown experimentally to be unstable to the formation of longitudinal vortices observed as lateral stripes [10]. In this pattern, the downstream velocity and the layer height vary periodically across the flow consisting of higher-slower and lowerfaster regions. The development is attributed to a mechanism analogous to the Rayleigh-Be´nard instability in heated liquid layers [12]. The average packing fraction 0031-9007=09=103(17)=178302(4)

av of this flow was below 0.3 corresponding to a relative density r ¼ av =s 0:5, where s is the static packing fraction (s is in between the random-closed-packed and random-loose-packed packing fractions [13]). This low density state, however, is hard to find either numerically or for some materials experimentally. Instead, we show that when increasing the plane inclination angle, the stripe state that emerges naturally is an instability of a dense uniform flow state, that this stripe state is robust and easy to find, and that the maxima of the downstream surface velocity correspond to the highest regions of the modulated height profile in qualitative agreement with the flow rule for the uniform state [14]. The flow was analyzed for a wide range of granular materials including different sized sand and glass beads and various copper samples with different particle shape. We demonstrate, using the apparatus illustrated in Fig. 1, that longitudinal stripes are robustly observed over a broad range of flow conditions with relative densities in the range 0:2 < r < 0:95 (corresponding to about 0:12 < av < 0:57), separated into the dense flow state for r * 0:6 (i.e., av * 0:36) and the dilute stripes for lower r as described below. In the present study, two experimental setups were used. The first apparatus, described in detail elsewhere [15], was used to characterize flow regimes over a wide range of . hopper camera 2 camera 1 laser sheet camera 3

laser

FIG. 1 (color online). Illustration of the experimental setup.

178302-1

Ó 2009 The American Physical Society

week ending 23 OCTOBER 2009

PHYSICAL REVIEW LETTERS

Because the system was enclosed in a cylindrical tube, precise measurements of the height profile hðyÞ and determination of the flow structure were performed in a second setup. It consisted of a glass plate with dimensions 2:27 m 0:4 m inclined at an angle ¼ 41:3 , see Fig. 1. In order to provide a random surface roughness, one layer of the d ¼ 0:4 mm sand grains was glued onto the glass plate or it was covered by sandpaper with similar roughness. The surface velocity of the flow, with downstream and transverse components us and vs , and the height profile hðyÞ, determined by a laser sheet, were obtained simultaneously using two cameras (Fig. 1). Using camera 3, the velocity at the bottom of the layer (ub and vb ) was taken at a location where the glass plate was clean. The flow velocities were determined using particle image velocimetry on image sequences taken at 2000 frames=s. The relative density r was measured using a method described in detail elsewhere [15]. The majority of the data presented in this Letter were obtained with sorted sand with d ¼ 0:4 0:05 mm and d ¼ 0:2 0:05 mm. The stripe state was also detected for glass beads with d ¼ 0:18 0:05 mm, d ¼ 0:36 0:05 mm, and for four sets of copper particles with similar size (d ¼ 0:16 0:03 mm) but various shapes. The angle of repose r varied between 20:9 r 33:8 for the materials tested. DEM simulations were also performed to investigate the instability. A soft particle model was used with a damped linear spring for the normal force (coefficient of restitution 0.8) and Coulomb friction for the tangential force (coefficient of friction 0.5). Particle stiffness was chosen so that the maximum overlap was less than 1%. Contact stresses were calculated according to the method of [16] using a Dirac delta function as the weight function. The time step was 1=10 of the binary collision time. The base was made of identical particles held at fixed positions taken from another simulation where a thick layer was allowed to form randomly. All quantities were nondimensionalized using the particle diameters and gravity. The instability was found over a range of parameter values: slope angle 34–39 , restitution 0.80–0.95 and width greater than 50. Below, we present results from one typical simulation with a slope angle 37. If the slope angle in the simulations is reduced below r , then the flows come to rest with a packing fraction of 0.6. We use this packing fraction to normalize the results. The system was periodic in the x direction (downslope length 24.3) and the y direction (cross-slope width 120.15). The number of particles simulated was 55 761, so the volume of the particles over the xy area corresponded to a height of 10 (height at 100% packing fraction). The system was run for several months until steady state was achieved. In the experiments, the stripe state for all the materials tested has qualitatively similar characteristics, and there is no sharp transition between the states observed in the dense and dilute flow regimes. Nevertheless, the flow structure of the dense flow state is quite different from the dilute flow case reported earlier [10]. The structure of the dense stripe

state consists of relatively narrow, dense, fast-moving regions that are also the highest. In the dilute regime, the fastmoving region corresponds to a height minimum, as schematically illustrated in Figs. 2(b) and 2(c). The continuous transition between the two regimes is characterized by increasing height of the slow-moving region as the plane inclination is increased as illustrated in Figs. 3(a)–3(c) or in movies taken for various materials at [17]. The density decreases with increasing in a similar way for all materials as shown in Fig. 2(a), where r is plotted as a function ~ ¼ H=d and normalof the normalized hopper opening H ized plane inclination tan= tanr . Generally, stripes are only observed for tan=tanr > 1:25. Stripes with the structure typical for the dense regime are observed for 0:6 < r < 0:95 (corresponding to 0:36 < av < 0:57), whereas the dilute regime, when it is observed, exists in the range for 0:2 < r < 0:7 (i.e., 0:12 < av < 0:42). In the following, we characterize the stripe structure in the dense regime using experimental data shown in Fig. 3 obtained for sand with d ¼ 0:4 mm and d ¼ 0:2 mm and numerical results shown in Fig. 4 obtained from DEM simulations. The flow pattern [17] has a downstream surface velocity us ðyÞ and a lateral surface velocity vs ðyÞ. The downstream velocity has a relatively large modulation of ðusmax usav Þ=usav 0:2, whereas the lateral velocity is very slow with maximal value vsmax 0:04usav where usav denotes the average downstream surface velocity. Cross sections of the velocity of the fully developed state in simulations are shown in Figs. 4(a) and 4(b) showing very similar downstream velocity modulation, but somewhat smaller lateral

a hom no oge flow neo stripes us f dense sim low u lat transition io zone n

stripes

dilute

wa ves

PRL 103, 178302 (2009)

gas phase

dense regime

dilute regime

c b FIG. 2 (color online). (a) The phase diagram of the system presenting the mean flow density r as a function of tan= tanr ~ ¼ H=d based on data taken for sand, glass beads, and and H various copper samples. Various flow regimes are indicated; the bullet in the dense stripe domain corresponds to the simulation data presented. Illustration of the flow structure in (b) dense and (c) dilute regimes, gray levels indicate local density.

178302-2

PHYSICAL REVIEW LETTERS

PRL 103, 178302 (2009)

a

25

d

20

u~s

tanθ/tanθr= =1.56

h

~

50 40 30

15 10 5

b

0

tanθ/tanθr=

0

u~s

e

60 50

100 y~

50

λ

=1.92

150

f

c

0

tanθ/tanθr=

20

40

λ / hav

=2.13

~ us

week ending 23 OCTOBER 2009

70 60 0

50

100

y~

150

200

4 3 2 0

80

60

~y

g

100

d = 0.2mm d = 0.4mm

5

10

15

20

25

30

~ hav

FIG. 3 (color online). Images of the flow taken in reflected light (illumination from the right) and normalized downstream pﬃﬃﬃﬃﬃﬃ surface velocity u~s ¼ us = gd as a function of the normalized transverse coordinate y~ ¼ y=d for sand with d ¼ 0:2 mm and downstream distance from the outlet x ¼ 1:55 m at plane inclinations (a) ¼ 42:6 ; (b) 48.5 and (c) 52.2 , corresponding to tan= tanr ¼ 1:56, 1.92, and 2.19, respectively. Data obtained at x ¼ 2:13 m for sand with d ¼ 0:4 mm at ¼ 41:3 : (d) height profiles hðyÞ taken at various hopper openings H, (e) laserline intensity (exposure time 4 ms), (f) transmitted light intensity, and (g) dimensionless wavelength =hav of the pattern as a function of the normalized mean flow thickness h~av ¼ hav =d. To adjust hav , the hopper opening H was varied.

velocities. In the experiments, the lateral velocity measured at the surface vs ðyÞ and at the bottom of the layer vb ðyÞ is of opposite direction in accord with the flow structure obtained from the simulations [see Fig. 4(b)]. The experimentally observed height variation in Fig. 3(d) agrees nicely with the simulation results [Figs. 4(a)–4(d)]. The fast stripes correspond to a higher narrow maximum of the hðyÞ curve, but another set of less pronounced maxima is present between them. Thus, instead of a sinusoidal hðyÞ profile observed for the dilute regime, a more complex hðyÞ is seen where a higher, global maximum corresponds to the fast flowing region and a lower, local maximum corresponds to the slower region. The double peak is also seen in the transmitted light intensity in Fig. 3(f), but the wavelength of other important measures of the pattern, e.g., the downstream velocity or velocity in the yz plane [see Figs. 3(a)–3(c), 4(a), and 4(b)], do not change during the dilute-dense transition so the emergence of a nonsinusoidal height profile does not correspond to a full frequency doubling. The pattern is observed only in a finite range of the flow thickness with the strongest amplitude at 10 < h~ < 18. The wavelength of the pattern is related to the mean flow thickness hav as 2:8 < =hav < 4:5 as shown in Fig. 3(g), so the cross section of a roll is elongated as compared to the traditional RayleighBe´nard rolls with nearly circular cross section [9]. Numerical simulations enable us to visualize spatial variations of the relative density, Fig. 4(c), and of the

FIG. 4 (color online). Results of the DEM simulations: (a) cross section of the downstream velocity u~, (b) speed in pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ﬃ ~ 2 ) with streamlines, (c) relative denthe yz plane (~ s ¼ v~2 þ w sity r ¼ =0:6 or packing fraction , (d) the inertial number I, (e) dependence of relative density on I and best quadratic fit, and (f) dependence of effective friction on the inertial number I. Solid line is best cubic fit. Dashed line is best fit to Pouliquen model (Eq. 2 in [7]).

inertial number I in Fig. 4(d). The inertial number, usually defined for incompressible flows [7], can be extended for pﬃﬃﬃﬃ pﬃﬃﬃﬃ the compressible case. We define I ¼ d jD0 j= p, where D0 is the deviatoric strain tensor [D0 ¼ D TrðDÞ=3], the density, ﬃ p the pressure, and we use the norm jAj ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ TrðAAT Þ=2. We do not consider normal pressure differences and define p ¼ TrðÞ=3, where is the stress tensor and 0 ¼ þ p. This is equivalent to the Pouliquen definition when the flow is incompressible [TrðDÞ ¼ 0]. The inertial number is proportional to the shear rate and to the ratio of the collisional stress to the total stress. The flow has the highest density in the fastmoving region where I (and shear rate) is lowest. We identify this region as a ‘‘plug’’ sliding fast on top of a ‘‘boiling’’ region with very low relative density and high inertial number and shear rate, a configuration reminiscent of the Leidenfrost effect [18] when a droplet of liquid lifted by its vapor is hovering above a hot surface. Experimental data visualizing the level of fluidization at the surface [Fig. 3(e)] and the profile of the transmitted light intensity [Fig. 3(f)] fully agree with this picture. The relative density determined from the simulation data decreases monotonically with increasing shear rate, see Fig. 4(e), and shows an amazing collapse over a wide range of densities and values of the inertial number I. This suggests that, at least for fast chute flows, a simple equation

178302-3

PRL 103, 178302 (2009)

PHYSICAL REVIEW LETTERS

of state giving the pressure as a function of shear rate and density is possible. To test the rheology, we calculate the effective friction coefficient by ¼ Trð0 D0 Þ=pjD0 j. This is the that minimizes the residual error j0 þ D0 jD 0 j j. Simpler calculations of , e.g., ¼ xz =zz or ¼ xz =p, produce poor results [no collapse would be seen in Fig. 4(f)], due to the complicated strain field. This definition is an extension of the Pouliquen model [7] to include compressible flows. Figure 4(f) shows as a function of I and demonstrates a reasonably good collapse. The data do not fit well with Pouliquen rheology and is much more strongly curved and appears to even decrease for large I [Fig. 4(f)]. At low shear rates, increases with increasing I, but at I ¼ 0:7, there appears to be a turnover above which decreases with increasing I. We believe that this behavior of the system is a key feature leading to the instability. Namely, by increasing the flow thickness above a certain value, the inertial number near the plane reaches a threshold above which the effective friction starts decreasing. As a consequence, the local inertial number (shear rate) grows even more, leading to stronger fluidization near the plane. At the same time, the fluidization drops in the upper layer (a plug develops), and this plug slides even faster on top of the expanded fluidized region. This mechanism agrees with other results for chute flows [5], but in that case the simulations were too narrow for the lateral instability to develop. On the surface, the plug absorbs material from the two sides as the surface fluidization is larger in that region, so the plug grows. In the simulations, as the instability develops, the mean flow speed decreases until a new lower velocity is reached at which point the instability has saturated. Thus, the instability may play an important role on steeper slopes, where no simple, steady state is expected, by increasing the effective viscosity. Though we have framed our discussion in terms of , it is difficult to draw too many conclusions from these results for a Pouliquen rheology. Our simulation data shows density differences, normal stress difference, and that the deviatoric stress tensor is not aligned with the strain tensor. These all disagree with the assumptions of the Pouliquen rheology indicating that a considerably more complicated rheology is necessary along with an equation of state to describe these flows. This system provides a very interesting case for studying granular rheologies because there is a complicated strain and stress field but very simple boundary conditions. Since the flow is steady, accurate measurements of all the flow variables in a simulation are possible, the only constraint being computer time. This system is therefore ideal for developing and validating granular theories in new ways. An intriguing possibility is that the lateral ridges and furrows observed in large rock avalanches [19] may be the result of the same instability. This work was funded by the US Department of Energy, Contract No. DE-AC52-06NA25396. Authors benefited from discussions with I. S. Aranson, W. B. Daniel, O.

week ending 23 OCTOBER 2009

Pouliquen, M. K. Rivera, E. Somfai, and M. van Hecke. T. B. was supported by the Bolyai Ja´nos research program and the Hungarian Scientific Research Fund Grant No. OTKA-F60157. J. N. M. was supported by the Engineering and Physical Sciences Research Council (UK).

*[email protected]i.hu [1] R. A. Bagnold, Proc. R. Soc. A 225, 49 (1954); 295, 219 (1966). [2] G. D. R. MiDi, Eur. Phys. J. E 14, 341 (2004); C. S. Campbell, Powder Technol. 162, 208 (2006); I. S. Aranson and L. S. Tsimring, Rev. Mod. Phys. 78, 641 (2006). [3] S. B. Savage, J. Fluid Mech. 92, 53 (1979). [4] R. Delannay, M. Louge, P. Richard, N. Taberlet, and A. Valance, Nature Mater. 6, 99 (2007). [5] N. Taberlet, P. Richard, J. T. Jenkins, and R. Delannay, Eur. Phys. J. E 22, 17 (2007). [6] L. E. Silbert, D. Ertas, G. S. Grest, T. C. Halsey, D. Levine, and S. J. Plimpton, Phys. Rev. E 64, 051302 (2001); D. Ertas, and T. C. Halsey, Europhys. Lett. 60, 931 (2002); L. E. Silbert, J. W. Landry, and G. S. Grest, Phys. Fluids 15, 1 (2003); D. M. Hanes and O. R. Walton, Powder Technol. 109, 133 (2000). [7] P. Jop, Y. Forterre, and O. Pouliquen, Nature (London) 441, 727 (2006); J. Fluid Mech. 541, 167 (2005). [8] I. S. Aranson, L. S. Tsimring, F. Malloggi, and E. Cle´ment, Phys. Rev. E 78, 031303 (2008). [9] M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys. 65, 851 (1993). [10] Y. Forterre and O. Pouliquen, Phys. Rev. Lett. 86, 5886 (2001); J. Fluid Mech. 467, 361 (2002). [11] S. N. Prasad, D. Pal, and J. M. Ro¨mkens, J. Fluid Mech. 413, 89 (2000); M. Y. Louge and S. C. Keast, Phys. Fluids 13, 1213 (2001); Y. Forterre and O. Pouliquen, J. Fluid Mech. 486, 21 (2003); S. L. Conway, D. J. Goldfarb, T. Shinbrot, and B. J. Glasser, Phys. Rev. Lett. 90, 074301 (2003). [12] E. M. Sparrow and R. B. Husar, J. Fluid Mech. 37, 251 (1969). [13] J. D. Bernal and J. Mason, Nature (London) 188, 910 (1960); G. Y. Onoda and E. G. Liniger, Phys. Rev. Lett. 64, 2727 (1990). [14] O. Pouliquen, Phys. Fluids 11, 542 (1999). [15] T. Bo¨rzso¨nyi and R. E. Ecke, Phys. Rev. E 74, 061301 (2006). [16] I. Goldhirsch and C. Goldenberg, Eur. Phys. J. E 9, 245 (2002). [17] See EPAPS Document No. E-PRLTAO-103-008943. Movies and images of the pattern can be downloaded from http://www.szfki.hu/~btamas/gran/stripes.html. For more information on EPAPS, see http://www.aip.org/pubservs/epaps.html. [18] J. G. Leidenfrost, De Aquae Communis Nonnullis Qualitatibus Tractatus (University of Duisburg, Duisburg, Germany, 1756); Int. J. Heat Mass Transf. 9, 1153 (1966); A.-L. Biance, C. Clanet, and D. Que´re´, Phys. Fluids 15, 1632 (2003). [19] K. Kelfoun, T. Druitt, B. van Wyk de Vries, and M.-N. Guilbaud, Bull. Volcanol. 70, 1169 (2008).

178302-4