est destinée au dépôt

grown polymer particles. To account for this PSD in a numerically trac table manner, we choose to describe particles of two characteristic sizes: (i) large particles (1,600 μm) which are representative of polyeth ylene particles in the reactor and will yield a realistic fluidized bed hy drodynamics and (ii) fine particles (80 μm) which are characteristic of injected catalyst particles. The growing aspect of particle sizes is not accounted for, as the simulated time is short compared to particles res idence time in the reactor. Therefore, the three present phases are the gaseous phase and the two particulate phases.

One hundred tons of large particles are present from the start of the simulation whilst fine particles are injected through the four catalyst in jectors. Both categories of particles are subjected to an exothermic reac tion (see Section 2).

The powder and gas material properties are gathered in Table 1. Note that the simulation does not account for the variation of these physical properties according to temperature and pressure. For the gas phase, the reference pressure is 24 bars.

Operating, boundary and initial conditions

We start with an homogeneous distribution of solid volume fraction in the bed part: α f 0 Initial gas and particle temperatures are set to 100 °C. Operating conditions are defined by taking into account the exothermic poly merization reaction. Thermal source terms are used to implement this exothermic reaction. Fine particles are 30 times more reactive than large ones per mass unit. Their reaction dissipated powers are Q f = 6.0kW/kg f and Q l = 0.2kW/kg l . Enthalpy source terms for each phase (see Eq. ( 4) and Fig. 2) are then given by:

S g 0:0 kW=kg g ð5Þ
S f α f ρ p Q f ð6Þ
S l α l ρ p Q l ð7Þ
For the simulation, the total mass of large particles will be M l = 100,000 kg l and the total mass of fine particles will be negligible for the time period simulated: about 5.2 kg f after 25 s. Therefore, the total reaction dissipated power will be about Q l M l = 20MW which is of the same order of magnitude than in an actual reactor.

The geometry has four kinds of boundary conditions summarized by Fig. 6 (right). The fluidization plate, through which air was injected at a constant mass flow rate, is modeled by a gas inlet with uniform superfi cial velocity equal to V g = 0.6m/s and set temperature. This boundary was seen as a wall by solid phases which induces a null flux for solid phase enthalpy. The detail of boundary conditions applied for each phase is given in Table 2.

The wall boundary condition was a no slip condition for particles and a friction condition for the gas. The no slip wall condition for parti cles correspond to a zero flux for random kinetic energy [START_REF] Fede | 3D numerical simulation of a lab-scale pressurized dense fluidized bed focussing on the effect of the particle-particle restitution coefficient and particle-wall boundary conditions[END_REF]. A Neumann type boundary condition is imposed for gas and solid phase enthalpies on walls.

The top of the reactor is a free outlet with set gauge pressure. Should a depression appear just below this surface, the outlet boundary condi tion acts as a gas inlet with a set temperature of 100 °C.

To take into account the polymerization reaction, we consider an exothermic reaction through heat source terms on particles. Enthalpy transfers between gas and large particles and between gas and fine par ticles are taken into account.

V g = 4U mf α inj, g = 0.99976 k = 1.0 × 10 4 m 2 ⋅ s 2 k = 1.0 × 10 4 m 2 ⋅ s 2 ε = 1.0 × 10 3 m 2 ⋅ s 3 ε = 1.0 × 10 3 m 2 ⋅ s 3
Free outlet for enthalpy. Zero flux of enthalpy. T g = 50 °C

T inj, g = 50 °C Fine part.

No-slip. Q inj f = 0.18 t/h/injector Free flow with imposed gauge pressure. α inj, f = 2.4 × 10 5 No-slip. q p 2 = 5.0 × 10 5 m 2 ⋅ s 2 q fp = 1.0 × 10

Mesh construction and mesh partitioning

In this study, the mesh refinement is a key point. We first detail how to generate a billion cells mesh for the industrial scale geometry, which is not feasible using a classical one step mesh generator. This requires more memory than available on most workstations or supercomputers compute nodes. Moreover, only few meshing software are parallelized. The use of supercomputers large memory nodes could be considered, but this approach would be prohibitively long due to the limit in com putational power of these nodes, and this approach would not scale to even bigger meshes. The use of a distributed approach is thus required. We then designed and applied a non conventional methodology using a multi step, multi software approach detailed hereafter.

To generate the large mesh used in following simulations, we relied on multiple useful features of Code Saturne [START_REF]Code Saturne Version 5.0.4 User Guide, EDF R&D; Fluid Dynamics, Power Generation and Environment Department Multi-Phase and Reactive Flow Group, 6 quai Watier[END_REF] in terms of mesh modification and extrusion, namely (i) its surface mesh extrude capability; (ii) its mesh merging feature (either matching or non matching) and (iii) its compatibility with user defined mesh transfor mations. The strategy has been the following:

