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 g s = 1.6 m/s). The authors found that at 10 s, complete segregation occurred at a superficial gas velocity of U g s = 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 m f ), pressure drops were predicted accurately when the superficial gas velocities were higher than U m f . 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 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 m f ), superficial gas velocity ( U g s ), 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 g s = 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 g s = 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 g s ), 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 g s = 0.008 m/s), whereas for other velocities (of approximately U g s = 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 in-house 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.

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.

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.

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.

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

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.

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₁, 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.

“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.

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^-1s^-1 for glass and 0.000017 kgm^-1s^-1 for air. The material properties used for simulations were as indicated by Jayarathna and Halvorsen (2011).

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.

The simulations were carried out using commercial fluid software Ansys Fluent 18.1. A first order upwind scheme was used to solve momentum, volume fraction, turbulent kinetic energy, and turbulent dissipation rate equations. A phase-coupled 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⁻⁴ 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.

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 (Lim et al., 1995) 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. R e g = ρ g U g d p s i , s j μ g (46) A r g = ρ g ρ s i − ρ g d p s i , s j 3 μ g (47) U * = R e g A r g 1 3 (48) d p s i , s j = ∑ i = 1 2 x s i d p , s i (49) where x s i is the initial weight fraction of the solid phase s_i.

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 g s = 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 g s = 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 g s = 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 g s = 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). 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 g s = 0.3 m/s, U g s = 0.45 m/s, and U g s = 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 g s = 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 g s = 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 g s = 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 gas velocity increased ( U g s = 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 g s = 0.45 m/s), and complete segregation did not take place. With further increase in superficial gas velocity ( U g s = 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 g s = 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 g s = 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 g s = 0.3 m/s, where a jet of fluid rose and caused smaller bubbles to rise at the bed surface. With increased superficial gas velocity ( U g s = 0.45 m/s), however, the bed tended to be segregated, while with further increased superficial gas velocity ( U g s = 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 g s = 0.45 m/s) to a turbulent regime ( U g s = 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 g s = 0.3 m/s, the bed reached minimum fluidization, while for U g s = 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 g s = 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 g s = 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 g s = 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 g s = 0.45 m/s and U g s = 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 g s = 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 time-averaged, 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 g s = 0.3 m/s and U g s = 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.