CFD simulations to study bed characteristics in gas–Solid fluidized beds with binary mixtures of Geldart-B particles: A qualitative analysis

The bed dynamics of unary and binary fluidized beds play a key role in understanding the pressure drop and hence provides an opportunity for performance improvement of the beds. In the present work, characteristics of fluidized beds with binary mixtures of Geldart-B particles were investigated using CFD simulations. The phenomena of segregation and mixing using simulations were studied, both qualitatively and quantitatively, at a range of superficial gas velocities (0.3–0.6 m/s) and two different bed heights. The study was divided into two parts. In Part I, the current study, a qualitative analysis of flow patterns for seven different binary mixtures, is presented. The quantitative analysis, including particle and gas velocity profiles, particle volume fraction profiles, and correlations for minimum fluidization velocity and pressure drop, will be presented in Part II of this work. A mathematical model consisting of an Eulerian-Eulerian model with RNG k-ε model and KTGF model to capture the bubble dynamics was used. The standardized values of coefficients and plastic stresses have been used for all simulations. The CFD model was validated using experimental data from the literature. Qualitative predictions of volume fraction profiles of small-sized particles showed that, for mixtures within a range of 40%–60% Geldart-B type large particles, the bubble and solid particle dynamics were different from those of single particles of the superficial gas velocities considered. In contrast to the single particles in the given superficial gas velocity range that were in bubbling regime, the binary particles showed a transition from bubbling to slugging to turbulent regime, as demonstrated by qualitative analysis. A homogeneous regime was observed for lower superficial gas velocities for mixtures consisting of 0%–20% large particles.


Introduction
Gas-solid (GS) fluidized beds have been important in various applications, such as drying (Yohana et al., 2020), granulation (Behzadi et al., 2009), blending, combustion, gasification (Roy et al., 2021), and conversion of methanol to olefins (Chang et al., 2019), for more than seven decades. One of the several advantages of fluidized beds over fixed beds is their ability to be operated isothermally, with minimal axial temperature gradients (Menéndez et al., 2019). Analytical models of fluidized beds are complex, involving dynamics and transport phenomena of two or more phases; for example, gas and solids/particles and bed phases such as bubble, cloud, and emulsion. The first pioneering work in analytical modeling of fluidized beds was carried out by Harris et al. (2002) and Yoshida et al. (1969), and involved taking mass balances. However, modeling the dynamics of fluidized beds poses a challenge due to particle size distribution (PSD), influence of geometric parameters, such as column diameter, height of bed to column diameter ratio, properties of gas and particles or mixture of particles, operating parameters such as superficial gas velocity on bed pressure drop, and minimum fluidization velocity. Due to the influence of the aforementioned parameters, the GS flow in a fluidized bed changes and encounters different regimes, namely homogeneous, bubbling, turbulent, fast fluidization, and pneumatic regimes (Bi and Grace, 1995;Lim et al., 1995;Bi, 2011). Geldart (1973) defined four particle size groups for bed expansion, known as Geldart-A, -B, -C, and -D. Based on these particle sizes, it can be determined whether the bed can be fluidized, how much the bed can be fluidized, and the type of fluidization. With the advent of the 21st century, researchers emphasized the need for studies on the hydrodynamics of Geldart-B-type particles and binary mixtures (Zhang et al., 2006). In fluidized beds with binary systems, the primary fraction that forms the top layer of the bed, or the one that floats, is called the flotsam and the one at the bottom layer, or the one that sinks, is called the jetsam. A major challenge in binary systems is that the beds reach equilibrium with either mixing or segregation of particles as two extremes. These dynamics are studied using advanced experimental techniques that include non-intrusive methods like tomography, radioactive particle tracking (RPT) (Roy et al., 2021), and intrusive techniques that include pressure probes for measuring pressure drops and quality of fluidization, and optical probes for measuring particle diameter, particle velocity, etc. Dynamics and transport phenomena in fluidized beds have also been studied both qualitatively and quantitatively using computational fluid dynamics (CFD) (mostly Eulerian-Eulerian approaches)/mathematical modeling (Cooper and Coronella, 2005;Du et al., 2006;Gao et al., 2009;Pei et al., 2010;Zaabout et al., 2010;Chang et al., 2012;Mostafazadeh et al., 2013;Benzarti et al., 2014;Sahoo and Sahoo, 2016;Bakshi et al., 2017;Agrawal et al., 2018;Chang et al., 2019;Daryus et al., 2019;Khezri et al., 2019;Shrestha et al., 2019;Kotoky et al., 2020).
In the current study, the experimental lab scale fluidized bed was operated in transition regime for a binary mixture. Hence, relevant literature works on bubbling, slugging, turbulent, and fast fluidizing regimes are summarized in Table 1. In this part of the work, numerical studies on qualitative analysis of the low patterns were elaborated, showing prominent experimental and numerical works highlighting the operating regimes, operating parameters, major findings, limitations/opportunities, geometrical details of equipment (diameter and height), and particle characteristics (size, shape, etc.).

