# Steady granular flow in a rotating drum: A theoretical nonlocal model for characterizing stress, velocity, and packing fraction profiles, encompassing grain shape effects from convex to highly concave

Weiyi Wang, Jonathan Barés, Mathieu Renouf, Emilien Azéma

► To cite this version:

Weiyi Wang, Jonathan Barés, Mathieu Renouf, Emilien Azéma. Steady granular flow in a rotating drum: A theoretical nonlocal model for characterizing stress, velocity, and packing fraction profiles, encompassing grain shape effects from convex to highly concave. Physical Review Research, 2024, 6 (4), pp.043310. 10.1103/PhysRevResearch.6.043310. hal-04871478

HAL Id: hal-04871478

https://hal.science/hal-04871478v1

Submitted on 7 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Steady granular flow in a rotating drum: A theoretical nonlocal model for characterizing stress, velocity, and packing fraction profiles, encompassing grain shape effects from convex to highly concave

Weiyi Wang<sup>1,*</sup> Jonathan Barés,<sup>1,†</sup> Mathieu Renouf,<sup>1,‡</sup> and Emilien Azéma<sup>1,2,3,§</sup>

<sup>1</sup>LMGC, Université de Montpellier, CNRS, Montpellier 34090, France

<sup>2</sup>Department of Civil, Geological, and Mining Engineering, Polytechnique Montréal, Montréal 2500, Canada

<sup>3</sup>Institut Universitaire de France (IUF), Paris 75005, France

<figure><img src="image_2.png" /></figure>

(Received 4 January 2024; accepted 25 October 2024; published 24 December 2024)

In this work, we investigate the effect of particle shape, ranging from spherical to highly concave, on steady flows in a rotating drum, a system that facilitates a continuous transition from a jamming state at greater depths to an unjamming state at shallower regions. We develop an analytical model to elucidate granular behavior within the rotating drum: (i) First, by decomposing the shear stress, we reconcile the discrepancy between simulation data and theoretical predictions, establishing a relationship with the angle of repose. (ii) Second, we extend the generalized Bagnold scaling, coupled with a nonlocal fluidity relation based on packing fraction, providing a framework for a correlation between shear stress, shear rate, and packing fraction. Additionally, we introduce a characteristic length to quantify the influence of particle shape and drum speed. This analytical model offers explicit functional forms for physical quantity profiles, which are validated experimentally in a thin rotating drum and numerically in a two-dimensional rotating drum. Our results demonstrate that this model accurately describes the change in velocity from the bottom of the drum to the free surface. Moreover, for different shapes of particle and drum speeds, the characteristic length captures the interplay between shear stress, shear rate, and the variation of packing fraction.

DOI: 10.1103/PhysRevResearch.6.043310

## I. INTRODUCTION: FLOW WITHIN A ROTATING DRUM

Granular matter is ubiquitous in both natural and industrial environments, featuring diverse complexity of shapes. Systems composed of these materials are commonly observed and studied both at the macroscopic and microscopic scales. One of their main characteristics is their ability to exhibit both solidlike (jammed) and fluidlike (unjammed) behavior [1–3]. Still, in this latter state, they do not behave as Newtonian fluid. Indeed, in contrast to the well-established constitutive relationships like the Navier-Stokes equations that aptly describe liquid flow, no single analytical model captures the evolution of flowing granular systems. In other words, the rheological behavior in granular flows differs from that of traditional fluids, representing an active area of intense investigation across several fields for decades [2–6].

One of the most complicated situations for this type of flow is when at least one surface of the system is not constrained

and is free to move, as the one observed near the surface of the flow in a rotating drum. In this latter case, there exist distinct flow modes that can simultaneously cohabit within the system. More specifically, in the so-called “rolling regime” (i.e., a flow governed by inertial effects with a well-defined dynamics [2,7–9]), we observe a solidlike behavior close to the walls and deep in bulk, where the particle motion is slow and considered as static. In contrast, in the upper part, a liquidlike behavior develops, characterized by layers of grains with larger inertia. In these layers, the flow evolves from a quasistatic regime just above the solidlike zone to inertial in the layers close to the free surface with an almost exponential increase in grain velocity. In this last layer, the flow is characterized by a nearly linear increase of grain velocity [2,8–12]. Finally, at the free surface particles mainly interact via binary collisions [8,9]. The transitions between each of these regimes, within a single system, are still a matter of debate [2,8,13,14]: At what depth do they occur? How thick they are? How do these thicknesses vary as a function of the drum speed and size? What about the properties of the grains?

This set of seemingly simple questions reveals a complex and still open problem: In the simplest and most common flow geometry, that of the rotating drum, there is no yet well-defined constitutive equations accounting for (i) the phase-flow transition and (ii) allowing to predict the stress and velocity fields. To address these issues, and more generally to provide a more scientific description of multidirectional granular flows in rotating drums, several paths may be explored.

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

<footnote>*Contact author: weiyi.wang@umontpellier.fr</footnote>

<footnote>†Contact author: jonathan.bares@umontpellier.fr</footnote>

<footnote>‡Contact author: mathieu.renouf@umontpellier.fr</footnote>

<footnote>§Contact author: emilien.azema@umontpellier.fr</footnote>
This paper is structured as follows: In Sec. II, we provide a comprehensive review of the state-of-the-art strategies for modeling granular flows, with a particular focus on drum geometry. This includes an examination of scaling laws, and local and nonlocal approaches. Section III delves into recent research regarding the influence of grain shape on flow behavior and how accounting for grain shape effects poses challenges to current rheological models. In Sec. IV, we introduce a theoretical model based on the concept of “fluidity” to predict stress and velocity fields in drum flow, considering the effect of grain shape. Sections V and VI detail the development of two-dimensional (2D) discrete element modeling and three-dimensional (3D) laboratory experiments to validate the new theories. Finally, in Sec. VII, we present our conclusions and outline future research directions.

## II. AN OVERVIEW OF MODELING STRATEGIES

### A. Scaling law approach

Scaling laws have been formulated to relate some macroscopic observables to microscopic features of the flows. These laws are most of the time established heuristically or empirically. For instance, in Govender’s comprehensive review [15], it is shown that the mean flow velocity $\langle v \rangle$ scales with the liquidlike flow thickness $h$ as $\langle v \rangle \propto h^m$, where $m$ varies from 1 to 3 [2, 10, 16–18]. The dynamic angle of repose $\langle \theta \rangle$ is also found to vary as a power law of angular speed $\Omega$ with exponent between 1 and 2.6 depending on the system properties [18–20]. Richer scaling laws have also been proposed, combining the Froude number ($F_r = R_d \Omega^2/g$, where $R_d$, $\Omega$, and $g$ are the cylinder radius, rotating velocity, and gravity acceleration, respectively) with other system parameters to describe the flow rates [21–23] or the dynamic angle of repose [23].

In fact, there is still no agreement on a set of equations to govern these systems. Crafting such equations is challenging due to the range of system parameters and flow scenarios involved.

### B. The $\mu(I)$ rheology and its limitation in multiphase flows

To analytically describe the evolution of a granular flow, a second way actively pursued consisted of building constitutive laws that better capture the main features of the flow rheology. After decades of efforts, compiling experimental and numerical data on flows in multiple configurations, a consensus has emerged in the form of the so-called “$\mu(I)$ rheology” [2]. This law empirically stated that the apparent macroscopic friction coefficient $\mu = \tau/P$ and the packing fraction $\phi$ both depend on the so-called “inertial number” $I$ [24–26], with $\tau$ being the shear stress and $P$ the confining pressure. The state parameter $I$ is defined as the ratio of the shear time $\dot{\gamma}^{-1}$ imposed by the flow rate $\dot{\gamma}$ over the particle relaxation time $(\rho d/P)^{1/2}$ for a particle of density $\rho$ and equivalent diameter $d$ [2]. It is shown that the shear stress is proportional to the pressure, through the effective friction coefficient as $\tau = \mu(I)P$, and the volume fraction $\phi$ is a function of $I$ as $\phi = \phi(I)$. In general, a linear dependence is obtained for small values of $I$ both for $\mu(I)$ and $\phi(I)$ [24]. More advanced functions have also been proposed to consider larger values of $I$ [3, 16, 17] or values depending on the geometry [13, 25, 27].

It is remarkable that this model agrees with Bagnold’s theory. Indeed, assuming that the momentum transfer between particles in adjacent layers results from binary collisions during flow [6, 28, 29] proposed

$$ \tau = f_1(\phi) \rho d^2 \dot{\gamma}^2, \quad (1a) $$

$$ P = f_2(\phi) \rho d^2 \dot{\gamma}^2, \quad (1b) $$

where $f_1$ and $f_2$ are functions depending only on $\phi$ such that $f_1(\phi) = \mu[I(\phi)]/I^2(\phi)$ and $f_2 = I^{-2}(\phi)$. Equation (1) has been recovered using dimensional analysis [28] and shown to be valid for all shear rates [30].

The $\mu(I)$ model when combined with a continuous approach (for the stress field) correctly predicts the velocity fields in various homogeneous flow geometries [3, 26]. Nevertheless, a number of limitations still remain [2, 3, 31]. Indeed, for instance, this model manages to characterize the rheology of the flow in the bulk material but fails near the walls, or for finite-size flows where boundary effects dominate [32, 33]. In the inclined plane flow, for example, the bottom and free surfaces are close and the angle at which the flowing layer comes to rest is determined by the layer thickness independently from the $\mu(I)$ prediction [2, 34]. Similarly, this model cannot accurately describe the transitions, from the quasistatic to the static regime and from the inertial to the collisional one [3].

This latter point is particularly critical when dealing with the flow in a rotating drum since these different regimes coexist. In this flowing geometry, the collinearity between the deviatoric shear stress and strain tensor is not observed, preventing the establishment of a general tensorial form of the constitutive law as in unidirectional flow [3, 35]. Nevertheless, the scalar form of the $\mu(I)$ relationship holds throughout the drum [13, 36, 37]. However, in cases such as these, or under conditions of high and low inertial numbers, the $\mu(I)$ law becomes “ill posed” due to its reliance on incompressible fluidlike behavior [38]. It is possible to modify the original $\mu(I)$ function to maximize the range of well posedness [39]. This approach has been successfully used to employ to simulate two-dimensional flow in square and triangular rotating drums containing mixtures of grains [40, 41]. In these simulations, the $\mu(I)$ rheology is active only in a surface layer, transitioning to a highly viscous Newtonian fluid deeper within the drum. This introduces a creep state, allowing grains deep within the drum to experience almost solid-body rotation. However, the theory in [39] can still become ill posed at high inertial numbers. New compressible theories have been developed [42–44]. However, as numerous studies have shown [13, 22, 36, 37], this does not imply that stress and velocity fields can be predicted using a continuum approach with a Bagnold-type profile, as is typically done for free-surface flow.

Another approach is still possible. By assuming a permanent flow, momentum balance equations along and normal to the flow direction yield a linear stress distribution across the depth as [2]

$$ \tau(z) = \rho g (H - z) \sin\langle\theta\rangle, \quad (2a) $$

$$ P(z) = \rho g (H - z) \cos\langle\theta\rangle, \quad (2b) $$
<figure><img src="image_3.png" /><figcaption>FIG. 1. (a) Typical shear stress profile $\tau(z)$ as a function of depth measured in drum geometry (red line). The gray zone represents the deviation of $\tau(z)$ from $\tau_0$ (gray dashed line) the shear stress given by Eq. (2). The blue line represents the additional shear stress from the drum geometry (see discussion in Sec. IV B). (b) Typical velocity profile $v_x(z)$ as a function of depth measured in drum geometry (red line). The gray zone represents the deviation of $v_x(z)$ from $v_0$ (gray dashed line) the typical Bagnold-type profile for a free-surface flow as in inclined plane geometry. The blue line represents the velocity profile imposed by the drum $v_{\text{drum}}(z)$ (see also Sec. IV B).</figcaption></figure>

