DEM/CFD Simulations of a Pseudo-2D Fluidized Bed: Comparison with Experiments

: The present work investigates the performance of a mesoscopic Lagrangian approach for the prediction of gas-particle ﬂows under the inﬂuence of different physical and numerical parameters. To this end, Geldart D particles with 1 mm diameter and density of 2500 kg/m 3 are simulated in a pseudo-2D ﬂuidized bed using a Discrete Element Method (DEM)/Large-Eddy Simulation (LES) solver called YALES2. Time-averaged quantities are computed and compared with experimental results reported in the literature. A mesh sensitivity analysis showed that better predictions regarding the particulate phase are achieved when the mesh is ﬁner. This is due to a better description of the local and instantaneous gas-particle interactions, leading to an accurate prediction of the particle dynamics. Slip and no slip wall conditions regarding the gas phase were tested and their effect was found negligible for the simulated regimes. Additional simulations showed that increasing either the particle-particle or the particle-wall friction coefﬁcients tends to reduce bed expansion and to initiate bubble formation. A set of friction coefﬁcients was retained for which the predictions were in good agreement with the experiments. Simulations for other Reynolds number and bed weight conditions are then carried out and satisfactory results were obtained.


Introduction
Based on their effectiveness regarding gas-solid heat and mass transfers, fluidized beds are among the best options for developing economically and environmentally viable techniques for fossil-fuel-based energy generation. Such systems involve complex physical mechanisms such as momentum, heat and mass exchanges between the gas and the particulate phases. In addition, fluidized beds exhibit an unsteady and inhomogeneous behavior leading to wide characteristic length and time scales. Particle-particle and gas-particle interactions at the micro-scale (1 to 5d p , where d p is the characteristic length of a particle) result in meso-scale structures, such as bubbles (10 to 100d p ), which can affect the macro-scale gas-particle flow [1]. Further, different local behaviors can be observed depending on the local particle-particle, gas-particle and particle-wall interactions, and may profoundly modify the bed hydrodynamics. Therefore, various numerical approaches have been developed over the past decades to simulate those flows at microscopic, mesoscopic and macroscopic scales with the aim of elucidating the mechanisms underlying the origin and the evolution of the heterogeneous gas-particle flow pattern. The accurate prediction of the underlying physics makes possible to improve existing processes and to design more efficient new facilities. In this context, the development of reliable numerical approaches is an essential prerequisite.

Modeling Strategy
All the numerical simulations presented in this work are performed using the finite-volume code YALES2 [9], a Large-Eddy Simulation (LES) and Direct Numerical Simulation (DNS) solver based on unstructured meshes. This code solves the low-Mach number Navier-Stokes equations for turbulent reactive flows using a time-staggered projection method for constant [10] or variable density flows [11]. YALES2 is specifically tailored for solving these low-Mach number equations on massively parallel machines with billion-cell meshes thanks to a highly optimized linear solvers [12].
Recently, a meso-scale four-way coupling approach for the modeling of solid particles has been implemented in the YALES2 solver. This approach combines DEM approach to represent the solid phase with LES equations solved on an Eulerian unstructured grid for the fluid phase. The CFD/DEM solver has been thoroughly optimized for massively parallel computing. It features a dynamic collision detection grid for unstructured meshes and packing/unpacking of the halo data for non-blocking Message Passing Interface (MPI) exchanges.

Gas Phase Modeling
The governing equations for the averaged fluids flow are obtained from the filtering of the unsteady, low-Mach number Navier-Stokes equations, taking the local fluid and solid fractions into account. If G is a filtering kernel (see for instance [7]), the local fluid fraction ε is defined as: where V f is the volume occupied by the fluid. Defining Φ(x, t) as a function of position and time, the volume filtered field < Φ > (x, t) refers to the regular spatial average and is computed by taking the convolution product with the filtering kernel G, giving: Further details concerning the volume filtering operations can be found in [8]. The governing equations finally read: ∂ε ∂t + ∇ · (ε < u >) = 0, ρ ∂ ∂t (ε < u >) + ρ∇ · (ε < u > ⊗ < u >) = −∇ < P > +∇ · (ε < τ >) + ερg + 1 ∆V < F > p→ f , (4) u, ρ, P, and τ are the gas velocity, density, dynamic pressure respectively and viscous stress tensor. The closure used to calculate the turbulent viscosity µ t is the localized dynamic Smagorinsky model [13] proposed by [14,15]. The termF p→ f is the momentum source term due to particle displacement. ∆V is the local control volume. Details concerning the computation of this term can be found in Section 2.3. For the sake of clarity, the notation <> will be dropped for the averaged quantities in the following. The use of the LES model for CFD/DEM simulations in case of dense granular flows is questionable. Indeed, there are very few publications dealing with the effect of inertial particles on a turbulent flow in a confined domain. It is thus very difficult to conclude whether the actual flow is turbulent or not. Nevertheless, previous publications (see, for instance, [8]) used this kind of model in fluidized bed simulations. In our case the bulk Reynolds number based on the bed width is relatively high (5600 ≤ Re ≤ 8000) so that, a single gas-phase flow in the same condition, is likely to be turbulent or at least in a transition regime. We performed a series of additional computations in order to assess the effect of the turbulent viscosity on the bed hydrodynamics. Results showed that, both qualitatively and quantitatively, such a model did not change significantly the hydrodynamics of the fluidized bed (see Appendix A). This tends to confirm that in the present configuration, the global behaviour of the flow is driven by the largest scales and by particle-particle contact.

Discrete Particle Modeling
The forces acting on a particle in motion can be divided into two categories, volume and surface forces. The unique volume force in the present situation is the weight whereas, the surface forces consist of hydrodynamic and contact forces. While the hydrodynamic forces arise from fluid-particle interactions, such as drag and lift forces, the contact forces are due to particle-particle and particle-wall interactions. Such forces can be further classified in collision and adhesive forces. In the present work, due to the high solid/gas density ratio, as reported in Table 1, buoyancy, lift and added mass forces are neglected. To date, theoretical expressions for the history force in dense regimes are still missing. Therefore, such a force is not taken into account in our modeling. Concerning the contact forces, only collision forces are accounted for. Adhesive effects that may originate from electrostatic or Van Der Vaals forces can safely be here neglected since, as it can be seen in Table 1, the simulated particles are relatively large (Geldart D particles).
Discrete particle models, in which each particle is tracked is tracked in a Lagrangian fashion, have clearly shown their ability to simulate the behavior of granular flows, and originate from molecular dynamics methods initiated in the 1950s [3]. A soft-sphere model [16] is employed to compute contact between each particle. They are allowed to overlap other particles or walls in a controlled manner. A resulting contact force accounting for particle-particle and particle-wall repulsion is thus added in the momentum balance of each particle. Particle movement is then given by Newton's second law for translation assuming high solid/gas density ratio: where m p , u p , and x p are the particle mass, velocity, position respectively. F D is the drag force, F G is the gravity force, F P is the pressure gradient force and F C is the collision force. The relation between F D , F C in Equation (5) and F p→ f in Equation (4) is detailed in Section 2.3. Preliminary simulations including particle rotation provided similar results as those without rotation. As a consequence, for computational reason, the particle rotation is not accounted for in this paper.

