Julien Chauchat

Marc Médale

A three-dimensional numerical model for dense granular flows based on the µ(I) rheology

Keywords: Dense granular flows, Visco-plastic flows, Regularisation techniques, Finite Element Method

This paper presents a three-dimensional implementation of the so-called µ(I) rheology to accurately and efficiently compute steady state dense granular flows. The tricky pressure dependent visco-plastic behaviour within an incompressible flow solver has been overcome using a regularisation technique along with a complete derivation of the incremental formulation associated with the Newton-Raphson algorithm. The computational accuracy and efficiency of the proposed numerical model have been assessed on two representative problems that have an analytical solution.

Then, two application examples dealing with actual lab experiments have also been considered: the first one concerns a granular flow on a heap and the second one deals with the granular flow around a cylinder. In both configurations the obtained computational results are in good agreement with available experimental data.

Introduction

Dense granular flows are ubiquitous in geophysics, e.g. landslides, avalanches, debris flows, but also in many industrial processes, e.g. pharmaceutical production, chemical industry, food and agricultural industries, energy production. Over the last two decades, an important research effort [START_REF] Jaeger | Granular solids, liquids, and gases[END_REF][START_REF] De Gennes | Reflections on the mechanics of granular matter[END_REF][START_REF] Pouliquen | Scaling laws in granular flows down rough inclined planes[END_REF][START_REF] Gdrmidi | On dense granular flows[END_REF][START_REF] Da Cruz | Rheophysics of dense granular materials: Discrete simulation of plane shear flows[END_REF] has led to the formulation of a dense granular flow rheology (µ(I)) [START_REF] Forterre | Flows of dense granular media[END_REF] that allows to describe most of the observations from different flow configurations [START_REF] Pouliquen | Scaling laws in granular flows down rough inclined planes[END_REF][START_REF] Gdrmidi | On dense granular flows[END_REF][START_REF] Jop | Crucial role of side walls for granular surface flows: consequences for the rheology[END_REF][START_REF] Forterre | Kapitza waves gives support for three-dimensional granular flow rheology[END_REF]. More recently, Jop et al. [START_REF] Jop | A constitutive law for dense granular flows[END_REF] have proposed a tensorial formulation of this rheology that opens the way to three dimensional simulations of dense granular flows in complex configurations. The literature survey shows only few publications concerning the implementation of the dense granular flows rheology in a multidimensional numerical model [START_REF] Chauchat | A 3D numerical model for incompressible two-phase flow of a granular bed submitted to a laminar shearing flow[END_REF][START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes model with a mu(I)-rheology[END_REF][START_REF] Staron | The granular silo as a continuum plastic flow: The hour-glass vs the clepsydra[END_REF]. However, such tools are awaited by industrial and geophysical communities working with granular matter.

The dense granular flow rheology is of frictional nature meaning that the tangential stresses are proportional to the normal ones, i.e. pressure, by an effective friction coefficient denoted as µ. From a rheological point of view it exhibits a visco-plastic behaviour with a non-conventional shear thinning character. It has been demonstrated from the dimensional analysis of plane shear for infinitely rigid particles that the only dimensionless parameter that controls the Preprint submitted to Journal of Computational Physics September 11, 2013 stress/shear rate relationship is the number I = d γ / p/ρ, designated as the inertial parameter [START_REF] Da Cruz | Rheophysics of dense granular materials: Discrete simulation of plane shear flows[END_REF], in which γ is the modulus of the shear rate, p is the pressure exerted on the granular media, d and ρ represent the particles diameter and density, respectively. From the analysis of discrete numerical simulations [START_REF] Da Cruz | Rheophysics of dense granular materials: Discrete simulation of plane shear flows[END_REF] and experimental measurements [START_REF] Pouliquen | Scaling laws in granular flows down rough inclined planes[END_REF] a master curve for the evolution of the friction coefficient µ against the inertial parameter I emerges, showing that this dimensionless number is the relevant parameter controlling the mechanical behaviour of the dense granular flow in the liquid regime [START_REF] Forterre | Flows of dense granular media[END_REF]. This master curve has been fitted using a rational functional form (equation 4 in section 2) that corresponds to the dense granular flows rheology µ(I). This local rheology has since been applied with some success to different flow configurations, e.g. inclined plane, flow on a pile, rotating drums, however it should be kept in mind that some limitations still remain [START_REF] Forterre | Flows of dense granular media[END_REF]. They are mainly linked to non-local effects such as correlated motions with typical length scales greater than the particles diameter or in the gaseous regime in which particles velocity fluctuations can be transported by the mean flow.

Concerning the numerical implementation of models related to dense granular flows, Chauchat and Médale [START_REF] Chauchat | A 3D numerical model for incompressible two-phase flow of a granular bed submitted to a laminar shearing flow[END_REF] have considered a Coulombian rheology (µ=constant) to study the configurations of an immersed granular media submitted to a laminar shearing flow. Their three-dimensional Finite Element steady-state model was built using a penalisation method and a simple regularisation technique, leading to quite slow convergence rates, despite using an incremental formulation in the Newton-Raphson algorithm. Moreover, two quite restrictive assumptions were also made in the previous work, considering fixed free surface location and a constant particle volume fraction inside the bed layer. On the other hand, Lagrée et al. [START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes model with a mu(I)-rheology[END_REF] have used a two-dimensional explicit in time Finite Volume solver (projection method) and a simple viscosity clipping technique to deal with the µ(I) rheology. Using a two-phase approach coupled with a Volume Of Fluid technique to capture the mesoscopic air-granular interface evolution the authors have successfully applied their model to study the collapse of a dry granular column under gravity effect.

Both models have been validated against analytical solutions demonstrating the good implementation of the dense granular flow rheology. However each model has its own main weaknesses: fixed free-surface and constant granular concentration assumptions for Chauchat and Médale model [START_REF] Chauchat | A 3D numerical model for incompressible two-phase flow of a granular bed submitted to a laminar shearing flow[END_REF] and weak computational efficiency to reach steady state solutions involving rigid zones for Lagrée et al. model [START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes model with a mu(I)-rheology[END_REF]. It should also be pointed out that the latter numerical model would experience computational issues to simulate flow configurations involving re-circulation zones, since as like most standard incompressible flow solver it will generate negative dynamical pressure in these locations, meanwhile using an absolute value of the pressure in the inertial number definition. This point will be further discussed in the present paper where a physically and computationally consistent solution is proposed.

The main objective of the present work is to accurately and efficiently compute three-dimensional dense granular flows based on the µ(I) rheology. To reach this goal, a consistent and efficient regularisation technique has been implemented in our finite element framework (section 2). Extensive validations and spatial convergence analyses have been performed on two idealised test cases: the vertical chute flow and the flow down a rough inclined plane (section 3). The capabilities of the present model are then illustrated on two application examples (section 4). The granular flow on a heap with sidewalls is presented and compared with existing experimental measurements [START_REF] Jop | Crucial role of side walls for granular surface flows: consequences for the rheology[END_REF] and numerical simulations [START_REF] Jop | A constitutive law for dense granular flows[END_REF]. Finally, the granular flow around a cylinder is presented and the evolution of the drag coefficient against the flow velocity are discussed.

Model formulations

Mathematical model

In the present model an idealised granular material constituted with mono-disperse spherical beads of diameter d and density ρ is considered. Given a cartesian coordinate system (O, x, y, z) where x represents the horizontal direction, y the spanwise direction and z the upward direction (x and z being likely tilted with an angle α from the horizontal and vertical direction, respectively), the velocity vector and its cartesian components are denoted by -→ u = (u, v, w), respectively. The continuity and momentum equations read:

- → ∇ • - → u = 0, (1) ρφ
∂ - → u ∂t + - → ∇ • - → u ⊗ - → u = - --→ ∇p + - → ∇ • τ + ρφ - → g , (2)
where p and τ represent the pressure and the deviatoric part of the stress tensor (interganular stresses due to particleparticle interactions), -→ g is the gravity acceleration vector and φ denotes the volume fraction of particles. The latter is assumed to be uniform and constant (φ = 0.6) in space and time in the present model. These governing equations ( 1) and ( 2) are supplemented with initial and boundary conditions that will be presented for each case in section 3 and 4.

Following Jop et al. [START_REF] Jop | A constitutive law for dense granular flows[END_REF] and Forterre and Pouliquen [START_REF] Forterre | Flows of dense granular media[END_REF] a frictional rheology based on the µ(I) law is assumed for the modelling of the deviatoric part of the stress tensor:

τ = µ(I) p γ γ , (3)
where the rate of strain tensor γ is defined as γ = ∇ -→ u + (∇ -→ u ) T and its magnitude is given by the square root of its second invariant γ = 1 2 T r γ2 . In equation (3), µ(I) represents the phenomenological friction coefficient that depends on the inertial parameter I according to the following relationship:

µ(I) = µ s + ∆µ I 0 /I + 1 and I = γ d p/ρ , (4)
and ∆µ = µ 2µ s . In relationship (4), µ s stands for the static friction coefficient or the so-called tangent of the repose angle, µ 2 is a dynamical friction coefficient and I 0 is an empirical coefficient of the rheology. These two last parameters have been determined from experiments and discrete numerical simulations. Their respective values will be given for each case in the following sections.

Following our previous work on dense granular flows [START_REF] Chauchat | A 3D numerical model for incompressible two-phase flow of a granular bed submitted to a laminar shearing flow[END_REF], one can define a particulate viscosity η p as:

η p = µ(I) p γ , (5)
such that the shear stress tensor can be rewritten in an usual way for fluid mechanics problems:

τ = η p γ. ( 6
)
In the following, the model equations are made dimensionless using d as the length scale, d/g as the time scale and φρgd as the stress scale:

             - → ∇ • - → u = 0 ∂ - → u ∂t + - → ∇ • - → u ⊗ - → u = - - → ∇ p + - → ∇ • η p ∇ - → u + ∇ - → u T -cos α - → e z + sin α - → e x , (7)
where the bar notations represent dimensionless variables and η p = µ(I) p/ γ with I = γ / φp. We point out here that the final expression of I depends on the chosen dimensionless form of the equations. For the sake of clarity, the bar notations are dropped in the following and all the variables correspond to their dimensionless counterparts.

Finite Element Model

In this section, a numerical model designed to solve for the coupled set of equations [START_REF] Carey | Finite Elements[END_REF] in arbitrary three dimensional computational domains is presented. The Finite Element Method (FEM) is used to build up the variational formulation and the related spatial discretisation in view of a tractable algebraic system. Let us denote the computational domain by Ω and its boundary by ∂Ω. From the set of partial differential equations [START_REF] Carey | Finite Elements[END_REF] the weak formulation of the coupled mass conservation and momentum equations reads [START_REF] Carey | Finite Elements[END_REF], [START_REF] Pironneau | Méthodes des éléments finis pour les fluides[END_REF]:

Find - → u h ∈ U h and p h ∈ Q h satisfying (8) ∀ -→ δu ∈ V h and ∀ δp ∈ Q h where: U h = { - → u h ∈ H 1 (Ω) | - → u = - → u Dirichlet on ∂Ω Dirichlet } Q h = {δp ∈ L 2 (Ω)} V h = { -→ δu ∈ H 1 0 (Ω) | -→ δu = - → 0 on ∂Ω Dirichlet }                            Ω δp - → ∇ • - → u h -λ p p h dΩ = 0 , Ω -→ δu • D - → u h Dt dΩ - Ω p h - → ∇ • -→ δu dΩ + Ω η p 2 δγ : γh dΩ + Ω -→ δu • cos α - → e z -sin α - → e x dΩ - ∂Ω T -→ δu • - → T M, - → n ds - ∂Ω OBC -→ δu • -p h 1 + η p γh - → n ds = 0 , (8)
where

λ p << 1 is the penalisation parameter, δγ = ∇ -→ δu + ∇ -→ δu T , γh = ∇ - → u h + ∇ - → u h T ,
(∂Ω OBC = ∂Ω -∂Ω Dirichlet -∂Ω T )
. This later condition is encountered when the computational domain is arbitrarily truncated with respect to the physical one, so that an appropriate numerical or open boundary condition has to be provided implicitly [START_REF] Papanastasiou | A new outflow boundary condition[END_REF] to behave as passively as possible.

Regularisation technique

As reported in [START_REF] Chauchat | A 3D numerical model for incompressible two-phase flow of a granular bed submitted to a laminar shearing flow[END_REF] the visco-plastic characteristic of the Coulomb's granular rheology can be straightforwardly implemented using a simple regularisation technique. In the present paper four different regularisation techniques have been investigated to appropriately implement the µ(I) rheology in a computationally suitable, accurate and efficient way. In the following, the particulate viscosity ( 5) is re-written in different regularised ways where λ r << 1 stands for the regularisation parameter:

• The simple regularisation presented in Frigaard and Nouar [START_REF] Frigaard | On the usage of viscosity regularisation methods for visco-plastic fluid flow computation[END_REF], transposed to the µ(I) rheology (denoted as simple regularisation in the following):

η s p = µ s p γ +λ r + ∆µ p I 0 √ φp+ γ +λ r (9)
• The mixed Bercovier-Engelman [START_REF] Bercovier | A finite-element method for incompressible non-newtonian flows[END_REF] -simple regularisation, transposed to the µ(I) rheology (denoted as Bercovier-Engelman regularisation in the following):

η be p = µ s p γ 2 +λ 2 r 1/2 + ∆µ p I 0 √ φp+ γ +λ r . (10)
• The mixed Papanastasiou's [START_REF] Papanastasiou | Flows of materials with yield[END_REF] -simple regularisation, transposed to the µ(I) rheology (denoted as Papanastasiou regularisation in the following):

η papa p = µ s p 1 -e -γ /λ r γ + ∆µ p I 0 √ φp+ γ +λ r . (11)
• The proposed regularisation, inspired from the Bercovier-Engelman's regularisation [START_REF] Bercovier | A finite-element method for incompressible non-newtonian flows[END_REF] and transposed to the µ(I) rheology:

η mc p =         µ s + ∆µ γ I 0 √ φp + γ         p γ 2 +λ 2 r 1/2 . ( 12
)
The first term of the latter relationship represents the yield stress that determines the transition between rigid and flowing regions, whereas the second one can be viewed as a Non-Newtonian viscous contribution to the stress. This last term diverges only if p and γ vanish together. This condition can be met, for example, at the free surface of a granular flow where the pressure vanishes. It should also be noted that both the original Papanastasiou's regularisation [START_REF] Papanastasiou | Flows of materials with yield[END_REF] and its transposition to the µ(I) rheology [START_REF] Da Cruz | Rheophysics of dense granular materials: Discrete simulation of plane shear flows[END_REF] are not definite when γ is strictly equal to zero.

Non-linear solution algorithm

As we are mainly interested in steady state solutions the advective terms and visco-plastic stress ones can induce highly non-linear behaviour in the governing equations. For this purpose, a linearising stage based on the Newton-Raphson algorithm has been introduced. Its incremental form is obtained by taking the first variation of the variational formulation [START_REF] Chauchat | A 3D numerical model for incompressible two-phase flow of a granular bed submitted to a laminar shearing flow[END_REF]. Owing to the presented regularisation techniques, a consistent tangent matrix can be derived for each of them, leading to a high convergence rate of the steady state solution algorithm. The incremental form associated to the proposed regularisation technique (12) reads:

I ′ = Ω -→ δu • - → u ′ h ∇ - → u h + - → u h ∇ - → u ′ h dΩ - Ω p ′ h - → ∇ • -→ δu dΩ + Ω 1 2 δγ : η mc p γ ′ + η ′ mc p γ dΩ - ∂Ω OBC -→ δu • -p ′ h 1 + η mc p γ ′ + η ′ mc p γ - → n ds ( 13
)
where

γ ′ = ∇ - → u ′ + ∇ - → u ′T ; γ ′ = γ : γ ′ 2 γ (14)
and

η ′ mc p = η mc p p - 1 2 A γ p ′ -          η mc p γ γ 2 +λ 2 r -A p          γ ′ (15)
with A = ∆µI 0 √ φp

I 0 √ φp + γ 2 γ 2 +λ 2 r 1/2
In all the previous relationships the prime superscript stands for an incremental quantity computed at every iteration within the Newton-Raphson solution algorithm.

It is also noteworthy that the pressure term that appears in both η mc p and η ′ mc p [START_REF] Forterre | Kapitza waves gives support for three-dimensional granular flow rheology[END_REF] and ( 15) is the total pressure p tot = p hst + p dyn , i.e., the sum of hydrostatic p hst and dynamic p dyn pressures. As in the present model the granular media is assumed to be a continuum and the flow to be incompressible, the computed dynamical pressure is the one that satisfies the incompressibility constraint. Therefore, as in almost all incompressible flow solvers, it is not positive definite. On the other hand, the particulate viscosity must remain positive from basic physical consideration. Consequently, the total pressure appearing in the definition of the particulate viscosity ( 13) and [START_REF] Frigaard | On the usage of viscosity regularisation methods for visco-plastic fluid flow computation[END_REF] has to be regularised as well. In the present work we have introduced the following regularised expressions for that purpose:

p regu = 1 2 p tot + p 2 tot + λ 2 r and p ′ regu = p ′ dyn 2             1 + p tot p 2 tot + λ 2 r             . (16)
Therefore, the regularised pressure entering the µ(I) rheology is restricted by below to λ r when the total pressure becomes negative.

Applying this approach all along with the regularising and linearising stages, round-off errors, if any, take place consistently in both Jacobian (tangent matrix) and residual computations. Using this method enables to achieve stable and accurate computations even from such unfavourably ill-conditioned algebraic systems.

Implementation

The spatial discretisation of both the weak integral form [START_REF] Chauchat | A 3D numerical model for incompressible two-phase flow of a granular bed submitted to a laminar shearing flow[END_REF] and its first variation ( 13) is performed with piecewise quadratic polynomial approximation for the velocity and piecewise linear discontinuous approximation for the pressure. They have been implemented with 27-nodes iso-parametric hexahedra Lagrange finite elements built on tri-quadratic Lagrange polynomials (H27) [START_REF] Medale | A parallel-computer implementation of the Asymptotic Numerical Method to study thermal convection instabilities[END_REF][START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of Asymptotic Numerical Method in the vicinity of bifurcations[END_REF].

The in-house research code has been implemented using the PETSc library [START_REF] Balay | PETSc Web page[END_REF][START_REF] Balay | PETSc users manual[END_REF], which provides various parallel iterative solvers and several interfaces to third party parallel direct solvers. This last feature is particularly welcomed in the present problem as, on the one hand the incompressibility constraint is dealt with a penalisation method and on the other hand the visco-plastic character of the dense granular flows rheology is dealt with regularisation techniques.

As a consequence, all the algebraic systems are highly ill-conditioned, so they have been solved by the MUMPS parallel direct solver [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF][START_REF] Amestoy | Multifrontal parallel distributed symmetricand unsymmetric solvers[END_REF][START_REF] Amestoy | Hybrid scheduling for the parallel solution of linear systems[END_REF]. Moreover, the pressure term that appears in ( 12) and ( 15) is implemented according to [START_REF] Gdrmidi | On dense granular flows[END_REF].

Validations

Several configurations representative of the three-dimensional µ(I) dense granular flow rheology can be considered for the numerical verifications and validations of the four implemented regularisation techniques, but the two that have been retained hereafter are in addition described by analytical solutions. The first one consists of the vertical chute flow problem for which an original analytical solution is derived and used to perform a spatial convergence analysis.

The second test case corresponds to the dry granular flow over a rough inclined plane that corresponds to the so-called Bagnold's profile.

Dry granular vertical-chute flow

For the purpose of testing the numerical model the streamwise uniform flow between two rough parallel plates is considered (figure 1). This flow is a simplified configuration representing the flow in a silo. An analytical solution for the dry granular vertical-chute flow based on the µ(I) rheology [START_REF] Jop | Ecoulements granulaires sur fond meuble[END_REF][START_REF] Forterre | Flows of dense granular media[END_REF] for the frictional shear stress can be obtained. All the details concerning the derivation of the analytical solution are presented in appendix A and the resulting velocity profile expression reads as follows:

w(y) =                                      -I 0 √ φp                 y - b 2 + ∆µ p ln                 ∆µ p -(y -y 0 ) ∆µ p - b 2 -y 0                                 i f y ∈ [y 0 ; b/2] -I 0 √ φp                 y 0 - b 2 + ∆µ p ln                 ∆µ p ∆µ p - b 2 -y 0                                 i f y ∈ [0; y 0 ] (17)
where y 0 = µ s p. It should be mentionned that this analytical solution exhibits some physical limitations in the quasi-static limit (I → 0) [START_REF] Forterre | Flows of dense granular media[END_REF] since it predicts zero thickness of the shear bands whereas they have finite size in experimental observations. Therefore, in the following test case we have chosen parameter values to ensure the granular flow to be in the liquid regime.

The spatial convergence analysis of the numerical solution is presented hereafter for this vertical-chute flow with a particulate pressure fixed to p = 28 and a dimensionless domain width set to b = 28. The governing equations are solved with a penalisation parameter λ p = 10 -6 when not otherwise mentioned and the rheological parameters are fixed to µ s = 0.383, ∆µ = 0.26 and I 0 = 0.279, following the values reported by Jop et al. [START_REF] Jop | Crucial role of side walls for granular surface flows: consequences for the rheology[END_REF].

In the first part of this analysis the regularisation parameter is set to λ r = 10 -4 and the corresponding numerical results have been computed for the four regularisation techniques presented in section 2. Figure 2 shows the comparison of the numerical results with the analytical solution for the velocity, inertial parameter, shear rate and tangential stress profiles for a uniform mesh of size h = b/80. The very good agreement obtained for this flow configuration enables us to validate the implementation of the µ(I) rheology with the four regularisation methods presented in section 2.

Figure 3a) plots the Root Mean Square (RMS) numerical error with respect of the analytical solution versus the mesh size h for the four regularisation techniques. The simple regularisation method appears to be the worst one with the highest errors and the worst spatial convergence. On the other hand, the Bercovier-Engelman's regularisation, the Papanastasiou's one and the one introduced in this paper exhibit similar convergence properties for coarse and medium mesh sizes, but for the finest meshes the proposed regularisation produces an error one order of magnitude lower than the three others. zone, where a numerical creeping flow appears instead of a perfectly rigid plug. The magnitude of this artificial creeping flow is of order of the regularisation parameter, which has been verified on velocity profiles at different iteration numbers during the computations (not shown here).

In order to understand the deviation of the global error between the Bercovier-Engelman, the Papanastasiou and the proposed regularisation, the same profiles as in figure 2 are presented in figure 4 for the finest mesh (h = b/320). It appears that, for this mesh size, except the proposed regularisation technique all the others present spurious oscillations on the velocity gradient and a significant deviation of the computed velocity profile compared with the analytical solution. The error increase observed in figure 3a) is linked to these oscillations. This result supports the remark of subsection 2.4 advocating that the proposed methodology enables to achieve stable and accurate computations.

In the second part of this analysis, the convergence of the proposed regularisation method with respect to the regularisation parameter λ r is analysed. Figure 5a) shows the evolution of the RMS error versus the mesh size h. For sufficiently fine meshes a decrease of λ r reduces substantially the error. This is better illustrated in figure 5b) where the evolution of the numerical solution error against the regularisation parameter is presented. It clearly appears from this figure that for low spatial resolution the global error is dominated by the discretisation error. When the spatial resolution is high enough the global error is dominated by the regularisation of the dense granular rheology and consequently the regularisation parameter strongly influence the error. For h = b/320 it is observed that the error decreases linearly with the regularisation parameter. The simulations carried out with a penalisation parameter λ p one order of magnitude smaller than the previous one does not induce any improvement on the error. This results from the fact that the incompressibility constraint is not very restrictive in such flow configurations where the pressure field is quite uniform. However, it is expected that the penalisation parameter will have a comparable influence on the error as the regularisation parameter for complex flow configurations with non negligible dynamical pressure.