1. Generate 4 mono layer O Grid type meshes representing hori zontal slices of the geometry: This step was performed using the software Simail [START_REF] Sas | SIMAIL version 8.0 manuel utilisateur, DISTENE, Campus Teratec, 2 rue de la Piquetterie, F-91680[END_REF]. Three of these meshes consist in 476,928 hexahedral cells each (see Fig. 3. Using user defined mesh transformations, first dilate the mesh radi ally to obtain the expected diameters, then generate the bulb shape (disengagement zone) by shifting 3D coordinates of all vertices in the three top cylinders. 4. Using the matching mesh merging feature of Code Saturne [START_REF]Code Saturne Version 5.0.4 User Guide, EDF R&D; Fluid Dynamics, Power Generation and Environment Department Multi-Phase and Reactive Flow Group, 6 quai Watier[END_REF], merge all 3D meshes into a final single mesh.

These 3D steps are summarized by Fig. 8. The final mesh contains 1,002,355,456 hexahedrons whose mean volume is about 0.125 cm 3 , i.e. 5 mm of characteristic size (see Fig. 9). The whole mesh connectivity is conformal. The refinement of this mesh is such that each injector of diameter d inj = 0.21m is defined by more than 350 cells. The mesh file, saved in a double precision binary format is about 209 GB in size.

All operations performed using Code Saturne [START_REF]Code Saturne Version 5.0.4 User Guide, EDF R&D; Fluid Dynamics, Power Generation and Environment Department Multi-Phase and Reactive Flow Group, 6 quai Watier[END_REF] benefited from the parallelization of the code and ran on CALMIP supercomputer (see Section 5.1 hereafter). This mesh creation step requires 90 nodes to have enough memory. The mesh construction step is critical for the sim ulation quality. It is a quick step in terms of CPU hours, note however that it was actually a long step in terms of dedicated human time: about 4 weeks.

Main characteristics of this mesh are detailed in Table 3. Quality criteria of NEPTUNE_CFD and Code Saturne are defined in their user's manual [START_REF]Code Saturne Version 5.0.4 User Guide, EDF R&D; Fluid Dynamics, Power Generation and Environment Department Multi-Phase and Reactive Flow Group, 6 quai Watier[END_REF]. The few cells marked as bad by respective criterion are located over the bed, in the top of the bulb. They will not impact the quality of the simulation.

To benefit from HPC capabilities of NEPTUNE_CFD [START_REF]NEPTUNE CFD Version 4.0.1 User Guide, EDF R&D; Fluid Dynamics, Power Generation and Environment Department Multi-Phase and Reactive Flow Group, 6 quai Watier, 78401[END_REF], partitioning of the mesh is required. Each CPU core will handle computations for one domain and will exchange information to CPU cores handling neighbors domains through MPI communications for ghost cell data. We look for a balanced partitioning (same number of cells in each domain) which also minimizes the number of neighbor domains for each domain, to limit inter domain communications to the bare minimum. The number of ghost cells (halo) is also a highly important metric of partitioning, and should be as small as possible. Other criteria to select a good partitioning algorithm or library are their stability (repeatable and stable) and the computation time required to perform this partitioning. Fig. 10 presents a decomposition in 64 domains of the mesh for two partionner, PT SCOTCH and Morton method. Each color corresponds to a domain solved by a CPU core.

Multiple partitioning libraries are available and have been assessed (see Table 4). Sequential partitioners are not able to handle such a large mesh. We then evaluated five parallel partitioners: PT SCOTCH [START_REF] Pellegrini | PT-Scotch version 6.0.4 user guide[END_REF], ParMETIS [START_REF] Karypis | ParMETIS Version 4.0.2 Manual[END_REF], Morton Curve, Hilbert Curve and Block methods. PT SCOTCH and ParMETIS make use of eponymous libraries, using de fault options, and are based on graph partitioning. Morton and Hilbert curve approaches are algorithm that use the corresponding space filling curves (see Fig. 11) to sort mesh cells. The obtained curve is then split into the required number of domains which induces perfect domain balance.

With this billion cell mesh, we wanted to create mesh partitioning from 1260 up to 36,000 cores, i.e. to create 1260 up to 36,000 domains.

ParMETIS (Parallel version of METIS version 4.0.2 [START_REF] Karypis | ParMETIS Version 4.0.2 Manual[END_REF]) was not able to handle such a large mesh over such a high number of cores, it was not stable and simply crashed and is thus discarded from the comparison with other methods. PT SCOTCH could generate a high number of domains, but could not run on many CPU cores either, otherwise the partionning step would fail. This induces a very slow domain partitio ning due to the low computational power of a limited number of cores. Nevertheless, this library generated the best results in terms of neighbor and ghost cells balancing. The domain unbalance was about 20%. Overall, this disequilibrium is compensated by the low number of ghost cells and induced more favorable computation times. Even though this library yields interesting results, the fact that it does not scale robustly for high domain counts (N4000) is very limiting. More over, this algorithm is based on a random seed and is thus non reproducible. We then saved the best partitioning among multiple calls to the library and re used it for multiple simulations.

The Block partitioning uses a dummy approach (simply partitioning by slabs based on initial ordering) which may be useful in some cases, but here, the domain unbalance was over 10%, the neighboring unbal ance was bad as well as the ghost cell unbalance.

Morton curve and Hilbert curve algorithms behave very similarly: stable, efficient, perfect balance for the chosen criteria (same number of local cells per domain ±1 cell). The Morton curve led to slightly better For large meshes, the space filling curve partitioning is thus pre ferred, with Morton being the best choice here. This is consistent with observations on many other large cases handled with Code Saturne.

NEPTUNE_CFD allows creating all needed mesh partitioning before simulation runs, using a pre processing mode. This is useful mainly to allow partitioning from a different number of ranks, and is often faster than recomputing the partitioning at each restart (depending on the performance of the parallel I/O subsystem). The different partitioning rank mappings are exported and can be re used according to CPU core number wished for simulation. Using this functionality we created all mesh partition maps before performing any simulation. Depending on the number of CPU cores for a set simulation, we read the adequate partition map during the simulation preprocessing step.

Results and discussion

Hardware and software resources

Numerical simulations of the present work have been made on the supercomputer Olympe at the Toulouse University Computing Center2 and on the supercomputer Gaia of EDF R&D. Main characteristics of these supercomputers are summa rized in Table 5.

The supercomputer Olympe is a BullSequana X1000, last genera tion of the European vendor ATOS aimed to Exascale computations. This cluster is made of 374 compute nodes interconnected with an infiniband fabric (EDR 100 GB/s) in a fat tree topology (blocking fac tor 2:1). Processors are Intel® "Skylake" 6140 (18 cores per socket@ 2.3 GHz), 13,392 cores as a whole in the cluster. All nodes have a bi socket architecture, i.e. two processors of 18 cores each. There are three kinds of compute nodes: 360 classical bi socket nodes, 12 nodes with GPU accelerator and 2 large memory nodes. All simula tions ran on classical bi socket partition, except the 362 nodes simu lations which used some GPU nodes but without using GPU accelerators. For this compute partition, nodes memory reaches 192 GB of RAM, for a total amount of 69 TBytes of distributed RAM. The overall peak performance is 1.365 Pflop/s. For storage a Lustre file system with 40 GB/s of bandwidth and 1.5 PB of space is con nected to nodes through an infiniband fabric. The supercomputer Gaia is a Atos Bull Cluster, also last generation of the European vendor ATOS aimed to Exascale computations. Gaia's ar chitecture is close to that of Olympe, but is three times bigger in terms of nodes and peak performances. Gaia overall peak performance is about 3.05 Pflop/s which ranks it 127th of TOP500 citeptop500 [START_REF]TOP500 Supercomputer Sites[END_REF] as of june 2019. The local storage is based on IBM Spectrum Scale (GPFS) and DDN SFA14KE with 200 GB/s of bandwidth and 16 PB of space. The parallel file system is connected to nodes through Intel OPA v1 and the total core number is 42,912.

Presented simulation benefited from the opportunity offered by pre production operations associated with the installation of Olympe and Gaia supercomputers. These installation phases allowed using these systems in an exclusive and dedicated manner. During the set up phase, strong interactions with CALMIP and EDF R&D Scientific Informa tion System architects and HPC engineers were necessary to circumvent all technical issues due to the fact that supercomputers were at this time newly installed. The size of the simulation made it an excellent bench mark for these two supercomputers and allowed fine tuning and vali dating them for IO, MPI, filesystems and node stability.

Major issues were two folds: Input/Output for one, and runtime for more than 10,000 processes for the other. It is noteworthy that memory footprint of the simulation was very low for such a scale, and was not an issue. Indeed, less than 10 TB of distributed RAM was necessary, com pare to the capacity of Olympe mentioned above (69 TB of distributed RAM). Due to the expected duration of the simulation, we needed to fre quently write checkpoint file to restart the simulation and to avoid loos ing CPU computation time in case of software or hardware failure. The size of the simulation lead to a check point restart file of 1.34 TBytes. On both supercomputers, libraries and software have been compiled using Intel Compiler 18.2 and relied on the MPI library IntelMPI 18.2. We used the capability of the software NEPTUNE_CFD [START_REF]NEPTUNE CFD Version 4.0.1 User Guide, EDF R&D; Fluid Dynamics, Power Generation and Environment Department Multi-Phase and Reactive Flow Group, 6 quai Watier, 78401[END_REF] to write and read in parallel (MPI IO) accordingly to the domain decomposition (partitioning). As we used the whole Olympe machine, this would have put a lot of pressure on the underlying file system (i.e. Lustre, character istics mentioned above). Nonetheless, Olympe and Gaia exhibited a great stability for reading and writing operations and with good perfor mances. For the runtime performance we experienced more difficulties. On Olympe, we encountered locks in communications with simulation above a threshold, namely over 5760 MPI processes (160 nodes).

To circumvent the problem we had to finetune the Intel MPI library that we linked against NEPTUNE_CFD to handle inter processes com munications. This tuning is quickly detailed in section 5.3. As we found the right parameterization on Olympe, we were able to run suc cessfully above 5760 MPI processes and reach the upper bound of 13,032 MPI processes (362 compute nodes). To sum up Olympe simula tions consumed 5 millions of core hours, and generated sets of data of 120 TB.

On Gaia, we have been able to reach the upper bound of 36,000 MPI processes (1000 nodes). Simulations which ran on Gaia consumed 10 millions of core hours and generated sets of data of 200 TB. Gaia benefited from previous experience on Olympe and no MPI issues were faced at large scales, though the same finetuning as been applied. Gaia being in installation phase, some issues related to OPA connections were also detected during NEPTUNE_CFD [START_REF]NEPTUNE CFD Version 4.0.1 User Guide, EDF R&D; Fluid Dynamics, Power Generation and Environment Department Multi-Phase and Reactive Flow Group, 6 quai Watier, 78401[END_REF] installation. Using OPA, for performance reasons, we choose to use unpopulated simulations using 35 cores out of 36 on each node.

Run organization and data management

Numerous successive simulations using checkpoint/restart were required to reach the total simulated time of 25 s due to some issues related to the scale of this simulation. These aspects are further discussed in Section 5.3.

A critical aspect that has been considered before running the simula tion is that of data export and postprocessing. On top of the raw chal lenge of running such large simulations, one goal of this study is to obtain a dataset characteristic of the functioning of an industrial fluid ized bed at the smallest scales currently achievable. The way these re sults may be used is discussed in Section 5.5 but generating these results is a whole challenge by itself.

In terms of data volume: NEPTUNE_CFD solves 22 PDEs. Exporting all variables in EnSight Gold Format (single precision) would take 200 GiB per time step. The simulation mean time step is about 2 × 10 -4 s, which corresponds to 125,000 timesteps for the whole 25 s simulation. It would then be impossible to save results either over the whole do main and/or for all timesteps. In the selected approach, 10 variables are exported at a rate of 40 exports per simulated second in 13 volumic and surfaces zones. Selected export areas are:

• 4 thick planes along the radial direction; • 7 thick planes at specific heights; • a whole cylindrical volume encompassing the catalyst injectors;

• the reactor external surface.

Despite reducing significantly the size of exports, the generated vol ume is still significant and very difficult to handle, manage and post process.

The Lustre file system is scalable, has high performances and allows parallel I/O operations. When using additional library for I/O (i.e. MPI IO in this case, but also HDF5 or parallel netCDF), reading and writing can be done in parallel from several nodes into single shared file. We used a technique called file striping to increase I/O performance. The Lustre file system is made up of an underlying set of I/O servers and disks called Object Storage Servers (OSSs) and Object Storage Targets (OSTs) respectively. A file is said to be striped when I/O operations access multiple OSTs concurrently. We enabled 4 possible concurrent accesses to each file, which allowed increasing the available I/O bandwidth significantly.

Computational performances (HPC)

The first learning from this billion mesh simulation is about NEPTUNE_CFD computational performances. NEPTUNE_CFD is powered by Code Saturne which was already known for its peta scale perfor mances [START_REF] Fournier | Optimizing code Saturne computations on petascale systems[END_REF][START_REF] Fournier | Recent and upcoming changes in code saturne: computational fluid dynamics HPC tools oriented features[END_REF]. During this simulation, we have been able to evaluate NEPTUNE_CFD HPC from 35 nodes up to 1000 nodes, i.e. from 1224 cores up to 36,000 cores on two supercomputers of last generation. The following analysis is based on the averaged CPU time value of 2 runs for each core number and we considered the mean value. No signif icant difference of CPU time is observed. NEPTUNE_CFD scalability has been evaluated on a restart simulation of 200 iterations after a transient step. To obtain an evaluation of computing performances, undisturbed by I/O operations, the first and last iterations are discarded from the analysis. During these remaining 198 iterations, no disk access occurs except for printing parallel performances (negligible overhead). Table 6 summarizes the effective CPU time required per timestep (aver aged over the 198 iterations). Computation elapsed time is very close to this time. On Olympe, we evaluated scalability from 35 up to 362 nodes and on Gaia, we started at 80 nodes up to 1000 nodes. As the two super computers use same processors, their performances can be directly compared despite slightly different architectures.

This study corresponds to a strong scaling study. Strong scaling is de fined as how the solution time varies with the number of cores for a fixed total problem size. It is hard to achieve a good strong scaling at larger process counts since the communication overhead increases in proportion to the number of processes used. The speedup is defined as the ratio between the elapsed time to execute a program on reference node number (here 1260 cores on Olympe) and on a set of concurrent n nodes (n × 36 cores) with n N 35. As architectures of both supercom puters are similar, the reference time from Olympe is also used to ana lyse speed up on Gaia. The efficiency is defined as the ratio between speedup and n the number of nodes. 30 s. Depending on their size, particle clusters have a life time ranging between 0.5 s up to 5 s which is in accordance with fre quency of apparition of bubbles ≤5 Hz as usually observed in the literature.

As stated previously, the simulation ran "only" for 25 s of physical time. In the first 5 to 10 s of transient regime, the bed is destabilized and the simulation loses track of initial conditions. Therefore, we have at best 15 s of useful data with the behavior of a bubbling fluidized bed reactor. This is not enough to access converged mean values and ex plains why we only focus on instantaneous values fields and unsteady values. Accessing meaningful statistics would require at least 200 s. Note however that, as far as clusters formation prediction is concerned, they are fully resolved as the simulation timestep (δt ≈ 2 × 10 -4 s) is significantly lower than their life time. Considering previously de scribed characteristic times, over the short simulated period, the total mass of particles can be considered stationnary. Nevertheless, the simu lated last 15 s are representative of hydrodynamics unsteadyness.

Fig. 15 presents the large particle solid fraction after 15 s. The poly ethylene volume fraction is high with a mean value inside the bed around 〈α l 〉 ≈ 0.45. Close to the wall, this value is higher between 0.56 and 0.6. In contrast, at the injector locations, the large particle vol ume fraction is null due to the gas and fine particles injection. At the top of the fluidized bed, we can see very small clusters and thread like structures. Only such a highly detailed computation could capture these small scale structures.

Fig. 16 also shows the large particle volume fraction but on a cut plane at the center of the reactor with two recursive zooms. Note the high local gradients of solid volume fraction α l . Along only a few cell, α l ranges between 0.0 and 0.6 which is close to the maximum solid vol ume fraction 0.64. We captured very small size structures such as bub bles and clusters (see Fig. 16) whose size is significantly larger than the mesh characteristic size. Though we may not claim that this simulation is independent from the mesh refinement, we got as close as currently possible from a fully resolved industrial scale reactor simulation.

If we focus on catalyst particles (see Fig. 17), the initial catalyst vol ume fraction is null and catalyst is injected through the four injectors. At the beginning the catalyst particles rise in the reactor and penetrate inside the dense fluidized bed, then after 10 s, catalyst particles are present everywhere in the reactor. Note that their volume fraction reaches up to α f = 10 -3 whereas they are injected with a lower fraction of α inj,f = 2.4 × 10 -5 which illustrates cluster formation (see Fig. 18).

One can observe a good dispersion of fine particles in the reactor with (i) an accumulation close to injectors, (ii) intermediate concentra tions over these injectors and up to the top of the fluidized bed and (iii) low values at the bottom and in the bulb. If the simulation ran a little while longer, we would expect to see fine particles going out of the re actor by elutriation even if the catalyst velocity in the disengagement zone is lower than in the fluidized bed. As illustrated by Fig. 19, large particles are mainly segregated in high density clusters (high peak on the right hand side at α l ∈ [0.5;/0.6]). A significant part of the volume is almost completely depleted from large particles (α l ∈ [0; 0.05]). In between (α l ∈ [0.05;0.50]), a flat distri bution of local volume fraction is observed.

One very interesting result is the thermal aspect. We took into ac count an exothermic reaction to be representative of the polymerization reaction through enthalpy source terms on particulate phases and gas heats up through particle gas enthalpy transfers. Fig. 20 plots the gas temperature (on the left) and the fine particle temperature. Fig. 21 shows axial temperature profiles for all three phases, at the center of the reactor, and near the wall on a x = 0 plane. Near the fluidiza tion inlet, gas has the set temperature of 50 °C and particles have tem perature in between 50 and 98 °C. Close to injectors at height z = 6.0 m, both gas and fine particles have a temperature of 50 °C and large particle temperature drops around 94 °C. Nevertheless, quickly after entering the reactor, all phase temperatures reach a narrow range between 97 °C and 101.6 °C. Being careful with the temperature scale, we can observe on Fig. 20 that a fluidized bed is an excellent mixer which yields a good temperature homogeneisation. Despite a total produced power of 20 MW, and catalyst particles that are very re active, the range between min and max temperatures in the bulk of the reactor is about a few degrees (≈4 °C).

Maximum gas and fine particles temperatures are observed right above injectors where α f is the highest, as fine particles are the most re active (see section 2 and Eqs. ( 6) and ( 7)) and will have had a long enough residence time to heat up. It is important to notice that inside the dense fluidized bed, high local temperature gradient may appear which would not be noticeable on coarser meshes. Nevertheless, the time scale for temperature evolution is significantly higher than hydro dynamics time scales. Reaching a thermic equilibrium would require simulating several hundreds of seconds.

All these transient fields seem to represent the expected physics of an actual fluidized bed. We are able to simulate at industrial scale a re active polydispersed pressurized fluidized bed reactor with an highly detailed mesh.

To be more quantitative, we can focus on wall pressure profiles, measured at t = 16.5 s as shown by Fig. 22. These profiles exhibit a linear slope with a maximum pressure drop of 63,000 Pa. A simple force bal ance, neglecting the friction force of particle on the wall, allows checking the total mass of particle:

M p P max S R g 102; 100 kg ð10Þ
with S R the section of the reactor cylindrical part. This total mass is con sistent with the initial mass of large particles (M l,0 = 100t). The influ ence of injectors in planes at 0°and 90°does not disturb these profiles significantly. The height of the bed can be deduced by extending the bed pressure profile up to the height axis, thus obtaining a bed height of ≈16.5 m.

From highly resolved simulations to sub grid modeling approaches

In the present paper we show that, nowadays, it is possible to per form a numerical simulation of an industrial scale geometry using a very fine mesh. However, such a highly detailed numerical simulation can only be performed with a very well parallized efficient code and it requires a huge amount of computational ressources. For engineers and researchers, computational ressources are limited and the number of accessible cores as well.

Since several years researchers have identified the effect of the small scale solid structures present in particulate flows which are not accounted for when the numerical simulation is performed with coarse mesh [START_REF] Schneiderbauer | CFD modeling and simulation of industrial scale olefin polymerization fluidized bed reactors[END_REF][START_REF] Cloete | Grid independence behaviour of fluidized bed reactor simulations using the two fluid model: detailed parametric study[END_REF]. Agrawal et al. [START_REF] Agrawal | The role of mesoscale structures in rapid gas-solid flows[END_REF], Igci et al. [START_REF] Igci | Filtered two-fluid models for fluidized gas-particle suspensions[END_REF], Parmentier et al. [START_REF] Parmentier | A functional subgrid drift velocity model for filtered drag prediction in dense fluidized bed[END_REF] and Özel et al. [START_REF] Özel | Development of filtered Euler-Euler two-phase model for circulating fluidised bed: high resolution simulation, formulation and a priori analyses[END_REF] show that, in dense and in circulating fluidized beds, when small scale solid structures are neglected, the first order effect appears on the drag force term (here in the term I g→p,i from Eq. ( 3)). Basically, the drag force is overestimated with a coarse mesh. As a consequence, in dense fluid ized beds, the use of a coarse mesh leads to an over prediction of bed height and in circulating fluidized beds the solid mass flux is overestimated. To overcome this effect and then to correctly predict the dynamics of particulate flows, researchers develop a new ap proach called Filtered Two Fluid Model (FTFM) [38 40]. Such an ap proach is widely inspired from the Large Eddy Simulation (LES) for single phase flows.

The use of a coarse mesh is similar to a filter applied on PDE equations. For a given function f(r) let us define the filtered value as f ðxÞ R f ðrÞGðx-rÞdr where G is the filter kernel. From this the filtered solid volume fraction reads

α p x ð Þ Z α p r ð ÞG x-r ð Þdr ð11Þ
and the gas and particle velocities

α g x ð ÞU g;i x ð Þ Z α g r ð ÞU g;i r ð ÞG x-r ð Þdr ð12Þ
α p x ð ÞU p;i x ð Þ Z α p r ð ÞU p;i r ð ÞG x-r ð Þdr: ð13Þ
In the following α k and U k;i are called the computed variables be cause those are the ones computed when a coarse mesh is employed. When such a filtering procedure is performed on the mathematical model given in Section 2, it leads to additional terms. As in single phase flow, the particle Reynolds stress can be decomposed as

g α p U p;i U p; j α p U p;i U p; j þ τ sgs ij ð14Þ
where the first term on the right hand side is the particle Reynolds stress expressed with the computed variables and τ ij sgs is the subgrid stress tensor that requires a model. However from a highly resolved nu merical simulation it is possible to measure τ ij sgs . Indeed, in the cylinder part of the reactor the mesh is nearly uniform and cartesian. Hence, we have extracted the data from a slice and a discrete spatial filter has been applied. Several filter widths Δ f have been applied up to 64Δ DNS where Δ DNS is the cell size of the highly resolved simulation. The results are shown by Fig. 23 where it can be seen that the instantaneous spatial average of shear components are nearly null. In contrast, the vertical

α p À Á ≈ g Δ f ; α p À Á R αmax 0 g Δ f ; α p À Á dα p : ð22Þ
Fig. 27 shows the evolution of hðα p Þ for several values of filter width. The different curves merge into a unique function. In the literature two models have been proposed. First, Parmentier et al. [START_REF] Parmentier | A functional subgrid drift velocity model for filtered drag prediction in dense fluidized bed[END_REF] from 2D highly resolved numerical simulation of dense fluidized proposed

h α p À Á - α p α max s 1-u ð Þ 2 1-C h;2 α p α max þ C h;3 α p α max 2 " # ð23Þ
where u α p =α max and the constant are set to C h,2 = 1.88 and C h,3 = 5.16. On another hand, Özel et al. [START_REF] Özel | Development of filtered Euler-Euler two-phase model for circulating fluidised bed: high resolution simulation, formulation and a priori analyses[END_REF] used 3D highly resolved numeri cal simulations of a periodical circulating fluidized bed and proposed

h α p À Á -tanh α p C h;1 α p α max s 1-u ð Þ 2 1-C h;2 α p α max þ C h;3 α p α max 2 " # ð24Þ
with C h,1 = 0.1. Fig. 27 shows that gðΔ f ; α p Þ= R αmax 0 gðΔ f ; α p Þdα p is inde pendent from the filter width. Also Fig. 27 exhibits that the function pro posed by Parmentier et al. [START_REF] Parmentier | A functional subgrid drift velocity model for filtered drag prediction in dense fluidized bed[END_REF] is in good accordance with measures from the highly resolved industrial scale numerical simulation.

Conclusion and outlook

To describe all flow scales in a polymerization fluidized bed reactor, very fine simulations are required. In the present study, one demon strates that the simulation of a fluidized bed at industrial scale with a billion mesh cells using to 36,000 cores is now achievable. This has been possible thanks to recent advances in the computational power of supercomputers and the developpement of efficient massively paral lel CFD solvers. This is the first time that such a detailed insight view of an industrial reactor has been obtained. Our previous largest simulation used a meshing ten times smaller than the current one [START_REF] Hamidouche | Numerical simulation of multiphase reactive flows[END_REF]. Now, thanks to the achieved size of mesh cells, we may consider simulating complex geometries of industrial reactors. Due to the novelty of performing that large simulations, new challenges appear, both in terms of HPC, and in the exploitation of the simulation results. Tackling these challenges required a close collaboration between research labo ratories, supercomputing centers and an industrial partner which de velops a massively parallel solver up to date in terms of available models and numerical methods.

HPC challenges occured at all steps of this project. First, the preprocessing of the simulation was tricky. To generate and partition an unstructured mesh with 1,002,355,456 hexahedral cells, a multi software procedure has been specifically designed. Multiple partitio ning methods have been compared in terms of stability, performance and quality of the domain decomposition. Second, the simulation ran on up to 36,000 cores on newly installed supercomputers which had to be fine tuned to handle that massive computations. The simulation 23) and the dashed line to Eq. ( 24) both computed with α max = 0.592 and normalized by their own integral.

