Laboratory experiment and discrete-element-method simulation of granular-heap flows under vertical vibration

TL;DR

This study combines laboratory experiments and discrete-element-method simulations to analyze granular heap flows under vertical vibration, confirming the applicability of a nonlinear diffusion transport model and revealing the shear band structure's characteristics.

Contribution

It demonstrates that a nonlinear diffusion transport model accurately describes surface velocities and internal flow profiles in vibrated granular heaps, supported by both experiments and simulations.

Findings

Velocity decreases exponentially with depth from the surface.

Shear band thickness scales with local pile height.

Flow profile integration matches the nonlinear diffusion model.

Abstract

Granular flow dynamics on a vertically vibrated pile is studied by means of both laboratory experiments and numerical simulations. As already revealed, the depth-averaged velocity of a fully-fluidized granular pile under strong vibration, which is measured by a high-speed laser profiler in the experiment, can be explained by the nonlinear diffusion transport model proposed by our previous paper (Tsuji et al, Phys. Rev. Lett. 120, 128001 (2018)). In this paper, we report that a similar transport model can be applied to the relation between the surface velocity and slope in the experiment. These facts are also reproduced by particle-scale numerical simulations based on the discrete element method. In addition, using these numerical results, the velocity profile inside the fluidized pile is measured. As a result, we show that the flow velocity decreases exponentially with depth from the surface of the pile, which means that a clearly fluidized region, also known as shear band structure, is localized around the surface. However, its thickness grows proportionally with the local height of the pile, i.e., the shear band does not consist of a fluidized layer with a constant thickness. From these features, we finally demonstrate that the integration of this exponentially-decreasing velocity profile is consistent with the depth-averaged velocity predicted by the nonlinear diffusion transport model.

Tables

Table 1. Table 1: Granular media used in the experiments. A.O. and JIS represent AS ONE (supplier of materials) and Japanese Industrial Standards.

Material	$d$ (mm)	$\rho$ (g/cm3)	tan${\theta }_c$	Note
1. Alumina ball	$0.5\pm 0.1$	3.9	$0.45$	A.O.
2. Alumina ball	$1.0\pm 0.1$	3.9	$0.45$	A.O.
3. Zirconia ball	$0.5\pm 0.1$	5.9	$0.46$	A.O.
4. Rough sand	$1.0\pm 0.3$	2.6	$0.65$	JIS

Table 2. Table 2: List of Δ ​ t 0 Δ subscript 𝑡 0 \Delta t_{0} for various analysis positions r 𝑟 r .

$r$ (mm)	$5$	$15$	$25$	$35$
$\Delta t_0$ (s)	$0.12$	$0.08$	$0.04$	$0.00$

Equations

$${\bar{v}}_{\mathrm{t}}=\frac{c{v}_{\mathrm{vib}}}{{\mu }^2}\frac{|\nabla h|}{1-(|\nabla h|/\mu {)}^2},$$

$${\delta }_{ij}=\left(d_i+d_j\right)/2-\left|{\mathbf{r}}_{ij}\right|,$$

$$v_{ij}^{(\mathrm{n})}={\mathbf{v}}_{ij}\cdot {\mathbf{n}}_{ij},$$

$${\mathbf{f}}_{ij}^{(\mathrm{n})}=\Theta \left({\delta }_{ij}\right)\left(k^{(\mathrm{n})}{\delta }_{ij}-{\eta }^{(\mathrm{n})}{v}_{ij}^{(\mathrm{n})}\right){\mathbf{n}}_{ij},$$

$${\mathbf{f}}_{ij}^{(\mathrm{t})}=-k^{(\mathrm{t})}{\mathbf{u}}_{ij},$$

$${\mathbf{u}}_{ij}={\int }_{\mathrm{stick}}\left({\mathbf{v}}_{ij}^{(\mathrm{t})}-\frac{({\mathbf{u}}_{ij}\cdot {\mathbf{v}}_{ij}){\mathbf{r}}_{ij}}{{\left|{\mathbf{r}}_{ij}\right|}^2}\right)dt,$$

$${\mathbf{v}}_{ij}^{(\mathrm{t})}=\left({\mathbf{v}}_{ij}\cdot {\mathbf{t}}_{ij}\right){\mathbf{t}}_{ij}-\frac{1}{2}\left(d_i{\mathbf{\omega }}_i+d_j{\mathbf{\omega }}_j\right)\times {\mathbf{n}}_{ij},$$

$${\mathbf{T}}_i=-\sum _j\frac{d_i}{2}{\mathbf{n}}_{ij}\times {\mathbf{f}}_{ij}^{(\mathrm{t})}.$$

$$m_i\ddot{{\mathbf{r}}_i}=\sum _j{\mathbf{f}}_{ij}+m_i\mathbf{g},$$

$$I_i\dot{{\mathbf{\omega }}_i}={\mathbf{T}}_i,$$

$$e=\exp \left(-{\eta }^{(\mathrm{n})}{t}_{\mathrm{col}}/m\right),$$

$$t_{\mathrm{col}}=\frac{\pi }{\sqrt{2k^{(\mathrm{n})}/m-{{\eta }^{(\mathrm{n})}}^2/m^2}}$$

$$h(r)=\frac{1}{5}\sum _{n=r/d-2}^{r/d+2}{h}_{\mathrm{act}}(nd).$$

$${\bar{v}}_{\mathrm{t}}=\frac{1}{h}{\int }_0^h{v}_{\mathrm{t}}(z)dz,$$

$$q(r,t)=\frac{1}{rdt}{\int }_0^r|h(r^{\prime },t+dt)-h(r^{\prime },t)|r^{\prime }d{r}^{\prime }.$$

$$v_{t0}=\frac{c_0{v}_{\mathrm{vib}}}{{\mu }^2}\frac{|\nabla h|}{1-(|\nabla h|/\mu {)}^2},$$

$$c_0\approx 2c,$$

$$v_{\mathrm{t}}=v_{t0}\exp \left(-\frac{z^{\prime }}{h_{\mathrm{s}}}\right),$$

$$\frac{h_{\mathrm{s}}}{h}=a|\nabla h|,$$

$${\bar{v}}_{\mathrm{t}}$$

$$f(|\nabla h|)=a|\nabla h|\left\{1-\exp \left(-\frac{1}{a|\nabla h|}\right)\right\}.$$

