L Bennani
email: lokman.bennani@imft.fr

H Neau

C Baudry
email: cyril.baudry@edf.fr

J Laviéville

P Fede

O Simonin

Numerical simulation of unsteady dense granular flows with rotating geometries

Keywords:

In chemical engineering applications, it is not uncommon to encounter reactors featuring rotating parts. As these rotating parts are present in order to enhance processes such as chemical reactions and/or ensure homogeneity, it is essential to take them into account to perform predictive numerical simulations. This aspect can be particularly challenging, even more so when complex industrial geometries are to be treated. In this paper a numerical methodology for simulating unsteady granular flow in rotating geometries is presented. The method is based on splitting the domain into static and rotating parts. The information between rotating and static parts is passed by a non-conformal mesh matching technique. The presented methodology is validated numerically by comparing its results with other conventional methods. The method is then applied to an industrial scale problem. The applicability of the method and the way it may be used to investigate complex flow is demonstrated.

Therefore this approach enables to consider the full geometry of complex reactors. It opens the door to further investigation, optimization and design of industrial scale chemical processes.

Introduction

Flow of granular materials in rotating geometries is encountered in a wide range of chemical engineering applications, such as the mixing process in axially stirred chemical reactors or grain drying in rotary kilns [START_REF] Boateng | Rotary Kilns: Transport Phenomena and Transport Processes[END_REF][START_REF] Soares | Polyolefin Reaction Engineering[END_REF]. These rotating parts are often present so as to enhance the chemical process. The gain in performance can be investigated experimentally. However, experimental investigation of large scale chemical reactors involving such a motion is often complicated. Hence there is a need for efficient numerical simulation strategies which are able to capture the flow of granular materials in rotating geometries.

A possible numerical approach is the Discrete Element Method (DEM). In effect, this method has been used to inves- Another method, which is the one adopted in this work, is to solve the flow equations in an Eulerian framework. Although it is more difficult to account for complex phenomena with this approach, it does not require the tracking of every particle. Hence more suited for the simulation of large scale reactors [START_REF] Zeren | Numerical study of solid particle axial mixing in a fixed cylindrical drum with rotating paddles[END_REF]. Nevertheless, taking rotating parts into account is a complex task. Indeed, if the geometry of interest features rotational symmetries a possible solution is to solve the governing equations in a rotating frame of reference or by using sliding wall boundary conditions (well adapted to rotating drums for instance). For example, [START_REF] Huang | Numerical studies of particle segregation in a rotating drum based on Eulerian continuum approach[END_REF] and [START_REF] Santos | Experimental and CFD study of the hydrodynamic behavior in a rotating drum[END_REF] have investigated the granular flow in a rotating drum using an Eulerian framework and by imposing a no-slip boundary condition at the rotating wall. In addition, [START_REF] Ahmadzadeh | Population balance equations' application in rotating fluidized bed polymerization reactor[END_REF] studied a rotating fluidized bed reactor using a two dimensional rotating frame approach. However, in most problems such simplifications are not possible. This can be due, for example, to the presence of chimneys or extraction pipes. Hence, the simulation of industrial reactors containing rotating parts remains a great challenge and a more elaborate numerical strategy is needed.

In this paper, a rotating mesh method for simulating unsteady granular flows in complex large scale geometries is presented. It is a first step towards solving the aformentionned problems. In addition, the unsteadiness is an aspect of prime importance for practical applications, in particular when chemical reactions occur. The main idea is to split the domain into static and rotating parts and to use a nonconformal mesh matching technique to connect the domains [START_REF] Edf | Code Saturne 4.0 Theory Guide[END_REF].

The main goal of this work is to validate the presented methodology by comparing it to other well known methods. Another approach could have been to use the method of manufactured solutions [START_REF] Blais | On the use of the method of manufactured solutions for the verification of CFD codes for the volume-averaged Navier-Stokes equations[END_REF][START_REF] Oberkampf | Verification and Validation in Scientific Computing[END_REF]. The reader should bear in mind that this work is of a numerical nature. The specific models or boundary conditions can be changed/implemented without any difficulty related to the rotating mesh method.

The first part of the paper is focused on presenting the physical modelling approach and the rotating mesh method.

The numerical validation of the method is then addressed by considering a simple rotating drum test case. Indeed, as this geometry presents rotational symmetry, a simulation using a sliding wall boundary condition is also possible. A first assessment is performed by comparing the results obtained with both methods.

Thirdly, a flat bladed stirrer configuration investigated by [START_REF] Stewart | Granular flow over a flat-bladed stirrer[END_REF] is simulated. This configuration is more complex and features interesting flow physics. The case also presents a rotational symmetry making it possible to simulate in a rotating frame of reference (by adding the inertial forces). Comparison of the obtained results allows to further validate the method.