s -1 / ! p =10 -6 ! r =10 -5 s -1 / ! p =10 -6 ! r =10 -6 s -1 / ! p =10 -6 ! r =10 -5 s -1 / ! p =10 -7
" h 2

Dry granular flow on an inclined plane: Bagnold's profile

In this subsection, the numerical model is used to compute the dry granular flow over a rough inclined plane.

The analytical solution is the so-called Bagnold's profile where the frictional granular shear stress is modelled using the µ(I) rheology [START_REF] Gdrmidi | On dense granular flows[END_REF][START_REF] Jop | A constitutive law for dense granular flows[END_REF][START_REF] Forterre | Flows of dense granular media[END_REF]. The analytical solution is recalled first and the spatial convergence analysis of the numerical solution is presented afterwards.

Considering the steady flow of a granular material over an infinite rough inclined plane only the longitudinal velocity component u does not vanish and depends only on the vertical upward direction (z = 0 at the bottom and z = H at the free surface). This problem has an analytical solution, the so-called Bagnold's profile, for which the detailed derivation is given in appendix B. The resulting velocity profile expression is only recalled here:

u(z) = 2 3 I α φ 0 cosα H 3/2 -(H -z) 3/2 , ( 18
)
where I α is given by the following relationship:

I α = I 0 tanα -µ s µ s + ∆µ -tanα . ( 19
)
The spatial convergence analysis of the numerical solution against the Bagnold's profile is presented hereafter for a tilt angle set to α = 0.406 rad (≈ 23 degree) and a dimensionless granular layer thickness set to H = 14. The rheological parameters are: µ s = 0.383, ∆µ = 0.26 and I 0 = 0.279, as in Jop et al. [START_REF] Jop | Crucial role of side walls for granular surface flows: consequences for the rheology[END_REF]. The governing equations are solved with a penalisation parameter λ p = 10 -6 and the µ(I) rheology is regularised using the proposed method [START_REF] De Gennes | Reflections on the mechanics of granular matter[END_REF] with a regularisation parameter set to λ r = 10 -8 .