had impressive metrics: more the 200 TB of data have been generated for 25 s of simulated physical time, each checkpoint restart file was about 1.3 TB and reading it on the whole of the supercomputer Olympe required up to 13 min. Postprocessing was to be considered before the simulation ran to limit the volume of generated data. Finally, the trans fer and post processing of such a crushingly high volume of data was a challenge by itself.

Though being feasible, this kind of simulation still requires tremen dous resources and is not easily performed. Each step (simulation setup, running and post processing) required about a month each, and expensive ressources. For this reason and for now, the goal of these large scale simulations must be to generate reference databases. Ideally, such a reference simulation should be longer than the one we performed, to ensure the convergence of averaged results. This was however the best that we could obtain for now, and supplementary advances will be required to reach converged fully resolved simula tions. Once reference results are generated, their exploitation may serve to improve the quality of lower cost simulations. This is achievable by exploiting fully resolved simulation results to measure sub grid statistics such as fluid particle drag, velocity or volume fraction vari ance. Sub grid models [17,29 32] are then tested a priori and a posteriori against the fully resolved simulation. This allows accounting for the effect of small and meso scale structures even when a coarse mesh is used.

As stated previously, one significant challenge has been the manage ment of significantly large datasets. In terms of actual storage space, I/O operations dedicated time, site to site transfer rates and associated human working time, we again reached the upper limit of what was fea sible. A strategy to significantly reduce these issues would be to switch the whole post processing to in situ visualization and co processing; i.e. post processing performed continuously during the simulation with tools like Catalyst ParaView [START_REF] Bauer | In SituMethods, infrastructures, and applications on high performance computing platforms[END_REF]. By switching toward these co processing approaches, storage footprint would significantly decrease which goes along with good practices for "green" supercomputing practices [START_REF] Feng | The Green500 list: encouraging sustainable supercomputing[END_REF]. For the past few decades, we analyzed computational performances for the simulation of the considered fluidized bed reactor mainly in terms of raw computational power (s simulated /h CPU ). However, in line with concerns for green practices, a better metrics than the one we presented previously would be the consider energetic performances (s simulated /J electrical ).