Finally, the method is pushed into the realm of industrial problems and is applied to a prototype of a horizontally stirred reactor. This proves that the presented approach enables to consider the full geometry of complex reactors. Hence, it provides a step further in the investigation and design of industrial scale chemical processes.

Physical model

The physical modelling is based on a particulate Eulerian approach. It is derived from a joint fluid-particle PDF equation allowing to obtain the transport equations for the mass, momentum and fluctuating kinetic energy of particle phases in rapid granular flows [START_REF] Balzer | Eulerian gas-solid flow modelling of dense fluidized bed[END_REF][START_REF] Simonin | Combustion and turbulence in two-phase flows[END_REF][START_REF] Balzer | Gas solid flow modelling based on the kinetic theory of granular media: validation, applications and limitations[END_REF]. These equations are standard in two-fluid modelling. Briefly, the set of equations consists of mass, momentum and fluctuating kinetic energy equations coupled by interfacial transfer terms. More details on this approach are given in Appendix A. Moreover, the numerical methodology presented here is general enough to be applied to other approaches.

The reader should bear in mind that the goal of this work is to assess the numerical performance of the rotating mesh method. That being said, given the high degree of compaction of the particle phase and the important presence of walls in the cases which are presented in this paper, the particle phase stress tensor and wall boundary conditions are briefly discussed.

Particle phase stress tensor

The flows under consideration in this study present a high degree of compaction. This means that frictional interaction between solid particles is an important aspect and has to be taken into account. This is done through the introduction of a frictional part in the stress tensor in the momentum equation.

The particle phase stress tensor p is modelled as the sum of a kinetic-collisional part derived in the frame of the kinetic theory of rapid granular flow modified by the drag force [START_REF] Boëlle | Second-order prediction of the particle-phase stress tensor of inelastic spheres in simple shear dense suspensions[END_REF], g p and a frictional part, f p . The collisional and kinetic part of the particle phase stress tensor is written using a viscosity assumption:

g p,ij = P p -p ∂U p,l ∂x l ı ij -p S p,ij , (1)
where P p is the granular pressure, p the bulk granular viscosity and p the kinetic-collisional viscosity (see Appendix A). Moreover, S p is defined in the following way:

S p,ij = ∂U p,i ∂x j + ∂U p,j ∂x i - 2 3 ∂U p,l ∂x l ı ij . (2)
Concerning the frictional part of the stress tensor, the approach proposed by [START_REF] Srivastava | Analysis of a frictional-kinetic model for gas-particle flow[END_REF] is retained. It is based on rheological laws used in soil mechanics which assume that, in the high compaction regime, the granular material follows a rigid-plastic behaviour. In this context, Srivastava and Sundaresan propose the following model for the frictional stress: where d p is the particle diameter, is the angle of internal friction, q 2 p is the particle kinetic agitation and 8q 2 p /3d 2 p is an estimate for strain rate fluctuations. p c is the critical state pressure defined by:

f p = p c I - √ 2 sin( )p c S p S p : S p + (8q 2 p /3d 2 p ) , (3)
p c (˛p) = ⎧ ⎪ ⎨ ⎪ ⎩ F ˛p -˛p,min r (˛p ,max -˛p) s if ˛p > ˛p,min 0 if ˛p ≤ ˛p,min (4)
where F, r, s are constants. The different parameters used here to define the frictional model are summed up in Table 1.

It should be noted that the values were taken as given by [START_REF] Srivastava | Analysis of a frictional-kinetic model for gas-particle flow[END_REF] except for ˛p,max . The modelling of frictional stresses in an Eulerian framework is still an active area of research [START_REF] Schneiderbauer | A comprehensive frictional-kinetic model for gas-particle flows: analysis of fluidized and moving bed regimes[END_REF][START_REF] Chialvo | Bridging the rheology of granular flows in three regimes[END_REF]. However, note that the model retained for this work suffices to reflect the basic physics of frictional regimes. The implementation of more sophisticated frictional models is part of current work and should present no difficulty associated to the rotating mesh method.

Particle phase wall boundary conditions

Throughout this work, the particle phase wall boundary condition that is used is a no-slip condition. Therefore, a zero flux is imposed for scalar variables, in particular for the solid random kinetic energy [START_REF] Fede | 3D numerical simulation of a lab-scale pressurized dense fluidized bed focussing on the effect of the particle-particle restitution coefficient and particle-wall boundary conditions[END_REF]. As for velocities, a no-slip condition is applied. Note this type of boundary condition corresponds to the maximum of friction at the wall. Other wall boundary conditions exist and have been investigated [START_REF] Soleimani | A comparison for different wall-boundary conditions for kinetic theory based two-fluid models[END_REF][START_REF] Fede | 3D numerical simulation of a lab-scale pressurized dense fluidized bed focussing on the effect of the particle-particle restitution coefficient and particle-wall boundary conditions[END_REF]. However, their implementation should present no difficulty associated to the rotating mesh method.

