Coupling non-local rheology and volume of ﬂuid (VOF) method: a ﬁnite volume method (FVM) implementation

. Additional to a behavior switching between solid-like and liquid-like, dense granular ﬂows also present propagating grain size-dependent e ﬀ ects also called non-local e ﬀ ects. Such behaviors cannot be ef-ﬁciently modeled by standard rheologies such as µ ( I )-rheology but have to be dealt with advanced non-local models. Unfortunately, these models are still new and cannot be used easily nor be used for various conﬁg-urations. We propose in this work a FVM implementation of the recently popular NGF model coupled with the VOF method in order to both make non-local modeling more accessible to everyone and suitable not only for single-phase ﬂows but also for two-phase ﬂows. The proposed implementation has the advantage to be extremely straightforward and to only require a supplementary stabilization loop compared to the theoretical equations. We then applied our new framework to both single and two-phase ﬂows for validation.


Introduction
Granular materials and more specifically granular flows are extremely ambiguous by nature in the sense that they display characteristics from both solids and liquids (even gas if we consider fast dilute flows). This duality has posed numerous problems to scientists who had to model and predict their behaviors. This changing nature is even more challenging near the transition between the flowing state and the jamming state, i.e. slow granular flows, where the propagation of microscopic stress within the media can influence the macroscopic flow phenomenology. We usually refer to this class of microscopic effects as non-local effects. Non-local effects are grain size dependent effects which propagate over a cooperative length of a few times the grain diameter and can manifest themselves in a variety of ways: creep flow developing underneath surface flows, grain size dependency of the tilt and stop angles for shallow layers down an inclined plane as well as the shear band width in various geometries [1], mechanicallyinduced creep flow far from the sheared area [2] or the wave-number dependency in inhomogeneous kolmogorov flows [3].
Capturing simultaneous slow creep flows and rapid surface flows is still an open question. While the Discrete Element Method (DEM) [4] could take on this challenge by individually tracking the grains, its calculation cost makes it very difficult to use without powerful computers. On the other hand, continuum models are more alike to Computational Fluid Dynamics (CFD) in which they are based on governing and constitutive relations. However, even the now famous µ(I)-rheology [5,6] which combines a rate dependent constitutive equation with an appropriate yield criterion specific to plasticity theory fails to capture the grain size dependencies within the flow and the smooth transition between the different states.
In the recent years, many spatially dependent non-local models have been proposed to address these shortcomings [7,8]. The one we will focus our work on is called steadystate Non-local Granular Fluidity (NGF) model and has been developed by Kamrin and Koval [9]. This model is inspired from the Kinetic Elasto-plastic (KEP) theory applied to emulsion flows [10,11] and already benefits from a Finite Element Method (FEM) [12] realization but has yet to be implemented in different frameworks.
Non-local modeling is still in its first steps, hence both its range of application and its accessibility are limited. In order to answer these issues, we propose a Finite Volume Method (FVM) [13] implementation of the steady-state NGF model enhanced with the Volume-of-Fluid (VOF) [14] method. The FVM discretization is one of the most popular framework used in CFD and is already being used by various commercial and open-source software. As for the VOF, it will make the modeling of non-local two-phase flows, i.e granular material coexisting with a second fluid, possible and extend the field of application of non-local modeling.

Non-local granular fluidity model
While the classical µ(I) rheology assumes that the stress tensor at a given location is a function of the shear rate at the same place [5,6], the NGF model is weakly nonlocal in the sense that the constitutive equation involves an additional state variable whose evolution is controlled by an independent equation sensitive to spatial gradients [9]. This new state variable is denoted as granular fluidity g and satisfies the following flow rule: whereγ is the strain rate and µ is the stress ratio µ = τ/P with τ the shear stress and P the pressure. The granular fluidity is a positive scalar state variable which may be interpreted as a pressure-weighted inverse viscosity and had initially been chosen for practical reasons leading to concerns about its physical meaning [7]. However, Zang and Kamrin showed that the granular fluidity could be described as a kinematic state variable through the relation g = (δv/d)F(φ), with δv the microscopic velocity fluctuations, d the grain size and φ the solid fraction [15]. Practically, the granular fluidity g can be seen as "the susceptibility to flow", i.e., a low value corresponds to slow creep flows while a high value refers to rapid flows. Similar to what Jop et al. did to generalize the µ(I) rheology [6], the flow rule (1) can be extended to higher dimensions [16]. We define the strain-rate tensor as Γ i j = (∂v i /∂x j + ∂v j /∂x i )/2, with v i the velocity field components and x i the spatial coordinates. We now assume that steady fully developed flows are isochoric which results in Γ kk = 0. The equivalent strain rate can then be written aṡ γ = 2Γ i j Γ i j . On the other hand, the Cauchy stress tensor and deviatoric tensor are defined through σ i j = σ ji and σ i j = σ i j − (1/3)(σ kk )δ i j whilst the shear stress, pressure and stress ratio are defined using the stress tensor invariants as τ = σ i j σ i j /2, P = −σ kk /3 and µ = τ/P. Finally, assuming the colinearity between the strain rate tensor and deviatoric stress tensor, we obtain with I the identity tensor. Since our main topic of interest in this paper is fully developed steady flows, we only consider the steady state form of the NGF rheology [17] where g loc and ξ are the local fluidity and the non-local diffusion length scale called cooperativity length, respectively, and are defined as follow: where d is the grain diameter, ρ s is the grain material density, µ s , µ 2 and I 0 are constant dimensionless material parameters whose values can be borrowed from the local rheology as µ s = µ (I → 0), µ 2 = µ (I → ∞), I 0 a dimensionless parameter characterizing the nonlinear rate-dependent response [18], ∆µ = µ 2 − µ s and A is a constant dimensionless material parameter which is also called non-local amplitude [17]. The model reduces to the local µ(I)-rheology when the flow is homogeneous and the Laplacian term vanishes: The NGF model has two contributions, one local with g loc and one non-local which is directly proportionate to the grain size d. This Laplacian term implies that the stress tensor is no longer a function of only the shear rate at the same location but is instead dependent on the neighbouring values. A straightforward way to compare both the local and non-local models is to consider an inhomogeneous flow where µ < µ s . In this situation g loc = 0 while the Laplacian term is non-zero. Consequently, the granular fluidity g is non-zero and is predicting a slow creeping flow over a few grain diameters where the local theory predicted a solid-like behavior. Thanks to this non-local term, we do not observe the sharp front between the static and flowing region caused by the Drucker-Prager criterion [19] in Eq. (6) anymore.