Overall, this work is a premiere which shows that fully resolved simulations of industrial scale fluidized bed reactors are now possible. They are of academic and industrial interest and should become more and more common in the next few years. Some of the challenges that we encountered should then become more easily tackled with solutions that starts becoming more accessible. With the expected evolution of HPC ressources, the numerical cost of this scale of simulation should de crease and reach a reasonnable level.

Massively parallel numerical simulation using up to 36 ,

 36 000 CPU cores of an industrial-scale polydispersed reactive pressurized fluidized bed with a mesh of one billion cells Hervé Neau, Maxime Pigou, Pascal Fede, Renaud Ansart, Cyril Baudry, Nicolas Mérigoux, Jérome Laviéville, Yvan Fournier, Nicolas Renon, Olivier Simonin To cite this version: Hervé Neau, Maxime Pigou, Pascal Fede, Renaud Ansart, Cyril Baudry, et al.. Massively parallel numerical simulation using up to 36,000 CPU cores of an industrial-scale polydispersed reactive pressurized fluidized bed with a mesh of one billion cells. Powder Technology, 2020, 366, pp.906-924. 10.1016/j.powtec.2020.03.010. hal-02735376

Fig. 7 .

 7 Fig. 7. Mono-layer O-Grid mesh generated with Simail [67] -476,928 hexahedral cells -850 cells along a diameter.

  7), and the top one contains 253,184 hexahedron. The three first meshes have 850 cells along a diameter. 2. Using the extrude functionality of Code Saturne [66], form 4 cylin ders from the four mono layer meshes (see Fig. 8: (1a): 1451 layers up to z = 16m, (1b): 471 layers up to z = 23m, (1c): 99 layers up to z = 25m and (1d): 152 layers up to z = 28m).

Fig. 9 .

 9 Fig. 9. Illustration of the final billion cell mesh. The greyed zoomed circle area represents a single injector.

Fig. 8 .

 8 Fig. 8. 3D meshes transformations with Code_Saturne. (1) mono-layer meshes extrusion; (2) Radial dilation; (3) and (4): bulb formation through vertices coordinates shifting; (5): submesh merging.

Fig. 10 .

 10 Fig. 10. Illustration of mesh partitioning of the geometry into 64 domains. (a) PT-SCOTCHexterior view. (b) PT-SCOTCH -centered cut plane. (d) Morton -exterior view. (d) Morton -centered cut plane.

Table 4

 4 Partitioner comparison. Stability indicates whether the partitioning step succeeded. Partitioning time, and the duration of a reference simulation that ran on obtained partionnings (CFD CPU time) are normalized by the results for Morton method. with 50 nodes 54,271-115,996 31,200-63,732 53,213-118,784 -Hilbert curve in terms of neighboring unbalance, ghost cell unbalance and simulation CPU time.

Fig. 11 .

 11 Fig. 11. (left) Illustration of the Morton (or Z) space-filling curve, used to sort mesh cells by the Morton partitioner. (right) Illustration of the Hilbert space-filling curve, used to sort mesh cells by the Hilbert partitioner.

Fig. 16 .Fig. 17 .

 1617 Fig. 16. Instantaneous large particle volume fraction on a cut plane with successive ×4 and ×10 zooms. White lines one the right sub-figure are the actual mesh contours. Note that the color has not been interpolated between mesh cells, hence the pixelated aspect of the right-most sub-plot.

Fig. 18 .

 18 Fig. 18. Instantaneous catalyst particle volume fraction after 15 s.

Fig. 19 .

 19 Fig. 19. Probability Density Function of large particles volume fraction. The data have been extracted from one vertical cut plane at time t = 16.5 s.

Fig. 20 .

 20 Fig. 20. Gas and fine particle temperatures after 15 s. Note the color-scale that is only linear between 98.00 °C and 101.23 °C.

Fig. 21 .ig. 22 .

 2122 Fig.[START_REF] Bennani | Numerical simulation of unsteady dense granular flows with rotating geometries[END_REF]. Axial temperature profiles for all three phases at the center of the reactor and near the wall on a x = 0 plane.

  h

Fig. 25 .

 25 Fig. 25. Correlation diagram of the drift velocity in z-direction with respect to the computed gas-particle relative velocity. The filter width is Δ f /Δ DNS = 24. Only 10% of points are represented.

Fig. 26 .

 26 Fig. 26. Function gðΔ f ; αpÞ computed with Eq. (21).

Fig. 27 .

 27 Fig. 27. Function hðαpÞ with respect to filtered solid volume fraction. The solid line corresponds to Eq. (23) and the dashed line to Eq. (24) both computed with α max = 0.592 and normalized by their own integral.

Table 2

 2 Boundary conditions.

		Fluidization inlet	Catalyst injectors	Outlet	Wall
	Gas	Q g = 756 t/h	Q inj,g = 19.3 t/h/injector	Free flow with imposed gauge pressure.	Friction model.

Table 1

 1 Powder and gas physical properties.

		Gas	Large particles	Fine particles
	Density at 24 bars (kg ⋅ m 3 )	22	850	850
	Viscosity×10 5 at 24 bars (Pa ⋅ s)	1.54	-	-
	Mean diameter (μm)	-	1600	80
	Restitution coefficient	-	0.9	0.9
	Specific heat (J ⋅ kg 1 ⋅ K 1 )	1728	2000	2000
	Thermal conductivity (W ⋅ m 1 ⋅ K 1 )	0.04	-	-

Table 3

 3 Mesh quality parameters.

	Billion cell mesh
	Cell Number	1,002,355,456
	Interior face number	3,005,213,856
	Boundary face number	3,705,024
	Vertices number	1,004,210,878
	Minimum control volume (m 3 )	6 × 1 0 8
	Maximum control volume (m 3 )	7 × 1 0 6
	Total domain volume (m 3 )	633.0
	Orthogonality quality	0 bad cell (0%)
	Offset	344 bad cells (≪0.001%)
	Least-Squares Gradient Quality	11,776 bad cells (≪0.001%)
	Cells Volume Ratio	0 bad cell (0%)

Table 5

 5 Supercomputer main characteristics.

	Olympe -peak performance 1.365 Pflop/s
	13,392 cores (2.3GHz)
	374 compute nodes	Cores / node Processor -GPU	RAM/node	1 node
	360 bi-socket	2 × 18 -Intel® Xeon® Gold Skylake 6140	192 GB	2.65 TF
	12 bi-socket + GPU	2 × 18 -Skylake 6140-4 GPU NVIDIA V100	384 GB	33.8 TF
	2 bi-socket memory	2 × 18 cores SKYLAKE 6140	1536 GB	2.65 TF
	Gaia -peak performance 3.05 Pflop/s
	127th TOP500 (06/2019) -42,912 cores (2.3GHz)
	1224 compute nodes	Cores / node Processor -GPU	RAM/node	1 node
	1192 bi-socket	2 × 18 -Intel® Xeon® Gold Skylake 6140	192 GB	2.65TF
	32 bi-socket + GPU	2 × 4 Skylake 5122-2 GPU NVIDIA V100	384 GB	18 TF
	48 bi-socket memory	2 × 18 -Intel® Xeon® Gold 6140	384 GB
	32 bi-socket memory	2 × 4 -Xeon Gold 5122	384 GB

www.calmip.univ-toulouse.fr

Acknowledgments

This work was granted access to the HPC resources of CALMIP supercomputing center under the allocation 2018 mesochallenge and EDF Gaia cluster during EDF Grands Challenges.

Declaration of Competing Interest

None.


Florian Sabatier

Renaud Ansart
email: renaud.ansart@ensiacet.fr

Huili Zhang

Jan Baeyens

Olivier Simonin

Experiments support simulations by the NEPTUNE_CFD code in an Upflow Bubbling Fluidized Bed reactor

• Fluidization phenomena differ with the tube length (from freely bubbling to slugging).

• Results are used to validate a 3D si- mulation approach by the NEPTUNE_CFD code.

• The effect of mesh refinement on the slugging behavior prediction is discussed.

• It is demonstrated that experiments and simulation predictions are in good agreement.

G R A P H I C A L A B S T R A C T A B S T R A C T

Long tubes with small internal diameter find increasing applications in indirect concentrated solar receivers using an Upflow Bubbling Fluidized Bed of Geldart-A powders as heat carrier. Although successfully demonstrated for tubes of 0.5 to 1 m length, longer tubes are required to increase the solar energy capture efficiency and capacity. The fluidization phenomena differ with the tube length, and freely bubbling fluidization will be transformed into slugging, thus hampering the heat transfer. The behavior of Geldart-A powders in tall tubes of small I.D. has not been extensively studied. The research experimentally investigated the different fluidization modes in a 4 m long tube, and demonstrated the occurrence of freely bubbling at the bottom section of the bed, and slugging from a bed depth in excess of about 1.25 m. Slug characteristics (frequency, length, velocity) were measured and correlated. The results were used to validate 3D numerical simulations based on an Euler-Euler approach in the NEPTUNE_CFD code applied to a fine mesh of 15,000,000 cells. A positive match between experimental and simulation results concerning frequency and velocity of large bubble structures was obtained. The effect of mesh refinement on the slugging behavior prediction was discussed.

Introduction