$$f(ϵ)=a(\mu -ϵ)\left\{1-\exp \left(-\frac{1}{a(\mu -ϵ)}\right)\right\}.$$

$$f(ϵ)$$

$${\bar{v}}_{\mathrm{t}}\sim \frac{0.5c_0{v}_{\mathrm{vib}}}{{\mu }^2}\frac{|\nabla h|}{1-(|\nabla h|/\mu {)}^2}.$$

$$h^{\prime }(r)=h(r)-g(r).$$

$$f_{cc}(\delta r)=\frac{1}{H^{\prime }(t)}{\int }_{r_s}^{r_e}{h}^{\prime }(r,t)h^{\prime }(r+\delta r,t+\Delta t)dr,$$

$$H^{\prime }(t)={\int }_{r_s}^{r_e}{h}^{\prime }(r,t{)}^{2}dr.$$

$$v_{t0}=\sum _{i=1}^6\left(\frac{\Delta r_i}{\Delta t_i}\cdot \frac{f_{cc}(\Delta r_i{)}^2}{f_{cc}^{tot}}\right),$$

Full text

I Introduction

Granular avalanches present peculiar liquid-like behavior, which has led to active discussions so far. A dense surface flow, which is also called a heap flow, is one of the avalanche types that we frequently encounter in daily life. The steady state of heap flows can be observed mainly in two situations: when granular media are continuously supplied onto the top of a static pile Lemieux and Durian (2000); Komatsu et al. (2001); GDR MiDi (2004); Jop et al. (2005, 2006); Katsuragi et al. (2010) or filled in a rotating drum GDR MiDi (2004); GRAY (2001); Bonamy (2002); Courrech du Pont et al. (2005); Yang et al. (2008); Kleinhans et al. (2011); Amon et al. (2013); Swisher and Utter (2014). In most of the experimental observations, the typical fluidized layer consists of $10^{0}-10^{1}$ grain diameters, which is spontaneously determined by the system itself. The profile of the flow velocity shows an exponential decay with depth from the surface layer Lemieux and Durian (2000); Courrech du Pont et al. (2005); Katsuragi et al. (2010). Heap flows are not in general produced unless slopes exceed more or less the angle of repose $\theta_{\rm c}$. This is however not the case when a system is subjected to perturbation.

Mechanical vibration is, above all, frequently used as a way to cause failure of frictional contacts among grains, which changes the rheological properties of granular media in a dramatic way. In a vibrating system, we can observe various phenomena peculiar to granular matter, such as convection Yamada and Katsuragi (2014), segregation Breu et al. (2003), compaction Iikawa et al. (2015), and friction weakening Caballero-Robledo and Clément (2009). Of course, heap flows are also induced by vibration. It is a well known fact that granular heap structure shows relaxation to a horizontally-flat surface when subjected to strong vibration, even if the slope is clearly less than $\theta_{\rm c}$ Jaeger et al. (1989); Sánchez et al. (2007); Swisher and Utter (2014); Khefif et al. (2018); Gaudel and Kiesgen de Richter (2018). This type of relaxation is ubiquitous in many natural phenomena. For instance, it has been claimed in planetary science that the granular-heap relaxation governs the terrain development on astronomical objects covered with granular beds called regolith (e.g., asteroid Itokawa) Richardson Jr. et al. (2005); Michel et al. (2009). Since the stability of a granular bed is dominated by a competition between gravity and vibration in many cases, a dimensionless parameter $\Gamma=A(2\pi f)^{2}/g$ Evesque and Rajchenbach (1989) is often used to characterize the occurrence condition of the above phenomena, where $A$ and $f$ are the amplitude and frequency of imposed vibration and $g$ is gravitational acceleration.

Our previous study Tsuji et al. (2018) has proposed the nonlinear diffusion transport (NDT) model in order to describe the behavior of granular particles driven by vertical vibration. The NDT model is derived on a basis of simple laboratory experiments and the energy equipartition model proposed by Roering et al. Roering et al. (1999). In Ref. Tsuji et al. (2018), we agitated a granular pile with relatively strong vertical vibration ($\Gamma\geq 2$) and measured the heap flow property under a vibrating system. By assuming the uniform fluidization of the entire pile (cf. Sánchez et al. (2007)), the NDT model for the depth-averaged velocity of a vibro-fluidized granular bed $\bar{v}_{\rm t}$ is given by

[TABLE]

where $|\nabla h|$, $\mu$, and $v_{\rm{vib}}$ are the slope of a pile, bulk friction coefficient, and maximum vibration velocity $A(2\pi f)$. This model also introduces a parameter $c$, which indicates the conversion efficiency from vertical vibration energy into horizontal granular transport energy. Furthermore, it turned out that this conversion efficiency $c$ is roughly constant independent of any experimental condition, implying the existence of universality.

Although the bulk flow property can be described from a macroscopic point of view in Ref. Tsuji et al. (2018), the NDT model cannot predict particle-scale behavior, such as velocity profile inside a fluidized pile. This limitation is also closely related to the following questions: Why is the depth-averaged velocity (Eq. (1)) determined by only the slope $|\nabla h|$ when $\mu$ and $v_{\rm vib}$ are fixed? On the analogy of conventional granular flows down an inclined plate (e.g., Pouliquen (1999); Silbert et al. (2001); Andreotti et al. (2013); Gray and Edwards (2014); Gaudel et al. (2016)), it might be more natural that the depth-averaged velocity also depends on the height of the pile. Besides, is the condition assumed in the derivation of the NDT model true that the entire granular bed is uniformly fluidized under vertical vibration? In a usual gravity-driven flow on a pile, actively flowing region is limited in the vicinity of the surface Lemieux and Durian (2000); Komatsu et al. (2001); Katsuragi et al. (2010); Courrech du Pont et al. (2005). There are also some experimental reports that a creep (very slowly moving) region exists even in a dense granular system subjected to vibration with $\Gamma\geq 2$ (e.g., Yamada and Katsuragi (2014)).