Modeling of Collisions
The total collision force F C acting on particle a is computed as the sum of all forces f col b→a exerted by the N p particles and N w walls in contact. As particles and walls are treated similarly during collisions, the b index refers to both: f col b→a = f col n,b→a + f col t,b→a .
Depending on the desired compromise between accuracy and numerical cost, a lot of soft sphere models can be found in the literature. Here, a linear-spring/dashpot [16] is used along with a simple Coulomb sliding model accounting for the normal f col n,b→a and tangential f col t,b→a components of the contact force, respectively, as in the work of Capecelatro [8]. For one particle (or wall) b acting on a particle a: f col n,b→a = −k n δ ab n ab − 2γ n M ab u ab,n if δ ab > 0, f col t,b→a = −µ tan ||f col n,b→a ||t ab .
The term δ ab is defined as the overlap between the a and b entities expressed using each particle radius r p and center coordinates x p such as: The system effective mass M ab is expressed using each particle mass m p such as: In case of a particle-wall collision, the wall is considered as a particle with infinite mass and null radius. This model requires three user-defined parameters; k n , γ n , and µ tan respectively accounting for the spring stiffness, normal damping, and friction coefficient of the a − b binary system. n ab and t ab respectively account for the unit normal vector from particle a towards entity b and the unit tangential vector. n ab is calculated as follows: with u ab being the relative velocity of the colliding system, its normal and tangential components are then given by: u ab,n = ((u a − u b ).n ab )n ab , (13) u ab,t = u ab − u ab,n .
Finally, t ab is given by: Using Newton's third law and a projection on n ab yields the following ordinary differential equation for the overlap evolution of an a − b binary system undergoing collision without taking any other force into account: where ω 0 stands for the system's natural frequency and is defined as: The damping parameter γ n accounting for the energy dissipation occurring during contact is calculated by the mean of a normal restitution coefficient e n verifying 0 < e n ≤ 1 such as: A contact time T C can also be analytically extracted from Equation (16) corresponding to the time during which the particle a and the entity b are overlapping:

Closure for Drag
The drag force F D acting on a particle p is written: where u @p is the local gas velocity interpolated at the center of the particle p and τ p is the drag relaxation time expressed as follows: where ρ p is the particle density, d p its diameter, µ the dynamic viscosity, ε @p is the local fluid fraction interpolated at the center of the particle, and C D is the drag coefficient. In order to compute C D , the classical closures proposed by Ergun (C D,ER ) [17] for high ε @p values and Wen & Yu (C D,WY ) [18] for low ε @p values are used along with the smoothing function φ gs introduced by Huilin & Gidaspow [19] to avoid discontinuities when switching models: in such a way that: and

Other Forces
The gravity force F G acting on a particle p is written: The pressure gradient force F P in Equation (5) reads as: where V p is the particle's volume and ∇P @p is the local pressure gradient interpolated at the center of the particle.

Phase Coupling
The coupling between the particle and fluid phases is a key point for the modeling of particle-laden flows, especially when the particle size approaches the Eulerian cell size. Many Eulerian fields have to be interpolated at the center of the particles for the numerous closures, as shown in Section 2.2. In the YALES2 solver, particles are located in a unique mesh cell (C) using the position of their center. For any Eulerian scalar or vector field Φ(x, t), its value taken at the particle p center Φ @p (t) obeys: Here i is a node index so that 'i ∈ C' means 'all nodes composing the mesh cell C in which the particle p is located'. ω p,i is the interpolation weight of the particle p on cell node i and is calculated using a trilinear interpolation on tetrahedra and on hexahedra. The same interpolation weights are used for data transfer from grid to particles (interpolation step) and from particle onto the grid (projection step). The conservative projection operator needed to compute F p→ f (see Equation (4)) is thus written on each node i as: and the fluid fraction at node i is written: Here ∆V i denotes the control volume of node i and 'p ∈ SC i ' means 'all particles belonging to any surrounding cell (SC) of node i'. Referring to Figure 1, it means that any particle belonging to one of the cells will be accounted for when computing ε (as well as F p→ f ) at node 3. Figure 1. 2D representation of a particle p located in a triangular cell ( ) moving towards a neighboring cell ( ). : mesh nodes. : contour of node 3 control volume. ω p,i : interpolation weight of particle p on node i. This method consisting in distributing particle quantities only in the cell where its center resides, referred to as Particle Centroid Method (PCM), can lead to large calculation errors in particular regarding the fluid fraction, as pointed out in [20]. This is partly due to the fact that many CFD/DEM codes feature a staggered grid where the fluid fraction is defined at cell centers, causing dramatic discontinuities in time and space derivatives when a particle enters or leaves a cell. On the contrary, in the YALES2 code the fluid fraction, as all the Eulerian fields, is computed at the grid nodes. As depicted in Figure 1, it is then straightforward that the particle crossing from the green cell to the red one won't cause any discontinuity on the computation of neither ω p,1 nor ω p,3 involved in their interface. Moreover ω p,2 won't be much affected neither during the crossing because as the particle approaches the cell's interface, ω p,2 tends towards 0. This is still true in 3D cases and on Cartesian meshes. Nevertheless, it is well known that the PCM method can induce inaccuracies and lead to numerical instabilities because it cannot prevent the fluid fraction from reaching unrealizable values, in particular when dealing with close to unity particle diameter/mesh cell size ratios. In such cases, the fluid fraction value can locally decrease below the theoretical packing limit. To cope with this limitation, a filtering operator well suited for distributed memory machines is used. Taking a 2D case as shown on the left in Figure 2, this filtering operator is built for any Eulerian scalar or vector field noted Φ i , its filtered value beingΦ i . At node i 1 ,Φ i reads: where ∆V i 1 is the control volume associated to node i 1 and S mn is the part of ∆V i 1 contained in the face delimited by nodes i 1 , i m and i n , as shown on the left in Figure 2. If all the control volumes are equal, on a structured mesh for instance, Equation (31) becomes: Figure 2. On the left: 2D representation of an unstructured mesh. : mesh nodes. The control volumes of nodes i 1 ( ), i 2 ( ), i 3 ( ), i 4 ( ), i 5 ( ), i 6 ( ), i 7 ( ) are shown. The control volume of node i 1 is composed of surfaces S 23 , S 34 , S 45 , S 56 , S 67 and S 72 . On the right: 3D representation of a regular Cartesian mesh part. : node i 1 . : nodes at faces' center. : nodes at edges' center. : nodes at corners.
The same type of filter can be derived in 3D. The following equation gives the value ofΦ i 1 in a 3D structured case with all equal control volumes as shown on the right in Figure 2, where RN is the abbreviation of Red Nodes, BN of Blue Nodes and GN of Green Nodes: This fully conservative operation being performed on all volumes at the same instant provides a fully filtered field, and can be repeated several times to increase the filter width. In all the simulations presented in this paper, only one filtering step was used. It should be underlined that for the computation ofε, the filtering step is applied before dividing by the local control volume in order to conserve total solid mass over the whole computational domain volume V: The properties of such a filtering operator, i.e., its moments, are not straightforward to determine on unstructured meshes but it can be noticed that it is based on direct neighbors and thus doesn't need distant nodes, hence its attractiveness regarding parallelism. The main drawback is that the filter width can't be directly obtained because it depends on the local mesh size. Thus, when using this filtering operator, the user cannot prescribe the filter width. Other projection methods have been developed to circumvent this drawback while ensuring interesting mathematical properties as moment conservation for instance. One can cite the work of Capecelatro [8], who used a Gaussian kernel filter in a method called mollification for the same kind of application as presented here, and the work of Mendez [21] in which a projection operator based on high-order moments conservation was built for deformable red cell membrane modeling purposes.