Gas/solid fluidization powders can involve C (cohesive), A (aeratable), B (bubbling) and D (coarse) particles [START_REF] Kong | Bubbling and slugging of Geldart Group A powders in small diameter columns[END_REF][START_REF] Baeyens | An investigation into slugging fluidized beds[END_REF][START_REF] Geldart | Gas Fluidization Technology[END_REF]. The different powder groups can be expressed in terms of the relevant particle size, d p , and gas/solid densities, g and s [START_REF] Kong | Bubbling and slugging of Geldart Group A powders in small diameter columns[END_REF][START_REF] Baeyens | Elutriation of fines from gas fluidized beds of Geldart A-type powders ' effect of adding superfines[END_REF].

For cristobalite sand, used in the present research ( ~2640 kg m s 3 ), C powders will be smaller than 33 µm; the aeratable group A will cover the range of 33 to 86 μm; group B powders will range from 86 to 446 μm; and group D powders will have a particle size in excess of 446 µm. Although aeratable group A powders are widely used in com- mercial fluidized bed catalytic reactors [START_REF] Kunii | Fluidization Engineering[END_REF], their recent application in Upflow Bubbling Fluidized Bed (UBFB) or Cross-flow Bubbling Fluidized Bed (CfBFB) indirect solar receivers opens new application perspectives [START_REF] Kong | Solids mixing in a shallow crossflow bubbling fluidized bed[END_REF][START_REF] Kang | Energy analysis of a particle suspension solar combined cycle power plant[END_REF][START_REF] Zhang | Experiments support an improved model for particle transport in fluidized beds[END_REF][START_REF] Zhang | High-efficiency solar power towers using particle suspensions as heat carrier and thermal energy storage[END_REF][START_REF] Zhang | The use of particle heat carriers in the Stirling engine concept[END_REF][START_REF] Zhang | Particle circulation loops in solar energy capture and storage: gas-solid flow and heat transfer considerations[END_REF] due to their high mobility [START_REF] Baeyens | Bed expansion and the visible bubble flow rate in gas fluidized beds[END_REF], high wall-to-bed heat transfer rates [START_REF] Baeyens | Modelling approach to the effect of equipment scale on fluidized bed heat transfer data[END_REF], and their easy circulation around the concentrated solar energy loops. Gas bubbles in group A powders moreover split and recoalesce frequently, resulting in a maximum stable bubble size [START_REF] Kong | Bubbling and slugging of Geldart Group A powders in small diameter columns[END_REF][START_REF] Valverde | Fluidization of fine and ultrafine particles using nitrogen and neon as fluidizing gases[END_REF], as further discussed in the paper. The bubble size is independent of the average particle size in Geldart-A powders. The bubble rise velocities appear to be ~0.3 to 0.4 m s -1 , irrespective of the bubble size, suggesting that the gross powder circulation controls the rise velocity [START_REF] Kong | Bubbling and slugging of Geldart Group A powders in small diameter columns[END_REF][START_REF] Valverde | Fluidization of fine and ultrafine particles using nitrogen and neon as fluidizing gases[END_REF][START_REF] Ye | Two-fluid modeling of Geldart A particles in gas-fluidized beds[END_REF]. Bubbles will transform into slugs at sufficiently high superficial gas velocity or in deep beds, as illustrated in Fig. 1. Slugs reduce particle mixing and minimize the gas/solid contacting [START_REF] Baeyens | Modelling approach to the effect of equipment scale on fluidized bed heat transfer data[END_REF] and reduce the wall-to-bed heat transfer coefficient [START_REF] Baeyens | Bed expansion and the visible bubble flow rate in gas fluidized beds[END_REF]. Bubbling and slugging are distinct hydrodynamic modes in a fluidized bed and should therefore be jointly considered in predicting chemical conversions and heat transfer. At velocities in excess of 0.5 m/s, depending on the particle size, a turbulent or ultimately circulating fluidized bed mode will occur and particle recycling is required to maintain the bed solids' inventory

Nomenclature

mm 1 m U U mf

The novelties of the present paper and the objectives are hence clearly defined, and result from the use of extensive experimental data in support of a dedicated NEPTUNE_CFD code simulation. The objectives involve a stepwise approach that:

(i) fully assesses properties of the system under scrutiny, with special emphasis on secondary particle properties such as size distribution and sphericity; (ii) studies the hydrodynamic behavior of selected particles experimentally in a small I.D. tube of 4 m height; (iii) introduces relevant simulation conditions into the NEPTUNE_CFD code, to define simulation geometry and mesh, fix boundary conditions (inlet, outlet, tube walls), and assesses the influence of the relevant particle and air properties; (iv) applies an Euler-Euler model, with governing equations described by Ansart et al. [START_REF] Ansart | Dense gas-particle suspension upward flow used as heat transfer fluid in solar receiver: PEPT experiments and 3D numerical simulations[END_REF]; (v) validates the approach through pressure gradients and solid volume fraction; (vi) compares experimental and simulation results towards slug formation, slug frequency, velocity, and occupancy profiles along the tube diameter and height; (vii) discusses the effect of refining the mesh characteristics on the slugging behavior prediction for Geldart-A powders; (viii) formulates fluid-dynamic conclusions and recommendations.

Experimental setup and procedures

Equipment and powders used

Within the concept of a "particle-in-tube" UBFB, a small I.D. cylindrical column was used of 50 mm I.D. and of a total height of 4 m, hence complementing the experimental research of Kong et al. [START_REF] Kong | Bubbling and slugging of Geldart Group A powders in small diameter columns[END_REF], limited to a 2 m high tube. Experimental procedures were detailed by Kong et al. [START_REF] Kong | Bubbling and slugging of Geldart Group A powders in small diameter columns[END_REF] and involve an operating air flow rate in the range of 0 to 5 × 10 -4 kg s -1 . The fluidized bed is fed by compressed air at ambient temperature and the bed exhaust is connected to atmosphere. The distributor plate had a pressure drop of 1560 Pa at a superficial air velocity of 0.1 m s -1 , significantly exceeding the common design recommendations [START_REF] Geldart | The design of distributors for gas-fluidized beds[END_REF]. The Geldart-A powder used was cristobalite silica. Particle properties were defined by the methods described by Kong et al. [START_REF] Kong | Bubbling and slugging of Geldart Group A powders in small diameter columns[END_REF], with mean particle size of 55 µm, absolute density of 2640 kg m 3 and tapped bulk density [START_REF] Geldart | Gas Fluidization Technology[END_REF] of 1380 kg m 3 . U mf was Fig. 2. Schematics of a solar tower combined cycle generation process [START_REF] Kang | Energy analysis of a particle suspension solar combined cycle power plant[END_REF].

