Multiscale Study of Reactive Dense Fluidized Bed for FCC Regenerator.

la particules catalyseurs e´coulements re´action global particule Abstract — Multiscale Study of Reactive Dense Fluidized Bed for FCC Regenerator — This study deals with reactive gas particle ﬂows like the coke combustion during the regeneration of FCC particles. In this kind of reactive ﬂow, the global reaction rate is usually bad predicted. In a ﬁrst approximation, the chemical scheme can be the reason because of the limitation of its modeling. It is usually based on macroscopic experimental results. The link between these macroscopic measurements and a local kinetics of the heterogeneous reaction occuring at the gas-particle interface is not conﬁrmed. Results of kinetics coming from experimental measurements are used and we try to highlight the problems that appear when the same kinetics are used at different of the particle volume fraction on the reaction rate and the necessity to develop a subgrid model for the reaction rate due to the natural clustering that appears in gas-particle ﬂows and non-linear additional terms appearing in ﬁltered equations.


INTRODUCTION
The oil cracking is an industrial process for the conversion of heavy molecules of hydrocarbons into lighter molecules for the use of fuel in the transport industry. Basically the oil cracking is performed in a riser with selective catalytic particles (FCC) and high temperature (400-500°C). During the reaction a thin layer of coke appears at the surface of the FCC particles. This thin coke layer inhibits the particle catalytic properties. In order to recover this property, particles are inserted into a regenerator where the presence of oxygen leads to removal of the deposited coke by a gas-solid combustion. The good mixing properties of a fluidized bed led engineers to use such a process for the regeneration of FCC particles.
The numerical simulation of the hydrodynamics of FCC regenerator is challenging because FCC particles are A-type in the Geldart's classification [1]. Agrawal [2] has shown that the main problem lies in the bad prediction of the small solid structures, smaller than the mesh size. These heterogeneities induced by the presence of clusters are not accurately solved when a too coarse grid is used in multi-fluid computations. When these small structures are neglected, it has been observed that the height of the bed is overestimated. Recent studies [3,4] proposed subgrid models as a correction of the drag force term when a coarse grid is used. A priori and a posteriori tests have shown that such models allow to predict the correct bed height.
Just as the drag force needs a subgrid model, the macroscopic gas production rate due to the heterogeneous reaction has to be modified as well. The local reaction rate, depending on both oxidant partial pressure and coke mass, is usually applied to estimate the total reaction rate in a mesh cell but it is based on the mean value of each component. These mean values are not representative of the real values at the position of the particle and can lead to errors in the prediction of the reaction rate.
In most FCC regenerators, air is injected through a bed of catalytic particles at a very low velocity so that the particulate flow is a dense fluized bed [5]. The gas phase finds a preferential way through the bed and usually at the center of the reactor where the walls do not have a significative impact. The consequence is that cavities of gas containing few particles appear in the bulk of the bed. The solid volume fraction around the bubbles ranges between 40 and 60% and is usually called the emulsion phase. Small bubbles appear in the lower part of the fluidized bed and they grew up by coalescence during ascension through the bed. Eventually, gas bubbles vanish above the bed when reaching the gas phase. It is also point out that the hydrodynamic properties of both dense and circulating fluidized beds are alike. The larger values of solid concentration are located near the walls of the reactor. Most of the gas injected through the fluidization grid goes up in the centre of the reactor. This effect leads to upward mean vertical velocities in the centre of the reactor for each phase and to negative velocities near the wall.
Specific studies have been performed previously by Bi et al. (2000) [6], Gobin et al. (2003) [7] or Fede et al. (2009) [8] about the hydrodynamics of gas-solid fluidized beds. Up to now very few works can be found about the kinetic scheme of the reaction in the gas phase or between two different phases in multi-fluid simulations. Eaton et al. (1999) [9] detailed the numerical approach for the study of reactive multiphase flows, especially for coal pyrolysis. Konan et al. [10] performed numerical simulations of a reactive fluidized bed with a N-phase CFD code in the framework of uranium fluorination. In 2008, Cao et al. [11] performed 2D numerical simulation of a FCC regenerator taking into account mass transfers between particle and gas due to the heterogeneous reaction and thermal transfers. The validity of numerical results is also an issue in reactive fluidized bed studies. The only data base available on the species mass fractions are usually collected at the outlet of the reactors. Local measurements inside the reactor are difficult and expensive and are often not really representative of the global flow.
The present paper reports the multiscale approach of a reactive fluidized bed for the simulation of FCC regenerators. Two upscaling effects will be analysed: from the particle to a cloud of particles, from a cloud of particles to the filtering size of resolution of continuity equations when using a multifluid approach.
Practically, the first upscaling is performed resolving Navier-Stokes equations of an air flow through a fixed array of particles. This allows to estimate the influence of some parameters on the gas mixture properties such as the reaction rate and the gas velocity field in the real conditions of FCC regeneration. The main objective is to develop volume source terms for the heterogeneous reaction based on a parametric study (impact of particle volume fraction and temperature, gas mixture). Only the impact of the particle volume fraction is highlighted here. The second upscaling is based on the computation of a low scale reactor which allows comparisons with experimental measurements. A 0D modeling permits to highlight the main issues. Multi-fluid numerical simulations are then performed in order to show the impact of the mesh on the reaction rate. Numerical results are also compared to experimental measurements.