Numerical and experimental studies on fluidized beds involving binary systems
The major challenges in CFD modeling of fluidized beds involve the modeling of solid-solid and solid-fluid interactions with the help of the kinetic theory of granular flow (KTGF) for Eulerian-Eulerian models. Furthermore, recent studies suggest that 3D models capture the dynamics better than 2D models. The numerical studies explained in this section focus on bubbling fluidized beds (BFBs) and turbulent fluidized beds (TFBs) operating with binary mixtures, along with a few studies on unary beds. The following is a discussion of numerical studies carried out over the last few decades, along with combined experimental and numerical studies. Huilin et al. (2003) performed simulations with binary mixtures using a KTGF model and the Euler-Euler approach (a multifluid model). The authors investigated the segregation of GS fluidized beds for binary mixtures for a bed height of 0.4 m and column diameter of 0.3 m. Finer particles tended to go up in the bed, while larger particles settled at the bottom, at lower superficial gas velocities (U gs = 1.6 m/s). The authors found that at 10 s, complete segregation occurred at a superficial gas velocity of U gs = 1.6 m/s. Furthermore, the authors observed that, with further increase in superficial gas velocity, solid volume fraction was more uniform in the axial direction. The authors concluded that the correct dynamics depended on distribution of particle size and energy dissipation due to solid-solid interactions.
Additionally, Philippsen et al. (2015) investigated the effect of various drag models to be used in fluidization and found that the Syamlal-O'Brien model (Syamlal and O'Brien, 1987) was the best drag model through which to assess the dynamics of fluidized beds. Daryus et al. (2019) compared two turbulence models, namely the standard k-ε and the RNG k-ε models to understand the effects of turbulence on CFD simulations of fluidized beds. The authors concluded that, while neither model could accurately predict the pressure drop for superficial gas velocities of less than minimum fluidization velocity (U mf ), pressure drops were predicted accurately when the superficial gas velocities were higher than U mf . Furthermore, the RNG k-ε model was found to predict the regimes and the static pressure distribution more accurately than the standard k-ε models. Cooper and Coronella, (2005) carried out numerical simulations for a bubbling fluidized bed reactor in the titanium refining industry with rutile (small size and high-density) and coke (large size and low density) particles. Outcomes signifying the importance of numerical simulations and bed characteristics included: 1. prediction of accurate dynamic similarity in flow patterns using mixing and segregation during scaleup; 2. prediction of bubble wake formations directly below the gas bubble and dynamics of the wake below the bubbles as depicted in their solid volume fraction contours; 3. eruption of the bubbles causing deposition of solids at the bed surface; and 4. downward movement of those bubbles that did not travel in the bubble wake; 5. flotsam and jetsam had similar but distinct velocity trajectories; 6. a minor difference in apparent slip velocity of bubbles and its influence on bed dynamics over     Objectives and findings: 1. Criterion for a mixed and segregated bed for binary systems was derived. 2. Comprehensive plots of pressure drop v/s superficial gas velocities for operating circulating fluidized beds were presented. 3. Analytical model to predict segregation in a FB, having binary mixture of different materials. 4. CFD simulations to study the effect of particle shape on bubble dynamics in bubbling FB. The bubble dynamics are significantly different for different shapes. 5. Correlation development for minimum fluidization velocity. Effect of unary and binary particle size distributions for different temperature and pressure. 6. CFD simulations for BFB of binary mixtures considering effect of particle size distributions and energy dissipation due to non-ideal particle-particle interactions. Importance of the KTGF model was highlighted. 7. Experimental study of binary mixtures of three different types of same sized and different density particles. Correlations were developed for minimum fluidization velocities. 8. Experimental investigations on segregation and mixing characteristics of BFBs containing Geldart-B particles. 9. CFD studies for mixing and segregation of binary mixtures in BFBs. The model developed was able to predict the characteristics for different operating conditions of binary mixtures. 10. Experimental investigation of segregation characteristics for binary mixtures in dual fluidized beds for change in operating parameters like fluidization velocity, particle mixture properties, and solid holdup. 11. Hydrodynamic characteristics of binary beds are significantly different than unary beds. 12. Comparative study of cluster formations, mass flux variation, and segregations in turbulent fluidization and fast fluidization regimes. 13. Cluster formation probability was higher in turbulent fluidized beds, while segregation extents were the same. 14. Mass fluxes were more dependent on particle properties in turbulent regime than fast fluidized regimes. 15. Hydrodynamics in a turbulent FB with binary mixture of polydisperse particles were studied using CFD with population balance. 16. The model gave important insights into the dynamics of particles with small and large differences in particle size distributions. 17. Study of dynamics of turbulent FB for different column diameters and development of a correlation for minimum fluidization velocity. 18. Experimental study of mixing and segregation of binary mixtures consisting of different percentage of biomass (5%-20%). 19. Amount of mixing increased with increased superficial gas velocity up to biomass concentration of 20%, after which it decreased. 20. Most of the materials used were non-spherical in shape and size and greater than 1 mm. 21. Two-dimensional CFD simulations using Geldart-B particles comparing standard k-ε and RNG k-ε models were carried out. The RNG k-ε model was found to be better than the standard k-ε model. 22. Experimental study of the effect of higher proportion of large particle sizes on bubble rise velocities. 23. An important finding was that increased higher fraction decreased bubble rise velocities, and hence particle velocities. 24. Data for velocity distribution were not available for binary and polydisperse beds. 25. Correlation development for minimum fluidization velocity and pressure drop using the artificial neural network. 26. Experimental investigations of pressure drop with the superficial gas velocity profile for binary mixture and compositions of particles at different places. 27. Experimental measurements and correlation development for minimum fluidization velocity for binary mixtures. 28. Predictions showed good match with experimental and published data. 29. Experimental investigations and correlation development for minimum fluidization velocity and pressure drop for binary mixtures. Predictions showed good match with experimental and published data. 30. Experimental investigations and correlation development for minimum fluidization velocity and pressure drop. Predictions showed good match with experimental data. 31. Experimental investigations and correlation development for minimum fluidization velocity and pressure drop. 32. CFD investigations were carried out for Geldart-A and -C type particles. 33. Sensitivity analyses of various geometric, operating, particle shape, size, and density were performed, and gas and solid fractions were analyzed. 34. The results will aid in design of fluidized bed reactors. 35. Experimental investigations were carried out in the dilute region of the riser in a CFB. 36. The shape of the axial and transverse profiles were dependent on the bed height and superficial gas velocities. 37. For lower bed heights, the movement is toward the center; with increased bed height, the movement is toward the wall. 38. Experimental investigation of gas mixing and axial dispersion in a bubbling fluidized bed using the RTD approach for linear low-density polyethylene was carried out, and a correlation for the dimensionless dispersion coefficient relating Re and aspect ratio was developed. 39. Experimental investigations to measure axial velocities of rising and falling particles were carried out for circulating fluidized bed riser. 40. The axial and transverse particle velocities were affected by superficial gas velocities and solids circulation rate. 41. Numerical investigations were carried out to understand the effect of particle diameter on bubbling GS fluidized beds. 42. Particle velocities decreased with increased particle diameter, which increased particle volume fractions. 43. The increase in particle velocity in the fluidization zone was higher for smaller particles and decreased with particle diameter. 44. CFD investigations of a turbulent fluidized bed with 2D and 3D simulations. Three-dimensional simulations were found to be more sensitive to specularity and restitution coefficients. 45. Two-dimensional simulations over-estimated particle volume fractions in the middle and top of the bed. 46. Twodimensional CFD simulations for FBs were carried out. 47. Increased bed height led to increased bed height but decreased average diameter of particles in the bed. 48. Significance of restitution coefficient in understanding the false segregation in beds in numerical simulations. 49. CFD simulations to study the effect of mesh size on transition from homogeneous to bubbling regime using Eulerian-Eulerian models. 50. The presence of a dilute region was dependent on selection of drag law, with Gidaspow and Syamlal-O'Brien models showing good predictions that omitted frictional stress and improper wall boundary conditions and showed appropriate minimum bubbling velocities. 51. Investigation of mixing dynamics and their dependence on operating conditions using CFD simulations for fluidized bed biomass gasification. 52. Bubble-induced solid micro-mixing induced solids up flow in nose and wake regions, and down flow along bubble walls. 53. Development of an analytical model for the fluidized bed. 54. Solid mixing was adversely affected in the presence of gas bypass, particularly in cases of heavier particles. 55. Threedimensional CFD investigations to analyze the capabilities of different drag models to predict the dynamics of turbulent fluidized beds filled with Geldart-B particles. 56. The Gidaspow model was found to be the best to predict drag coefficients per this investigation. 57. Development of a new experimental method to test the reactivity of particles in a GS fluidized bed. 58. CFD modeling to study heat transfer between particles in a fluidized bed. 59. Heat transfer coefficient increased with large particle size and superficial gas velocity. Limitations: 1. Investigations of flow patterns depicting volume fractions of solid particles, bubble dynamics, and mixing were not performed. 2. Axial velocity profiles for gas and particles across radial distance for different axial positions were not performed. 3. Particle volume fraction profiles across vertical centerline were not shown. 4. Correlations for minimum fluidization velocity and pressure drop were not developed. 5. The work was limited to binary mixture of particles of same size. 6. Densities of the particles were the same. 7. The particle sizes were larger than the Geldart classification of sizes. 8. Experimental investigations/CFD simulations were carried out at the same bed height. 9. Sensitivity analyses in terms of superficial gas velocities and different combinations of particle diameters were not performed. 10. Two-dimensional simulations were carried out, which did not show good predictions in the middle and top parts of the beds. 11. The work was limited to single particles of different sizes.

