3.1. The VOF method coupled with $\mu (I,\varTheta )$ rheology

Granular flows in both natural and industrial contexts often involve complex free-surface dynamics, such as column collapse, heap formation and landslide runout, where interfaces undergo large deformations, fragmentation, and topological transitions. To resolve these processes within a continuum solver, the VOF method (Versteeg Reference Versteeg2007) is employed. In this approach, the granular–air interface is described by a scalar volume-fraction field $\alpha$ , advected with the local velocity field and discretised within the FVM. This ensures strict phase conservation, accurate reconstruction of interfaces, and robustness in handling topological transitions, thereby providing a consistent foundation for coupling free-surface evolution with rheological closures.

The governing equations of the VOF framework comprise the continuity equation, the transport equation for $\alpha$ , and the momentum balance. Together, they provide a complete description of multiphase dynamics. For dense granular flows at low to moderate velocities, incompressibility is assumed, yielding the continuity condition

(3.1) \begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0, \end{equation}

where $\boldsymbol{u}$ denotes the velocity vector field. The volume-fraction transport equation is expressed as

(3.2) \begin{equation} \frac {\partial \alpha }{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }\left (\boldsymbol{u}\alpha \right )+\boldsymbol{\nabla }\boldsymbol{\cdot }\left [\boldsymbol{u}_r\alpha \left (1-\alpha \right )\right ]=0, \end{equation}

where $\alpha \in [0,1]$ represents the granular-phase volume fraction. Specifically, $\alpha =1$ indicates that the cell is fully occupied by the granular phase, $\alpha =0$ corresponds to a cell fully occupied by the gas phase, and $0\lt \alpha \lt 1$ denotes interfacial cells. The third term is an artificial compression term, with $\boldsymbol{u}_r$ denoting the relative interphase velocity. This contribution acts primarily in the interface region (typically around $\alpha$ approaching $\alpha \approx 0.5$ ) to counteract numerical diffusion, thus sharpening the interface and improving stability (Jasak Reference Jasak1996).

The momentum conservation equation is written as

(3.3) \begin{equation} \frac {\partial (\rho \boldsymbol{u})}{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }\left (\rho \boldsymbol{uu}\right )=\boldsymbol{\nabla }\boldsymbol{\cdot }\left (\eta\boldsymbol{\nabla }\boldsymbol{u}\right )+\left (\boldsymbol{\nabla }\boldsymbol{u}\right )\boldsymbol{\cdot }\left (\boldsymbol{\nabla }\eta \right )-\boldsymbol{\nabla }\!p_{\textit{eff}}- \textit{Gh}\,\boldsymbol{\nabla }\!\rho . \end{equation}