Figure 6 shows the validation of the numerical model against the Bagnold's analytical solution in terms of streamwise velocity component, vertical velocity gradient and shear stress profiles with h = H/40. The very good agreement obtained with respect to the analytical solution validate qualitatively and quantitatively the numerical model. Figure 7 shows the spatial convergence analysis obtained on six successively refined meshes h ∈ {H/5, H/10, H/20, H/40, H/80, H/160} expressed in RMS error and infinity norm. The best fit obtained on the four coarsest meshes gives an effective order of convergence of 2.46 for the RMS error and 1.5 for the infinity norm. Despite, this could be thought to be far from the optimal theoretical values for a smooth C 0 problem (3 and 2, respectively), it is however not so bad for a penalised and regularised problem. The slightly lower values obtained with the model on the present test case are almost certainly due to the regularisation of the µ(I) rheology and agree perfectly with previous studies [START_REF] Chauchat | A 3D numerical model for incompressible two-phase flow of a granular bed submitted to a laminar shearing flow[END_REF][START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes model with a mu(I)-rheology[END_REF].

The numerical results presented in figure 7 have been computed with a requested residual per degree of freedom of 10 -12 . The corresponding Newton-Raphson iteration counts and CPU time per iteration are summarised in Table 1.

Figure 8 shows the evolution of the velocity at the surface for both numerical and analytical solutions as a function of the tilt angle α where the analytical solution is given by (31) and (32) with z = H. The numerical model is able to reproduce fairly well the evolution of the surface velocity for a wide range of slope angle and more importantly the numerical model is able to capture the threshold of motion at α c ≈ 0.37 rad.