DIRECT NUMERICAL SIMUMATION OF A REACTIVE GAS FLOW AROUND FIXED PARTICLES ARRAY
In order to understand accurately what happens localy in the flow during gas-solid reaction, direct numerical simulations are performed at the scale of a cloud of particles. The direct numerical simulation of a reactive gas flow going through fixed particle arrays deals with Euler/Lagrange type simulations (Euler for the gas phase and Lagrange for the particle). The particles are described by their real geometry (this study only takes into account spherical particles) and the fluid flow is exactly resolved around the particles. Direct numerical simulation is used because the order-of-magnitude of the size of the array is lower than the smallest scales of the studied turbulent flow. Therefore unsteady Navier-Stokes equations are solved with an explicit and compressible CFD code. Computations are performed using Lax-Wendroff numerical scheme which is of second order of both time and space. The numerical simulation ends when every statistical values of the conservative variables converge in time (then the reference values for the mean gas mass flux is reached).

Geometry and Meshing
Particles are arranged in a cubic centered array at each corner and at the centre of a cube.
Four different volume fractions (corresponding to four different volume fractions encountered in the dense bed) are studied. The higher dilutions (1% and 5%) correspond to the concentrations of particles in the bubbles of gas flowing through the bed whereas the smaller dilution corresponds to the emulsion. The particle diameter d p ¼ 70 lm corresponds approximately to the mean particle diameter of FCC particle size distribution. Eventually, the characteristic size L corresponding to the size of the edge of the pattern of the network can be deduced (Tab. 1), depending on the particle volume fraction. Table 1 also reports the total volume of the domain (including particles) and the volume of the gas. The order-of-magnitude of the smallest cells is about 1 lm.
As the numerical scheme is stable for a CFL number (1) equal to 0:5, the resulting time step is about Át ¼ 10 À9 s: where c is the sound velocity for an ideal gas, depending only on the gas temperature. The next step consists in determining the yielding velocity of the gas going through the cloud of particles. To get an order-of-magnitude of this velocity, the method used consists in considering an equilibrium state of the homogeneous fluidization for each volume fraction taking into account three volume forces: the pressure gradient, the drag force and the gravity. Figure 1 shows the values of the relative mean velocity in term of the particle volume fraction depending on the pressure P and the temperature T .  Values of the gas-particle relative velocity in a fluidized bed for two operating conditions in term of the particle volume fraction.
The resulting particulate Reynolds number (2) is extremely small so that the flow can be consider as laminar: where U r is the mean gas-particle relative velocity and m g is the gas kinematic viscosity. Figure 2 shows a global overview of the meshing elementary pattern for a p ¼ 0:05. Figure 3 shows different slices of this grid. It contains 126 740 tetraedrons. The design of the mesh requiers particular vigilance on the number of cells meshing the surface of each particle in order to take into account the heterogeneities of the local mass transfer going through the gas-particle interface due to the heterogeneous reaction. But the size of the smallest cell must be large enough in order to obtain a time step as long as possible because the convergence of the reactive case needs a long physical time. A good prediction of the mass transfer between the particle and the gas leads to a good prediction of the species mass fraction around the particle. That is why the thickness of the first cell near the wall must be small enough. This is the most restrictive criteria in the design of this mesh because the Reynolds number is very small. Table 2 gives the mesh properties for each considered particle volume fraction. The characteristic mesh parameter is defined as the ratio between the particle perimeter and the reference lengh of the smallest cell located at the gas-particle interface (f ¼ pd p =V 1=3 cell ). The retained mesh parameter criteria in order to solve the reactive problem varies in a range of 40 to 44 cells per particle perimeter which leads to an accurate prediction of the reactive gas flow.
x z y Figure 2 Global overview of the tetraedric meshing of the cubic centered array of particles for a p ¼ 5%. Slices of the tetraedric meshing of the cubic centered array of particles for a p ¼ 5% in the planes of normal vectors (0,0,1) and (0,1,1).