Fluid Advancement Procedure
This section presents some numerical features of the YALES2 code. This code solves the filtered low-mach number Navier-Stokes equations presented in Section 2.1 with an explicit time advancement. Among the various implemented numerical schemes, a fourth-order central scheme was used for the spatial integration, and a fourth-order scheme called TFV4A [22] combining Runge-Kutta and Lax-Wendroff methods was used for the explicit time integration.
The time advancement uses a time-staggered projection method for variable density flows [11] described below in which the n superscript refers to discrete times:

Lagrangian phase advancement
First, the particles are advanced. The full description of this step is available in Section 2.4.2. After being relocated on the grid, ε n+3/2 can be computed using Equation (30).

Density prediction for scalar advancement
The density predictor (ερ) is then determined by the mean of the mass conservation equation (Equation (3)): 3. Velocity prediction Once (ερ) n+1 is known, the velocity can be predicted reusing the dynamic pressure gradient of the previous time step: 4. Velocity correction Velocity correction is performed by updating the pressure gradient: The Poisson equation aiming at calculating P n+1/2 2 is obtained by taking the divergence of Equation (37) and inserting the condition imposed by the following equation of mass conservation written for u n+1 : The Poisson equation finally reads: This linear system requires an efficient and accurate iterative solver. For all our simulations, a Deflated Preconditioned Conjugate Gradient (DPCG) algorithm [23] is used.
The resulting time advancement is fully mass and momentum conserving. The time step for the fluid phase ∆t is calculated at each solver iteration by enforcing a maximum value of 0.2 for the Courant-Friedrichs-Lewy number (CFL) criterion and 0.15 for the Fourier number criterion for all the simulations.

Particle Advancement Procedure
A second-order explicit Runge-Kutta (RK2) algorithm is used to advance the position x p , the velocity u p (see Equation (5)) of the particles in time: RK2 -1st step: RK2 -2nd step: where ∑ F refers to the summation of all forces acting on particle p (see Section 2.2). At the particle scale, several phenomena need to be integrated properly by the mean of an associated characteristic time, among these: drag, gravity, etc. Thus, several stability criteria have to be computed on each particle to determine the smallest time step needed for the most constraining characteristic time. In dense fluidized beds simulations the collision time step is generally the limiting one, so it will be noted ∆t p from now on.
∆t p must be inferior to the contact time T C described in Equation (19) in order to be able to solve collisions properly. Furthermore, it must be small enough to ensure numerical stability depending on the selected numerical scheme, without compromising the performances of the code. It can be noted that there is a unique value of T C because all the particles are identical. In our simulations, the following criterion was used: This criterion has been tested against the following non-dimensional analytical solution of Equation (16) that can be found in [24] for a normal collision of a single particle on a wall with given collision parameters as shown in Figure 3: with  The time evolution of both overlap and particle velocity during contact has been plotted in Figure 4a,b, respectively.  Figure 4. Non-dimensional overlap (a) and velocity (b) as a function of the non-dimensional time during the contact presented in Figure 3. Comparison between analytical method ( ), Euler method with t p = T C /80 ( ) and RK2 method with t p = T C /6 ( ).
The agreement with the analytical solution is not perfect. However the scheme is stable and in this case the absolute error on the restitution velocity was found to lie below 5%.

Performances and Parallelism
A critical point in the coupling between the gas and solid phase is the time advancement, as they can have very different time scales. The choice was made to sub-step the solid phase, of which time advancement is generally limited by the collision time step in dense fluidized beds, during a gas phase time step, which is itself governed by convective and diffusive time scales. The mean number of particle sub-steps observed during our simulations was approximately 3 for the refined mesh, 7 for the intermediate mesh and 16 for the coarse mesh (mesh details are described in the following section).
As the computation of the collision force requires a distance estimation for each particle pair a priori, solving Equations (40) and (41) can become critically time consuming. The work of Lubachevsky [25] provides an analysis of the linked-cell method to tackle this problem, which is the most commonly used method. The solution consists in the definition of a Cartesian grid superimposed on the grid mesh allowing each particle to search for its potential collision partners, as shown in Figure 5. First, each particle is located in one of the Cartesian grid cell. Then looping over the 8 surrounding cells plus the one where the particle stands (26 + 1 cells in 3D) provides a list of the closest particles, that are stored as potential collision partners. Eventually, during the computation of the collision force, only these particles are checked. : superimposed Cartesian grid. The particle of interest is noted p in the central cell ( ). The particles of which center belongs to this cell or one of the surrounding ones ( ) are stored as potential collision partners of particle p.
The linked-cell method is used even if the grid is Cartesian. This choice is justified by the fact that the cells of the detection grid should probably be (i) at least as big as the particle diameter (ii) as small as possible, in order to prevent too many neighboring particles from being detected. These requirements may be hard to guarantee if the Cartesian mesh is too coarse, or locally refined (which is one of the future aims of the code). For these reasons, a Cartesian grid is superimposed in any case.
Parallelism also requires special treatment for particles, as collision might occur between some of them although they don't belong to the same processor domain. To cope with this requirement, a ghost particle method is used, which is also classical in CFD/DEM. Ghost particles are identified using a cell halo surrounding each processor domain as shown in Figure 6 for processor ranked #1 in a cylindrical geometry discretized with an unstructured mesh. The particles belonging to the cell halo are exchanged between involved processors with the necessary data to compute collision force only by the mean of non-blocking Message Passing Interface (MPI) exchanges. The size of the cell halo is locally determined by the cell mesh size and particle diameter to avoid unwanted distant particles from being identified as ghost particles, which is not straightforward on unstructured meshes. On each processor, ghost particles are treated to be located as well on the afore mentioned Cartesian detection grid (see Figure 5).  Figure 6. Ghost particle method principle shown for a part of a cylindrical domain. The different processor domains are colored accordingly. The processor of interest is ranked #1 ( ) and its closest neighbors are ranked #2, #3, #4, #5 and #6. A particle entering the white cell halo around #1 will be sent to #1 as a ghost particle by the processor it belongs to.

