The performance of the μ(I)-rheology model on flat bottom silos discharge

The aim of this work is to explore the capability of the μ(I)-rheology model and its numerical implementation in addressing a silo discharge problem by computational simulation. In order to do so, the model was implemented in the general structure of an Eulerian multiphase solver based on the Volume-Of-Fluid (VOF) method of the OpenFOAM(R) suite. First, the implementation is validated against the results of another Lagrangian and Eulerian codes in a two-dimensional discharge problem. After that, the model is tested against the experimental results of a lab-scale and industrial-scale discharge problem. While the results of the first one were satisfactory in terms of discharge rate, for the latter one, the model exhibits disagreements in the flow patterns inside the silo. The study shows the limits of applicability of the standard formulation of the model for real scale silos and sets the ground for further discussion and improvements.


Introduction
Silos are largely used in the agricultural industry for handling and storage of granular materials. Among the many variables involved, the rate of discharge is one of the most important characteristics to consider, in search of an optimal design and performance. The behaviour of the grains during discharge has been of interest for many decades due to the wide range of regimes and phenomena involved [1]. In order to study these regimes, computational simulations come as a non-expensive and non-intrusive highly detailed tool that can complement (and sometimes even substitute) the experimental approaches. There are two main branches available to study these problems by computational means: Eulerian and Lagrangian methods. The latter seeks to represent the multiple interactions between particles in detail. However, the modelling of individual particles in systems of the scale of industrial silos makes this approach expensive to be tractable for full systems simulations [2]. On the other hand, the Eulerian approach relies on addressing averaged conservation equations on the granular media, which results in a continuum set of equations. In this approach, the scales of the system do not represent a big obstacle as for the Lagrangian techniques but accurate modelling of the granular rheology is crucial.
The flow in granular materials is complex and highly dependent on the concentration and properties of the par- * e-mail: cesarvenier@gmail.com ticles. In this sense, the rheology can be modelled by considering the behaviour of the flow [3]. First, a solid-like behaviour may be attributed to quasi-static regimes (such as piles of grains), where the soil mechanics and frictional theories are suitable for describing it [4]. On the opposite side, when the granular media is in a highly diluted state, the kinetic theory of granular flow comes into place [5]. This regime resembles the behaviour of a gas and can be found in applications such as a pneumatic transport of grains. However, the cope of both of these states does not cover the whole spectrum of regimes. During the discharge of a silo, whether the grains move in funnel or mass flow, the granular media presents a behaviour that resembles the flow of a liquid phase. This type of regime cannot be characterized by the aforementioned theories and needs of a constitutive law on its own. In this aspect, Da Cruz et al. and Jop et al. [6,7] developed the µ(I)-rheology model that seeks to model this intermediate state.
The performance of the µ(I)-rheology model and its implementation have been largely discussed for several applications in the literature. Staron et al. [8] showed the suitability of the model through a qualitative comparison against a Lagrangian approach on a two-dimensional silo. Barker et al. [9] discuss the ill-posed nature of the model for certain ranges of the inertial number, Franci et al. [10] showed the drawbacks of the standard implementation of the model and proposed two regularizations of the µ(I) function. The group led by Prof. Kamrin [11,12] worked on coupling the visco-plastic flow behavior of the granu-lar media with an elasticity law to account for the stagnant behavior when the shear stress is below the yield criteria. However, there is still a lack of discussion regarding the performance of the model for different scales of silos during discharge. It is the purpose of this work to consider these aspects and discuss the applicability of the standard model for industrial-scale silos.

Modeling and implementation
The development of the rheological model is based on the observations of Da Cruz et al. [6], where it is shown that the shear stress is proportional to the normal stress by means of a coefficient proportional to the so-called inertial number I. Later on, Forterre et al. [3] and Jop et al. [7] formulate the model for three dimensions, which results in: where η, p,γ, d P and ρ s are the granular phase viscosity, the pressure field, the second invariant of the shear stress tensor, the particle diameter and the solids bulk viscosity, respectively. For this model, the µ(I) function is defined as : where µ 1 , µ 2 and I 0 are material parameters of the model. The implementation of the model has been done on the OpenFOAM(R) platform [13] based on the general structure of the viscosity models, and the simulations are performed using an Eulerian multiphase solver based on the Volume-Of-Fluid (VOF) method. As mentioned by Jop et al. [7], the yield criterion is reached when |γ| tends to zero. However, since the implementation is made on a flow solver structure, η cannot be allowed to diverge. To adress this issue, some authors proposed to use elasto-plastic models to fulfill the yield criteria and to correctly model static conditions [11,12].
In this work, a maximum viscosity threshold is included to avoid singularities as done by Popinet et al. [14]. In this way, the granular material would always behave as a fluid and stagnant zones can only be replicated as regions with a very high effective viscosity.