This study aims to develop a better understanding of heap flow dynamics by the mutually complementary analysis of laboratory experiments and particle-scale simulations. In the former analysis, after reviewing the work of Ref. Tsuji et al. (2018) briefly, the surface flow property is analyzed by the pattern matching of profiles using the same data, and then compared to the depth-averaged flow property. In the latter analysis, numerical simulations are utilized to obtain the information that cannot be accessed from experimental data. Particularly, the velocity profile inside the pile is extensively studied from a microscopic point of view in contrast to Ref. Tsuji et al. (2018). Finally, in order to prove the consistency between the NDT model (Eq. (1)) and particle-scale dynamics investigated in this paper, the depth-averaged velocity is deduced by integrating the internal velocity profile.

II Setup

This section firstly introduces the experimental configuration and the process of data acquisition. Then, we explain the setup of the numerical simulation, which attempts to reproduce the experiment, and the analysis method of those numerical data.

II.1 Laboratory experiment

The experimental method and data used in this study are the same as Ref. Tsuji et al. (2018). A schematic illustration of the experimental system is shown in Fig. 1(a). A conical granular pile with the angle of repose is created on a disk with radius $R=40~{}{\rm mm}$, which is horizontally mounted onto an electromechanical vibrator (EMIC 513-B/A). The granular materials used in the experiment are listed in Table 1. The surface of the disk is pasted with the same type of grains so that the heap structure can be sustained. After the creation of a pile, sinusoidal vertical vibration is continuously applied to the disk. Here, the amplitude is gradually increased during the initial short term $T_{\rm a}=0.5~{}{\rm s}$ to calmly reach stable vibration conditions without burst signals. The time when the vibration achieves a stable state is defined as $t=0~{}{\rm s}$.

The amplitude $A$ and frequency $f$ are varied in the range of $10^{-4}-10^{0}~{}{\rm mm}$ and $50-500~{}{\rm Hz}$ such that granular layers are fluidized. Actually, the onset criterion of fluidization occurrence is difficult to determine. Although $\Gamma$ seems to be the most relevant, the critical condition $\Gamma_{\rm c}$ fluctuates around $1<\Gamma_{\rm c}<2$ depending on $f$ Tennakoon and Behringer (1998); King et al. (2000). To avoid such complexity and focus only on clearly fluidized regimes, in this study a pile is subjected to relatively strong vibration of $\Gamma=2-10$.

Once granular media begin to flow to the outside of the disk, the shape of the pile gets relaxed (Fig. 1(b)). Outflowed grains are captured by an acrylic container surrounding the disk. The advantage of this experimental configuration is that the sidewall effect, which must be taken into consideration in the case of quasi-2D flows Jop et al. (2005, 2006), does not appear at all.

In order to measure the flow properties during the vibration, surface profiles of the relaxing pile are recorded by a high-speed laser profiler (KEYENCE LJ-V7080). Figure 1(c) shows profiles taken in the experimental condition of $A=0.04~{}{\rm mm}$ and $f=200~{}{\rm Hz}$ $(\Gamma=6)$ with Material 1 in Table 1. The sampling rate is 50 Hz, and the horizontal spacial resolution is 50 $\mu$m/pix. The number of pixels is $800$, which can just cover all of profiles at $r=0-40{\rm~{}mm}$, where $r$ is the radial distance from the center of the disk. The origin of the height coordinate $z=0$ is calibrated to the surface of the disk. The measurement error of the laser along the $z$ direction is less than $50~{}\mu{\rm m}$. These high resolution and accuracy of the laser measurement enable us to precisely capture the grain-scale movement as shown in the magnified view of profiles (Fig. 1(d)), where we can see that the surface layer moves almost keeping its profile pattern. Experiments were performed three times for each set of conditions in order to check the reproducibility. Henceforth, the experimental conditions shown in Figs. 1(b)$-$(d) are used for subsequent plots unless otherwise noted.

II.2 Numerical simulation

II.2.1 Contact model

Particle-scale simulations are conducted in both two dimensions (2D) and three dimensions (3D) by means of the discrete element method (DEM) (Cundall and Strack, 1979). Modeled in this study are $N$ polydisperse disks/spheres of constant material density with the maximum diameter $d$ and mass $m$. A small polydispersity equally ranging from $d/\sqrt{2}$ to $d$ is given in order to prevent the crystallization of a system (Iikawa et al., 2016). The position, velocity, and angular velocity of $i$th particle are denoted by $\bm{r}_{i}$, $\bm{v}_{i}$, and $\bm{\omega}_{i}$, respectively. The total force applied on $i$th particle is determined by a combination of gravity and the contact force with $j$th particle $\bm{f}_{ij}$, which consists of the normal part $\bm{f}_{ij}^{\rm(n)}$ and tangential part $\bm{f}_{ij}^{\rm(t)}$, i.e., $\bm{f}_{ij}=\bm{f}_{ij}^{\rm(n)}+\bm{f}_{ij}^{\rm(t)}$.

For a pair of two contacting $i$th and $j$th particles with diameters $d_{i}$ and $d_{j}$, the normal compression $\delta_{ij}$ and relative velocity ${v}_{ij}^{\rm(n)}$ are given by

[TABLE]

[TABLE]

where $\bm{r}_{ij}=\bm{r}_{i}-\bm{r}_{j}$, $\bm{v}_{ij}=\bm{v}_{i}-\bm{v}_{j}$, $\bm{n}_{ij}=\bm{r}_{ij}/\left|\bm{r}_{ij}\right|$. The normal contact force $\bm{f}_{ij}^{\rm(n)}$ is modeled using a normal spring constant $k^{\rm(n)}$ and viscous constant $\eta^{\rm(n)}$ as

[TABLE]

where $\Theta(x)$ is the Heaviside step function, i.e., $\Theta(x)=1$ for $x\geqslant 0$ and $\Theta(x)=0$ otherwise.

On the other hand, the tangential contact force $\bm{f}_{ij}^{\rm(t)}$ is modeled with a tangential spring constant $k^{\rm(t)}$ as

[TABLE]

where the tangential displacement $\bm{u}_{ij}$ is obtained by integrating the tangential relative velocities $\bm{v}_{ij}^{\rm(t)}$, which are expressed as

[TABLE]

[TABLE]

