Identifying spatial transitions in heterogenous granular flow

Abstract

It is well known that heterogeneous granular flows exhibit collisional, dense and creep regimes that can coexist in space. How to correctly predict and control such complex phenomena has many applications in both mitigation of natural hazards and optimization of industrial processes. However, it still remains a challenge to establish a predictive granular rheology model due to the lack of understanding of the internal structure variation across different regimes and its interaction with the boundary. In this work, we use DEM simulations to investigate the internal structure of heterogeneous granular flow developed at the center of rotating drum systems. By systematically varying the side wall conditions, we are able to generate various heterogeneous flow fields under different levels of boundary effects. Our extensive simulation results reveal a highly relevant micro-structural quantity $\delta \theta = |\theta _c - \theta _f|$, where $\theta _c$ and $\theta _f$ are the preferred direction of inter-particle contacts and the preferred direction of inter-particle force transmissions, respectively. We show that $\delta \theta$ can characterize the internal structure of granular flow in collisional, dense and creep regimes, and its variation can identify the transition between them. In particular, in dense and collisional regimes, the classical rheological relation between bulk friction $\mu$ and inertia number I holds, while in the creep regime, such relation breaks down and $\mu$ instead depends on $\delta \theta$. Our findings hold for all investigated flow fields regardless of the level of boundary effect imposed, and regardless of the amount of shear experienced. $\delta \theta$ thus provides a unified micro-structural characterization for heterogeneous granular flow in different regimes, and lays the foundation of establishing microstructure-informed granular rheology models.

Appendices

Appendix 1: Discrete element model calibration and validation

In this section we discuss how we determine the model parameters used in our simulations. In general, the interaction between rigid particles can be modeled by either solving a linear-complementarity problem (implicit dynamics, the NSCD) [55] or by penalizing the inter-particle penetration (explicit dynamics, the classical DEM). Despite the different underlying principle, they give consistent results within the scope of rigid particle dynamics [56]. For the classical DEM, there are various inter-particle contact laws with different level of sophistication. Following the discussion in [57], we choose the linear Hookean contact law and pick $k_n = 2\times 10^5 \bar{m}g/\bar{d}$ (large enough to ensure rigid particle limit for the gravity-driven surface flows considered in our study) with $k_t = 2/7k_n$, $\gamma _t = 0$, and $\gamma _n = -2\text {ln}e\sqrt{\bar{m}k_n/(\pi ^2+\text {ln}^2e)}$. We are left to determine the coefficient of restitution e, the inter-particle friction $\mu _p$, the particle-wall friction and possibly the addition of rolling friction $\mu _r$.

1.1 First stage calibration via column collapse tests

We use the column collapse test to perform preliminary model calibration by measuring the angle of repose $\theta _r$.We first glue glass beads to the center area of an aluminum sheet ($300\times 300$ mm) and the inner surface of two identical iron angle bars with height 50 mm, length 25 mm and width 16 mm (Fig. 6a). A typical column collapse test can be divided into three steps: (1) filling with glass beads the hollow rectangular tube formed by the two iron angle bars placed over the aluminum sheet center, (2) rapidly removing the bars, and (3) taking picture of the formed pile to measure $\theta _r$. We repeat the procedure for 50 times and $\theta _r$ is measured to have a mean of $13.87^\circ$ and a standard deviation of $0.576^\circ$.

Fig. 6

a Setup of the column collapse test, b variation of $\theta _r$ according to the change of $\mu _p$ with a fixed $e = 0.82$, c variation of $\theta _r$ according to the change of e with a fixed $\mu _p = 0.4$, and c variation of $\theta _r$ according to the change of $\mu _r$ with fixed $\mu _p = 0.4, e = 0.82$