The rotating mesh numerical method

The numerical simulations presented in this paper have been carried out using a Eulerian n-fluid modeling approach for gas-solid turbulent polydisperse flows developed and implemented by IMFT (Institut de Mécanique des Fluides de Toulouse) in NEPTUNE CFD version. NEPTUNE CFD is a multiphase flow software developed in the framework of the NEPTUNE project, financially supported by CEA (Commissariat l'Energie Atomique), EDF (Electricité de France), IRSN (Institut de Radioprotection et de Sûreté Nucléaire) and AREVA-NP. The numerical solver has been developed for High Performance Computing [START_REF] Neau | Neptune CFD high parallel computing performances for particle-laden reactive flows[END_REF][START_REF] Neau | High performance computing (HPC) for the fluidization of particle-laden reactive flows[END_REF].

Methodology

The rotating mesh method considered here operates by taking as inputs several meshes, some of them will be fixed (stator), others will be rotating (rotor) [START_REF] Audebert | Code saturne, mise en place et validation de la fonctionalité couplage rotor/stator pour la modélisation des pompes[END_REF]. The equations of motion are solved in each mesh in the absolute frame of reference.

For the rotating mesh, the governing equations must be formulated so as to take into account its motion. To do so, an ALE-like (Arbitrary Lagrangian-Eulerian) approach is used. As the ALE method is well established [START_REF] Hirt | An arbitrary Lagrangian-Eulerian computing method for all flow speeds[END_REF][START_REF] Hughes | Lagrangian-Eulerian finite element formulation for incompressible viscous flows[END_REF][START_REF] Donea | An arbitrary lagrangian-Eulerian finite element method for transient dynamic fluid-structure interactions[END_REF][START_REF] Donea | Arbitrary Lagrangian-Eulerian Methods[END_REF], only the basic general principles of the approach are presented. As shown in Fig. 1, let x be the Eulerian coordinates and be the so called reference coordinates. Moreover, ( , t) denotes the transformation from to x. Consider a quantity of interest, which can be viewed equivalently in both configurations through the relation: g(x, t) = g( ( , t), t) = g( , t). Hence, g( , t) represents the function as attached to the moving mesh. In this framework, one may relate the time derivative of g to that of g through the formula:

∂Jg ∂t = J ∂g ∂t + ∇ • (gW) (5)
where W = ∂ /∂t is the mesh velocity and J is the Jacobian of the transformation . In the present case, the mesh movement is a rotation and thus J = 1. The governing equations, in a rotating mesh, for a phase k hence take the following form:

∂ ∂t ( ˜k ˜ k ) + ∇ • (˛k k (U k -W)) = 0 (6) ∂ ∂t ˜k ˜ k Ũk + ∇ • (˛k k (U k -W) ⊗ U k ) = -˛k∇P g + ˛k k g +I k -∇ • (˙k) (7)
Eqs. ( 6) and ( 7) mean that in the finite volume procedure the mesh rotation will essentially impact the mass flux and convection terms. The general method for computing these terms is briefly recalled below. Fig. 2 shows the geometrical definitions used for the numerical scheme. L and R denote the left and right cell respectively, n LR is the unit normal to the face adjacent to L and R and S LR is its area. Below, XY denotes the vector joining a given point X to a given point Y.

The mass flux at a given interior face is computed using the following scheme:

[˛k k U k ] F .S LR n LR = ˛k,F (p k [ k U k ] A + (1 -p k )[ k U k ] B ) • S LR n LR (8) where [ k U k ] A is defined as: [ k U k ] A = k,A U k,A + 1 2 min(˛k ,A , ˛k,B ) ([∇U k ] A • AA + [∇U k ] B • BB )) (9)
and the ponderation coefficient p is defined as:

p k = c˛k ,A c˛k ,A + (1 -c)˛k ,B (10)
where c = |FB |/|A B | and ˛k,A is defined as:

˛k,A = ˛k,A + [∇˛] A • AA (11)
Note that if the cell belongs to a rotating part the velocity used to compute the mass flux is U k -W in accordance with Eq. ( 6).

As for the convective flux between cells L and R, noted F LR (˛k k U k , U k ), it is computed as follows:

F LR (˛k k U k , U k ) = ([˛k k U k ] F • S LR n LR ) • U k,F (12)
where U k,F is defined depending on the chosen numerical scheme. The one being used in this work is the Second-Order Linear Upwind Scheme (SOLU) which is given by:

U k,F = { U k,A + (∇U k ) A • AF if [˛k k U k ] F • S LR n LR ≥ 0 U k,B + (∇U k ) J • BF if [˛k k U k ] F • S LR n LR < 0 (13)
In Eq. ( 12) the computation of the convective flux is based on the mass flux. Therefore, if the face belongs to a rotating cell, no further modification is needed as the rotation velocity will already have been accounted for in the mass flux.

In addition to taking the mesh rotation velocity into account when computing the mass and convective fluxes, the actual rotation of the mesh has to be integrated into the code's time marching scheme. The time marching scheme is based on an elliptic oriented fractional step method [START_REF] Méchitoua | An unstructured finite volume solver for two-phase water-vapour flows based on an elliptic oriented fractional step method[END_REF]. The method can be seen as a prediction/correction scheme. The actual rotation of the mesh is performed between the prediction and the correction steps. Therefore, the algorithm can be summed up as follows:

1. Velocity prediction step:

U n k -→ U pr k
2. Rotate rotor mesh by an angle ω t (ω the angular velocity and t the current time step). Then, match to stator mesh using a non-conformal matching algorithm 3. Correction step: compute all variables at iteration n + 1 through the pressure-mass-momentum-energy coupling step As stated in step 2 a non-conformal mesh matching algorithm is used after rotation. Indeed, even if the matching between meshes is conformal at the first iteration, it will not be the case at the subsequent ones. Indeed, the rotor mesh will have rotated by an angle which does not necessarily coincide with the angular discretization of the meshes. Hence, in order to solve the governing equations, information has to be passed at the non-conforming interface between the stator and the rotor. This is done thanks to a mesh joining algorithm and has to be performed at each time step. The main idea of this method is to intersect the faces of both interfaces and split them so as to build conforming meshes. An illustration of this method is shown in Fig. 3. More information on this procedure can be found in EDF (2015).

It should also be noted that the rotation imposes a restriction on the time step. Indeed, so as to remain consistent with respect to the flow physics, the curvilinear distance traveled by the rotating mesh at the joining interface during a given time step should not exceed that of a mesh cell. In other words: ωR j t ≤ s j , where R j is the joining radius between rotor and stator and s j is the cell size in the tangential direction.

The rotating drum (without gravity)

The goal of this section is to numerically validate the rotating mesh method with a simple test case. Indeed, consider a rotating drum as illustrated in Fig. 4. The drum is rotated at ω = 1.21 rad s -1 and initially uniformly filled with gas and particles with volume fractions respectively ˛g = 0.56 and ˛p = 0.44. Moreover, as this case is purely numerical, gravity is not present. In this case, the gas flow is considered laminar while transport equations are solved for the particle agitation kinetic energy q 2 p (see Appendix A for further details). The problem is treated in two dimensions. Hence, the lateral faces are treated as symmetry boundary conditions. This implies a homogeneous Neumann condition for scalar variables. For the velocities, the symmetry boundary condition imposes a homogeneous Dirichlet condition on the normal component and a homogeneous Neumann on the tangential one. The outer cylinder boundary is treated as a wall, which means a noslip boundary condition for velocities and zero flux for scalar variables, in particular for the solid random kinetic energy [START_REF] Fede | 3D numerical simulation of a lab-scale pressurized dense fluidized bed focussing on the effect of the particle-particle restitution coefficient and particle-wall boundary conditions[END_REF].

As the geometry features rotational symmetry, the case can also be simulated by using a sliding wall boundary condition. This enables a purely numerical validation by comparing the results obtained with both numerical approaches. The sliding wall approach consists in imposing a tangential velocity at the cylinder's boundary. In this way, the rotating motion of the cylinder is effectively taken into account. On the other hand, as explained in Section 3, the rotating mesh method consists in considering two meshes. An inner fixed mesh and an outer rotating crown mesh. Numerical validation is then assessed by comparing the results obtained with both methods. Fig. 5 illustrates the differences between the two approaches. The mesh (Fig. 5) is composed of 1280 cells in the crown and 3840 cells in the fixed area, making a total of 5120 cells.

The density of the gas is set to g = 1.2 kg m -3 and its viscosity to g = 1.85 10 -5 Pa s. Table 2 gives the physical parameters of particle phase (glass beads).

As the flow exhibits a rotational symmetry, the results are extracted along the radial direction. In Fig. 6 the profiles of solid volume fraction obtained at different times (t = 10 s, 20 s, 30 s and 40 s) with both approaches are plotted and compared. A very good agreement between both methods is achieved. Fig. 7 shows a comparison between tangential velocities for both approaches and at various times. A very good agreement is observed. Starting out from the center of the drum, the velocity profile first presents a low but almost linear slope before transitionning to higher values and a higher slope. Steady state is reached for this case after approximately 90 s. As shown in Fig. 8, in this state, all the particles are packed near the wall and behave as a rigid body.

The above results show that the rotating mesh method is fully able of correctly reproducing the numerical results obtained with a conventional sliding wall approach. The simulations were performed with an adaptive time step. For both cases the time steps at every iteration were very close. At steady state, the time step was 0.018 s and limited by a maximum CFL of one. However, the rotating mesh method is more computationally expensive due to the non-conformal mesh matching which is required at each time iteration. Table 3 provides an idea of this additional cost.

Also, it is interesting to note that changing locations for the mesh matching interface does not change the results. Indeed different locations were tested but not reported here for conciseness.

The flat bladed stirrer

With the rotating mesh method now validated on a simple case, the road is open to the investigation of more complicated cases. In order to further assess the method, a very interesting case is that of the flat bladed stirrer, investigated experimentally by [START_REF] Stewart | Granular flow over a flat-bladed stirrer[END_REF]. In this section, the flat bladed stirrer case is simulated with both the rotating mesh method and the rotating frame method. Indeed, the geometry of the flat bladed stirrer (described below) has a rotational symmetry with respect to its vertical axis. One may hence also simulate this case in a rotating frame of reference, which requires the addition of the inertial forces (Coriolis and centrifugal).

Geometry and mesh

The geometry of the stirrer consists of a vertical shaft to which are attached two opposed flat blades. This assembly is con- It should be noted that a gap of 2 mm was left between the blade tips and the vessel walls in order to be able to use the rotating mesh method. Fig. 10(a) shows the skin mesh on the shaft and the blades. The stator mesh consisted in 240 000 elements while the rotor was discretized with 125 000 elements.

For the present study the filling is of 2.8 kg of particles (approximately 200 000). The shaft is rotated at 20 rpm. The physical properties of the gas and particle phases are the same as in section 4 except for the particle diameter which is set to 2.2 mm (see Table 2 for particle physical properties and 1 for the frictional viscosity parameters). For this case, all boundaries are walls.

Results and discussion

In this section, the results obtained using the rotating mesh method are compared to those obtained with the rotating frame method. For this purpose, all results are presented in the rotating frame of reference. It should be noted that although some aspects linked to the dynamics of the flow are pointed out, the main focus is the comparison of both approaches so as to assess the rotating mesh method. It is however interesting to note that all the computed flow patterns are consistent with the observations of [START_REF] Stewart | Granular flow over a flat-bladed stirrer[END_REF].

Starting with the flow pattern in a radial cut, Fig. 11 shows a comparison of the particle bed heights. They are plotted as a function of the curvilinear abscissa (s) and along the full circumference at r = 0.08 m. The predicted bed heights compare very well. One may also note that a heap is captured, which is a primary feature of the flow.

Fig. 12 shows the tangential velocity along the z-axis at r = 0.08 m and two different angles ( = 37.5 • and  = 60 • ). The agreement between the predicted velocities is very satisfactory. Moreover, the change in sign shows that in the higher part of the heap the particles are pushed back. Coherently, a recirculation pattern in the heap can be observed in Fig. 13. It shows the tangential projection of the particle velocity field on a radial cut at r = 0.08 m.

Concerning the flow in vertical cuts, Fig. 14 shows the tangential velocity along the r-axis at two different angles and heights. The agreement between the predicted velocities is excellent. The profiles in Fig. 14 show that the particle motion is towards the blade near the walls and away from the blade in the near shaft area. This pattern becomes clearer in Fig. 15, which shows the tangential projection of the particle velocity field on a vertical cut at z = 0.026 m. The global flow of the particles in this vertical cut is towards to blades. However, in the near shaft region, a small recirculation develops.

As for the previous case, the time step was adaptive and remained very close for both cases at every iteration. At steady state, the time step was 0.001 s and limited by a maximum CFL of one. Table 4 provides an idea of the additional cost when using the rotating mesh method.

Towards industrial scale applications: prototype of a horizontally stirred reactor

In this final section, the applicability of the rotating mesh method to full scale industrial problems is demonstrated. Indeed, it is for these kinds of complex applications that the method proves to be the most usefull. In this section, the hydrodynamics of a gas-particles flow inside a numerical prototype of a horizontally stirred reactor is simulated. As shown in Fig. 16, the numerical prototype has a chimney. This means that the use of methods based on rotational symmetries would not have been possible.

Geometry and mesh

The geometry of the reactor is strongly inspired from the illustrations that can be found in [START_REF] Soares | Polyolefin Reaction Engineering[END_REF]. It is composed of a horizontal stirrer and a chimney, as shown in Fig. 16(a). As stated by [START_REF] Soares | Polyolefin Reaction Engineering[END_REF], the horizontal vessel is usually 10-15 m long and 1.5-4 m in diameter. Therefore, in the case of our numerical prototype, the vessel is 11 m long and 4 m in diameter. Also, the chimney is 6 m high. The rotating part is the axial stirrer contained within the horizontal vessel whereas the static (stator) part is the rest of the vessel and the chimney (see Fig. 16(b)). It should be noted that this is not a reproduction of an actual reactor. It is a numerical prototype designed only to demonstrate the applicability of the method to industrial scale problems.

The mesh was composed of 1.3 million hexahedral elements. The rotating part was meshed with 600 000 elements and the static part with 700 000 elements. The resulting mesh is shown in Fig. 17.

Results and discussion

In this case, the reactor is filled with 25 tonnes of particles (located in the horizontal vessel). Their density is p = 500 kg m -3 and they have a diameter of d p = 0.5 mm. This is equivalent to more than 10 10 particles in the reactor. The physical properties of the gas are chosen so as to be more representative of industrial processes. Values from the literature where used as an order of magnitude (see for example [START_REF] Kaneko | DEM simulation of fluidized beds for gas-phase olefin polymerization[END_REF]). For this case, the density is g = 45 kg m -3 and the viscosity g = 1.1 × 10 -5 Pa s. Concerning the boundary conditions, the top of the chimney is an outlet whereas all other boundaries are walls.

It should be noted that for this case, the flow of the gas phase can no longer be assumed laminar. Hence, a kmodel which includes additional terms accounting for the influence of the particles on the fluid turbulence is used [START_REF] Vermorel | Numerical study and modelling of turbulence modulation in a particle laden slab flow[END_REF].

In terms of computational cost, on 100 cores, 4.3 h are needed to simulate 20 s of physical time.

The hydrodynamics may first be investigated by looking at the solid volume fraction. Fig. 18 shows solid volume fractions at the skin of the reactor. A kind of wavy structure is observed which is due to the blades passing and projecting particles towards the boundaries. Moreover, the peaks of the waves are shifted coherently with the blades that are passing by. It is also of interest to look at the particle velocity fields. Let's first consider a cut in the YZ plane, as shown in Fig. 20(a). The cut is done just after a set of blades has passed through the plane. In the upper half, the velocity vectors show a complex pattern. Where the blade has just passed the particles are pushed upwards. On the other hand, the particles in the most upper part of the vessel are moving back down, creating the kind of circular pattern which is observed.

Secondly, let's consider a cut in the XY plane (Fig. 20(b)). The cut is done at the same time as the one in the YZ plane. However, note that the cut is performed so as to intersect with a blade that is almost at 0 • . Here also, complex patterns can be observed, with particles being pushed upwards, others falling back down and the bulk of the bed being stirred at the bottom.

Conclusions and perspectives

The goal of this paper was to present a rotating mesh method for simulating dense gas-particle flows in rotating geometries. This work focuses on the validation of the method by com- paring it with more conventional approaches (sliding wall, rotating frame). Within the framework of an Eulerian n-fluid modeling approach, this method enables to treat problems featuring complex geometries. After having presented the physical modelling and the numerical method, a rotating drum case was considered. For the purpose of this case, gravity was not present and the drum was initially filled with a uniform gas-particle mixture. The problem was also treated using a sliding wall approach. The numerical results obtained with both methods were compared and showed excellent agreement. From this step, it appeared that the rotating mesh method is equivalent to the conventional sliding wall approach and is fully able of capturing the main features of the flow.

Thirdly, the dense granular flow in a flat bladed stirrer was simulated using the rotating mesh method and a rotating frame approach. The results from both methods compared very well. In addition, the method captured recirculation regions, which is an important feature of this flow. This shows that the numerical method is able to treat flows with more complicated patterns.

Finally, the applicability of the method to industrial size complex problems was demonstrated. Indeed, the method was used to investigate the hydrodynamics of the gas-particle flow inside a prototype of a horizontally stirred reactor. This enabled to identify several flow patterns and hence opens the doors to interesting investigations of gas-particle flows in complex geometries. Hence, the method has strong potential for applications to industrial reactor design and the investigation of chemically reacting particulate flows in rotating geometries.

Overall, the rotating mesh method therefore presents a promising perspective for the Eulerian simulation of dense granular flows in rotating geometries. This will enable the assessment and validation of physical models by comparing simulations to experiments making use of a rotating apparatus (such as rotating drums). Moreover, the gain in using this method is that it allows to simulate complex problems of industrial size. For example, a configuration of interest is the simulation of a full scale chemical reactor featuring chimneys and other complex parts, and including the modelling of chemical reactions. These aspects are part of ongoing work. q gp -2k g + V d,i V r,i (A.28) kg p→g represents the effect of particles. Note that the production of agitation by the wake of a particle is assumed negligible compared to the large scale turbulent agitation [START_REF] Vermorel | Numerical study and modelling of turbulence modulation in a particle laden slab flow[END_REF].

The Boussinesq hypothesis is used to model the Reynolds stress tensor, yielding: u g,i u g,i = - In Eq. (A.34), the first term on the right-hand-side is the turbulent transport of the particle agitation. That term is written by introducing the diffusivity coefficients: The second term on the right-hand-side of Eq. (A.34) represents the production of particle agitation by the gradients of the mean solid velocity. The third term is the interaction with the gas. Finally the fourth term is the particle agitation dissipation by inelastic collisions.

The transport equation for the gas-particle covariance reads:

˛p

  Fig. 1 -Reference and transformed meshes and their associated coordinates.

Fig. 3 -Fig. 4 -

 34 Fig. 3 -Illustration of interface joining: the red faces and black faces are joined by intersecting their edges and building new faces (dotted lines). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.) Adapted from EDF (2015).