To improve numerical conditioning and stabilise the pressure solution, we employ the dynamic pressure $p_{\textit{eff}}=p-\rho \textit{Gh}$ , which separates the dynamic and hydrostatic contributions of the pressure gradient. The mixture density and viscosity are obtained via volume-fraction weighting, $\rho =\alpha \rho _{{g}}+(1-\alpha )\rho _{\textit{air}}$ and $\eta =\alpha \eta _{{g}}+(1-\alpha )\eta _{\textit{air}}$ , where the subscripts ‘g’ and ‘air’ denote the granular and gas phases, respectively. Because the viscosity field $\eta$ varies spatially across the two-phase interface, the viscous term in (3.3) comprises both the standard diffusion contribution $\boldsymbol{\nabla }\boldsymbol{\cdot }(\eta\boldsymbol{\nabla }\boldsymbol{u})$ and the additional term $(\boldsymbol{\nabla }\boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{\nabla }\eta$ , which accounts for the effects of viscosity gradients across the two-phase interface.

A distinctive feature of the present framework is the determination of the effective granular viscosity $\eta _g$ through the $\mu (I,\varTheta )$ rheology. It is expressed as

(3.4) \begin{equation} \eta _g=\frac {\mu (I,\varTheta )\,P}{\left |\dot {\gamma }\right |}, \end{equation}

where $|\dot {\gamma }|$ is the shear-rate magnitude, and $\mu (I,\varTheta )$ represents the constitutive relation relating the stress ratio to both the inertial number $I$ and granular temperature $\varTheta$ , as defined in (2.4). The granular temperature $\varTheta$ is explicitly solved through a transport equation proposed by Irmer et al. (Reference Irmer, Brodsky and Clark2025):

(3.5) \begin{equation} t_0\dot {\varTheta }=Dd^2\,{\nabla}^2\varTheta -B\varTheta +AI^a, \end{equation}

where $t_0$ is the characteristic elastic collision time, $D$ represents the diffusivity of granular temperature, $B$ quantifies collisional dissipation, and $A$ prescribes the local production of velocity fluctuations. These coefficients appearing in (3.5) retain an explicit parametric dependence on the coefficient of restitution through order-unity prefactors, although their overall sensitivity to restitution is expected to remain weak in the dense granular flows considered here. The exponent $a$ retains its kinetic theory interpretation. The diffusion term allows fluctuation energy to propagate over distances of several particle diameters, thereby incorporating non-local stress transmission, while the source term models production from shear or vibrating boundaries. Equation (3.5) is intended for dense, contact-dominated granular flows, where velocity fluctuations are governed by frictional multi-body interactions. In this regime, it provides closure to the momentum balance via the $\mu (I,\varTheta )$ rheology, whereas dilute, collision-dominated flows are more appropriately described by kinetic-theory-based models. Discrete-element simulations of oscillatory shear flows further demonstrate that boundary-driven granular temperature scales with the square of the vibration velocity amplitude (Irmer et al. Reference Irmer, Brodsky and Clark2025), yielding a dimensionless wall temperature proportional to $A^2(2\pi f)^2$ . When only local effects are considered, the unsteady, convective and diffusive contributions in (3.5) vanish, and the granular temperature reduces to the local expression

(3.6) \begin{equation} \varTheta _{\textit{loc}}=\frac {A}{B}I^{3/2}. \end{equation}

Under this simplification, the two-phase VOF framework reduces to a single-phase formulation applicable to dense granular flows. In this case, the interface-tracking equation is omitted, and in the momentum balance, the buoyancy contribution $-\textit{Gh} \boldsymbol{\nabla }\!\rho$ vanishes due to constant density. The viscosity gradient term $(\boldsymbol{\nabla }\boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{\nabla }\eta$ , however, must still be retained because $\eta _g$ remains spatially variable through its dependence on $\mu (I,\varTheta )$ . The continuity equation remains unchanged.

3.2. Discretisation and iterative schemes

The numerical implementation is carried out within the OpenFOAM $\circledR$ platform using a finite-volume formulation. All governing equations are discretised in a cell-centred framework, with surface and volume integrals evaluated through Gauss’ theorem (Rusche Reference Rusche2002; Versteeg Reference Versteeg2007). Temporal derivatives are advanced using a first-order implicit Euler scheme, which provides unconditional stability while preserving computational efficiency. Convective terms are discretised with a second-order linear upwind scheme combined with gradient reconstruction (Beam & Warming Reference Beam and Warming1976), whereas diffusive terms employ central differencing. For non-orthogonal meshes, explicit correction terms are applied to maintain second-order accuracy (Ferziger, Perić & Street Reference Ferziger, Perić and Street2019). Pressure-gradient and surface-tension contributions are interpolated at cell faces to enforce strict momentum conservation (Rhie & Chow Reference Rhie and Chow1983).

The volume-fraction transport equation is solved using the multidimensional universal limiter for explicit solution (MULES) algorithm (Weller et al. Reference Weller, Tabor, Jasak and Fureby1998), which applies flux-corrected transport (Zalesak Reference Zalesak1979; Deshpande, Anumolu & Trujillo Reference Deshpande, Anumolu and Trujillo2012) to ensure boundedness of the volume fraction while minimising interface smearing. Coupling between pressure and velocity fields is achieved through the PIMPLE algorithm (Issa Reference Issa1986; Holzmann Reference Holzmann2016), which integrates SIMPLE-type outer iterations with PISO-type inner corrections. Within each time step, the algorithm proceeds sequentially through momentum prediction, pressure correction via a Poisson equation derived from continuity, and velocity correction to restore incompressibility. Multiple PISO iterations are adopted to accelerate convergence.

Viscosity and friction coefficients are updated explicitly after each pressure–velocity loop. To maintain numerical stability, the effective viscosity is bounded according to $\eta _{\textit{min}}\le \eta \le \eta _{\textit{max}}$ , where $\eta _{\textit{min}}$ prevents divergence, and $\eta _{\textit{max}}$ accommodates the elevated viscosities characteristic of quasi-static regimes. Previous $\mu (I)$ implementations often set $\eta _{\textit{max}}$ consistent with the Drucker–Prager yield condition. In the present work, where non-local quasi-static effects are emphasised, excessively small upper limits truncate the relevant range of viscosities, whereas overly large values promote spurious oscillations. Numerical tests indicate that adopting $\eta _{\textit{max}}$ two orders of magnitude larger than conventional $\mu (I)$ choices yields stable and physically consistent results (Lagrée et al. Reference Lagrée, Staron and Popinet2011). Similar constraints are applied to the effective friction coefficient, with the lower bound prescribed as $\mu _s=0.1$ , and the upper bound limited to suppress numerical instabilities (Kim & Kamrin Reference Kim and Kamrin2020).

The granular-temperature transport equation is discretised using the FVM. Time integration is performed with a first-order implicit Euler scheme, while convective fluxes are approximated using a first-order upwind discretisation. Within the present VOF formulation, the granular-temperature equation is formally advanced everywhere, although granular temperature is physically meaningful only within the granular phase. To localise the dynamics of the granular temperature $\varTheta$ to prevent spurious numerical spreading into the pure-fluid region, a phase-shielding strategy based on the granular volume fraction $\alpha$ is introduced. The resulting fully discrete implicit form of the transport equation reads

(3.7) \begin{equation} \begin{aligned} & \underbrace { \frac {\alpha _P (\varTheta _P^{n+1} - \varTheta _P^n)}{\Delta t} V_P }_{\text{time derivative term}} + \underbrace { \sum _f \alpha _f F_f \varTheta _f^{{upwind}} }_{\text{convection term}} - \underbrace { \sum _f \alpha _f D_{\!f} (\boldsymbol{\nabla }\varTheta )_f \boldsymbol{\cdot }\boldsymbol{S}_f }_{\text{diffusion term}} \\ & =\underbrace { \alpha _P V_P (-B t_0^{-1} \varTheta _P^{n+1} + A t_0^{-1} I_P^{3/2}) }_{\text{source term}} - \underbrace { (1 - \alpha _P) \chi _{\textit{air}} \varTheta _P^{n+1} V_P }_{\text{air-phase damping}}, \end{aligned} \end{equation}

where superscripts $n$ and $n+1$ denote successive time levels. The subscript $P$ refers to the control volume, and $f$ denotes the cell face. The quantities $\alpha _P$ and $\alpha _f$ represent the granular volume fraction at the cell centre and the cell face, respectively. The volumetric flux across face $f$ is defined as $F_f=\boldsymbol{u}_f\boldsymbol{\cdot }\boldsymbol{S}_f$ , where $\boldsymbol{S}_f$ is the face area vector. The upwind value $\varTheta _f^{{upwind}}$ is taken as $\varTheta _P^{n+1}$ if $F_f\geq 0$ , and $\varTheta _N^{n+1}$ otherwise, with $N$ denoting the neighbouring cell. The diffusion coefficient is defined as $D_{\!f}=Dd_g^2t_0^{-1}$ . For numerical stability, the dissipation term $-Bt_0^{-1}\varTheta _P^{n+1}$ is treated implicitly, whereas the production term $At_0^{-1}I_P^{3/2}$ is evaluated explicitly.

From a physical perspective, the volume-fraction weighting confines the transport, production and dissipation of granular temperature to regions occupied by grains. The additional damping term proportional to $(1-\alpha _P)$ enforces a rapid decay of $\varTheta$ at the phase interface, and ensures that granular temperature remains negligible in the pure-fluid phase. From a numerical perspective, this phase-shielding strategy suppresses spurious diffusion of $\varTheta$ across the interface, and prevents unphysical activation of the $\mu (I,\varTheta )$ rheology in fluid-only regions, thereby enhancing the robustness of the coupled VOF–granular rheology framework.

Model parameters are prescribed as follows. The elastic collision time scale is given by $t_0=d\sqrt {\rho /E}$ , with representative particle parameter $d=0.5\times 10^{-3}$ m, $\rho =2450$ $\mathrm{kg\,m^{-3}}$ , and $E=1.3\times 10^7$ $\mathrm{Pa}$ . Following Irmer et al. (Reference Irmer, Brodsky and Clark2025), the granular-temperature decay length $\sqrt {D/B}$ is of order unity under strong vibration ( $\varGamma \gt 10$ ), and varies between 1 and 4 under weaker forcing. In the present work, we adopt $D=1$ and $B=1$ , consistent with kinetic theory scaling. Using the scaling relation $\varTheta _{\textit{loc}}\propto I^{3/2}$ , we choose $A/B=0.15$ , yielding $A=0.15$ . The granular material considered in the computational examples corresponds to glass beads, for which the rheological parameters are specified as $\mu _s=0.38$ , $\mu _2=0.64$ and $I_0=0.279$ .

The overall solution procedure follows a segregated strategy. Within each time step, the solver advances sequentially through the volume-fraction equation, momentum prediction, pressure Poisson equation and velocity correction. Outer iterations are performed to ensure global convergence. A schematic representation of the iterative process is shown in figure 1.

Figure 1.

Algorithm flowchart. The two-phase flow solution in this work is based on interFoam, while the single-phase flow solution is based on pimpleFoam. Both employ an alternating inner–outer iteration scheme, with the distinction that two-phase flow solution requires prior solution of the phase-fraction equation.

3.3. Implementation of boundary conditions

Accurate representation of boundary conditions is essential for capturing velocity and granular-temperature distributions in creeping regions adjacent to confining walls. Since the FVM is based on flux variables, direct implementation of Coulomb-type frictional boundary conditions is not straightforward. Techniques such as Lagrange multipliers or penalty methods, commonly used in finite-element approaches (Simo & Laursen Reference Simo and Laursen1992; Huněk Reference Huněk1993; Ferrás et al. Reference Ferrás, Nóbrega and Pinho2013), are not directly applicable in the present framework. To address this challenge, a Navier-type partial-slip condition is employed. The tangential velocity at the boundary is expressed as

(3.8) \begin{equation} \boldsymbol{u}_b=f\boldsymbol{\cdot }\boldsymbol{u}_{\textit{ref}}+\left (1-f\right )\boldsymbol{\cdot }\left (\boldsymbol{I}-\hat {\boldsymbol{n}}\otimes \hat {\boldsymbol{n}}\right )\boldsymbol{\cdot }\boldsymbol{u}_o, \end{equation}

where $f$ denotes the slip coefficient controlling the degree of slip, $\boldsymbol{u}_{\textit{ref}}$ is the reference velocity (zero in the present context), $\boldsymbol{I}$ is the identity tensor, $\hat {\boldsymbol{n}}$ is the outward unit normal at the boundary, and $\boldsymbol{u}_o$ represents the velocity of the adjacent control volume. This formulation interpolates between the no-slip limit ( $f=1$ ) and perfect slip ( $f=0$ ), thereby providing a flexible representation of wall–particle interactions. In the present simulations, the coefficient $f$ is calibrated against experimental flow-rate scaling to ensure consistency with observed bulk behaviour.

For the granular-temperature field, zero-gradient Neumann boundary conditions are imposed to avoid physically inconsistent constraints. As granular temperature is generated by particle vibration and quantifies velocity fluctuations, it is particularly sensitive to boundary specification. When vibration acts as an external energy source, Dirichlet-type conditions can be prescribed in terms of excitation amplitude and frequency (Irmer et al. Reference Irmer, Brodsky and Clark2025). However, the assumption of vanishing granular temperature at a rough, fixed base is valid only in creeping regimes where the velocity profile decays exponentially. Numerous experimental and numerical studies have shown that in dense inertial flows with Bagnold-type profiles, the granular temperature remains finite near the wall; indeed, for frozen-particle base boundaries, its magnitude can remain comparable to bulk values (Bharathraj & Kumaran Reference Bharathraj and Kumaran2019). These observations indicate that enforcing $\varTheta =0$ at the boundary would amount to imposing an artificial dissipation sink and is generally inappropriate for dense flows. Instead, the zero-gradient Neumann condition provides a consistent and robust approximation, ensuring compatibility with both creeping and inertial flow regimes.

Figure 2.

Schematic diagram of steady dense granular flow calculation.

4.1. Flow configuration and model validation

Steady heap flows constitute a canonical benchmark for assessing granular rheological models, as they exhibit two distinct regimes within a single configuration. At the free surface, a rapid flowing layer develops whose thickness $h$ scales with the imposed volumetric flux $Q$ and is well described by local inertial rheology. Beneath this layer, however, the material undergoes slow creeping motion characterised by exponentially decaying velocities with depth. The decay length $\lambda$ , which governs the quasi-static region, depends only on system geometry and remains independent of $Q$ , thereby reflecting the intrinsic role of non-local effects.

Heap-flow experiments are often favoured because they generate stable, quasi-steady configurations suitable for statistical analysis and the derivation of scaling relations between flow rate and surface properties (Jop et al. Reference Jop, Forterre and Pouliquen2005). In numerical studies, however, channel-flow geometries are frequently adopted, primarily for computational efficiency: plane shear assumptions reduce three-dimensional problems to two dimensions (Liu & Henann Reference Liu and Henann2017), and boundary conditions are simplified. Despite these differences, both heap and channel flows display the same essential phenomenology, namely the coexistence of a rapid surface layer and a creeping bulk.

In the present study, we consider the classical configuration of gravity-driven channel flow (figure 2). A granular layer of depth $H$ is confined between two glass plates separated by width $W$ and inclined at angle $\theta$ . Boundary conditions are specified as follows. The basal surface is modelled as rough and no-slip. The sidewalls are treated as frictional boundaries using partial-slip boundary conditions. At the top, a free-slip boundary condition is imposed, corresponding to zero shear stress and zero normal velocity at the free surface. For the pressure field, a zero gauge pressure is prescribed at the free surface, while zero-gradient (Neumann) boundary conditions are applied at the base and the sidewalls.

Steady-state simulations are performed for inclined channel flows. Convergence is deemed to be achieved when two criteria are satisfied simultaneously. First, the normalised residuals of velocity and pressure obtained from the solid-phase governing equations, together with the normalised residual of the granular temperature obtained from the granular-temperature transport equation, all decrease below $10^{-4}$ . Second, the relative variation of the volumetric flux between successive iterations is less than $1\,\%$ . To minimise basal boundary effects, the channel depth is fixed at $H/d=20$ . All velocities are normalised by $\sqrt {Gd}$ . The resulting velocity fields reproduce the characteristic features of heap flows, including a reduction of velocity adjacent to the sidewalls, and enhanced transport in the channel core.

Figure 3.

Flow field velocity contours and wall velocity distribution calculated in narrow channels with $W/d=19$ . Normalised velocity contours obtained using the $\mu (I,\varTheta )$ model at (a) $\theta =26.1^\circ$ and (b) $\theta =28.0^\circ$ . Distribution of normalised wall velocity versus normalised depth in (c) linear coordinates and (d) semi-logarithmic coordinates. Solid lines represent the $\mu (I,\varTheta )$ model, dashed lines represent the $\mu (I)$ model, dash-dotted lines represent the NGF model with results extracted from Liu & Henann (Reference Liu and Henann2017), and symbols represent experimental data from Jop et al. (Reference Jop, Forterre and Pouliquen2007).

Figure 3 summarises the simulation results for narrow channels with $W/d=19$ . Figures 3(a,b) show normalised velocity contours at inclination angles $\theta =26.1^\circ$ and $\theta =28.0^\circ$ , respectively. In both cases, the $\mu (I,\varTheta )$ model captures the key flow features observed in experiments, namely reduced wall velocities and accelerated transport along the channel centreline. Figures 3(c,d) present depth-dependent wall velocity profiles in linear and semi-logarithmic coordinates. Both the $\mu (I,\varTheta )$ and NGF models successfully reproduce the exponential velocity decay characteristic of creeping regions, in close agreement with experimental measurements of Liu & Henann (Reference Liu and Henann2017). At higher inclination angles, a slight underestimation of velocity in the deep creeping region is observed. This deviation is likely related to the calibration of the wall boundary parameters against the bulk flow rate, which primarily reflects inertial-region dynamics and is therefore less sensitive to creep. In contrast, the classical $\mu (I)$ rheology predicts an abrupt cessation of motion once the yield criterion is violated, thereby failing to capture the creeping velocity tails.

Figure 4.

Velocity fields in a wide channel configuration ( $W/d=142$ ). Normalised velocity contours $v/\sqrt {Gd}$ displayed in the $y$ – $z$ cross-section, obtained from the $\mu (I,\varTheta )$ non-local rheological model for two volumetric flux conditions: (a) $Q/Wd\sqrt {Gd}=4.4$ and (b) $Q/Wd\sqrt {Gd}=91$ . (c) Transverse distribution of the dimensionless surface velocity $v_{\textit{surf}}/\sqrt {Gd}$ across the channel width $z$ . Experimental measurements from Jop (Reference Jop2006) are indicated by symbols, and predictions from the local $\mu (I)$ rheology are shown as dashed lines, while solid lines represent results from the non-local $\mu (I,\varTheta )$ formulation. The dash-dotted lines denote the NGF model (Liu & Henann Reference Liu and Henann2017). (d) Depth profiles of the dimensionless wall velocity $v_{wall}/\sqrt {Gd}$ as a function of distance $y$ beneath the free surface, shown for various flow rates $Q/Wd\sqrt {Gd}$ .

To examine the influence of lateral confinement, additional simulations were conducted in wide channels with $W/d=142$ , as shown in figure 4, and compared against the inclined-plane hopper–heap experiments of Jop and others (Jop et al. Reference Jop, Forterre and Pouliquen2005, Reference Jop, Forterre and Pouliquen2006, Reference Jop, Forterre and Pouliquen2007; Jop Reference Jop2006). The flow rate $Q$ considered here is defined as the total volumetric flux, consistent with the definition in Liu & Henann (Reference Liu and Henann2017), but differing from the unit-width flux adopted in Jop et al. (Reference Jop, Forterre and Pouliquen2005, Reference Jop, Forterre and Pouliquen2007). The dimensionless discharge is expressed as $Q/Wd\sqrt {Gd}$ . The numerical set-up does not precisely replicate experimental flow-control conditions; a proportional–integral–derivative (PID) feedback scheme was employed to regulate the imposed flux (see Appendix A).

Normalised velocity contours for $Q/Wd\sqrt {Gd}=4.4$ and $91$ (figures 4 a,b) reproduce the expected qualitative features, with accelerated motion in the central region, and reduced velocities near the sidewalls. The comparison of surface velocity profiles $v_{\textit{surf}}(x,y=0)$ with experimental measurements (Jop Reference Jop2006) across four volume flow rates (figure 4 c) demonstrates that the similar predictive capability of both the $\mu (I)$ and $\mu (I,\varTheta )$ models indicates that the surface velocity distribution is largely independent of non-local effects. By contrast, wall velocity profiles in figure 4(d) provide a more discriminating test. The $\mu (I,\varTheta )$ model exhibits significantly improved agreement with experimental data, whereas the classical $\mu (I)$ rheology systematically underestimates velocities in deeper wall-adjacent regions. Compared with the NGF model (Liu & Henann Reference Liu and Henann2017), the $\mu (I,\varTheta )$ rheology shows slightly reduced accuracy in the regime where normalised velocities fall below $10^{-3}$ . This discrepancy is attributed primarily to artefacts of the partial-slip boundary conditions employed here, rather than to inherent deficiencies of the $\mu (I,\varTheta )$ rheology.

Figure 5.

Relationships in steady, dense granular heap flow: (a) free-surface inclination angle ( $\theta$ ) versus dimensionless volume flow rate $Q^\prime$ , and (b) dimensionless maximum velocity $v_{\textit{max}}^\prime$ versus dimensionless flow rate $Q^\prime$ . Symbols represent experimental data from Jop et al. (Reference Jop, Forterre and Pouliquen2005). Dash-dotted lines represent theoretical predictions from the NGF model reported in Liu & Henann (Reference Liu and Henann2017). Dashed lines represent theoretical predictions from the $\mu (I)$ model extracted from Jop et al. (Reference Jop, Forterre and Pouliquen2005). Solid lines represent theoretical predictions from the $\mu (I,\varTheta )$ model obtained in the present work.

4.2. Macroscopic scaling relationships

To further evaluate the effectiveness of non-local rheological models in regimes where non-local effects play only a secondary role, we examine the scaling behaviour of macroscopic flow quantities. Unlike localised measurements in creeping regions, bulk observables – such as the total volumetric flux $Q$ , the centreline surface velocity $v_{\textit{max}}=v(y=W/2,z=0)$ , and the steady inclination angle $\theta$ – average the influence of non-local cooperativity across the entire flow domain. Consequently, their consistency with established scaling relations provides a stringent test of whether a non-local model correctly degenerates to local rheology when appropriate. Jop et al. (Reference Jop, Forterre and Pouliquen2005) demonstrated that for heap flows with different stable discharge rates, the relationships between the bulk flow rate, surface velocity and slope inclination collapse onto universal curves when expressed in the dimensionless form

(4.1a,b) \begin{align} Q^\prime =Qd/G^{1/2}W^{7/2}, \quad v^\prime _{\textit{max}}=v_{\textit{max}}d/GW^{3/2}, \end{align}

where $Q'$ and $v'_{\textit{max}}$ represent normalised discharge and surface velocity, respectively. Figure 5 compares simulation results from the $\mu (I,\varTheta )$ framework with these scaling laws. The partial-slip wall coefficient was selected consistently with the finite-element implementation of the NGF model, and multiple flow datasets were employed to enhance robustness (figure 5 a). Both the $\mu (I,\varTheta )$ and NGF models closely reproduce the experimental scaling relations, showing only minor deviations from measurements. Interestingly, the classical $\mu (I)$ rheology yields an almost perfect collapse of the normalised quantities. Nevertheless, the differences among all three models remain small, confirming that creeping zones do not influence global scaling. Thus the $\mu (I,\varTheta )$ and NGF models successfully capture creeping behaviour while preserving the predictive accuracy of macroscopic observables governed primarily by local inertial rheology.

Notably, the collapse observed in figure 5 highlights a fundamental property of granular rheology: within the rapidly flowing surface layer, where the inertial number is sufficiently large, non-local contributions become negligible, and global observables continue to obey simple scaling laws. This finding supports the argument of Kim & Kamrin (Reference Kim and Kamrin2020), who proposed that constitutive forms of the type $\mu \varTheta ^p=f(I)$ provide a unified description of granular flows across diverse geometrical configurations while preserving the validity of local rheology in its applicable regime.

A simple scaling analysis of the $\mu (I,\varTheta )$ model provides further insight. Introducing a characteristic length $L$ to represent the typical scale over which the granular temperature $\varTheta$ varies significantly, the diffusion term can be approximated as ${\nabla} ^2\varTheta \sim \varTheta /L^2$ . In the absence of appreciable non-locality, $\varTheta \approx \varTheta _{\textit{loc}}$ . Under steady-state conditions, $Dd^2\varTheta /L^2\ll B\varTheta$ ; with $B/D=1$ under strong vibration, this criterion reduces to $L\gg d$ . Thus when the characteristic variation length of $\varTheta$ greatly exceeds the particle diameter, diffusion effects become negligible. The flowing surface layer satisfies precisely this condition, explaining why non-local contributions vanish and why global rheological scaling relations remain valid.

4.3. Distinct length scales in layered flows

To elucidate the qualitative differences between the rapidly flowing surface layer and the deep creeping region, two classical length scales are employed (Jop et al. Reference Jop, Forterre and Pouliquen2007). The surface flowing layer is characterised by the mean flow depth $h=Q/(Wv_{\textit{max}})$ , where $v_{\textit{max}}$ denotes the centreline surface velocity. The creeping bulk is quantified by the decay length $\lambda$ , defined through the exponential form of the longitudinal velocity profile, $v(z=z_0,y)\sim \exp {(-y/\lambda )}$ , with $z_0$ denoting the channel-centre position. In experiments, owing to limited velocity field resolution, $\lambda$ is typically extracted from wall-based velocity profiles (Jop et al. Reference Jop, Forterre and Pouliquen2007). In the present simulation with $H/d=20$ , the creeping region extends up to approximately ten particle diameters above the base. Figure 6(a) presents the variation of normalised length scales $h/d$ and $\lambda /d$ with the dimensionless discharge $Q^\prime$ for both the $\mu (I,\varTheta )$ and NGF models in a narrow channel with $W/d=19$ . The results reproduce a central observation from Jop et al. (Reference Jop, Forterre and Pouliquen2007): the mean flow depth $h$ increases with discharge, whereas the decay length $\lambda$ remains essentially independent of flow rate. Both models achieve good quantitative agreement with experimental data, although the $\mu (I,\varTheta )$ model exhibits a slight residual dependence of $\lambda$ on velocity. Figure 6(b) further compares the scaling relation between $\lambda /d$ and channel width $W/d$ . The $\mu (I,\varTheta )$ predictions are consistent with NGF results, and both align with the $1/4$ power-law conjecture proposed by Liu & Henann (Reference Liu and Henann2017), shown as the dashed line.

Figure 6.

Comparison of the $\mu (I,\varTheta )$ model and NGF model at $W/d=19$ . (a) Relationship between flowing zone size $h/d$ and creeping zone decay length $\lambda /d$ versus dimensionless volume flow rate $Q^\prime$ (see text for definition). Solid lines represent predictions of the $\mu (I,\varTheta )$ model obtained in the present work, whereas dash-dotted lines represent NGF model results extracted from Liu & Henann (Reference Liu and Henann2017). (b) Relationship between creeping zone decay length $\lambda /d$ and channel width $W/d$ . The dashed line represents the reference $1/4$ power law.

These trends reflect fundamental differences in rheological formulation. In the NGF framework (Kamrin & Koval Reference Kamrin and Koval2012), the flowing surface layer follows the local constitutive law when $\mu \gt \mu _s$ , producing a rate-dependent flow thickness $h$ . In the creeping bulk, where $\mu \lt \mu _s$ , the model prescribes a vanishing local fluidity, $g_{loc}=0$ , thereby completely decoupling the local rheology from the driving rate, and producing a strictly rate-independent decay length $\lambda$ . In contrast, within the $\mu (I,\varTheta )$ framework, the surface layer also behaves locally when the inertial number is large, as $\varTheta _{\textit{loc}}/\varTheta \approx 1$ . In the creeping regime, however, as $I\rightarrow 0$ , the function $f(I)$ tends to zero and $\varTheta _{\textit{loc}}/\varTheta \ll 1$ , allowing a small amount of fluctuation energy to diffuse from the inertial layer into the creep region. This residual coupling produces a weak but finite rate dependence of $\lambda$ , rather than the strict rate independence enforced in NGF. Such behaviour is physically plausible and consistent with independent theoretical analyses (Martinez et al. Reference Martinez, González-Lezcano, Batista-Leyva and Altshuler2016). Thus the divergence between the two models ultimately reflects whether granular rheology becomes completely decoupled from inertial effects in the quasi-static limit.

Overall, the simulations demonstrate that the $\mu (I,\varTheta )$ model consistently captures the layered structure of dense heap flows, characterised by two coexisting length scales: the mean flow depth $h$ , reflecting the flux-dependent inertial dynamics of the surface layer, and the decay length $\lambda$ , encoding the geometry-controlled creep in the bulk. Across both narrow and wide channels, the $\mu (I,\varTheta )$ formulation reproduces these behaviours with quantitative accuracy. In narrow channels, it recovers the exponential velocity decay observed experimentally and predicted by NGF, overcoming the limitations of the classical $\mu (I)$ rheology. In wide channels, it yields wall velocity distributions in close agreement with experiments, and performs comparably to NGF, while retaining the computational efficiency of a local constitutive law. Collectively, these findings confirm that incorporating granular temperature substantially extends the applicability of inertial rheology, enabling a unified description of dense granular flows that encompasses both local inertial dynamics and non-local creeping phenomena.

5.1. Simulation configuration and parameters

Shear localisation is a defining feature of granular flows, in which deformation is concentrated within narrow shear bands while the surrounding bulk remains quasi-rigid. Predicting the width and morphology of these bands remains a central challenge in granular rheology. Local constitutive models, including the classical inertial rheology, typically yield shear zones spanning only a few grain diameters with unrealistically sharp gradients, in marked contrast to experimental observations (Jop Reference Jop2008). Measurement is further complicated by the frequent anchoring of shear bands to sidewalls. Split-bottom shear cells provide a particularly useful geometry in this regard: the shear band originates from the basal discontinuity and broadens upwards into the bulk, offering a clear manifestation of intrinsic localisation. Experimental studies using this configuration have consistently highlighted the importance of non-local effects (Fenistein & van Hecke Reference Fenistein and van Hecke2003; Fenistein et al. Reference Fenistein, van de Meent and van Hecke2004).

In the present work, shallow split-bottom Couette flows ( $H\le 50\,\mathrm{mm}$ ) are simulated using a two-phase VOF framework coupled with the $\mu (I,\varTheta )$ rheology. The apparatus configuration is illustrated in figure 7 and consists of an inner stationary cylinder of radius $R_i=65$ mm, a fixed central disk of radius $R_s=85\,\mathrm{mm}$ , and an outer rotating cylinder of radius $R_o=105\,\mathrm{mm}$ . Both the outer cylinder and basal disk rotate at constant angular velocity $\varOmega =0.16\,\mathrm{rad\,s^{-1}}$ , generating shear along the basal discontinuity between fixed and moving boundaries. To avoid spurious pressure constraints arising from centrifugal effects in annular geometries, the granular free surface is evolved dynamically through the two-phase VOF representation. No-slip velocity conditions are imposed on all solid boundaries, while pressure and granular temperature fields satisfy zero-gradient Neumann conditions. The initial granular column height is prescribed through the volume-fraction field.

Classical local models invariably predict shear bands of vanishing width. By contrast, the NGF model successfully reproduces the experimentally observed scaling of shear-band width with system size $W(d,H)$ , as well as the dome-shaped transition in deeper layers, by introducing a cooperativity length $\xi \propto d$ that reflects the scale of collective particle rearrangements (Henann & Kamrin Reference Henann and Kamrin2013). More recently, Faroux et al. (Reference Faroux, Washino, Tsuji and Tanaka2022) demonstrated that finite-volume implementations of non-local formulations can quantitatively capture these behaviours. Building on these advances, we combine the two-phase VOF method with the $\mu (I,\varTheta )$ rheology to assess its predictive capability for shear-band width and morphology, benchmarking the results against both theoretical expectations and experimental measurements.

Figure 7.

Schematic diagram of the split-bottom Couette cell.

5.2. Non-local prediction of shear-band characteristics

The normalised angular velocity field, defined as $\omega =v_\theta /(r\varOmega )$ , provides direct insight into shear localisation in split-bottom Couette flows. As shown in figure 8, shallow granular columns ( $H=10\,\mathrm{mm}$ ) exhibit significantly narrower shear bands than taller columns ( $H=30\,\mathrm{mm}$ ), while systems composed of smaller particles ( $d=0.35\,\mathrm{mm}$ ) develop narrower bands than those with larger particles ( $d=1.2\,\mathrm{mm}$ ). These observations confirm that shear-band width is jointly governed by the confinement geometry and the particle size. The $\mu (I,\varTheta )$ formulation predicts shear bands that broaden progressively with depth and display smooth velocity gradients, consistent with experimental reports of band widths exceeding 50 grain diameters (Fenistein & van Hecke Reference Fenistein and van Hecke2003). In contrast, the classical $\mu (I)$ rheology produces sharply localised deformation zones only a few diameters wide, reflecting the inherent limitations of local laws in the quasi-static limit.

Figure 8.

Contour plots of the normalised angular velocity field in the $r$ – $z$ plane of a split-bottom cell, obtained from simulations using different models at column heights 10 and 30 mm: (a–d) predictions from the $\mu (I,\varTheta )$ model; (e)–(h) predictions from the $\mu (I)$ model. Panels (a,b,e, f) represent particle diameter $d = 0.35\,\mathrm{mm}$ , whereas (c,d,g,h) correspond to $d = 1.2\,\mathrm{mm}$ .

Particle-size dependence further highlights the underlying mechanism. While both rheologies predict wider bands for larger particles, the effect is significantly stronger in the $\mu (I,\varTheta )$ framework. This outcome highlights grain diameter as an intrinsic length scale for shear localisation. In the NGF model, this role is captured by the cooperativity length $\xi \propto d$ , representing the extent of correlated particle rearrangements (Henann & Kamrin Reference Henann and Kamrin2013). In the $\mu (I,\varTheta )$ formulation, it arises naturally from the diffusion of granular temperature, which transmits fluctuation energy over multiple grain diameters. Accordingly, the $\mu (I,\varTheta )$ predictions agree qualitatively with experimental evidence of wide shear zones, whereas the $\mu (I)$ rheology yields unrealistically sharp gradients akin to narrow-band variational predictions (Unger et al. Reference Unger, Török, Kertész and Wolf2004). Quantitative comparisons reinforce these findings, as shown in figure 9. For particles with $d=0.35\,\mathrm{mm}$ at heights $H=10$ , $20$ and $30\,\mathrm{mm}$ , $\mu (I,\varTheta )$ simulations reproduce the experimentally observed broadening of shear bands with depth and the progressive inward migration of the band centre (figures 9 a,b), in close agreement with Fenistein & van Hecke (Reference Fenistein and van Hecke2003). Importantly, shallow columns ( $H\le 30\,\mathrm{mm}$ ) exhibit vertical bands that do not intersect the inner wall, again consistent with experimental reports. When fitted to the universal error function form

(5.1) \begin{equation} \omega \left (\lambda \right )=1/2+1/2\,\text{erf}(\lambda ), \end{equation}

the $\mu (I,\varTheta )$ results for $d=1.2\,\mathrm{mm}$ confirm excellent agreement across multiple filling heights. In contrast, the $\mu (I)$ rheology fails to reproduce the universality, predicting band widths that remain nearly independent of particle size (figure 9 c).

Figure 9.

Radial profiles of the normalised angular velocity at the flow surface, obtained from simulations using different rheological models and compared with experimental and theoretical references. (a) Comparison of the $\mu (I, \varTheta )$ and $\mu (I)$ models with experimental data from Fenistein & van Hecke (Reference Fenistein and van Hecke2003) for particle diameter $d = 0.35\,\mathrm{mm}$ . (b) Comparison of the $\mu (I, \varTheta )$ and $\mu (I)$ models with the velocity-scaling predictions of Fenistein & van Hecke (Reference Fenistein and van Hecke2003) for $d = 1.2\,\mathrm{mm}$ . (c) Compilation of all $\mu (I)$ model simulation results for varying column heights and particle diameters.

Additional simulations were carried out for particle diameters $d = 0.35{-}2.2\,\mathrm{mm}$ with filling heights $H = 5{-}35\,\mathrm{mm}$ , and the results are summarised in figure 10. The predicted dependence of shear-band width on layer depth across different particle sizes exhibits excellent agreement with both experimental measurements and the NGF model predictions (Henann & Kamrin Reference Henann and Kamrin2013). Specifically, the $\mu (I,\varTheta )$ predictions collapse onto the power-law relation

(5.2) \begin{equation} {W}/{d}\propto \left ({H}/{d}\right )^{2/3} \end{equation}

as given in the inset of figure 10, consistent with the scaling proposed by Henann & Kamrin (Reference Henann and Kamrin2013). This scaling highlights that the band width is not governed by stress gradients alone, but instead reflects cooperative particle rearrangements over a length scale proportional to grain size. By contrast, the $\mu (I)$ model systematically underestimates band width and fails to capture the experimentally observed broadening. It predicts bands of nearly constant and unrealistically small thickness, essentially independent of particle size, owing to the vanishing inertial number in the quasi-static regime.

Figure 10.

Relationship between the shear-band width $W$ and the packing height $H$ of the granular column for different particle diameters ( $d = 0.35$ , $0.8$ , $1.2$ and $2.2\,\mathrm{mm}$ , from bottom to top). Upon normalising both $W$ and $H$ by the particle diameter $d$ , the simulation data collapse onto a $2/3$ power-law scaling, $W/d \propto (H/d)^{2/3}$ , as illustrated in the inset, consistent with the observations of Henann & Kamrin (Reference Henann and Kamrin2013).

The split-bottom Couette simulations demonstrate that the $\mu (I,\varTheta )$ rheology provides a robust and physically consistent description of shear localisation in quasi-static granular flows. Unlike the classical $\mu (I)$ model, which predicts unrealistically narrow and size-independent shear zones, the $\mu (I,\varTheta )$ framework captures the progressive broadening of shear bands with depth, their inward migration with increasing filling height, and the particle-size dependence of shear-band width. The resulting velocity profiles collapse onto universal error function forms, and the width scaling $W\propto d^{2/3}H^{1/3}$ is reproduced in close agreement with both experimental measurements and NGF predictions. These results confirm that granular-temperature-based cooperativity furnishes a mechanistic explanation for non-local stress redistribution, reconciling the finite shear-band widths observed experimentally with continuum modelling. Collectively, the findings establish the $\mu (I,\varTheta )$ formulation as a unified rheological framework capable of recovering inertial behaviour in rapid flows while extending predictive accuracy to slow, quasi-static regimes dominated by non-local effects.