Via DEM simulations, we then perform numerical column collapse tests with the same configuration as in the experiments. We carry out two sets of simulations: (1) fixing $\mu _p = 0.4$ (a common choice for glass beads) and varying e from 0.1 to 0.82, and (2) fixing $e = 0.82$ [58] and varying $\mu _p$ from 0.1 to 0.8. From (1) we find that e has negligible effect on $\theta _r$ (Fig. 6c), and from (2) that $\theta _r$ first increases but later saturates with the increase of $\mu _p$ (Fig. 6b). In summary the above results suggest the necessity to incorporate rolling friction $\mu _r$, a parameter that imposes rotation hinderance [59] to model the interaction between non-spherical particles. Accordingly, we fix $e = 0.82, \mu _p = 0.4$ and vary $\mu _r$ from 0.01 to 0.2. Figure 6d shows the variation of $\theta _r$ against $\mu _r$. The numerical results indicate $\mu _r$ to be around 0.07 which is slightly larger than 0.01 in [60] where smooth glass spheres were used and slightly smaller than 0.1 in [61] where plastic spheres were used.

1.2 Second stage calibration and validation via rotating drum experiments

In the rotating drum experiments (Fig. 7a), after the surface flow becomes stationary under rotation speed $\omega = 2.59^\circ /\hbox {s}, 5.73^\circ /\hbox {s}$ and $11.23^\circ /\hbox {s}$, we use a high speed camera positioned against the glass plate to take images (288 px $\times$ 288 px corresponding to a 0.1389 mm/px resolution) for a time period of 10 s with a frame rate of 1000 fps. Accordingly the images cover an area of about 4 cm $\times 4$ cm at the drum center (Fig. 7b). From the sequence of images, we measure the dynamical angle of repose $\theta _d$ and the down-stream velocity $v_{yw}(z)$ near the glass plate. In terms of the former, we first binarize each image, then identify the pixels that represents the slope surface, and lastly use the identified pixels to perform a linear fit whose slope gives $\theta _d$ (Fig. 7d); as to the latter, we first use the open-source Particle Image Velocimetry (PIV) code [62] to compute the 2D velocity field $(v_1,v_2)$ by correlating boxes with dimension 8 px $\times 8$ px (corresponding to roughly $\bar{d}\times \bar{d}$), we then compute the velocity field under the frame rotating with the drum located at the drum center to get $(v_y ,v_z)$, lastly we compute $v_{yw}(z)$ by averaging $v_y$ within bands ($L_y = 20\bar{d}, L_z = 2\bar{d}$) positioned in parallel to y (slope surface) over a set of points that are placed every $\bar{d}$ distance along z (perpendicular to the slope surface) with $y=0$. Note that $|v_z| \ll |v_y|$ as the flow is nearly unidirectional. In the simulations, we generate images located at exactly the same location with exactly the same size (and resolution) as the ones taken from experiment (Fig. 7c), from which we follow the same image analysis procedure to find $\theta _d$. For $v_{yw}$, under the frame rotating with the drum at the drum center, we first pick particles located within $2\bar{d}$ away from the front-side flat wall, we then compute $v_{yw}(z)$ by averaging the particle velocity following the same procedure used in the experiments. Note that the width $2\bar{d}$ is picked to best represent the glass beads that are captured by the high speed camera.

Fig. 7

a The half-filled rotating drum with rear-side wall and inner-cylinder wall being glued with glass beads, b an image taken by the high speed camera at the drum center, c an image generated by numerical simulation with exactly the same location and size (resolution) as (b), d the binarized image of c for $\theta _d$ estimation, e time-averaged $\theta _d$ data for three different rotating speed $\omega$ estimated from experiments (red) and simulations (blue), f the down-stream velocity profile $y_{yw}(z)$ against the depth at the drum center calculated respectively for $\omega = 11.23^\circ /\hbox {s}$ from experiment (red triangle) and simulation (blue triangle), for $\omega = 5.73^\circ /\hbox {s}$ from experiment (red square) and simulation (blue square), and for $\omega = 2.59^\circ /\hbox {s}$ from experiment (red circle) and simulation (blue circle), and g the corresponding semi-log plot of (f). For e–g, the error bars represent the standard deviation associated with each time-averaged quantity

