(will be inserted by the editor)
Closure Relations for Shallow Granular Flows from Particle Simulations
Thomas Weinhart1,2· Anthony Thornton1,2,† · Stefan Luding1 · Onno Bokhove2
Received: date / Accepted: date
Abstract The Discrete Particle Method (DPM) is used to model granular flows down an inclined chute. We observe three major regimes: static piles, steady uniform flows and accelerating flows. For flows over a smooth base, other (quasi-steady) regimes are observed where the flow is either highly energetic and strongly layered in depth for small inclina-tions, or non-uniform and oscillating for larger inclinations. For steady uniform flows, depth profiles of density, ve-locity and stress have been obtained using an improved coarse-graining method, which allows accurate statistics even at the base of the flow. A shallow-layer model for granular flows is completed with macro-scale closure relations obtained from micro-scale DPM simulations of steady flows. We thus ob-tain relations for the effective basal friction, shape factor, mean density, and the normal stress anisotropy as functions of layer thickness, flow velocity and basal roughness. For collisional flows, the functional dependencies are well de-termined and have been obtained.
Keywords Discrete Particle Method · Coarse graining · Granular chute flow · Depth-averaging · Shallow-layer equations
1 General introduction 1.1 Background
Granular avalanche flows are common to natural environ-ments and industry. They occur across many orders of mag-nitude. Examples range from rock slides, containing upwards
1Dept. of Mechanical Engineering, Univ. of Twente, The Netherlands 2Dept. of Applied Mathematics, Univ. of Twente, The Netherlands †Numerical Analysis and Computational Mechanics group and
Multi-Scale Mechanics group, University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands, Tel.: +31 53 489 3301, Fax: +31 53 489 4833, E-mail: a.r.thornton@utwente.nl
of 1000 m3of material; to the flow of sinter, pellets and coke into a blast furnace for iron-ore melting; down to the flow of sand in an hour-glass. The dynamics of these flows are in-fluenced by many factors such as: polydisperity; variations in density; non-uniform shape; complex basal topography; surface contact properties; coexistence of static, steady and accelerating material; and, flow obstacles and constrictions. Discrete Particle Methods (DPMs) are an extremely pow-erful way to investigate the effects of these and other factors. With the rapid recent improvement in computational power the full simulation of the flow in a small hour glass of mil-lions of particles is now feasible. However, complete DPM simulations of large-scale geophysical mass flow will, prob-ably, never be possible.
One primary goal of the present research is to simu-late large scale and complex industrial flows using granular shallow-layer equations. In this paper we will take the first step of using the DPM [CS79, SEG+01, SGPL02, SLG03, Lud08] to simulate small granular flows of mono-dispersed spherical particles in steady flow situations. We will use a refined and novel analysis to investigate three particular as-pects of shallow chute flows: i) how to obtain meaningful macro-scale fields from the DPM simulation, ii) how to asses the flow dependence on the basal roughness, and iii) how to validate the assumptions made in depth-averaged theory.
The DPM simulations presented here will enable the con-struction of the mapping between the micro-scale and macro-scale variables and functions, thus enabling construction of a closed set of continuum equations. These mappings (closure relations) can then be used in continuum shallow-layer mod-els and compared with full DPM simulations (DPMs). For certain situations, precomputed closures should work; but, in very complicated situations the pre-established relations may fail. Heterogeneous, multi-scale modelling (HMM) is then an alternative [WEL+07], in which the local consiti-tute relations are directly used in the continuum model. In
HMM, continuum and micro-scale models are dynamically coupled with a two-way communication between the differ-ent models in selective regions in both space and time, thus reducing computational expense and allowing simulation of complex granular flows.
1.2 Shallow-layer models
Shallow granular continuum models are often used to sim-ulate geophysical mass flows, including snow avalanches [CGJ07], dense pyroclastic flows, debris flows [DI01], block and ash flows [DPP+08] and lahars [WSS08]. Such shallow-layer models involve approximations reducing the proper-ties of a huge number of individual particles to a handful of averaged quantities. Originally these models were de-rived from the general continuum incompressible mass and momentum equations, using the long-wave approximation [SH89, HSSN93, GWH99, WGH99, GTN03, BT12] for shal-low variations in the fshal-low height and basal topography. De-spite the massive reduction in degrees of freedom made, shallow-layer models tend to be surprisingly accurate, and are thus an effective tool in modelling geophysical flows. Moreover, they are now used as a geological risk assessment and hazard planning tool [DPP+08]. In addition to these geological applications, model applications involve small-scale laboratory chute flows containing obstacles [GTN03], wedges [HH05, GC07] and contractions [VATK+07], show-ing good quantitative agreement between theory and exper-iment.
In fluid dynamics, the Navier-Stokes equations are es-tablished with full constitutive equations. Nonetheless, the shallow-layer equations or Saint-Venant equations are often used in large scale situations where it is unpractical to solve the full Navier-Stokes equations. Our present aim is to di-rectly investigate the validity of the assumptions of granu-lar shallow-layer models first from discrete particle simula-tions, before obtaining fully 3D ‘kinetic theory’-style consti-tutive relations and simplifying these via the depth-integration process. A discussion of the full three-dimensional prop-erties of our particle simulations will be undertaken later. Here we restrict attention to the closures required for two-dimensional shallow granular flow equations.
A key difference between shallow-layer fluid models and granular ones is the appearance of a basal friction coeffi-cient,µ, being the ratio of the shear to normal traction at the base. In early models, a dry Coulomb-like friction law was used [SH89]. It impliesµto be constant, given by the tan-gent of the friction angle between the material and the base,
δ, i.e.,µ= tanδ. As a consequence constant uniform flow is only possible in such a model at that angleδ, independent of height. There is a considerable amount of experimental evidence, e.g., [DD99, GDR04], that suggests that such a simple Coulomb law does not hold on rough beds or for
moderate inclination angles. Furthermore, detailed experi-mental investigations using glass beads [Pou99] lead to an improved empirical ‘Pouliquen’ friction law characterised by two angles: the angle at which the material comes to rest,
δ1, where friction dominates over gravity and the angle,δ2, above which gravity dominates over friction and the material accelerates. Between these two angles steady flow is possi-ble, and in the limitδ1→δ2=δ the original Coulomb style model is recovered.
Since its formulation a lot of work has been performed on extending and understanding this Pouliquen law. The orig-inal law was obtained by retarding flowing material and mea-suring the angle at which the material stopped as a function of height hstop(θ), or equivalently, by inverting this relation,
θstop(h). For most materials, granular included, a greater angle is required to initiate stationary than to retard flow-ing material. Pouliquen and Forterre [PF02], by measurflow-ing the angle required to start motion, measuredθstart(h), i.e., the friction law for initially stationary material. As expected
θstartwas greater thanθstopand this information was used to extend the friction law to all values of the height and veloc-ity within the steady regime. Borzsonyi & Ecke performed a series of additional experiments: firstly, in [BE06] they looked at higher angles were the mean flow rates are close to the terminal velocity of a single particle, and found regions were the Pouliquen law is not valid and the Froude num-ber becomes inversely proportional to the height, as opposed to the linear relationship observed for steady flow. Borz-sonyi and Ecke, and Pouliquen and Forterre [BE07, FP03] have all worked on extending the original law to be valid for more complicated non-spherical materials like sand and metallic materials. Also, the effect of basal surface rough-ness has been systemically studied in [GTDD03] by vary-ing the size of both the free flow and fixed basal particles. For convenience, we define λ to be the size ratio of the fixed and the free particles. They observed a peak in rough-ness at a certain diameter ratio,λc, which depends on the compactness of the basal layer. Measured values of λc in [GTDD03] ranged between 1 and 3 for a monolayer of fixed particles. For fixed particles with smaller size andλ <λc, the range of angles where steady flow was observed de-creased, and eventually the steady flow regime completely vanished, i.e.,δ2−δ1→ 0 asλ→ 0 (yielding Coulomb type behaviour). For smaller flow particle diameters, i.e., with
λ >λc, there was also a reduction in friction, but weaker than in the smallλ case. For much largerλ, the friction sat-urated to a constant value, which they contributed to free particles that filled the holes in the basal surface and effec-tively created a stable basal surface of free particles. In a later publication [GDDT07], they extended this investiga-tion to flows containing two particle sizes.
Louge and Keast [LK01] modified the kinetic theory pre-sented in [Jen93] by modelling enduring contacts via a
fric-tional rate-independent stress component in order to obtain steady flow on flat frictional inclines. This work was later extended to bumpy inclines [Lou03]. Jenkins [Jen06] took a different approach and theoretically formulated a pheno-menological modification of granular kinetic theory to ac-count for enduring particle contacts. His idea is that endur-ing contacts between grains, forced by the shearendur-ing, reduce the collision rate of dissipation. Therefore a modification to the dissipation is introduced, which does not affect the stress. It leads to a law very similar to the one experimentally obtained by Pouliquen. He extended the theory in [Jen07] to very dissipative frictional particles, with a coefficient of restitution less than 0.7. Later, a detailed comparison with new experiments was performed, showing agreement with flows of low inclination [JB10].
Silbert et al. [SEG+01] used DPMs to simulated chute flow of cohesionless particles. They found that a steady-state flow regime exists over a wide range of inclination angles, heights and interaction parameters, in confirmation of the experiments of Pouliquen [Pou99]. They found for steady-state flows that the volume fraction is constant throughout the flow, in agreement with the assumptions of shallow-layer theory [SH89]. They also observed that the shear stress is proportional to the square of the shear and the flow veloc-ity scales with the height to the power 3/2. This result co-incides with Bagnold’s analysis of dilute binary collisions flows [Bag54]. They also observed small systematic devia-tions from isotropic stress, which shows a deviation from fluid-like behaviour. However, normal stresses do not ap-proach a Coulomb-yield criterion structure at the angle of re-pose except near the surface, hinting that the failure of flow starts near the surface. In [SGPL02], they looked at the effect of different basal types and found that for an ordered chute base the steady state regime splits into three distinct flow regimes: at smaller angles, the flowing system self-organizes into a state of low-dissipation flow consisting of in-plane ordering in the bulk; at higher angles, a high-dissipation regime similar to that for a rough base but with considerable slip at the bottom is observed; and, between these two sub-regions they observe a transitional flow regime characterized by large oscillations in the bulk averaged kinetic energy due to the spontaneous ordering and disordering of the system as a function of time. Finally, [SLG03] investigated the initia-tion and cessainitia-tion of granular chute flow more carefully and computed bothθstop andθstart. For inclinationsθ ≫θstop they observed a Bagnold rheology, forθ>
∼θstopa linear pro-file, and forθ≈θstopavalanching behaviour.
1.3 Overview of this study
Our present research is novel on the following three counts: Firstly, we compute more meaningful macro-scale fields from the DPM simulations than before by carefully
choos-ing the coarse grainchoos-ing function. In order to homogenize the DPM data, the micro-scale fields need to be coarse-grained to obtain macroscopic fields. Coarse-grained micro-scale fields of density, momentum and stress were derived directly from the mass and momentum balance equations by Goldhirsch [Gol10]. The quality of the statistics involved depends on the coarse graining width w, which defines the amount of spatial smoothing. For small coarse-graining width w, micro-scale variations remain visible, while for large w these smooth out in the macro-scale gradients. Since one of the objectives is to obtain the value ofµat the base, we use a novel adaptation of Goldhirsch’ statistics [Gol10] near boundaries. This is a non-trivial issue neglected in the literature; often continuum fields are simply not computed within a distance w of the boundary. Our new approach [WTLB11] is consistent with the continuum equations everywhere, enabling the construc-tion of continuum fields within one course-graining width of the boundary.
Secondly, we follow the approach of [GTDD03] and vary the basal particle diameter to achieve different basal condi-tions. For particles with smaller basal than flowing diam-eter,λ < 1, the flow becomes more energetic and oscilla-tory behaviour is observed. This phenomena has previously been reported in [SGPL02], but was achieved by changing the basal particles to a more regular, grid-like configuration. By investigating flow over fixed particles of different size than the free ones, we are able to quantify the roughness and numerically investigate the transition from rough to smooth surfaces. For smoother surfaces, we show that the parameter space can be split into to two types of steady flow, and we obtain a general friction law.
Finally, we test the assumptions made in depth-averaged theory and determine the required closure laws. For shallow granular flows, the flow can be described by depth-averaged mass and momentum-balance equations [GTN03]. Solving the depth-averaged equations requires a constitutive relation for the basal friction, a way to account for mean density variations, the shape of the velocity profile and the pressure anisotropy. We extract such data from DPMs obtained for steady uniform flows, and establish a novel, extended set of closure equations. Also, the depth-averaged equations are obtained under the assumptions that a) the density is con-stant in space and time and does not vary through the flow; b) the ratio between mean squared velocity and the squared mean velocity is known; c) the downward normal stress is lithostatic, i.e., balances the gravitational forces acting on the flow; and, d) the ratio between the normal stresses is known. Gray et al. [GTN03] assumed the latter ratio to be one. The depth profiles of these quantities are discussed by Silbert et al. in [SEG+01, SGPL02, SLG03] for steady flow. We extend these measurements to validate our DPM simula-tions, using our superior statistical procedure. Hence, we es-tablish the range in which the shallow-layer model is valid,
and the closure relations required for shallow-layer contin-uum simulations, for a much wider range of flow regimes than had been considered before.
1.4 Outline
We introduce the force model used in the DPM in Sect. 2, and the statistical method used to obtain macroscopic den-sity, velocity and stress profiles in Sect. 3. In Sect. 4, we dis-cuss the continuum shallow-layer equations for modelling granular flow including some macro-scale closures. The set up of the simulations is discussed in Sect. 5, and the steady-state regime is mapped for flows over a rough basal sur-face in Sect. 6. Depth profiles of the flow are introduced in Sect. 7, which are then used to characterize the steady flow over smoother surfaces in Sect. 8. Finally, the closure rela-tions for the shallow-layer model are established in Sect. 9, before we conclude in Sect. 10.
2 Contact law description
A Discrete Particle Method (DPM) is used to perform the simulation of a collection of identical granular particles. Bound-aries are created by special fixed particles, which generally will have different properties than the flow particles. Parti-cles interact by the standard spring-dashpot interaction model [CS79, WB86, Lud08], in which it is assumed that particles are spherical and soft, and that pairs have a single contact point.
Each particle i has a diameter di, densityρi, position ri, velocity viand angular velocity ωi. For pairs of two particles {i, j}, we define the relative distance vector ri j = ri− rj, their separation ri j= |ri j|, the unit normal ˆni j= ri j/ri jand the relative velocity vi j= vi− vj. Two particles are in con-tact if their overlap,
δn
i j= max(0, (di+ dj)/2 − ri j), (1)
is positive. A single contact point c at the centre of the over-lap is assumed, which is a valid assumption as long as the overlap is small. For our simulations the overlap between two particles is always below 1% of the particle radius; the simple contact model we employ thus captures the key pro-cess as there are no multiple contact points.
The force acting on particle i is a combination of the pairwise interaction of two particles. The force fi jrepresents the force on particle i from the interaction with particle j and can be decomposed into a normal and a tangential compo-nent,
fi j= fi jn+ fi jt. (2)
We assume that the particles are viscoelastic; therefore, they experience elastic as well as dissipative forces in both
normal and tangential directions. The normal force is mod-elled as a spring-dashpot with a linear elastic and a linear dissipative contribution,
fi jn= knδn
i jˆni j−γnvni j, (3)
with spring constant kn, damping coefficientγnand the nor-mal relative velocity component,
vni j= (vi j· ˆni j) ˆni j. (4)
For a central collision, no tangential forces are present, and the collision time tcbetween two particles can be calculated as tc=π/ s kn mi j− γn 2mi j 2 , (5)
with the reduced mass mi j= mimj/(mi+ mj). The normal restitution coefficient rc(ratio of relative normal speed after and before collision) is calculated as
rc= exp(−tcγn/(2mi j)). (6)
We also assume a linear elastic and a linear dissipative force in the tangential direction,
fi jt = −ktδδδti j−γtvti j, (7) with spring constant kt, damping coefficientγt, elastic tan-gential displacementδδδt
i j (which is explained later), and the relative tangential velocity at the contact point,
vti j= vi j− vn
i j− ωi× bi j+ ωj× bji, (8)
with bi j= − (di−δi jn)/2 ˆni j the branch vector from point
i to the contact point; for equal size particles bi j= −ri j/2. The elastic tangential force is modelled and derives from the roughness of the particle surface. Near the contact point, small bumps on a real particle would stick to each other, due to the normal force pressing them together, and elongate in the tangential direction resulting in an elastic force propor-tional to the elastic tangential displacement. It is defined to be zero at the initial time of contact, and its rate of change is given by dδδδti j dt = v t i j− (δδδt i j· vi j) ˆni j ri j , (9)
where the first term is the relative tangential velocity at the contact point, and the second term ensures thatδδδti j remains normal to ˆni j. The second term is always orthogonal to the spring direction and, hence, does not affect the rate of change of the spring length: it simply rotates it, thus keeping it tan-gential.
When the tangential to normal force ratio becomes larger than the particle contact friction coefficient, µc, for a real particle the bumps would slip against each other and their
elongation is shortened until the bumps can stick to each other again. This is modelled by a static yield criterion, trun-cating the magnitude ofδδδt
i j as necessary to satisfy|fti j| ≤
µc|f,ni j|. Thus, the contact surfaces are treated as stuck while |fi j| <t µc|fi j| and as slipping otherwise, when the yield cri-n terion is satisfied1.
The total force on particle i is a combination of contact forces fi j with other particles and external forces such as gravity g. The resulting force fiand torque qiacting on par-ticle i are fi= g + N
∑
j=1, j6=i fi j, and qi= N∑
j=1, j6=i bi j× fi j. (10)Finally, using these expressions we arrive at Newton’s equa-tions of motion for the translational and rotational degrees of freedom, mi d2ri dt2 = fi, and Ii d dtωi= qi, (11)
with mithe mass and Iithe inertia of particle i. We integrate (11) forward using Velocity-Verlet [AT93], formally second order in time, with an adequate time step of∆t= tc/50. The collision time tcis given by (5), while (9) is integrated using first-order forward Euler.
Hereafter, we distinguish between identical free flowing and identical fixed basal particles. Base particles are mod-elled as having an infinite mass and are unaffected by body forces: they do not move. This leaves two distinct types of collision: flow-flow, and flow-base. Model parameters for each of these collision types are changed independently.
3 Statistics 3.1 Coarse-graining
The main aims of this paper are to use discrete particle sim-ulations to both confirm the assumptions of and provide clo-sure rules required for the depth-averaged shallow-water equa-tions. Hence, continuum fields have to be extracted from the discrete particle data. There are many papers in the lit-erature on how to go from the discrete to the continuum, binning micro-scale fields into small volumes [IK50, SH82, Lud04, Lud09, LLV+01], averaging along planes [TED95], or coarse graining spatially and temporally [Bab97, SA04, Gol10]. Here, this will be achieved using a coarse-graining approach first described in [Bab97] based on a later deriva-tion of Goldhirsch [Gol10], extended further by us [WTLB11]
1 Meant for review stage only. It should be noted that in the absence
of dissipative forces and slipping, the system can be described as an Hamiltonian system: see Appendix A. Appendix B contains details on the tangential displacement. A pseudocode of the tangential force cal-culation is provided in Appendix C.
to account for boundary forces due to the presence of the base.
The method has several advantages over other methods because: (i) the fields produced automatically satisfy the equa-tions of continuum mechanics, even near the flow base; (ii) it is neither assumed that the particles are rigid nor spheri-cal; and, (iii) the results are even valid for single particles as no averaging over groups of particles is required. The only assumptions are that each particle pair has a single point of contact (i.e., the particle shapes are convex), the contact area can be replaced by a contact point (i.e., the particles are not too soft), and that collisions are not instantaneous.
3.2 Mass and momentum balance
3.2.1 Notation and basic ideas
Vectorial and tensorial components are denoted by Greek letters in order to distinguish them from the Latin particle indices i, j. Bold vector notation will be used when conve-nient.
Assume a system given by Npflowing particles and Nb fixed particles. Since we are interested in the flow, we will calculate macroscopic fields pertaining to the flowing parti-cles only. From statistical mechanics, the microscopic mass density of the flow,ρmic, at a point r at time t is defined by
ρmic(r,t) =
∑
Np i=1miδ(r − ri(t)) , (12)
whereδ(r) is the Dirac delta function and mi is the mass of particle i. The following definition of the macroscopic density of the flow is used
ρ(r,t) = Np
∑
i=1miW (r − ri(t)) , (13)
thus replacing the Dirac delta function in (12) by an in-tegrable ‘coarse-graining’ function W whose integral over space is unity. We will take the coarse-graining function to be a Gaussian W (r − ri(t)) = 1 (√2πw)3exp −|r − ri(t)| 2 2w2 (14)
with width or variance w. Other choices of the coarse-graining function are possible, but the Gaussian has the advantage that it produces smooth fields and the required integrals can be analyzed exactly. According to [Gol10], the answer de-pends only weakly on the choice of function, and the width w is the key parameter.
It is clear that as w→ 0 the macroscopic density defined in (14) reduces to the one in (13). The coarse-graining func-tion can also be seen as a convolufunc-tion integral between the
micro and macro definitions,
ρ(r,t) = Z
W(r − r′)ρmic(r′,t) dr′. (15)
3.2.2 Mass balance
Next we will consider how to obtain the other fields of in-terest: the momentum vector field and the stress tensor. As stated in Sect. 3.1 the macroscopic variables will be defined in a way compatible with the continuum conservation laws. The coarse grained momentum density p(r,t) is defined by pα(r,t) = Np
∑
i=1 miviαW(r − ri), (16)where the viα’s are the velocity components of particle i. The macroscopic velocity field V(r,t) is then defined as the ratio of momentum and density fields,
Vα(r,t) = pα(r,t)/ρ(r,t). (17)
It is straightforward to confirm that equations (13) and (16) lead to the continuity equation
∂ρ ∂t +
∂pα
∂rα = 0, (18)
with the Einstein summation convention for Greek letters.
3.2.3 Momentum balance
Finally, we will consider the momentum conservation equa-tion with the aim of establishing the macroscopic stress field. In general, the desired momentum balance equations are writ-ten as, ∂pα ∂t = − ∂ ∂rβ ρ VαVβ +∂σ∂αβ rβ +ρgα, (19) whereσαβ is the stress tensor, and gα a component of the acceleration vector of gravity.
Expressions (16) and (17) for the momentum p and the velocity V have already been defined. The next step is to compute their temporal and spatial derivatives, respectively, and reach closure. Taking the time derivative of (16) gives
∂pα ∂t = ∂ ∂t Np
∑
i=1 miviα(t)W (r − ri(t)) = Np∑
i=1 miv˙iαW(r − ri) + Np∑
i=1 miviα∂∂ tW(r − ri). (20)Using (11), the first term in (20) can be expressed as
Aα≡ Np
∑
i=1 miv˙iαW(r − ri) = Np∑
i=1 fiαW(r − ri). (21)In the simulations presented later the force on each par-ticle contains three contributions: parpar-ticle-parpar-ticle interac-tions, particle-base interacinterac-tions, and the gravitational body force. Hence, fiα= Np
∑
j=1, j6=i fi jα+ Nb∑
k=1 fikbα+ migα, (22)where fi j is the interaction force between particle i and j, and fikbthe interaction between particle i and base particle k, or base wall if the base is flat. Therefore, we rework (21) as
Aα= Np
∑
i=1 Np∑
j=1, j6=i fi jαWi+ Np∑
i=1 Nb∑
k=1 fikbαWi+∑
i=1 miWigα, (23)where Wi= W (r − ri). The last term in (23) can be simpli-fied toρgαby using (13). From Newton’s third law, the con-tact forces are equal and opposite, such that fi j= − fji.Hence,
Np
∑
i=1 Np∑
j=1, j6=i fi jαWi= Np∑
i=1 Np∑
j=1,i6= j fjiαWj= − Np∑
i=1 Np∑
j=1,i6= j fi jαWj. (24)where in the first step we interchanged the dummy summa-tion indices. It follows from (24) that (23) can be written as Aα = 1 2 Np
∑
i=1 Np∑
j=1, j6=i fi jα(Wi− Wj) + Np∑
i=1 Nb∑
k=1 fikbαWi+ρgα = Np∑
i=1 Np∑
j=i+1 fi jα(Wi− Wj) + Np∑
i=1 Nb∑
k=1 fikbαWi+ρgα. (25) Next, we will write Aα as the divergence of a tensor in order to obtain a formula for the stress tensor. The following identity holds for any smooth function WWj− Wi= Z 1 0 ∂ ∂sW(r − ri+ sri j) ds = ri jβ ∂ ∂rβ Z 1 0 W(r − ri+ sri j) ds, (26)
where ri j= ri− rj; we used the chain rule and differentia-tion to the full argument of W(·) per component.
The next step extends Goldhirsch’ analysis near a bound-ary. To obtain a similar expression for the interaction with base particles, we write
−Wi= Z ∞ 0 ∂ ∂sW(r − ri+ srik) ds = rikβ ∂ ∂rβ Z ∞ 0 W(r − ri+ srik) ds, (27)
which holds because Widecays towards infinity. Substitut-ing identities (26), (27) and (13) into (25) leads to
Aα= − ∂ ∂rβ Np
∑
i=1 Np∑
j=i+1 fi jαri jβ Z 1 0 W(r − ri+ sri j) ds −∂∂r β Np∑
i=1 Nb∑
k=1 fikbαrikβ Z ∞ 0 W(r − ri+ srik) ds +ρgα. (28) From [Gol10], it follows that the second term in (20) can be expressed as follows∑
i miviα∂∂ tW(r − ri) = − ∂ ∂rβ " ρVαVβ+ Np∑
i miv′iαv′iβWi # , (29) where v′i is the fluctuating velocity of particle i, with com-ponents given byv′iα(t, r) = viα(t) −Vα(r,t). (30)
Substituting (28) and (29) into momentum balance (19) yields
∂σαβ ∂rβ = ∂ ∂rβ " − Np
∑
i=1 Np∑
j=i+1 fi jαri jβ Z 1 0 W(r − ri+ sri j) ds − Np∑
i=1 Nb∑
k=1 fikbαrikβ Z ∞ 0 W(r −ri+sri j) ds − Np∑
i miv′iαv′iβWi # . (31) Therefore the stress is given byσαβ = − Np
∑
i=1 Np∑
j=i+1 fi jαri jβ Z 1 0 W(r − ri+ sri j) ds − Np∑
i=1 Nb∑
k=1 fikbαrikβ Z ∞ 0 W(r −ri+ srik) ds− Np∑
i miv′iαv′iβWi. (32) The kinetic component of the stress tensor as well as the contact stress coming from normal forces are symmetric; only the contact stress from tangential forces can contribute to the antisymmetric part of the stress tensor. In our simula-tions the tangential forces contribute less than 5% to the total stress in the system, such that the stress is almost symmetric. Equation (32) differs from the results of [Gol10] by an additional term that accounts for the stress created by the presence of the base. The contribution to the stress from the interaction of two flow particles i, j is spatially distributed along the contact line from ri to rj, while the contribution from the interaction of particles i with a fixed particle k is distributed along the line from ri to rk, extending further beyond rk. We explain the situation as follows, see Fig. 1. Stress and density profiles are calculated using (15) and (32)x ρ σ ∂xσ -2 -1 0 1 2 -2 0 2 x ρ σ ∂xσ -2 -1 0 1 2 -2 0 2
Fig. 1: Stress and density profiles are shown for two one-dimensional two-particle systems, each with two particles of unit mass at positions
x= ±1, and repelling each other (so with d > 2 for our granular case).
In the top figure, both particles are flowing, while in the bottom figure the left particle is fixed and the right one flowing.
for two one-dimensional two-particle systems, each with two particles of unit mass at positions x= ±1, repelling each other with a force| f | = 1 and with w = 0.2. In the top figure, both particles belong to the flowing species, so the density is distributed around the particles’ center of mass and the stress along the contact line. In the bottom figure, the left particle is a fixed base particle and the right particle is a free flowing one, so density is distributed around the flowing particle’s center of mass and the stress along the line extending from the flowing particle to negative infinity.
The strength of this statistical method is that the spa-tial coarse graining satisfies the mass and momentum bal-ance equations exactly at any given time, irrespective of the choice of the coarse graining function. Further details about the accuracy of the stress definition (32) are discussed in [WTLB11]. The expression for the energy is also not treated in this publication, see [Bab97].
4 Mathematical background
In this section, we briefly outline the existing knowledge to develop a continuum shallow-layer theory for granular flow.
4.1 Shallow-layer model
Shallow-layer models have been shown to be an effective tool in modelling many geophysical mass flows. Early ava-lanche models were formulated by adding gravitational ac-celeration and Coulomb basal friction to shallow-layer mod-els [GEY67, KE73]. Similar dry granular modmod-els have been derived using the long-wave approximation [SH89, HSSN93, Ive97, GWH99, WGH99] for shallow variations in the flow height and slope topography and included a Mohr-Coulomb rheology by the use of an earth pressure coefficient. The key
to these theories is to depth-integrate general three-dimensi-onal equations in the shallow direction, resulting in a system of two-dimensional equations which still retains some infor-mation about variations in thickness.
Let Oxyz be a coordinate system with the x-axis down-slope and the z-axis normal to a channel with mean down-slope
θ. For simplicity, we further consider boundaries, flows, and external forcing to be (statistically) uniform in y. The contin-uum macro-scale fields are thus independent of y, while the DPM simulations remain three-dimensional and will be pe-riodic in y. The free-surface and base location are z= s(x,t) and z= b(x), respectively. The thickness of the flow is thus
h(x,t) = s(x,t) − b(x), and the bulk density and velocity components are ρ and u= (u, v = 0, w)t, respectively, as functions of x, y, z and t.
The three-dimensional flow viewed as continuum is de-scribed by the mass and momentum balance equations (18) and (19). At the top and bottom surface, kinetic boundary conditions are satisfied: D(z −s)/Dt = 0 and D(z−b)/Dt = 0 at their respective surfaces, and with material time deriva-tive
D(·)/Dt =∂(·)/∂t+ u∂(·)/∂x+ w∂(·)/∂z
(since we assumed v= 0). Furthermore, the top surface is traction-free, while the traction at the basal surface is essen-tially Coulomb-like. We decompose the traction t= tt+ tnˆn in tangential and normal components, with normal traction tn= −ˆn · (σˆn), where ˆn is the outward normal at the fixed base. The Coulomb Ansatz implies that tt = −µ|tn|u/|u| with friction factorµ> 0. Note that µ generally can be a function of the local depth and the flow velocity. Its deter-mination is essential to find a closed system of shallow-layer equations.
We consider flows that are shallow, such that a typi-cal aspect ratio ε between flow height and length, normal and alongslope velocity, or normal and alongslope varia-tions in basal topography, is small, of order O(ε). Further-more, the typical friction factorµis small enough to satisfy
µ= O(εγ) withγ∈ (0,1). We follow the derivation of the depth-averaged swallow layer equations for granular flow by [GTN03] without assuming that the flow is incompress-ible. Instead we start the asymptotic analysis from the di-mensionless form of the mass and momentum conservation equations (18) and (19), assuming only that the density is independent of depth at leading order. Density, velocity, and stress are depth averaged as follows
¯ () =1 h Z s b () dz. (33)
In the end, one retains the normal stress ratio K= ¯σxx/ ¯σzz, the shape factorα= u2/ ¯u2, and the frictionµas unknowns. The goal is to investigate whether these unknowns can be expressed as either constants or functions of the remaining
shallow flow variables, to leading order in O(ε). The latter variables are the flow thickness h= h(x,t) and the depth-averaged velocity ¯u= ¯u(x,t). At leading order, the momen-tum equation normal to the base yields that the downward normal stress is lithostatic,σzz(z) = ¯ρg cosθ(s − z) + O(ε). Depth-averaging the remaining equations, while retaining only terms of order O(εµ), yields the dimensional depth-averaged shallow-layer equations, cf., [VATK+07, BT12],
∂( ¯ρh) ∂t⋆ + ∂ ∂x⋆( ¯ρhu) = 0, (34a) ∂ ∂t⋆(h ¯ρu¯) + ∂ ∂x⋆ h ¯ραu2+K 2gh 2ρ¯cosθ = gh ¯ρS, (34b) with S= sinθ−µ u¯ | ¯u|cosθ− ∂b ∂x⋆cosθ. (34c)
To demarcate the dimensional time and spatial scales, we have used starred coordinates. These scales differ from the ones used before in the particle dynamics and the dimen-sionless ones used later in the DPM simulations. The shallow-layer equations (34) consist of the continuity equation (34a) and the downslope momentum equation (34b). The system arises also via a straightforward control volume analysis of a column of granular material viewed as continuum from base to the free surface, using Reynold-stress averaging and a leading order closure with depth averages.
While the mean density ¯ρ can be modelled as a system variable by considering the energy balance equation, we will assume that it can be expressed as a function of height and velocity ¯ρ(h, ¯u). Thus, the closure to equations (34) is de-termined when we can find the functions ¯ρ(h, ¯u), K(h, ¯u),
α(h, ¯u), and µ(h, ¯u). In Section 9.2, we will analyze if and when DPM simulations can determine these functions.
4.2 Granular friction laws for a rough basal surface
The friction coefficient,µ, was originally [HM84] taken to be a simple Coulomb typeµ= tanδ, whereδ is a fixed fric-tion angle. Note that in steady state for a flat base with b= 0, the shallow-layer momentum equation (34b) then yields
µ= tanθ. Pure Coulomb friction implies that there is only one inclination, θ =δ, at which steady flow of constant height and flow velocity exists. That turns out to be unre-alistic. Three parameterizations forµhave been proposed in the literature.
Firstly, Forterre and Pouliquen [FP03] found steady flow in laboratory investigations for a range of inclinations con-cering flow over rough basal surfaces. They measured the thickness hstop of stationary material, left behind when a flowing layer was brought to rest, with the following fit
hstop(θ)
Ad =
tan(δ2) − tan(θ) tan(θ) − tan(δ1)
where δ1 is the minimum angle required for flow, δ2 the maximum angle at which steady uniform flow is possible, d the particle diameter, and A a characteristic dimensionless length scale over which the friction varies. Note that hstop diverges for θ=δ1 and is zero forθ=δ2. For h> hstop, steady flow exists in which the Froude number, the aspect ratio between flow speed and surface gravity-wave speed (F= ¯u/√g cosθh), is a linear function of the height, F=β h
hstop(θ)−γ, for for δ1<θ<δ2, (36) whereβ andγ are constants independent of chute inclina-tion and particle size. Provided one uses the steady state as-sumptionµ= tanθ to hold (approximately) in the dynamic case as well, one can combine it with (35) and (36) to find an improved empirical friction law
µ=µ⋆(h, F) = tan(δ1) +
tan(δ2) − tan(δ1)
βh/(Ad(F +γ)) + 1. (37)
It is a closure for µ in terms of the flow variables, which has been shown its practical value. Note that in the limit
δ1→δ2=δ the Coulomb model is recovered.
Secondly, in an earlier version [Pou99], another, expo-nential fitting was proposed for hstop, as follows
h′stop(θ) A′d = ln tan(θ) − tan(δ′ 1) tan(δ2′) − tan(δ1′) , for δ1′<θ<δ2′ (38) with the same limiting behaviour, and primes used to denote the difference in the fit. It yields the friction factor
µ=µ′(h, F) = tanδ1′+ tanδ2′− tanδ1′ e n −β′h
A′d(F+γ′) o
. (39)
Equation (35) did, however, prove to be a better fit to exper-iments and is computationally easier to evaluate.
Finally, an additional correction to (36) was made in [Jen06] to include dependence on the inclination accounting for enduring particle contacts in very dissipative frictional flows, F=βJ h hstop(θ) tan2(θ) tan2(δ 1)−γ J, (40)
for which we can use any appropriate fit for hstop. It leads subsequently to a more complicated evaluation of the fric-tion law forµ. We omit further details. We will compare our DPM simulations against these rules and fits for the rough basal surface, and extend it to smoother surfaces in the up-coming sections.
z g
y
x Fig. 2: DPM simulation for approximated height H= 17.5, inclination θ= 24◦and the diameter ratio of free and fixed particles,λ= 1, at
time t= 2000; gravity direction g as indicated. The domain is periodic
in x- and y-directions. In the z-direction, fixed particles (black) form a rough base while the surface is unconstrained. Colours indicate speed, increasing from blue via green to orange.
5 Simulation description
In this section, DPM simulations are used to simulate monodis-persed granular flows. Parameters have been nondimension-alised such that the flow particle diameter d= 1, mass m = 1 and the magnitude of gravity g= 1. The normal spring and damping constants are kn= 2 · 105 andγn= 50; thus the contact duration is tc= 0.005 and the coefficient of resti-tution isε= 0.88. The tangential spring and damping con-stants are kt= (2/7)knandγt=γn, such that the frequency of normal and tangential contact oscillation and the normal and tangential dissipation are equal. The microscopic fric-tion coefficient was taken to beµc= 1/2.
The interaction parameters are chosen as in [SEG+01] to simulate glass particles of 0.1 mm size; this corresponds to a dimensional time scale ofpd/g = 3.1 ms and dimensional velocity scale√dg= 0.031 ms−1. The above parameters are identical to the simulations of Silbert et al., except that dis-sipation in tangential direction,γt, was added to damp rota-tional degrees of freedom in arresting flow. Adding of such tangential damping removes all vibrational energy for flows otherwise arrested. The differences in the simulations with
γt =γn reported here are minor relative to the case with
γt = 0. Silbert et al. also investigated the sensitivity of the results to the particle interaction parameters tc,ε, the ra-tio kn/kt, andµc; they found that while the density of the bulk material is not sensitive to these interaction
parame-t Ek in /h Eel a i θ= 60◦ θ= 50◦ θ= 40◦ θ= 30◦ θ= 29◦ θ= 28◦ θ= 27◦ θ= 26◦ θ= 25◦ θ= 24◦ θ= 23◦ θ= 22◦ θ= 21◦ θ= 20.5◦ θ= 20◦ 0 500 1000 1500 2000 100 101 102 103 104 105
Fig. 3: Shown is the ratio of kinetic to mean elastic energy over time for H= 20, basal roughnessλ= 1, and varying chute anglesθ. Flow stops for inclinationsθ≤ 20.5◦, remains steady for 21◦≤θ≤ 29◦and accelerates forθ≥ 30◦(dashed lines).
ters, the flow velocity increased with decreasing frictionµc. Nonetheless, the qualitative behaviour of the velocity pro-files did not change.
The chute is periodic and of size 20× 10 in the x- and y-directions and has a layer of fixed particles as a base. The bottom particles are monodispersed with (nondimensional) diameter λ. Various basal roughnesses are investigated by takingλ= 0, 1/2 to 2 in turn, withλ = 0 as flat base. This bottom particle layer is obtained by performing a simula-tion on a horizontal, smooth-bottom chute. It is filled with a randomly distributed set of particles and we simulate un-til a static layer about 12 particles thick is produced. Then the layers of particles at height z∈ [9.3,11]λ are selected, remaining particles deleted, and the selected ones moved downwards by 11λ. The layer is thick enough to ensure that no flowing particles can fall through the base. Their posi-tions are fixed.
Defining a ‘filling height’ H, an integer amount of about 200 H particles are inserted into the chute. To insert a par-ticle, a random location(x, y, z) ∈ [0,20] × [0,10] × [0,I] is chosen, where I= H. If the particle at this position over-laps other particles, the insertion is rejected, and the inser-tion domain is enlarged by increasing I to I+ 0.01 to en-sure that there is enough space for all particles. The initial packing fraction is aboutρ/ρp= 0.3. Thus the particles ini-tially compact to an approximated height H, giving the chute enough kinetic energy to initialize flow. Dimensionless time is integrated between t∈ [0,2000] to allow the system to reach a steady state. A screen shot of the system in steady state is given in Fig. 2.
θ h H= 10 H= 20 H= 30 H= 40 19 20 21 22 23 24 25 26 27 28 29 30 5 10 15 20 25 30 35 40 45 50
Fig. 4: Overview of DPMs forλ= 1, with markers denoting the flow
state at t= 2000: arrested •, steady ◦, and accelerating ∗ flows. Grey dash-dotted lines mark thickness h for fixed H from H= 10, 20, 30, 40
upward. The demarcation line is fitted to hstop(θ) in equation (35)
(solid line) and h′stop(θ) in (38) (dotted line). Error bars mark intervals
establishing the demarcation line. Red crosses denote the demarcation between arrested and accelerating flow found in [SEG+01].
To ensure that the size of the periodic box does not in-fluence the result, we compared density and velocity profiles of the flow at an angleθ= 24◦and height H= 17.5 for do-main sizes Lx= 10, 20, 40, Ly= 10 and Lx= 10, 20, 40,
Ly= Lx/2, and found no significant changes.
6 Arrested, steady, and accelerating flow
From the experiments of Pouliquen [Pou99], granular flow over a rough base is known to exist for a range of heights and inclinations. DPMs by [SEG+01] also showed that steady flows arose for a range of flow heights and (depth-averaged) velocities or Froude numbers. Their simulations did, how-ever, provide relatively few data points near the boundary of arrested and steady flow to allow a more adequate fit of the stopping height. In this section, we therefore per-form numerous DPMs at heights and angles near the sep-aratrix between the steady flow regime and the regime with static piles. To study the full range of steady flow regimes, simulations were performed for inclinationsθ varying be-tween 20◦ and 60◦ and approximated or ‘filling heights’
H= 10, 20, 30 and 40. In Section 8, we will repeat (some
of) these simulations for varying base roughness.
Whether a steady flow regime has been reached, is quan-tified here by assessing the time whereafter the ratio of ki-netic energy normalized by the mean elastic potential en-ergy becomes time independent. This is shown in Fig. 3, where we plot such an energy ratio for a rough base, con-stant height, and varying chute angle. The elastic potential energy is averaged over t∈ [1000,2000] to minimize fluctu-ations after start-up, but any interval larger than 100 appears
sufficient. For chute angles at most 20.5◦the kinetic energy vanishes after a short time, thus the flow arrests; for chute angles between 21◦− 29◦, a constant value is reached, indi-cating steady flow; and, for inclinations above 29◦the en-ergy keeps increasing: thus flow steadily accelerates. If the energy ratio remained constant within t∈ [1800,2000], the flow was deemed steady, otherwise the flow was deemed to be either accelerating or stopping.
Unlike fluids, the free surface of granular flows, and thus the flow height, are not well defined. In [SEG+01], the height of the flow was estimated by H, which is equivalent to as-suming a constant packing fraction ofρ/ρp=π/6. How-ever, the exact height h= s − b of the flow varies from the approximated height H due to compaction of the flow and H is typically an overestimate. In [VATK+07], the surface of the flow was defined by the time-average of the maxi-mum vertical position of all flow particles. One could also define the free surface of the flow as the height where the density vanishes. The latter two methods, however, have the disadvantage that saltating particles can lead to slightly over-estimated flow heights.
Instead, we will define the height via the downward nor-mal stress. For steady uniform flows the downward nornor-mal stress is lithostatic, i.e., balances the gravitational weight, such that
σzz(z) = Z ∞
z ρ
(z′)g cosθdz′. (41)
This is a direct consequence of the momentum balance equa-tions. Thus,σzz(z) has to decrease monotonically; the base and free surface are the heights at whichσzz(z) reaches its maximum and minimum value, respectively. However, in or-der to avoid effects of coarse graining or single particles near the boundary, we cut off the stressσzz(z) on either boundary by defining threshold heights
z1= min{z : σzz< (1 −κ) max
z∈Rσzz} and (42a)
z2= max{z : σzz>κmax
z∈Rσzz} (42b)
withκ= 1%. We subsequently linearly extrapolate the stress profile in the interval(z1, z2) to define the base b and sur-face height s as the height at which the linear extrapoation reaches the maximum and minimum values ofσzz, respec-tively,
b= z1− κ
1− 2κ(z2− z1), s = z2+
κ
1− 2κ(z2− z1). (42c) The variable most sensitive to these height choices is ¯ρ. It shows well-defined functional behaviour for our definition of height, shown later. This is not the case if we define height by the density or the method in [VATK+07]. The threshold
κ= 1% was chosen because the results in Fig. 12 were rel-atively insensitive to the choice ofκ at or above 1%.
To determine the demarcation line hs(θ;λ) between ar-rested and steady flow with good accuracy, we performed a set of simulations with initial conditions determined by the following algorithm. Starting from an initial ‘filling height’
H= 4 and inclination θ = 21◦, the angle was increased
in steps of 1◦until eventually a flowing state was reached. From this initial flowing state, the height was increased by 2 particle heights, if the flow arrested, or else the angle de-creased by 1/2◦, assuming the curve is monotonically de-creasing. Flow was defined to be arrested when the energy ratio Ekin/hEelai = 10−5 fell below within 500 time units, otherwise the flow was classified as flowing. In simulations in which such arrested flows were continued after t= 500, a further decrease of kinetic energy was observed, thus val-idating the approach. This procedure yields intervals of the inclination angle for each height and, vice versa, height in-tervals for each angle, between which the demarcation line lies. The values presented in [SEG+01] deviate at most 0.5◦ from our observations, perhaps due to the preparation of the chute bottom, or the slightly different dissipation used. A demarcating curve between steady and arrested flow was fit-ted to equations (35) and (38) by minimizing the horizon-tal, respectively vertical, distance of the fit to these intervals, see Fig. 4. Fitting hstop(θ) yields better results than h′
stop(θ) for all roughnesses and only the fit (35) will be used here-after. Similar fits will be made in Section 8 for varying basal roughness. It leads us to a study of the depth profiles for steady state flow in the following section.
7 Statistics for uniform steady flow
To obtain detailed information about steady flows, we use the statistics defined in Sect. 3. Since the flows of interest are steady and uniform in x and y, density, velocity and stress will be averaged over x, y and t. The resulting depth pro-files will depend strongly on the coarse-graining width w, which needs to be carefully selected. Representative depth profiles for particular heights, inclinations and basal rough-nesses will also be analyzed.
Depth profiles for steady uniform flow are averaged us-ing a coarse grainus-ing width w over x∈ (0,20], y ∈ (0,10] and t∈ [2000,2000 + T]. The profile of a variableχ is thus defined as hχiT w(z) = 1 200T Z 2000+T 2000 Z 10 0 Z 20 0 χw (t, x, y, z) dx dy dt, (43) withχwin turn the macroscopic field(s) of density, momen-tum and stress, as defined in (13), (16) and (32). We average in time with time snapshots taken every tc/2 units.
To determine an appropriate time averaging interval T , we calculate the rate of change in momentum from the
den-1 1 T r 100 101 102 10−2 10−1
Fig. 5: Depth-averaged norm of the momentum rate of change, r=
Rs
b|∂t(ρu)|2dz, with∂t(ρu) determined by (44) for varying time
aver-aging interval T . Steady flow at height H= 20 and inclinationθ= 24◦
was used. Temporal fluctuations decrease inversely proportional to the length of the time averaging interval.
z− b ρ / ρp w= 0.1 w= 0.25 w= 0.5 w= 1 0 5 10 15 20 25 30 0.45 0.5 0.55 0.6 0.65 0.7
Fig. 6: Particle volume fractionρ/ρpfor H= 30,θ= 24◦, andλ= 1
for varying coarse graining widths w. Microscopic layering effects are visible for w< 0.5. The density is constant in the bulk and decreases
slightly in the basal layer.
sity, velocity and stress fields by
∂(ρu)
∂t =∇· σ −ρg− u ·∇(ρu). (44)
For steady flow, the temporal variations in mass and mo-mentum should approach zero when averaged over a long enough time interval T . This is shown in Fig. 5, where we plot the depth-averaged norm of the momentum rate of change for varying time averaging interval. For T≥ 100, the tempo-ral fluctuations decrease to less than 2% of the largest term,
¯
ρg, in the momentum equation. In the remainder, we choose T = 100 as the averaging interval.
The effect of varying coarse-graining width w is shown in Fig. 6, which shows the z-profile of particle volume frac-tionρ/ρp, whereρpis the particle density. For small w we observe strong oscillations of about 0.9 particle diameters width, particularly at the base. The microscopic oscillations are increasingly smoothed out and finally vanish as we ap-proach w= 0.5. For larger w, such as w ≥ 1, the macroscopic
z− b ←σxx σzz→ σyy→ −σxz→ 0 5 10 15 20 25 5 10 15 20 25
Fig. 7: Normal and shear stresses for H= 30,θ= 28◦, andλ= 1. Shear σxz and downward normal stress σzz are balanced by gravitational
forces, see equation (41). The other normal stresses show anisotropic behaviour. (z − b)/h u / u λ= 0,θ= 24◦ λ= 1/2,θ= 22◦ λ= 1/2,θ= 26◦ λ= 2/3,θ= 24◦ λ= 5/6,θ= 24◦ λ= 1,θ= 24◦ λ= 2,θ= 24◦ Bagnold 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 @ @ @ @ @ I λ
Fig. 8: Flow velocity profile of thick flow for H= 30, θ= 24◦, λ= 0, 2/3, 5/6, 1, 2 and for H = 30, θ= 22◦, 26◦,λ = 1/2. For a
rough base withλ≥ 5/6, we see a Bagnold velocity profile (dashed line), except near the surface. For a smooth bases withλ≤ 2/3, the profile becomes more convex. Forλ= 1/2,θ< 24◦, the flow velocity
shows layering while still observing the Bagnold profile. Forλ= 0,
a considerable slip velocity is observed. Forλ= 2, the basal shear is
small due to flow particles trapped between basal particles so that the definition of the base b(x) is rather fuzzy.
gradients at the base and surface are smoothed out, an un-wanted effect of the coarse-graining. The same behaviour is observed in the stress and velocity fields. Smoothing over the microscopic effects makes it impossible to observe mi-croscopic layering in the density, which we still wish to identify in our averaged fields. Hence, we choose w= 0.25 as the coarse-graining width, such that layering effects re-main visible along with macroscopic gradients.
θ h λ= 1/2 λ= 2/3 λ= 5/6 λ= 1 λ= 2 19 20 21 22 23 24 25 5 10 15 20 25 30 35 40 45 50
Fig. 9: Demarcation line hs(θ;λ) for varying basal roughness. Markers
denote the midpoint of the intervals around which the curve was fit-ted. Steady flow is observed at smaller inclinations for smoother bases. While the smaller angleδ1,λvaries only slightly, the larger angleδ2,λ
decreases rapidly with the smoothness. Forλ= 0, the demarcation line
is vertical atθ= 12.5◦.
The microscopic oscillations at the base indicate a strong layering effect of particles near the boundary, despite the rough bottom surface. The macroscopic density throughout the flow is almost constant in the bulk and decreases slightly at the base. An approximately constant density profile is a feature of all steady flows and is an important assumption of depth-averaging.
Non-zero stress components are plotted in Fig. 7. We observe that the stress components are nearly symmetric. Shear stresses σyx andσyz are zero since the flow veloc-ity in y-direction vanishes. For steady flow, the downward normal stressσzz(z) is lithostatic and satisfies equation (41) with a maximum error of 0.5%. Since the density is nearly constant, we obtain a linear stress profile, another assump-tion of depth-averaged theory. Applying the momentum bal-ance (19) to steady uniform flow further yields that the shear stress satisfies σxz =
R∞
z ρ(z′)g sinθdz′. Thus, the macro-scale frictionµsatisfiesµ=σxz/σzz= −gx/gz= tanθ. This relation is locally satisfied for all steady flow cases to an accuracy of|θ− tan−1(µ)| < 0.4◦. The remaining normal stress components,σxx andσyy, are not constrained by this mass balance. We thus see in Fig. 7 significant anisotropy in the amplitude of the normal stresses, in particular inσyy, showing that the confining stress is largest in the flow direc-tion, except for very small inclinations. It is always weak-est in the lateral or y–direction with fluctuations at the base that are in phase with the fluctuations of the density. Gen-erally, the anisotropy increases with higher inclinations and smoother bases, as will be analyzed in future work.
8 Transition from rough to smooth base
Next, we study the effect of smoother bases on the range of steady flows by decreasing the diameterλ of the base parti-cles, with the limiting case of a flat bottom wall forλ= 0. Such an extensive numerical study of the effects of chang-ing bottom roughness appears to be novel. To that effect, the DPM simulations from Section 6 were extended such that results for basal roughnessesλ = 0, 1/2, 2/3, 5/6, 1, and 2 can be compared. Before we show the hstop–curves for these simulations, we investigate the extent to which changes in basal roughness lead to more complex density and velocity profiles.
We summarize the density profiles seen without explic-itly showing the results. For decreasing basal roughnessλ, we observe that the microscopic oscillations and the dip in density at the base increase, while the bulk density remains constant. Forλ = 1/2 andλ= 2/3 and small inclinations, we see steady flow that is strongly layered throughout the flow. In contrast, forλ = 2 there is a low flow density in the basal region, since some of the free particles are small enough to sink a little into the base, forming a mixed layer of fixed and free particles.
Velocity profiles for H= 30 andθ= 24◦are shown for varying basal roughness in Fig. 8. Forλ = 1, we observe the Bagnold profile [Bag54] for thick collisional flows, dif-fering only at the surface. For very thin flows at H= 10 or inclinations near the arresting flow regime, the profile dif-fers strongly from the Bagnold profile and becomes linear. For smoother bases, the flow velocity increases, and the pro-file becomes more concave. Weak to stronger slip velocities are observed forλ< 2/3. Forλ= 0, thicker flows have con-stant velocity throughout the depth, almost corresponding to plug flow.
A family of demarcation curves hstop(θ;λ) between steady and arrested flow is shown in Fig. 9. The curve fits are based on
hstop(θ;λ) = Aλdtan(δ2,λ) − tan(θ) tan(θ) − tan(δ1,λ)
, δ1,λ <θ<δ2,λ, (45)
in which the dependencies onλ are explicitly denoted. The fitting parametersδ1,λ,δ2,λ, Aλ appearing in (45) are found in Table 1. Again, a fit based on the original equation (35) (or (45)) rather than Pouliquen’s early fit (38) yields the best results.
For a flat bottom, such thatλ = 0, steady flow initiates and resides at or very tightly around an inclinationθ= 12.5◦ for all heights, see Fig. 9. It is in agreement with the angle found in the laboratory experiments of [GTDD03]. Hence, for a smooth base the flow is steady only at a single inclina-tion, arrests for lower inclinations and accelerates for larger
θ h δ3,1/2→ δacc,1/2→ 19 20 21 22 23 24 25 26 27 28 29 30 5 10 15 20 25 30 35 40 45 50 z/h ρ / ρp 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 t Ek in /h Eel a i 0 500 1000 1500 2000 0.5 1 1.5
Fig. 10: Top: Overview of DPM simulations forλ= 1/2, with markers
denoting the flow state at t= 2000: arrested •, layered ×, oscillating
⋄, steady ◦, and accelerating ∗ flows. Demarcation line hstop(θ; 1/2) is
fitted according to (35). Bottom left panel: Profile of particle volume fraction of layered flow at H= 20,θ= 22◦. Bottom right panel: Ratio of kinetic over mean elastic energy for oscillating flow at H= 30,θ=
24◦. An example of oscillating flow is seen Fig. 6.
inclinations. Such behaviour is in line with laboratory obser-vations and DPM simulations, e.g. [VATK+07]. For 1/2 <
λ ≤ 2, we observe Pouliquen-style behaviour in Fig. 9. The angleδ1,λof flow initiation is nearly constant with respect to
λ. In contrast, the range of angles at which both steady and arrested flow is possible,δ2,λ−δ1,λ, is maximal forλ= 1 and decreases for smoother chutes with λ < 1, as follows from Table 1. This has been reported in [GTDD03] for lab-oratory experiments, who also observed a slight decrease of the intervalδ2,λ−δ1,λ forλ >λc≈ 2. However, theirλc was measured for basal particles fixed at the same height and depended on the compactness of the base. We observe a slight decrease ofδ2,λ forλ= 2; however, the fitting curves in Fig. 9 do mildly overlap forλ≥ 1.
We recall thatδ1,λ andδ2,λare fitting parameters for the hstop-curve (45) which does not necessarily imply, though it is expected, that the flow accelerates for angles greater than
δ2,λ. Surprisingly, while steady flow is observed exclusively forθ∈ (δ1,1,δ2,1) whenλ= 1, the range of angles associ-ated with steady flow for smoother chutes (i.e., whenλ< 1) extends to greater inclinations withθ>δ2,λ. For these latter cases,δacc,λ>δ2,λ is defined as the smallest angle at which accelerating flow is observed; the DPM simulations show
that
δacc,λ= (
25◦± 1◦ ifλ= 0,
29◦± 1◦ otherwise. (46)
Note that for theλ= 0 case, between angles 12.5◦and 25◦, the flow is steady and layered, because the friction factor is nonzero.
The above is illustrated in Fig. 10 forλ = 1/2, where one observes the two different steady state regimes. At higher angles,δ3,1/2<θ<δacc,1/2, a disordered regime similar to that for a rough base is observed. At smaller angles,δ1,1/2<
θ<δ3,1/2, the flowing system self-organizes into a state of layered flow consisting of ordering in the x–y–plane for the bulk (bottom left panel of Fig. 10), except for a small inter-mediate region,θ≈δ3,1/2, where a transitional flow regime can be found. It is characterized by large oscillations in the ratio of bulk averaged kinetic to elastic energy due to a spon-taneous ordering and disordering, or stop-and-go flow, of the system as a function of time (lower right panel). The same flow regimes have been observed in [SGPL02], where the smoother bottoms were achieved by arranging the base par-ticles in a grid-like fashion. In contrast, we always use of a fully disorded base and vary the roughness by changing the basal particle size.
9 Closure relations for the depth-averaged model The goal of this section is to close the shallow-layer equa-tions (34) by a determination of the basal friction µ, the mean density ¯ρ, the stress ratio K, and the velocity profile
α, using our DPMs. A demarcation will be made of the flow regimes in which such a determination is possible.
9.1 Frictionµof shallow-layer model
For the rough base several friction laws have been proposed, as detailed in Section 4.2. In the following, we will compare these friction laws for the rough base,λ = 1, as well as for varying basal ratiosλ.
λ δ1,λ δ2,λ Aλ βλ γλ err 0 12.25 12.25 − 1.500 −4.065 0.886 1/2 17.913 21.357 11.959 0.300 −0.236 0.229 2/3 17.217 24.302 11.088 0.215 −0.143 0.104 5/6 17.425 26.828 7.465 0.203 0.039 0.113 1 17.708 32.780 3.355 0.196 0.040 0.121 2 17.518 29.712 5.290 0.189 0.080 0.137
Table 1: Table of fitting parameters δ1,λ, δ2,λ, Aλ for the curve hstop(θ;λ) andβλ,γλ for the flow rule (47), including the variance