where $g$ is the gravitational acceleration and $H$ is the height of the free surface along the $z$ direction (see Fig. 1). These relations are verified for flows down an inclined plane [2,3] and even for heap flows [2]. However Eq. (2a) fails in the case of flow within a rotating drum [13,14], while Eq. (2b) continues to be valid.

In Fig. 1(a) we present a sketch of the shear stress profile measured in the stationary rolling regime (full red curve) and its prediction from Eq. (2a) (gray dashed line). This latter, although acceptable in the first inertial layers, differs from the measurements deep in the quasistatic and static layers. In these layers, the shear stress tends to become independent of depth. It is often argued that the discrepancy between measurement and predictions comes from the fact that the flow is not perfectly homogeneous. As a major consequence, the effective coefficient of friction $\mu$ also varies with depth, and thus cannot be directly related to free-surface angle $\langle\theta\rangle$ [3,13,14]. Then, a constitutive law for the shear stress within a rotating drum is still missing.

Moreover, the streamwise velocity profile sketched in Fig. 1(b) shows a logarithmic decrease as it extends further from the inertial-flow zone deep into the static bulk zone. It moves away from the Bagnolian-type profile deduced from the $\mu(I)$ rheology [3,14]. Consequently, the $\mu(I)$ law does not apply to the drum case, requiring a nonlocal extension of the Bagnold scaling.

### C. Nonlocal rheology and high-order fluidity parameter

#### 1. Generalities

For some flow geometries, the previous observations emphasize the need to consider the nonlocal nature of the granular flow. What occurs at a given point in the flow is influenced by what diffuses from its immediate vicinity. To implement this in a coherent way, we introduce a characteristic length scale $l_m$ related to the particle diameter. This defines a scale of interaction that is distinct from the actual

size of the individual particles. Thus, a coefficient $\alpha$ is introduced such that $l_m = \alpha d$. This is assumed to account for the occurrence of collective motion of some clusters at different scales [2,12,31,45–47] and location within the flow. Then, following [2] Bagnold equation (1) is modified as

$$ \tau(z) = l_m^2(z) \rho \dot{\gamma}^2(z). \quad (3) $$

In the case of the rotating drum, it has been shown that this coherent length scale $l_m$ goes to zero at the free surface. Then it progressively increases as it approaches the transition into the static phase leading to a divergence close to the walls. This behavior drastically differs from other homogeneous like flows where $l_m$ remains most of the time relatively constant in bulk [2]. The challenge therefore lies in finding the “right” physical model to best describe how this coherent length scale originates and evolves as a function of the flow conditions.

As explained in the extensive reviews by Bouzid et al. [48] and Kamrin [49], a fruitful way to overcome this challenge has consisted of building an extra parameter describing the local state of the system. It relies on the introduction of an order parameter $f$ describing the dynamical transition from solidlike to fluidlike behavior. This order parameter quantifies the degree of fluidity displayed by a given region within a granular system. This parameter $f$ satisfies a partial Ginzburg-Landau-type phenomenological equation based on a diffusion process of the form [50–52]

$$ T \underbrace{\frac{Df}{Dt}}_{\text{evolution, material derivative}} = \underbrace{l^2 \nabla^2 f}_{\text{diffusive term}} + \underbrace{\mathcal{I}(f)}_{\text{source term}}, \quad (4) $$

where $T$ and $l$ are characteristic time and length scales, respectively. $\mathcal{I}(f)$ is a “source term” designed to switch the stability of the phase $f$ from solidlike to liquidlike behavior. Under this framework, Eq. (3) is often rewritten in a more general form as $f = \mathcal{F}(\tau, \dot{\gamma}, P, I, \dots)$. Then, the main issue consists of finding the control parameters and the source function $\mathcal{I}(f)$ that fit with the system dynamics.

One of the earliest approaches is the so-called “partially fluidized theory” (PFT), where $f$ is assumed to vary from 1, for solidlike behavior, to 0 for fluidlike behavior. The characteristic length scale $l$ is taken as the mean grain diameter [50,53,54]. The PFT manages to go further than the traditional $\mu(I)$ model and captures the solid-to-fluid phase transition, particularly concerning the initial and stopping heights along an inclined plane [50,53] or some specific flow patterns [55,56].

Then, the “nonlocal granular fluidity” (NGF) approach has emerged as a powerful tool for quantitatively predicting various nonlocal phenomena in various loading geometries [5,52,57]. The NGF can be seen as an extension of the work of Goyon et al. [4] and Bocquet et al. [58] where the fluidity is assumed to be the inverse of the viscosity: $f = \dot{\gamma}/\mu$. Also, the source term $\mathcal{I}(f)$ is chosen to vanish when $f$ tends to a local fluidity $f_{\text{loc}}$ that depends on $\mu$. The characteristic length scale $l$ represents long-range interactions diverging as the flow approaches the solid phase and going to 0 in the fluidized zone. More recently, Zhang and Kamrin [59] have formulated a microscopic interpretation that has then been verified experimentally [60] and found to be equivalent to the macroscopic form (under given boundary conditions) by Poon
and colleagues [61]. It turns out that the fluidity originates in the spatial fluctuations of the particle velocities. Nevertheless, the NGF is based on an explicit formulation of the $\mu(I)$ law which makes its applicability quite tricky in the drum geometry where stress information is missing. Moreover, both $\mu(z)$ and $I(z)$ are nonuniformly distributed in the drum.

Other strategies have been considered, for instance, by assuming that the fluidity is only defined from the inertial number [i.e., $f=I$ together with an expansion of the $\mu(I)$ law] [48,51].

## 2. Discussion: The case of the rotating drum

Nonlocal models for the flow in the rotating drum geometry (in particular to predict the stress and velocity fields) is still largely an open question. Nevertheless, some successes have been achieved (at least for the “rolling” regime) with the PFT approach. For instance, Horpe and Khakhar showed that the PFT (as developed initially by Arranson and Tsimring [50,53]) reproduces the velocity field within a rotating drum [14], although they do not perfectly fit the numerical and experimental data. Nonetheless, their model relies on four *ad hoc* adjustable parameters, together with a preliminary fit of the shear stress $\tau$, in the absence of a predefined stress model. However, their velocity model does not follow a Bagnolian profile.

At the same time as Aranson and Tsimring [50], Bonamy and Mills [62] developed their own nonlocal model. In this latter, the flowing granular medium is assumed to behave like a network of granular chains embedded in a Newtonian viscouslike fluid. By decomposing the total stress as a linear combination of Coulombic friction stresses, viscous stresses, and stresses due to immersed chains, they propose a simple equation based on solely two parameters (the shear rate and a characteristic length linked to the particle diameter). Their equation fits the data well in the liquid flow region, but does not describe the region near the free surface and fails to reproduce the transition to the solid region, particularly at high rotation speeds. Bonamy revisited the PFT approach, and modified the original model by forcing the viscosity to follow a Bagnold’s law [63]. He obtained a model that fits the experimental data much better than the earlier version, but the justification for maintaining a Bagnolian approach is missing.

# III. THE IMPORTANCE OF GRAIN SHAPE AND OVERALL OBJECTIVE OF THE PAPER

## A. Grain shape effects

Most of the results discussed in the previous sections have been obtained for model granular materials, i.e., materials composed of disks (2D) or spheres (3D). Whether they are minimalist like the $\mu(I)$ model, or more sophisticated as the “fluidity” concept, most of the models are all able to explain complex collective phenomena, but they mainly ignore the fact that the complexity can arise from the materials itself. The intrinsic properties of the grains composing the system can dramatically change its rheology [64–71].

For example, in the quasistatic limit (i.e., $I \ll 1$), the shear strength $\mu$ increases with the grain angularity. But for larger grain angularities it may saturate towards a maximum

value [72,73] or even decrease toward values close to that of an assembly of disks, depending on the contact friction value [74]. Numerous systematic studies have also highlighted nonlinear variations in the solid fraction with grain elongation [64,65] or grain shape nonconvexity [65,66,68,75,76]. In the case of homogeneous flows, the $\mu(I)$ rheology is still well observed, at least for convex [77] or slightly concave [75,78] grains.

However, for the same thickness, the avalanche initiation angle is greater for nonspherical grain packing [77,78], demonstrating that the shape of the grains modifies the local properties. In silolike geometries, faster discharges were observed for elongated particles than for rounded particles [79]. In contrast, grain discharges slow down with increasingly angular grains [80]. This again reveals nonlocal effects where elongated grains manifest nematic ordering [81,82], facilitating their ejection, whereas increasingly pronounced arching effects occur with angular grains induced by face-to-face contact [72]. In drum geometry two distinct flow regimes were observed depending on rotation speed and grain shape [83]: (i) at high rotation speeds irregularly shaped grains led to a higher dynamic angle of repose caused by the granular packing expansion, (ii) at low rotation speeds, on the contrary, the angle is more related to the packing fraction. For elongated and flat particles, differences in the velocity profile measured in the active layer were reported as the aspect ratio increases [84]. Similarly, nonspherical particles showed less axial dispersion than glass beads [85].

Among all the possible grain shapes, concave grains are of particular interest. Indeed, highly concave particles demonstrate steep transitions between different dynamic states [69]. They can even maintain stability and loading without any external confinement through their interlocking capability and long-range correlation [86–90].

At the same time, they exhibit very low packing fraction [66,76]. For instance, the column stability of concave U-shaped particles shows that the decrease in particle packing fraction offsets the increase in entanglement with concavity, these two trends conspire to generate a maximum of resistance to separation in collections of nonconvex particles of intermediate aspect ratio [67]. Very recently, it has also been reported that sheared hexapods develop a secondary flow profile that is completely opposite to that of convex grains in the same geometry [91]. This makes the systems composed of these particles an ideal system for exploring the fundamental principles governing transitions arising from grain shape properties.

The task is very challenging. In addition to theoretical aspects aimed at generalizing rheological and nonlocal models, systematic studies need to be developed to continuously assess the effect of grain shape change (from convex to highly concave) on rheological properties. From an experimental point of view, a fundamental challenge is to design a sufficient number of particles while systematically controlling their shapes and mechanical properties. In addition, accounting for particle shape in numerical simulations using discrete element methods presents a number of technical hurdles that are both geometrically and computationally complex. An example is detection of contacts and the calculation of forces between particles of arbitrary shape, particularly concave ones, which can have a large number of contacts.
<figure><img src="image_4.png" /><figcaption>FIG. 2. A granular system is flowing in a half-filled rotating drum of radius R and lateral thickness W. The unit vector $e_x$ oriented with an angle $\theta$ from the horizontal is parallel to the free surface and $e_z$ is perpendicular to it. $\Sigma$ (in red) is a narrow area in which the flow is considered invariant along $x$ in the steady-state regime. $q(z)$ represents the value of any system observable, averaged in space and time in the green area.</figcaption></figure>

## B. Objectives and outline of the paper

In this work, we investigate the effect of grain shape non-convexity on the intrinsic rheology of a “rolling” steady-state flow within a rotating drum. To this end, we design 2D numerical studies and 3D experiments in which the grain shape is systematically varied from circular (and spherical, respectively) to cross-shaped particles with very thin arms. The drum speed is also systematically varied to explore different rolling regimes. Ultimately, we aim to develop a nonlocal model for predicting the stress and velocity fields. This model has to be sufficiently general to capture the effects that grain shape can have on the different flow regimes (depth and thickness) observed in the drum geometry.