Solver implementation 3.1 Governing equations
Which one from FVM or FEM is the most suitable for granular modeling is still an open question and there is no evidence to substantiate either claim. Hence, our choice of using FVM instead of FEM is based on its ease of implementation and its wide use in CFD simulations. For the implementation of the NGF viscosity model, we have chosen to use the incompressible InterFoam [20] solver included in the OpenFoam software. InterFoam is a solver which makes use of the FVM discretization as well as the PIMPLE algorithm, a combination of the SIMPLE (Semi-Implicit Method for Pressure Linked Equations) and PISO (Pressure-Implicit with Splitting of Operators) algorithms [21]. Figure 1. A typical control volume (mesh cell) used for the discretization of the solution domain using FVM. The control volume represents a portion of the continuum and not a reference particle. Its size is defined by the numerical stability criterion of the system and the mesh independence of the results which may lead to volumes smaller than a particle.
Besides FVM, the choice of InterFoam has been motivated by its use of the volume-of-fluid method which changes it into a two-phase flow solver and makes it applicable to a wider range of applications than a standard single phase solver. The governing equations including the continuity equation, the momentum equation and the VOF equations are as follow: where x is the position vector, G is the gravity, p d = P − ρG · x is a modified pressure, η is the mixture dynamic viscosity, ρ is the mixture density, ρ 1 = φρ s is the bulk density of the granular phase, U r = U 1 − U 2 is the compression velocity and ∇ · U r γ(1 − γ) is the interface compression term which does not affect the solution but helps to sharpen the interface (Weller-VOF implementation [20]). It should be noted that the designation "compression" only relates to the VOF-interface smearing and does not relate to compressible flow. We take the solid fraction φ = 0.62 near the random close packing limit of quasi-monodisperse spherical grains.

Viscosity model
In addition to the governing and VOF equations, the NGF model has to be implemented in the InterFoam solver. Using the FVM method [13], the discretized form of the NGF fluidity equation (3) naturally becomes where V p is the control volume (see Fig. 1), the subscript f designates the value of the variable at the center of the face, S is the outward-pointing face area vector and (∇g) f is calculated using a gaussian linear correction scheme. Unlike (1) which is fully explicit, both g terms are treated implicitly. However, while Eq. (1) and (12) are extremely simple and straightforward they are also subject to numerical instabilities. Especially Eq. (1) which presents a ratio between the strain rate and the fluidity g which can potentially lead to an alternation of extreme values between g and µ. In order to ease the finding of both the friction coefficient and the non-local fluidity we set up a relaxation loop which allows us to calculate both fields simultaneously similar to the pressure-velocity coupling procedure.

Results
As a first validation case for our model implementation we consider an annular shear cell of height H = 10 mm as shown in Fig. 2(A). The inner cylinder of radius R i = 51 mm is rotated at a rate Ω and the outer cylinder is chosen as R o = 63 mm so that R i − R 0 = 16d with d = 0.75 mm the grain diameter, that way we minimize the influence of the outer wall on the experiment results. We consider glass beads of density ρ s = 2450 kg/m 3 which leads to µ s = 0.38, µ 2 = 0.64, I 0 = 0.279 [18] and A = 0.48 [22]. We use no slip boundary conditions for both the inner and outer walls, perfect slip for the bottom and a free surface for the top. We take n i (∂g/∂x i ) = 0 with n i the surface normal for the fluidity g while P is specified at the outer wall. The fluidity boundary condition has been chosen for its simplicity as well as its low intrusiveness as discussed in the supporting information in [22]. We compare our results using the FVM-NGF model with the FEM-NGF implementation [22] and experiments [23] in Fig. 2(B). We plot the surface tangential velocity V θ normalized by the inner wall velocity V Ri = ΩR i as a function of the normalized distance from the inner wall (r − R i )/d. The symbols represent different wall speeds and highlight the near rateindependent response in the annular shear geometry. We can observe the shear band developing from the inner wall and quickly decaying into the bulk. Our implementation is showing a good agreement with comparative data. However, our solver is not limited to single-phase granular flows but can also be used for more challenging applications such as 3D frictional heap flows as described in Fig. 3(A)&(B). Similarly to [16], we simulate a granular layer going down an incline surrounded by friction walls  with µ wall = 0.32 topped with a newtonian fluid whose density and viscosity are defined as ρ f /ρ s = 0.00625 and ν f / gH 3 = 0.1 using VOF. All the material parameters are kept identical as before while d = 0.5 mm, H = 20d, W = 19d and θ = 28 degree. Fig. 4(A) shows the two layers within the computational domain while Fig. 4(B) shows the granular layer flow field and its characteristic shear band near the free surface with the bow shape resulting from the frictional walls.

Conclusion
In conclusion, we proposed a new straightforward and user-friendly continuum framework capable of modeling non-local phenomena in single-phase and two-phase granular flows. It is based on the NGF model, FVM discretization and VOF method and has shown good agreement against experiments and FEM implementation results.