Fig. 5 -Fig. 6 -

 56 Fig. 5 -Illustration of the sliding wall and rotor/stator mesh approaches.

Fig. 7 -Fig. 8 -

 78 Fig. 7 -Comparison of the tangential velocities obtained at different physical times.

Fig. 9 -

 9 Fig. 9 -Stirrer geometry and dimensions.

Fig. 10 -

 10 Fig. 10 -Stirrer mesh and rotor/stator.

Fig. 13 -

 13 Fig. 13 -Tangential projection of the particle velocity field on a radial cut at r = 0.08 m (rotating mesh).

Fig. 14 -

 14 Fig. 14 -Tangential velocity along the r-axis at two different angles and heights.

Fig. 15 -

 15 Fig. 15 -Tangential projection of the particle velocity field on a vertical cut at z = 0.026 m (rotating mesh).

Fig. 16 -

 16 Fig. 16 -View of the geometry.

Fig. 18 -

 18 Fig. 18 -Solid volume fraction at the boundaries after two consecutive blade passes.

Fig. 18

 18 Fig. 18 shows solid volume fractions in an YZ cut. Here again, structures are formed by the blade motion which sends the particles upwards.It is also of interest to look at the particle velocity fields. Let's first consider a cut in the YZ plane, as shown in Fig.20(a). The cut is done just after a set of blades has passed through the plane. In the upper half, the velocity vectors show a complex pattern. Where the blade has just passed the particles are pushed upwards. On the other hand, the particles in the most upper part of the vessel are moving back down, creating the kind of circular pattern which is observed.