[ [START_REF] Zhang | Operation diagram of circulating fluidized beds (CFBs)[END_REF][START_REF] Pitie | Circulating fluidized bed heat recovery/storage and its potential to use coated phase-change-material (PCM) particles[END_REF][START_REF] Chan | Particle velocities and their residence time distribution in the riser of a CFB[END_REF]. Small Internal Diameter (I.D.) beds are normally used in laboratory scale experiments, or represent the constrained bed cells in between vertical heat exchanger tubes in a catalytic reactor [START_REF] Kunii | Fluidization Engineering[END_REF]. The current development of particle suspension solar heat receivers use such small I.D. beds with Geldart-A powders as heat carriers [START_REF] Zhang | Experiments support an improved model for particle transport in fluidized beds[END_REF][START_REF] Zhang | High-efficiency solar power towers using particle suspensions as heat carrier and thermal energy storage[END_REF]. Fig. 2 depicts the layout of a MW th hybrid combined cycle solar tower power plant, to be tested on pilot-scale at the CNRS-Themis solar tower (France), without the steam block [START_REF] Kong | Solids mixing in a shallow crossflow bubbling fluidized bed[END_REF][START_REF] Kang | Energy analysis of a particle suspension solar combined cycle power plant[END_REF][START_REF] Zhang | High-efficiency solar power towers using particle suspensions as heat carrier and thermal energy storage[END_REF]. The pilot-scale set-up is described elsewhere and includes the particle suspension receiver, hot particle storage [START_REF] Zhang | High-efficiency concentrated solar power plants need appropriate materials for high-temperature heat capture, conveying and storage[END_REF], and crossflow fluidized bed air heater to preheat the turbine compressed air [START_REF] Kong | Solids mixing in a shallow crossflow bubbling fluidized bed[END_REF][START_REF] Kang | Energy analysis of a particle suspension solar combined cycle power plant[END_REF]. Exhausts from the turbine are fed to a steam block. The design of the receiver is of paramount importance in order to be able to operate at as high a temperature as possible (750 to 850 °C).

The present work considers the scale-up of the Upflow Bubbling Fluidized Bed (UBFB) receiver concept. Previous papers have provided information concerning the on-sun testing of a single tube (36 mm I.D., 0.5 m long) and of a 150 kW th pilot module consisting of 16 parallel tubes of I.D., each long, using group A powders fluidized at superficial gas velocities of ~0.03 to 0.20 m s -1 [START_REF] Zhang | High-efficiency solar power towers using particle suspensions as heat carrier and thermal energy storage[END_REF][START_REF] Zhang | Particle circulation loops in solar energy capture and storage: gas-solid flow and heat transfer considerations[END_REF][START_REF] Flamant | Dense suspension of solid particles as a new heat transfer fluid for concentrated solar thermal plants: on-sun proof of concept[END_REF][START_REF] Perez-Lopèz | On-sun operation of a 150kWth pilot solar receiver using dense particle suspension as heat transfer fluid[END_REF]. A further scale-up to multi-MW capacity will involve the use of the parallel tubeconcept with larger I.D. tubes and of taller height towards an optimum use of the concentrated heliostat beams and provide sufficient heat capture surface area. Although receiver heights of 10 to 18 m are now common practice in molten salt applications [22], the gas-solid hydrodynamics in the UBFB will limit the maximum applicable height in function of the tube geometry, with special emphasis on its height, due to the occurrence of slugging (Fig. 1). The onset of slugging mostly depends upon the tube diameter and the excess gas velocity, i.e. [START_REF] Kong | Bubbling and slugging of Geldart Group A powders in small diameter columns[END_REF]. Since slugs have a limited bubble-induced particle mixing, they reduce the wall-to-bed heat transfer [START_REF] Baeyens | Modelling approach to the effect of equipment scale on fluidized bed heat transfer data[END_REF]. It is hence important to investigate the bubbling/slugging effects i n a g eometry s imilar t o the currently scaled-up particle-in-tube concept. P 0.34

U mf U mb U mf U mb U t

Measuring techniques

The behavior of fluidized beds at the onset of slugging and full slugging mode was assessed by three different methods, as previously used by Baeyens and Geldart [START_REF] Baeyens | An investigation into slugging fluidized beds[END_REF] and Kong et al. [START_REF] Kong | Bubbling and slugging of Geldart Group A powders in small diameter columns[END_REF]. Visual observations, pressure drop and its fluctuations and high speed camera imaging were jointly used and assessed. Whereas visual observations cannot be used in the bottom (bubbling) section of the bed due to overlapping bubbles, and pressure drop measurements show significant fluctuations, the combined analyses of the three methods provided a satisfactory insight into the full slugging mode. Photographic illustration of the observed wall and axi-symmetric slugs is given in Fig. 3.

Experimental procedure

The bed heights at U mf and U mb were measured. The gas velocity was then increased to well into the slugging regime. Flow rate readings, visual and P observations, and images were made for given periods of time (between 2 and 3 min). All data points presented in the figures are averages of 3 repeat experiments. The deviation among 3 repeat results and the averages was below 5%. The evolution of expanded bed height (H) in relationship with the bed height at minimum bubbling (H mb ) with excess gas velocity is illustrated in Fig. 4. The time average pressure drop across the whole bed P, and the pressure fluctuations, presented as standard deviation to the average pressure drop P ( ), were recorded and presented as ratio in Fig. 5. Leva et al. [START_REF] Leva | Fluid Flow Through Packed and Fluidized Systems[END_REF] tentatively predicted the onset of slugging when pressure fluctuations exceed 5 per cent of the pressure drop through the bed at the given gas flow rate. This rule of thumb was in fair agreement with our observations. At U U ( ) mf values beyond the onset of slugging, big bubbles rapidly transformed into wall-slugs, which prevailed throughout the bed height for about 1.25 m. They are subsequently transformed into axi- symmetric slugs in deeper beds. The axi-symmetric slugging mode hampers the heat transfer from the receiver wall to the bed (poor solids mixing and long contact of an air slug with the wall).Further experimental results will be given in Section 4 below.

Numerical simulation description

Mathematical model 3.1.1. Eulerian modeling approach

To access to the particle behavior inside the receiver tubes, two CFD approaches (Eulerian-Lagrangian or Eulerian-Eulerian) can be used. The Eulerian-Lagrangian approach treats particles as discrete elements when the fluid flow is characterized by Eulerian equations. The motion for each individual particle is solved, with inclusion of the effects of determined from a common versus U graph, and was cm s -1 . The minimum bubbling velocity of 0.53 cm s -1 corresponded with the superficial air velocity when first bubbles could be seen erupting at the bed surface. The accuracy on was 0.01 cm s -1 while being 0.02 cm s -1 for . Values of and are clearly different. The particle shape was determined by OPTEC DV 320 microscope imaging and the sphericity ( ) varying between 0.64 and 0.71 was calculated from Cavarretta et al. [START_REF] Cavarretta | Applying 2D shape analysis techniques to granular materials with 3D particle geometries[END_REF] and Cho et al. [START_REF] Cho | Particle shape effects on packing density, stiffness, and strength: natural and crushed sands[END_REF].

was calculated as 0.203 m s -1 according to the method described by Geldart [START_REF] Geldart | Gas Fluidization Technology[END_REF]. For the minimum fluidization height of 230 cm, the column was filled with 5 kg of solid. = 0.55

s min , = 0.64 s max , q p 2 3.1.

Consideration of powder PSD and non-sphericity

As shown in Figs. 6 and7, the powder used presents a very broad particle size distribution. Since particles follow a Gaussian-like monomodal size distribution, the complete size distribution was simplified to a single particle class with two separate diameters: (See Fig. 8)

= d d d drag [3; 0] 3 [1; 0] (1)
d drag for the drag coefficient term, as defined in Eq. 1 and = d d p 50 , the median diameter, for every other equation (momentum transfer, random kinetic energy, collisions) as previously applied by Sabatier et al. [START_REF] Sabatier | Bubbling of Geldart Group A Powders in Small Diameter Columns: Experiments and 3D Numerical Simulations focussing on the particle size distribution and non sphericity[END_REF]. The drag term coefficient of Gobin et al. [START_REF] Gobin | Fluid dynamic numerical simulation of a gas phase polymerization reactor[END_REF] was modified to include effects of non-sphericity based on the work of Bagheri and Bonadonna [START_REF] Bagheri | On the drag of freely falling non-spherical particles[END_REF] as expressed in Eqs. [START_REF] Baeyens | An investigation into slugging fluidized beds[END_REF]. Two coefficients k N and k S were introduced, correcting the drag term in the Stokes and Newton regime respectively, and determined by the measurement of two shape factors defined by Bagheri and Bonadonna [START_REF] Bagheri | On the drag of freely falling non-spherical particles[END_REF].

= < C Re C C C ( ) if 0.7 min[ , ] else 0.7 d p d WY g d WY d Erg g , , , nsph nsph nsph
(2.a)

= + + + C k Re Re k k k 24 1 0.125 0.46 1 d WY S p p N S N Re k k g , 2/3 5300 / 1.7 nsph p N S (2.b) = + = C k Re k Re v d µ 200 (1 ) 7 3 , | | d Erg S g p N p g g r p p g , (2.c) = < > = V v U U rp i r i p p i g i , , , , (2.d
)

Geometry and mesh

Geometric model

The geometric model used for simulation was the experimental 4 m long pipe of 5 cm I.D., followed by a convergent head. To prevent the occurence of a local exit pressure drop that could result in a numerical error for the code, the convergent head (although not used experimentally) accelerates the particle suspension before the tube exhaust, in order to prevent occasional backflow and thus numerical artifacts. particle collisions and forces acting on the particles due to the gas [START_REF] Xie | Eulerian-Lagrangian method for threedimensional simulation of fluidized bed coal gasification[END_REF]. For this method to be very accurate, every particle has to be taken into account, likewise every particle-particle and particle-wall collisions. Since the calculation time is directly linked to the number of particles simulated, in the specific case studied here, with a particle size less than , more than 10 11 particles need to be accounted for. This approach is clearly impossible to solve with current calculation resources. The Eulerian-Eulerian approach treats particle and fluid p hases a s interpenetrating yet separated and continua; although coupled Eulerian transport equations for mass, momentum, and energy are written for both phases. This modeling approach allows the simulation of large fluid-particle s ystems, i ncluding r eactive fl ow, he at tr ansfers [START_REF] Hamidouche | Numerical Simulation of Multiphase Reactive Flows[END_REF][START_REF] Benoit | Threedimensional numerical simulation of upflow bubbling fluidized bed in opaque tube under high flux solar heating[END_REF] and polydisperse systems [START_REF] Fotovat | Sand-assisted fluidization of large cylindrical and spherical biomass particles: experiments and simulation[END_REF][START_REF] Fede | 3D Numerical Simulation of Polydisperse Pressurized Gas-Solid Fluidized Bed[END_REF]. In this approach, particles size and inertia are linked to the mesh size: with very small and light particles, calculation time may still be a difficulty to overcome.

Three-dimensional numerical simulations of the fluidized bed were performed using the NEPTUNE_CFD code. This Eulerian n-fluid unstructured parallelized multiphase flow s oftware [START_REF] Méchitoua | An unstructured finite volume solver for two phase water/vapour flows based on an elliptic oriented fractional step method[END_REF] i s d eveloped in the framework of the NEPTUNE project, rallying CEA (Commissariat á l'Energie Atomique), EDF , IRSN (Institut de Radioprotection et de Sûreté Nucléaire), and Framatome. The modeling approach for turbulent polydisperse gas-solid flows i s i mplemented b y t he I nstitut de Mécanique des Fluides de Toulouse (IMFT) and the Laboratoire de Génie Chimique (LGC) and has been developed for High Performance Computing [START_REF] Hamidouche | Numerical Simulation of Multiphase Reactive Flows[END_REF].

Mathematical models implemented

The mathematical approach used is a hybrid method [START_REF] Morioka | Modeling of gas and solid particles two-phase flow and application to fluidized bed[END_REF] deriving transport equations by phase ensemble average for the fluid phase and including the effects of the interstitial fluid towards the dispersed phase by the kinetic theory of granular flows [START_REF] Boelle | Second-order prediction of the particle-phase stress tensor of inelastic spheres in simple shear dense suspensions[END_REF]. The frictional stress tensor is taken into account in the particle stress tensor using the model proposed by Srivastava and Sundaresan [START_REF] Srivastava | Analysis of a frictional-kinetic model for gas-particle flow[END_REF], with a threshold value for the frictional viscosity activation and a maximal particles compaction of . Transport equations (mass, momentum and fluctuating k inetic e nergy) a re s olved f or e ach p hase a nd coupled through interphase transfer terms. The drag law used is based on two different correlations (Wen and Yu [START_REF] Wen | Mechanics of fluidization[END_REF] and Ergun [START_REF] Ergun | Fluid flow through packed columns[END_REF]) depending on the value of the solid volume fraction [START_REF] Gobin | Fluid dynamic numerical simulation of a gas phase polymerization reactor[END_REF][START_REF] Simonin | Drag force modelling in dilute to dense particle-laden flows with mono-disperse or binary mixture of solid particles[END_REF]. The gas flow i s considered as laminar because the variance effects o f g as l ocal v elocity in the equation of momentum are negligible likened it to the drag term. Finally, , the particle random kinetic energy, is integrated in a transport equation through the solid phase. All mathematical equations are described by Ansart et al. [START_REF] Ansart | Dense gas-particle suspension upward flow used as heat transfer fluid in solar receiver: PEPT experiments and 3D numerical simulations[END_REF].

Geldart-A particles are highly sensitive to the mesh cells size, but a mesh-independent solution based on the bed expansion is possible [START_REF] Agrawal | The role of meso-scale structures in rapid gas-solid flows[END_REF][START_REF] Parmentier | A functional subgrid drift velocity model for filtered drag prediction in dense fluidized bed[END_REF]. A sensitivity study was performed on columns of the same diameter but only of 1 m height. This study was based on the convergence of the spatially averaged solid fraction s , depending on the mesh nodes number on a diameter. For each case, = r D (nodes on a diameter 1)/ and = z r 2 . Table 1 shows the evolution of the solid volume fraction s , spatially averaged over the bed. As the mesh cells size decreases, the spatially averaged solid volume fraction increases, thus reflecting a bed height decrease. Beyond 80 nodes along the bed diameter, the average solid volume fraction is almost constant (less than 1% variation). Due to computational resources, the rest of the study will hence focus on a 80 nodes along the bed diameter mesh. Finally the mesh used for the simulation is built from about 15,000,000 hexahedrons of approximately = r 0.63 mm and = z 0.89 mm.

Phase properties

The particles used for the simulation were cristobalite silica, fluidized by air at atmospheric pressure (101325 Pa) and ambient temperature (298 K). Due to the bed height, the pressure gradient along the tube height cannot be neglected and the gas density cannot be considered as constant but needed to be computed using the perfect gas law. It is noteworthy that cinematic viscosity is automatically computed by the code from density and dynamic viscosity. The initial solid mass in the column for numerical simulations was 5 kg. Phase physical properties used for the simulation are summarized in Table 2.

Boundary conditions

Due to the mesh size, numerical simulations were highly CPU-time consuming and therefore only one intermediary velocity case has been studied. The industrial gas velocity operating range for the UBFB in solar receivers ranges from 7 to 56× U mf , with an average fluidization velocity of 25 × U mf . The distributor was considered as a horizontal wall for the particles and an inlet with an imposed uniform superficial velocity for the gas, equals to

= × = U U 25 0.085 mf m s -1 .
Due to the expansion of the gas flow towards the top of the bed, the superficial gas velocity will increase and close to 30 × = U 0.102 mf m s -1 . At the top of the fluidized bed, a free outlet for both gas and particles was defined. The same wall-type boundary condition was applied for both the particles at the distributor, and for both phases at the side walls: this condition involves no-slip for velocity, and zero-flux for random kinetic energy for the solid phase, corresponding to a large wall roughness or a high non-sphericity of the particles [START_REF] Fede | 3D numerical simulation of a lab-scale pressurized dense fluidized bed focussing on the effect of the particle-particle restitution coefficient and particle-wall boundary conditions[END_REF]. It is noteworthy that no-slip boundary condition was assumed over partial-slip boundary condition due to the lack of tube roughness measurement and thus the impossibility to accurately justify for a value of the models friction and bouncing parameters. In addition, the no-slip boundary condition, considering elastic frictionless bouncing at the walls, and thus with a particle-to-wall restitution coefficient = e 1 w , was found to in good agreement with PEPT experimental measurements for dense fluidized beds [START_REF] Fede | 3D numerical simulation of a lab-scale pressurized dense fluidized bed focussing on the effect of the particle-particle restitution coefficient and particle-wall boundary conditions[END_REF].

Simulation procedure

Numerical simulations have been carried out on 700 parallelized cores. Classically, to post-treat and analyze results from instantaneous numerical simulation, calculations on long physical times (> 100 s) are needed in order to extract converged temporal averages of simulated physical quantities. But for this simulation, jointly due to the very high fluctuation induced by the crossing of slugs, the tube height and the slug velocity, physical times of more than 500 s would be needed to achieve converged statistic time average, which would represent months of calculation. Consequently, new post-treatments were developed, focusing on their convergence speed. Every calculation hence followed a three-step scheme, where the first 20 seconds account for a transitory regime, including the initial bed destabilization; a 30 seconds of hydrodynamic transition to complete fluidization; and finally, a 50 seconds calculation allows to extract instantaneous variables for each mesh cell at a 20 Hz frequency. This time delay is certainly sufficient, since the axi-symetric slug velocity is around 0.4 m s -1 and with a frequency of ~0.5 to 1 Hz according to Kong et al. [START_REF] Kong | Bubbling and slugging of Geldart Group A powders in small diameter columns[END_REF]. The time step was adaptive and calculated in order to respect a Courant-Friedrish-Levy and a Fourier condition inferior to 1 and 10 respectively. The convergence criterion was calculated through the sum of every volume Fig. 8. 3D mesh for the numerical simulation 15,000,000 cells.

< i max ( 1 ( ) ) i ncel k s
vol , where i is a given cell within the whole domain ncel and vol is set to 10 8 . The partial differ- ential equations are discretized with a second-order centered scheme and the solution is time-advanced by a first order scheme. The numerical methods used by NEPTUNE_CFD are further described by Hamidouche et al. [START_REF] Hamidouche | Numerical Simulation of Multiphase Reactive Flows[END_REF] and Neau et al. [START_REF] Neau | NEPTUNE_CFD High Parallel Computing Performances for Particle-Laden Reactive Flows[END_REF].

Results and discussion

Experimental results

Slug description

The slug frequency, as determined by high speed camera, varies with the considered level in a bed of given total height, as illustrated in Fig. 9 using the average values obtained by imaging and pressure drop measurements. Combining earlier literature findings [START_REF] Werther | Hydrodynamics and mass transfer between the bubble and emulsion phases in fluidized beds of sand and cracking catalyst[END_REF][START_REF] Yamazaki | Fluidization'85 Science and Technology[END_REF][START_REF] Weimer | On dense phase voidage and bubble size in high pressure fluidized beds of fine powders[END_REF], Kong et al. [START_REF] Kong | Bubbling and slugging of Geldart Group A powders in small diameter columns[END_REF] proposed equations for the slug frequency, f s , the slug length, L s , and the slug intervals, W s . Essential equations are repeated hereafter.

s 1/2 1/2 (5)
The slug frequency is independent of the gas velocity, as confirmed by the experiments, and the constant k must be experimentally determined. Average slug frequencies were determined and were illustrated in Fig. 9.

Frequencies level off for bed levels exceeding 100 to 150 cm, and axi-symmetric slugs start to form at heights exceeding approximately 150 cm. For bed levels above ~1 m, when wall slugs are completely formed, f s is independent of L and about 1 Hz. The data obtained are ~30% higher than proposed earlier by Baeyens and Geldart [START_REF] Baeyens | An investigation into slugging fluidized beds[END_REF]. For axi- symmetric slugs, the equilibrium slug frequency is 0.5 Hz. The slug length is a nearly linear function of U U mf , as illustrated for the ex- perimental measurements in Fig. 10, and corresponds with the prediction of Eq. ( 4). Wall slug lengths of 2 to ~15 cm were imaged at specific levels in the bed and for a given excess air flow rate. Axisymmetric slugs with length of over 15 cm were imaged at bed levels in excess of 1.5 m and higher values of U U mf . The level in the bed clearly determines the slug length at a given U U mf , and slug growth continues along the level in the bed by further coalescence. From the image analysis, the slug velocity was also determined and plotted in Fig. 11.

The bubble or slug velocities increase slightly with excess gas velocity, irrespective of the level in the bed when exceeding 0.8 m. The measured slug velocities correspond with the theoretical prediction of the wall-slugging regime, as indicated in the Figure caption. At lower levels in the bed, wall slugs transform into axi-symmetric slugs with velocities according to the appropriate equation, at higher U U mf values.

When the level in the bed exceeds 1.5 m, slug velocities dropped to ~0.3 to 0.4 m s -1 , corresponding to axi-symmetric slugs, as calculated from:

= + U U U gD ( ) 0.35 s mf (6)

Solids fraction

The bubble fraction , the bed voidage the solids fraction = ) p were calculated from the experimental results and are presented in Figs. 12 and13, respectively. The bubble fraction was calculated from the bed expansion data ( = H 1 /

B mf

). The bed solids fraction was determined from the measured pressure drop across the bed ( = P [ ( 1) ] s

). The bubble fraction in the bed is considerably higher than that for group B powders (expected < 0.1), and the expanded bed is characterized by a fairly low solids fraction [START_REF] Kong | Bubbling and slugging of Geldart Group A powders in small diameter columns[END_REF][START_REF] Baeyens | Bed expansion and the visible bubble flow rate in gas fluidized beds[END_REF]. Both values are important in calculations of e.g. the heat transfer coefficient where a surface renewal model strongly depends upon both p and . Heat transfer measurements and its modeling are the subject of ongoing research. mf cm s -1 . Regimes: freely bubbling, wall slugging and transition to axi-symmetric slugging, axi-symmetric slugging; f ws and f as : expected frequency limits for wall slugging and for axi-symmetric slugging, respectively.

Pressure gradient along the bed height

In order to validate the simulation results, the pressure drop across the bed was compared with the experimental results at U U mf 0.099 m s -1 . It was found to average to solid fraction of = 0.29 s , in fair agreement with experimental results of Fig. 13 ( 0.33 s ), considering the inaccuracy of measuring P across the bed when con- siderable fluctuations occur due to slugging. The slug frequency is shown to be independent of the excess air velocity. Since the volumetric excess gas flow rate equals the product of slug frequency and slug volume, the slug length for a constant frequency and slug diameter ( ~tube diameter), will slightly increase from the onset height to the top of the bed: the hydrostatic pressure indeed decreases with bed height, and the air expands from 1.47 kg m 3 at the bottom of the bed to 1.18 kg m 3 at the top. The air superficial velocity hence increases from × U 25 mf at the distributor to × U 30 mf at the top of the bed, meaning that the slug length will increase by maximum 20% towards the top of the bed.

Bed height fluctuations P P ( )/

In the same way as for experiments, the pressure fluctuations were monitored just above the distributor and the ratio of pressure fluctuation to total pressure drop was determined. This ratio was predicted to be 0.06 against the experimentally measured 0.065 at U 0.099 mf m s -1 (Fig. 5). Pressure fluctuations above 5% meet the criterion of Leva et al. [START_REF] Leva | Fluid Flow Through Packed and Fluidized Systems[END_REF] for the onset of slugging.

Spatially averaged instantaneous solid volume fraction along the bed, z t ( , ) s

As explained before, new post-treatments had to be developed to avoid convergence issues as found with standard physical quantities' horizontal profiles. Thus, the instantaneous solid volume fraction r z t ( , , , ) s was spatially averaged over cylindrical slices defined along the vertical z axis, with their diameter equal to the tube I.D. and with a vertical span of about 9.55 mm. Fig. 14 illustrates the evolution of this spatially averaged solid fraction z t ( , ) s along the tube height, at a given time = t 71.9 s. Similar results have been obtained for several time steps but for clarity only this single simulation time is presented. From the bottom to the top of the tube, two patterns are distinguishable, an initial unstructured small fluctuation corresponding to bubbling flow, followed by stronger isolated peaks highlighting the presence of large bubbles or slugs. In order to filter out bubbling flow, a threshold value, based on a simplified bubble-emulsion model and = d D 0.66 B was set as a criterion for the onset of slugging [START_REF] Kong | Bubbling and slugging of Geldart Group A powders in small diameter columns[END_REF], and defined = 0.22 sslugs as a maximum solid volume fraction for slug detection. The calculation leading to the value of sslugs is detailed in Eq. ( 7), where V slice and h slice account for the volume and the height of the considered slice, respectively.

= + = V h d h d V ( /4) 0.3 /4 0.22 s s slice slice B slice B slice 2 2 slug ( 7
)
Several peaks under the threshold value highlight the presence of slugs after 50 cm bed height and are in good agreement with experi- mental results. A clear fringe pattern evolving periodically along the bed is visible. Considering slugs as a progressive wave propagation along the bed, their vertical velocity was measured graphically with a value of 0.4 m s -1 , in good agreement with experimental results for axi-symmetric slugging. In addition, the fringe pattern only appears after a first random region of 0.3 to 0.5 m at the bottom of the bed corresponding to the experimental height for the onset of slugging.

Temporal Fast Fourier Transform of the spatially averaged solid volume fraction at different heights, t z FFT[ ( )]( ) s

Since slugs can be considered as a periodic wave propagation in the • At the bed bottom (0.1 m), the signal is very random with no dis- tinctive frequency, due to the freely bubbling bed behavior and a wide span of bubble sizes and velocities due to bubble coalescence.

• At 1.0 m height, the amplitude of frequencies below 1.8 Hz is higher than lower in the bed, with a similar superimposed random signal as before: bubbles have coalesced but smaller bubbles continue to coincide with larger ones.

• At 3.0 m height, the high amplitude frequencies zone is now under 1 Hz, with stronger peaks, however with a background signal of higher frequencies due to the presence of small bubbles in between slugs.

After the onset on slugging, a higher amplitude zone emerges due to slugs crossing the studied bed level. Bubbles however remain to form a background signal. In line with the experimental results the frequency tends to decrease between 1 m and 3 m height. Furthermore, the tem- poral and spatial 2D Fast Fourier Transform of the spatially averaged solid fraction was performed leading to similar results as aforementioned. In addition, considering slugging as a progressive wave, their vertical velocity was computed multiplying their frequency by their wave number leading once again to a velocity close to 0.4 m s -1 , in good agreement with experimental results for axi-symmetric slugging.

Final assessment of the NEPTUNE_CFD simulation

NEPTUNE_CFD was found to be able to reproduce accurately the bed expansion (pressure drop along the bed and thus the bed height). In addition, large diluted structures related to slugging were successfully predicted, and in good agreement with experimental results concerning their onset height within the bed. Their rise velocity and frequency was found to be of the same order of magnitude as experimental results for axi-symmetric slugs. No wall slugging, with a higher occurrence frequency and higher rise velocity, was distinguishable. The discrepancy on P P ( )/ (0.06 against 0.065) is also directly linked to the slug shape. Bigger slugs, as found during experiments, would induce larger pressure drop fluctuations when bursting at the top of the bed. The differences and agreements between experiments and simulation results are summarized in Table 3.

Although experimental and simulation results are in fair agreement, additional improvements are required and in order to enhance the slug description in numerical simulation, several parameters were studied. Fede et al. [START_REF] Fede | 3D numerical simulation of a lab-scale pressurized dense fluidized bed focussing on the effect of the particle-particle restitution coefficient and particle-wall boundary conditions[END_REF], Goldschmidt et al. [START_REF] Goldschmidt | Hydrodynamic modelling of dense gas-fluidised beds using the kinetic theory of granular flow: effect of coefficient of restitution on bed dynamics[END_REF] have shown that a decrease of particle-particle coefficient e c trends to increase the for- mation of bubbles and the suspension heterogeneity. In the present study, a decrease of e c from 0.9 to 0.7 allowed to keep the same slug velocity and frequency, but did not significantly improve the slugs characteristics, and trends to decrease slightly the total bed height. solid fraction value along the vertical axis z, the frequency analysis of the spatially averaged solid fraction value was performed at different levels in the bed by computing its Fast Fourier Transform (FFT). Fig. 17 highlights these FFT results at three different levels:

Conclusions

Experimental results and simulations implementing an Eulerian nfluid approach were compared in order to better understand the flow behavior of Geldart group A particles in small I.D. tube within the frame of an Upflow Bubbling Fluidized Bed. Experiments have indicated that the onset height for wall slugs is around 1 m, which is then followed for bed height exceeding 1.5 m by axi-symmetric slugs. The slugging fre- quency was found to be independent of the gas velocity. Due to both the low velocity and length of axi-symmetric slugs, it is of paramount importance to be able to predict, and thus avoid their formation. Indeed, the presence of axi-symmetric slug has a strong impact on the suspension heat transfer coefficient, in particular within the frame of the solar receiver particle-in-tube concept. 3D simulations were performed using ultra fine mesh of more than 15,000,000 cells. Because of characteristic time of slugging phenomenon, classical post-treatments based time averaged physical quantities have been replaced in favor of new posttreatments necessitating less calculation time to converge. Simulation results were found to be in good agreement with experiments concerning the bed average solid volume fraction and successfully reflect the formation of large bubbles structures, hence highlighting the capability of the NEPTUNE_CFD code to simulate such particles flow.

The onset height of slugging was in fairly good agreement with experimental results but simulation was unable to distinct wall slugging from axi-symmetric slugging. In addition, the slug description in numerical simulation results is not as obvious as with experimental results.

Preliminary simulations, carried out with in 2D, highligted the effect of a strong mesh refinement (beyond the minimum cells size for a converged bed height) over slugs prediction. Slugs with size and shape corresponding to experimental observations were predicted for mesh cell size around × d 2 drag for the studied powder. However, this cells size cannot be suitable for industrial geometries and thus the development of subgrid models is essential for an accurate representation of slugging regimes of Geldart group A particles.

Finally, Kong et al. [START_REF] Kong | Bubbling and slugging of Geldart Group A powders in small diameter columns[END_REF] tentatively predicted a delay of slugging to deeper bed heights when operating at higher temperatures. Experiments and simulations at higher temperatures up to 475 °C are in the process of being examined.

Then, various authors have shown that the particles-wall boundary conditions may have an effect o n b ubbles s tructures, b ed h eight and flow i nternal b ehavior [START_REF] Benyahia | Evaluation of boundary conditions used to model dilute, turbulent gas/solids flows in a pipe[END_REF]. R esults d epicted i n t his s tudy h ave been based on a no-slip boundary condition (corresponding to maximum of friction) and elastic particle-wall bouncing. In contrast, the other extreme boundary condition, the free-slip wall boundary conditions (pure elastic frictionless rebounds of spherical particles on a flat w all) has been tested, resulting in a strong diminution of bed height as observed by Fede et al. [START_REF] Fede | 3D numerical simulation of a lab-scale pressurized dense fluidized bed focussing on the effect of the particle-particle restitution coefficient and particle-wall boundary conditions[END_REF], and without any improvement on slug description.

A comprehensive study was carried out on extremely refined meshes, with cells size from . Since 3D simulations with such a fine mesh cell size are currently impossible to be achieved for in an industrial geometry, these tests were performed in 2D over a 2 m height high fluidized bed of 0.03 m I.D., with multiple cartesian grids. The 3D simulation previously performed corresponds to a cell size of about . Within the studied cell size range, the bed height remained constant but the local flow structure of the particle suspension was substantially modified, according to Fig. 18. Only the bottom 1.3 m of the beds has been represented for clarity reasons. From a cell size of , the suspension heterogeneity was emphasised, which allowed for the formation of increasingly larger bubbles structures, but without clearly showing the presence of axi-symmetric slugging. For cell size between however, the suspension heterogeneity kept increasing, showing the formation of slugs with size and shape corresponding to the experimental observations. These extremely refined m eshes e nabled t o c apture m esostructures r esponsible f or slug formation. Finally, the transposition of this mesh cell size to the 3D, 4 m high, experimentally studied geometry, would require a ~4,300,000,000 cell mesh. It is noteworthy that more than just increasing the cell number, the cell size reduction would also imply a time-step duration reduction in order to respect equations resolution stability criteria.

zFig. 1 .

 1 Fig. 1. Bubbling and slugging modes.

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. Photographic illustrations of Axi-symmetric slugs in the powder at the level of 150 cm in the bed: a) = U U 9 mf