Frontiers in Energy Research
frontiersin.org passage of both time and additional bubbles. The authors have substantiated the aforementioned outcomes for a wide range of particle sizes and superficial gas velocities. Mazzei et al. (2010) carried out numerical simulations for a binary mixture of particles to understand two cases. Case 1: Investigation of minimum fluidization velocities at which the mixture no longer remains fixed, but starts segregating, and transient fluidization takes place; and Case 2: The mixture becomes steadily fluidized and fully mixed. The authors assessed the following: 1. numerical stability of models in fast segregating beds; 2. mutual effects of plastic viscosity and granular temperature; 3. the role played by plastic solid stress; and 4. selection of an appropriate time-step to ensure invariance of numerical results. The authors emphasized the role of plastic stress in the modeling of collapsing monodisperse fluidized beds. The authors found that, in the case of collapsing monodisperse fluidized beds, plastic stress and plastic solid viscosity are important, whereas plastic solid pressure can be ignored. The authors further used the standardized parameter to find the bed characteristics; for instance, the minimum fluidization velocity (U mf ), superficial gas velocity (U gs ), necessary for complete mixing, and velocity for oscillating pressure drop. The authors used the multifluid model or KTGF model as specified via Ansys Fluent software.
Jayarathna and Halvorsen, 2011 carried out both experimental and numerical investigations with different binary mixtures of glass particles and studied the pressure drop and volume fraction changes for lab scale fluidized beds. The authors conducted experiments for two different bed heights, each for a range of superficial gas velocities (U gs = 0.3-1 m/s). The authors validated their numerical model through experimental measurements of pressure drop with CFD predictions and found moderate agreement due to lack of distributor availability. Furthermore, the authors observed that, at slugging conditions, bubbles were moving upward using a zigzag path. Mostafazadeh et al. (2013), with the help of their in-house code, carried out numerical investigations for mixtures of 1 mm and 2 mm particles with densities of 2,400 kgm -3 and 2,500 kgm -3 , respectively, for a superficial gas velocity range of U gs = 0.5-2.5 m/s. The authors observed that increased mass fraction of small particles from 49% to 59% led to increased bed height and decreased average diameter of particles in the bed. The authors also observed that differences in restitution coefficient can cause segregation, even among particles of the same size and density. Hence, an appropriate restitution coefficient value is needed for suitable bed characteristics. Benzarti et al. (2014) examined the ability of the mathematical/CFD models to predict dynamics of TFBs filled with Geldart-B particles. The authors investigated the significance of drag coefficient models and restitution coefficient values on the prediction of dynamics of fluidized beds in CFD. The authors concluded that the restitution coefficient, which accounts for the inelasticity of the particle-particle collisions, needs to be considered, especially when the superficial gas velocity is less than the minimum fluidization velocity. The authors concluded that, for Geldart-B particles, the Gidaspow model (Huilin et al., 2003) gave the most reasonable results, both in terms of qualitative and quantitative predictions. Furthermore, with a specularity coefficient value of 1 and a restitution coefficient of 0.9, the model gave near accurate predictions. While analyzing the effect of superficial gas velocity (U gs ), the authors also found that increased superficial gas velocity caused particles to be entrained into the dilute region of a turbulent fluidized bed. Sande and Ray (2014) carried out numerical studies of transition from a homogeneous to bubbling regime for Geldart-A particles and concluded that the drag laws played an important role in the identification of the dilute region of fluidization. The authors, in their qualitative analysis using CFD, also found that inappropriate selection of wall boundary conditions and inclusion of frictional stress led to inappropriate predictions of minimum fluidization velocity. Such studies have not been carried out for Geldart-A, -B, -C, and -D particles. Both the Gidaspow (Huilin et al., 2003) and Syamlal-O'Brien models (Syamlal and O'Brien, 1987) gave good results for moderate superficial gas velocities (of approximately U gs = 0.008 m/s), whereas for other velocities (of approximately U gs = 0.01 m/s), the Wen Yu drag law model gave good results. Sahoo and Sahoo (2016) carried out CFD simulations for Geldart-C and -A fine particles (monodisperse particles) in a cylindrical fluidized column. The effect of parameters, such as static bed height, particle density, size of particle, and superficial velocity of fluidizing medium were studied and compared with experimental results. The bed expansion and pressure drop variation with increased superficial velocity of the gas was found to be similar to that of conventional fluidized beds. The authors were able to simulate and confirm that fluidization under normal conditions is a challenge for Geldart particles due to action of strong cohesive forces. Bakshi et al. (2017) carried out CFD simulations to study the effects of solids mixing on thermal and concentration gradients, and on the performance of fluidized bed reactors. The authors found that the bubble-induced solids were responsible for the micro-mixing during the up flow of the solids. This included the wake region during the up flow of solids. Furthermore, the mixing of solids was affected by gas bypass or through flow, particularly during fluidization of heavier particles. The authors also investigated the dynamics of the motion of gas and solids, and their interaction, under specific operating conditions. Chang et al. (2019) studied dynamics in fluidized beds with Geldart-B particles (binary systems) with 2D and 3D simulation approaches. An important aspect in Eulerian-Eulerian modeling is the restitution coefficient. Hence, analysis of sensitivity of the restitution coefficient was carried out and it was determined that a value of 0.9-1 for the restitution coefficient predicted realistic results for Geldart-B particles. Furthermore, the authors found that 2D simulations predicted the dynamics of the dense phase (bottom layer) well, whereas they over-estimated aspects of the dynamics of the middle and upper regions. The 2D simulations also over-estimated the bubble sizes and bed expansion, solid concentration, and solid velocities compared to experimental results. Hence, the authors suggested that 3D simulations should be carried out to obtain realistic results in studies of the dynamics of fluidized beds with Geldart-B particles. Kotoky et al. (2020) carried out CFD simulations using an inhouse code for Geldart-B particles to analyze the bed dynamics of unary fluidized bed reactors. The authors concluded that, with increased particle diameter, the particle velocity at any section in a fluidized zone decreased while the particle volume fraction increased, i.e., particle velocities were higher for smaller sized particles, especially in the dilute region of the bed, whereas velocities were lower for larger particle sizes. Hence, maximum value of time-averaged volume fractions was found for larger particles at the bottom of the reactor.