Discussion

The verification and validation of our numerical model that implements the µ(I) rheology have been performed on two representative problems that possess an analytical solution. These two test cases are relevant for benchmarking the µ(I) rheology implementation as they contain most of the specific features exhibited by dense granular flows: onset of motion, plastic-rigid transition, shear-dependent rheology, etc. Among the four regularisation techniques presented in section 2.3 the simple one always performs the worst, both for the global error at convergence and the effective order of convergence. The Papanastasiou, the Bercovier-Engelmann and the proposed regularisation techniques produce, approximately, all the same solution error at convergence for appropriate mesh size and regularisation parameter.

However, contrarily to the others the Papanastasiou regularisation does not bear neither an initial null velocity field nor any rigid motion flow, so we arbitrarily initialised it with the solution provided by the simple regularisation technique.

The proposed regularisation method [START_REF] De Gennes | Reflections on the mechanics of granular matter[END_REF] performs the best in achieving small global errors at the lowest computational cost (Newton-Raphson iteration counts at convergence). The results obtained in these two test cases evidence that the proposed regularisation method is the most consistent way to regularise the µ(I) rheology in our Finite Element framework, as it enables to achieve the most stable and accurate computations from such unfavourably ill-conditioned algebraic systems. The various sensitivity analyses performed in this paper reveal that it exists an optimal value of the combination between the regularisation parameter and mesh size (∝ λ r h 3 ) below which any subsequent reduction leads to intractable ill-conditioned algebraic systems that deteriorates the overall solution quality, even with a robust direct solver such as MUMPS [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF][START_REF] Amestoy | Multifrontal parallel distributed symmetricand unsymmetric solvers[END_REF][START_REF] Amestoy | Hybrid scheduling for the parallel solution of linear systems[END_REF]. This explains why the actual order of spatial convergence observed in these two test cases always deteriorates for the smallest mesh sizes for a given value of the regularisation parameter.

Unfortunately, the two analytic solutions used herein are only one-dimensional, as only one velocity component is not null and only varies along one spatial direction (the transverse one). However, it should be reminded that three-dimensional analytic solutions are extremely few for fluid flows, even for the Newtonian constitutive law, which is linear. So, one can figure out how cumbersome could be the work of designing a full three-dimensional analytic solution that integrates the highly non-linear µ(I) rheology. Nevertheless, based on the previous analysis more complex three-dimensional flow configurations can be computed using the proposed regularisation technique for which the solution accuracy is controlled by the numerical constant: ε h r ∝ λ r h 3 .

Applications

In this section, the numerical model presented in section 2 is applied to two configurations: the surface granular flow on a heap and the granular flow around an immersed cylinder. These two applications illustrates the capabilities of the proposed model to compute dense granular flows in rather complex configurations.

Dry granular flow on a heap