Two-dimensional silo discharge
First, the implemented model is validated by solving a two-dimensional silo discharge as proposed by Starton et al. [8]. The results are presented in a non-dimensional form for different discharge conditions over a silo with a flat bottom, where W represents the non-dimensional discharge rate (as defined by Staron et al.) and d represents the outlet size.    Fig. 2 shows the pressure field using both Eulerian softwares and DEM. A qualitative agreement is observed between both Eulerian tools with slight differences near the lateral walls which might be related to the boundary condition adopted. Fig. 3 shows the velocity profiles on a horizontal line (at y = 0.15H, where H is the initial filling). Here, a good agreement is found between both Eulerian software, with results that are closer between OpenFOAM and DEM for the horizontal components of the velocity in the transition region between discharge and the quasi-static region.

Three-dimensional silo discharge
This section focuses on three-dimensional cylindrical (flatbottomed) silos discharge. Two experiments with different dimension are considered to assess the performance of the rheology model. The first one is based on the work of Uñac et al. [15], which will be referred to as Test 1, and the second one is based on the experimental results of Chen et al. [16], which will be referred to as Test 2. Unlike Test 1, Test 2 consists of an industrial-scale silo with an excentric discharge (ε) where the outlet flow is controlled by a discharge feeder. A scheme of the problem can be observed in Fig. 4 and the dimensions and physical parameters for Test 1 and 2 are presented on Table 1. The simulations were performed using 64.800 cells on an O-grid discretization for both tests. The time step was adjusted to guarantee stability and the PISO-SIMPLE algorithm was set up to perform 10 SIMPLE iterations and 5 PISO iterations which ensured a decrease of 6 orders of magnitude of the pressure residuals per time step.
For Test 1, a set of simulations, each one with the same conditions but with a different outlet diameter, were performed (d = 3, 4.5, 5.9, 7.5, 9, 10.5 and 12 cm). Once a steady outflow was obtained (after the initial transient state where the fields change from stagnant to discharge conditions), the rate of discharge is compared with the experimental results. This comparison is shown in Fig. 5 along with a linear fitting of the discrete points.
In these conditions, it is expected that the discharge rate takes the form of the Beverloo scaling for flatbottomed silos: W = C (d − f ) 2.5 , where C and f are adjustment parameters of the correlation. The results shown in Fig. 5, indicate that a linear correlation matches the simulated results with high accuracy (with a coefficient of determination R 2 > 0.99). It is worth to mention that the standard form of the µ(I) model does not take into account particle-scale effects (e.g. blockage during discharge). Therefore, as the diameter of the outlet decreases, the discharge rate should tend to zero. This suggest that the CFD results should lose the linear trend given by the fitting in Fig. 5 as it reaches the zero diameter of the outlet. That being said, the experimental and CFD results are in very good agreement for the studied conditions, with a tendency of the CFD to overpredict the rate of discharge in comparison to the experimental results as the outlet diameter tends to increase. Fig. 6 shows the main fields distribution on a mid vertical plane (showing only half of the plane up to the axis line due to the high symmetry of the results) at t = 10 seconds after the discharge begins. The qualitative evolution of all the fields obtained through CFD simulations agrees with the expected behaviour, showing a transition from a mass flow to a funnel flow discharge (as represented through the solids fraction α s ). Also, stagnant zones in the circular corner near the base, where the pressure and the solids viscosity reach their maximum values, and high gradients of the velocity field close to the orifice are observed. On the other hand, for Test 2, Fig. 7 shows that the velocity magnitude inside the silo predicted by CFD differs from the experimental measurements. Although there the general behaviour is correctly reproduced, there are a few aspects that deserve further attention. The excentric discharge emphasizes the need of accurate modelling of the wall boundary conditions to account for the grains dynamic in contact with the wall (Fig. 7 shows the results for a partial-slip boundary condition with a blending factor of 0.5). This becomes clear by comparing the velocity contours near the cylindrical wall. Moreover, the experimental data indicates that there is a larger stagnant region moving away from the discharge in com-parison to the CFD prediction. These differences become relevant when the interest of the study is focused on, for example, the mixing patterns of grains inside a silo during discharge. In this situation, an accurate representation of the velocity field inside the silo becomes crucial.Another way to modify the flow patterns is to change the thresholds of maximum allowable viscosity, which acts as an artificial adjustment on the capability of the material to flow or to act more rigid. An elastic law inclusion for solidlike behaviors below the yield criterion [11,12] could also narrow the size of the non-stagnant regions in the simulation. At the current state of development of the rheological model, preliminary studies, consisting of parameters adjustment and validation against experiments, cannot be omitted for the simulation of industrial-scale silos.

Conclusions
The performance of the µ(I)-rheology model to address a flat-bottom silo discharge by computational simulation was studied. The model was implemented in the framework of a VOF solver on the OpenFOAM(R) platform. The model was assessed for a two-dimensional silo discharge obtaining a good agreement against the reported results using DEM and Eulerian approaches in terms of pressure distribution and velocity fields. Then, two threedimensional cases of different dimensions were tested and the discharge rates and flow distribution were compared against the experimental results. For the lab-scale test, the CFD results show a field distribution that correlates well with the expected behaviour and a good agreement with the experiments. On the other hand, the velocity profiles predicted by CFD inside the silo do not agree with the experimental observations on an industrial-scale silo with a controlled discharge rate. The reasons behind this may be related to the boundary condition adopted and the rheological parameters involved. In any case, the results indicate that, at this stage of development of the model, industrialscale experimental data is essential to assess and complement the computational approach.