Conservation of mass and momentum is governed by the continuity and momentum equations for both solid and gas phases. For the two-fluid model (TFM), which is also known as the Eulerian-Eulerian model, the averaging approach is used as a multiphase model. By definition, the sum of volume fractions of all the phases is equal to one: α g + α s = 1 , (1) where α g and α s are volume fraction of gas and solid phase in the mixture respectively.

The continuity equations for gas and solid phases are solved separately with time and space and are expressed as follows: ∂ ∂ t ( α g ρ g ) + ∇ . ( α g ρ g V g → ) = 0 (2) ∂ ∂ t ( α s ρ s ) + ∇ . ( α s ρ s V s → ) = 0 (3)

The following set of equations represents the conservation of momentum for the gas phase and solid phase: ∂ ∂ t ( α g ρ g V g → ) + ∇ . ( α g ρ g V g → V g → ) = − α g ∇ p g + ∇ . τ g ¯ ¯ + α g ρ g g → + β g s ( V g → − V s → ) (4) ∂ ∂ t ( α s ρ s V s → ) + ∇ . ( α s ρ s V s → V s → ) = − α s ∇ p g + ∇ . τ s ¯ ¯ + α s ρ s g → + β g s ( V s → − V g → ) , (5) where V g → and V s → are gas and solid-phase velocities, respectively. Moreover, τ g ¯ ¯ is the gas-phase strain tensor and is expressed as τ g ¯ ¯ = α g μ g ( ∇ V g → + ∇ V g τ → + α g ( λ g − 2 3 μ g ) ) ∇ . V g → I ¯ ¯ , (6) where I ¯ ¯ is the unity tensor.

Similarly, solid-phase strain tensor is expressed as τ s ¯ ¯ = α s μ s ( ∇ V s → + ∇ V s τ → ) + α s ( λ s − 2 3 μ s ) ∇ . V s → I ¯ ¯ . (7)

The above Eqs 4, 5 are the interphase momentum transfer between gas and solid phase, which is an important term in gas-solid fluidization modeling. Solid particle fluidization results from the pressurized air (gas) drag exerted on the particulate phase. For the gas-solid interphase exchange coefficient, several drag models have been used so far. Syamlal–O'Brien drag function is more promising and accurate when compared with the other models, so it is more appropriate for predicting the hydrodynamics of gas-solid flows (Kaneko et al., 1999; Taghipour et al., 2005). Terminal velocities of particles are essential parameters in fluidization phenomena and drag law is based on its measured values (Syamlal and O'Brien, 1989; Syamlal and O'Brien, 2003). The empirical correlation for the exchange coefficient in terms of the volume fraction and relative Reynolds number is expressed as follows: β g s = 3 4 α s α g ρ g V r , s 2 d s   C D ( Re s V r . s ) | V s → − V g → | . (8)

For Syamlal–O'Brien drag function, the fluid-solid drag function is described by Dalla Valle (1948) according to the following relation: C D = ( 0.63 V r . s + 4.8 V r , s Re s ) 2 , (9) where Re s is the Reynolds number of solid particles described by Richardson and Zaki (1954) according to the following relation: Re s = ρ g d s | V s − V g → → | μ g . (10)

Terminal velocity of the solid-phase particles is described by Garside and Al-Dibouni (1977) using the following relation: V r , s = 0.5 ( A − 0.06 Re s + ( 0.06 Re s ) 2 + 0.12 Re s ( 2 B − A ) + A 2 ) , (11) where A = α g 4.14 and B = 0.8 α g 1.28 for α g ≤ 0.85 and B = α g 2.65 for α g > 0.85

Many studies have been carried out on the kinetic theory of fluidized bed behavior’s granular flows (KTGF). The detailed derivation has been explained by Gidaspow (1994) and Peirano and Leckner (1998). Similar to thermodynamic temperature in gas, a granular temperature in the solid phase is also a measure of velocity fluctuations and is defined as follows: T = 1 3 U ′ s 2 . (12)