The surface flow of a dry granular material on a heap has been investigated experimentally by Jop et al. [START_REF] Jop | Crucial role of side walls for granular surface flows: consequences for the rheology[END_REF] (see figure 9). In their experiments the authors measured four quantities: the mean flow rate per unit of width Q (m 2 .s -1 ), the inclination of the free surface with the horizontal α, the thickness of the flowing layer h(y), and the free-surface The proposed numerical model is used to simulate this configuration using 70x100 H27 finite elements in the y and z directions respectively, for symmetry reasons only one half of the domain is discretised. The total height of the domain is fixed to H/d = 50. Boundary conditions are imposed as follows: no-slip at the bottom, no-slip at the lateral wall, symmetry condition at the vertical median plane, and inflow-outflow conditions at the inlet and outlet respectively. The penalisation parameter is fixed to λ p = 10 -7 and the regularisation parameter is fixed to λ r = 10 -5 .

Following the methodology used by Jop et al. [START_REF] Jop | A constitutive law for dense granular flows[END_REF], the inclination angle α is imposed and the flow rate is computed a posteriori. The tilt angle is chosen from trial and errors to get a flow rate comparable with the measurements. By doing so a typical error of less than 2% is achieved. Figure 10 shows comparisons of the free-surface velocity profiles and of the thickness of the flowing layer between experiments from Jop et al. [START_REF] Jop | Crucial role of side walls for granular surface flows: consequences for the rheology[END_REF], numerical simulations from Jop et al. [START_REF] Jop | A constitutive law for dense granular flows[END_REF] and (mesh of 1x70x100 H27 FE on the half domain, λ p = 10 -7 , λ r = 10 -5 , dashed line ---), numerical solutions from Jop et al. [START_REF] Jop | A constitutive law for dense granular flows[END_REF] (grid size of 71x80 on the full domain, solid line -) and experimental results from Jop et al. [START_REF] Jop | Crucial role of side walls for granular surface flows: consequences for the rheology[END_REF] (circles •).

* =90.1 Q * =90.1 Q * =40.2 Q * =14.9 Q * =14.9 Q * =40.2 b) a)

Dry granular chute flow over a cylinder

The second application example concerns a dry granular chute flow over a horizontal cylinder. This study aims at investigating the influence of the applied external force and confining pressure onto the granular flow rate and the induced drag force acting on the cylindrical obstacle. To conduct this parametric study the cylinder diameter denoted as D is used as the reference length scale and the particles are assumed to have a dimensionless diameter of d * = d/D = 0.1 and a uniform volume fraction φ = 0.6. The computational domain is L x /D = 20 large in the crossflow direction, L y /D = 0.1 thick in the spanwise direction (quasi two-dimensional problem) and L z /D = 30 long in the vertical streamwise direction. The cylinder axis is aligned in the spanwise direction and located at the centre of the computational domain. The granular media is submitted to a driving body force acting downwards in the vertical direction, together with a confining pressure (p ext ). A no slip boundary condition ( -→ u = -→ 0 ) is imposed both along the cylinder wall and the two lateral confining walls, an outflow boundary condition is applied (implicit Neumann boundary condition [START_REF] Papanastasiou | A new outflow boundary condition[END_REF]) at the outlet, and finally a symmetry boundary condition (v = 0) is prescribed over the two other walls.

The parametric study is conducted with the proposed regularisation model (( 12), ( 14) and ( 16) with µ s = 0.4, ∆µ = 0.25, I 0 = 0.3, λ p = 10 -8 and λ r = 10 -6 ) for five confining pressures. The mesh is made up of 42200 H27 finite elements organised in a single layer (only one element) in the spanwise direction. The mesh size is graded toward the cylinder around which 160 elements are disposed with a thickness of D/50. The mesh is built onto 509220 nodes so the resulting algebraic system contains 1527660 velocity unknowns.

For every confining pressure the steady state problem solutions are computed from an uniformly null initial solu-tion. Then, the dimensionless vertical body force (F v adim = ρφgL x p ext ) is incrementally applied in a continuation procedure up to the unsteadiness limit. In order to get a deeper insight of this dense granular flows the relevant primary variable fields are plotted for a dimensionless external loading of F v adim = 0.054 in figures 13 (whole computational domain) and in figures 14 (obstacle vicinity). The regularised pressure field plotted in figure 13(a) displays the highest pressures at the upstream obstacle stagnation region and a very low pressure zone downstream the obstacle. The fields of streamwise and crossflow velocity components are plotted in figures 13(b) and 13(c), respectively. Close-up views in the obstacle vicinity (figures 14) reveal that very thin boundary layers develop in regions where the tangential velocity is the highest. Moreover, wiggles also appear in the cylinder wake, as the granular flow approaches the steady-state limit for these parameter values. Indeed, owing to the assumption of incompressible granular flow at constant volume fraction the dynamical pressure is negative and justifies the use of the regularised pressure relationship [START_REF] Gdrmidi | On dense granular flows[END_REF] in the µ(I) rheology implementation. This can be clearly seen in figure 15 where both the dynamical pressure, the total pressure and the regularized one are plotted along the left half cylinder perimeter (-π ≤ θ ≤ 0) versus the azimuthal angle (θ = 0 at the downstream stagnation point), along with the normal and tangential stresses and their regularised counterparts which are plotted on the right half perimeter (0 ≤ θ ≤ π). It turns out that within the present assumptions the dynamical pressure behind the cylinder can no longer be directly used to compute the granular viscosity in the µ(I) rheology.

# t # t (regu) # n # n (regu)

21

On the other hand, the regularised particulate viscosity [START_REF] De Gennes | Reflections on the mechanics of granular matter[END_REF] is plotted in figure 16(c) along with its static and dynamic components, plotted in figures 16(a) and 16(b), respectively, for two dimensionless external loadings of F v adim = 0.044 and 0.054. The static component spans over seven orders of magnitude between the rigid body motion far from the motionless rigid walls and the highly sheared zones close to them. On the other hand the dynamical component is maximum in the cylinder wake and varies only within one order of magnitude. Therefore the regularised particulate viscosity (figure 16(c)) tends to a very tiny value in the cylinder wake as the pressure tends to its regularised value (close to zero) in this region.

Concluding remarks

This paper presents a full three-dimensional implementation to accurately and efficiently compute steady state dense granular flows. The physical model considered in the present work is based on a continuum and incompressible approach supplied with the so-called µ(I) rheology to represent dense granular flows. The key points to get a successful implementation enabling accurate and efficient computations for this kind of models are twofolds: i) the way to deal with the tricky pressure dependent visco-plastic rheology within an incompressible flow solver; ii) the capability to solve efficiently such very stiff non-linear algebraic systems. To achieve these goals four regularised µ(I) relationships have been proposed, implemented and tested against various representative problems with analytical solutions (vertical chute flow and flow over an inclined plane). It turned out that, in our finite element framework, the one of eq. ( 12) is not only the most accurate and efficient one, but also the most stable and robust as the problem becomes extremely stiff.

Moreover, the pressure term arising in all these regularised µ(I) relationships needs also to be regularised as it is not necessarily positive definite while computed within an incompressible flow solver. Finally, the computational accuracy and efficiency of the proposed steady state solution algorithm has been obtained thanks to complete derivations of the incremental formulations associated with the Newton-Raphson algorithm.

Further to comprehensive validations on academic configurations, two application examples that exhibit the peculiarities of dense granular flows have also been considered. The first one is the granular flow on a heap which has been previously investigated experimentally and numerically elsewhere. Our computational results agree qualitatively well with both experiments and previous numerical ones, assessing the present model capability to reproduce pronounced 3D effects. The second application example deals with the vertical chute flow over a cylinder. This configuration points out the necessity to perform a pressure regularisation to obtain a computationally reliable µ(I) implementation when dealing with dense granular flows over or around obstacles. The obtained results demonstrate that the implemented model is able to qualitatively reproduce the main trends of this particular configuration (flow rate and cylinder drag force) along with some more intricate ones such as the upstream-douwnstream asymmetry induced by the strong coupling between the dynamical pressure and the granular effective viscosity.