where $\bm{t}_{ij}=\left(-r_{ij,y}/\left|\bm{r}_{ij}\right|,r_{ij,x}/\left|\bm{r}_{ij}\right|\right)$, and the second term in the integral of Eq. (6) insures that $\bm{u}_{ij}$ always lies on the tangent plane of the contact point. In Eq. (6), “stick” means that the integral is performed only when $\left|\bm{f}_{ij}^{\rm(t)}\right|<\mu_{m}\left|\bm{f}_{ij}^{\rm(n)}\right|$ is satisfied, where $\mu_{m}$ is a microscopic friction coefficient. This condition indicates that the Coulomb friction criterion holds in quasistatic motion: $\left|\bm{f}_{ij}^{\rm(t)}\right|=k^{\rm(t)}|\bm{u}_{ij}|$ in the “stick” region of $\left|\bm{u}_{ij}\right|<\mu_{m}\left|\bm{f}_{ij}^{\rm(n)}\right|/k^{\rm(t)}$, while $\left|\bm{f}_{ij}^{\rm(t)}\right|$ remains $\mu_{m}\left|\bm{f}_{ij}^{\rm(n)}\right|$ in the “slip” region of $\left|\bm{u}_{ij}\right|\geq\mu_{m}\left|\bm{f}_{ij}^{\rm(n)}\right|/k^{\rm(t)}$. In addition, the torque $\bm{T_{i}}$ of $i$th particle is given by

[TABLE]

Using the contact forces introduced above, the translational and rotational accelerations of $i$th particle are determined by Newton’s second law:

[TABLE]

[TABLE]

where $m_{i}$ and $I_{i}$ are the mass and moment of inertia of $i$th grain, and $\bm{g}=(0,0,-g)$ is gravity. In our simulation, $d$, $m$, and $g$ are set to be unity, and all of the quantities are computed in dimensionless forms, where the unit time is $\sqrt{d/g}$. After the computation, we give dimensions to all of quantities using the same units as the experiment $(d=0.5~{}{\rm mm},m=~{}2.5\times 10^{-4}~{}{\rm g},~{}{\rm and}~{}g=9800~{}{\rm mm/s^{2}})$.

Another important point in DEM simulations is how to choose mechanical parameters. This study adopts sufficiently hard spheres/disks ($k^{\rm(n)}=10^{4}~{}[mg/d]$) so that this choice does not have an influence on the simulation result. In fact, the simulation result changes mostly only within error even if $k^{\rm(n)}=10^{5}~{}{\rm or}~{}10^{6}~{}[mg/d]$ is used. The ratio $k^{\rm(t)}/k^{\rm(n)}$ depends on the material property of particles, which is typically set to be $k^{\rm(t)}/k^{\rm(n)}=2/7-2/3$ in DEM simulations. However, we found that the result is not sensitive to this ratio in this study. Although all of the results shown below are obtained with $k^{\rm(t)}/k^{\rm(n)}=2/7$, the data using $k^{\rm(t)}/k^{\rm(n)}=2/5~{}{\rm or}~{}2/3$ change only within error as well. In contrast, a microscopic friction coefficient $\mu_{m}$ has a little influence on the result. In this study, relatively large friction $\mu_{m}=0.8$ is used, the reason of which is explained in Sec. II.2.2.

The energy loss due to inelastic collisions is characterized by the coefficient of restitution $e$, which is defined as the ratio of the post-collisional to pre-collisional normal relative velocity and also can be analytically written as

[TABLE]

where $t_{\rm col}$ is the collision time:

[TABLE]

To investigate the dependence on the degree of inelasticity, we vary $e$ from $0.7$ to $0.9$ by arranging $\eta^{\rm(n)}$, which corresponds to typical granular matter.

Last but not least, care must be taken for the time step of the calculation $\delta t$. In our simulation, the bulk flow property is not sensitive to $\delta t$ once it becomes less than $t_{\rm col}/20$. Therefore, all of the simulations in this study are conducted with the time step $\delta t<t_{\rm col}/20$.

II.2.2 Simulation procedure

Here, the simulation procedure from the creation of a pile to vibration is explained. Firstly, particles are randomly filled into a triangle (2D) or cone (3D) space on a fixed base. The base is composed of the same grains with diameter $d$, which are placed without gaps, and the system size is, unless otherwise noted, the same as the experiment ($R=40$ mm). The initial packing fraction $\phi_{\rm init}$ is set to be slightly smaller than the jamming point Otsuki and Hayakawa (2009): $\phi_{\rm init}=0.75~{}{\rm(2D)}~{}{\rm and}~{}0.55~{}{\rm(3D)}$. The slope of a filled triangle/cone is always larger than the angle of repose $\theta_{\rm c}$. The number of simulated particles is $\sim 5\times 10^{3}~{}{\rm(2D)}~{}{\rm and}~{}\sim 5\times 10^{5}~{}{\rm(3D)}$.

Next, by imposing gravity to all of particles simultaneously, a pile with the angle of repose is spontaneously created. Since simulated particles are completely spherical, the angle of repose $\theta_{\rm c}$ produced in DEM simulations is smaller than that of real grains. In addition, $\theta_{\rm c}$ slightly depends on a microscopic friction coefficient $\mu_{m}$ Zhou et al. (2001, 2002). In this study, $\theta_{\rm c}$ increases with $\mu_{m}$, and almost saturates in the range of $\mu_{m}\geq 0.8$. We have confirmed that the flow property does not change once $\mu_{m}$ exceeds $0.8$ as the bulk frictional property does not change either. In order to reproduce as realistic a pile as possible, $\mu_{m}=0.8$ is employed in this study, where ${\tan\theta_{\rm c}}=0.35~{}({\rm 2D})~{}{\rm and}~{}0.40~{}({\rm 3D})$.

Finally, particles at the base are vertically vibrated with $A\sin(2\pi f)t$ in the same way as the experiment. Note that we use a smaller amplification period ($T_{\rm a}=0.25$ s) than laboratory experiments so that the flow property can be measured in broad $|\nabla h|$ ranges. The time evolution of a simulated vibrating pile is shown in Fig. 2 for both 2D and 3D cases. The vibration condition is same as Fig. 1 and grains with $e=0.8$ are used, which are also used for subsequent plots unless otherwise specified. More than $10$ and $5$ simulation runs for each set of conditions were conducted for the 2D and 3D cases, respectively.

II.2.3 Production of surface profiles