Frontiers in Energy Research
frontiersin.org

Discussion
From Section 1.1: 1. Turbulence models like the RNG k-ε model are better than the standard k-ε model for both unary and binary mixtures. 2. KTGF can capture bubble dynamics for bubbling beds, including bubble movement in the bed, bubble wake, and bubble eruption. 3. Roles of plastic solid stress, plastic viscosity, granular temperature, plastic solid pressure in unary/monodispersed beds using commercial software Ansys Fluent have been standardized and found to predict bed characteristics well via comparison with experimental measurements. 4. The standardized values of the restitution coefficient and specularity coefficient should be used in predicting correct bed characteristics using CFD models. The values reported in the literature are in the range of 0.9-1. 5. For predicting the suitable drag coefficient, the Gidaspow (Huilin et al., 2003) and Syamlal-O'Brien (Syamlal and O'Brien, 1987) models were found to be most appropriate when focusing on bed dynamics. However, if thermal and concentration gradients are coupled with bed dynamics, the Wen and Yu drag law provides better results.

Objective of the present work
In the present work, the flow patterns (both steady and transient) of binary mixture particles with the same densities were investigated using CFD simulations. Geldart-B particles were used at different operating conditions. For this purpose, geometry available in the current literature (Jayarathna and Halvorsen, 2011) was considered. The CFD model considers the standard values for different parameters, such as friction pressure, plastic viscosity, plastic pressure, specularity coefficient, and restitution coefficient, as reported in the literature, and the drag and other laws used in KTGF modeling. In future work, a sensitivity analysis will be carried out for different combinations wherever suitable. The model will then be validated with experimental data from the literature. In the case of good agreement, seven different binary mixtures will be taken, and simulations for three different superficial velocities, each for two different bed heights, will be carried out. The quantitative analysis will be carried out in Part II of the study and reported in a subsequent article.
The originality of this manuscript lies in: 1. the comprehensive combination of the particle size of binary mixtures and the operating parameters considered; 2. the CFD model that considered all the current best practices; and 3. investigation of whether unusual bed characteristics were present in any of the cases considered.
2 Mathematical modeling 2.1 Assumptions 1. No mass transfer between the phases is taking place in the system. 2. Two different solid phases of the same density, but containing particles of different sizes, are simultaneously interacting with each other and with the gas phase. 3. All the solid particles are spherical. 4. The gas fluid phase is a Newtonian fluid. 5. No other force or energy, other than gravity, is affecting the fluidized bed system in any manner.