However, a closer quantitative agreement with experiments should probably require some fitting of the numerical parameters that enter the µ(I) rheology and its implementation along with more appropriated boundary conditions at solid walls instead of the trivial no slip condition. Moreover, several other improvements will be undertaken in the near future to broaden the validity domain of the proposed model. A very interesting feature would be to account for variable particle volume fraction within the computational domain, assuming a constant particle volume fraction is no longer relevant in many actual situations. Finally, another attractive direction for industrial or geophysical applications is to model free surface flows, but the way to go is still long.

The mechanical equilibrium is simple, the shear stress balances the vertical component of the gravity and the pressure is constant across the section. The y-component of the momentum equation gives a constant distribution of the pressure: p = constant. The pressure is a free parameter that controls the velocity and the flow rate. Integrating the z-momentum equation over y gives:

τ yz = y + C ( 21
)
where C is an integration constant that should be zero for symmetry reasons. Considering the linear evolution of the shear stress [START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes model with a mu(I)-rheology[END_REF] its maximum arises at the wall y = b/2. Since the granular media exhibits a plastic threshold, given by µ s p, two cases has to be distinguished, if µ s p ≥ b/2, the shear stress τ yz is always lower than the threshold value and the granular media is not sheared accordingly with the viscoplastic granular rheology. The solution is trivial in

Introducing the expression µ(I) = µ s + ∆µ/(I 0 /I + 1) in the previous equation gives:

µ s p + ∆µp I 0 I + 1 = y. ( 23
)
Using the definition y 0 = µ s p and after some algebra the analytical profile for the inertial number I reads:

I = I 0 y -y 0 ∆µ p -(y -y 0 ) . ( 24
-I 0 φp                         y 0 -

Figure 1 :

 1 Figure 1: Sketch of the vertical-chute flow.

  Figure3b) shows the evolution of the error against the analytical solution versus Newton-Raphson iteration counts for the four regularisation techniques and h = b/160. The initial condition was set to a uniform null velocity in the whole computational domain for three out of four regularisation methods, meanwhile the Papanastasiou's regularisation was initialised by the simple regularisation solution at convergence. So, it should be pointed out that the iteration counts of the Papanastasiou's regularisation has been shifted by the iteration counts of the simple regularisation. The proposed method (12) is shown to converge faster than the Bercovier-Engelman regularisation[START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of Asymptotic Numerical Method in the vicinity of bifurcations[END_REF] and the Papanastasiou one[START_REF] Da Cruz | Rheophysics of dense granular materials: Discrete simulation of plane shear flows[END_REF] with similar global error. The simple regularisation (9) converges at the same rate as the proposed regularisation but with a global error at convergence approximately one order of magnitude higher. It is also noteworthy that the error reaches a minimum three to four times lower than the final converged value. This latter error is a modelling one introduced by all the regularisation techniques in the central

Figure 2 :

 2 Figure 2: Vertical-chute flow test case with p = 28 and b = 28. Comparison of computed (h = b/80) and analytical solution profiles: vertical velocity component w (upper left), transverse velocity gradient dw/dy (lower left), inertial parameter I (upper right) and shear stress τ (lower right).

Figure 3 :

 3 Figure 3: Vertical-chute flow test case with p = 28 and b = 28. a) Spatial convergence analysis of the four regularisation methods (Simple, Bercovier-Engelman, Papanastasiou and Médale-Chauchat) with λ r = 10 -4 , expressed in RMS error: e = 1/N N (w iw ex i ) 2 . b) RMS error versus Newton-Raphson iteration count for the four regularisation techniques (h = b/160).

Figure 4 :

 4 Figure 4: Vertical-chute flow test case with p = 28 and b = 28. Comparison of computed (h = b/320) and analytical solution profiles: vertical velocity component w (upper left), transverse velocity gradient dw/dy (lower left), inertial parmeter I (upper right) and shear stress τ (lower right).

Figure 5 :

 5 Figure 5: Vertical-chute flow test case with p = 28 and b = 28, computed with the proposed regularisation method (12). a) Spatial convergence analysis for different values of the regularisation parameter λ r = {10 -4 ; 10 -5 ; 10 -6 } and the penalisation parameter λ p = {10 -6 ; 10 -7 } expressed in RMS error: e = 1/N N (w iw ex i ) 2 . b) RMS error versus regularisation parameter λ r for five meshes h ∈ {b/10; b/20; b/40; b/80; b/160; b/320}, a penalisation parameter of λ p = 10 -6 (blue lines) and λ p = 10 -7 (red circles), the black dashed line corresponds to a linear evolution.

Figure 6 :

 6 Figure 6: Granular flow over an inclined rough plane for α = 0.406 rad and H = 14. Comparison of computed (h = H/40) and analytical solution profiles in the layer: u-velocity component (left), vertical velocity gradient du/dz (center) and shear stress τ (right).

Figure 7 :

 7 Figure7: Granular flow over an inclined rough plane for α = 0.406 rad and H = 14. Spatial convergence analysis of the numerical solution in terms of RMS error and infinity norm. The RMS error is defined as e rms = 1/N N (u iu ex i )2 and the infinity norm is defined as e ∞ = max N u iu ex i .

Figure 8 :

 8 Figure 8: Granular flow over an inclined rough plane for H = 14. Comparison between numerical (h = H/80, ) and analytical (solid line) free surface velocity u(z = H) as a function of the tilt angle α.

�Figure 9 :

 9 Figure9: Sketch of the heap-flow experimental configuration[START_REF] Jop | Crucial role of side walls for granular surface flows: consequences for the rheology[END_REF]

  the present model solutions. The proposed model reproduces quantitatively the experiments for both quantities even if the thickness of the flowing layer is slightly overestimated in the present numerical results. Nevertheless, one can observe differences in the velocity profiles close to the wall between computations from Jop et al. (2006) and those from the present ones. They could be due to: i) different spatial resolutions (70 grid points for Jop et al. and 281 nodes in our computations), ii) different spatial discretisation schemes of the governing equations (Finite Differenceversus Finite Element methods) and iii) different numerical approximation of the particulate hydrostatic pressure that enters the rheology. As stated by the authors, these results show that the proposed rheology gives quantitative agreement for this complex three-dimensional flow. Moreover, from a computational point of view, this comparison allows to validate the choice of a regularisation technique to deal with complex granular flows involving moving-static transition. As an illustration, figure11shows the longitudinal velocity profile in a cross section of the channel for the three dimensionless flow rate Q * ∈ {14.85 ; 39.71 ; 88.72}.

Figure 10 :

 10 Figure 10: Granular flow on a heap. a) Surface velocity profiles; b) Depths of the flowing layer across the channel width. Present numerical results

