CFD Hydrodynamics Investigations for Optimum Biomass Gasifier Design

Biomass gasification is nowadays considered a viable option for clean energy production. Furthermore, still more efforts need to be spent to make this technology fully available at commercial scale. Drawbacks that greatly limit the full-time plant availability—and so its economically feasibility—mainly concerns syngas purification by contaminants such as tars. Different technological approaches were investigated over last two decades with the aim to increase both the plant availability and the overall efficiency by keeping, at the same time, CAPEX and OPEX low. Among technologies, fluidized beds are surely the most promising architectures for power production at thermal scale above 1 MWth. Gasifier can be surely considered the key component of the whole power plant and its proper design, the main engineering effort. This process involves different engineering aspects: thermo-structural, heat, and mass transfer, and chemical and fluid-dynamic concerns being the most important. In this study, with the aim to reach an optimal reaction chamber design, the hydrodynamics of a bubbling fluidized bed reactor was investigated by using a CFD approach. A Eulerian–Eulerian multiphase model, supported by experimental data, was implemented to describe the interactions between the solid and fluid phases inside the reactor while a discrete dense phase model (DDPM) model was considered to investigate momentum exchange among continuous phases and solid particles simulating char. Different process parameters, such as the bed recirculation rate and the particles circulation time inside the bed, were at least analyzed to characterize the hydrodynamics of the reactor. Results indicate that the recirculation time of bed material is in the order of 6–7 s at bench scale and, respectively, of 15–20 s at full scale. Information about solid particles inside the bed that should be used to avoid elutriation and agglomeration phenomenon, suggest that the dimension of the mother fuel particles should not exceed the value of 5–10 mm.