From the obtained DEM data, the surface profile is produced to confirm the consistency with the experiment. The resolution along the radial direction is set to be $1d$, in which the position of the highest-located particle is given as the actual height $h_{\rm act}$. However, this raw profile scatters due to the saltation motion of particles, which could lead to a large error on the result. Therefore, to stabilize the data, we take the moving-average with a $5d$ window in the following way:

[TABLE]

As an example, the surface profiles produced by Eq. (13) are drawn onto 2D piles in Fig. 2(a), which are in good agreement with surface particles.

III Experimental result

III.1 Depth-averaged velocity

First, the consistency between the NDT model (Eq. (1)) and experimental data, which has been confirmed by Ref. Tsuji et al. (2018), is briefly reviewed. According to the NDT model, the depth-averaged velocity along the horizontal direction $\bar{v}_{\rm t}$, defined as

[TABLE]

where $v_{\rm t}(z)$ is flow velocity at position $z$, is a function of only slope $|\nabla h|$ when the other experimental conditions such as $\mu$ and $v_{\rm vib}$ are fixed. To check this, $\bar{v}_{\rm t}$ and $|\nabla h|$ are measured at four different points $(r=5,15,25,35~{}{\rm mm})$ and various time ($t=0.1\times\sqrt{2}^{n}~{}{\rm s}~{}(n=0,1,2,\cdots)$). $\bar{v}_{\rm t}$ can be computed by the volume flux $q$ divided by the height $h$ for each position $r$ and time $t$, i.e., $\bar{v}_{\rm t}(r,t)=q(r,t)/h(r,t)$, where $q(r,t)$ is defined as the volume of granular media flowing across a unit length per unit time and can be calculated from the volume change of the relaxing pile as

[TABLE]

$|\nabla h|$ is locally measured by the linear least-squares method using the neighbor profiles of $r\pm 5~{}{\rm mm}$. The detailed measurement method of these quantities is explained in Ref. Tsuji et al. (2018) and its supplementary material.

In Fig. 3(a), $\bar{v}_{\rm t}$ is plotted against $|\nabla h|$, where all of the data are collapsed into a single curve and can be fitted by Eq. (1). As shown in the inset of Fig. 3(a), this scaling can also be observed for experimental data obtained in different experimental conditions. Since Ref. Tsuji et al. (2018) assumes $\mu=\tan\theta_{\rm c}$ for the sake of simplicity, the fitting parameter is only $c$ here. The dependence of $c$ on various experimental conditions is also shown in Fig. 3(b). Interestingly, $c$ is almost not sensitive to any experimental condition, such as the vibration condition and type of granular material. This means that the conversion efficiency from vertical vibration energy into horizontal granular transport energy could be a universal constant.

III.2 Surface velocity

Next, we investigate whether the NDT model holds in particle scale. By the pattern matching of subsequent profile images, the velocity of a surface flow is measured at the same spatial and temporal resolutions as the depth-averaged velocity (Fig. 3(a)). The applied algorithm is similar to the particle image velocimetry (PIV) method Lueptow et al. (2000), although the profile data of this study is one dimensional. Small spacial windows for the profile pattern matching are given by $\pm 5~{}{\rm mm}$ at each position $(r=5,15,25,35~{}{\rm mm})$ as well as computation of slope $|\nabla h|$. Then, adequate and appropriate time interval $\Delta t$ is chosen so that the distance of displacement can be clearly identified as can be seen in Fig. 1(d). This distance is determined by finding the position where one-dimensional cross-correlation function between two profiles shows the maximum value. Here, we define $\Delta r$ as the $r$ component of this displacement (not along the surface) as the depth-averaged velocity $\bar{v}_{\rm t}$ is measured along the $r$ direction as well. Consequently, the surface velocity along the $r$ direction $v_{\rm t0}$ can be estimated by $\Delta r/\Delta t$. The range of $\Delta t$ is properly chosen in order for the typical $\Delta r$ to be a few grain diameter. The detail of this pattern matching is summarized in Appendix A.

Figure 3(c) shows the comparison between $v_{\rm t0}$ and $|\nabla h|$, where the data are scaled by a single curve as well as $\bar{v}_{\rm t}$ in Fig. 3(a). The data can also be fitted by the NDT model:

[TABLE]

where $c_{\rm 0}$ is a constant, physically corresponding to the conversion efficiency of the surface flow. The dependence of $c_{\rm 0}$ on experimental conditions is shown in Fig. 3(d). Note that the data in the range of $v_{\rm vib}>70~{}{\rm mm/s}$ is not plotted because the variation of surface profiles is too large to apply the pattern matching algorithm in these strong vibration conditions. Although the data range is slightly limited, the same flat trend as Fig. 3(b) is observed. However, the value of $c_{\rm 0}=0.11$ is different from $c=0.068$:

[TABLE]

which suggests that the flow velocity on the top of the relaxing pile is approximately twice as large as the depth-averaged velocity; and the flow has a structure that the velocity decreases as going deeper from the surface.

This tendency is qualitatively consistent with shear band structure of conventional heap flows Katsuragi et al. (2010), i.e., a clearly fluidized region is localized around the surface. However, the derivation of the NDT model assumes a uniformly-fluidized granular pile Tsuji et al. (2018). To verify the consistency between the observed results and the NDT model, the internal velocity profile of granular flows has to be investigated. To solve this matter, numerical simulations are much more convenient than experimental approaches for the setup of this study (Fig. 1(a)). In the next session, we will go into the discussion on numerical results.

IV Simulation result

IV.1 Depth-averaged velocity

The consistency with the experimental result is firstly confirmed for simulation data. Figure 4(a) shows the relation between the depth-averaged velocity $\bar{v}_{\rm t}$ and slope $|\nabla h|$ obtained in 2D and 3D simulations. While $\bar{v}_{\rm t}$ and $|\nabla h|$ are measured at three different points $(r=10,20,30~{}{\rm mm})$, other measurement methods are the same as the experiment (Fig. 3(a)). Both data are collapsed into single curves, and can be fitted by Eq. (1). It should be noted, however, that both $\mu$ and $c$ are left as free fitting parameters. Although $\mu$ was fixed by $\tan\theta_{\rm c}$ in the experiment Tsuji et al. (2018), $\mu$ can be in general different from $\tan\theta_{\rm c}$ Ghazavi and Mollanouri (2008) and also tends to decrease in the presence of vibration Caballero-Robledo and Clément (2009). In fact, $\mu=0.28~{}({\rm 2D})~{}{\rm and}~{}0.35~{}({\rm 3D})$ result in better agreements with the data of Fig. 4(a) than $\mu={\tan\theta_{\rm c}}=0.35~{}({\rm 2D})~{}{\rm and}~{}0.40~{}({\rm 3D})$. Besides, these fitting-based $\mu$ values are almost independent of simulation conditions (only within $\pm 0.02$) in the range of $v_{\rm vib}=30-80$ mm/s and $e=0.7-0.9$. Thus, the fixed values $\mu=0.28~{}({\rm 2D})~{}{\rm and}~{}0.35~{}({\rm 3D})$ are employed for all of the analyses of DEM simulations.