We use $\theta _d$ and $v_{yw}(z)$ measured with $\omega = 11.23^\circ /\hbox {s}$ for model calibration and the rest two for model validation. Prior to calibration, according to the column collapse test results, we fix $\mu _{\mathrm{p}} = 0.4$, and $e = 0.82$ that best represents the property of glass beads, although the latter has negligible effect for simulating steady granular flow [57]. Further, as the front-side plate is also made from glass, we fix the associated wall friction to be 0.4. The only left parameter to calibrate is the inter-particle rolling friction $\mu _{\mathrm{r}}$. Observing Fig. 6d, we vary $\mu _{\mathrm{r}}$ to be 0, 0.03, 0.05 and 0.07 and take the particle-wall rolling friction to be zero. As the front-side wall is flat, zero wall rolling friction is a reasonable choice. By solely using $\theta _d$ we identify that $\mu _{\mathrm{r}} = 0.03$ gives the best estimation ( $\mu _{\mathrm{r}} = 0$ underestimates $\theta _d$ while the rest two lead to overestimation). What’s more, when $\mu _{\mathrm{r}} = 0.03$, simulation and experiment show excellent agreement (Fig. 7f) in terms of $v_{yw}(z)$. The choice of $\mu _{\mathrm{r}}$ = 0.03 is then validated (Fig. 7e–g) by comparing both $\theta _d$ and $v_{yw}(z)$ obtained from simulations to those measured from experiments under $\omega = 2.59^\circ /\hbox {s}$ and $\omega = 5.73^\circ /\hbox {s}$.

Appendix 2: Additional simulation results

2.1 “S-shape” surface profile

For a direct comparison, Fig. 8 shows the surface shape profile for simulations performed under respectively periodic boundary condition, frictional side walls with $W/\bar{d} \simeq 21$ and wall friction 0.4 and that with $W/\bar{d} \simeq 10$ and wall friction 0.4. It can be observed that as the effective drum width is decreased from infinite (periodic boundaries) to $W/\bar{d} \simeq 21$ and to $W/\bar{d}\simeq 21$, the “S-shape” profile becomes more obvious under more significant side wall friction effect, especially when the rotation speed is large such as when $\omega = 67.38^\circ /\hbox {s}$ (Fig. 8j, o) in our case.

Fig. 8

Surface shape profiles. Images in first row from a–e are those under periodic lateral boundaries, in second row from f–j represents those under frictional side walls with $W/\bar{d} \simeq 21$, and in third row from k–o represents those under frictional side walls too but with $W/\bar{d} \simeq 10$. Wall friction is 0.4

2.2 Effect of lateral boundary condition on the co-directionality between s and $\dot{\gamma }^d$

Figure 9 showcases the spatial variation of the misalignment angle $\alpha$ for all considered drum configurations under $\omega = 33.69^\circ /\hbox {s}$, where $\alpha$ is defined as the angle between the principle directions of s and those of $\dot{\gamma }^d$. It can be observed that the presence of frictional side wall has a great impact on the value of $\alpha$, and the more frictional the side walls and the narrower the drum, the less co-directional are s and $\dot{\gamma }^d$. $\alpha$ is generally large (up to $60^\circ$) near the side walls (especially on the bumpy side), and beneath $z = z_\text {th}$ in the creep flow region. The green solid line indicates locations with local deformation $|\dot{\gamma }^d|/2\Omega$ equalling one [21] where $\Omega = \omega /360^\circ$; locations above this line has $|\dot{\gamma }^d|/2\Omega > 1$ while $|\dot{\gamma }^d|/2\Omega < 1$ for locations below this line. It thus may be understood that the large $\alpha$ in the creep flow region is due to the lack of shear. For locations above $z = z_\text {th}$, shear deformation is sufficient, and the misalignment $\alpha$ can be attributed to side wall perturbation.

Fig. 9

Spatial variation of $\alpha$ under $\omega = 33.69^\circ /\hbox {s}$ for all drum configuration considered; from left to right: $W/\bar{d}\simeq 10$ with wall friction of 0.4, periodic boundary, $W/\bar{d}\simeq 21$ with wall friction of 0.2, 0.4, 0.6 and 0.8, respectively. The grey solid lines indicate $z = z_\text {th}$ and the green ones represent locations with local deformation $|\dot{\gamma }^d|/2\Omega = 1$, where $\Omega = 360^\circ /\omega$

2.3 Tests of several non-local models