Numerical Simulation Organization
The main strategy for the direct numerical simulations is reminded hereafter. The numerical simulation must represent an air flow going through an infinite network of particles in the x-direction. The particularity of this computation is that upstream from x ¼ 0, the combustion is frozen (no coke deposition is considered for x < 0) whereas downstream from this x-coordinate, particles are loaded with 1% (mass fraction) of coke. The upstream part (cold flow) is solved performing a periodic numerical simulation in each direction setting the gas mean mass flux in the x-direction.
Periodic computations can be performed implementing dynamic source terms [12] in the momentum equation projected on the x-direction (S 0;x ) and in the energie equation (S 0;E ). They are evaluated with an iterative process, solution of Equations (3) and (4): with: where q ref and u x ref are the reference values of the density and the velocity respectively. These reference values have to be reached when the computation converges. Figure 4 shows the gas axial velocity field resulting from the periodic computation inside the cubic centered network for two different particle volume fractions (5% and 30%). A slice of the elementary pattern has been performed in order to show the isolines of velocity between the central sphere and the sphere located at the corners. This velocity is here expressed as a percentage of the reference velocity: u y ¼ u=u ref Â 100. As the particulate Reynolds number is very small, an attached flow is obtained. The gas flow fixes to the sphere after the particle. The presence of lateral sphere induces a containment effect. That is why the gas velocity increases between two spheres. This phenomena also delays the critical value of the Reynolds number which allows the gas to detach the particle. In the streamwise direction (x-direction), the interaction with both upstream and downstream particle leads to stagnation zones.
Then, the gas flow is reactive for x > 0. The coke deposited on the particles reacts with the dioxygen following the reaction: Due to the non periocicity of the reactive computations, boundary conditions have to be set for the downstream part of the network. Thus, at the inlet of the particle array (x = 0), the velocity field is deduced from the periodic cold computation. Uniform species mass fraction and temperature are set to constant values ( Fig. 5). At the outlet, a constant pressure is assumed.
The resolution of the periodic simulations allows to perform the reactive cases. The velocity profile obtained  Normalized velocity field for a) a p ¼ 0:05 and b) a p ¼ 0:30.
during the periodic simulation is set as a constant and non-reflecting boundary condition [13] at the inlet of the domain for the reactive computations. The pressure P ¼ 2 bar is imposed at the outlet of the network corresponding to the pressure of the gas in FCC regenerators. The heterogeneous reactions are taken into account using a dedicated boundary condition where the normal mass flux at the gas-particle interface is computed using the previous kinetics [14]. In order to obtain a stationary solution of the problem, the carbon quantity transported by the particle is assumed to remain constant and the coke mass fraction is chosen equal to 1%. In the gas phase, the carbon monoxide also reacts with the dioxygen and leads to an undirect production of carbon dioxide: These two reactions are supposed to reproduce accurately enough the real kinetic scheme. Particularly every reactions involving water (H 2 O) and dihygrogen (H 2 ) in the gas phase are neglected. They are assumed to be present in the reactor in a so small concentration to have a negligible impact on the global burning rate.
The kinetic law of the heterogeneous reaction is provided by measurements on a fixed bed and is given by: where A 1 is a pre-exponential factor, P O 2 the dioxygen partial pressure, n 1 the exponent law, B 1 the activation temperature and T the temperature. The resulting local mass transfer can then be written: where m C is the total catalyst coke content and S p its total surface. The molar ratio / is also depending on the dioxygen partial pressure: The retained values of each constant parameter are given in Table 3.
Therefore, the heterogeneous reaction is only described as a local mass flow going through the fixed particle interface taking into account reaction occuring inside the particle due to the particle porosity and depending on both values of temperature and dioxygen mass fraction at the particle interface.
The physical properties of the gas going through the cloud of particles are provided by the operating conditions encountered in FCC regenerators. The air is injected at the inlet at T ¼ 1 000 K and P ¼ 2 bar. All other properties are reported in Table 4. It is also assumed that the temperature at the surface of the particle is uniformly equal to T ¼ 1 000 K. The low variation of temperature inside the reactor (50 K) leads to low errors on the reaction rates and allows to make this assumption. Therefore, the resolution of the energy equation inside the particle is not performed. Some properties are imposed as an inlet boundary condition Sketch of the configuration for the DNS.   Table 5 summarized the numerical values of the velocity obtained for each particle volume fraction.
Ideally, it is appropriate to study the reactive flow along a large number of elementary patterns so that a position x can exist inside the computational domain where all the dioxygen has reacted. Basically, it depends on the mean velocity of the gas. Indeed, the faster the domain is supplied in dioxygen, the higher is the position of dioxygen total consumption. Then, it depends on the compactness of the network because the lower is the compactness, the higher is the distance the gas needs to cover in order to meet the same coke quantity. Eventually, the lower the compactness, the higher is the gasparticle relative velocity.
Because of these two additional effects it should be appropriate to perform numerical simulations through very higher distance for low compactnesses than for dense cases. Here, in order to compare the numerical simulations of each particle volume fraction, the same number of particles is imposed in each case. This number is arbitrary chosen to 12. It represents six elementary patterns of the cubic centered network.
The stationnary solution of the reactive case for each dilution is computed. The computation is performed when a stationnary state is reached. Figure 6 shows the species mass fraction field on a slice of the cubic centered network (plane z = 0). Dioxygen mass fraction is injected with a value Y O 2 ¼ 0:23 at the inlet of the network, corresponding to the maximum value encountered in industrial FCC regenerators. The computation provides informations on the consumption rate of dioxygen in the most dilute case. In the spanwise direction, concentration gradients are very low. Because the temperature is not high, the combustion rate is not significant compared to the convection effect. Moreover, this computation shows that, even for very dilute particle volume fraction, the isolines of species mass fractions near the particles are not like in several studies undertaken before by Xu and Fu (1997) [15], Blake (2002) [16] or Wichman and McMaster (2007) [17] concerning a single reactive particle. Due to particle-particle interactions, the isolines of species mass fractions are significantly different.