The $c$ values computed from the DEM simulation are plotted onto those obtained from experimental data in Fig. 4 (b). It seems that there is no significant difference between simulations and experiments in terms of the $c$ values. This result supports the fact that our DEM simulation can reproduce laboratory experiments well. In addition, $c$ does not depend on spatial dimensions (Figs. 4(a) and (b)). Therefore, 2D simulations are used to investigate the parameter dependence and analyze the velocity profile in Sec. IV.2.

The influence of the restitution coefficient $e$ is also investigated here. Figure 4(c) shows the relation between $c$ and $e$, which suggests that $c$ is not sensitive to $e$. This can be understood in the following sense: da Cruz et al. da Cruz et al. (2005) report that $e$ has no influence on the bulk frictional property in plane shear flows. In fact, the value of $\mu$ determined by the fitting of Eq. (1) is little influenced by $e$ as mentioned above. Given that the frictional property does not depend on $e$ in our system, the forces considered in the derivation of the NDT model are not influenced Tsuji et al. (2018). It is thus natural that the bulk flow property characterized by $\bar{v}_{\rm t}$ does not change depending on $e$.

IV.2 Velocity profile

The velocity profile inside the relaxing pile, which was technically challenging to address from experimental data, can be measured from the DEM simulation data. Here, we introduce the $z^{\prime}$ coordinate defined by the depth from the surface ($z^{\prime}=h-z$): $z^{\prime}=0~{}{\rm and}~{}h$ correspond to the surface and bottom, respectively. The velocity profiles are measured as a function of $z^{\prime}$ at various time and $r=10,20,30$ mm as Fig. 4(a) does. As an example, profiles obtained at $r=20~{}{\rm mm}$ and various time are shown in Fig. 4(d), where velocity profiles depend on not only $z^{\prime}$ but also $t$ as $h$ and $|\nabla h|$ vary with $t$.

On the analogy of steady heap flows Komatsu et al. (2001); Katsuragi et al. (2010), we found that the flow velocity $v_{\rm t}(z^{\prime},h,|\nabla h|)$ decreases exponentially with the depth $z^{\prime}$, implying the presence of shear band structure. The data are fitted by the following function:

[TABLE]

where $v_{\rm t0}$ and $h_{\rm s}$ represent the surface velocity and the characteristic thickness of a clearly fluidized layer. The dependence of $v_{\rm t}$ on $h$ and $|\nabla h|$ should be included in $v_{\rm t0}$ and $h_{\rm s}$.

An example of the fit by Eq. (18) is shown in the inset of Fig. 4(d). Despite being able to fit the data around the surface $(z^{\prime}=0)$ in a good quality, misfit to the exponential curve comes to remarkably appear as going deeper. This might be caused by the boundary effect as particles are trapped at the bottom, i.e., $v_{\rm t}\approx 0~{}{\rm at}~{}z^{\prime}=h$. Besides, when $h$ is too small, an exponential profile cannot be distinguished from another form of function such as a linear profile. To avoid these matters, only the upper halves ($z^{\prime}>h/2$) of sufficiently thick profile data ($h>5d$) are used for the fit by Eq. (18). Figure 4(e) shows all of the profile data, where $v_{\rm t}$ and $z^{\prime}$ are scaled by $v_{\rm t0}$ and $h_{\rm s}$. All of the profile data are collapsed into a single curve, which support the validity of an exponential profile (Eq. (18)).

Next, according to the experimental result (Fig. 3(c)), $v_{\rm t0}$ is a function of $|\nabla h|$, and can be fitted by Eq. (16). Figure 4(f) shows the relation between $v_{\rm t0}$ and $|\nabla h|$ obtained by DEM simulations. Although the data particularly at early $t$ scatter compared to Fig. 4(a), the NDT model for the surface flow (Eq. (16)) almost holds, suggesting the consistency with the experimental result. The best fitted $c_{\rm 0}$ values are also plotted in Fig. 4(c), which are approximately twice as large as $c$ (Eq. (17)) and do not depend on $e$. This relation also agrees with the experimental result (Eq. (17)).

Finally, the relation between $h_{\rm s}$ and $|\nabla h|$ is shown in Fig. 4(g), which suggests complex dependence. $h_{\rm s}$ is a function of not only $|\nabla h|$ but also the analysis position $r$ that can be read by $h(r)$. However, we empirically found that $h_{\rm s}/h$ is almost scaled by only $|\nabla h|$ as shown in Fig. 4(h), where a linear dependence can be observed:

[TABLE]

where $a$ is a constant, which is estimated as $a=2.41$ by means of the linear least-squares method. This trend means that the shear band gets enhanced as increasing the inclination. When $|\nabla h|$ is fixed, Figs. 4(g) and (h) suggest that $h_{\rm s}$ increases with $h$, but the ratio $h_{\rm s}/h$ is almost constant independent of $h$. The inset of Fig. 4(h) shows the analysis results obtained in various vibration conditions, where $h_{\rm s}/h$ is little influenced.

IV.3 Dependence on system size

From these characteristics of the shear band structure (Eq. (19) and Fig. 4(h)), it can be anticipated that the result is not influenced by the system size $R$ or the ratio $R/d$. In fact, the NDT model holds with the same value of $c$ when $R/d$ is varied by changing grain diameters ($d=0.5$ and $1.0$ mm) under the fixed system size ($R=40$ mm) in the experiment Tsuji et al. (2018). In contrast to this, here DEM simulations attempt to expand the understanding of the dependence on the system size by directly changing $R$.