• The velocity fluctuation model Based on 2D numerical simulation results on granular flow in annular shear cells, the model proposed in [44] show that an improved relation between $\mu$ and I by adding the effect of velocity fluctuation is able to capture the failure of the one-to-one $\mu (I)$ relation in the creep flow region. The key prerequisite of this model is that the variation of $\mu -I$ and that of $\Delta -I$ are both not one-to-one and are mutually correlated, where $\Delta = |\delta v|\sqrt{\rho _s/P}$ is called the fluctuation number. However, Fig. 10a shows that for all flows considered, there is a global collapse between $\Delta$ and I. It appears that this model does not apply to the 3D flows considered in our study.

• The gradient expansion model In the gradient expansion model [17], the value of $\mu$ in the $\mu -I$ relation is modified by considering an additional contribution $\bar{d}^2\nabla ^2I/I$which is scaled by a phenomenological constant $\nu > 0$, assuming short range correlation between particle motions in different locations. Physically, the laplacian term captures on average, how does the I value at a certain point compares to its surrounding area. For instance, if a point is surrounded by a more fluid-like area (higher I), the laplacian term is positive and leads to the decrease of $\mu$ at that point. Since $\nu > 0$, essentially as long as $\nabla ^2I\ne 0$, the model will report a modification on $\mu$. Figure 10b shows the typical drum-width-averaged I and $|v|/\sqrt{g\bar{d}}$ against z with the case under $\omega = 33.69^\circ /\hbox {s}, W/d \simeq 21$ and wall friction of 0.4. The variation of I against z roughly follows the same trend as that of $|v|/\sqrt{g\bar{d}}$: it linearly decays in the collisional region and exponentially decays in the dense region. Thus in the collisional region $\nabla ^2 I = 0$ and the model reports no modification on $\mu$, which is consistent with our observation. However, in the dense region, $\mu$ will be modified according to the model since $\nabla ^2I \ne 0$—this contradicts our computations that confirm the applicability of the $\mu -I$ relation in the dense region. Again, it appears that this model does not apply to the 3D surface flows considered in our study.

• The fluidity model The fluidity model [16] implicitly modifies the value of $\mu$ in the $\mu -I$ relation by considering the nearby region contribution through the Laplacian of the fluidity parameter $g = \dot{\gamma }/\mu$ expressed as $\xi ^2\nabla ^2g$, where $\xi$ is defined as the cooperative length that diverges as approaching the jamming transition. Although the mathematical expression looks similar to that of the gradient expansion model, the assumed underlying physical mechanism is different. g is found to obey the following microscopic relation [29]: $\dot{\gamma }^{d}\bar{d}/\mu = \delta v F(\phi )$. Figure 10c shows how the normalized fluidity $\dot{\gamma }^{d}\bar{d}/(\mu \delta v)$, varies with volume fraction $\phi$. It can be observed for all data computed from $W/\bar{d} \simeq 21$ with varying side wall friction, they collapse well onto a single curve who shape resembles the one identified in [29]. However, this curve is clearly drum-width-dependent, as when the effective drum width is respectively infinite (periodic boundary, colored in blue) and $10\bar{d}$ (colored in yellow), no collapse can be observed, even in the dense and collisional region where the $\mu (I)$ rheology in its invariant form holds. Memory effect (insufficient shear) observed for $I < I_\text {th}$ (see Fig. 9) may explain the failure of the model in the creep flow region; while its break down in the fast flow regime ($I >I_\text {th}$ with sufficient shear deformation) reveals the non-trivial effects of side wall friction that have not been considered in the granular fluidity model [29]—indeed, even though investigated in 3D configuration, the considered flow fields have shear only in z direction as boundaries along both x and y direction are treated as periodic.

Fig. 10

a$\Delta -I$ data for all cases considered collapse onto a single master curve. Data are shown as the drum-with-average values with error bars representing the associated variation. b Variation of drum-width-averaged $|v|/\sqrt{g\bar{d}}$ and I against z with the case under $\omega = 33.69^\circ /\hbox {s}, W/d \simeq 21$ and wall friction of 0.4. c Relation between the normalized fluidity $|\dot{\gamma }^d|\bar{d}/\mu \delta v$ and volume fraction $\phi$. Data are shown as the drum-with-average values for clarity

2.4 Additional results from micro-scale analysis