Assessment of the Solver Performances
Assessments of the code global performances are first extracted from a canonical isothermal case, disregarding its physical relevance. This case consists of a cubic box meshed with tetrahedra. Particles are randomly seeded in the box, with a mean porosity of about 0.54, and each tetrahedron contains roughly 11 particles. A fluid phase is present to account for the computational cost related to interpolation and projection steps. The particle timestep is chosen such that, ten particle timesteps are performed for each fluid timestep, which corresponds to a usual sub-stepping configuration. Code assessments are obtained from a single fluid timestep. All the tests were carried out on the regional supercomputer Myria of the CRIANN center (Centre Régional Informatique et d'Applications Numériques de Normandie, France), featuring a Intel Omni-Path interconnect. The processors used are Intel Broadwell 14 cores running at 2.4 GHz with 128 GB RAM (about 4 GB memory per core) for total peak power of 400 TFlop/s. The speed-up is first obtained by running the same simulation on different numbers of cores, ranging from 532 (reference case) to 4144. Each simulation roughly gathers 38M tetrahedra and 410M particles. The reference CPU time t re f CPU being associated with the temporal loop of the solver on N procs re f = 532 cores, the speed-up for a CPU time t CPU on N procs is calculated by: and is illustrated on Figure 7a. The solver exhibits a good scalability up to 4144 cores, with a speed-up reaching 80% of the ideal scaling. Secondly, the scale-up of the solver is quantified by an evalutation of the performances at constant load per core on different numbers of cores, ranging from 252 (reference case) to 4144. The number of particles per core is about 99k, and the number of tetrahedra is about 9.1k per core. The scale-up is given by: where N re f c and N c are the number of cells in the reference case and in the current case, respectively. The scale-up curve is represented on Figure 7b, which shows an excellent scaling.
B n l r M + 2 0 a S m d t 1 b 3 8 T f T K Z m 9 Z 7 l u R n e 9 S 1 p w N 7 P c c 6 D 1 p H j u Y 5 3 X q 3 U n H z U R e x h H 4 c 0 z x P U U E c D T f I e 4 h F P e L b q V m R l 1 t 1 n q l X I N b v 4 t q y H D 6 d k j 3 k = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 B n l r M + 2 0 a S m d t 1 b 3 8 T f T K Z m 9 Z 7 l u R n e 9 S 1 p w N 7 P c c 6 D 1 p H j u Y 5 3 X q 3 U n H z U R e x h H 4 c 0 z x P U U E c D T f I e 4 h F P e L b q V m R l 1 t 1 n q l X I N b v 4 t q y H D 6 d k j 3 k = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 B n l r M + 2 0 a S m d t 1 b 3 8 T f T K Z m 9 Z 7 l u R n e 9 S 1 p w N 7 P c c 6 D 1 p H j u Y 5 3 X q 3 U n H z U R e x h H 4 c 0 z x P U U E c D T f I e 4 h F P e L b q V m R l 1 t 1 n q l X I N b v 4 t q y H D 6 d k j 3 k = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 W X e f q V Y h 1 + z i 2 7 I e P g C p x 4 9 6 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " w B P P R U B 0 y i D + 5 y P P W X e f q V Y h 1 + z i 2 7 I e P g C p x 4 9 6 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " w B P P R U B 0 y i D + 5 y P P W X e f q V Y h 1 + z i 2 7 I e P g C p x 4 9 6 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " w B P P R U B 0 y i D + 5 y P P

= " > A A A C t X i c j V L L S g M x F D 0 d X 7 V W r W s 3 g 0 V w V T J u d C n o w m U F + 4 B a Z C Z N a + y 8 T D J C K f 6 A W z 9 O / A P 9 C 2 / i C G o R z T A z J + f e c 5 K b m y i P p T a M v V S 8 p e W V 1 b X q e m 2 j X t v c 2 m 7 U u z o r F B c d n s W Z 6 k e h F r F M R c d I E 4 t + r k S Y R L H o R d N T G + / d C 6 V l l l 6 a W S 6 G S T h J 5 V j y 0 B D V v m 4 0 W Y u 5 4 S + C o A R N l C N r P O M K I 2 T g K J B A I I U h H C O E p m e A A A w 5 c U P M i V O E p I s L P K B G 2 o K y B G W E x E 7 p O 6 H Z o G R T m l t P 7 d S c V
o n p V a T 0 s U + a j P I U Y b u a 7 + K F c 7 b s b 9 5 z 5 2 n 3 N q N / V H o l x B r c E P u X 7 j P z v z p b i 8 E Y x 6 4 G S T X l j r H V 8 d K l c K d i d + 5 / q c q Q Q 0 6 c x S O K K 8 L c K T / P 2 X c a 7 W q 3 Z x u 6 + K v L t K y d 8 z K 3 w J v d J f U 3 + N n N R d A 9 b A W s F V w w V L G L P R x Q G 4 9 w g n O 0 0 S H L E R 7 x 5 J 1 5 t 9 7 d x z 3 w K u W

O a g i H w 1 Z e c E t + p B g y D A B R w x N O E K A l J 4 u f H h I i O t h T p w i J G y c 4 w E l 0 m a U x S k j I H Z M 3 y H N
u j k b 0 9 x 4 p l b N a J W I X k V K F 0 e k k Z S n C J v V X B v P r L N h f / O e W 0 + z t x n 9 w 9 x r Q q z G i N i / d I v M / + p M L R o D n N k a B N W U W M Z U x 3 K X z J 6 K 2 b n 7 p S p N D g l x B v c p r g g z q 1 y c s 2 s 1 q a 3 d n G 1 g 4 2 8 2 0 7 B m z v L c D O 9 m l 9 R g / 2 c 7 l 0 G r X v O 9 m n / t o Y g D H O K Y 2 n i K c 1 y i g S Z Z j v C I J z w 7 V 4 5 0 p p 9 X w S n k d 2 I f 3 4 Z z / w E G 0 o 6 S < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " H 9 E N / z V X o 5 t 8 G y Q x s / Z f F M Y p e l Y = " > A A A C x 3 i c j V H L S s N A F D 2 N r 1 p f V Z d u g k V w F Z K C 6 L L g R n c V 7 A N q k S S d t k P T T J h M i q W 4 8 A f c 6 p + J f 6 B / 4 Z 0 x B b W I T k h y 5 t x 7 z s y 9 N 0 g i n i r X f S 1 Y S 8 s r q 2 v F 9 d L G 5 t b 2 T n l 3 r 5 Performances of the code are also assessed on an isothermal pressurized gas-fluidized bed carried out at the university of Birmingham [26]. The investigated case gathers 10M particles in a cylindrical domain meshed with 3.7M tetrahedra. Here, the performances are extracted over 1 s of physical time after the stability of the bed fluidization has been assessed by monitoring both the bed height and the pressure loss across the system. Thus, this case deals with realistic conditions where the execution speed highly depends on the local physics (presence of dense or void zones), and where, the fluid and particle timesteps are recomputed throughout the simulation. The tests were carried out on the Curie supercomputer of the TGCC center (Trés Grand Centre de Calcul, France), featuring a InfiniBand QDR Full Fat Tree interconnect. The used nodes comprise two Intel Sandy Bridge octo-core processors running at 2.7 GHz with 64 GB RAM (about 4 GB per core).

m K T I a s E Y p I y H b g p y z i M W s o r i L W T i T z x 0 H E W s H o X M d b E y Z T L u J r N U 1 Y d + w P Y t 7 n o a 8 0 V T r x q r f l i u u 4 Z t m L w M t B B f m q i / I L b t C D Q I g M Y z D E U I Q j + E j p 6 c C D i 4 S 4 L m b E S U L c x B n u U S J t R l m M M n x i R / Q d 0 K 6 T s z H t t W d q 1 C G d E t E r S W n j i D S C 8 i R h f Z p t 4 p l x 1 u x v 3 j P j q e 8 2 p X + Q e 4 2 J V R g S + 5 d u n v l f n a 5 F o Y 8 z U w O n m h L D 6 O r C 3 C U z X d E 3 t 7 9 U p c g h I U 7 j H s U l 4 d A o 5 3 2 2 j S Y 1 t e v e + i b + Z j I 1 q / d h n p v h X d + S B u z 9 H O c i a F Y d z 3 W 8 K 7 d S c / J R F 3 G A Q x z T P E 9 R w w X q a J D 3 E I 9 4 w r N 1 a Q l r Y t 1 9 p l q F X L O P b 8 t 6 + A A 5 p o + t < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " H 4 5 d I s v S o P 1 K a O E 3 b h p O Q H n F 2 u 8 = " > A A A C x 3 i c j V H L S s N
The speed-up is obtained by running the simulation on various numbers of cores, ranging from 64 (reference case) to 1024 cores, and is illustrated on Figure 8. On the whole, the solver reaches 55% of its ideal scaling value, which is acceptable given the high dispersion of the particles amongst the cores, causing their de-synchronizing.

Configuration and Meshes
In this work, CFD/DEM simulations of a fluidized bed, similar to that used in the experiments reported by Patil et al. [27], are performed. A sketch of the simulated configuration of the fluidized bed is shown in Figure 9 (left) including its dimensions. The bed contains inert glass particles fluidized by fresh nitrogen injected through the bottom of the bed at different gas flowrates. Physical properties of both gas and particles are summarized in Table 1. In the experiments, the bottom area of the bed was equipped with a small circular nozzle of 1.2 cm diameter, through which no gas was supplied during the experiments. In order to reproduce the effect of the switched nozzle as faithfully as possible, a zero gas flow was set through an area located at the center of the lower horizontal section of the bed. This area was considered as a wall for both gas and particles and is referred to as the bottom wall in the simulations. In the simulations, three different grid refinements were used and their characteristics are reported in Table 2, including that of the bottom wall retained for each mesh. As an example, a sketch of the intermediate mesh is illustrated in Figure 9 (right).  Table 1. Gas and particle properties [27].

Description of the Simulation Cases
Several simulations based on different grid refinements, combined with different physical parameters, are investigated. In our modeling strategy, a friction coefficient, previously denoted as µ tan in Section 2.2.1, is required to account for the tangential contact force during binary particle and particle-wall collisions. Since such a coefficient is not provided in the work of Patil et al. [27], we examined available values from the literature. As reported by Lorenz et al. [28], its value may depend on the state and history of the particle surfaces. The authors experimentally evaluated the friction coefficient, of glass particles having experienced a considerable time in granular flows, to 0.177 and 0.141 for inter-particle and particle-spent aluminum plate collisions, respectively. In the work of Goldschmidt et al. [29], values of 0.10 and 0.09, as measured by Gorham et al. [30], were respectively used in the modeling of the particle-particle and the particle-wall frictional collisions in a dense gas-fluidized bed configuration. In Gorham et al., both the particles and the pseudo-2D bed were made of glass. In the following sections, we denote the particle-particle and the particle-wall friction coefficients as µ p and µ w , respectively. All the simulations are summarized in Table 3. A mesh sensitivity analysis, drawn from simulations C1, C2 and C3, for which µ p = µ w = 0.2 is assumed, is first presented. In a second step, various values of the friction coefficients, µ p and µ w , are investigated.
This study refers to the sets of simulations C4 and C5, as shown in Table 3. The C4 simulations are performed using a constant value of µ p , with different values of µ w , whereas in the C5 simulations, the effect of µ p is investigated for a constant value of µ w . From these investigations, values for the friction coefficients are chosen by comparing numerical results of the volume fraction and the axial flux of the particulate phase with the available experimental measurements. Finally, additional Reynolds number and bed weight conditions are used as variable parameters for the simulations C6 and C7. It has to be mentioned that the inlet gas velocity, U f , given in Table 3 is the superficial velocity through the horizontal cross-sectional area of the bed (so, this area also includes the bottom wall). For each tested gas velocity, the corresponding particle Reynolds number, given as Re p = ρU f d p /µ, is reported in the table, together with the Stokes number St = ρd p U f /(18µ). It has to be noted that, the bulk Reynolds number based on the bed width is relatively high (5600 ≤ Re ≤ 8000). Therefore, for all the performed simulations, a turbulent viscosity based on Smagorinsky model is used. The additional parameter of interest is the restitution coefficient during particle-wall collisions, e w , which is considered equal to that of the particle-particle collisions, e n .

Results
Non-dimensional results are presented in this section. Since the lengths of the simulated configuration are different in all three directions a different length scale is defined in each direction. As a consequence, x−, y− and z−distances are divided by the bed width (L = 0.08 m), the bed depth (D = 0.015 m) and the bed height (H = 0.25 m), respectively. The time scale is defined as t/(d p /U f ). The velocities are divided by the inlet gas velocity (U f ) and the mass flowrates by the inlet gas mass flowrate (ρU f ). In order to characterize the bed hydrodynamics in an Eulerian fashion, the instantaneous fields of any Lagrangian variable, ψ p , are obtained by performing a spatial average over the grid cell. In addition, a 2D-behavior is considered for this fluidized bed configuration since the hydrodynamic variations in the depth direction are negligible. This can be observed in Figure 10, which displays the instantaneous particle volume fraction, α p , on different slices chosen at y/D = 0.2 (close to the back side), y/D = 0.5 (in the middle of the bed) and at y/D = 0.8 (close to the front side).
As an example, in Figure 11, profiles along the depth direction at the height z/H = 0.16 are also given for the time-averaged vertical velocity of the solid phase normalized by the inlet gas velocity (U p,3 /U f ). It can be observed that variations along the depth direction of such a quantity are not strong. Thus, two-dimensional variables are computed by averaging over the bed thickness (i.e., in the y-direction). In the present work, the 2D-variables corresponding to the computed field, < ψ p >, are denoted as < ψ p > xz . All the simulations are run for 20 s in order to reach a steady regime. In such a regime, particles and fluid hydrodynamic fields are statistically stationary and time-averaged quantities, < ψ k >, for which k denotes the gas (g) or the particulate (p) phase, may be computed for comparison with the experimental measurements.

Effect of the Grid Refinement
In this section, the simulations C1, C2 and C3 are presented and discussed. The three simulations differ from each other by their mesh size, as shown in Tables 2 and 3, and consequently by their filtering kernel length ( 4∆x). However, because the kernel length is always larger than the particle diameter, the volume-averaged equations are well posed in the sense of Jackson [7]. Snapshots of the time-averaged fields concerning the volume fraction, < α p > xz , and the mass flux, < α p ρ p U p > xz , of the particulate phase are depicted in Figure 12 for the different grid refinements. Experimental results are also included in the right column. The simulations match quite well the experiments. The particle volume fraction exhibits high values (dense regions) close to the side and bottom walls and small values (dilute regions) in the center of the bed. However, the average height of the fluidized bed seems to be slightly overestimated by the simulations. The solid flux exhibits a macroscopic double mixing loop and, on average, the particles move upward at the center of the bed and downward close to the side walls. Comparison to experiments reveals that the center of the mixing loops, as predicted by the simulations, is located slightly higher than that experimentally observed. The three simulations indicate clearly that the mesh size affects the bed hydrodynamics prediction. The reason is that particles continuously interact with the gas flow through the drag and the pressure gradient forces, which then contributes to the particle dynamics evolution. Since the mesh resolution affects the hydrodynamics prediction of the gas phase, this consequently modifies the particle dynamics predictions. The upper panel of Figure 12, shows that the finest grid leads to the best predictions, as a result of an accurate resolution of the gas hydrodynamics. The high solid concentrations close to the side walls, as reported by the experiments, are reproduced with a relatively good accordance by the fine-grid simulation. Nevertheless, refining the mesh does not seem to have a significant effect on the height of these dense regions from the bottom of the bed, nor on the mean bed height. Indeed, refining the mesh "locally" improves the bed hydrodynamic predictions by improving the fluid velocity field around the solid particles, but it does not influence the predicted overall bed height. This latter is instead directly affected by the selected drag force model. This is the main reason why one-dimensional models are accurate regarding the prediction of the mean bed height, given a sufficiently appropriate drag model. The analysis of the time-averaged fields of the solid flux reveals, apart from the progressive intensification of these fields due to the increasing mesh resolution, that this latter did not influence the position of the observed double mixing loop.  Figure 13 shows profiles of the time-averaged vertical solid mass flux normalized by the gas mass flowrate, < α p ρ p U p ,3 > xz /ρU f . The profiles are taken at the height z/H = 0.092 (z = 2.3 cm) above the bottom of the bed, for which experimental data are available. We notice that the predictions are improved when the finest mesh is used, especially close to the side walls. However, at the center of the bed, the effect of the bottom wall on the axial solid flux is still poorly reproduced even when a finer mesh is used.
With the attempt to improve the numerical results, the wall boundary condition type on the gas phase was also investigated. Within the framework of Navier-Stokes equations, the wall boundary condition for the gas is no-slip. However, such a condition is questionable for averaged equations, as those employed in this work. In addition, the presence of the particles may profoundly modify the flow hydrodynamics. In the present work, slip and no-slip fine-grid simulation cases were run and compared to each other. The time-averaged results (not shown here) relative to the particulate phase were only very slightly affected by the selected wall condition of the gas phase. Such a behavior is inherent to inertial particles, as those simulated in the present work, for which the Stokes number is much greater than unity, as shown in Table 3. Figure 13. Time-averaged normalized vertical mass flux of the solid phase, < α p ρ p U p ,3 > xz /ρU f , at the height z/H = 0.092 above the bottom bed, for different grid refinements. µ w = µ p = 0.2, Re p = 70 and bed weight 75 g.

Investigation of the Frictional Effects
Effects of the inter-particle and the particle-wall friction coefficients on the bed hydrodynamics are here discussed. In the previous study, the computational time required by the fine-grid simulation was nearly twice higher than that of the intermediate-grid simulation. Therefore, the present study is carried out using the intermediate mesh.
The time-averaged fields concerning the particle volume fraction, as obtained from the sets of simulations C4 and C5 are shown in Figure 14. As previously mentioned, a dense region is observed close to the side and bottom walls, whereas, the center of the bed is characterized by a dilute solid concentration. In addition, the average height of the fluidized bed seems to decrease as the frictional effects become stronger. However, more subtle differences between the situations can be observed. The upper panel as obtained by the simulations C4, reveals that, for the same value of the particle friction coefficient, when the wall frictional effects become stronger, the dense regions become denser and larger near the side walls, while they are sligthly narrowing above the bottom wall. Results also show that such regions move slightly upwards as µ w increases. It has to be reminded that the influence of the wall friction coefficient is also exerted at the front and the back walls, but should be stronger in the confined parts of the bed, namely at the side walls, where the denser regions of the bed are observed. The bottom panel, obtained by the simulations C5, shows that, when µ p increases up to 0.2, the bed hydrodynamics exhibit similar behavior as that observed when µ w varies and µ p is kept constant (simulations C4). Such results are in agreement with the work of Yang et al. [31], who showed the same trend of particle volume fraction distribution when µ p was increased from 0.05 to 0.15, in their comparison study between two fluid model (TFM) and discrete particle model (DPM) simulations. However, unlike simulations C4, continuing to increase the particle-particle friction coefficient above the value of 0.2 in the simulations C5 resulted in decreasing the solid concentration near the side walls and increasing it above the bottom wall. The inter-particle friction acts on the whole bed volume and, as a consequence, high values of the inter-particle friction coefficient may considerably modify the bulk-bed hydrodynamic behavior. Finally, when the particle friction is not considered (µ p = 0 for µ w = 0.1), a homogeneous fluidized bed is obtained, which is not the case for µ p = 0 and µ w = 0 (simulations C4). These various bed behaviors inform about the effect of the particle friction on the bubble formation. As an exemple, snapshots of the instantaneous bed porosity (gas volume fraction), obtained at different instants by the simulations C5 for two values of µ p , are depicted in Figure 15. In the case of frictionless particles (µ p = 0), no bubble formation occurs, whereas in the case of relatively frictional particles (µ p = 0.2), a realistic bubble growth is observed. This is consistent with the work of Hoomans et al. [32], who showed the same bubble patterns in 2D-fluidized bed DPM simulations. Figure 14. Time-averaged particle volume fraction, < α p > xz , for different values of the particle-wall µ w and the particle-particle µ p friction coefficients. Upper panel: constant µ p and different µ w . Bottom panel: constant µ w and different µ p . Re p = 70 and bed weight 75 g. Time-averaged particle agitation energy, q 2 p , has been computed as q 2 p = u p,i u p,i , where u p,i = u p,i − U p,i . In these expressions, u p,i and U p,i are the i th components of the fluctuating and the mean particle velocities, respectively, with the mean particle velocity computed as U p,i =< u p,i >. The normalized results are shown in Figure 16 for the two sets of simulations. The first panel obtained from the set C4 shows that most of the particle agitation is produced at the approximate height z/h = 0.25. This height corresponds to the average height of the dense bed (α p > 0.2). Furthermore, it increases when friction at wall becomes stronger. Results of the simulations C5, as depicted in the bottom panel, show a decline, on average, of the particle agitation energy when the particle friction increases. In addition, some particle agitation is produced at the height z/h = 0.25, similarly to what is observed in the C4 simulations. On the contrary, fewer amount of particle agitation is produced as the particle frictional effects become more significant. These results show that, for low friction between particles, the fluidized bed may be considerably heterogeneous when friction at wall is high and become homogeneous when friction between particles is strong. It has to be noted that, in all these simulations, the computed particle agitation energy in the upper part of the bed (z/H ≥ 0.375), has a very low significance since, in this part, the particle volume fraction is nearly zero. Figure 16. Time-averaged dimensionless particle agitation energy, q 2 p /U 2 f for different values of the particle-wall µ w and the particle-particle µ p friction coefficients. Upper panel: constant µ p and different µ w . Bottom panel: constant µ w and different µ p . Re p = 70 and bed weight 75 g. Figure 17 shows the time-averaged field of the solid mass flux for the two sets of simulations. From the visualizations, it appears that friction affects significantly the magnitude of the solid flux and the extent of the mixing loops. In some cases, it also affects their position from the bottom of the bed. As previously observed for the particle volume fraction distribution, the position of the centre of the double mixing loop moves upwards when µ w increases (upper pannel). This point is not straightforward when µ p is increased (bottom panel). However, qualitatively speaking, the magnitude of the solid flux exhibits a maximum for µ p = 0.2 in both sets of simulations. Figure 17. Time-averaged solid mass flux, < α p ρ p U p > xz , for different values of the particle-wall µ w and the particle-particle µ p friction coefficients. Upper panel: constant µ p and different µ w . Bottom panel: constant µ w and different µ p . Re p = 70 and bed weight 75 g.
A quantitative analysis of the solid flux may be performed via its vertical profiles at a given height of the bed. Profiles of the time-averaged normalized vertical solid mass flux are depicted in Figure 18a,b for the two sets of simulations. The profiles are taken at the height z/H = 0.092 (z = 2.3 cm) above the bottom of the bed, for which experimental data are available. In all the simulations, two distinct parts for each profile can be observed. The first part corresponds to a nearly flat profile located in the region spanning the center of the bed and the second one is a sloped profile starting from the side walls. The expansion of one part to the detriment of the other depends on the values of the friction coefficients µ w and µ p . These coefficients also determine the maximum values of the upward and the downward mass fluxes in the flat and the sloped parts, respectively. The figures show that when µ w or µ p increases up to 0.2, the downwarding flux increases and the upwarding flux, close to the center of the bed, decreases. This leads, at the same time, to a steeper slope of the solid mass flux profile in the sloped part and to a less broad flat part at the center of the bed. For a stronger friction, either at wall or between particles, larger sloped parts are still observed, but with less steep slopes, leading to the decrease of both the downwarding and the upwarding solid fluxes. Therefore, despite the differences observed concerning the volume fraction and the agitation energy of the particulate phase, simulations C4 and C5 predict similar characteristics (field and magnitude) of the solid mass flux. It is conjectured that, in this confined configuration, when the wall friction coefficient increases, the front and back walls may considerably contribute together with the side walls to propagate stronger wall frictional effects towards the bulk bed. This has already been shown in a previous CFD/DEM study of Li et al. [33] for bed thicknesses of 10 and 20 particle diameters in a bubbling fluidized bed. As a consequence, stronger wall friction may produce almost the same effect on the solid flux as that induced by the friction between particles.
Finally, comparison with the experimental measurements shows that, the case with µ p = µ w = 0.1, is particularly interesting as it allows to reproduce the effect of the bottom wall on the axial solid flux ( Figure 18) and to better predict the location of the solid mixing loops from the bottom of the bed (Figure 17 vs. experimental data from Figure 12). Nevertheless, the solid concentration close to the side walls seems to be slightly underestimated ( Figure 14 vs. experimental data from Figure 12) and the solid loops more elongated than that experimentally reported. However, the predicted value of 0.1 for both the particle-particle and the particle-wall friction coefficients is noteworthy. In fact, the friction coefficient is a physical parameter that depends on the mechanical properties and the surface morphology of the solid particles, mainly Young's modulus and the surface roughness. Its experimental measurement is also very sensitive to some environmental factors, such as wet or dry conditions and the exerted normal stress. For this reason, a wide range of values are reported in the literature for glass beads under various conditions. Values ranging from 0.139 to 0.464 have been determined by Ishibashi et al. [34] for 1.6 mm-sized particles under dry conditions and various normal forces. Sandeep and Senetakis [35] reported a unique value of 0.2 for glass particles of 2 mm diameter. Additionally, other values may also be found in the literature as this is previously discussed in Section 3.2. The predicted value of 0.1 in our simulations is deemed to be acceptable and it is in very good agreement with that used by Goldschmidt et al. [29] in their simulation of glass particles in a pseudo-2D dense fluidized bed. Figure 18. Time-averaged vertical solid mass flux normalized by the gas mass flow rate, < α p ρ p U p ,3 > xz /ρU f , for different values of the particle-wall and particle-particle friction coefficients, and comparison with experiments. Re p = 70 and bed weight 75 g.

Application of the Results to Other Physical Configurations
In order to generalize the previous results to different physical configurations, other Reynolds number and bed weight conditions are investigated using the intermediate mesh. The same value of µ w and µ p and equal to 0.1 is used in the simulations since, as shown in Section 4.2, this leads to minor deviations from the experiments conducted for Re p = 70. A gas no-slip condition at the wall is also retained. In the simulation C6, a superficial gas velocity of 1.71 m/s corresponding to Re p = 100, as reported in Table 3, is imposed at the bottom of the bed. Predictions obtained by this simulation concerning the time-averaged fields of volume fraction and mass flux of the solid phase are reported in Figure 19a,b, respectively. Experimental measurements are also included for comparison purpose. One can see that the simulation reproduces very well the experimental solid distribution and the solid mass flux field. However, slight mismatches can be observed concerning the mean bed height and the position of the solid mixing loops. At the height z/H = 0.092, numerical and experimental vertical mass fluxes are displayed on Figure 20. Globally, results are conveniently reproduced by the simulation and confirm the accuracy of the selection µ w = µ p = 0.1.
(a) < α p > xz (b) < α p ρ p U p > xz Figure 19. Numerical (a) against experimental (b) time-averaged fields of the volume fraction and the mass flux of the solid phase. Re p = 100 and bed weight 75 g. In the simulation, µ w = µ p = 0.1. Figure 20. Numerical against experimental time-averaged vertical mass flux, < α p ρ p U p ,3 > xz /ρU f , of the solid at the height z/H = 0.092 above the bottom of the bed. Re p = 100 and bed weight 75 g. In the simulation, µ w = µ p = 0.1.
Results of a fluidization experiment of 125 g of particles by nitogen injected at a velocity of 1.54 m/s, corresponding to Re p = 90, have also been reported in Patil et al. [27]. This operating point is simulated in the present work (simulation C7) and results displayed in Figure 21a,b, together with the experimental data. Based on the investigations reported in Section 4.2, a value of the friction coefficient that allows to match as best as possible the experiments is retained. However, for all the simulations investigated so far, common deviations from the experiments are observed, namely concerning the mean bed height that appears overestimated, and the lateral solid recirculation zones which seem to be more elongated than that experimentally reported. The present simulation, using a higher bed weight, also does not sufficiently match the experiments, although the general pattern of the bed is globally well reproduced. Improvement of the numerical results should be achieved by, first, selecting an accurate drag law, which would allow to decrease the mean bed height, and then, by refining the mesh, which, as shown in Section 4.1, allows to accurately reproduce the solid concentration near the side walls.
(a) < α p > xz (b) < α p ρ p U p > xz Figure 21. Numerical (a) againt experimental (b) time-averaged fields of the volume fraction and the mass flux of the solid phase. Re p = 90 and bed weight 125 g. In the simulation, µ w = µ p = 0.1.

Conclusions
In this work, a DEM/LES approach is used to simulate a confined fluidized bed configuration for several physical and numerical conditions. Results are compared with experimental data of [27]. This study shows that, compared with the experimental measurements, predictions regarding the volume fraction and the mass flux of the particulate phase, obtained from the fine grid, are slightly better than that of the coarse and the intermediate grids. This slight improvement is attributed to a better prediction of the gas-particle interactions through the drag and the pressure gradient forces. The present work is also dedicated to frictional effects between particles and between particles and walls on the bed hydrodynamics behavior. It is shown that, increasing either the coefficient of the inter-particle friction or that of the particle-wall friction, the average bed height decreases and the bubble formation is enhanced. At low friction between particles, increasing the wall friction coefficient leads to similar solid circulation patterns (fields and magnitudes) as in the simulations with increasing friction between particles. Thus, it is conjectured that, in such confined 2D-configuration, when the wall friction coefficient increases, the friction of the front and back walls may have a significant impact on the inner-flow hydrodynamics of the fluidized bed. This finally may lead to the same bed hydrodynamics behavior as that observed when the friction coefficient between particles increases. Nevertheless, this point has to be further investigated. Indeed, to clearly understand the contribution of the friction at the front and back walls, simulations with various bed thicknesses combined with different values of the friction coefficients have to be performed. Additional simulations, performed at higher Reynolds numbers and/or other bed weights, showed globally good agreement with the experiments, reproducing almost similar particle volume fraction and solid circulation patterns as the experiments. However, deviations from the experiments concerning the average height of the dense bed is still exhibited in all the simulations. This could be further enhanced by studying other modeling features, such as the drag law. Funding: Computations were performed on the IFPEN supercomputer ENER110 and on the supercomputer CURIE under the project GENCI A0032B07345. The ENER110 and GENCI supercomputers teams are gratefully acknowledged. This work is done in the frame of the MORE4LESS project funded by the French National Research Agency.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Influence of LES Model on the Numerical Results
In order to clearly assess the influence of the LES model on the numerical simulations we performed additional simulations with and without LES model. The physical and numerical parameters corresponds to case C6 in Table 3. Figure A1 shows that there almost no difference between the case using a LES model and the case with no model. The closure law for small scale structures have a weak impact on the flow which means that the dissipation is mainly governed by the large scale structure and by particle-particle friction.