Models
Different models were used for modeling the interaction between the solid phase and the gas phase. Table 2 shows the models used for quantities.

Mathematical modeling with equations
The Eulerian model, or two-fluid model, considers each phase as a continuum, where the phases are interacting and interpenetrating in nature. The solid phase may be assumed to be a pseudo-fluid. For the given study, the Eulerian model is used for the modeling of the fluidized bed system.

Continuity equations
The continuity equation of a phase i is given by: where ε represents the volume fraction of the phase i, ρ represents the density of the phase i, and U represents the velocity of the phase i.
In the system, there are three phases interacting with each other. They are given by: 1. g representing the fluid gas phase; 2. s 1 representing solid phase with smaller particle size; and 3. s 2 representing solid phase with larger particle size.

Momentum equations
The momentum equation for the gas fluid phase is given as: z zt ε g ρ g U g + ∇ · ε g ρ g U g U g −∇ · P g + ∇ · Π g + ε g ρ g g where P g represents the total fluid pressure, Π g represents the stress tensor of the gas phase, g represents gravitational acceleration, and ϕ g,si represents the drag interaction coefficient between the gas phase and the solid phase s i . The stress tensor of the gas phase is given by: where μ g represents the dynamic viscosity of the gas phase, ∇ · U g represents the divergence of U g , and ∇ · U T g represents the divergence of the transpose of U g .
The solid phase momentum equation is given by: where Π si represents the total stress tensor for the phase s i and ϕ si,sj represents the drag interaction coefficient interacting between the solid phase s i and the solid phase s j . The total solid phase tensor is given by: where P si represents the total solid phase pressure of the phase s i , I si represents the moment of inertia of particles of the phase s i , ξ si represents the granular bulk phase viscosity of the phase s i , and μ si represents the solid phase granular viscosity for the phase s i . The solid phase granular viscosity given by Syamlal and O'Brien (1987) is: where μ si(col) represents the collisional viscosity, μ si(kin) represents the kinetic viscosity, and μ si(fr) represents the frictional viscosity. The collisional viscosity is given as: where d p,si represents the particle size diameter of the phase s i , which is the same for all, G si,sj represents the radial distribution of the solid-solid particle interaction between the solid phases s i and s j , Θ si represents the granular temperature of the phase s i , and e si represents the total coefficient of restitution for the phase s i . The kinetic viscosity is given by: The total solid phase pressure is given by the Ma-Ahmadi model (Ahmadi and Ma, 1990) as: The granular bulk phase viscosity is given by the mathematical model of Lun et al. (1984), as shown: The total coefficient of restitution is given as: e si e si,si + e si,sj 2 where e si,si represents the coefficient of restitution between the similar particles of the phase s i and e si,sj represents the coefficient of restitution between the dissimilar particles of the phases s i and s j . The radial distribution function is given by the Ma-Ahmadi model (Ahmadi and Ma, 1990): where ε si(max ) represents the maximum possible volume fraction for the solid phase s i . The granular temperature is calculated using the algebraic model: where k 1 , k 2 , k 3 , and k 4 are equation constants given by:

Turbulence governing equations
The turbulence-based modeling of the system was carried out using the Renormalization Group RNG k-ϵ model for turbulent viscosity, since the previous sensitivity analysis provided good results. The model equations are similar to the standard k-ϵ model, with the constant C μ in turbulent viscosity modeled by a differential equation. A constant value of 0.0845 can also be derived from the differential equations. In the current study, the constant value was provided. The model uses the following equations: where r and z represent the directions, μ tur represents the turbulent viscosity, k represents the turbulence kinetic energy, ϵ represents the dissipation rate of turbulence kinetic energy, U r represents the component of U g in the direction of r, Y k represents the turbulence kinetic generation due to mean velocity gradients, Y b represents the buoyancy turbulence kinetic energy generation, Y c
In calculation of the aforementioned quantities, S represents the modulus of mean rate of strain tensor, β represents the coefficient of thermal expansion, g r represents the component of gravity in the direction of r, Pr tur represents the turbulent Prandtl number given as 0.85, M tur represents the turbulent Mach number, and Z represents the compressibility of the fluid gas.

Kinetic energy equations
KTGF is used for kinetic based modeling of the fluidized bed system and is the extended version of the kinetic theory of gases. The model assumes unequal granular temperature for different phases and uses collisions as a potential source of energy transfer and a variable affecting the granular temperature. The model equation for granular temperature of a solid phase is as follows: where q si represents the collisional heat flux for solid phase s i and γ si represents the dissipation of the turbulent kinetic energy due to particle collisions. The collisional heat flux is given by: where P si(col) represents the collisional pressure generated by particle collisions, m s represents the combined mass of the solid phases s i and s j , and d si,sj represents the average particle size of the solid phases s i and s j . The collisional pressure given by Gidaspow and Huilin (1996) is: where Δ is an equation constant, n si and n sj represent the total number of particles of the solid phases s i and s j , respectively, and m si and m sj represent the single particle masses of the solid phases s i and s j , respectively. The equation constant Δ is given as: The average particle size of two solid phases is given by: where d si represents the particle size of the solid phase s i and d sj represents the particle size of the solid phase s j . The combined mass for two solid phases is defined as: The single particle mass is calculated as: The total number of particles is defined as: The turbulent kinetic energy dissipation by particle collisions is given by the Gidaspow and Huilin model (Gidaspow and Huilin, 1996) as:

Drag equations
The Syamlal-O'Brien model (Syamlal and O'Brien, 1987) was used for the drag modeling of the fluidized bed system. The model equation for the gas-solid particle drag interaction is as follows: where C dr represents the drag coefficient of the gas-solid system represented by Eq. 38, ϑ r is the equation constant, and where Re si is the Reynolds number of the solid phase s i . The Reynolds number for the solid phase is given as: The equation constant ϑ r is represented by the following equation: where A and B are given as: The particle-particle drag interaction coefficient is governed by the Syamlal and O'Brien model (Syamlal and O'Brien, 1987) as: ϕ si,sj 3 1 + e si,sj π 2 + π 2 f 8 ε si ε sj ρ si ρ sj d si + d sj where f represents the coefficient of friction from interaction between the solid phases s i and s j .

Frictional equations
The frictional pressure is derived from KTGF and is given by Syamlal et al. (1993) as: The frictional viscosity given by Schaffer is as follows: where P fr represents the frictional pressure, α is the angle of internal friction taken as 30°, and Π si(2D) represents the second invariant of the deviatoric stress tensor.

Geometry and mesh details
A 3D cylindrical geometry of radius 0.036 m with a height of 1.4 m was created in Ansys Workbench 18.1. The bottom was designated as the air inlet, while the upper circular geometry was used as the outlet. Figure 1A shows the 2D plane from 3D geometry with bed height and axial positions at t = 0 s, and Figure 1B shows the 3D geometry representation at t = 0s. The diameter of the cylinder is given by D 1 , and H represents the total height of the cylinder. At t = 0 s, h represents the initial bed height of the glass particles. Figure 2 shows different views of the mesh used for simulations, including the axial and radial zoomed views of the 3D geometry.

Frontiers in Energy Research
frontiersin.org

Boundary conditions
"Velocity-inlet" is used as a boundary for the gas inlet condition. "Pressure-outlet" is used as a boundary condition for outlet. The packing limit for glass particles was taken as 0.63. No-slip condition was applied on the walls. The time-step used in simulation was 0.001 s. The total simulation was run for 7 s.

Material properties
Glass particles were used as the solid phase and air was used as the fluid phase. Glass particles of two different sizes, 154 μm and 488 μm, were used, and average particle sizes varied between groups. The densities that were used are as follows: 2485 kgm -3 for glass and 1.22 kgm -3 for air. The viscosities used were 0.00082 kgm -1 s -1 for glass and 0.000017 kgm -1 s -1 for air. The material properties used for simulations were as indicated by Jayarathna and Halvorsen (2011).

Grid sensitivity
Three different meshes were used for simulations: Mesh 1 with 173,040 elements, Mesh 2 with 267,786 elements, and Mesh 3 with 497,568 elements. The mesh elements are hexahedral and more refined near the wall, with a near-wall yplus of around 30. The initial volume fractions in the 2D plane are shown for both bed heights (h s1 = 0.335 m and h s1 = 0.635 m) in Figure 3. The axial velocity magnitudes in the radial direction for a bed height of 0.335 m were plotted. Figure 4 shows the radial velocity profile at Position 2 of Figure 1 for different meshes. The deviation between Mesh 2 and Mesh 3 was 2%, while maximum deviation for Mesh 1 and Mesh 2 was 10%. Hence, Mesh 2 was used. Figure 4 shows axial and radial views of Mesh 2.

Method of solution
The simulations were carried out using commercial fluid software Ansys Fluent 18.1. A first order upwind scheme was

Frontiers in Energy Research
frontiersin.org used to solve momentum, volume fraction, turbulent kinetic energy, and turbulent dissipation rate equations. A phasecoupled SIMPLE scheme was used to solve pressure-velocity coupling. For transient formulation, a first order implicit scheme was used. Convergence criterion for continuity was 0.001; it was 10 −4 for other equations. The parametric data (for initial bed heights, as shown in Figure 3, including superficial velocity, binary mixtures, and both individual and average particle sizes) used for the simulations are shown in Table 3.

Results and discussion
First, the standard values and models for the parameters, as discussed in Section 1.2, were chosen from those available in the literature. Similarly, the drag and turbulence model was chosen per the literature. The model was validated with experimental data available from published studies that used these standard settings. The transient solid particle dynamics in the bed was then presented in the form of qualitative solid volume fraction contours to understand the segregation and mixing characteristics for different particle size mixtures considered in the study. In this section, 100 × x% mixture represents the percentage of large particles and 100 × (1 − x)% represents the percentage of the small particles. All the contours presented in the figures are steady-state time-averaged volume fraction contours of particles. Herein, x represents the weight fraction of large particles. Lim et al. (1995) have emphasized the importance of particle size, particle composition, and baffles, which lead to transition from bubbling or slugging regime to turbulent regime. The criteria for the dimensionless velocity that characterizes the regimes are given by Eq. 3. An effort has been made to identify the regime in which the present work was carried out, per analysis demonstrated by Lim et al. (1995). Eqs 46-49 represent the dimensionless numbers and velocities, as well as the average particle diameter for a binary mixture. Figure 5 shows the plot of the regime analysis  for the superficial velocities considered in the present work. The analysis shows that the entire zone is in bubbling regime. However, it must be noted that the analysis derived by Lim et al. (1995) was based on experimental data from unary beds.

Regime analysis
where x si is the initial weight fraction of the solid phase s i .