Fig. 19 -Fig. 20 -

 1920 Fig. 19 -YZ cut: solid volume fraction after two consecutive blade passes.

F

  gp (1 + ˛pg 0 c ) / 1 +

Table 1 -Frictional model parameters.

 1

	˛p,min	˛p,max	F[Pa]	r	s
	0.5	0.64	0.05	2	5

Table 3 -CPU times to simulate 100 s of physical time (on one core)

 3

	Computation time

Table 4 -CPU times to simulate 2.5 s of physical time (on forty cores)

 4

		Computation time
	Rotating mesh	2 h
	Rotating frame	1.4 h

4. Solid Random Kinetic Energy Transport Equation.

							t g		∂U g,i ∂x j	+	∂U g,i ∂x j	+	2 3	k g + t g	∂U g,k ∂x k	ı ij	(A.29)
	where the turbulent viscosity is given by:
	t g =	2 3	k g		t g 1 + C 12	␣ p p ␣ g g	t gp F gp	1 -	q gp 2k g	-1	(A.30)
	with C 12 = 0.34.
	Moreover t g is given by:
	t g = C	3 2	k g					(A.31)
	where C = 0.09.
	The transport equation for the turbulent dissipation reads:
	˛g g		∂ ∂t	+ U g,j	∂ ∂x j			= -˛g g k g	C ,1 u g,i u g,i	∂U g,i ∂x j	+ C ,2
										+	∂ ∂x j	˛g g	∂x j g ∂ t	+ p→g	(A.32)
	with = 1.3 and where
	p→g = C ,3 k g	kg p→g	(A.33)
	with C ε,3 = 1.2.
	The solid random kinetic energy transport equation writes:
	˛p p		∂q 2 p ∂t	+ U p,j	∂q 2 p ∂x j	= -	∂ ∂x j	˛p p K kin p + K col p	∂q 2 p ∂x j
											-	g p,ij	∂U p,i ∂x j
											-	˛p p gp F	2q 2 p -q gp