The micro-structure can be investigated by three kinds of quantities that range from lower-order (L) to higher-order (H) within each kind: (1) geometry-associated ones range from volume fraction $\phi$ (L), coordination number Z (L) to the angular distribution of contact orientation (H) [47, 48]; (2) inter-particle-force-associated ones range from normal (and tangential) force p.d.f. distribution (L) [47, 48] to their angular distributions (H) [47, 48]; and (3) kinematics-associated ones where the lower order quantity can be the velocity fluctuation $\delta v$ [29, 44]. Inspired by [48, 63], we may regard $a_\chi ^v = (\lambda _3 - \lambda _1)/(\lambda _3 + \lambda _1)$ as the higher order quantity where $\lambda _3(\lambda _1)$ is the maximum (minimum) eigenvalue of the tensor $\chi = \langle |\delta v| n_i^v n_j^v \rangle /\langle |\delta v|\rangle$. Here “$\langle \cdot \rangle$” denotes the average over all particles and $n^v$ is the direction of velocity fluctuation associated with each particle. Accordingly $a_\chi ^v$ has range from 0 to 1 and reflects how “cooperative” the particle motions are: a value close to 1 in the quasi-static flow region implies the formation of “granular eddy” [24, 64].

Fig. 11

Spatial variation of the coordination number Z and $a_\chi ^v$ against the depth z. The symbol shape represents different rotation speed, while the symbol color represents different drum configurations in terms of drum width W and side wall friction. Error bars represent the variation of each shown quantity across the drum width

In principle, we are trying to find certain micro-structural quantities whose variations against I, (R1) exhibit a clear transition as passing through $z = z_{\text {th}}$, and more importantly, (R2) show consistence with that of $\mu$ against I for rheological considerations: when $I > z_{\text {th}}$ they are free of drum configuration effect and rotation speed effect while they become lateral boundary dependent and rotation speed dependent as soon as $I \le I_{\text {th}}$. Many of the aforementioned quantities, as we discover, only satisfy R1, such as Z and $a_\chi ^v$. Figure 11a shows the variation of Z against the depth z. Z varies case by cases when $z \le z_\text {th}$ and remains almost constant when $z>z_0$. However, on the contrary, the stress ratio $\mu$ shows a global collapse when $z \le z_\text {th}$ but varies case by case when $z>z_\text {th}$. $a_\chi ^v$ shows a slightly different spatial variation (Fig. 11b): it remains constant when $z\le z_\text {th}$ in a similar way to that of $\delta \theta$. However, $a_\chi ^v$ seems to be boundary condition dependent instead. When $z >z_\text {th}$, it rapidly increases independently with respect to the boundary condition and rotation speed, which is inconsistent with the variation of $\mu$ either. Following the momentum transferring argument, we have also investigated the portion of formed granular clusters based on either the velocity fluctuation following [5, 64] or the local volume fraction fluctuation (achievable via the computed Voronoi diagram) following [5, 14]. The variation of such portion along the depth show similar trend as that of $a_\chi ^v$ that does not satisfy (R2).

2.5 Determination of $z_\text {th}$ based on the drum-width-averaged velocity profile

Figure 12 shows the velocity magnitude profile against depth at the drum center |v|(z) and its variation across the drum width, for the same set of simulations considered here. When the lateral boundaries are periodic, |v|(z) shows negligible differences across the drum width. For the cases with frictional side walls ($W/\bar{d} \simeq$ 21 and 10), |v|(z) varies majorly as vertical translation without shape alteration. Therefore without sacrificing much the accuracy we determine $z_\text {th}$ based on the drum-width-averaged |v|(z) profile.

Fig. 12

|v|(z) profile variation across the drum width for different rotation speeds with each color representing a certain location x between $-8\bar{d}$ and $9\bar{d}$ ($W/\bar{d}\simeq 21$) and between $-3\bar{d}$ and $4\bar{d}$ ($W/\bar{d}\simeq 10$). First row from a–e: data extracted from simulations with periodic lateral boundaries; second row from f–j: data extracted from simulations with frictional side walls ($W/\bar{d}\simeq 21$); third row from k–o: data extracted from simulations with frictional side walls ($W/\bar{d}\simeq 10$); wall friction is 0.4