Figure 11 :

 11 Figure 11: Velocity profiles obtained with the present model (Grid size: 1x70x100, half domain, penalisation parameter: λ p = 10 -7 and regularisation parameter: λ r = 10 -5 ) for three different dimensionless flow rate (Q * ∈ {14.85 ; 39.71 ; 88.72}).

  Plotting the dimensionless drag force (F D = Cyl -→ T • -→ e z ds) versus the Froude number in figure 12 (left panel) (Fr = U/ √ gD , where U = Q v √ gD/(D 2 L x L y ) and Q v is the computed volumetric granular flow rate) reveals that for the highest Froude numbers considered the dynamical contribution to the drag force reaches the same magnitude as the static one. A slope break is observed on the curve at Froude numbers around 0.005 that corresponds to the limit between the artificial creeping flow regime induced by the regularisation and the actual flowing one as revealed in the close-up view. Finally, the drag coefficient, defined as C D = F D /( 1 /2 ρφU 2 ), is plotted versus the Froude number based on the effective cylinder diameter in figure 12 (right panel) for values higher than 0.005 that corresponds to the validity domain of the present numerical solution. In this regime, the numerical results are in quite good agreement with the experimental ones reported in Chehata et al. [9]. The best results are obtained with the highest external pressure (p ext = 10 ρφgD). Even if the numerical results are not strictly quantitative the dependency of the drag coefficient to the Froude number is quite well reproduced. As shown by Chehata et al. [9], in such slow granular flow, the drag force is almost independent of the flow velocity corresponding to drag coefficient proportional to U -2 . This behavior is also confirmed on the drag force versus Froude number evolution where the drag force tends to an asymptotic value as the Froude number increases. The nice collapse of the drag coefficient divided by the external pressure versus the Froude number observed in the close-up view of the right panel shows the linear dependency of the drag force to the external pressure. A better agreement with the experimental data could certainly be obtained by adjusting the empirical parameters of the rheology and/or modifying the geometry and theboundary conditions (such as the no-slip boundary conditions at the side walls). However, it is not our objective here to reproduce exactly this configuration. Our point is more to demonstrate the reliability of the proposed numerical model to compute dense granular flows in such complex configurations.