Appendix A.

Acknowledgements

This work was granted access to the HPC resources of CALMIP supercomputing center under the allocation 2015-0111. This work was performed using HPC resources from GENCI-CINES (Grant 2015-c20152b6012).

Appendix A. Mathematical Model

This appendix gives the set of equations of the n-fluid Eulerian model. The Eulerian n-fluid approach used is a hybrid method, in which the trans-port equations are derived by using phase ensemble average for the fluid phase and by using kinetic theory of granular ows supplemented by fluid effects for the dispersed phase, thanks to joint fluid-particle probability density function (PDF) approach [START_REF] Simonin | Combustion and turbulence in two-phase flows[END_REF]. In the following equations,subscript k = g refers to the gas phase and k = p to the particle phase. ˛p p in the particle transport equations represents n p m p where n p is the number density of particle centers and m p is the mass of a single particle: ˛p = n p m p / p is an approximation of the local volume fraction of particle p. Hence, gas and particle volume fractions, ˛g and ˛p have to satisfy:

Volume Fraction Transport Equation

The mass balance equation writes

where ˛k is the volume fraction of the phase k, k the material density and U k,i the i th component of the mean velocity. In Eq. (A.2) the right-hand-side is equal to zero because no mass transfer takes place.

Appendix A.2. Mean Momentum Transport Equation