Introduction
Amongst all energy sources for power generation, renewables had in the last decade the highest rate of growth. IEA (International Energy Agency) projections in its SDS (Sustainable Development Scenario), project bioenergy share by 2030 from by the actual 590 TWh (2019) to about 1170 TWh (2030), with an annual increase of 6% (electricity generation from bioenergy) [1]. Moreover, it is expected that several policies and market initiatives will encourage this trend. In some countries such as China, where use biomass fueled power plant is actually promoted by feed-in tariff, development of biomass fueled power plant, is expected to great extent. The overall goal is to improve air quality to meet environmental policies supporting the phase out of coal-fired boilers and the uncontrolled burning of agricultural residues in the field it a self. As China other countries, such as Brazil and India are implementing similar policies to encourage the development of bioenergy as a way to support that these methods are not yet well understood and their application at commercial scale is actually limited [3,14,15].
Other than the selection of the right operational conditions and the injection of additives (to prevent also the formation of ashes agglomerates) or catalysts during gasification, the approach used in the gasifier design plays a key role in tar formation mechanisms. In the last two decades efforts have been made to investigate the hydrodynamic of fluidized bed reactors, both experimentally and numerically. The main aim is to well understand phenomenon involved during gas-solid interactions as well as to know how these can influence the reactor behavior. It was because of the complexities involved in the process of biomass gasification (fluid dynamic concerns, heat and mass transfer, kinetic of reactions, and so on) that the gasifier design was mainly based on empirical correlations from experimental derivation on bench scale reactors or pilot plant. Furthermore, such empirical formulas had limited applications under certain conditions [16]. Recent developments in numerical techniques and the availability of computational infrastructures, make the application of CFD affordable to investigate design of fluidized bed reactors and their evaluation. In this field two types of methods can mainly be applied: the Eulerian-Eulerian (E-E) approach and the Eulerian-Lagrange (E-L) approach [17]. The E-E approach considers both fluid and particles in the continuous interpenetrating phases. In such a manner granular solid phase is treated as a pseudo-fluid. The E-L approach still considers fluid phase (i.e., gas phase) as continuous but the granular one, is treated by tracking a large number of particles through the computational domain. The E-E approach is less computational expensive compared to E-L approach and is surely the most widely applied in the study of the hydrodynamic of fluidized beds. In the following text, a brief and not exhaustive summary of some contributions in the field is selected with reference to that considered in this study. Kuramoto et al. (1985) [18], developed a new system for circulating fluidized bed within a single vessel where two distinct sections of the reactor were separated by a vertical plate. The two regions were fluidized by a different flow rate so that a circulation of bed material between them was observed. Furthermore, the authors, experimentally investigated the behavior of the system under different operational conditions demonstrating that the proposed configuration potentially had the same advantages of conventional dual-bed system, and that it can be applied to biomass gasification. Kuramoto et al. (1986) [19], experimentally studied the circulation of dense fluidized particles in a two-dimensional bed which was divided in two portions by a vertical plate with an opening. The tests demonstrated how the circulation rate of solid through the opening was controlled by the superficial gas velocity of the up-flowing chamber and by the ratio between the opening area and the cross-sectional area of the down-flowing chamber. Empirical correlations that coupled opening ratio and pressure drop across the opening were also developed. Di Felice et al. (1991) [20], experimentally evaluated the validity of scaling relationships between a pair of dynamically similar gas fluidized beds confirming their application. The study was extended to a wide range of particles: from fine (Geldart group A) to large particles (Geldart group D). The authors highlighted as scaling relationships fail to explain observed slugging regime characteristics. Nieuwland et al. (1994) [21], developed a one-dimensional model for the riser section of a circulating fluidized bed, which describes the steady-state hydrodynamic inside the chamber. Gas and solid phases were considered as two continuous media fully penetrating each other. Zimmermann et al. (2005) [22], investigated the hydrodynamics and kinetics of gas-solid fluidized bed highlighting how the use of both Gidaspow and Syamlal O'Brien drag law overestimate bed expansion in comparison to experimental data. Moreover, they evaluated as the modified Syamlal O'Brien drag law based on the minimum fluidization conditions well predict both the expected bubbling fluidization behavior and bed expansion if compared to experimental evidences. Foscolo et al. (2007) [23], experimentally investigated the behavior of cold model reactor to confirm the feasibility of the pilot scale reactor configuration consisting of two interconnected bubbling fluidized beds. The configuration studied by Foscolo is the landmark of our reactor development. Actually, 1 MW th biomass gasification pilot plant as result of this study, is working, for research purposes, at the ENEA (Italian National Agency for New Technologies, Energy and Sustainable Economic Development) research center Trisaia at Rotondella. Papadikis et al. (2008Papadikis et al. ( -2010 [24][25][26], studied the momentum transport on solid particles, modelling solid fuel, inside bubbling fluidized bed of sand, as continuous. E-E two phases model was implemented to describe the behavior of bubbling bed. Furthermore, kinetic model was implemented to describe gasification reactions inside the reactor. Deza et al. (2009) [27], used a multifluid E-E model to represent the gas and solid phases as continuous. Simulations were carried out to investigate about the behavior of a cold-flow glass bead fluidized bed using as comparison two different drag laws. Numerical data were compared with experimental ones by X-ray flow visualization computed tomography highlighting as Gidaspow model can be used to describe the hydrodynamic of a biomass fluidized bed. Armstrong et al. (2010) [28], used the E-E two-fluid model (TFM) which assumes the gas-solid phases as continuous and fully interpenetrating to evaluate the hydrodynamic behavior of CFB riser during transition from bubbling to circulating fluidized bed. 2D and 3D model were both considered with the implementation of modified drag laws (EMMS). In its study also bubbles diameter was evaluated and compared with the Davinson's model. All evaluations were reported to be in good agreement with experimental data. Gerber et al. (2010) [29], used an E-E multiphase approach for modeling unsteady wood gasification in fluidized bed reactor. Three interpenetrating continuous phases were used to describe the behavior of dispersed phases inside the reactor: one representing wood and two char with different particles diameter. No inert material was considered inside the bed. 2D simulations that also implemented wood pyrolysis, char gasification, and homogeneous gas phase reactions were carried out by using MFIX code by NETL (National Energy Technology Laboratory). Liu et al. (2013) [30], used a comprehensive 3D model by using E-E approach to simulate biomass gasification in a CFB reactor. Turbulent k − ε model was coupled with the kinetic theory of granular flow to simulate the hydrodynamic of the gas-solid system. Kinetics of homogeneous and heterogeneous reactions were also implemented in the whole model to evaluate the effects of same operation parameters on the performances of the CFB reactor. Bidwe et al. (2014) [22], experimentally investigated cold model of dual fluidized bed system (BFB gasifier and CFB regenerator) for solid fuels gasification by using Glicksman's simplified scaling ratios. They showed as parameters under investigation can be useful applied (with respect to scaling laws) in the design of the pilot plant. Canneto et al. (2015) [31], numerically studied bubbles formation and their paths to investigate the behavior of the ENEA's Internally Circulating Bubbling Fluidized Bed (ICBFB) under different fluidization conditions. CFD model by using E-E approach was developed and experimental data from cold model were compared to theoretical ones to study the fluidization quality and to estimate the circulation of solid particles inside the bed. Pecate et al. (2019) [27], experimentally investigate the hydrodynamic of fast internally circulating fluidized bed (FICFB) used for biomass gasification. In its study the influence of bed temperature (between 20 and 950 • C) and the fluidizing agent (air or steam) on the hydrodynamic of a dense fluidized bed of olivine was analyzed. Numerical correlations were used to fit experimental data.
Most of the numerical studies analyzed in the literature are referred to the simple cylindrical geometry of the reactor at bench scale, and many of these are focused only on the evaluation of the main phenomenon involved in the fluidizing process. This work proposes the use of CFD techniques as a new approach to the real rector design by using information achievable by the scaled cold model reactor. In this study the hydrodynamic of the new concept of the ENEA's ICBFB (Internally Circulating Bubbling Fluidized Bed), as revision of the 1 MW th pilot plant [32,33], is under investigation. This configuration merges both advantages of bubbling and internally circulating fluidized beds by keeping all reaction zones in only one vessel. A CFD E-E based model, referred to the real geometry of the scaled cold reactor, was implemented to investigate the behavior of the bed. To describe solid-fluid interactions a modified Syamlal O'Brien drag law is considered while DDPM (discrete dense phase model) was used to discretize solid particles modelling char to evaluate both the interactions with solid and fluid phases.

Materials and Methods
This work deals about the hydrodynamic optimization of the new configuration of the ENEA's ICBFB pilot scale reactor, Figure 1. Within only one vessel, two bubbling fluidized beds of different cross-sectional area are delimitated by a vertical plate at the lower edge of which an interconnection opening between them (fluidized beds) is realized. Different flow rate of fluidizing agents (air) are introduced in the two chambers so that different regimes of fluidization are induced. These conditions create, inside the reactor, two areas at different bed density, which promote the circulation of the inert bed material through the opening from the denser chamber (down-flowing bed, DFB) to the lesser one (up-flowing bed, UFB). By varying the difference between the two superficial velocity inside the two reaction chambers, a different circulation rate of bed material is induced. opening between them (fluidized beds) is realized. Different flow rate of fluidizing agents (air) are introduced in the two chambers so that different regimes of fluidization are induced. These conditions create, inside the reactor, two areas at different bed density, which promote the circulation of the inert bed material through the opening from the denser chamber (down-flowing bed, DFB) to the lesser one (up-flowing bed, UFB). By varying the difference between the two superficial velocity inside the two reaction chambers, a different circulation rate of bed material is induced.

Cold Model Definition
Numerical investigations were carried out on bench scale cold model reactor in fluid dynamic similar to the real ICBFB pilot plant. For this purpose, scale transposition laws were assumed to be in agreement with some dimensionless parameters (nine dimensional groups) as proposed by Glicksman [34]. In such a manner phenomenon involved and the hydrodynamic behavior observed on bench scaled cold model can directly be transposed to real reactor. Different applications of this criterion as carried out by Kolar and Leckner [35], where a 12 MWth CFB scale transportation is involved, demonstrate how it is very difficult to meet at the same time, all nine dimensionless groups. In this study, the methodology proposed by Foscolo et al. [23], where only four of these parameters must be satisfied, was considered: (1) Flow number:

Cold Model Definition
Numerical investigations were carried out on bench scale cold model reactor in fluid dynamic similar to the real ICBFB pilot plant. For this purpose, scale transposition laws were assumed to be in agreement with some dimensionless parameters (nine dimensional groups) as proposed by Glicksman [34]. In such a manner phenomenon involved and the hydrodynamic behavior observed on bench scaled cold model can directly be transposed to real reactor. Different applications of this criterion as carried out by Kolar and Leckner [35], where a 12 MW th CFB scale transportation is involved, demonstrate how it is very difficult to meet at the same time, all nine dimensionless groups. In this study, the methodology proposed by Foscolo et al. [23], where only four of these parameters must be satisfied, was considered: Archimedes number : Ar = d p 3 ·ρ g · ρ p − ρ g ·g (1) Flow number : Length number : with respect to bench scale tests, with known fluidizing flow (i.e., air), well defined temperature, pressure (i.e., reference conditions), and geometry of the real reactor, thus satisfying the above written dimensionless parameters, the geometry of the cold model is fully determined, as well. In addition to such parameters, also the Reynolds and Froude numbers (even if their combination defines the Archimedes number) were verified as, additional conditions: Froude number : Main results of the scaling procedure for the cold model definition are summarized in Table 1. As shown, bed material with particle density of the order of ρ p 8900 kg/m 3 and with particle diameter of about d p 0.14 mm must be chosen to meet all scaling parameters. In addition, the satisfaction of the density number criterion imposes to use as fluidizing flow a medium with density of ρ g 0.99 kg/m 3 at 25 • C. All conditions were satisfied by considering copper powder (ρ p = 8930 kg/m 3 , d p mean = 0.137 mm) as bed material and a mixture of air and helium at 3% wt. at 25 • C as fluidizing flow. In conclusion, bench scale reactor, operated in cold flow conditions (25 • C, 1 bar), must be in length scale ratio with the real gasifier of 3.65 (scale 1:3.65) and should use copper powder as bed material and a mixture of air and helium as fluidizing flow.

Mathematical Model
Eulerian-Eulerian multiphase approach was considered to describe gaseous flow (fluidizing flow as primary phase) and granular solid material (bed material as secondary phase) behavior. As defined earlier in the introduction, all interacting phase are considered, by E-E approach, i.e., as continuous and interpenetrating. The volume fraction occupied by each one phase cannot be occupied by another one: where n is the number of interacting phases, ε q the volume fraction of phase q. In addition, each phase must satisfy conservation equations: both continuity and momentum. The first of these, for a generic phase q, is expressed by m ij the mass transfer rate between phases i and j, S q the source term for phase q (generally equal to zero). If close system, with no mass transfer between phases is considered, the previous equation is simplified as Under the same assumptions and additionally when volume, lift and virtual mass induced forces are negligible, the momentum conservation can be expressed, for gas phase, by where m is the number of solid phases that interact with the gaseous one. The first term at left-hand side accounts for unsteady acceleration terms while the second one for the acceleration related to convective flow. At right-hand side, the term ε g ∇p accounts for pressure variation, while ∇· = τ g represents viscous forces. In this latter, the stress tensor is given by where µ g and λ g are respectively the cinematic and bulk viscosity of gas flow and = I the identity tensor.
The term m s=1 → R sg in Equation (10) accounts for gas-solid interacting forces and is given by where K sg is the interphase exchange coefficient, while → v s − → v g is the relative solid to gas velocity. Conservation of momentum must be satisfied also by solid phase: where the term ∇p s accounts for solid pressure forces, while the n p=1 → R ps represents the interaction forces between the gas or solid phase p and the solid phase s. The interphase exchange coefficient is given by where the parameter f provide for the drag coefficient C D as function of the relative Reynolds number Re s : Different expressions of the drag coefficient are available in the technical literature. The most suitable for fluidizing systems are those as proposed by Wen-Yu [36], Syamlal-O'Brien [37], and Gidaspow [38] whose results seem to be in good agreement with experimental ones with flow conditions exceeding the minimum fluidization [39]. Furthermore, the Gidaspow model, that is a combination of f Wen-Yu and Ergun, is the most widely used in the field of fluidizing flow concerning bubbling fluidized bed [25]. In this study the parametrized model of Syamlal-O'Brien was used to describe the interactions between gas and solid phases [40]: The terminal velocity for the solid phase is determined as log 10 c 1 log 10 ·0.85 (18) Here, particles diameter is defined and minimum fluidization velocity is known experimentally. Parameters c 1 and d 1 are evaluated by solving iteratively Equation (18) by following the numerical procedure proposed by Syamlal et al. [40].
Solution of Equation (13) imposes the use of a closure model for the term p s . Different approaches are proposed in the literature though the widely used is that of Lun et al. [41]: where Θ s is the granular temperature, e ss the coefficient of restitution and g 0,ss is the radial distribution function that accounts for the probability of collision between solid particles. The stress tensor = τ s in Equation (13), has the same formulation of that one defined for gas phase = τ g by Equation (11). It contains the granular viscosity µ s and the granular bulk viscosity λ s Terms originated by momentum exchanging in solid particles collisions and translations. The first of these can be evaluated as sum of three viscosity terms: The collisional µ s,col , the kinetic µ s,kin and the frictional µ s, f r by Instead, the granular bulk viscosity accounts for frictional phenomenon originated by the compression and the expansion of the solid granular phase: Numerical coupling between continuous phases and discrete one which model single solid particles injected into the bed, was realized by considering a DDPM (dense discrete particle model) which solves the trajectory of a single particle by mean the equation: where f is the drag factor and τ u a response time parameter: Different numerical implementations are available for the drag factor evaluation. In this study that of Putnman as proposed by Papadikis et al. [24], as given below, are used: The third term, F vm , at right hand side of Equation (22) considers unsteady effects related to virtual mass forces: Exhaustive details about mathematical model can be found in the technical literature [42].

Computational Methodology
Ansys Fluent code was used to formulate the governing equations of the problem adopting a finite volume approach. PC-SIMPLE (phase-coupled semi-implicit method for pressure linked equation) algorithm with a first order formulation of convective terms, was selected to carry out transient simulations by setting 10 −3 s as time step. Overall simulation time of 20 s was imposed at each run for a total number of 20.000-time steps (each one with 100 iterations). Two UDFs (user defined function) were implemented in the original code to introduce a modified form of both the interphase exchange coefficient K sg , Equation (14), and the drag factor f p , Equation (24). The first one parameter, was determined by implementing the numerical procedure proposed by Syamlal et al. [40] when minimum experimental fluidization velocity was known. The second term was instead defined as proposed by Papadikis et al. [24]. Both 2D and 3D numerical grids with different refinement grade were considered during calculations. Initial solid volume fraction as well as initial field velocity, were imposed both for freeboard and granular bed region. Constant velocity profiles were imposed at inlet sections (velocity inlet), while atmospheric pressure was imposed at outlet section (pressure outlet) of the numerical grid. Gaseous phase velocity at walls was imposed to be equal to zero (no-slip boundary conditions) while no tangential stress conditions (free-slip boundary conditions) was imposed for the solid phase.

Results and Discussion
Validation of mathematical model is of crucial importance in order to achieve good accuracy in results with respect to physical observations. The importance is even greater if obtained results are to be used as guideline in the designing of equipment. The effective hydrodynamic conditions inside a fluidized bed reactor, is very difficult to capture. In view of the fact that experimental measure methodologies do not have to perturbate fluid dynamic conditions inside the reaction chamber, in last two decades or so, well known non-invasive techniques, are assuming a great importance. The tomography and fluoroscopy with X-ray are well known techniques widely used as allowing measurements of the phases distribution inside the reactor, by means of the implementation of image analysis algorithm (DIA-digital image analysis).

Numerical Benchmark
In this study, mathematical model was experimentally validated by using data available by Deza et al. [43] which conducted tests on bench scale cold model by the implementation of X-ray methodology. A schematic representation of the apparatus used during tests is shown in Figure 2. In the same figure, characterized properties of particles used as bed material, are also shown. Experiments were carried out by using an acrylic cylindric shaped reactor whose diameter is equal to 9.5 cm. Glass beads were used as bed material (with an initial bed height of 10 cm) and air as fluidizing agent by imposing a superficial gas velocity of U g = 1.3 ·U m f . It must be highlighted that bed material used by Deza et al. in their study was different from that identified for the cold model considered here, where copper powder was selected. This methodological choice is justifiable in view of the lack of available experimental data for copper powder. Furthermore, it must be stated that both glass and copper particles are classifiable as Geldart Group B, so that the phenomenological nature of the observations should be the same. This means that mathematical model, parametrized with respect to the different particles used and describing the two processes, should have the same response.  Preliminary experiments were conducted on the ENEA's ICBFB pseudo-2D cold model, Figure  3; PIV (Particle Image Velocimetry) methodology was adopted to investigate the velocity field structure, seem to confirm these assumptions. A schematic representation of the ICBFB pseudo-2D cold model is shown in Figure 7. Reactor, 10 mm depth in perpendicular direction to the plane, was equipped with a front panel realized in polycarbonate in order to investigate solid motion inside the two chambers. A black back panel was instead used in order to create the right light contrast for the acquisition. The acquisition area was lighted by using two light led sources while a high speed camera, provided by La Vision, was used for frame acquisitions. Cross-correlation method was used to correlate frames at different time in order to solve the solid velocity field. A comparison of two different images of tests conducted, under the same fluidizing conditions, with 0.14 mm copper powder and 0.55 mm glass beads, are shown in Figure 3. The comparison should be regarded as qualitative. Of course, the two experiments were carried out in the same fluidizing conditions but Preliminary experiments were conducted on the ENEA's ICBFB pseudo-2D cold model, Figure 3; PIV (Particle Image Velocimetry) methodology was adopted to investigate the velocity field structure, seem to confirm these assumptions. A schematic representation of the ICBFB pseudo-2D cold model is shown in Figure 7. Reactor, 10 mm depth in perpendicular direction to the plane, was equipped with a front panel realized in polycarbonate in order to investigate solid motion inside the two chambers. A black back panel was instead used in order to create the right light contrast for the acquisition. The acquisition area was lighted by using two light led sources while a high speed camera, provided by La Vision, was used for frame acquisitions. Cross-correlation method was used to correlate frames at different time in order to solve the solid velocity field. A comparison of two different images of tests conducted, under the same fluidizing conditions, with 0.14 mm copper powder and 0.55 mm glass beads, are shown in Figure 3. The comparison should be regarded as qualitative. Of course, the two experiments were carried out in the same fluidizing conditions but with different density number, in two cases. Nevertheless, the behavior of the two systems appear to be the same. Dimension of the bubbles is in the same order of magnitude. Moreover, bed expansion and up-flowing chamber is visibly less dense than down-flowing one. Furthermore, circulation of bed material through the opening between the two chambers is well visible in the two cases. Under these assumptions, three different grid refinement, by considering as limiting maximum cell size s c < 10·d p , were considered in the ratio 1:0.5:0.125 (9.10·d p , 4.55·d p and 2.27·d p ) with respect to the coarsest mesh. In such a manner 19 × 80, 38 × 160 and 76 × 320 grid patterns were obtained. Simulations were carried out by imposing for each run, a global simulation time of 40 s with a temporal discretization of 10 −3 s. Only results between 5 s and 40 s were considered by averaging values every 0.01 s. With the aim to validate computational code, the same mathematical model was also solved by using Mfix code by NETL. A User Defined Function (UDF) was used to modify the original computational code and to implement the parametric model of Syamlal-O'Brien for the calculation of the interphase exchange coefficient. A comparison between theoretical (Ergun equation), experimental and numerical (CFD) pressure drop trough the bed as function of the gas superficial velocity, is shown in Figure 4. Additionally, collected numerical parameters for the UDF definition c 1 and d 1 , calculated by using the procedure proposed by Syamlal et al. [40], are shown in Figure 4.     As it is possible to see, numerical case labelled as Fluent CFD 76 × 320 Mesh-UDF, implementing UDF, fits very well theoretical (red dashed line) and experimental data (red circle). The implemented model provides a good estimation of both minimum fluidization velocity and the correlated pressure drop. Only a slight underestimation and overestimation respect to experimental data (below and above the minimum fluidization velocity value), were registered. Other numerical cases (Fluent without UDF and MFix) only well predict pressure drop above the minimum fluidization conditions. A comparison between CFD solid phase volume fraction obtained by calculations carried out by using different grid refinement, is shown in Figure 5. Images are referred to calculations at simulation time t = 20 s by using Fluent code implementing UDF. As noticeable, coarsest mesh (19 × 80) does not allows to describe with sufficient details, spatial distribution of bubbles as well as their dimensions and borders. Only bed expansion is sufficiently predicted. On the contrary finer meshes (both 38 × 160 and 76 × 320), well describe, with about the same level of detail, bubbles distribution inside the bed, their borders and dimensions: bubble size at bed surface was estimated to be in the order of 2-4 cm whose values are in good agreement with theoretical 0-D solutions (calculations carried out by using correlations whose results are independent from model dimensions). Furthermore, in the case of 76 × 320 grid refinement, visual accordance with X-ray radiography images, obtained by Deza et al. [44] at the same time, is very good.
A comparison between experimental and numerical (CFD) profiles of the void of fraction at different height above the flow distribution plate (z = 4 cm and z = 8 cm), as function of column diameter and column height, is proposed in Figure 6. In the diagrams, experimental data are collected with reference to two perpendicular directions (X-slice and Y-slice). As is evident, numerical and experimental trends are comparable provided finest meshes are considered for both codes here used (Fluent and Mfix). Different proposed profiles in Figure 6a are different in both trend and range. Most probably, the differences between the two groups of curves (theoretical and experimental) are imputable to the real reactor set-up. In particular, the large oscillation in values at z = 4 cm is related to the presence of the flow distribution plate that, in the real reactor, has a finite number of holes through which fluidizing flow is introduced. Furthermore, a channeling phenomenon near the reactor walls is evident so that a greater range of values in the pressure profile at the base of reactor, is registered. 160 and 76 × 320), well describe, with about the same level of detail, bubbles distribution inside the bed, their borders and dimensions: bubble size at bed surface was estimated to be in the order of 2-4 cm whose values are in good agreement with theoretical 0-D solutions (calculations carried out by using correlations whose results are independent from model dimensions). Furthermore, in the case of 76 × 320 grid refinement, visual accordance with X-ray radiography images, obtained by Deza et al. [44] at the same time, is very good.

× 80
38 × 160 76 × 320 (a) (b) (c) A comparison between experimental and numerical (CFD) profiles of the void of fraction at different height above the flow distribution plate (z = 4 cm and z = 8 cm), as function of column diameter and column height, is proposed in Figure 6. In the diagrams, experimental data are collected with reference to two perpendicular directions (X-slice and Y-slice). As is evident, numerical and experimental trends are comparable provided finest meshes are considered for both codes here used (Fluent and Mfix). Different proposed profiles in Figure 6a are different in both trend and range. Most probably, the differences between the two groups of curves (theoretical and experimental) are imputable to the real reactor set-up. In particular, the large oscillation in values at z = 4 cm is related to the presence of the flow distribution plate that, in the real reactor, has a finite number of holes through which fluidizing flow is introduced. Furthermore, a channeling phenomenon near the  All this promote a not uniform distribution of the flow inside the chamber with the formation of large bubbles where local conditions are favorable. This is not evident in numerical case (CFD) where a uniform distribution of flow at the reactor base was imposed. The observed phenomenon is less evident at z = 8 cm, Figure 6b, because of the greater distance from the flow distributor. Noticeable is, in all cases considered, the slightly differences with respect to experimental X-slice curves even if 3-D solutions are considered, Figure 6c.
Nevertheless, the maximum percentage error is less than 3.5% (coarsest Mfix mesh). The independence of results toward grid refinement was evaluated following ASME (American Society of Mechanical Engineers) procedure for estimation of the uncertainty due to discretization in CFD applications [45]. GCI procedure (grid convergence method), based on the Richardson extrapolation method, was used to evaluate the selected variables (for the purpose: pressure drop through the bed and mean void fraction at z = 4 cm and z = 8 cm) do not exhibit a monotone dependence toward grid refinement.

ENEA's ICBFB Cold Model Analysis
The validated numerical code was then applied to study the behavior of the cold model describing the ENEA's ICBFB. A schematic representation of the cold model considered here is shown in Figure 7. Furthermore, the main properties of the material used (copper powder) for describing the inert fluidizing bed and fluidizing conditions imposed at the inlet section of the two flowing chambers, are shown in the same figure.
Fluid dynamics similarity conditions (Table 1) impose to use as bed material, copper powder with mean particle diameter equal to d p = 137 µm. Of such particles, no fluidizing properties are experimentally known. For the purpose, minimum fluidization conditions were determined by using the correlation proposed by Coltters et al. [46]: (26) whose values are in good agreement with experimental determinations (R = 0.990) and suitable for particles with diameter 3 µm ≤ d p ≤ 900 µm and density 2.7 g cm 3 ≤ ρ p ≤ 11.3 g cm 3 .
The above correlation was successfully tested with experimental data by Chiba et al. [47] considering reference values for copper shot of diameter d p = 163 µm. CFD calibration curve implementing Syamlal-O'Brien parameters for the considered particles is shown in Figure 8. As is possible to see, CFD numerical results are in good accordance with theoretical ones and accurately predict both minimum fluidization velocity (9.67 cm/s) and pressure drop at incipient fluidization conditions (5448.57 Pa). and mean void fraction at z = 4 cm and z = 8 cm) do not exhibit a monotone dependence toward grid refinement.

ENEA's ICBFB Cold Model Analysis
The validated numerical code was then applied to study the behavior of the cold model describing the ENEA's ICBFB. A schematic representation of the cold model considered here is shown in Figure 7. Furthermore, the main properties of the material used (copper powder) for describing the inert fluidizing bed and fluidizing conditions imposed at the inlet section of the two flowing chambers, are shown in the same figure.  Fluid dynamics similarity conditions (Table 1) impose to use as bed material, copper powder with mean particle diameter equal to  Main CFD numerical results for the analyzed case are reported in Table 2 where other than main bed properties, also calculated mass flow rate through the opening is proposed. The proposed numerical results are related to the best grid refinement at solution convergence. This latter was evaluated by following the ASME procedure [45] based on the Richardson extrapolation method. In the same table, also theoretical mass flow rate through the opening, evaluated as proposed by Kuramoto et al. [19], were considered for comparison purpose. In this way, two different correlations were used to evaluate numerically the mass flux through the opening connecting the two different regions of the bubbling fluidized bed. The first one, proposed by Jones and Devidson [19], is suitable Main CFD numerical results for the analyzed case are reported in Table 2 where other than main bed properties, also calculated mass flow rate through the opening is proposed. The proposed numerical results are related to the best grid refinement at solution convergence. This latter was evaluated by following the ASME procedure [45] based on the Richardson extrapolation method.

Copper Powder
In the same table, also theoretical mass flow rate through the opening, evaluated as proposed by Kuramoto et al. [19], were considered for comparison purpose. In this way, two different correlations were used to evaluate numerically the mass flux through the opening connecting the two different regions of the bubbling fluidized bed. The first one, proposed by Jones and Devidson [19], is suitable for small opening ratio where S 0 is the opening section, S d the downflowing chamber cross sectional area, ρ d the DFB bed density, ∆P 0 the pressure drop across the opening, and C s a particle discharge coefficient ranging between 0.5 and 0.65. The second correlation here used, was that proposed by Kuramoto et al. [19]: where C s0 is the discharge coefficient evaluated at ε u = 1 and equal to 0.5. It must be highlighted that in the above correlations, W d is the mass flux [kg/m 2 s] calculated at the mid-section of the downflowing chamber, so that its value is greatly influenced by the real geometry of the reactor. This condition explains the slightly difference between the calculated values and CFD one. By analyzing numerical results, it is clearly evident as the two regions inside the reactor (UFB and DFB), are characterized by a different behavior. This is well noticeable by evaluating for each of these the bed expansions ratio, the bubbles fraction and the mean void fraction terms. DFB chamber is denser than UFB one, that also exhibits a greater bed expansion respect to the initial condition. Furthermore, a greater fraction of bubbles was observed in the UFB region, where also a positive value (in direction + y) of the mean solid velocity was registered.
On the contrary DFB region is characterized by a negative (in direction-y) value of the same term of velocity. Again, the flux vector evaluated on the opening plane, exhibits as its main component that one in perpendicular direction to the plane along the positive x axes. All these conditions clearly are explained by a motion of the bed material between the two regions inside the reactor from the denser to less dense chamber across the opening at the bottom side of the vertical plate. All the above considerations are well visible in Figure 9, where a comparison between results obtained with the coarsest and finest mesh is proposed for different simulation time. Images show, on the background, the solid phase volume fraction and, in overlapping, the solid phase velocity field. Each of these, clearly shows the bed motion between the two reactor regions. Two different macro-motions are well discernable: the first one (MM1) in counter clockwise direction, that developing along the reactor walls, involves both chambers and, the second one (MM2) concerning only DFB chamber. This latter is composed by two sub-motions: the first one (SM1) identifiable at the bottom left side of the DFB chamber near the flow distributor and, the second (SM2) well visible at the upper right side near the vertical plate. SM1 develops in counter clockwise direction and involves the left side region of the DFB chamber. The closed pathways outlined by this could be responsible of the formation of stagnation regions where reacting material (e.g., char particles) or residual ashes could be trapped causing hot spots or agglomerates. When this motion involves the upper bed surface (left upper area), finer particles like ashes mixed to fine char, could be elutriated and trapped near the reactor wall. Both these are dangerous conditions that must be avoided in order to not compromise the operating functionality of the reactor. SM2, instead, involves the right region of the DFB chamber in clockwise direction. It develops from the bottom surface of the chamber toward the upper bed surface where its motion is inverted. A descending solid flow near the vertical plate wall, is originated that drags solid material toward the opening at the bottom side of the plate. Here solid materials are mixed to that from MM1, that crossing the opening, enters the UFB chamber.
Entering materials are then vigorously transported by the ascending gas flow toward the surface of the bed where they are reintroduced in the down flowing chamber. Comparing bed structures resulting from the solution of the mathematical model with the two different grid refinements, it is noticeable as the coarsest mesh can't capture with good level of detail smallest formation inside the bed. Moreover, bubbles distribution and their dimensions cannot be evaluated with sufficient grade of precision. Only bed expansion and main solid flow macro-motions are sufficiently predicted by coarse mesh. A comparison between results obtained by using the two levels of refinement are shown in Figure 9. In the proposed diagrams, velocity and void fraction profiles are compared. Circulation of the bed material between the two zones of the reactor is also confirmed by the velocity profiles, Figure 10. As shown UFB profiles exhibit mean positive values of velocity at the contrary of DFB ones whose mean values are negative. By the same diagram also the inversion of solid flows near the reactor walls are visible. Furthermore, UFB velocity profiles are shifted toward the wall of the vertical plate where the formation of bubbles of greater dimensions is promoted. This latter condition is well confirmed by the analysis of the distribution of the void fraction along the chamber, Figure 11a,c. At least, comparing void fraction profiles as function of the bed height, Figure 11b,c, it can be observed that UFB chamber exhibits a greater bed expansion than DFB one and how this latter is denser than the first one. In brief, comparing results, it can be stated that coarsest mesh over predict velocity of about 45% and void fraction of about 4% are well compared with fine mesh (values at convergence). A visual representation of 3D-CFD solution at simulation time t = 10 s, is proposed in Figure 12, where solid phase volume fraction and solid velocity field on main section planes of the reactor are shown. X-Y slicing plane, shows well the distribution of the field velocity in the two chambers, Figure 12a, while granular solid phase iso-surfaces depicted in 12b, delimitate bubbles formations inside the reactor. This latter image confirms general phenomenon discussed above. Bubbles of dimensions in the order of 4-6 cm, mainly develop in UFB chamber near the wall surfaces of the vertical plate by-passing, in this region, reacting flows. The evaluated mass flow rate across the opening, in 3D-CFD simulation, equal to W x dA O = 6.65 kg/s, is comparable both with 2D and theoretical results proposed in Table 2. In order to investigate the recirculation time of solids inside the reactor, char particles by biomass gasification, were considered. Experimentally measured density of char particles was evaluated to be equal to ρ p = 975 kg/m 3 . Five values of diameter ranging from 0.2 mm to 1.0 mm were considered as average diameter of char particles in the bed, as result of primary fragmentation and particle shrinkage. These values can be correlated to the original fuel size (mother fuel particle) as proposed by Gómez et al. [48].  Figure 9. Solid phase volume fraction and field velocity comparison between CFD solutions with different grid refinement and at various simulation time. Figure 9. Solid phase volume fraction and field velocity comparison between CFD solutions with different grid refinement and at various simulation time.
in Figure 9. In the proposed diagrams, velocity and void fraction profiles are compared. Circulation of the bed material between the two zones of the reactor is also confirmed by the velocity profiles, Figure 10. As shown UFB profiles exhibit mean positive values of velocity at the contrary of DFB ones whose mean values are negative. By the same diagram also the inversion of solid flows near the reactor walls are visible. Furthermore, UFB velocity profiles are shifted toward the wall of the vertical plate where the formation of bubbles of greater dimensions is promoted. This latter condition is well confirmed by the analysis of the distribution of the void fraction along the chamber, Figure 11a,c. At least, comparing void fraction profiles as function of the bed height, Figure 11b,c, it can be observed that UFB chamber exhibits a greater bed expansion than DFB one and how this latter is denser than the first one. In brief, comparing results, it can be stated that coarsest mesh over predict velocity of about 45% and void fraction of about 4% are well compared with fine mesh (values at convergence). A visual representation of 3D-CFD solution at simulation time t = 10 s, is proposed in Figure 12, where solid phase volume fraction and solid velocity field on main section planes of the reactor are shown. X-Y slicing plane, shows well the distribution of the field velocity in the two chambers, Figure 12a, while granular solid phase iso-surfaces depicted in 12b, delimitate bubbles formations inside the reactor. This latter image confirms general phenomenon discussed above. Bubbles of dimensions in the order of 4-6 cm, mainly develop in UFB chamber near the wall surfaces of the vertical plate by-passing, in this region, reacting flows. The evaluated mass flow rate across the opening, in 3D-CFD simulation, equal to ∮ W = 6.65 / , is comparable both with 2D and theoretical results proposed in Table 2. In order to investigate the recirculation time of solids inside the reactor, char particles by biomass gasification, were considered. Experimentally measured density of char particles was evaluated to be equal to = 975 / . Five values of diameter ranging from 0.2 mm to 1.0 mm were considered as average diameter of char particles in the bed, as result of primary fragmentation and particle shrinkage. These values can be correlated to the original fuel size (mother fuel particle) as proposed by Gómez et al. [48].  As proposed by the authors, the initial fuel size , is related to the average size of char particle inside the bed, , , by the following correlation: As proposed by the authors, the initial fuel size d f , is related to the average size of char particle inside the bed, d ch,b , by the following correlation: where ϕ is the shrinkage factor, n 1 the number of fragments, n 2,m and σ multiplication factors determined, experimentally. Experimental determinations proposed by the same authors indicate for wood chips a value of d ch,b d f = 0.34, which was assumed here as reference for numerical evaluations. Table 3 summarizes the scaling parameters used to study the behavior of char particles inside the cold model reactor. The table present data related to char particles of size 0.80 mm and 1.00 mm, respectively. The only ones whose results were considered of interest. 3D numerical simulations were carried out by setting up a DDPM as above discussed. Five equivalent particles were injected at simulation time t = 10 s when inert bed material was well fluidized. Simulations were terminated at simulation time t = 20 s and data about particles (e.g., position, velocity, and momentum) were collected every 0.01 s. Results showed that particles whose diameter is lesser than 0.22 mm (equivalent to 0.80 mm char particles), segregate on bed surface in the upper left region of the bed, near the reactor walls, where stagnant fluid dynamic conditions are encountered. Instead, particles whose diameter is in the range 0.22 mm to 0.27 mm (equivalent to 1.00 mm char particles) well follow bed motion, circulating between the two reaction chambers.
The calculated recirculation time of average particles is in the order of 5 s for 0.22 mm particle, and 6-7 s for 0.27 mm particle, well comparable with the average bed recirculation time estimated to be in the order of 6-7 s. 2D-images of particle motion inside the reactor bed are shown in Figure 13, where white and black dots are used to identify respectively 0.22 mm and 0.27 mm particles. As shown particle of averaged size 0.22 mm, follows the main bed macro-motion, circulating continuously between the chambers, while 0.27 mm particle, seems to be trapped in the bottom section of the UFB chamber where closed pathway solid flow lines are encountered. All these conditions suggest, in last instance, the dimension of fuel particles that must be used in the real reactor to avoid both segregation and agglomeration phenomenon. As discussed above, by using the correlations proposed by Gómez et al. [48], these shouldn't exceed the estimated value of 5-10 mm.  Furthermore, information about solid particles inside the bed that should be used to avoid elutriation and agglomeration phenomenon are achieved. Particles whose dimensions are less than Figure 13. Visual representation of equivalent char particles position inside cold model reactor at different simulation time: white dot for particle of diameter d p = 0.22 mm; black dot for particle of diameter d p = 0.27 mm. Images also show the instantaneous granular solid phase volume fraction and the solid velocity field.

Conclusions
In this paper, the hydrodynamic of cold model reactor was investigated numerically in order to get guidelines in the design procedure of a real reactor. CFD method was applied by treating continuous both granular solid phase and gaseous phase. Eulerian-Eulerian multiphase approach was then implemented to describe the two interacting phases while a DDPM was considered to evaluate interactions between continuous phases and solid particles modelling char.
Mathematical model was validated first by using experimental data available in the technical literature and then applied to the numerical case under investigation. The geometry of the ICBFB cold model reactor and the fluidizing conditions were fully defined by applying scaling laws to preserve fluid dynamic similarities between the real pilot reactor at full scale (1 MW th ), operating at high temperature, and that one at bench scale operating at ambient conditions. Behavior of the inert bed material was investigated under the design fluidizing conditions and critical conditions (e.g., elutriation and agglomeration phenomenon) were identified.
Recirculation rate of bed material between the two reaction chambers was also evaluated and compared with theoretical one calculated by using semi-empirical correlations by different authors. The behavior of solid particles modelling char was also investigated by injecting solid particles of different size inside the well fluidized bed and by keeping for each one all information (e.g., position, velocity, and momentum) at different simulation time.
Results indicate that the recirculation time of bed material is in the order of 6-7 s at bench scale and, respectively, of 15-20 s at full scale.
Furthermore, information about solid particles inside the bed that should be used to avoid elutriation and agglomeration phenomenon are achieved. Particles whose dimensions are less than 0.22 mm (0.80 mm at full scale reactor) should segregate on bed surface agglomerating near the reactor walls, while particles greater than 0.27 mm (1.00 mm at full scale reactor) should migrate to the bottom section of the UFB chamber were stagnant fluid dynamic conditions were detected. Particles of dimensions ranging between 0.22 mm and 0.27 mm, instead, well follow the bed motion circulating between the two reaction chambers.
All these information suggest that the dimension of the mother fuel particles, to feed to the real pilot reactor, should not exceed the evaluated value of 5-10 mm.
Results obtained should be considered as preliminary. Further investigations need to be done in order to validate numerical CFD results, experimentally. Developments of the work are devoted to experimental investigations on the ENEA's pseuo-2D cold model where PIV methodology is implemented to solve the solid velocity field. Moreover, the effects of fluid dynamic implications related to the real bed mixture when both solid fuel particles and inert bed material are considered, need to be investigated further. From numerical point of view, implementation of the chemistry with related kinetics of reaction is the next step.
Funding: Research activities were funded by H2020 Spring G2E project supported by the Italian ministry of economic development.