# IV. ANALYTICAL MODEL

## A. Notations and main assumptions

In the rotating drum shown in Fig. 2, we define the reference frame $\mathcal{R}_0 = (e_x, e_z)$, where $e_x$ (respectively $e_z$) is parallel (respectively perpendicular) to the free surface. This latter is inclined of an angle $\theta$ with the horizontal. In the reference frame, the flow can be considered as quasihomogeneous at the center of the drum, within a thin elementary slice $\Sigma$ (red area in Fig. 2), parallel to $e_z$ and located around $x = 0$. The slice can be divided into elementary horizontal layers stacked on top of each other and parallel to the flow, one such layer is indicated as green in Fig. 2. The value of a given continuum quantity $q(z)$ (packing fraction, velocity, stress, etc.) at depth $z$ is then defined as the average of the corresponding quantity given at the grain scale for all the particles, whose center of mass is inside the layer. We assume that the flow is sufficiently steady and quasihomogeneous across the flowing layer within the slice $\Sigma$, so that

$$ 
\begin{aligned}
\frac{\partial q(z)}{\partial t} &\approx 0, \\
\frac{\partial q(z)}{\partial x} &\approx 0.
\end{aligned}
\qquad (5) $$

We also assume the packing fraction to be a global variable $\phi$ in the granular flow. We define it as the ratio of the volume of the particles to the volume of the packing $\phi = V_{\text{grain}}/V_{\text{packing}}$.

Finally, the lateral thickness $W$ is assumed to be wide enough with respect to the apparent diameter of the grains $d$ to prevent side effects induced by lateral walls.

## B. Stress profile

Starting from the Cauchy momentum equation that describes the momentum transport in any continuum, we get [92]

$$ \underbrace{\frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla) \mathbf{v}}_{\text{material derivative, } Dv/Dt} = \frac{1}{\rho} \underbrace{(\nabla \cdot \boldsymbol{\tau} - \nabla P)}_{\nabla \cdot \boldsymbol{\sigma}} + \mathbf{g}, \qquad (6) $$

where $\mathbf{v}$ is the velocity vector field of the flow [depending on time $t$ and space $(x, z)$], $\rho$ is the material density, $g$ is the gravitational acceleration, and $\boldsymbol{\sigma}$ is the stress tensor that we can decompose into $P$ the hydrostatic pressure and $\boldsymbol{\tau}$ the deviatoric stress tensor.

When studying the drum geometry, it is classical to decompose the velocity $\mathbf{v}$ into two distinct sources using the superposition principle [2,13,14,35]: one linear component following the same profile as the rotating drum $\mathbf{v}_d$ [blue curve in Fig. 1(b)], and another, nonlinear one, representing the velocity relative to the drum $\mathbf{v}_r$. This gives

$$ \mathbf{v} = \mathbf{v}_r + \mathbf{v}_d. \qquad (7) $$

According to the first principle of special relativity, stating that all physical laws take their simplest form in an inertial frame and that there exist multiple inertial frames interrelated by uniform translation [93], we postulate that the same holds true for the stress tensor leading to the following decomposition:

$$ \boldsymbol{\sigma} = \boldsymbol{\sigma}_r + \boldsymbol{\sigma}_d, \qquad (8) $$

where $\boldsymbol{\sigma}_r$ is the stress tensor relative to the drum, and $\boldsymbol{\sigma}_d$ the stresses resulting from the rotation induced by the drum. However, while in Eq. (7) at least $\mathbf{v}_d$ is well defined, in Eq. (8) neither $\boldsymbol{\sigma}_d$ nor $\boldsymbol{\sigma}_r$ are known a priori.

As a first approximation, in the drum reference frame ($\mathcal{R}_0$), the flow can be considered similar to the flow on a heap (assuming no lateral wall effects [16,94]). This analogy can be invoked thanks to the so-called self-similarity of the velocity profiles in these two configurations [2,94,95]. Thus, according to our main assumptions [see Eq. (5)], substituting $\boldsymbol{\sigma}_r$ and $\mathbf{v}_r$ within Eq. (6) in the reference frame $\mathcal{R}_0$ leads naturally to the following equations [see also Eq. (2) in Sec. II B]:

$$ \tau_r(z) = \rho g(H - z) \sin\langle\theta\rangle, \qquad (9a) $$

$$ P_r(z) = \rho g(H - z) \cos\langle\theta\rangle, \qquad (9b) $$

where $\tau_r(z)$ and $P_r(z)$ are the average shear and normal stresses components of $\boldsymbol{\sigma}_r$, respectively. It should be noted that Eq. (9) is always verified in heap flow geometry from the free surface deep into the static bulk [96]. Thus, since the repose angle $\langle\theta\rangle$ remains constant, on average, within $\mathcal{R}_0$, the relative effective coefficient of friction $\mu_r$ remains constant too and
can be written as

$$ \mu_r = \frac{\tau_r(z)}{P_r(z)} = \tan\langle\theta\rangle. \qquad (10) $$

Moreover, it is also reported in the literature that Eq. (9b) fits correctly the normal stress distribution in the rotating drum geometry [13]. Thus, it can be stated that $P(z) = P_r(z)$, where $P(z)$ is the hydrostatic component of $\sigma$. This implies that $\sigma_d$ carries only the shear stress $\tau_d$.

Therefore, we can define a “basal” friction coefficient $\mu_d$. This latter is transmitted only by the drum and results from wall effects that propagate through the material [97]. Accordingly, it can be written as

$$ \tau_d(z) = -\mu_d P(z). \qquad (11) $$

The negative sign comes from the fact that the orientation of $\tau_{\text{drum}}(z)$ must correspond with the outer layer drag speed. Identifying the basal friction law is a complex task and has led to different strategies depending on the studied geometry [97–99].