The same DEM simulations as explained in Sec. II.2 are conducted under a larger plate $R=80$ mm with grains of $d=0.5$ mm. The experimental conditions are $A=0.04$ mm and $f=200$ Hz $(\Gamma=6)$, which are the same as those used in Fig. 4(h). The velocity profile ${v}_{\rm t}(z^{\prime})$ is measured at three positions ($r=20,40~{}{\rm and}~{}60~{}{\rm mm}$) and various time $t$, and the characteristic thickness of the shear band $h_{s}$ is estimated for each data as well. The result is shown in Fig. 5, which is compared to the data obtained with $R=40$ mm (same data as Fig. 4(h)). As expected, identical shear band structure expressed by Eq. (19) can be observed even if we change $R$ or $R/d$. In addition, as shown in the inset of Fig. 5, it has been confirmed that the bulk flow property characterized by $\bar{v}_{\rm t}$ is independent of $R$ or $R/d$, which is also consistent with Ref. Tsuji et al. (2018).

V Discussion

V.1 Consistency with the NDT model

Since DEM simulations have revealed the particle behavior inside the pile so far, the NDT model can be derived by integrating the velocity profile along the vertical direction. Using Eqs. (14), (18) and (19), the depth-averaged velocity $\bar{v}_{\rm t}$ can be computed as

[TABLE]

where

[TABLE]

Although the dependence on $h$ does not appear in $\bar{v}_{\rm t}$ as predicted by the NDT model, Eq. (20) is not equal to Eq. (1). However, as reported by Ref. Tsuji et al. (2018), the NDT model does not hold in smaller $|\nabla h|$ ranges, where the relaxation almost halts leaving finite slopes. In fact, Figs. 3(a) and 4(a) show that the model curves exhibit misfits to the data as the slopes approach zero. With respect to this reason, it can be speculated that the inertial energy supplied to grains is insufficient to overcome potential barriers of their neighbors in small $|\nabla h|$ ranges Jaeger et al. (1989); Roering (2004). On the other hand, this fact suggests that the NDT model is mainly suitable to predict the bulk flow property in a large $|\nabla h|$ range, where non-linearity clearly appears and the relaxation dynamics is practically dominated. We therefore focus on how Eq. (21) behaves in the vicinity of a divergence point ($|\nabla h|\simeq\mu$).

To this end, we introduce a variable $\epsilon=\mu-|\nabla h|$, and Eq. (21) is subjected to variable conversion:

[TABLE]

In order to observe the behavior of this function around $\epsilon\to 0~{}(|\nabla h|\to\mu)$, Taylor expansion is conducted as follows:

[TABLE]

where $a=2.41$ and $\mu=0.28$ are substituted. Therefore, in the limit of $\epsilon\to 0$, which is important for the fitting of the NDT model, Eq. (20) can be read as

[TABLE]

Since Eq. (1) is identical to Eq. (24) with Eq. (17), the NDT model has been consistently reproduced from the integration of the velocity profile obtained by DEM simulations. For small $|\nabla h|$ regimes, $c$ should decrease linearly as shown in Fig. 4(h) and Eq. (23). However, its effect is rather limited to discuss the practical relaxation of the pile.

Another remarkable point revealed by DEM simulations is that heap flows create shear band structure. In other words, clearly fluidized regimes are localized around the surface with thickness $h_{\rm s}$, below which creeping flows exhibit. This fact is in contrast to the assumption used in the derivation of the NDT model Tsuji et al. (2018) that the whole pile is uniformly fluidized as illustrated in Fig. 6(a). Actually, however, $h_{\rm s}$ is proportional to $h$ (Eq. (19)) when $|\nabla h|$ is fixed, which leads to a true image of the flow drawn in Fig. 6(b). Beside, Figs. 4(e) and 5 imply that the velocity profile is similar at any position and time independent of heap size as long as the vibration is strong enough to mobilize the whole granular pile. This characteristic differs from a conventional shear band structure with a constant thickness everywhere as observed in heap flows in the absence of vibration (e.g., Katsuragi et al. (2010)). Although the detailed pictures of Figs. 6(a) and (b) are different from each other, $h$ determines the fluidized thickness $h_{\rm s}$ in both cases. This is why $\bar{v}_{\rm t}=q/h(\propto q/h_{\rm s})$ can be described as a function of only $|\nabla h|$, and the NDT model derived on a basis of Fig. 6(a) is still valid for the granular flow consisting of peculiar shear band structure as shown in Fig. 6(b).

V.2 Potential applicability

It is also noteworthy that sediment transport from soil-mantled hillslopes shows a similar nonlinear property. The relation between sediment flux and hillslope gradient exhibits nonlinearity like Eq. (1). In other words, the flux increases divergently as the slope approaches a certain critical slope, which is reported by field observations Roering et al. (1999) and field measurements Gabet (2003), where environmental disturbance (e.g., earthquakes, rainsplash and biogenic activity) is considered to mobilize regolith particles. This natural process is mimicked by granular flows with acoustic noise in laboratory experiments Roering et al. (2001); Roering (2004); Furbish et al. (2008); and with random perturbation in DEM simulations BenDror and Goren (2018); Ferdowsi et al. (2018). Although the applied perturbation types of these studies are different from mechanically-controlled vibration used in our study, similar nonlinear transport properties are reported. From this similarity, it can be expected that the framework of our modeling for heap flows in the presence of vibration will be potentially applicable to other experimental configurations with different disturbance types.

V.3 Limitations and future works

Finally, several limitations of the model are discussed here. Although these limitations introduced below are beyond the scope of the present paper, they are important open issues left for future works.

The first limitation is that the vibration range where the NDT model can be applied is limited. This study focuses on the vibration conditions of $v_{\rm vib}=10-200$ mm/s and $\Gamma=2-10$. When increasing the vibration strength above this range, the transition into a granular-gas phase Jaeger et al. (1996) will occur, where the NDT model is no longer suitable. Conversely, as approaching a critical fluidization condition $\Gamma=\Gamma_{c}$, the NDT model will break down at some point, where the whole layer is not fluidized, i.e., the characteristic shear band drawn in Fig. 6(b) will not be created. The critical conditions to distinguish these multiple regimes will need to be investigated.