Model validation
For model validation, two mixtures of 0% and 40% were simulated at 0.235 m of initial bed height and superficial gas velocities varying from U gs = 0.184 m/s to 0.225 m/s each. Figure 6 shows a deviation of around 5%-7% between experimental data and numerical predictions for a binary mixture with 0% large sized particles (or 100% small particles), whereas there was less than 3% deviation for a binary mixture with 40% large particles. The deviation is attributed to the absence of distributor details from the published literature. Figure 7i shows the steady-state time-averaged solid volume fraction contours for different superficial gas velocities and 0% mixture (100% fine particles) for a bed height of 0.635 m. For superficial gas velocity of U gs = 0.3 m/s, a well-mixed pattern can be observed. However, slugs of particles seemed to deposit at different axial locations at the walls. A large bubble with solid particles was observed at an axial location of z/H = 0.5 when the maximum height of fluidization was 1.25 m at steady state. A similar flow pattern was observed for a higher superficial gas velocity of U gs = 0.45 m/s. However, the solids were deposited on the right-hand wall of the bed between the dimensionless heights of z/H = 0.35 and z/H = 0.8. A bubble formed on the left-hand wall with particles moving in the space between the bubble and wall. For a superficial gas velocity of U gs = 0.6 m/s, a large bubble was seen at the outlet with a wake below and followed by another bubble. A prominent zigzag pattern was observed from bottom to top, with slugs of solid particles alternating on the right and left wall. These results confirm the bubble wake and bubble formation, as has been reported in published literature (Cooper and Coronella, 2005).

FIGURE 6
Variation of pressure drop as function of gas superficial velocity (U gs ) for bed height = 0.235 m. Numbers 1 and 2 denote CFD simulations for increasing order of particle sizes, respectively, while symbols denote experimental measurements. ■-0% large particles; ▲-40% large particles.

Frontiers in Energy Research frontiersin.org
Furthermore, bubble formation and dynamics were also captured by the model. Figure 7ii shows steady-state, time-averaged volume fraction contours for two different bed heights of 0.335 m and 0.635 m for the three superficial gas velocities. The binary mixture contained 20% large particles and 80% small particles. Figure 7iiA shows that, for initial bed height of 0.335 m, the fluidized bed steady-state heights were 0.45 m, 0.5 m, and 0.65 m for superficial gas velocities of U gs = 0.3 m/s, U gs = 0.45 m/s, and U gs = 0.6 m/s, respectively. A well-mixed pattern was observed for the three superficial gas velocities considered at this initial bed height. With a superficial gas velocity of U gs = 0.45 m/s, a small layer of dense solid particles accumulated at the top, indicating that most of the finer particles go to the top, resulting in segregation. On the other hand, at a superficial gas velocity of U gs = 0.6 m/s, slugs of particles were formed, as in the previous case shown in Figure 7i. Furthermore, a large bubble was formed at the top of the fluidized bed at this superficial gas velocity. These patterns are similar to those observed by Lan et al. (2014), where partial segregation was predicted. Figure 7iiB shows the solid phase volume fractions for a higher initial bed height (z = 0.635 m) for the same set of conditions as in Figure 7iiA. When the superficial gas velocity was lower (U gs = 0.3 m/s), the finer particles formed larger slugs in the upper half of the bed, while the lower half had mostly coarser particles. The fluidized bed height was approximately 0.85 m. As the superficial

Frontiers in Energy Research
frontiersin.org gas velocity increased (U gs = 0.45 m/s), bubbles formed in the bed, while the slug sizes of finer particles decreased and became thinner and covered greater length at the top of the bed. A bubble formed at the midpoint of the bed (z = 0.55 m). A further increase in superficial gas velocity caused the top portion to be occupied by finer particles and the bottom portion to consist of coarser particles, with some area in the middle covered by a large bubble. Thus, complete segregation was observed at the highest velocity. This also corresponded with results reported by Lan et al. (2014), who found similar patterns where the top bed was well-mixed, while the bottom part was stagnant with coarser particles. Figure 8 shows the time-averaged steady state solid phase volume fraction contours for a binary mixture of 40% large particles and 60% small particles. Figure 8A shows that, for a lower bed height, there was complete segregation of flotsam and jetsam. However, mixing was observed when the superficial gas velocity was increased (U gs = 0.45 m/s), and complete segregation did not take place. With further increase in superficial gas velocity (U gs = 0.6 m/s), the following characteristics were observed: bubbles occupied the top area, while finer particles were restrained to the middle of the bed, and the top of the bed consisted of mixed particle sizes. The bottom of the bed consisted mostly of jetsam, which denotes intermediate  mixing. Figure 8B shows that, for a bed with an initial height of 0.635 m, volume fractions in jetsam were higher for lower superficial gas velocities (U gs = 0.3 m/s), and the fluidized bed height was 0.75 m, indicating that, due to the presence of large particles, there was less mixing and greater segregation. The scenario changed with an increase in superficial gas velocity. For U gs = 0.45 m/s, the bed was still segregated, but some mixing occurred. For 0.6 m/s, the bed was well mixed. Figure 9 shows the time-averaged steady-state solid-phase volume fraction for a binary mixture of 60% large particles and 40% small particles. An interesting observation can be made from Figure 9A for a superficial gas velocity of U gs = 0.3 m/s, where a jet of fluid rose and