The energy transport equation, which incorporates the granular temperature, is needed to describe this phenomenon. The following equation has been proposed by Ding and Gidaspow (1990): 3 2 [ ∂ ( α s ρ s T ) ∂ T + ∇ . ( α s ρ s T ) V → s ] = ( − p s I ¯ ¯ + τ s ¯ ¯ ) : ∇ V → s − ∇ . ( k T ∇ T ) − γ T − J T , (13) where k T is the granular conductivity, γ T is the dissipation due to inelastic particle-particle collisions, and J T is the dissipation or generation of granular energy caused by the exchange of fluctuating energy between the two phases.

In a fluidized bed simulation, the bulk viscosity of the mixture should be accounted for by Newtonian fluids. The solid shear and normal stresses occur due to granular particles’ random motion and strongly depend on the velocity gradients. There are usually three solid bulk viscosity parts; one is collision viscosity and the others are kinetic viscosity and frictional stress part. The collision viscosity is described by the following relation (Gidaspow et al., 1991): μ s , c o l = 4 5 α s ρ s d s g 0 , s s ( 1 + e s s ) ( T π ) 1 2 α s , (14) where g_0,ss is the radial distribution function, which is interpreted as a measure of particle contact probability–particle contact and the coefficient of restitution. g 0 , s s = ( 1 − α s α s , max ) − 2.5 α s , max , where α s , max is the maximum packing limit.

The kinetic part of the bulk viscosity is described by the following relation (Syamlal et al., 1993): μ s , k i n = α s d s ρ s T π 6 ( 3 − e s s ) [ 1 + 2 5 ( 1 + e s s ) ( 3 e s s − 1 ) α s g 0 , s s ] . (15)