The second limitation is that the modeling in this study completely neglects the contribution of velocity fluctuation. Since the energy is transfered through a vibrating disk, the boundary conditions, such as those proposed in Ref. Richman (1993), should be satisfied, which also enables us to calculate a profile of granular temperature in the system. According to a kinetic theory for dense fluidized flows (e.g., Jenkins and Berzi (2012)), it is predicted that the viscosity, which connects the shear stress with shear rate, changes as a function of local granular temperature. Therefore, taking velocity fluctuation into consideration would be important to theoretically explain the velocity profile which should be governed by viscosity.

The third limitation is that the bulk frictional coefficient $\mu$, which is defined as the ratio of the shear stress to the pressure, is assumed to be constant. The dependence of $\mu$ on vibration conditions has been investigated by experiments Sánchez et al. (2007) and DEM simulations Khefif et al. (2018). These studies observe a spreading granular droplet under horizontal vibration, and report that $\mu$ varies as a function of an inertial dimensionless parameter $I=v_{\rm vib}/\sqrt{gd}$, or the square root of a shaking parameter $S={v_{\rm vib}}^{2}/gd$ Pak and Behringer (1993). From this, one might speculate that the dynamics of heap flows under vibration can be discussed on the analogy of local rheology Jop et al. (2005, 2006).

The last important open question is “what underlying nature determines $c$?” Since the value of $c$ is much less than $1$, most of the inputted energy is not used for the bulk granular transport. Moreover, strictly speaking, the value of $c$ (Fig. 3) shows a slight upward trend with $v_{\rm vib}$ when $f$ is fixed, although $c$ can be regarded as a constant approximately. The reason for this could be related to some missing factors described above. In any case, to solve this issue it is necessary to fully evaluate the energy partition among dissipation (inelasticity and friction), random motion (granular temperature), and collective motion (mean flow determining the value of $c$).

VI Conclusion

For the purpose of understanding the granular heap flow on a pile fully-fluidized by relatively strong vibration, we have experimentally and numerically studied the relaxation dynamics of a granular pile on a vertically-vibrating plate. To explain the relation between the depth-averaged velocity and local slope, the NDT model (Eq. (1)) has been proposed in Ref. Tsuji et al. (2018), which turns out to be applicable to the surface velocity (Eq. (16)) as well. These results are also satisfied in the DEM simulations, which support the universality of the modeling. The comparison of the model fitting parameters $c$ and $c_{0}$, which are constant independent of both experimental and numerical conditions, suggests that the surface velocity is approximately twice as large as the depth-averaged velocity (Eq. (17)). This result predicts that the flow velocity decreases as going deeper from the surface, which has been confirmed by measuring the internal velocity profiles obtained by DEM simulations. Moreover, it has been revealed that the relaxing pile creates shear band around the surface with exponentially-decreasing velocity profile (Eq. (18)). Its characteristic thickness, however, is not constant but proportional to the local height of the pile. We have also confirmed that these flow properties are independent of the system size. Finally, by integrating the exponential velocity function with this peculiar shear band structure from the base to the surface, the depth-averaged velocity described by the NDT model can be successfully deduced. Although these results are mostly based on empirical findings for now, the bulk transport picture proposed by Ref. Tsuji et al. (2018) has been consistently bridged to the particle-scale detailed picture revealed by this study.

Acknowledgments

This study has been supported by JSPS KAKENHI Grants No. 15H03707, No. 16H04025, No. 17H05420, No. 17J05552, and No. 18H03679.

Appendix A Algorithm of the profile pattern matching

Here the way of estimating $\Delta r$ is explained. Since $\Delta r$ is defined as the $r$ component of the distance that a profile moves during $\Delta t$, firstly a profile, which is actually inclined along slope, is projected onto the horizontal axis. To this end, as illustrated in Fig. 7, a profile $h(r)$ is fitted by a linear function ($g(r)=a_{0}+a_{1}r$, where $a_{0}$ and $a_{1}$ are fitting parameters), and then the profile along the $r$ axis $h^{\prime}(r)$ is created as

[TABLE]

$\Delta r$ can be determined by finding the position where one dimensional cross correlation function $f_{cc}(\delta r)$ between two profiles $h^{\prime}(r,t)$ and $h^{\prime}(r,t+\Delta t)$ shows the maximum value. The cross correlation function is defined as

[TABLE]

where $r_{s}$ and $r_{e}$ represent the start- and end-points of the spacial window considered for the analysis, and the normalization term $H^{\prime}(t)$ is given by

[TABLE]

Here, $r_{e}-r_{s}=10~{}{\rm mm}$ is employed as explained in Sec. III.2. Since the correlation function is normalized by ${H^{\prime}(t)}$ in Eq. (26), $f_{cc}(\delta r)=1$ corresponds to the complete match of two profiles. By changing $\delta r$ systematically, $\Delta r$ can be estimated, where $f_{cc}(\Delta r)$ shows the peak value of the correlation function. In this study, six different $\Delta t$ are applied to obtain reliable data, which increase consecutively at intervals of the temporal resolution of the profile measurement by the laser ($0.02$ s). Hereafter, six different $\Delta t$ are denoted by $\Delta t_{i}$, and their corresponding $\Delta r$ are expressed as $\Delta r_{i}$, where $i=1,2,\cdots,6$. $\Delta t_{i}$ is given by $\Delta t_{0}+0.02i~{}{\rm s}$, where a constant $\Delta t_{0}$ depends on the analysis position $r$. In general, the surface velocity increases with $r$ as a local slope gets steeper. Therefore, it is preferable for the accurate measurement that smaller $\Delta t_{0}$ is chosen for larger $r$ ranges. Table 2 shows the list of $\Delta t_{0}$ for various analysis positions $(r=5,15,25,35~{}{\rm mm})$. These values are chosen in order for the typical value of $\Delta r$ to be a few grain diameter.

Finally, the surface velocity $v_{\rm t0}$ is estimated using the weighted average:

[TABLE]

where $f_{cc}^{tot}=\sum_{i=1}^{6}f_{cc}(\Delta r_{i})^{2}$. Note that only the velocity with $f_{cc}(\Delta r_{i})>0.4$ is used to calculate $v_{\rm t0}$ so that the data obtained by low-correlated matching, which could be inaccurate, can be removed. We have confirmed that the result does not change once this threshold value becomes larger than $0.4$.