DEM simulations of quasi-two-dimensional flow of spherical particles on a heap without sidewalls

Surface flows of granular materials find several important applications in both nature as well as industry. The effect of sidewalls on such flows is known to be large. Here, we study the rheology of such flows on a quasi two-dimensional heap without sidewalls, at different mass flow rates. It is seen that the surface angle of the heap, for all the mass flow rates, is the same and corresponds to the neutral angle. System variables such as the velocity, volume fraction and stresses are reported as a function of depth from the free surface of the heap. The friction coefficient and volume fraction are also studied as a function of the scaled local shear rate and these are also found to be independent of the mass flow rate. The behaviour observed in the present work is different from that reported in previous studies of surface flows with side walls.


Introduction
In granular mechanics, a surface flow encompasses a shallow layer of material flowing over a 'fixed bed' composed of the same material. These flows are of paramount importance for the description of several natural phenomena such as flow of sand dunes, avalanches, pyroclastic flows, etc. From an industrial viewpoint also, these special types of flow have applications in transportation, storage, and mixing of materials such as grains, powders, pellets, etc.
The erosion or deposition of material into the 'static' bed is postulated to be governed by the difference between the surface angle and neutral angle (surface angle of the heap when there is no erosion or deposition of materials ) of the heap [1,2]. However, it was later shown that the bed is actually not static. Rather, the bed region shows a 'creeping' motion with velocity profiles decaying exponentially within the heap [3]. Surface flows in channels with sidewalls have been extensively studied, and the stability of such flows has been attributed to the friction of the sidewalls [4]. In such confined systems, the surface angle has been shown to increase with increase in mass flow rate [5]. Similar flows have been studied in rotating cylinder systems, and the transition between the 'creeping' bed and 'flowing' layer has been stated to be analogous to a glass transition [6]. Development of constitutive relations to describe such flows is a challenge and the µ(I) rheology, which relates the ratio of the stresses to a rescaled shear rate, has been successful in characterising the behaviour in the flowing part of such systems [7,8].
The role of sidewalls in such flows has been investigated by varying the channel width from 20 particle diameters to about 600 particle diameters, and it has been observed that, with increase in the channel width, flows * e-mail: 184020003@iitb.ac.in * * e-mail: khakhar@iitb.ac.in become slower and thicker [9]. It has also been reported that frictional side-walls aid in order-disorder transition of granular flows of monodisperse spheres down an incline with a bumpy base, whereby both the friction coefficient and the packing fraction is altered [10]. Recently, it has also been reported that, for surface flow over an unconfined asymmetric conical heap, the surface angle is significantly lower than the angle of repose as well as independent of flow rates [11]. Hence, the presence of side-walls seems to affect the flow behaviour on heaps quite significantly.
Here, our aim is to study surface flow of granular materials in a system where the boundary effects in the spanwise direction are minimal. In the following sections, the simulation methodology and the rheological characterisation of the resulting system are discussed. The paper is organised into the following sections: section 2 covers the methodology used, significant results are highlighted in section 3 and section 4 summarises the findings.