Frontiers in Energy Research
frontiersin.org caused smaller bubbles to rise at the bed surface. With increased superficial gas velocity (U gs = 0.45 m/s), however, the bed tended to be segregated, while with further increased superficial gas velocity (U gs = 0.6 m/s), bubbles and slugs of solids started forming, indicating transition to a turbulent regime. Figure 9B shows a similar analysis for a higher bed height. Here, for a lower superficial gas velocity, the bed remained stagnant, while with increased superficial gas velocity there was transition from bubbling (U gs = 0.45 m/s) to a turbulent regime (U gs = 0.6 m/s). Figure 10 shows interesting results for a binary mixture with 80% large solids and 20% fine solids. Figure 10A shows that, for U gs = 0.3 m/s, the bed reached minimum fluidization, while for U gs = 0.45 m/s it had a fluid jet that entered the bed and a bubble that adhered to the wall. An extremely interesting flow pattern was observed for a velocity of a superficial gas velocity of U gs = 0.6 m/s. Alternate slugs of fine and dense mixtures were observed rising up the bed. Figure 10B shows similar patterns for a case of higher initial bed height. An interesting pattern was observed at a superficial gas velocity of U gs = 0.6 m/s, which showed bubble formation and its rise at the bottom of the bed similar to the one observed, both experimentally and numerically, by Cooper and Coronella (2005). Figure 11 shows time-averaged steady-state flow patterns for a binary mixture with 100% large solids and 0% fine solids. It was observed that for both bed heights of 0.335 m and 0.635 m and a low superficial gas velocity of U gs = 0.3 m/s, the bed remained as a fixed bed and no fluidization was possible, as shown in Figure 11A. For higher superficial gas velocities of U gs = 0.45 m/s and U gs = 0.6 m/s, a bubbling fluidized regime was observed. For a higher initial bed height of 0.635 m ( Figure 11B), bubbles formed at the bottom and adhered near the wall. No mixing was observed for the superficial gas velocities, but slug formation was observed for U gs = 0.6 m/s. Interestingly, until 40% large particle diameter, for most of the superficial gas velocities and bed heights, we observed small structures of bubbles and solid slugs which were representative of the slugging/ turbulent regime. However, the flow patterns of a binary mixture with 20%-80% large particles showed more mixing than the cases with 40% and 60% large particles.
Although both cases seemed to represent a turbulent fluidization regime, two distinct questions were posed: 1) Would mixtures of between 20% and 40% large particles be in turbulent range or transition range? and 2) what are the transient dynamics of this process?
Since all the cases from Figures 7-11 were steady-state timeaveraged, it was worth observing the transient flow patterns for an intermediate mixture composition and high superficial gas velocity for a different bed height (in between the bed heights already considered). A binary mixture of 25% small particles and 75% large particles was considered for analysis with an initial bed height of 0.48 m. Figure 12 shows the volume fraction contours for the same. Transient volume fraction contours show that, at the end of one second, a large bubble is formed at the top. After each subsequent second, the intermixing throughout the column was evident as turbulent fluidization (although the criterion for fluidization requires confirmation). A similar exercise with increasing superficial gas velocity showed that height of the bed increased as the superficial gas velocity increased, with bubbles forming at U gs = 0.3 m/s and U gs = 0.75 m/s. A detailed study on this can be carried out to determine whether the bed undergoes turbulent fluidization that indicates a turbulent regime, in contradiction to the bubbling regime predicted in Figure 5. This is, however, outside the scope of the present work.

Conclusion
The conclusions drawn from this study are: 1. Qualitative flow patterns and quantitative gas and particle velocity profiles indicate transition from bubbling and slugging regime to turbulent regime for some of the binary mixtures considered. These observations at dimensionless velocity, per the criteria of Lim et al. (1995), are different and may be attributed to the presence of particle size distribution, resulting in breakage of bubbles/slugs during bed expansion. Frontiers in Energy Research frontiersin.org 2. Low volume fractions of 0%-20% of large particles and low superficial gas velocity of U gs = 0.3 m/s with no large particles result in homogeneous regimes for both bed heights, while higher superficial gas velocities of U gs = 0.45 m/s and U gs = 0.6 m/s show intermixing at higher axial locations of the bed, and completely mixed steady-state profiles are observed. For mixtures with 20%-40% volume fraction range of large particles for both bed heights, gas bubbles were seen on the near-wall zone, and about 30% of the bed remained segregated at superficial velocities of 0.3 m/s and 0.45 m/s, with mixing restricted to the top part of the bed; the bed was well mixed under the 0.6 m/s condition. 3. For a 40% volume fraction of large particles, the bed remained 70% segregated at both bed heights and superficial velocity of 0.3 m/s. For a higher superficial velocity, the bed was well mixed. 4. For a 60% volume fraction of large particles, the bed was largely segregated for both bed heights and at lower superficial velocity of 0.3 m/s, while for a higher superficial velocity (0.45 m/s), the bed at lower height was well mixed. At the higher bed height, 80% of the bed was well mixed. Furthermore, for a higher superficial velocity of 0.6 m/s, the bed was well mixed at both bed heights. 5. For mixtures with 80% and 100% large particles at lower superficial velocity of 0.3 m/s, the bed did not fluidize, while fluidization of approximately 30%-35% was observed for a superficial velocity of 0.45 m/s at a bed height of 0.635 m. For a lower bed height, the amount of fluidization was around 65%. For the highest velocity considered, the bed was well mixed for lower bed height (0.335 m) and 70% mixed for a higher bed height (0.635 m).
5 Future work 1. In Part II of this two-part series, the present conclusions will be substantiated with comprehensive study of the gas and particle velocity profiles, as well as particle volume fraction profiles, for all the particle mixtures, bed heights, and superficial velocities considered in the present study. 2. Furthermore, the generated data will be used to develop correlations for minimum fluidization velocity and pressure drop for binary mixtures in Part II of the study. 3. Similar comprehensive studies will be taken up, focusing on simulations using the discrete element method (DEM) for labscale fluidized beds of Geldart-B and other Geldart group particles, which have practical applications depending on experimental data.
4. Simulation studies focusing on particles larger than 1 mm should be conducted, as has been reported in the experimental literature.

Data availability statement
The simulation data supporting the conclusion of this article will be made available by the authors, without undue reservation.