Figure 12 :

 12 Figure 12: Granular chute flow around a cylinder for different confining pressures: drag force divided by external pressure (F D /p ext ) versus Froude number U/ √ gD (left panel) and drag coefficient (C D = F D /( 1 /2 ρφU 2 ) versus Froude number based on the effective cylinder diameter U/ g(D + d) (right panel). In the close-up view the drag coefficient is divided by the external pressure whereas the horizontal axis is the modified Froude number. The symbols represent experimental data reported in Chehata et al. [9]

Figure 13 :

 13 Figure 13: Granular chute flow around a cylinder, iso-values of the primary variables for a dimensionless external loading of F v adim = 0.054.

Figure 14 :

 14 Figure 14: Granular chute flow around a cylinder, close-up in the vicinity of the cylinder of the iso-values of the primary variables for a dimensionless external loading of F v adim = 0.054.

Figure 15 :

 15 Figure 15: Granular chute flow around a cylinder for F v adim = 0.054. Plots along the cylinder perimeter versus the azimuthal angle of the dynamical pressure, total pressure, regularized one (plotted for -π ≤ θ ≤ 0), normal and tangential stresses and their regularised counterparts (plotted for 0 ≤ θ ≤ π).

  (a) Static component µ s /( γ 2 +λ 2 r ) 1/2 . (b) Dynamic component ∆µ/(I 0 √ φp + γ ). (c) Particulate viscosity (12).

Figure 16 :

 16 Figure 16: Granular chute flow around a cylinder, regularised particulate viscosity fields for or two dimensionless external loadings of F v adim = 0.044 (half left) and 0.054 (half right).

  direction between the two plates (y ∈ [-b/2; b/2]). z denote the vertical upward direction as shown in figure 1. The dimensionless momentum equations reduce to: 0 = ∂τ yz ∂y -1 ; 0 = -∂p ∂y

  this case: w(y) = 0 for y ∈ [-b/2 ; b/2]. On the contrary, if µ s p < b/2, for y ≥ µ s p = y 0 the material is sheared whereas it is not for values of y ≤ y 0 . The velocity profile exhibits a plug in the central region of the flow. Therefore, for p < b/(2µ s ) and assuming a µ(I) frictionnal rheology for the shear stress τ yz = µ(I) p the zcomponent of the momentum equation (20) reads: µ(I)p = y, (

) 25 )

 25 Introducing the expression for I = γ / φp in the previous equation gives an equation for γ = |dw/dy|: dw dy = I 0 φp yy 0 ∆µ p -(yy 0 ) . (For y ∈ [y 0 ; b/2], |dw/dy| and w are negative therefore integrating over y, gives: w(y) = -I 0 φp yy 0 + ∆µ p ln (∆µ p -(yy 0 )) + C * ∈ [y 0 ; b/2] ;

Table 1 :

 1 Granular flow over an inclined rough plane for α = 0.406 rad and H = 14. Newton-Raphson iteration counts and CPU time obtained with λ r = 10 -8

	h	H/5	H/10	H/20 H/40 H/80 H/160
	Degrees of freedom	297	567	1107 2187 4347 8667
	Newton-Raphson iteration counts 10000 10000	13	9	8	8
	CPU time (s)	109.56 149.16 2.68 2.28 2.63	3.6
	(s/iteration)	0.011	0.015	0.21 0.25 0.33	0.45

Acknowledgements

The authors thank Y. Forterre and O. Pouliquen for very fruitful discussions regarding the granular rheology and also acknowledge Y. Jobic for technical support. Funding from Agence Nationale de la Recherche (Project Dunes ANR-07-3-18-3892) is gratefully acknowledged and finally this work was also granted access to the HPC resources of IDRIS under the allocation 2012-96212 made by GENCI (Grand Equipement National de Calcul Intensif).

A. Analytical solution of the dry granular chute flow

The analytical solution of the dry granular vertical-chute flow based on the µ(I) rheology [START_REF] Gdrmidi | On dense granular flows[END_REF][START_REF] Jop | A constitutive law for dense granular flows[END_REF][START_REF] Forterre | Flows of dense granular media[END_REF] for the frictional shear stress is presented. Considering the steady flow of a granular material between two vertical parallel plates only the vertical downward velocity component w is not vanishing that only depends on y, the transverse and with no-slip boundary condition at the wall (y = b/2), the following velocity profile is found:

B. Analytical Bagnold profile

Considering the steady flow of a granular material on an infinite inclined plane only the longitudinal velocity component u is not vanishing and depends only on z the vertical upward direction (z = 0 at the bottom and z = H at the free surface). The mechanical equilibrium in this case is very simple, the shear stress balances the longitudinal component of the gravity and the pressure balances its vertical component. The momentum equations in dimensionless form reduce to:

The z-component of the momentum equation gives an hydostatic distribution of the pressure: p = Hz cosα (with p = 0 at z = H). Assuming a µ(I) frictionnal rheology for the shear stress τ = µ(I) p , the x-component of the momentum equation (28) simplifies as:

and at last µ(I) = tanα (constant with z). Using the definition of the inertial number I = du/dz / p , one can write the analytical expression for the vertical velocity gradient:

and by integrating on z with u(z = 0) = 0 as boundary condition, one obtains the Bagnold velocity profile:

Introducing the expression µ(I) = µ s + ∆µ/1 + I 0 /I, one finds the following expression for I α :


N El Gharbi

A Benzaoui

R Bennacer

Aerodynamics simulation of operating rooms

Keywords: Operating room, aerodynamics simulation, turbulent model, comfort, Airflow, Indoor air quality

The hospital is the place where we find simultaneously people whose health state is weakened, or vulnerable, and pathogenic micro-organisms able to worsen their health. Because these, the quality of the air in hospital must be in conformity with precise criteria in the buildings such everyday usage, and in particular in the buildings where some risk of specific pollution exist as in operating rooms. In this paper, we present a modelisation and three dimensional numerical study made in an operating room.

Results could be used when new operating rooms must be conceived or others existing will be modified. It could be used to avoid risks and to allow all controls of any risks. The air flow modelling aim, by analysing the stream coherence, is to control the contamination in the operating area by bringing out the contamination dispatching within the operating room. It will allow to clearly understanding the complex coupled phenomena thanks to animation and to virtual 3D reality. The use of real data of the studied operating room (geometry, volume, extracted and blown air flow, temperature and hygrometry) will allow the determination of the drainage of air, the distribution of temperature, and the zones of poor air distribution. To be able to find the exact distribution zones of the contaminants, the chosen turbulent model has to be accurate with a good ability to predict the recirculation of the air. The obtained results show that only one of the four tested models can correctly define the recirculation of the air.

Introduction

The air conditioning system of a hospital operating room must provide a comfortable and healthy environment for the patient and the surgical team. Thermal comfort can be achieved by controlling the temperature, the humidity, and the air flow. A healthy environment can be achieved by minimizing the risk of contamination through appropriate filtration and air distribution scheme. To ensure these optimal conditions, a study of the aerodynamics flow in a conditioned operating room must be made, by using the digital simulation. [START_REF] Chen | Comparison of different k-ε models for indoor air flow computations[END_REF] Choice of the turbulent model According to literature, Nielsen P. V (1974), was one of the first researchers having used the k-ε turbulent model in 2D, to study the air movement and heat transfer in conditioned room [START_REF] Nielsen | V Flow in air conditioned rooms Ph[END_REF]. Murakami et al (1994). From the University of Tokyo quoted by [START_REF] Weathers | A study of computational fluid dynamics applied to room air flow Master of Science[END_REF] added to the standard kε model a wall function to model the air flow, in a clean room where the number Air Change by Hour (ACH) is very high. [START_REF] Chow | Y Performance of ventilation system in a non-standard operating room[END_REF] used the standard k-ε model using eddy viscosity hypothesis fort their numerical predictions on the ventilation of an operating theatre in a hospital of Honk Kong, [START_REF] Chow | Y Performance of ventilation system in a non-standard operating room[END_REF]. Similarly the work made by Monika et al(2004), who used the same turbulent model and the flow near the boundaries were represented by using the standard logarithmic law [START_REF]Diagonal air-distribution system for operating rooms: experiment and modeling[END_REF]. They focused on the contaminant diffusion in an experimental operating room. [START_REF] Chen | Comparison of different k-ε models for indoor air flow computations[END_REF] [START_REF] Chen | Comparison of different k-ε models for indoor air flow computations[END_REF] studied five modified k-ε models and compares the obtained numerical results with existing experimental results. The RNG k-ε model was found slightly better than the standard kε model for simulating air flow displacement.

While [START_REF] Luo | Numerical study of three dimensional turbulent flows in a habitat with coupled heat and mass transfer Ph[END_REF] [START_REF] Luo | Numerical study of three dimensional turbulent flows in a habitat with coupled heat and mass transfer Ph[END_REF] compared three turbulent models, RNG k-ε , realizable k-ε and SST k-w , where an enhanced treatment near the walls was used for the two k-ε models. The study of the 3D ventilation shows that the SST k-w model predicts better the velocity profiles in the vicinity of floor. For our study the choice of the turbulent model is crucial to well predict the air recirculation. Such recirculation will be favourable for the accumulation of the contaminants. Four turbulent models (two equations) were studied:

• The standard k-ε model with walls function • The realizable k-ε model with non equilibrium walls function ;

• RNG k-ε model with non equilibrium walls function;

• SST k-w model.

The choice of the treatment associated to walls as each models was concludes from the results of the literature search [START_REF] Luo | Numerical study of three dimensional turbulent flows in a habitat with coupled heat and mass transfer Ph[END_REF].

Adopted methodologies

The first phase is the determination of our operating room characteristics: 2.1 Specific Area with very high risks for the patients (class 4 according to NF S 90-351 or ISO 5 according to ISO'S DIN 14644-1).

Aerodynamic

• The distribution of air is diagonal • The renewal of air is 40% of new air and 60% of recycled air, with 6.08ACH

• The air flow ventilation is equal to 1852m3 / h, with a velocity of 0.51m/s.

• The recycled flow is equal to1110m3 / h, with a velocity of 0.61m/s.

• The stale flow is equal to 495m3 / h, with a velocity of 0.27m/s.

• The blowing temperature is 17.5°C with 40% humidity.

Geometrical

The following parameters were measured:

• Volume of the room.

• Dimensions of the walls and of the doors.

• Dimensions and positions of the various objects in the room (Scialytic-operating table -medical equipment -staffs). • Dimensions and positions of the inlet and outlet diffuser.

The second phase consists in calculating the heat flows released by the occupants, lighting, by the machines and transmission, to be able to use them as boundary conditions.

Presentation of the operating room

The studied operating room is assigned to the service of the neurosurgery. The arriving air in the room is already filtered in the station of air treatment, using primary filters of 88% and 95% effectiveness, and then in a second time on the level of the inlet diffuser using HEPA filters. For our study, we consider the room empty without staff neither patient nor medical equipment, to define the air recirculation without presence of disturbance which can increased these recirculation. The following heat flows were taken into account: 1 The heat dissipation

The fluid is supposed incompressible and a second order scheme is used to discreditez the diffusion terms. The solution is obtained in primitive variable P-V and the pressure and velocity coupling is obtained by the SIMPLEC algorithm. The convergence criterion was based on the maximum error less than a prescribed value taken equal to 10 -6 for energy, and 10 -4 for the other equations plans were selected Each plan is localised in the middle of the inlet diffuser represented on figure 2 and localised at:

• z =1.13m in the medium of the first inlet diffuser.

• z = 2.54m in the medium of the second inlet diffuser.

• z = 3.95m in the medium of the third inlet diffuser. z = 4.47m in the medium of the fourth inlet diffuser. We notice about the time duration of the iterations made with these four turbulent models: standard K-ε, realizable k-ε, RNG k-ε and SST k-w, the first two models converge more quickly compared to the last two models. After the convergence of the four models, we obtain the following results: With the four turbulent models there are some recirculation areas of the air with a return of the velocity vectors in the opposite direction of the flow when it runs up against the wall. These recirculation's are visualized on the level of the floor with SST k-w model, still better with the standard model Kε , on the other hand RNG k-ε model reflects very well these recirculation's in the interior of the room.

X[m] Y[m]
• At the plan z = 2.54m With the realizable k-ε model and SST k-w model, one does not visualize the zones of recirculation, this recirculation are weak with the two other models.

X[m] Y(m]
• At the plan z = 3.54m Figure 5 The recirculation at z = 3.95m It is the plan which passes by the operating table, the realizable k-ε model does not envisage recirculation of air, with SST k-w model we observe a weak recirculation meadows of the diffuser inlet, on the other hand with the standard k-ε these recirculation's are very well visualized under the table, similar for the model RNG k-ε which, in addition to the recirculation under the table, appear other ones in the room. • At the plan z = 5.47m Figure 6 The recirculation at z = 5.47m

From this point of view one can also see the two outlet extraction of the air, always with the realizable K-ε model it does not have there recirculation of the air, one can visualize with SST k -w model, but visualized well with the two last models.

Conclusion

Aerodynamics modelling makes it possible to optimize the way of air flow, in order to decrease for the patients, any contamination risks inside the protected zone. The visualization of the zones of recirculation of air in the studied operating room, with the four models turbulent, makes it possible to conclude that model RNG k-ε is the best to be adapted to define these zones of recirculation in all the room. The results obtained during this study, tend to confirm the practical interest of modelling. This solution makes it possible to bring concrete answers for the control of contamination, by highlighting the zones of bad distribution of air potentially contaminated.

Figure 1

 1 Figure 1 Geometry of the operating room

Figure 2

 2 Figure 2 Plans of study

3•Figure 3

 3 Figure 3 The recirculation at z = 1.13m

3Figure 4

 4 Figure 4 The recirculation at z = 2.54m

3

  Realizable k-e at z = 5.47m