When the solid-phase volume fraction becomes close to the packing limit, the friction stress part of viscosity becomes more critical (Syamlal and O'Brien, 1989). μ s , f r = p s ⁡ sin ⁡ φ 2 I 2 D . (16)

The solid pressure, the angle of internal friction, and I 2 D are the second invariants of the deviatoric stress tensor. The solid pressure is composed of both kinetic and collision terms as described by the following: p s = α s ρ s T + 2 ρ s ( 1 + e s s ) α s 2 g 0 , s s T . (17)

Transient three-dimensional (3D) simulations of the fluidized bed were performed. The solver used for the simulation was a commercial software package, ANSYS FLUENT 16.0 (Fluent, 2012). For multiphase modeling, the Eulerian multiphase model was used for the analysis. The pressure-based transient solver was used for the analysis. Phase coupled SIMPLE algorithm was used for pressure–velocity coupling. Governing partial differential equations were discretized using the second-order upwind discretization scheme, whereas the high-resolution interface capturing (HRIC) was used for volume fractions of phases. The bounded second-order implicit transient formulation was used for more accurate solutions. Typically, a time step of 0.001°s with 20 inner iterations per time step was used. This number of inner iterations was found to be adequate to achieve convergence for most time steps. For all the equations, 10^–3 residual convergence criterion was set during all the simulations.

The geometry of the gasifier used for the simulation is adopted from the experimental work of Yudin et al. (Yudin et al., 2016), having a cylindrical column of height H = 260 mm and diameter D = 108 mm as mentioned in Table 1. The air distributor plates were of three different designs made up of aluminum with 115 mm diameter “D” and 8 mm thickness “t,” as illustrated in Table 2. The first among the three distributors was a perforated distributor plate (Figure 1A) consisting of N = 89 holes having 4 mm diameter “d₀” each and arranged in an evenly spaced triangular pitch pattern. Such air distributors are generally used in most fluidized bed reactors and are considered reference distributors for conventional fluidized beds. The second distributor used a circular edged distributor with N = 8 air slots with 6.5 mm in width “W” and 29 mm long “L.” The air angle of attack in this distributor is perpendicular to the fluidized bed material ( α = 90°), as shown in Figure 1B. The third distributor was the modified version of the 90° distributors, a swirling inclined distributor having similar specs as 90° distributors, but the air angle of attack is inclined to 45° to the bed material shown in Figure 1C. The area opening ratio “ γ ” for both 90° and 45° distributors was 13%, and for the perforated plate distributor was 11%.

Geometric details of all three types of distributors are explained in Table 3.

Geometry details of the gasifier along with distributors are shown in Figure 2.

Due to the symmetry of geometry for perforated and 90° distributor plates, the quarter region of the geometry was created compared to the 45° slotted plate configurations for simulation.

Geometries of all types of distributors were created using design modular and computational mesh was generated using a meshing tool of ANSYS workbench. Due to the symmetric geometry of both perforated and 90° plates, the mesh was generated on a quarter half of the domain using symmetry boundary conditions (Figure 3). However, a full 360° mesh was generated for 45° plates due to nonsymmetric geometry. Before selecting the final mesh for the analysis, a grid-independent study (explained in the next section) was carried out to get more accurate simulation results.

Different mesh sizes for porous plates were evaluated in the present study. The coarse, medium, and fine meshes were 0.11, 0.28, and 0.45 million, respectively, for the perforated distributor plate. The pressure drop across the bed was measured for each mesh at various superficial velocities. The pressure drop for coarse mesh gave crude results as compared to medium and fine meshes. The sequence of mesh independence strategy is shown in Figure 4.

The comparison of pressure drops for coarse medium and fine meshes for perforated plate distributor is tabulated in Table 3.

Medium and fine mesh results were in close agreement with the experimental results, as shown in Table 3. The medium and fine meshes showed inconsiderable pressure drop differences with each other, so medium mesh was chosen as an optimum mesh for the rest of the simulations.

The solution was initialized with all the specified field variables of the entire computational domain’s solid and gas phases, as listed in Table 1. In addition, simulations were initialized with only air in the domain with no solid particles and superficial air inlet velocity.

The bed material was patched into the lower part of the domain so that the initial patched region containing bed material had an aspect ratio of 0.2 and 0.4 subsequently. The initial solids volume fraction patched in was 0.55. The initial air volume fraction patched in the solid-phase zone is 0.45 and above that zone, such as free broad area, is patched as 1. At the bottom inlet boundary, the constant specified gas inlet velocity and volume fraction were specified for each case. Air inlet velocities range from 0 to 5 m/s with an interval of 0.5 m/s for each simulation run for every distributor plate. The pressure outlet boundary condition was assigned to the upper outer boundary of the computational domain. For the gas phase, no-slip boundary conditions were applied at the wall section of the domain. In contrast, the following boundary equations were applied for the particle’s tangential velocity on the wall and the wall's granular temperature (Lun and Savage, 1987; Patankar, 2018). V s , w → = − 6 μ s α s , max 3 π φ ρ s α s g 0 , s s T s ∂ V s , w → ∂ n (18) T s , w = − k s T s e s s , w ∂ T s , w ∂ n + 3 π ρ s α s V s → g 0 T s 3 / 2 6 α s , max e s s , max . (19)

For the validation of the proposed computational model, simulation results were compared with the experimental work of Yudin et al. (2016). Pressure drop predicted by CFD simulation was compared with experimental findings using different plate designs. The range of simulations was carried out at different superficial velocities, i.e., from initial low velocity up to minimum fluidization velocity and even higher velocities to predict pressure drop across the sand bed. The particle size for the simulation was taken as 0.75 mm for all the distributor runs. Figures 5–7 show the pressure drop and their comparison with the published data. Computational results showed that by increasing the superficial air velocity, the pressure drop also increased to minimum fluidization velocity. The pressure drop was not further increased by increasing the velocity of air through the distributors.

To observe the effect of distributor configuration on fluidized bed hydrodynamics, the comparison of pressure drops results was obtained using different distributor designs. For example, perforated plate (4 mm diameter, 89 holes), circular edged slotted type (angle of attack of air 90°), and the novel swirling type (air angle of attack 45°) were investigated, as seen in Figure 8, which is the comparison of pressure drop for three different types of distributor plates. Pressure drop was highest for perforated distributor plate followed by 45° plate and lowest pressure drop was observed for 90° distributor plate.

The pressure drop for the perforated plate distributor is on the higher side due to the flow resistance offered by the smaller circular opening holes located at a relatively smaller distance apart. So significant pressure loss occurred due to the higher kinetic head development while crossing the plate’s narrow holes, depicting the nozzle effect. The early fluidization occurrence of a 45° distributor plate (Figure 8) is due to the inlet airflow's axial and radial components. This axial and sideways flow allows the packed particles of sand to become loosely bound since chaotic motion starts right from the lower portion of the bed, permitting early fluidization. The pressure drop using a 45° plate is slightly higher because the flow has to travel a longer path, thus experiencing more frictional head loss than 90° slotted plate (Aworinde et al., 2015). In fluidized bed hydrodynamics literature, it is generally found that distributor pressure drop becomes significantly high (Sreenivasan and Raghavan, 2002). This effect is most likely found in fluidized beds by using perforated plate distributors due to the formation of smaller bubbles and clusters (Rahimpour et al., 2017; Yudin et al., 2020).

It can be observed qualitatively from Figure 9 that solid-phase distribution is different for different distributor plates.

In general, the central region of the channel, whether its pipe or duct, has a maximum velocity at its center. However, the gas-solid mixing pattern in a fluidized bed is significantly affected by the distributor plate’s orifice arrangements (triangular pitch, radial pitch, or square pitch). In Figure 9A-1, initially, at time 1 s, the bubbles start to grow near the central region of the plate and rise until they reach the surface of the bed, where they collapse, causing the central region to become vacant. This vacant space is filled by the solids from the annular region of the bed, causing the airflow to make channels through the annular region as time proceeds, i.e., time, t = 3 and 5 s. This flow behavior is due to the radial distribution of orifice, which enhances the coalescence of bubbles through the annular region of the fluidized bed. Afrooz et al. (2006) have examined the similar behavior of solid particle movement inside a fluidized bed using a perforated distributor plate having a radial distribution of orifice arrangements.

Moreover, it is also concluded that at low fluidization velocity (u_o ∼ 1.5*u_mf) and having a radial distribution of orifice arrangement, more channeling of air bubbles is usually formed through the annular region of the bed with a broader solid central region. At higher superficial air velocities, such as u_o ≥ 2.5*u_mf, these air bubbles start rising through the orifices by making an angular path to form coalescence with the neighboring bubbles. It results in a slug flow that reaches the bed’s surface, where the larger air bubbles collapse and the process continues as time proceeds. Such transition phenomenon from lower to higher fluidization velocities results in the shift of the annular region of bubbles toward the central core of the fluidized bed. Due to the larger coalescence of air bubbles toward the central zone, the solid region central core of the fluidized bed is narrower. This mixing behavior can be analyzed through qualitative results of solid-phase distribution contours in Figures 9A-2, A-3. This behavior is also depicted by Mu et al. (2020), where the solids have a higher velocity at the annular region when using low superficial gas velocities. However, as the fluidization velocity of air increased, the velocity profile of solids shifted toward the central core of the bed. By examining the contours of slotted plate distributors (b) and (c) of Figure 9, it is revealed that a larger bulk flow of air can quickly be passed through the slots due to a bigger open area ratio. It is showing a lower kinetic head as compared to the perforated plate. Therefore, less pressure was dropped across a 90° slotted plate distributor.

Moreover, the particle movement in a 90° slotted plate distributor comes out to be straight upward, having a maximum axial component. Therefore, the stagnant zone of solid concentration is high enough adjacent to the 90° slots describing the sufficient portion of dead zones at the lower portion of the bed during the fluidization phenomenon, as shown in Figure 9B. The pressure drop across the 45° slotted plate is slightly higher than the 90° slotted plate due to the restrictive angular flow through the distributor’s slots. Furthermore, it can be noticed from the contours of a solid phase (Figure 10) that there is a significant skewed movement of an airflow across the 45° slotted plate. It allows the flow path to be two-dimensional or swirling due to both axial and lateral flow components. Such flow distribution permits the better mixing of solid particles within the fluidized zone and eliminating dead zones within the fluidization regime.

Iso-surface of solid-phase distribution further explains the flow pattern inside the fluidization region. From Figure 11, it can be seen how the solid and air phases demonstrate different flow patterns by using different distributor plates. As shown in Figure 11C, a significant swirling flow pattern is observed for a 45° slotted plate distributor. This flow behavior helps to mix binary mixtures and enhance heat transfer and mass transfer rates in the hot model (McAuley et al., 1994; Yang et al., 2021).

Perforated plate distributor shows uniform flow distributor in Figure 9A. However, there are many chances of dead zones in the lower portion of the bed or near the distributor plate where the airflow direction is almost axial along the fluidization column. It can also be revealed from the qualitative analysis of CFD results using a 90° distributor plate. A large portion of the stagnant region near the lower part of the fluidization column is due to the maximum axial flow with the minimum radial current and less mixing. The only advantage for the 90° distributor plate is its lower pressure drop than that of the other two distributor plates.

Solid particle distribution along the height of the fluidization column appears to be diverse for different distributor plates. Such solid-phase distribution patterns describe the overall mixing scenario within the fluidization zone. Figure 12 explains the solid particles’ distribution to varying heights of a column using other distributor plates. It can be observed that for the perforated distributor plate, the solid-phase concentration near the distributor plate is maximum, which indicates the existence of numerous dead zones near the lower portion of the sand bed. The solid-phase concentration falls sharply as we go up the column and eventually becomes zero above the fluidization zone (Afrooz et al., 2006). Using a 90° distributor plate, the graph of Figure 12 shows the higher concentration of the solid-phase volume fraction near the plate but not as much when compared to the perforated distributor plate. This effect is due to the larger opening area of the flow to allow the particles to disperse more rapidly, causing a relatively lower concentration near the plate. Using a 90° plate distributor, the solid-phase concentration gradually falls with a more significant slope than the perforated distributor plate. Finally, the 45° distributor plate showed the minimum solid-phase concentration near the distributor plate region due to the swirling flow. This swirling flow is attributed to the slanted (45° inlet air) inlet boundary causing axial and radial flow components to develop in the lower portion of the sand bed.

As we go up the bed’s height, the solid volume fraction first increases and decreases to the fluidization column’s freeboard region. Quantitatively, the mixing behavior in perforated distributor plates exhibits an initial volume fraction of around 0.58. It falls rapidly as go up the riser (7.7% of column height); 90° slotted plate shows an initial lower volume fraction of around 0.5. Then, it exhibits an even broader sand volume fraction along with the column height (13.46% of column height). Finally, the 45° distributor plate reveals the highest range of volume fraction through the riser height (17.3% of column height), indicating the better mixing characteristics of the fluidized zone. This broad spectrum of solid-phase concentration for the column’s height indicates the better mixing characteristics of the flow inside the fluidized bed. Qualitatively, the solid-phase distribution along the fluidization column's height can be analyzed through volume fraction contours of solid-phase for each distributor plate. Figure 13 describes the particle distribution on the cut plane section at specified heights of the column. The blue color represents the pure gas phase and the red color the solid phase. By observing the perforated distributor plate solid volume fraction contours, it can be monitored that air passes through the perforated plate’s holes by observing the cut plane section near the plate. Still, as we go up the column, the blue color fades down, indicating the pressure drop of fluidizing air as it passes through the sand bed.

Moreover, the lateral component of air velocity and the axial component start to develop, which is depicted by the enlargement of air bubbles. Air passage through the 90° distributor plate slots clearly describes the larger air bubbles rising upward and expanding as they approach the bed surface. This phenomenon helps to mix the particles only on the upper portion of the bed. A significant portion of the solid phase remains almost stagnant, leading to dead zones' formation at the lower portion of the fluidized bed. Solid volume fraction contours of the 45° distributor plate in Figure 13 describe the two-dimensional flow phenomena, i.e., the axial flow, which is along with the column's height, and the lateral or side-wise flow parallel to the radial direction. In this way, as the air enters the 45° distributor plate, the velocity of air splits up into two main components: one is parallel to the riser height (vsin θ ) and the other component of velocity is perpendicular to the inlet flow (vcos θ ), where θ is the flow angle, which varies between zero and 90°. Due to this two-dimensional flow, the resultant flow becomes swirling as it goes up the fluidization column (Shukrie et al., 2016; Yudin et al., 2016). Such flow reduces the stagnant region within the fluidization regime, as shown in Figure 13, which shows the contours of a 45° distributor plate depicting the swirling flow pattern of solid-phase as going up the column. Batcha et al. (2013) have also obtained a similar swirling flow trend after carrying out the CFD simulations series using inclined blade distributors with different blade angles. The optimum swirling flow will combine the minimum pressure drop and the high uniformity of solid particles to ensure better mixing; such phenomena are also illustrated in our findings using a 45° distributor plate.