The mean momentum transport equation writes

where P g is the mean gas pressure and g i the gravity acceleration. k,ij the effective stress tensor which is composed of two parts:

where u k,i is the velocity fluctation u k,i u k,j are the turbulent or kinetic stresses and k , ij the viscous or collisional stresses for the gas phase or the particulate phase respectively. Gas-Particle Interphase Momentum Transfer. In Eq. (A.3), I p,i is the mean gas-particle interphase momentum transfer without the mean gas pressure contribution. According to a large particle to gas density ratio only the drag force must be accounted for in the mean gas-particle interphase momentum transfer term, which is written as:

and I g,i = -I p,i .

(A.5)

The mean fluid-particle relative velocity, V r,i , is given in terms of the mean gas and solid velocities:

where V d is the drift velocity, which can appear due to turbulence [START_REF] Gobin | Fluid dynamic numerical simulation of a gas phase polymerisation reactor[END_REF] or sub grid effect [START_REF] Parmentier | A functional subgrid drift velocity model for filtered drag prediction in dense fluidized bed[END_REF].

Following [START_REF] Gobin | Fluid dynamic numerical simulation of a gas phase polymerisation reactor[END_REF] the drift velocity due to turbulence is modelled by:

where t gp the characteristic timescale of turbulence 'seen' by the particles (including crossing trajectory effects) and q gp= u g,i u p,i the gas-particle covariance.

The particle relaxation time scale writes 1

where Cd is the drag coefficient and with

The variance of the fluid-particle relative veloticy is given by: v r,i v r,i = 2q 2 p + 2k g -2q gp (A.10)

To take into account the effect of large solid volume fraction [START_REF] Gobin | Fluid dynamic numerical simulation of a gas phase polymerisation reactor[END_REF] proposed the following correlation for the drag coefficient .11) where C d,Erg is the drag coefficient proposed by [START_REF] Ergun | Fluid flow through packed columns[END_REF]:

and C d,WY by [START_REF] Wen | Mechanics of fluidization[END_REF]:

The particle Reynolds number is given by

Moreover, t gp is related to the characteristic timescale of turbulence, t g , through the following relation:

with C ˇ= 1.8 and where t g is given by the turbulence model. Solid Stress Tensor. The solid stress tensor writes p,ij = g p,ij + f p,ij , where f p,ij is the frictional stress tensor described in section 2. In the frame of the kinetic theory of rapid granular ows modified by the drag force [START_REF] Boëlle | Second-order prediction of the particle-phase stress tensor of inelastic spheres in simple shear dense suspensions[END_REF] the granular stress tensor can be computed from separate transport equations. However, for numerical implementation we use the following approximated form for the granular stress [START_REF] Boëlle | Second-order prediction of the particle-phase stress tensor of inelastic spheres in simple shear dense suspensions[END_REF][START_REF] Balzer | Eulerian gas-solid flow modelling of dense fluidized bed[END_REF][START_REF] Balzer | Gas solid flow modelling based on the kinetic theory of granular media: validation, applications and limitations[END_REF] g p,ij = P pp ∂U p,k ∂x k ı ijp S p,ij (A.16) where the strain rate tensor S p is defined by Eq. ( 2). The granular pressure, viscosities and model coefficients are given by P p = 2 3 ˛p p q 2 p 1 + 2˛pg 0 (1 + e c ) (A.17

The collision time scale c is given by [START_REF] Simonin | On the spatial distribution of heavy-particle velocities in turbulent flow: from continuous field to particulate chaos[END_REF]:

where gp is a fluid-particle correlation coefficient that takes into account the effect of turbulence on collisions and is given by: gp = q 2 gp 4q 2 p k (A.25)

Moreover, the radial distribution function, g 0 , is computed according to [START_REF] Lun | The effects of an impact velocity dependent coefficient of restitution on stresses developed by sheared granular materials[END_REF] as g 0 (˛p) = 1 -˛p ˛max -2.5␣m a x (A.26) where ˛m a x = 0.64 is the closest random packing.

Appendix A.3. k -Model Adapted to Gas-Particle Flows

The transport equation for the gas turbulent kinetic energy reads :