Fig. 6 .

 6 Fig. 6. SEM photograph of powder at 150× magnification.

Fig. 7 .

 7 Fig. 7. Cristobalite powder PSD, analysis by Malvern MS2000.

Fig. 9 .

 9 Fig. 9. Bubble and/or slug frequency, f (s -1 ) versus level of measurement in the bed at = U U 12.7

Fig. 10 . 4 . 2 .

 1042 Fig. 10. Illustration of bubble diameters or slug lengths at different excess gas velocities for cristobalite, measured at different bed levels in a bed of = H 300 cm mf

Fig. 11 .Fig. 12 .Fig. 13 .Fig. 15

 11121315 Fig. 11. Slug velocity versus U U mf , measured at different bed levels in a bed of = H 300 cm mf

Fig. 14 .Fig. 15 .

 1415 Fig. 14. Instantaneous profile of the spatially averaged solid fraction along the vertical axis z, at = t 71.9 s.

Fig. 16 .

 16 Fig. 16. Instantaneous map of the spatially averaged solid fraction along the vertical axis function of time with t 0 offsets to 30 s of physical time.

Fig. 17 .

 17 Fig. 17. FFT of the spatially averaged solid fraction z t ( , ) s

Fig. 18 .

 18 Fig. 18. Study on the impact of the mesh refinement onto slugs description: 2D simulations of a 2 m height geometry, instantaneous field of the solid fraction at the same physical time for each case.