Results and Discussion
Finally, downstrean from the spheres the dioxygen mass fraction is lower than upstream because the interaction with the following particle decreases the convection effect. (Fig. 6 isoline Y O 2 ¼ 0:2279). As the dioxygen supplied by convection at this position is lower, the reaction has a more significative effect.
At the outlet of the network (x ¼ 6L a p ), the dioxygen mass fraction reaches a value about 0:2254.
Quantitatively, the production of carbon dioxide is very low for a p ¼ 0:01. The maximum mass fraction reaches 37:4 Â 10 À4 at the outlet of the network. In addition, carbon monoxide production is very low and even lower (24:1 Â 10 À4 ) than the carbon dioxide production. For this operating point, two reasons lead to the low production of carbon in gas phase: the high dilution and the high gas-particle relative velocity (also due to the high dilution).
Moreover, we can also point out that the convective boundary layers around all the particles lead to a species boundary layer (isolines In the emulsion, the solid mass fraction is about 50%. This case is interesting because the main part of the dense fluidized bed is occupied by this phase. Such a compacity leads to a relative gas-particle mean velocity equal to 0:67 Â 10 À2 m=s. Figure 7 shows the species mass fraction fields resulting from the computation. The dioxygen almost disappears at the position x ¼ 540 Â 10 À6 m. This simulation highlights the reaction thickness between the bubbles and the emulsion phases. Carbon species are here released in high quantities because all the dioxygen is converted. The carbon dioxide reaches the value of Y CO 2 ¼ 17 Â 10 À2 at the outlet of the domain and the carbon dioxide reaches the value of Y CO ¼ 11 Â 10 À2 . The mass production  cannot be neglected in this configuration because it has an impact on the velocity field around the particles. Finally, the direct numerical simulations of the reactive flow show that there is a high correlation between the particle volume fraction and the value of dioxygen mass fraction (Tab. 6). However, the reaction rate at the surface of the particle is not significant and no impact of the resulting surface mass flux on the velocity field appears. Direct numerical simulations also highlight that, in such a reactive gas-particle fluidized bed, the reaction mainly occurs at the interface between the emulsion and the bubbles of gas flowing through the bed. The upscaling from the particle to a cloud of particles has been addressed partially. However, the parametric study using direct numerical simulation has to continue. Especially, the impact of both gas and particle temperatures as well as coke mass quantity deposited on the particle must be highlighted. The complete study of the local particle combustion is necessary in order to find accurate correlations for the gas-solid mass transfer for FCC regenerators multifluid computations.
The next part of the paper will focus on the upscaling from the cloud of particles to the filtering size in Euler/ Euler computations, based on low-scale reactor numerical simulations.

EXPERIMENTAL SET-UP
A first order validation of a macroscopic numerical simulation usually requires the knowledge of experimental measurements on the macroscopic characteristics of the flow. For gas particle dense fluidized beds, the mean height of the bed is the first criteria for the validation of the numerical simulation. If one wants to take into account the chemical reactions occuring in the bed, the knowledge of the evolution of species concentrations at the outlet of the reactor can be precious. Experiments were performed at IFP Energies nouvelles on the lowscale reactor shown by Figure 8. The particle size distribution is given by Figure 9. The coke deposition at the surface of the particle reaches Y C ¼ 0:871%. Y C is the mass ratio between coke and FCC particles. In all this study, the density of the solid phase is supposed to be equal to q p ¼ 1 400 kg=m 3 even if one knows that the coke deposition usually increases FCC particles density due to the deposition in the micropores. Measurements are undertaken in two steps. The first one consists in injecting nitrogen (N 2 ) as a fluidization gas. When the fluidization is established, air is incorporated to the nitrogen. Then, a depleted oxygen air Y O 2 ¼ 9:2% is obtained and reacts slowly with the coke following reaction (6). The resulting fluidization velocity reaches V f ¼ 0:041 m=s. The gas pressure P g is equal to 1:76 bar and the temperature is regulated to about 945 K. The experiments are batch for the particle. Therefore the coming study is totally transient. Figure 10 shows the species mass fractions evolution measured at the outlet of the pilot.

0D MODELING OF THE REACTIVE FLUIDIZED BED
A simplified sketch of the reactor is given by Figure 11. The fluidized bed is supposed to be a perfectly mixed    reactor in the range z 2 ½0; h. Hence, the mass balance in the domain leads to: where q g and U g are respectively the gas density and velocity, a p the particle volume fraction and C g the gas-solid mass transfer. V f is the fluidization velocity and q inj the injected gas density. The unsteady momentum balance equation is not solved here but an equilibrium state is assumed leading to the balance Equation (12): The mass fraction balance equation of the species k (k ¼ O 2 , CO, CO 2 ), Y k , is given by: where / k is the source term of species k due to gas-solid and homogeneous combustions: Then, the mass fraction of N 2 is found using the definition Lastly, the catalyst coke content is given by Y C ¼ m c =m p and its evolution follows the balance equation: where m p is the total mass of solid in the bed. Practically, the total mass of the particles and the bed height of the homogeneous reactor are imposed and the temporal evolution of all mass fractions is computed.
The 0D modeling of such a reactor gives several informations. Does the prediction of the high of the bed have an effect on the species mass fractions evolution? Is it possible to have a good prediction of these evolutions knowing the height of the bed? Figure 12 shows the temporal evolution of the dioxygen mass fraction in the dense bed provided by a 0D modeling and by measurements at the outlet of the reac-  Comparison between dioxygen mass fraction evolutions provided by the 0D modeling for different values of h and experimentally. of U g found with the 0D computation, the time the gas needs to reach the upper part of the reactor is not significant. It means that the 0D modeling predicts the dioxygen appearance in the bed before it is obtained experimentally. The critical coke mass fraction on the particles is reached faster using the modeling. The reason is that the global reaction rate in the fluidized bed must be overestimated. The particle volume fraction affects also the characteristic time t c;b . The higher the fluidized bed is, the less is the particle volume fraction and the longer is the characteristic time. The particle concentration effect increases the reaction rate significatively. Then one can suppose that the physical fluctuations of the particle volume fraction in a dense fluidized bed due to the gas cavities going throw the bed (which are not taken into account with the 0D modeling) affect the global reaction rate. This shows the limitation of the 0D modeling in order to predict the evolution of the species mass fractions in the bed and anticipate the residence time of the particle in such a reactor to totally delete the coke. Figure 13 shows the evolutions of carbon dioxide obtained with the 0D modeling and experimentally. This figure also shows that the variations in carbon dioxide mass fraction are different numerically (with a bed height imposed to h ¼ 0:737 m) and experimentally (even if the same behaviours of the curves are observed).
During the simulation, this mass fraction takes quasiinstantaneously the value 6.6% at t ¼ 0 and then increases following an exponentiel law in order to reach a maximum value of about 10% at t ¼ 3 500 s. Then, it decreases at 0.1% in about 300 seconds. One can then suppose that all the coke has burnt at the physical time t b 2 3 500; 3 850 ½ . Figure 14 comes to validate this theory. It shows the catalyst coke content evolution during the recovering process.
Only the numerical data are available for this parameter because of the difficulties to obtain it experimentally. This graph shows that the mass fraction becomes actually lower than (1 Â 10 À4 % for t ¼ 3 800 s) for the same bed height h ¼ 0:737 m.
Experimentally, the production of carbon dioxide at the beginning of the air injection occurs slower than one can observe numerically. The carbon dioxide mass fraction increases linearly from 0 to about 8% in 400 seconds. Then one observes an inflexion point and an other linear evolution with a weaker slope. The carbon dioxide reaches the maximum value of about 10:7% at t ¼ 4 200 s. Finally, when all the coke has reacted, the carbon dioxide mass fraction decreases progressively. During the decreasing phase the slope of the curve is flatter experimentally than numerically. Moreover, it does not decrease immediately to Y CO 2 ¼ 0. When the value reaches 1:7%, a second inflexion point is observed and the slope becomes very small so that Y CO 2 still is equal to 0:6% at t ¼ 7 500 s.
This second inflexion point is not reproduced numerically. It can be due to the limitation of the kinetic law, considering that the coke stored in the pores of the particles reacts following the same law as the coke deposited at the surface. Actually, the production of carbon dioxide after t ¼ 5 000 s might be due to the gazification of the coke from the pores. Figure 15 shows also the evolution of the total carbon monoxide. Comparison between catalyst coke content evolutions provided by the 0D modeling for different values of h.
In order to improve the numerical prediction of the low scale reactor combustion and to account for the gas bubbles formations through the emulsion phase, the last approach deals with 3D numerical simulations using an Euler/Euler approach.

Mass Balance Equation
A summary of the Eulerian description of a multiphase flow has been addressed by Enwald et al. (1996) [18]. In 2009, Konan et al. [10] describe additional terms and equations in order to take into account a reactive flow. This paper deals with a monosolid fluidized bed. Only two phases are considered. If k takes the value p for the solid phase or g for the gas phase, the mass balance equation is given by: where a k and q k denote the volume fraction and the density of the phase k, respectively. U k;i is the ith component of the mean velocity of the phase k. C k represents the mean mass transfer of the phase k with all other phases and satisfies: As the mass transfer occurs from the solid to the gas, we will focus on the positive term C g for the modeling.
For the single solid phase p, the diameter is suppposed to be constant and equal to the mean diameter of the size distribution. The reaction law occuring at the particle surface is very slow and takes also into account the coke presence inside the micropores of the particle. The diameter of the particle is not modified by the heterogeneous reaction. It should be more appropriate to compute the density variation of the particle phase during the coke combustion. In a first approximation, a constant particle density q p is assumed.

Momentum Balance Equation
The Eulerian formulation of the solid phase allows to write both particle and gas momentum equation in a generic way. Assuming that the gravity g i is the only external contribution of the volume forces involved in this configuration, it can be written as follow: The kinetic stress tensor R k;ij is usually modeled using a Boussinesq approximation. It is the summation of the Reynolds tensor and the molecular collisons for the gas phase or the particle collision for the solid. P g is the mean gas pressure. The mean particle interface velocitỹ U r due to the gas mass production during the reaction is usually assumed equal to the mean solid velocity because: -of the sphericity of the particle; -the partial pressure of reactives around the particle is assumed to be homogeneous. We are going to show below that this second assumption is not really true in most of cases (particularly for dilute particles).
Eventually, the gas-solid momentum transfer I k;i occurs through the drag law: where s gp is the particle relaxation time defined bellow: According to Gobin et al. [7], C D is evaluated by coupling Wen and Yu and Ergun's laws: Comparison between carbon monoxide mass fraction evolutions provided by the 0D modeling for different values of h and experimentally.
where V d;i denotes the drift velocity due to the turbulence transport of the solid phase by the fluid turbulence. Its modeling was suggested by Deutsch and Simonin (1991) [19] and takes the following expression: where D t gp is the fluid-particle dispersion coefficient, proportional to the fluid-particle fluctuating velocities covariance q gp and the characteristic time of the fluid turbulence seen by the particle s t gp and is defined by: with: where r a is taken equal to the Schimdt or the Prandtl number of the gas phase and C b is the coefficient of crossing trajectories equal to 1:8. n r is the ratio between the mean component of the gaz-particle relative velocity and its fluctuating component: The fluid-particle fluctuating velocities covariance is the solution of the differential equation (Simonin et al., 1993 [20]): The dissipation and the productin rates are given respectively by Equations (30) and (31): Finally, since g r ¼ s t gp =s F gp Fe´vrier and Simonin [21] write R gp;ij as:

Modeling of the Particle Kinetic Stress Tensor
A closure modeling of additionnal terms in the momentum equation for the solid phase is needed due to the contributions of fluctuating velocities and collisional effects. The sum of the Reynolds stress tensor and the collisional terms in dense flows is modeled by the following approach (describing a Boussinesq approximation): k p denotes the bulk viscosity: l p takes into account the stress viscosity and its formulation is the summation of the turbulent and the collisional viscosities l p ¼ a p q p m t p þ m c p . According to Simonin (1995) [22], the turbulent viscosity of the particulate phase is written: This turbulent contribution is also involved in the computation of the collisional viscosity: Eventually, the particle kinetic agitation balance allows to close the problem: The diffusion term needs a model for the total diffusivity coefficient of the dispersed phase agitation K t p . It is provided by the summation of a kinetic part K kin p and a collisional part K coll p : The production of agitation due to the gas-particle mass transfer is usually assumed to be: where the two collisional coefficients are written: Þð49 À 33e c Þ ð 43Þ and the characteristic time of particle-particle collision is: P q p represents the effect of the gas phase on the particle agitation: when e p is the particle agitation dissipation due to the collisions:

Species Mass Balance Equation
The numerical simulation of a reactive flow requires one additional equation for each species of the gas phase such as the mass fraction balance equation: Therefore, assuming the gas to be ideal, its density evolution is determined by: The mean molecular weight is deduced from:

Filtered Mass Balance Equation
The analysis of the results of continuity equations for Geldart A particles provided by different levels of grid refinement permits to highlight the strong grid dependency. Therefore, the subgrid terms can be determined using a filtered approach. Assuming that the solid volume fraction field obtained with a very fine mesh can be written a p ðx; tÞ, then the filtered value of this solid volume fraction computed with a coarse grid is: where Gðx; xÞ is the weighted function of the filtered value centered inx [23].
The filtered values of the density and gas velocity are then given by: a g ðx; tÞq g ðx; tÞŨ g;i ðx; tÞ ¼ Z V Gðx; xÞa g ðx; tÞq g ðx; tÞ The resulting filtered gas mass balance equation is then given by: Eventually, the non-linearity of the term C g leads to additional terms, not only in the mass balance but in every continuity equations. These terms need to be evaluated. Therefore, the next part of this paper deals with the resolution of non-filtered continuity equations with a coarse and a fine mesh.

3D NUMERICAL SIMULATION OF THE REACTOR USING A MULTI-FLUID CFD CODE
Numerical predictions of the low-scale reactor are undertaken using multi-fluid CFD code NEP-TUNE_CFDV1.08Tlse. NEPTUNE_CFD is a multiphase flow software developed within the framework of the NEPTUNE project which is a french industrial consortium, financially supported by CEA, EDF, IRSN and AREVA NP. A mono-solid gas-particle flow is assumed and for each phases, two different wall boundary conditions are used. Basically, for the gas phase we used a classical friction law dedicated to disturbed flows described with k À e model. For the solid, a no-slip boundary condition appears to be appropriate [8].

Geometry and Meshing
The computational domain must be meshed so that the numerical solution of the discretized equations converges to the real solution at each iteration. It is often found that the numerical simulation depends on the size of the mesh. The non-dependency of the mesh is reached for a very low characteristic size. That leads to a high number of cells to mesh the domain and to a very small time step. The high number of cells requires to perform the computation with a lot of processors in order to improve the real time of computation. The small time step leads to a higher number of iterations. Computations are usually performed on a coarser mesh using subgrid scale modeling of physical phenomena in order to improve the CPU time.
The two O-Grid meshes used in order to performed this computations are shown in Figure 16.

Numerical Simulation Organization
The first step of the computation consists in fluidizing the fixed initial bed of particles with an inert gas.
Initiation of the computational domain is performed inserting homogeneously a particle total mass of m p ¼ 16:5 kg along an arbitrary height of the reactor. We assume that the heterogeneous reaction does not change the particle mass and diameter. That is why only two phases are considered for the computation with only one solid phase of constant diameter d p ¼ 70 lm. This assumption is quite valid because the mass fraction of coke deposited on a particle is usally between 10 À3 and 1%. Both gas and solid velocities are initialized to 0 in all the domain and the gas phase is only composed with nitrogen, so that no reaction occurs at t ¼ 0. This homogeneous and arbitrary initiation does not match with any physical behaviour. A physical solution of the non reactive problem must be generated before starting to study the reactive case. The first step of the numerical simulation consists in fluidizing the bed with an inert gas (N 2 ). Then, the bed reaches a so-called statistically stationnary state when the time averaged values of each variable does not change in each point of the domain. Especially, the mean velocity components and mean volume fraction of each phase remain constant in time. Finally, the bed is considered as a perfectly-mixed reactor. Thus, a physical initial solution for the reactive computation is used. For each mesh, during the five first seconds, the unsteady flow is computed in order to obtain a perfectly-mixed reactor, then, during the following 20 seconds, time average values are computed.
For Geldart A particles, the height of the fluidized bed is usually not well predicted numerically except if a very fine mesh is used. But when a very fine mesh is used, the time step collapses and only allows to compute few physical seconds. Numerical simulations undertaken here are reactive and unsteady. Experimentally, using the law scale reactor, the total combustion of the coke deposited on the solid phase takes a couple of hours.  Therefore, it is not possible to use a non dependent mesh in order to follow the temporal evolution of species mass fraction in such a reactive fluidized bed. That is why a subgrid model for the drift velocity proposed by Parmentier [3] is used in order to have a good prediction of the drag law with such a mesh. Thus, performing numerical simulations for 2 different meshes, the height of the bed should remain the same and the impact of the mesh size on the prediction of the global reaction rate can be observed.
The subgrid model of the drag law consists in modifying the gas-particle relative velocity using a drift velocitỹ V d;i .
Then the apparent relative velocity is the sum of the resolved relative velocityW r;j and this drift velocity. This leads to:Ṽ The K ij values are constant. They are supposed to be equal to 0 if i 6 ¼ j. Thus, only the values of K xx and K yy must be determined. The function g is determined by Parmentier [3] resolving the same flow (of FCC particles) and performing a mesh dependency study. He obtained that g can be approach by: It is defined as a product of a function of the mesh parameter Á H and a function of the mean particle volume fraction a p .
Moreover if u ¼ a p =a m , where a m ¼ 0:64 is the maximum geometric particle volume fraction for a random compactness of a monodispersed particulate flow, h is defined as: Eventually, f can be written as: where a ¼ 6:13 Â 10 À2 and Á H is given by: where the characteristic size of the filtering Á is assumed to be dependent on the mesh size Á 2 ¼ 2Á 2 x . Thus Figure 17 shows the mean pressure profil in the reactor for the two different grids.
In each computation, the subgrid drift velocity is used in order to take into account the small structures. The pressure profils represent two straight lines separated by a break up of the slope around P À P ref ¼ 0 at the hight h ¼ 0:760 m. The mean height of a dense fluidized bed is typically characterized by this break up. Experimentally, this height is measured at h ¼ 0:737 m. Finally, the subgrid modeling for the drag law permits to predict the appropriate mean height of the fluidized bed with a mistake only equal to ð0:760 À 0:737Þ=0:737 ¼ 3:12%. At the same time, simulations are undertaken without using the subgrid modeling on the same meshes. In this configuration, the height of the bed is overpredicted and the bed even reaches the divergent part of the reactor (Fig. 18). Quantitatively, the instantaneous solid volume fraction is very low compared to the subgrid model. We can recall that in the emulsion phase the order-ofmagnitude of the particle volume fraction is about 0:50. Figure 15 highlights that the prediction of the volume fraction in the emulsion phase is well predicted using the subgrid model. Without the model the particle volume fraction is only about 0:30. Figure 19 shows a comparison between the temporal evolution of coke total mass in the reactor provided by the 0D approach and the 3D numerical simulations given by: Mean pressure profile in the fluidized bed for two different levels of grid refinement.

Results and Discussion
For each case, it follows a linear evolution and the total coke consumption occurs for t ¼ 4 100 s.
In the 3D numerical simulation, the volume mean value of the species mass fraction k is given by: where V bed is the mean volume of the bed taking into account the mean heigh of the bed h bed found to be the height where the pressure gradient changes. Figure 20 shows the evolution of the carbon dioxide mass fraction in the bed. The results of the 3D numerical simulation are different from the 0D prediction. In the first step (when the dioxygen injection starts), the evolutions are similar. However the break up of the slope does not match. It happens between 1 and 2 minutes when enough carbon gas has been produced and totally filled the reactor. The 0D prediction of the critical value of carbon dioxide is about Y CO 2;S ' 6:6% whereas the 3D numerical prediction is about Y CO 2;S ' 8%. When this value is reached, the carbon dioxide mass fraction continues to increase weakly exponentially. After this second part of the flow, the maximum value of carbon dioxide is reached. The 0D modeling predicts a maximum value Y CO 2;max about Y CO 2;max ' 10% when the 3D numerical prediction gives Y CO 2;max ' 9%.
The last part of the reaction starts when the mean value of the carbon dioxide collapses fastly. During this third part, the coke quantity transported by the solid becomes very scarce and does not always react with the dioxygen injected through the bed. The dioxygen starts to appear inside the bed. This phenomena occurs earlier with the 0D numerical simulation than with the 3D approach. Again, this is due to the different probability for a particle to cross a bubble of dioxygen between the 3D numerical simulation and the 0D approach. So the delay in the 3D numerical simulation is due to the fact that the coke solid combustion goes faster with the 0D approach.
Moreover, Figures 20, 21 and 22 show that the numerical simulations with both fine and coarse mesh lead to the same species mass fraction evolutions in the bed. The impact of small structures cannot be pointed out using these two meshes. In order to go deeper in the Evolution of carbon dioxide mass fraction. study, very accurate computations have to be performed with an extremely fine mesh.

CONCLUSION AND PROSPECTS
The multiscale study of a reactive flow has been performed in order to understand the reasons of a difference between experimental investigations and current numerical simulations. The low scale analysis consists in a direct numerical simulation of the gas phase through an array of reactive particles. These simulations highlight that the particle volume fraction of the cloud and the reactant species mass fraction are correlated. As the macroscopic reaction rate depends on both particle volume fraction and dioxygen species mass fraction, filtering gas mass balance equation leads to a correlation between these two parameters. In order to get an evaluation of this correlation, multi-fluid numerical simulations at the pilot scale have been performed for two different meshes. However, we highlight that the impact of the size of the grid is not significant. In order to point out a mesh dependency, we should perform numerical simulations with a very fine mesh but the computational time would be very long. In addition, numerical simulations of an industrial case should be better to highlight the heterogneneities in the bed because at the lab scale, the reactor is not wide enough to see bubbles of dioxygen going up through the bed, so that the heterogeneous reaction only occurs at the bottom of the bed. This effects does not allow to identify any mesh dependency. Evolution of carbon monoxide mass fraction.