In our case, we formulate these statements. (i) First, by construction, $\mu + \mu_d = \mu_r = \tan\langle\theta\rangle$, where $\mu = \tau/P$, with $\tau$ the shear stress extracted from $\sigma$. Since it is known that $\mu = \mu(z)$ evolves linearly with $z$ in the drum geometry [36] (see also Appendix A), then $\mu_d = \mu_d(z)$. More accurately $\mu_d$ is linear with $z$. (ii) Second, $\mu_d$ is assumed to vanish at the free surface for a sufficiently large system, meaning that $\mu_d|_{z=H} = 0$. (iii) Third, near the bottom wall $\mu_d|_{z=0} = \mu_w$, with $\mu_w \in ]0, \tan\langle\theta\rangle[$ an effective coefficient of friction between the flowing layer and the bottom. All these conditions lead to the following:

$$ \mu_d(z) = \mu_w \frac{H-z}{H}. \qquad (12) $$

Consequently, substituting this expression of $\mu_d$ in Eq. (11) gives

$$ \tau_d(z) = -\mu_w P(z) \frac{H-z}{H}. \qquad (13) $$

Interestingly, Eq. (13) is reminiscent of the shear stress profiles induced by the walls in vertical-chute flows [2,100–102]. This implies a nontrivial analogy where, at least for stresses, the flow within a rotating drum can be viewed as an intricate combination of heap flow and pipe flow. Finally, noting that the basal conditions near the wall are necessarily the same for $\mu$ and $\mu_d$ (i.e., $\mu|_{z=0} = \mu_d|_{z=0}$), we get that $\mu_w = \frac{1}{2}\tan\langle\theta\rangle$. Thus, inserting this last expression into Eq. (13), together with Eq. (9b) [recalling that $P(z) = P_r(z)$], we get an explicit formulation of $\tau_{\text{drum}}$ of the form

$$ \tau_d(z) = -\rho g \frac{(H-z)^2}{2H} \cos\langle\theta\rangle. \qquad (14) $$

A sketch of $\tau_{\text{drum}}(z)$ is shown in blue in Fig. 1(a). Then, combining Eqs. (14) and (9) together with (8), the shear stress component $\tau$ within a rotating drum can be written as:

$$ \tau(z) = \rho g \frac{H^2 - z^2}{2H} \sin\langle\theta\rangle. \qquad (15) $$

Finally, we define the granular flow density at coordinate $z$, $\rho(z)$, using the packing fraction of the assembly $\phi(z)$ and the

grain density $\rho_0$: $\rho(z) = \rho_0\phi(z)$. We are then able to scale the shear and normal stresses within a rotating drum to get their theoretical expression as follows:

$$ \frac{\tau(z)}{\rho_0 g d} = \phi(z) \frac{H^2 - z^2}{2Hd} \sin\langle\theta\rangle, \qquad (16) $$

$$ \frac{P(z)}{\rho_0 g d} = \phi(z) \frac{H-z}{d} \cos\langle\theta\rangle. \qquad (17) $$

### C. Velocity profile: A nonlocal model based on the packing fraction

In this section, based on recent work using the PFT approach, we develop a nonlocal velocity model that is (i) compatible with Bagnold’s scaling law and (ii) accounts for the effect of the grain shape.

#### 1. Preliminary step: Equivalence between fluidity and packing fraction profiles

We start from the Ginzburg-Landau-type phenomenological equation [Eq. (4)] presented in Sec. IIC1. Following Aranson and Tsimring [50], the source term can be chosen as $\mathcal{I}(f) = f(1-f)(f-\delta)$, where $\delta$ is a function with value between 0 and 1 defined as $\delta = (\mu - \mu_{\text{dyn}})/(\mu_{\text{sta}} - \mu_{\text{dyn}})$. It is built from $\mu_{\text{sta}}$ and $\mu_{\text{dyn}}$ the dynamic and static friction coefficient respectively, i.e., the friction coefficient of a given granular material measured at small and high inertial number, respectively. This choice is motivated by the simple fact that the source terms must have extrema both at $f = 0$ (fluidlike) and $f = 1$ (solidlike). The parameter $\delta$ is introduced to control the range in which both static and dynamic phases coexist. In practice, $\delta$ is set to 0.5. Under the hypothesis of a steady and fully developed flow given in Eq. (5), and the definition of the source term given above, an explicit expression for $f$ can be derived with the form [50]

$$ f(z) = \frac{1}{2} \left\{ 1 - \tanh \left( \frac{z - z_0}{\sqrt{8}l} \right) \right\}, \qquad (18) $$

where $z_0$ is such that $f(z_0) = \frac{1}{2}$. This corresponds approximately with the depth at the transition between the quasistatic and the inertial flow regimes (see Fig. 3). Note that, unlike Aranson and Tsimring [50], we do not assume that the characteristic length $l$ is related to the mean diameter $d$ of the grains. On the contrary, we assume that $l$ is a critical length depending on the shape of the grains. Therefore, this length scale is considered as a free parameter in our approach.

It is also well documented in the literature that the evolution of the packing fraction $\phi(z)$ within a rotating drum follows a relatively simple trend as a function of $z$, as depicted in Fig. 3 [35,94]. In this figure, $\phi$ is constant and equals to $\phi_s$ in the solid phase of the drum. Then, it decreases linearly with $z$ in the liquid phase from $\phi_s$ at $z = z_s$ (i.e., at the transition from the static to the liquid regime) to $\phi_c$ at $z = z_c$ (i.e., at the transition from the liquid to the collisional regime). Finally, it tends towards zero at the surface, for $z > z_c$. Since $z_0$ is assumed to be the depth at the transition between the quasistatic and inertial regimes, we define $\lambda = z_0 - z_s$ the thickness of the quasistatic flow zone (see Fig. 3). From these definitions, it is easy to obtain the following expression
<figure><img src="image_5.png" /><figcaption>FIG. 3. Illustration of the packing fraction profile $\phi(z)$ as expressed by Eq. (19), where $z_s$, $z_0$, and $z_c$ are the depths at the transition from solid to quasistatic and finally collisional states, respectively. $\lambda = z_0 - z_s$ is the thickness of the quasistatic flow zone.</figcaption></figure>

that gives the evolution of the packing fraction with depth:

$$ \frac{\phi_0 - \phi(z)}{\phi_s - \phi_0} = \frac{z - z_0}{\lambda}, \qquad (19) $$

where $\phi_0 = \phi(z_0)$. Now, multiplying the left- and right-hand sides of Eq. (19) by $\xi = \lambda/(\sqrt{8l})$ and inserting it into Eq. (18) we obtain a definition of the fluidity $f$, where the dependence on $z$ appears only through $\phi(z)$. This permits to define a functional $F$ depending only on $\phi$:

$$ F(\phi) = \frac{1}{2} \left\{ 1 - \tanh \left( \xi \frac{\phi_0 - \phi}{\phi_s - \phi_0} \right) \right\} \equiv f(\phi(z)) \qquad (20) $$

Thus, we see that the fluidity parameter can be described in terms of the packing fraction $\phi$ only, by considering the functional $F$ through the global parameter $\xi$. This latter allows linking of two, a priori, unknowns parameters that both depend on the grain properties (i.e., shapes, sizes, ...): the thickness of the quasistatic flow $\lambda$ and the characteristic length parameter $l$. Considering that $F(\phi) \to 1$ when $\phi \to \phi_s$, together with the form of the tanh($\xi$) function, this suggests that $\xi$ can be chosen to be greater than at least 2, but at the same time cannot tend to infinity since $\lambda < z_c - z_s$. Thus, $\xi$ is a fitting parameter that has to be determined numerically by adjusting Eq. (20) on (18), and, contrary to $\lambda$ and $l$, it is not expected to depend on the grain properties. Thus, it can be called a “universal” parameter.

## 2. Velocity profile built from the fluidity

Now we aim to generalize the Bagnold scaling by coupling it with the fluidity function. Based on the works mentioned in Sec. II B, two remarks can be made.

On the one hand, it appears that in its general Bagnolian form, the shear stress is related to the flow rate through a function $f_1(\phi)$ that depends on the packing fraction [see Eq. (1a)]. Actually, several studies have evidenced that the function $f_1$ exhibits a divergence as the packing fraction $\phi$ approaches its maximum value $\phi_s$. More precisely, there is evidence that it is necessary for $f_1(\phi)$ to scale with $(\phi_s - \phi)^{-2}$ in the vicinity of $\phi_s$ [6,29]. Thus, incorporating the concept of fluidity into Eq. (1a) together with the aforementioned conditions on $f_1(\phi)$

a natural expression of $f_1(\phi)$ is

$$ f_1(\phi) = \frac{k^2}{[1 - F(\phi)]^2}, \qquad (21) $$

where $k$ is a dimensionless constant that does not depend, a priori, on the grain properties that are naturally captured by $\phi$. $k$ can be considered as a second “universal” parameter. The square form of $k$ in this equation will be justified.

On the other hand, and as already discussed, several authors have proposed to modify the Bagnolian equation by introducing a characteristic length scaling, through a parameter $\alpha$, with the particle diameter $l_m(z) = \alpha(z)d$. This scaling depends, a priori, on the flow depth $z$ [see Eq. (3)]. Thus, by assuming that this description relies essentially on the same physical mechanisms as those of Eq. (1a), we can assume that $\alpha^2(z) = f_1(\phi(z))$, and thus we get

$$ \alpha(z) = \frac{k}{1 - F(\phi(z))} = \frac{k}{1 - f(z)}. \qquad (22) $$

Therefore, we can generalize the Bagnold scaling from Eqs. (3) and (22), which elucidates the intricate interplay between shear stress, shear rate, and packing fraction in the form

$$ [1 - f(z)]^2 \tau(z) = \rho(kd)^2 \left( \frac{\partial v_x(z)}{\partial z} \right)^2. \qquad (23) $$

Finally, inserting the expression of the shear stress [Eq. (16)] and the expression of the fluidity [Eq. (18)] into Eq. (23), and remembering that $\lambda = z_0 - z_s = \xi\sqrt{8l}$, we obtain the following expression for the derivative of the velocity profile:

$$ \frac{\partial v_x(z)}{\partial z} = \frac{1}{2kd} \sqrt{g \frac{H^2 - z^2}{2H} \sin\langle\theta\rangle} \left[ 1 + \tanh\left(\xi \frac{z - z_s}{\lambda} - \xi \right) \right] - \Omega. \qquad (24) $$

Thus, we are able to give a theoretical expression for the velocity field of the flow within a rotating drum. The above equation is based on four parameters that must be adjusted. Two of them, $k$ and $\xi$, are presumed to be independent of the grain properties. On the contrary, the two others, $z_s$ and $\lambda$, depend on the grain shape and describe the different thicknesses and depths of the flow zones. For instance, $z_s$ can be easily determined from the packing fraction profile. To determine the three other parameters, we can rely on the relationship between $f(z)$ and $F(\phi)$ in Eq. (20) since Eq. (18) can be simply rewritten as a function of $\xi$ and $\lambda$ only.

# V. NUMERICAL VALIDATION

## A. Discrete-element modeling, system parameters, and steady state

Two-dimensional simulations were carried out using the contact dynamics (CD) method [103]. The CD method is a discrete-element method (DEM) in which small-scale effects are considered into nonregularized contact laws together with a nonsmooth formulation of the particle dynamics. In other words, contrary to the strategy adopted by the molecular dynamic approaches (also called soft-DEM) [104], no elastic repulsive potential nor smoothing of the Coulomb friction law
<figure><img src="image_6.png" /><figcaption>FIG. 4. (a) Example of 2D star-shaped particles with varying concavities (η). (b) Snapshot of the rotating drum in the permanent rolling flow regime. The red zone shows the area Σ in which the averages are calculated. (c) Zoom within the rotating sample for η = 0.9. The arrows show the velocity field.</figcaption></figure>

is used to determine contact forces. The unknown variables are the particle velocities and contact forces simultaneously found using a (parallelized [105]) iterative algorithm based on a nonlinear Gauss-Seidel scheme. Finally, the equations of motion are integrated by an implicit time-stepping scheme. This method is numerically unconditionally stable and particularly well adapted for simulations of a large assembly of frictional particles. This is particularly true for complex, nonconvex-shaped particles where there may be multiple points of contact between two particles. The method has been extensively employed for the simulations of granular materials in two and three dimensions. For a detailed description of the CD method, see [103]. For our simulations, we used the simulation platform LMGC90 developed in our laboratory [106,107].

In this study, we consider cross-shaped particles with rounded-cap ends on each branch extremities as illustrated in Fig. 4(a). Such shapes can be easily described using a “concavity” parameter [65] defined as

$$ \eta = \frac{r - r_0}{r}, \qquad (25) $$

where r is the radius of the circumscribed disk and 2r₀ is the branch thickness. For technical reasons, we consider two ways of modeling particles. In the cases where η ≤ 0.5, the grains are modeled as four overlapped disks of radius r₀ whose centers lay at the corners of a square of edge r₀√(2). In the cases where η > 0.5, the grains are built with two rectangles

of length L = 2(r - r₀) and four disks of radius r₀ at the ends of the rectangles. The contacts between two cross-shaped grains can be reduced to a combination of contacts between disks for η ≤ 0.5. In contrast, for η > 0.5 three situations may arise: cap-to-cap, cap-to-line, and line-to-line contacts. Cap-to-cap and cap-to-line contacts are considered as one contact point (i.e., disk-disk or disk-polygon contacts, respectively). In the framework of the CD method, it is common to represent line-to-line contacts as two contact points. This is what is done here by considering two cap-to-line contacts. The implementation of line-to-line contacts in the framework of the CD method is described in detail in Ref. [72]. In the following, the concavity parameter η is varied from 0 (disk) to 0.9 in steps of 0.1.

Nₚ randomly oriented grains of radius r are first laid under the action of the gravity within a drum of radius 150r [see Fig. 4(b)]. The number Nₚ of grains ranges from 9104 for η = 0 to 20 707 for η = 0.9 so that the drum is half-filled whatever η. A small particle size distribution is introduced around the mean radius r̄ to avoid crystallization. Typically, the smallest particles have a radius of r = 1.6 mm while the largest have r = 2.4 mm. The friction between grains is set to 0.2 while the friction between the grains and the drum is set to 0.9. This prevents slipping at the boundaries. Finally, in order to isolate the effect of grain shape, in all the simulations the grain mass is kept constant by adjusting the material density ρ₀ for each η.

Then, a constant angular velocity Ω is applied to the drum. Ω is varied in [2,4,6,8,10] rpm for each values of η. For every configuration, we make sure that a permanent steady-state regime is reached. Under rotation, all systems dilate slightly. The free surface is then no longer parallel to the reference frame but inclined with an angle θ as illustrated in Figs. 4(b) and 4(c). Figure 5(a) shows the evolution of θ for Ω = 2 rpm (Ω = 10 rpm in the inset) as a function of the number of revolutions and for different values of η. For each pair (η, Ω), on average, θ remains constant. This evidences that the flow remains in a steady-state regime, also called the “rolling” flow regime. For every simulation, this steady state is maintained for approximately five revolutions before stopping the simulations. Moreover, Fig. 5(b) shows that ⟨θ⟩ (i.e., the value of θ averaged between revolutions 1 and 5) is an increasing, but nonlinear, function of η that seems to saturate for η > 0.7, for all Ω.

In the following subsections, the data we present are averaged over time in the steady-state regime. These averaged profiles are computed within the Σ zone at the drum center whose thickness is fixed to 40r [see Fig. 4(b)]. This zone is then vertically divided into layers of thickness equal to δz = 1.5r. This choice is motivated by a preliminary sensitivity study of the profiles not presented here. Averages are taken over 480 states saved between the first and the last revolution in the steady-state regime.

## B. Stress profile prediction

In this section we aim to validate the predictions of σzz(z) and σxz(z) given by Eqs. (17) and (16), respectively. To do so, we first need to evaluate the evolution of the granular stress tensor and the packing fraction as a function of z.
<figure><img src="image_7.png" /><figcaption>FIG. 5. (a) Evolution of the free-surface angle $\theta$ as a function of the number of revolutions for a drum speed $\Omega = 2$ rpm ($\Omega = 10$ rpm inset) and different particle shapes. (b) $\langle \theta \rangle$, the value of $\theta$ averaged between revolutions 1 and 5, as a function of $\eta$ after one revolution at drum speeds $\Omega = 2, 6, 10$ rpm. The droplike symbols give (vertically) the histograms of $\theta$ values at drum speeds $\Omega = 2$ and $10$ rpm.</figcaption></figure>

## 1. Definitions

To compute the stress profiles within the assemblies, we first construct the granular stress tensor $\sigma$. For each particle $p$, we calculate the internal moment tensor $M_{ij}^p = \sum_{c \in p} F_i^c r_j^c$ [97], where $F_i^c$ represents the $i$th component of the force applied on particle $p$ at contact $c$, and $r_j^c$ denotes the $j$th component of the position vector of contact $c$. The sum runs over all contacts $c$ of particle $p$. Subsequently, a Voronoi tessellation, covering the volume (or area in 2D) occupied by the grains (see Appendix D), is employed to measure the effective volume $V_v^p$ occupied by particle $p$. Finally, the components of the granular stress tensor $\sigma(z)$ at an altitude $z$ are given by [97]

$$ \sigma^c(z) = \frac{1}{V_z} \sum_{p \in [z, z+\delta z]} M^p, \quad (26) $$

where the sum runs over all the particles $p$, with their center of mass within $[z, z + \delta z]$; and $V_z$ is the sum of the local volumes $V_v^p$ of the corresponding particles. Finally, along the same line, the packing fraction profile $\phi(z)$ is built using the Voronoi tessellation as

$$ \phi(z) = \frac{1}{V_z} \sum_{p \in [z, z+\delta z]} V^p, \quad (27) $$

<figure><img src="image_8.png" /><figcaption>FIG. 6. Packing fraction profile $\phi(z)$ for $\eta = 0.8$ and 0 (inset) for all rotation speed $\Omega$ (color level). The solid lines are linear approximations of the packing fraction in the static regime (vertical line) and in the flow regime (sloping line). The intersection of these two lines, marked with a black disk, corresponds to the height $z_s$, which we use to mark the transition height between the static and the flow regimes inside the drum.</figcaption></figure>

where $V^p$ is the volume (area in 2D) of a particle $p$ whose center of mass belongs to the slice $[z, z + \delta z]$.

## 2. Profiles

Figure 6 shows the packing fraction profiles $\phi(z)$ for $\eta = 0.8$ and 0 (inset) for all rotation speed $\Omega$ simulated in this study. $\phi(z)$ is in good agreement with the description given in Sec. IV C. Indeed, $\phi(z)$ is nearly uniform in the bulk and equal to $\phi_s$ until it declines nearly linearly in the flowing layers, from $\phi_s$ at $z_s$ to $\phi_c$ at $z_c$ and tends towards zero at the free surface, i.e., for $z > z_c$. Yet, close to the drum boundaries, the packing fraction is also found to significantly decrease within a small region one or two grain diameters thick. These profiles also evidence an expansion of the assemblies at large-$\eta$ values: $\phi$ declines as $\eta$ increases. An expansion is also noticeable with $\Omega$ but only in the flowing zone.

Figure 7(a) shows the variation of $z_s$ as a function of $\eta$ for the different rotation speed $\Omega$. This quantity is measured from a linear fit of the packing fraction profiles as explained in Sec. IV C and illustrated in Fig. 3. As a reminder, $z_s$ is a characteristic depth of the system that is the input to the velocity profile equation [Eq. (24)]. $z_s$ evolves nonlinearly with $\eta$: it first increases and, beyond $\eta = 0.4$, remains relatively constant. Meanwhile, $z_s$ declines as $\Omega$ increases but conserves its trend as a function of $\eta$. Figure 7(b) shows the variation of quantities $\phi_s$ and $\bar{\phi}$ as functions of $\eta$ for $\Omega = 2$ and $10$ rpm. We recall that $\phi_s$ is the packing fraction averaged in the solidlike zone, while $\bar{\phi}$ is the packing fraction averaged in the liquid zone between $z_s$ and $z_c$. These quantities coincide regardless of the rotation speed $\Omega$. Interestingly, these packing fractions slightly increase with $\eta$, reach a maximum for $\eta = 0.5$, and then quickly decline as $\eta$ further increases. A similar nonmonotonous behavior has been observed in previous works with elongated [65,108], angular [65,66,75], or slightly nonconvex [66,72,75] grains in both two and three dimensions.

Then, Fig. 8 displays the normal stress $\sigma_{zz}(z)$ (a) and shear stress $\sigma_{xz}(z)$ (b) profiles for different values of $\eta$, and for different loading speed $\Omega = 2$ and $10$ rpm. We observe that
<figure><img src="image_9.png" /><figcaption>FIG. 7. (a) Evolution of $z_s$ (the depth of the boundary between static and quasistatic zones) as a function of the shape parameter $\eta$. (b) Evolution of $\phi_s$ (the packing fraction averaged in the solidlike zone) and $\bar{\phi}$ (the packing fraction averaged in the liquid zone, between $z_s$ and $z_c$) as functions of the shape parameter $\eta$ for $\Omega = 2$ (dashed line) and $\Omega = 10$ (plain line).</figcaption></figure>

$\sigma_{zz}(z)$ decreases linearly with $z$ and decreases with $\eta$, from the top layer, $\sigma_{xz}(z)$ first decreases with depth, and then saturates close to the drum border at a distance between $0.2H$ and $0.3H$ depending on $\eta$. When increasing this latter parameter, the shear stress first increases (in absolute value), but beyond $\eta = 0.5$ it goes back to values slightly lower than that of disk assemblies, for $\eta = 0.9$.

The stress profiles $\sigma_{zz}(z)$ are well approximated by Eq. (17) for all pairs $(\eta, \Omega)$ used in this study. The prediction of $\sigma_{xz}(z)$ given by Eq. (16) also nicely reproduces the main variations of the shear stress profiles across the drum for all pairs $(\eta, \Omega)$, from the free surface up to deep in the solidlike phase. It correctly captures the constant shear stress profile observed close to the border. The minor mismatch with the numerical data mainly comes from the fact we impose $\sigma_{xx} = \sigma_{zz}$ to write Eq. (16). This is valid in bulk, but not necessarily close to the border, consistent with previous works [13,109], discussion is added in Appendix C on this specific aspect.

Finally, following of Eq. (16) the nonmonotonous variation of $\sigma_{xz}$ as a function $\eta$ can be better understood from the combined contributions of the evolution of $\bar{\phi}$ and $\langle\theta\rangle$ with this same parameter $\eta$. Indeed, for small-$\eta$ values both $\sin\langle\theta\rangle$ and $\bar{\phi}$ increase, which explains the increase of $\sigma_{xz}$. For $\eta > 0.5$ the rapid decrease of $\bar{\phi}$ with $\eta$ induces a decrease of $\sigma_{xz}$. As a result, the product $\bar{\phi} \sin\langle\theta\rangle$ increases first with $\eta$ but declines from $\eta = 0.5$ up to values lower than that of disk packings (see Fig. 9). The same holds for the variation of $\sigma_{zz}$

<figure><img src="image_10.png" /><figcaption>FIG. 8. Normal stress profile $\sigma_{zz}$ (a) and the shear stress profile $\sigma_{xz}$ (b) for different values of $\eta$ and for $\Omega = 2$ rpm and $\Omega = 10$ rpm (inset). The simulation data are presented as symbols and the model predictions [Eqs. (17) and (16), respectively] by dashed lines.</figcaption></figure>

with $\eta$. In this latter case, the small increases of $\bar{\phi}$ at small $\eta$ are compensated by the decrease of $\cos\langle\theta\rangle$. On the contrary, at larger $\eta$, both $\bar{\phi}$ and $\cos\langle\theta\rangle$ decline, so that their product declines with $\eta$ (see Fig. 9).

### C. Velocity profile prediction and fluidity

In this section, we test the predictions of velocity profiles given by Eq. (24). We discuss the effect that the grain shape $\eta$ and the rotation speed $\Omega$ have on the thickness $\lambda$ of the inertial flow zone.

<figure><img src="image_11.png" /><figcaption>FIG. 9. Evolution of the slopes $\bar{\phi} \cos\langle\theta\rangle$ (dashed line) and $\bar{\phi} \sin\langle\theta\rangle$ (full line) of the corresponding normal [Eq. (17)] and shear [Eq. (16)] stress profiles, as functions of $\eta$, for $\Omega = 2$ rpm (circle) and $\Omega = 10$ rpm (squares).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIG. 10. Velocity profiles measured in the numerical simulations along the $e_x$ direction at drum speed $\Omega = 2$ rpm (a) and $\Omega = 10$ rpm (b) for different $\eta$. In both cases (a) and (b), the plain line represents the velocity profile of the model [Eq. (24)], and, the blue dashed line represents the velocity profiles of the drum. The simulation data are shown as scatter plots.</figcaption></figure>

Figure 10 displays the $x$-velocity profiles $v_x$, averaged in the slice $\Sigma$, for $\eta \in [0, 0.5, 0.9]$ and at drum speeds $\Omega = 2$ rpm (a) and $\Omega = 10$ rpm (b). Symbols are direct measurements of the numerical simulations while the plain lines are the predictions given by Eq. (24). For every shape and drum speed, a solidlike behavior is displayed deep in the packing. This corresponds with the speed profile following the straight blue dashed lines in Fig. 10. Then, in the upper layers, the speed profiles increase and differ from the drum speed at different depths, depending on the values of $\eta$ and $\Omega$.

The prediction of the relative velocity profile, given by Eq. (24), is shown in black lines in Fig. 10 by fitting the free parameters ($k, \xi, \lambda$). $z_s$ is fixed from the packing fraction profiles curves (see Fig. 6). The theoretical prediction fits the data well for every $\eta$ and $\Omega$. The way the fitting is carried out is first by fitting the “universal” constant $k$ by minimizing the sum of residuals for all the fits, corresponding with all the simulations. We get $k = 10$. Then, the two unknowns $\xi$ and $\lambda$ are solved using the two equations (24) (for the velocity profile) and (20) (for the fluidity profile) simultaneously. Doing so, we find that $\xi = 2$, independently of $\eta$ and $\Omega$ (see also Fig. 17 in Appendix B). On the contrary, as shown in Fig. 11(a), $\lambda$ varies significantly with $\eta$ and in a lower proportion with $\Omega$. More precisely, $\lambda$ first decreases with $\eta$ and then reaches a plateau from $\eta = 0.4$ while it continuously increases with $\Omega$. From $\lambda$ we can then compute the transition depths between the quasistatic and inertial flow $z_0$. Interestingly, as shown in

<figure><img src="image_13.png" /><figcaption>FIG. 11. The scaled characteristic thickness of the inertial zone $\lambda/H$ (a), and of the scaled critical depth $z_0/H$ (b) as functions of the concavity $\eta$ for different drum speeds $\Omega$.</figcaption></figure>

Fig. 11(b), this latter remains independent of $\eta$ but decreases with $\Omega$.

The strength of the model we propose in this paper lies not only in the correct modeling of the velocity, density, and stress profiles for any grain shape and loading speed but also in the ability to predict the thickness of the different flow layers and their depth.

## VI. EXPERIMENTAL VALIDATION

In this section, we develop a series of grain flow experiments in a rotating drum. We used molded rigid cross-shaped 3D grains with different levels of convexity to test our theoretical velocity model experimentally.

### A. Experimental setup

The rotating drum is shown in Fig. 12(a). It is a homemade device that has an inside diameter of 28 cm and a lateral thickness of 5 cm. It is filled with monodisperse particles of shape varied from spherical to highly concave as presented in Fig. 12(b). All the particles are circumscribed in a sphere of diameter $d = 12$ mm. The concave ones consist of three spherocylinders that extend toward the direction of the faces of a regular cube. The radius of the spherocylinders $r_0$ gradually increases from 0.75 to 6 mm (sphere). These particles are made by injection molding of high-density polyethylene (HDPE). A dedicated mold producing clusters of 20 particles has been custom made [110]. This material's friction coefficient is very low, close to 0.1, and its Young modulus is quite high, around 1 GPa. This makes the particles rigid and slippery. The concavity parameter $\eta$ [see Eq. (25)] is varied such that $\eta \in [0, 0.33, 0.5, 0.58, 0.67, 0.71, 0.75, 0.79, 0.83, 0.875]$.

The drum has smooth transparent glass on both axial sides, and the radial inner surface has regularly raised band to prevent particle sliding. A constant volume of 900 mL of particles is loaded into the drum for each experiment. The drum is illuminated with LED light, and a camera [111] is positioned perpendicular to the front glass to image the system at a high
<figure><img src="image_14.png" /><figcaption>FIG. 12. (a) Illustration of the experimental setup: A camera, B drum, C stepper motors, D lighting. (b) Particles with varying shapes, ranging from spherical ($\eta = 0\%$) to highly concave ($\eta = 87.5\%$). When the concavity $\eta$ increases, the spherocylinder diameter gradually decreases from 12 to 1.5 mm. (c) Velocity field obtained from particle image velocimetry (PIV) superimposed on a photo of the flowing grains for particles with concavity $\eta = 79\%$. The dark blue dotted line represents the studied slice $\Sigma$, which has a width of six particle diameters (7.2 cm).</figcaption></figure>

frequency (60 frames per second). The rotation speed of the drum is controlled by tuning the speed of two synchronized stepper motors below it [see Fig. 12(a)]. We tested three angular velocities $\Omega \in [1.93, 2.91, 4.83]$ rpm. For each couple ($\eta, \Omega$) the drum is first rolled for 2 min to reach the steady state, while image recording lasted 5 min.

The mean behavior for each set of parameters $\eta$ and $\Omega$ is obtained by averaging over three independent data sets. Each experiment was repeated 3 times, so in total 33 experiments were performed. Image analysis shows that the flows stay in the rolling regime for these pairs of parameters: ($\eta, \Omega) \in \{([0, \dots, 0.58], 1.93) ; ([0, \dots, 0.67], 2.91) ; ([0, \dots, 0.79], 4.83)\}$. In the following, only these data sets are used. Figure 13 shows $\langle \theta \rangle$ as a function of $\eta$ for the three $\Omega$ rotation speeds. Obviously, the experimental and numerical values [see Fig. 5(b)] of $\langle \theta \rangle$ cannot be directly compared since these are 2D simulations vs 3D experiments with different grain shapes. However, in both cases the measured $\langle \theta \rangle$ is an increasing function of both $\eta$ and $\Omega$.

<figure><img src="image_15.png" /><figcaption>FIG. 13. The time-average free-surface angle $\langle \theta \rangle$ as a function of $\eta$ and for different drum speed $\Omega \in [1.93, 2.91, 4.83]$ rpm. The vertical histogram represents the distribution of the surface angles for a rotational speed $\Omega = 4.83$ rpm.</figcaption></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>FIG. 14. Velocity profiles measured experimentally along the $e_x$ direction at drum speed $\Omega = 1.93$ rpm (a) and $\Omega = 4.83$ rpm (b) for different $\eta$. In both cases (a) and (b), the black line represents the velocity profile of the model [Eq. (24)], and the blue dashed line represents the velocity profiles of the drum. The experimental data are shown as scatter plots.</figcaption></figure>

## B. Velocity profile prediction

From the images, we deduce the instantaneous velocity field $v$ using particle image velocimetry (PIV) as shown in Fig. 12(c). Then, we decompose this velocity field into directions parallel ($e_x$) and perpendicular ($e_z$) to the free surface. The so obtained $x$-velocity component $v_x$ is averaged in time within the slice $\Sigma$ to obtain the velocity profiles $v_x(z)$ (see Fig. 14). The blue dashed line shows the drum velocity as a function of the depth. The black line is the velocity profile fitted by the analytical model [see Eq. (24)].

The velocity profiles obtained experimentally follow the same trends as those obtained numerically. From experiments, we cannot measure the evolution of the packing fraction. So we cannot determine the parameter $z_s$ a priori, and Eq. (20) cannot be used based on the evolution of the solid fraction. Consequently, the parameter $\xi$ cannot be estimated precisely. Nevertheless, remembering that $\lambda = z_0 - z_s = \xi\sqrt{8l}$, Eq. (24) can be rewritten as

$$ \frac{\partial v_x(z)}{\partial z} = \frac{1}{2kd} \sqrt{g \frac{H^2 - z^2}{2H} \sin\langle\theta\rangle \left[1 + \tanh\left(\frac{z - z_0}{\sqrt{8l}}\right)\right]} - \Omega. \quad (28) $$

In this new form, Eq. (28) is based on three parameters, i.e., $k$, $z_0$, and $l$ that have to be determined. The same strategy as previously is used to determine these three parameters. First, $k$ is fixed to 4 by minimizing the sum of all residuals for all fits. Once $k$ is determined, the two last parameters are fitted for each experiment. As it is seen in Fig. 14, the prediction of
<figure><img src="image_18.png" /><figcaption>FIG. 15. The scaled characteristic thickness of the inertial zone $\lambda/H$ (a) and the scaled critical depth $z_0/H$ (b) as functions of the concavity $\eta$ for different drum speeds $\Omega$ in the experiments.</figcaption></figure>

the relative velocity profile, given by Eq. (28) (lines) by fitting the free parameters ($k, z_0, l$) is good for all $\eta$ and $\Omega$.

### C. Discussion and comparison with numerical results

In order to compare the 2D and 3D values of the fitting parameters, Fig. 15 displays the variations of the scaled characteristic thickness of the inertial zone $\lambda/H$ (a), and of the scaled critical depth $z_0/H$ (b) in the experiments. Both are plotted as functions of the concavity $\eta$ for different drum speeds $\Omega$. Two comments: (1) to calculate $\lambda$ it is necessary to have an estimate of the parameter $\xi$. As explained in Sec. IV C, $\xi$ is a purely numerical parameter used to ensure the equality of Eq. (20). In this sense, it does not depend on the dimensionality of the problem and consists simply in studying the convergence of the function $\tanh \xi$ to 1. For this reason, we can rely on the 2D results and assume $\xi = 2$ as a first approximation. (2) To support this choice, an inverse analysis can be performed, consisting of taking Eq. (24) and imposing $\xi = 2$, which allows to find values for $\lambda$ and $z_0$ similar to the previous method, giving confidence in the measurement of the fitting parameters.

Thus, consistently with the numerical results, the critical depth $z_0$ decreases when the rotating speed $\Omega$ increases. Yet, for a given $\Omega$, it is almost constant as long as $\eta$ is low enough, and increases (faster as in the 2D numerical simulations) with this latter parameter for values above 0.5. The thickness of the inertial flow zone $\lambda/H$ is quite independent of the rotation speed $\Omega$, as already observed numerically. On the other hand, its variation as a function of $\eta$ in the experiments slightly differs from the simulations. Specifically, $\lambda/H$ decreases less steeply than in the numerical simulations. We also observe an increase for the highest $\eta$ values, whereas it tends to plateau in the simulations. Finally, in addition to the very similar trends observed for these parameters as a function of $\eta$ and $\Omega$, we note that similar (or at least very close) values are obtained

for the $\eta = 0$ case. For example, for $\Omega = 2$ rpm, $\lambda/H$ is close to 0.3 in the 3D experiments and close to 0.35 in the 2D simulations. The same holds for $z_0/H$ close to 0.90 in both approaches.

### VII. CONCLUSION

In this paper, we have presented an analytical model to describe the evolution of the steady granular flow in a rotating drum. Our approach relates the shear stress $\tau$ to the angle of repose of the free surface $\theta$ and introduces the shear rate $\dot{\gamma}$ and packing fraction $\phi$ into the generalized Bagnold equation through the “fluidity” $f$. We have found $f$ to depend only on the packing fraction $\phi$.

To test our model, we designed a series of 2D numerical simulations (based on a discrete-element method) and 3D experiments involving grain flow in a rotating drum. In addition, to test the robustness of the theoretical model, we went as far as considering the case of granular systems made of very nonconvex grains. Grains with a very nonconvex shape have the peculiarity of being entangled and thus present flow patterns that can be radically different from those of convex grains.

Our numerical and experimental results confirm the validity of the model for steady granular flows and its ability to accurately describe the velocity profile, regardless of the grain shapes. Furthermore, our simulations confirm the model’s efficacy in predicting stress and packing fraction within the assembly. The combination of experiments and simulations compared to the model demonstrates the power of an approach introducing the concept of “fluidity” in the analytical study. Additionally, we provide an explicit form for this fluidity function $[F(\phi(z))]$ for cross validation.

In our analytical model, the “fluidity” can be viewed as an approximate solution to the diffusion equation based on Landau’s theory. To avoid the artificial building of a source term in the diffusion equation, we introduce a constant $k$ into the fluidity. This fluidity is then incorporated into the Bagnold scaling with the Prandtl mixing length scale to correct the mixing length. The successful combination of Landau theory and Prandtl mixing length in a macroscopic system primarily governed by inertia provides us with a more intuitive understanding of phase transitions and fluid mechanisms. This hyperbolic logistic function appears in many statistical problems related to phase transitions and provides a highly approximate solution to bridge the different states of the order parameter. It is widely used not only in physical phase transitions but also in ecological [112], epidemiologic [113], and chemical phase transitions. We apply a hyperbolic tangent function to the order parameter fluidity in granular flow. This approach may have applications in similar rheological systems of amorphous glassy materials, such as concentrated emulsions, pastes, and molecular glasses.

In addition, our analytical framework incorporates most of the key physical quantities in particle flow. In particular, we evidenced a critical height $z_s$, which marks the transition between the solid and fluid regimes, and a critical thickness $\lambda$ characterizing the inertial flow zone. Our approach, through the concept of fluidity, highlights a characteristic length scale specific to the grain shape $l$, which is explicitly linked to the
critical thickness by the relationship we establish $\lambda = \xi\sqrt{8l}$. In other words, our model explicitly incorporates the grain shape through this intrinsic grain length scale $l$.

In this paper, we have focused on the modeling of the macroscopic flow. However, a lot of work remains to be done to characterize the microscopic properties, particularly in the presence of very nonconvex grains. As mentioned above, these grains can become entangled and most likely form complex clusters connected by multiple contacts. Such a description remains an open research topic and will be the subject of a detailed study in forthcoming publications.

Finally, it would be interesting to reexamine the results presented in this paper with other grain shapes, in particular anisometrically elongated grain shapes, and different fill rates. In general, elongated grain shapes, such as spheroids and spherocylinders, tend to develop an orientational order that affects force transmission and friction behavior [82]. This “nematic” ordering occurs even when the grains interact only by contact and friction. It is hoped that the theory developed in this paper for predicting stress and velocity profiles will remain valid for more simple and convex grain shapes as well as for different filling. It would then be interesting to quantify the thickness of phase changes (solid to liquid) as a function of both grain elongation and filling, in order to characterize the underlying microstructures in each of the phases involved.

## ACKNOWLEDGMENTS

The authors acknowledge financial support from Agence nationale de recherche MICROGRAM (Grant No. ANR-20-CE92-0009). We also acknowledge the support of the High-Performance Computing Platform MESO@LR. Finally, we would like to thank G. Camp, G. Genevois, and Q. Chapelier for their technical help in setting up the experimental device and B. Gallard and S. Buonomo for their help in making particles.

## APPENDIX A: MACROSCOPIC FRICTION PROFILE

Figure 16 shows the variation of the macroscopic friction profile $\mu(z)$ as a function of $z$. As we can see, $\mu(z)$ evolves approximately linearly with $z$ in our numerical simulations (as also reported in [36]).

## APPENDIX B: NUMERICAL APPROXIMATION FOR SOLVING EQ. (20)

Figure 17 shows the evolution of $f(z)$ (lines) and $F(\phi(z))$ (symbols) for $\eta = 0$ (a) and $\eta = 0.9$ rpm (b). The data always collapse when $\xi$ is chosen equal to $\xi = 2$ and $\lambda$ given in Fig. 11(a).

## APPENDIX C: REWRITING EQ. (16) WITHOUT NEGLECTING EXTRA TERMS

When developing Eq. (16) we assumed that the stresses $\sigma_{xx}$ and $\sigma_{zz}$ were equal. This assumption is valid in the bulk but may be too strong near the walls, more precisely at the bottom of the drum. Although the approximation proposed

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>FIG. 16. Effective friction coefficient $\mu$. (a) $\mu$ is plotted for different values of $\eta$ and for $\Omega = 2$ rpm. (b) $\mu$ is plotted for different values of $\eta$ and for $\Omega = 10$ rpm.</figcaption></figure>

by Eq. (16) reproduces the shear stress profiles very well, for all shapes and drum angular speed, it can be significantly improved without the previous assumption. In this case, the

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /><figcaption>FIG. 17. Fluidity profile $f(z)$ (plain line) and the specific form of fluidity profile as a function of the packing fraction $F(\phi(z))$ (scatter plot) for $\eta = 0$ (a) and $\eta = 0.9$ (b).</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>FIG. 18. Corrected shear stress, for different values of η and for Ω = 2 rpm (a) and for different values of η and for Ω = 10 rpm (b).</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>FIG. 19. Voronoi tessellation in a snapshot of simulation. (a) For concave particles with η = 0.5. (b) For highly concave particles with η = 0.9.</figcaption></figure>

shear stress model reads as

$$ \frac{\tau(z)}{\rho_0 g d} = \phi(z) \frac{H^2 - z^2}{2Hd} \sin\langle\theta\rangle + \frac{\sqrt{|\sigma_{xx}^2(z) - \sigma_{zz}^2(z)|}}{\sigma_{zz}(z)}. \quad (C1) $$

Figure 18 shows the evolution of the normalized shear stress profile τ for all shapes at Ω = 2 rpm (a) and Ω = 10 rpm (b). It also shows the approximation proposed by Eq. (C1), where σxx(z) and σzz(z) are measured directly from the simulations. As we can see, the approximation is more accurate than that given by Eq. (16) in the region close to the drum border.

## APPENDIX D: VORONOI TESSELLATION

To compute the packing fraction and stress, we used Voronoi tessellation to partition and determine the packing volume of each particle. Figure 19 shows the Voronoi meshing for highly nonconvex grains.

[1] H. M. Jaeger, S. R. Nagel, and R. P. Behringer, Granular solids, liquids, and gases, Rev. Mod. Phys. 68, 1259 (1996).

[2] G. D. R. MiDi, On dense granular flows, Eur. Phys. J. E 14, 341 (2004).

[3] Y. Forterre and O. Pouliquen, Flows of dense granular media, Annu. Rev. Fluid Mech. 40, 1 (2008).

[4] J. Goyon, A. Colin, G. Ovarlez, A. Ajdari, and L. Bocquet, Spatial cooperativity in soft glassy flows, Nature (London) 454, 84 (2008).

[5] K. Kamrin and G. Koval, Nonlocal constitutive relation for steady granular flow, Phys. Rev. Lett. 108, 178301 (2012).

[6] F. Tapia, O. Pouliquen, and É. Guazzelli, Influence of surface roughness on the rheology of immersed and dry frictional spheres, Phys. Rev. Fluids 4, 104302 (2019).

[7] Y. L. Ding, J. P. K. Seville, R. Forster, and D. J. Parker, Solids motion in rolling mode rotating drums operated at low to medium rotational speeds, Chem. Eng. Sci. 56, 1769 (2001).

[8] D. Bonamy, F. Daviaud, and L. Laurent, Experimental study of granular surface flows via a fast camera: A continuous description, Phys. Fluids 14, 1666 (2002).

[9] N. Jain, J. M. Ottino, and R. M. Lueptow, An experimental study of the flowing granular layer in a rotating tumbler, Phys. Fluids 14, 572 (2002).

[10] J. Rajchenbach, Granular flows, Adv. Phys. 49, 229 (2000).

[11] J.-N. Roux and G. Combe, Quasistatic rheology and the origins of strain, C. R. Phys. 3, 131 (2002).
[12] A. V. Orpe and D. V. Khakhar, Scaling relations for granular flow in quasi-two-dimensional rotating cylinders, Phys. Rev. E **64**, 031302 (2001).

[13] M. Renouf, D. Bonamy, F. Dubois, and P. Alart, Numerical simulation of two-dimensional steady granular flows in rotating drum: On surface flow rheology, Phys. Fluids **17**, 103303 (2005).

[14] A. V. Orpe and D. V. Khakhar, Rheology of surface granular flows, J. Fluid Mech. **571**, 1 (2007).

[15] I. Govender, Granular flows in rotating drums: A rheological perspective, Minerals Eng. **92**, 168 (2016).

[16] P. Jop, Y. Forterre, and O. Pouliquen, Crucial role of sidewalls in granular surface flows: consequences for the rheology, J. Fluid Mech. **541**, 167 (2005).

[17] P. Jop, Y. Forterre, and O. Pouliquen, A constitutive law for dense granular flows, Nature (London) **441**, 727 (2006).

[18] G. Félix, V. Falk, and U. d'Ortona, Granular flows in a rotating drum: the scaling law between velocity and thickness of the flow, Eur. Phys. J. E **22**, 25 (2007).

[19] J. Rajchenbach, Flow in powders: From discrete avalanches to continuous regime, Phys. Rev. Lett. **65**, 2221 (1990).

[20] R. Khosropour, E. Valachovic, and B. Lincoln, Flow and pattern formation in a binary mixture of rotating granular materials, Phys. Rev. E **62**, 807 (2000).

[21] N. Taberlet, P. Richard, and E. J. Hinch, S shape of a granular pile in a rotating drum, Phys. Rev. E **73**, 050301(R) (2006).

[22] F. Pignatel, C. Asselin, L. Krieger, I. C. Christov, J. M. Ottino, and R. M. Lueptow, Parameters and scalings for dry and immersed granular flowing layers in rotating tumblers, Phys. Rev. E **86**, 011304 (2012).

[23] L. F. Orozco, J.-Y. Delenne, P. Sornay, and F. Radjai, Rheology and scaling behavior of cascading granular flows in rotating drums, J. Rheology **64**, 915 (2020).

[24] F. da Cruz, S. Emam, M. Prochnow, J.-N. Roux, and F. Chevoir, Rheophysics of dense granular materials: Discrete simulation of plane shear flows, Phys. Rev. E **72**, 021309 (2005).

[25] T. Hatano, Power-law friction in closely packed granular materials, Phys. Rev. E **75**, 060301(R) (2007).

[26] C. H. Rycroft, K. Kamrin, and M. Z. Bazant, Assessing continuum postulates in simulations of granular flow, J. Mech. Phys. Solids **57**, 828 (2009).

[27] D. Bonamy, P.-H. Chavanis, P.-P. Cortet, F. Daviaud, B. Dubrulle, and M. Renouf, Euler-like modelling of dense granular flows: application to a rotating drum, Eur. Phys. J. B **68**, 619 (2009).

[28] C. S. Campbell, Rapid granular flows, Annu. Rev. Fluid Mech. **22**, 57 (1990).

[29] Y. Forterre and O. Pouliquen, Physics of particulate flows: From sand avalanche to active suspensions in plants, C. R. Phys. **19**, 271 (2018).

[30] G. Lois, A. Lemaître, and J. M. Carlson, Numerical tests of constitutive laws for dense granular flows, Phys. Rev. E **72**, 051303 (2005).

[31] D. Ertaş and T. C. Halsey, Granular gravitational collapse and chute flow, Europhys. Lett. **60**, 931 (2002).

[32] T. Miller, P. Rognon, B. Metzger, and I. Einav, Eddy viscosity in dense granular flows, Phys. Rev. Lett. **111**, 058002 (2013).

[33] P. G. Rognon, T. Miller, B. Metzger, and I. Einav, Long-range wall perturbations in dense granular flows, J. Fluid Mech. **764**, 171 (2015).

[34] L. E. Silbert, D. Ertaş, G. S. Grest, T. C. Halsey, D. Levine, and S. J. Plimpton, Granular flow down an inclined plane: Bagnold scaling and rheology, Phys. Rev. E **64**, 051302 (2001).

[35] P.-P. Cortet, D. Bonamy, F. Daviaud, O. Dauchot, B. Dubrulle, and M. Renouf, Relevance of visco-plastic theory in a multi-directional inhomogeneous granular flow, Europhys. Lett. **88**, 14001 (2009).

[36] C.-C. Lin, M.-Z. Jiang, and F.-L. Yang, Developing a rheological relation for transient dense granular flows via discrete element simulation in a rotating drum, J. Mech. **36**, 707 (2020).

[37] T. M. Povall, I. Govender, and A. T. McBride, Dense granular flows in rotating drums: a computational investigation of constitutive equations, Powder Technol. **393**, 238 (2021).

[38] T. Barker, D. G. Schaeffer, P. Bohorquez, and J. M. N. T. Gray, Well-posed and ill-posed behaviour of the $\mu(i)$-rheology for granular flow, J. Fluid Mech. **779**, 794 (2015).

[39] T. Barker and J. M. N. T. Gray, Partial regularisation of the incompressible $\mu(I)$-rheology for granular flow, J. Fluid Mech. **828**, 5 (2017).

[40] T. Barker, M. Rauter, E. S. F. Maguire, C. G. Johnson, and J. M. N. T. Gray, Coupling rheology and segregation in granular flows, J. Fluid Mech. **909**, A22 (2021).

[41] E. S. F. Maguire, T. Barker, M. Rauter, C. G. Johnson, and J. M. N. T. Gray, Particle-size segregation patterns in a partially filled triangular rotating drum, J. Fluid Mech. **979**, A40 (2024).

[42] T. Barker, D. G. Schaeffer, M. Shearer, and J. M. N. T. Gray, Well-posed continuum equations for granular flow with compressibility and $\mu(I)$-rheology, Proc. R. Soc. London Ser. A **473**, 20160846 (2017).

[43] J. Heyman, R. Delannay, H. Tabuteau, and A. Valance, Compressibility regularizes the $\mu(I)$-rheology for dense granular flows, J. Fluid Mech. **830**, 553 (2017).

[44] D. G. Schaeffer, T. Barker, D. Tsuji, P. Gremaud, M. Shearer, and J. M. Gray, Constitutive relations for compressible granular flow in the inertial regime, J. Fluid Mech. **874**, 926 (2019).

[45] L. Prandtl, 7. bericht über untersuchungen zur ausgebildeten turbulenz, Z. Angew. Math. Mech. **5**, 136 (1925).

[46] E. DeGiuli, G. Düring, E. Lerner, and M. Wyart, Unified theory of inertial granular flows and non-brownian suspensions, Phys. Rev. E **91**, 062206 (2015).

[47] K. Saitoh and T. Kawasaki, Critical scaling of diffusion coefficients and size of rigid clusters of soft athermal particles under shear, Front. Phys. **8**, 99 (2020).

[48] M. Bouzid, A. Izzet, M. Trulsson, E. Clément, P. Claudin, and B. Andreotti, Non-local rheology in dense granular flows: Revisiting the concept of fluidity, Eur. Phys. J. E **38**, 125 (2015).

[49] K. Kamrin, Non-locality in granular flow: Phenomenology and modeling approaches, Front. Phys. **7**, 116 (2019).

[50] I. S. Aranson and L. S. Tsimring, Continuum theory of partially fluidized granular flows, Phys. Rev. E **65**, 061303 (2002).

[51] M. Bouzid, M. Trulsson, P. Claudin, E. Clément, and B. Andreotti, Nonlocal rheology of granular flows across yield conditions, Phys. Rev. Lett. **111**, 238301 (2013).
[52] K. Kamrin and D. L. Henann, Nonlocal modeling of granular flows down inclines, *Soft Matter* **11**, 179 (2015).

[53] I. S. Aranson and L. S. Tsimring, Continuum description of avalanches in granular media, *Phys. Rev. E* **64**, 020301(R) (2001).

[54] D. Wolfson, L. S. Tsimring, and I. S. Aranson, Order parameter description of stationary partially fluidized shear granular flows, *Phys. Rev. Lett.* **90**, 254301 (2003).

[55] I. S. Aranson and L. S. Tsimring, Patterns and collective behavior in granular media: Theoretical concepts, *Rev. Mod. Phys.* **78**, 641 (2006).

[56] F. Malloggi, J. Lanuza, B. Andreotti, and E. Clément, Erosion waves: Transverse instabilities and fingering, *Europhys. Lett.* **75**, 825 (2006).

[57] D. Liu and D. L. Henann, Size-dependence of the flow threshold in dense granular materials, *Soft Matter* **14**, 5294 (2018).

[58] L. Bocquet, A. Colin, and A. Ajdari, Kinetic theory of plastic flow in soft glassy materials, *Phys. Rev. Lett.* **103**, 036001 (2009).

[59] Q. Zhang and K. Kamrin, Microscopic description of the granular fluidity field in nonlocal flow modeling, *Phys. Rev. Lett.* **118**, 058001 (2017).

[60] F. Fazelpour, Z. Tang, and K. E. Daniels, The effect of grain shape and material on the nonlocal rheology of dense granular flows, *Soft Matter* **18**, 1435 (2022).

[61] R. N. Poon, A. L. Thomas, and N. M. Vriend, Microscopic origin of granular fluidity: An experimental investigation, *Phys. Rev. E* **108**, 064902 (2023).

[62] D. Bonamy and P. Mills, Diphasic non-local model for granular surface flows, *Europhys. Lett.* **63**, 42 (2003).

[63] D. Bonamy, Phénomenes collectifs dans les matériaux granulaires. ecolements de surface et réarrangements internes dans des empilements modeles, Ph.D. thesis, Université Paris Sud - Paris XI, 2001.

[64] P. W. Cleary, The effect of particle shape on simple shear flows, *Powder Technol.* **179**, 144 (2008).

[65] B. Saint-Cyr, K. Szarf, C. Voivret, E. Azéma, V. Richefeu, J.-Y. Delenne, G. Combe, C. Nouguier-Lehon, P. Villard, P. Sornay et al., Particle shape dependence in 2D granular media, *Europhys. Lett.* **98**, 44008 (2012).

[66] A. G. Athanassiadis, M. Z. Miskin, P. Kaplan, N. Rodenberg, S. H. Lee, J. Merritt, E. Brown, J. Amend, H. Lipson, and H. M. Jaeger, Particle shape effects on the stress response of granular packings, *Soft Matter* **10**, 48 (2014).

[67] N. Gravish and D. I. Goldman, *Entangled Granular Media* (Wiley, Hoboken, NJ, 2016), pp. 341–354.

[68] A. D. Rakotonirina, J.-Y. Delenne, F. Radjai, and A. Wachs, Grains3D, a flexible dem approach for particles of arbitrary convex shape part III: extension to non-convex particles modelled as glued convex particles, *Comput. Part. Mech.* **6**, 55 (2019).

[69] T. A. Marschall and S. Teitel, Athermal shearing of frictionless cross-shaped particles of varying aspect ratio, *Granular Matter* **22**, 4 (2020).

[70] M. Cox, D. Wang, J. Barés, and R. P. Behringer, Self-organized magnetic particles to tune the mechanical behavior of a granular system, *Europhys. Lett.* **115**, 64003 (2016).

[71] Y. Zhao, J. Barés, H. Zheng, C. Bester, Y. Xu, J. Socolar, and R. P. Behringer, Jamming transition in non-spherical

particle systems: pentagons versus disks, *Granular Matter* **21**, 90 (2019).

[72] E. Azéma, N. Estrada, and F. Radjai, Nonlinear effects of particle shape angularity in sheared granular media, *Phys. Rev. E* **86**, 041301 (2012).

[73] E. Azéma, F. Radjai, and F. Dubois, Packings of irregular polyhedral particles: Strength, structure, and effects of angularity, *Phys. Rev. E* **87**, 062203 (2013).

[74] T. Binaree, E. Azéma, N. Estrada, M. Renouf, and I. Preechawuttipong, Combined effects of contact friction and particle shape on strength properties and microstructure of sheared granular media, *Phys. Rev. E* **102**, 022901 (2020).

[75] E. Azéma, F. Radjai, B. Saint-Cyr, J.-Y. Delenne, and P. Sornay, Rheology of three-dimensional packings of aggregates: Microstructure and effects of nonconvexity, *Phys. Rev. E* **87**, 052205 (2013).

[76] N. A. Conzelmann, M. N. Partl, F. J. Clemens, C. R. Müller, and L. D. Poulikakos, Effect of artificial aggregate shapes on the porosity, tortuosity and permeability of their packings, *Powder Technol.* **397**, 117019 (2022).

[77] E. Azéma, Y. Descantes, N. Roquet, J.-N. Roux, and F. Chevoir, Discrete simulation of dense flows of polyhedral grains down a rough inclined plane, *Phys. Rev. E* **86**, 031303 (2012).

[78] W. Han, D. Wang, and Y. Wu, The study on flow properties of nonconvex particles down a rough inclined plane, *Phys. Fluids* **33**, 063314 (2021).

[79] J. Tang and R. P. Behringer, Orientation, flow, and clogging in a two-dimensional hopper: Ellipses vs. disks, *Europhys. Lett.* **114**, 34002 (2016).

[80] E. Goldberg, C. M. Carlevaro, and L. A. Pugnaloni, Clogging in two-dimensions: effect of particle shape, *J. Stat. Mech.: Theory Exp.* (2018) 113201.

[81] R. C. Hidalgo, I. Zuriguel, D. Maza, and I. Pagonabarraga, Role of particle shape on the stress propagation in granular packings, *Phys. Rev. Lett.* **103**, 118001 (2009).

[82] E. Azéma and F. Radjai, Force chains and contact network topology in sheared packings of elongated particles, *Phys. Rev. E* **85**, 031303 (2012).

[83] D. Höhner, S. Wirtz, and V. Scherer, A study on the influence of particle shape and shape approximation on particle mechanics in a rotating drum using the discrete element method, *Powder Technol.* **253**, 256 (2014).

[84] H. Ma and Y. Zhao, Investigating the flow of rod-like particles in a horizontal rotating drum using dem simulation, *Granular Matter* **20**, 41 (2018).

[85] O. Dubé, E. Alizadeh, J. Chaouki, and F. Bertrand, Dynamics of non-spherical particles in a rotating drum, *Chem. Eng. Sci.* **101**, 486 (2013).

[86] Y. Zhao, K. Liu, M. Zheng, J. Barés, K. Dierichs, A. Menges, and R. P. Behringer, Packings of 3D stars: stability and structure, *Granular Matter* **18**, 24 (2016).

[87] J. Barés, Y. Zhao, M. Renouf, K. Dierichs, and R. P. Behringer, Structure of hexapod 3D packings: Understanding the global stability from the local organization, *EPJ Web of Conferences* (EDP Sciences, Les Ulis, France, 2017), Vol. 140, p. 06021.

[88] K. Dierichs and A. Menges, Towards an aggregate architecture: designed granular systems as programmable matter in architecture, *Granular Matter* **18**, 25 (2016).
[89] K. Murphy, L. Roth, D. Peterman, and H. Jaeger, Aleatory construction based on jamming: Stability through self-confinement, *Architect. Design* **87**, 74 (2017).

[90] D. Aponte, N. Estrada, J. Barés, M. Renouf, and E. Azéma, Geometric cohesion in two-dimensional systems composed of star-shaped particles, *Phys. Rev. E* **109**, 044908 (2024).

[91] M. Mohammadi, D. Puzyrev, T. Trittel, and R. Stannarius, Secondary flow in ensembles of nonconvex granular particles under shear, *Phys. Rev. E* **106**, L052901 (2022).

[92] D. J. Acheson, *Elementary Fluid Dynamics* (Clarendon, Oxford, 1991).

[93] H. A. Lorentz, A. Einstein, H. Minkowski, H. Weyl, and A. Sommerfeld, *The Principle of Relativity: A Collection of Original Memoirs on the Special and General Theory of Relativity*, Dover Books on Physics and Mathematical Physics (Dover, New York, 1952).

[94] O. Pouliquen and F. Chevoir, Dense flows of dry granular material, *C. R. Phys.* **3**, 163 (2002).

[95] C.-Y. Hung, C. P. Stark, and H. Capart, Granular flow regimes in rotating drums from depth-integrated theory, *Phys. Rev. E* **93**, 030902(R) (2016).

[96] N. S. Ray and D. V. Khakhar, Dem simulations of quasi-two-dimensional flow of spherical particles on a heap without sidewalls, *EPJ Web of Conferences* (EDP Sciences, Les Ulis, France, 2021), Vol. 249, p. 03034.

[97] B. Andreotti, Y. Forterre, and O. Pouliquen, *Granular Media: Between Fluid and Solid* (Cambridge University Press, Cambridge, 2013).

[98] O. Pouliquen and Y. Forterre, Friction law for dense granular flows: application to the motion of a mass down a rough inclined plane, *J. Fluid Mech.* **453**, 133 (2002).

[99] E. C. P. Breard, J. Dufek, L. Fullard, and A. Carrara, The basal friction coefficient of granular flows with and without excess pore pressure: implications for pyroclastic density currents, water-rich debris flows, and rock and submarine avalanches, *J. Geophys. Res.: Solid Earth* **125**, e2020JB020203 (2020).

[100] R. M. Nedderman, *Statics and Kinematics of Granular Materials* (Cambridge University Press, Cambridge, 1992).

[101] R. Artoni, A. Santomaso, and P. Canu, Shear bands in granular flow through a mixing-length model, *Europhys. Lett.* **80**, 34004 (2007).

[102] R. Artoni, A. Santomaso, and P. Canu, Coupling between countercurrent gas and solid flows in a moving granular bed: The role of shear bands at the walls, *Int. J. Multiphase Flow* **37**, 1209 (2011).

[103] J. J. Moreau, Numerical aspects of the sweeping process, *Comput. Methods Appl. Mech. Eng.* **177**, 329 (1999).

[104] P. A. Cundall and O. D. L. Strack, A discrete numerical model for granular assemblies, *Géotechnique* **29**, 47 (1979).

[105] M. Renouf, F. Dubois, and P. Alart, A parallel version of the non smooth contact dynamics algorithm applied to the simulation of granular media, *J. Comput. Appl. Math.* **168**, 375 (2004).

[106] F. Dubois and R. Mozul, *LMGC90, 11th national structural calculation conference* (2013), https://hal.science/hal-01717115.

[107] P. Sánchez, M. Renouf, E. Azéma, R. Mozul, and F. Dubois, A contact dynamics code implementation for the simulation of asteroid evolution and regolith in the asteroid environment, *Icarus* **363**, 114441 (2021).

[108] D. B. Nagy, P. Claudin, T. Börzszönyi, and E. Somfai, Rheology of dense granular flows for elongated particles, *Phys. Rev. E* **96**, 062903 (2017).

[109] T. M. Povall, Dense granular flow in rotating drums: a computational investigation of constitutive equations, Ph.D. thesis, University of Cape Town, 2018.

[110] The mold has been made by cfo-outillage.fr. Injections have been done in IMT Mines Alès on a Krauss Maffei KM50/180CX. Batches of particles can be available upon request to the authors.

[111] USB3 camera Basler acA3088-57um, 6.4 Mpx, 60 frames per second.

[112] R. Pearl and L. J. Reed, On the rate of growth of the population of the united states since 1790 and its mathematical representation, *Proc. Natl. Acad. Sci. USA* **6**, 275 (1920).

[113] S. Y. Lee, B. Lei, and B. Mallick, Estimation of covid-19 spread curves integrating global data and borrowing information, *PLoS One* **15**, e0236860 (2020).