Table 1

 1 Instantaneous solid volume fraction

	t 60 s, spatially averaged within a s at = z 50 cm. z 0 to = vertical slice from = nodes number on the tube diameter 60 70 80 90
	spatially averaged solid fraction s	0.282	0.288	0.295	0.297

Table 2

 2 Phases properties used for the numerical simulation (air and powder).

	Gas density at	Gas	Absolute	Particle	Drag	k S	kN	Particle-particle	Frictional	Maximal
	P 101.3 kPa = T 298 K and =	dynamic viscosity	particle density	median volume diameter, d 50	diameter, ddrag			restitution coefficient, e c	viscosity activation threshold, s min ,	particles compaction, s max ,
	kg m -3	Pa s	m -3	μm	μm
	1.1847	1.85 × 10 -5	2640	73.3	45	1.5	4	0.9	0.55	0.64

Table 3

 3 Comparison of experimental and simulation results, at =

				U	25	×	U mf .
				Experimental	Simulation results
				results
	Average solid fraction	0.33	0.29
		P ( )/	P	0.065	0.06
	Onset	wall slugging	0.9 m	>1 m
	height
		axi-symmetric	1.5 m	Between 1 and 3 m
			slugging
	Frequency	wall slugging	1 Hz	No distinction
					possible,
					but < 1 Hz
		axi-symmetric	0.5 Hz
			slugging
	Rise	wall slugging	0.63 m s -1	Not predicted
	velocity	axi-symmetric	0.4 m s -1	0.4 m s -1
			slugging

Acknowledgments

This work was supported by the European Union's Horizon 2020 research and innovation program under Grant Agreement 727762, project Next-CSP. Additional funding was obtained from Beijing Advanced Innovation Center for Soft Matter Science and Engineering of the Beijing University of Chemical Technology. This work was granted access to the HPC resources of CALMIP under the allocation P1132 and of CINES under the allocation gct6938 made by GENCI.