Simulation Methodology
Flow of particles on the surface of a heap, composed of the same material as that of the flowing layer, is simulated by means of the soft-particle discrete element method (DEM), using the open source software LAMMPS (http://lammps.sandia.gov) [12]. The simulation domain (see fig. 1) consists of a rectangular box, having a rough base, formed by pouring particles of 1 mm diameter before the start of the heap formation process. Periodic boundary conditions are applied in the y direction. The length (L x ) of the simulation domain is 400d and the width (L y ) of the box is 15d where d is the diameter of the particles. Such a system mimics the behaviour of a real 3D heap. Particles having diameter d = 1 mm and density ρ p = 2.5 g/cm 3 are poured from a slit at the top left corner of the simulation box ( fig. 1).Results for four different mass flow rates viz. 57.2 g/s, 73.0 g/s, 88.3 g/s and 92.9 g/s are reported in section 3. Three different slit widths viz., 7.5d, 9.5d and 11.5d are used in order to achieve the mass flow rates of 57.2 g/s, 73.0 g/s and 88.34 g/s respectively. All components of the velocity of the particles being poured are kept zero for these cases. For the mass flow rate of 92.9 g/s, the slit width is kept fixed at 11.5 d while the z component of the velocity of the particles being poured is changed to -2 cm/s. In the simulations, the force between particles is modelled using the Hertzian model, i.e., the normal pushback force for two overlapping particles is proportional to the area of overlap of the two particles. The elastic constant for normal contact k n = 2570000 dynes/cm 2 , elastic constant for tangential contact k t = 2k n /7, viscoelastic damping constant for normal contact γ n = 5000 s −1 and the visco-elastic damping constant for tangential contact γ t =γ n /2 are used. The value for the coefficient of static friction is taken as 0.5. The values of all these parameters are the same as those for the H3 model in ref. [14]. The equations of motion are integrated by utilizing the velocity-Verlet integration scheme with a time step of ∆t = 10 −5 s. The quantities of interest (such as velocity in the flow direction, volume fraction, stresses etc.) are averaged over bins of dimension 100d × L y × 1d centred around x = 250d (shown by the green box in fig. 2). The shear rate (γ) is obtained by numerically differentiating the velocity profile.

Results and Discussion
For all four cases, the angle of inclination of the free surface was found to be 21.6 • and corresponds to the neutral angle [2], since the heap is steady and there is no erosion or deposition. A similar behaviour was observed in ref. [11], where the equilibrium neutral angle is nearly constant for almost a three-fold increase in the mass flow rate in an experimental study of surface flow over an asymmetric conical heap without side-walls. As can be seen from fig. 2, the flow in the heap without side-walls is comprised of a fluid-like zone, and a fixed bed, where the material flows relatively slowly. The bottom wall consists of a bumpy base formed by randomly fixing a layer of particles to it. The right wall is kept open so that the particles can flow out from it. Due to the absence of any resistance to flow at the exit as well as steeper surface angles, the particles accelerate near the end of the free surface. An inspection of the streamlines plotted in fig. 3 shows that the flow is a slightly converging flow, with the streamlines progressively coming closer to each other with length in the flow direction. However, in the central region (x = 20cm to x = 30cm), the flow may be approximated to be unidirectional and fully developed.  Fig. 4 gives the variation of the system variables with depth (z p ) in the flowing region for the range of mass flow rates studied. For all these plots, the origin is fixed at the free surface of the heap. The coordinates have been suitably transformed, so that z p is perpendicular to the free surface of the heap (as shown in fig. 2). We can also see that the profiles of the system parameters are qualitatively the same for all the flow rates. Hence, the inflow conditions do not affect the heap flow significantly. The velocity profiles are smooth and mostly linear in the flowing zone while showing an exponential decay deeper into the heap, as reported in previous studies [3]. The exponential decay of velocity with depth inside the heap can be clearly seen by plotting the velocity profile on a semi-log scale ( fig. 4  b). The local shear rate (γ) does not vary linearly with √ z p ( fig. 4 c), indicating that the velocity profiles do not follow a Bagnold scaling (v x ∝ z 3/2 ) [15] and hence do not follow the µ(I) model of ref. [7]. Although at first glance, the volume fraction seems to be constant (fig 4d), a closer inspection reveals that, with increasing depth inside the heap, the volume fraction actually increases slightly. The normal stress and the shear stress profiles show a linearly increasing trend with depth from the free surface into the heap. The Cauchy equations for stresses in 2-D predict the slope for a plot of normal stress (σ zz ) versus depth (z p ) to be ρ b gcosθ and that for shear stress (τ xz ) vs. z p to be ρ b gsinθ. These values have also been calculated assuming a φ value of 0.6 and are plotted as dashed lines in fig. 4(e) and (f), respectively. It can be clearly seen that the slopes of the profiles obtained from the simulations are same as those predicted from the Cauchy equations.
The ratio of the macroscopic time scale for bulk deformation to the microscopic time scale for grain rearrangements, also known as the inertial number (I), is a convenient way to characterize granular flows into either slow flows (I < 0.1) or rapid flows (I > 0.1) [20]. Here, I is calculated as: The ratio of the shear stress (τ xz ) to the normal stress (σ zz ) gives the effective friction coefficient (µ). Fig. 5 shows the variation of the effective coefficient of friction (µ) with the inertial number (I). It can be seen that, for most parts of the system, the value of µ is constant, around 0.36, with deviations occuring near the free surface of the heap shown in fig. 4 a, b, c, and d as a dotted line corresponding to I = 0.09. Data for layers about 10 particle diameters wide lie above the cut off at I = 0.09. Fig. 6 shows the spatial distribution of inertial number for different parts of the system. The span of I indicates that, for most part, the system is in a regime of slow or quasi-static flow, with I < 0.1. However, the velocities in the flowing region are quite high. Also, the data for different mass flow rates collapses to a single curve for I < 0.1, and hence,in this region, the µ(I) rheology is independent of the mass flow rate. For I > 0.03, the friction coefficient begins to decrease with increasing values of I. Such behaviour has been previously reported by [17] for I > 0.4. This decreasing trend was attributed to the rapid reduction in solid volume fraction at higher Inertial numbers and consequent transition to a kinetic regime. For flow in the quasi-static to inertial regime, taking place in a split-bottom shear cell, a modified µ(I) rheology, where, the effective friction coefficient varies as logarithmic of the inertial number, was suggested [16]. For the present system, when I < 0.1, the friction coefficient (µ) is proportional to log(I). We fit the data for I < 0.1 to the following relation: From fig. 5, the values of µ 0 and α are 0.38 and 0.0012, respectively. So, for all practical purposes, α is negligible compared to µ 0 and hence, the µ(I) rheology for the present system essentially reduces to the Coulomb Yield Condition (τ xz ≤ σ zz tan θ), for slow granular flows. Also, as shown recently in ref. [18] and [19], the granular temperature also plays a key role in determining the rheology of such systems and hence, along with µ and I, the granular temperature also needs to be considered,especially for the region where I > 0.03. Fig. 7 shows the variation of the solid volume fraction (φ) with the inertial number (I). As with the case of µ and I, for I < 0.1, there seems to be a reasonable collapse of all the data to a single curve. We fit the data to an equation of the form: From fig. 7, it can be seen that equation 3 fits the data below I < 0.1. The values of φ 0 and β are 0.598 and 0.0049, respectively.

Summary
In this work, the flow of granular materials in a quasi-two dimensional heap without sidewalls was characterised.
Results for four different mass flow rates, obtained by means of DEM simulations, were presented. For all the different mass flow rates studied, the surface angle of the heap was found to be 21.6 • , which appears to be the lowest possible angle that can support flow in the absence of side-wall friction. A similar behaviour has been noted for surface flow over an asymmetric conical heap without sidewalls, where, the equilibrium neutral angle is independent of the flow rate [11]. The velocity profiles decayed exponentially with depth inside the heap. In the flowing zone, the velocity profiles were nearly linear, as opposed to Bagnold [15] profiles. The data for friction coefficient µ as well as the volume fraction φ collapses to a single curve for all flow rates, and was found to be proportional to the logarithm of the inertial number I below I < 0.1.