Effects of Inlet Velocity Profile on the Bubble Dynamics in a Fluidized Bed Partially Filled with Geldart B Particles

Abstract

In this study, a two-dimensional computational domain featuring gas and solid phases is computationally studied for Geldart-B-type particles. In addition to the baseline case of a uniform gas-phase injection velocity, three different inlet velocity profiles were simulated, and their effects on the fluidized bed hydrodynamics and bubble dynamics have been studied. An in-house computer program was developed to track the bubbles and determine the temporal evolution of their size and position prior to their breakup. This program also provides information on the location of bubble coalescence and breakup. The gas-solid interactions were simulated using a Two-Fluid Model (TFM) with Gidaspow’s drag model. The results reveal that the bed hydrodynamics feature a counter-rotating vortex pair for the solid phase, and bubble dynamics, such as coalescence and breakup, can be correlated with the vortices’ outer periphery and the local gradients in the vorticity.

1. Introduction

Gas-solid fluidized beds have varied applications in the chemical, pharmaceutical, and food industries. For instance, fluidized beds find applications in fluidized catalytic cracking (FCC) in the chemical industry, wet granulation in the pharmaceutical industry, and food-processing techniques, such as fluidized bed freezing, drying, and fermentation [1,2,3]. They have gained traction in the clean energy sector as well especially in the biomass energy industry [4,5]. Fluidized beds are also implemented in concentrated solar power plants (CSP plants), where the solar thermal energy is transferred to particles inside a solar receiver which stores thermal energy. In order to enhance the heat transfer from the particles, gas-solid fluidization techniques have been implemented in the past [6,7]. Fluidization techniques are also employed in particle-to-supercritical carbon dioxide (sCO₂) heat exchangers in CSP plants to enhance the heat transfer between the moving packed bed of particles and the serpentine arrangement of circular tubes carrying sCO₂ [8].

Fluidization results in enhanced mixing between the two phases through complex bubble dynamics and gas-solid motion in its vicinity. The gas-solid mixing strongly depends on the size of the formed bubbles as well as their coalescence and/or breakup characteristics. The bubble dynamics depends on the particle size distribution, fluidization conditions, and geometry of the fluidized bed. For a superficial gas velocity which is slightly above the minimum fluidization velocity for the bed, Geldart B particles (150–1000 microns) exhibit bubbling [9]. These bubbles effectively mix the particles and enable better heat and mass transport. However, the initial fluidization regimes just above the minimum fluidization velocity display either a slugging motion of large bubbles or smaller bubbles coalescing to form larger bubbles as they find their way to the top of the bed. These large bubbles are detrimental to the desired mixing effects of fluidization, which is a phenomenon observed in [10,11,12,13]. The wake of a bubble enhances mixing as it generates small-scale vortices in the trailing regions of its path [14]. In a study by Askarishahi et al. [15], two counter-rotating vortex pairs (CRVPs) were found to be generated with one pair at the bottom and the other at the top. The size of the lower CRVP was reduced while that of the upper pair elongated axially (in the direction of the bulk motion of the fluid phase) with an increase in the fluidization velocity. With an increasing aspect ratio of the bed, the hydrodynamics in the lower bed remained largely unaffected; however, the upper bed vortices elongated axially.

Over the past few decades, several methods have been explored to mitigate the generation of large-scale bubbles, with a common approach involving the introduction of solid structures inside the bed, which results in the breakup of large bubbles into smaller bubbles [10,16]. For instance, Köksal and Vural [17] implemented a double-plate distributor at the bed inlet, which led to a pulsating flow. The frequency of the motion of plates resulted in a certain inlet air pulsation frequency, and this flow conditioning resulted in a decrease in bubble size. Kleijn van Willigen et al. [11] implemented a body-force augmentation technique by applying an electric field to control the size of the bubbles. The general trend observed in the literature is that bubbles coalesce toward the top of the bed [18].

Wang et al. [19] proposed a bubble coalescence and breakup model which could predict their frequency of occurrence along with the bubble size distribution in the bed. An investigation conducted by Wang et al. [20] agreed with the proposed correlation in which the dynamics of coalescence and breakup were studied by implementing three-dimensional bubble-tracking techniques to characterize the frequency of the occurrence of these unique phenomena particular to fluidized beds. It was observed that when the frequencies of bubble coalescence and breakup matched, an equilibrium time-averaged bubble diameter was achieved. In a similar attempt focused on the understanding of bubble dynamics, Choi et al. [21] proposed a correlation between the bubble size distribution and the frequency of coalescence and breakup; Wytrwat et al. [22] proposed a correlation between bubble holdup, velocity, and vertical bubble length. It was found that the shape of the bubbles had a strong correlation with the superficial gas velocity. A critical review by Karimipour and Pugsley [23] provides a list of correlations which relates the inlet velocity with the bubble size distribution, coalescence, and breakup frequencies.

Bubble dynamics in fluidized beds has also been studied by a combined DEM-CFD approach, for instance, by Zhou et al. [24]. DEM-CFD is a computationally expensive approach where the flow interaction with each particle is captured along with the particle-particle and particle-wall interactions. Typically, particle sizes greater than 500 microns (Geldart B particles) have been simulated through the DEM-CFD approach [25,26,27]. Another approach which is less computationally expensive involves the modeling of the two phases (solid and gas) as an interpenetrating continuum, where drag models are used to capture the gas-solid interactions; this method is commonly known as a Two-Fluid Model (TFM) [28]. The Syamlal-O’brien model [29] and the Gidaspow drag model [30] are the ones which are commonly implemented in TFM studies.

To accurately predict the gas-solid dynamics in a fluidized bed through TFM modeling, several different input parameters are required, such as the solid-solid coefficient of restitution [31,32], solid-wall specularity coefficient [33,34,35], and, finally, particle diameter [36,37]. Furthermore, an appropriate selection of a gas-solid drag model is required to accurately capture the interphase drag forces. The Gidaspow drag model and Syamlal-O’Brien drag models are some of the earliest models which have been proven to have accurate predictive capabilities in several investigations [38,39,40,41,42,43]. In an investigation by Jung et al. [44], a fluidized bed was experimentally analyzed to obtain the granular temperature, solid velocity field, and bed height. The same research group performed a numerical study in which the TFM results were validated with the experimental data [9]. Two codes, namely, MFix and IIT, were used in the study to compare their predictive abilities.

In bubbling fluidized beds, the mixing of bed elements is achieved through agitation as the bubbles find their way through the bed in the direction of a favorable net pressure gradient. In order to achieve the desired bed mixing results, it is imperative to understand the fundamental bubble dynamics in fluidized beds through probing into the fluidized bed domain, visualization, and analysis of the gas-solid flow field. With respect to a balance between numerical prediction accuracy and computational cost, the TFM is a preferred choice for modeling the gas-solid flow in fluidized beds. This modeling approach, in particular, is well suited for capturing the bubbles in bubbling fluidized beds [9]. As mentioned earlier, one of the methods used to enhance the mixing of fluidized beds is through the introduction of obstacles such as pins (or cylinders) inside the bed, which results in the reduction in bubble size, thus affecting the mixing characteristics. The introduction of obstacles also provides control over the resultant heat transfer, as they can be strategically arranged to achieve the desired performance of the fluidized beds. However, it has been found that such structures when placed in the fluidized bed have a low lifetime, as they are subjected to a volatile environment which includes bubble formation and their breakup, resulting in the random motion of solid particles and leading to erosion over time. Furthermore, the introduction of obstacles results in a net reduction in the available volume for packing the solid medium in the total fluidized bed volume [10,16,17]. The objective to achieve the desired mixing results in fluidized beds with the economical usage of fluidized gas requires the development of cost-effective solutions which have long lifetimes. The present study is motivated by this aim, where we propose inlet flow conditioning to achieve the desired fluidization properties. Different inlet velocity profiles can be achieved by installing variable porosity mesh at the entrance of the fluidized beds, which under the action of uniform plenum supply pressure would result in specific inlet velocity profiles. This approach has its engineering relevance, as it is an easy modification to the existing fluidized beds, while at the same time, it provides flow control options and has a superior lifetime. To this end, it is imperative to understand the fluidized bed behavior in response to different inlet velocity profiles, such that the fundamental understanding can be leveraged to provide design recommendations for superior mixing properties of fluidized beds.

There have been a few studies which were focused on understanding the effects of the superficial gas velocity magnitude and bed aspect ratio of the fluidization characteristics. These characteristics have been studied through the statistical analysis of bubble velocity and size distributions [23]; however, the fundamental understanding of the complex bubble dynamics is still limited, more so, from the perspective of bubble control. This study is focused on understanding the effects of different inlet velocity profiles on fluidization properties through a novel bubble tracking and analysis method for TFM-based modeling.

The following sections provide details of the mathematical modeling, computational domain, flow conditions, and bubble tracking methodology, which is followed by results and their discussion.

2. Mathematical Modeling

2.1. Governing Equations

The fluidization process has been studied via a Two-Fluid Model (TFM) following the Eulerian-Eulerian approach, where the granular phase is modeled as a fluid. The gas and solid phases are modeled as an interpenetrating continuum. The flow governing equations are presented below [40,45]:

We used the following continuity equations:

$$\frac{\partial }{\partial t}\left({\rho }_g{\alpha }_g\right)+\nabla .\left({\alpha }_g{\rho }_g\overrightarrow{v_g}\right)=0$$

(1)

$$\frac{\partial }{\partial t}\left({\rho }_s{\alpha }_s\right)+\nabla .\left({\alpha }_s{\rho }_s\overrightarrow{v_s}\right)=0$$

(2)

In the above equations, $\alpha$ represents the volume fraction, the subscripts ‘g’, ‘s’, represent the gas and solid phases, respectively, $\rho$ represents the density, $\overrightarrow{v}$ represents the velocity, and ‘t’ represents time.

In addition, we used the following momentum equations:

$$\frac{\partial }{\partial t}\left({\alpha }_g{\rho }_g\overrightarrow{v_g}\right)+\nabla .\left({\alpha }_g{\rho }_g\overrightarrow{v_g}\overrightarrow{v_g}\right)=-{\alpha }_g\nabla p+\nabla .{\overline{\overline{\tau }}}_g+{\alpha }_g{\rho }_g\overrightarrow{g}+K_{gs}\left(\overrightarrow{v_s}-\overrightarrow{v_g}\right)$$

(3)

$$\frac{\partial }{\partial t}\left({\alpha }_s{\rho }_s\overrightarrow{v_s}\right)+\nabla .\left({\alpha }_s{\rho }_s\overrightarrow{v_s}\overrightarrow{v_s}\right)=-{\alpha }_s\nabla p-\nabla p_s+\nabla .{\overline{\overline{\tau }}}_s+{\alpha }_s{\rho }_s\overrightarrow{g}+K_{gs}\left(\overrightarrow{v_g}-\overrightarrow{v_s}\right)$$

(4)

In the above equations, ‘p’ represents pressure, $\overline{\overline{\tau }}$ represents shear stress tensor, $\overrightarrow{g}$ represents acceleration due to gravity, and $K_{gs}$ represents the drag coefficient (discussed later).

The shear stress tensors are given as the following:

$${\overline{\tau }}_g={\alpha }_g{\mu }_g\left(\nabla \overrightarrow{v_g}+\nabla {\overrightarrow{v_g}}^T\right)-\frac{2}{3}({\alpha }_g{\mu }_g.\left(\nabla .\overrightarrow{v_g}\right)I)$$

(5)

$${\overline{\tau }}_s={\alpha }_s{\mu }_s\left(\nabla \overrightarrow{v_s}+\nabla {\overrightarrow{v_s}}^T\right)-{\alpha }_s\left({\lambda }_s-\frac{2}{3}{\mu }_s\right)\left(\nabla .\overrightarrow{v_s}\right)I$$

(6)

In the above equations, $\mu$ is the dynamic viscosity, and $\lambda$ is the bulk viscosity.

We also used the following granular temperature equation:

$$\frac{3}{2}\left(\frac{\partial }{\partial t}\left({\alpha }_s{\rho }_s{\theta }_s\right)+\nabla .\left({\alpha }_s{\rho }_s\overrightarrow{v_s}{\theta }_s\right)\right)=\left(-p_s\overline{\overline{I}}+{\overline{\overline{\tau }}}_s\right):\nabla \overrightarrow{v_s}+\nabla .\left(k_{{\theta }_s}\nabla {\theta }_s\right)-{\gamma }_{{\theta }_s}+{\varphi }_{gs}$$

(7)

In the above equation, ${\theta }_s$ is the granular temperature, ${\gamma }_{{\theta }_s}$ is the collisional dissipation of energy, and ${\varphi }_{gs}$ is the energy exchange between the gas and solid.

In Equation (5), $p_s$, ${\mu }_s$ and ${\lambda }_s$ can be determined as functions of the granular temperature, radial distribution function, solid volume fraction, and the thermo-physical properties of the solid phase. These relationships are comprehensively reviewed in [45]. Granular temperature ${\theta }_s$ is not the same as the conventional temperature term. ${\theta }_s$ arises from the kinetic theory of granular flows just as temperature does from the kinetic theory of gases. This quantity measures the average rate of change in velocities of the granular phase at a fixed point in time. Granular temperature is the quantification of random oscillations of particles. Mathematically, it is the average of three variances in three dimensions of the instantaneous velocities of the particles [44].

Since the gas-solid flow is modeled considering the solid phase as a continuum, the individual solid particles were not resolved and tracked. Further, the solid-solid interaction was modeled through the granular shear and bulk viscosity. This is a complex process, as the solid phase, in reality, is arranged in a discrete manner. Hence, an estimation of the net shear experienced by particles based on their local distribution is required. This phenomenon is captured through statistical methods by formulating a radial distribution function, which is the probability of finding other particle(s) in radially outboard directions in reference to a fixed particle. In the present study, the radial distribution function formulated by Ahmadi and Ma [46] was implemented.

In gas-solid flows, two types of flow regimes can be observed when attempting to close the momentum equations. These are plastic flow and viscous flow. In gas-solid flows, the solids collide with one another and exchange momentum. In plastic flows, momentum transfer due to friction dominates, while in viscous flow, the momentum transfer due to collisions dominates. More details are provided in [45]. Since the inlet flow velocity used in the present study was only slightly higher than the minimum fluidization velocity, the frictional stresses were ignored and only solid-phase pressure was used as the flow was viscous, where momentum transfer is dominated by collisions compared to friction. For the viscous stresses (momentum transfer via collisions), the formulation given by Lun et al. [47] was implemented.

2.1.1. Drag Models

The solid-solid interaction was modeled in MFix, following the methodology prescribed by Syamlal et al. [45]. As mentioned earlier, several gas-solid drag models have been developed in the past, among which the Syamlal-O’Brien and Gidaspow models are widely used [29,30]. The drag models are constructed to yield the drag coefficient, $K_{gs}$, which is included in the momentum conservation equations. The Syamlal-O’Brien drag model is given as the following:

$$K_{gs}=\frac{3}{4}\frac{{\alpha }_s{\alpha }_g{\rho }_g}{v_{r,s}^2{d}_s}{C}_d\left(\frac{R{e}_s}{v_{r,s}}\right)\left|\overrightarrow{v_s}-\overrightarrow{v_g}\right|$$

(8)

where $v_{r,s}$ is the terminal velocity and is given by

$$v_{r,s}=0.5(A-0.06R{e}_s+\sqrt{{\left(0.06R{e}_s\right)}^2+0.12R{e}_s\left(2B-A\right)+A^2})$$

(9)

where

$$A={\alpha }_g^{4.14}$$

and

$$B=0.8{\alpha }_g^{2.65}$$

for ${\alpha }_g$ $\le 0.85$ and $B=0.8{\alpha }_g^{2.65}$ for ${\alpha }_g$ > 0.85

The Gidaspow model is given as

$$K_{gs}=150\frac{{\alpha }_s^2{\mathsf{\mu }}_g}{{\alpha }_g^2{d}_s^2}+1.75\frac{{\alpha }_s{\rho }_g\left|\overrightarrow{v_s}-\overrightarrow{v_g}\right|}{{\alpha }_g{d}_s}\mathrm{ }\mathrm{ }For\mathrm{ }{\alpha }_g\le 0.8$$

(10)

$$K_{gs}=\frac{3}{4}{C}_d\frac{{\alpha }_s{\rho }_g\left|\overrightarrow{v_s}-\overrightarrow{v_g}\right|}{d_s}{\alpha }_g^{-2.65}\mathrm{ }\mathrm{ }For\mathrm{ }{\alpha }_g>0.8$$

(11)

where

$$C_d=\frac{24}{Re}\left[1+0.15{\left(R{e}_s\right)}^{0.687}\right]; R{e}_s\le 1000$$

(12a)

and

$$C_d=0.44; R{e}_s>1000.$$

(12b)

where

$$R{e}_s=\frac{{\alpha }_g{\rho }_g{d}_s|\overrightarrow{v_s}-\overrightarrow{v_g}|}{{\mu }_g}$$

(13)

2.1.2. Gas Turbulence Modeling

In MFiX, the gas turbulence is modeled as a combination of the gas viscosity and the eddy viscosity, which is the Boussinesq hypothesis [48] for turbulent flows. In MFiX, eddy viscosity was set to zero by default, implying that the stress tensors were purely a function of the gas molecular viscosity. In the present study, turbulence was not modeled as it is established that in gas-solid flows (and for conditions similar to those studied here), a dense regime of the solid phase exists, and the contribution of gas turbulence is negligible [9]. In an investigation by Loha et al. [40], the effect of turbulence modeling was studied in contrast with the laminar flow model, where no appreciable differences were found between the two modeling approaches.

3. Numerical Considerations

The computations were carried out in a 2D domain of size 15.5 cm (width) and 40 cm (height). The domain was discretized with quad elements with 31 divisions along the width (resulting in 32 cells) and 82 divisions (resulting in 83 cells) along the height of the bed. The discretized computational domain is shown in Figure 1. The recommendations on the grid size for the two-fluid model simulation as prescribed in [9,40,49] is adopted in this study, where the size of one cell was ~10 times the particle diameter.

The continuity and momentum residuals were set to 0.001 with the temporal discretization scheme of Implicit Euler and spatial discretization of MFiX-Superbee. The linear solver was of the bi-conjugate gradient stabilization type. The maximum time-step size was 0.01 s, and the minimum time-step size was 1 × 10⁻⁶ s. For the momentum formulation, the entire framework of equations and methods listed so far belongs to the class of Model A [50]. The simulations were conducted for a physical time of 40 s. Each of the simulations took nearly two hours to complete, and the simulations were conducted via the Dell Precision 7865 Tower, which has an AMD Ryzen Threadripper PRO 5945WX processor, 12 cores, 32 GB RAM.

4. Gas and Solid Properties, Initial Conditions, and Boundary Conditions

Air was modeled as an incompressible fluid. The solid phase was modeled as a granular phase and treated as a continuum along with the gas phase. The particle diameter was 530 microns, and the density of particles was 2500 kg/m³. The particle-particle coefficient of restitution was set to 0.99. These properties are summarized in

Table 1. Properties of the fluid medium.

Property	Value
Air viscosity	1.5 × 10−5 Pa.s
Air density	1.2 kg/m3
Solid-Wall-e	1
Solid-Solid-e	0.99
Particle diameter	530 microns
Particle density	2500 kg/m3
Operating pressure	101,325 Pa
Gravity	9.81 m/s2
Particle-particle mean shear	Granular temperature-based
Particle-particle bulk mean shear	Granular temperature based

. In addition to these parameters, the body force due to earth’s surface gravity (9.81 m/s²) was modeled.

To facilitate the convergence during the initial physical times, the initial shock to the entrance region of the bed had to be reduced. To this end, the initial bed height was set to 20 cm with a void fraction of 0.6. A no-slip wall boundary condition was applied for the gas phase. For the solid phase, the boundary treatment was based on the formulation provided by Johnson and Jackson [51] with a specularity coefficient of 0.6. The coefficient of restitution for the wall-particle pair was set as unity. For the baseline case, a constant velocity (0.587 m/s) inlet boundary condition was applied at the inlet, with a void fraction of unity, implying that only air entered the computational domain. The outlet boundary was set as Neumann boundary type with a fixed atmospheric pressure of 101,325 Pa.

The effect of the inlet velocity profile on the hydrodynamics of the fluidized bed is studied via a comprehensive scientific procedure in this paper. Firstly, the TFM is extensively validated with prior benchmark experimental studies. Second, a novel bubble identification and tracking algorithm was developed and validated with experimental data. The numerical solutions were then studied in conjunction with the bubble dynamics to understand the fundamental relationship between the fluidized bed hydrodynamics and the location of bubbles and their size distribution. The mixing characteristics of the fluidized bed for different inlet velocity profiles were studied through the above approach. The primary hypothesis behind this investigation is that the fluidized bed mixing characteristics can be controlled by conditioning the inlet gas velocity, as the resultant solid and gas phase motion in the fluidized bed as well as near-wall hydrodynamics will be affected by the inlet conditions. Further, through the prior knowledge of the relationship between the bubble dynamics and the localized mixing, the above-mentioned approach is adopted to test the hypothesis through a comprehensively validated numerical model.

For the other inlet velocity profiles investigated in this study, the profiles were entered via a subroutine written in FORTRAN and compiled in MFiX. In this study, four different velocity profiles were studied, as shown in Figure 2. For all profiles, the average velocity ($u_0$) was kept the same (0.587 m/s). The four velocity profiles are also summarized in

Table 2. Test matrix.

Configuration	Test Condition
Baseline	0.587 m/s uniform inlet velocity
Parabolic	Refer Figure 2 for velocity profile
Trapezoidal 1	Refer Figure 2 for velocity profile
Trapezoidal 2	Refer Figure 2 for velocity profile

.

5. Bubble Tracking Methodology

The present study is focused on understanding the fluidized bed properties through analyzing the bubble dynamics by tracking the void fraction results as predicted by the TFM. The post-processing of the predicted results on the gas-solid dynamics involved a bubble tracking approach, which was implemented through the development of an in-house computer program. The solid volume fraction contours were generated at different time steps during the course of physical run time, and these contours were exported from MFiX to a separate folder in jpeg format. The exported images were in grayscale mode where the solid phase was represented by pixels with higher intensity, and gas was represented by those with lower intensity (first row of Figure 3). Upon the visual examination of individual images of the solid-volume fraction in the grayscale mode with trial and error, a threshold value of 0.8 was considered appropriate to mark a clear difference between the gas and solid phases in the fluidized bed. The grayscale images exported from MFiX were then imported in MATLAB R2022.a for further processing based on the above criterion to facilitate the bubble dynamics study. In order to obtain a well-defined bubble geometry which had clear borders such that an exact coordinate of its centroid could be determined, the threshold criteria of 0.8-pixel intensity was applied uniformly across all the solid-volume fraction images exported in the grayscale format (second row of Figure 3). The next step in this process is illustrated for a sample bubble inscribed in a red box as shown in the first row of Figure 3. An in-built function in MATLAB R2022.a was used to identify the boundaries of a region marked by a change from 0 to 1 or 1 to 0. The region boundary information as obtained from the above step in pixel notation was then converted to physical length scales to account for the physical area corresponding to a pixel and its location in the 2D computational domain. Once the bubble boundary information was stored, the information on the location of the corresponding centroids of different bubbles was stored for a given binary scale image. Figure 3 shows the evolution of bubbles going from left to right, and this evolution is effectively captured by the above bubble tracking methodology, as illustrated in the third row of Figure 3.

The stored information on bubble boundaries and centroids was then used for tracking individual bubbles and to study crucial phenomena of coalescence and breakup.

To facilitate a continuous evolution of bubbles, the physical time step size of 0.01 s was used, which ensured that the shift in the location of centroid ($∆x,∆y)$ of a given bubble was small, aiding in the accurate tracking of bubbles in both space and time. This information was used to determine the average bubble size. The next phase of bubble tracking involved a frame-by-frame sweep where the above information was tracked to identify the locations of bubble coalescence and breakup. The bubble coalescence and split locations were also monitored and stored. The algorithm for the bubble dynamics study is shown through a flowchart in Figure 4.

6. Results and Discussion

This section presents the validation of the computational model, which is followed by results and discussion on the gas-solid interactions where fluidized bed hydrodynamics is presented along with the bubble dynamics study. The interphase Nusselt number as derived from the time-averaged flow-field is also presented.

6.1. Validation of the TFM Predictions with Experiments

The prediction results on granular temperature and axial velocity have been compared with the experimental results from Jung et al. [44], where the authors in [44] conducted fluidization experiments on a thin rectangular bed which was made from an optically clear material. A uniform fluidization was achieved through a gas distributor upstream of the inlet to the fluidized bed. The experiments were conducted on glass spheres with an average particle diameter of 530 μm. The instantaneous velocity of the particles was determined by capturing the images using a CCD camera. The particle streaks captured by the camera were used to determine the streamwise and spanwise components of the flow velocity. A detailed description of the experimental setup can be found in [44].

As part of the validation studies, the drag models proposed by Syamlal-O’Brien and Gidaspow were used to compare the respective predictions through the TFM with experimental data. For comparisons, important hydrodynamic properties of the fluidized bed such as the bed height, time-averaged axial velocity of the solid phase, and granular temperature were considered. The inlet profile for validation studies was set to 0.587 m/s acting uniformly across the inlet. The presented results were averaged over time, where data over 10–40 s of physical time were considered. The initial physical times were omitted from the time averaging. During the initial physical time (0–5 s), a symmetric fluidization pattern was observed, which was due to the initialization of the hydrodynamic variables, resulting in symmetric patterns in volume fraction distribution. Post the initial timeframe which featured symmetric flow behavior, the effects of the hydrodynamic disturbances started to mark their presence in the fluidized bed. The time-averaging choice of the 10–40 s timeframe was made such that the typical characteristics of the disturbed fluidized bed can be captured. A comparison of the time-averaged axial velocity profiles with the experimental profiles is shown in Figure 5. Both the drag models agreed well with the experimental data on the axial velocity profile sampled at a bed height of 0.14 m. The second comparison was based on the granular temperature as shown in Figure 5. In this case as well, both the drag models agreed well with the experiments with the Gidaspow model being relatively closer to the measurements.

Apart from the line variation of flow quantities, the model validation study includes a localized comparison of the time-averaged solid-phase volume fraction superimposed by the solid-phase velocity vectors, which is shown in Figure 6. The Gidaspow model showed a better agreement in terms of the fluidized bed height comparison with the experimental measurements. An axially stretched large counter-rotating vortex pair (CRVP) was observed from the solid-phase velocity vectors in both Gidaspow and Syamlal-O’Brien drag model-based predictions as well as in the results in Jung et al. [9]. This CRVP resulted in a predominantly vertically downward movement of particles along the walls and an upwash of solid particles in the vicinity of the bed centerline along the bed height.

Further, the grayscale images obtained from the solid-phase volume fraction predictions using Syamlal-O’Brien and Gidaspow models are compared at two different physical time instances with the experimental data (Figure 7). It can be observed that both drag models were able to capture the intricate features of the larger bubble observed in the experiments. In summary, the validation study reveals that both drag models, Syamlal-O’Brien and Gidaspow, have good predictive capabilities for the simulated condition in reference to the available experimental data. In this study, we have used Gidaspow’s drag model to study the inlet velocity profiles’ effect on the fluidized bed hydrodynamics.

6.2. Effects of Inlet Velocity Profile on the Bed Hydrodynamics

In reference to the uniform inlet velocity boundary condition, three different velocity profiles have been studied, as shown in Figure 2. The average velocity for all the velocity profiles was kept the same as the magnitude for the uniform inlet velocity case. The time-averaged axial velocity variation with the span of the fluidized bed at 0.14 m height, for the four different inlet velocity profiles, is shown in Figure 8.

For the parabolic inlet velocity profile for the gas phase, the solid-phase axial velocity observed at a bed height of 0.14 m was also similar to a parabolic shape, with a velocity peak occurring near the bed centerline, with a peak velocity magnitude of 20 cm/s. After comparing the trapezoidal velocity profiles with the uniform inlet velocity profile, it is observed that the trapezoidal profile resulted in an even closer agreement with the experiments by Jung et al. [44]. The fluidized bed characteristics corresponding to different inlet velocity profiles are shown through the time-averaged solid volume fraction superimposed by solid-phase velocity vectors (also time-averaged) in Figure 9.

The trapezoidal inlet velocity profiles resulted in a CRVP, which was stretched axially in the direction of favorable net pressure gradient. The shape of the individual vortices in the CRVP was more well defined in the trapezoidal 2 profile, where its shape can be characterized as rectangular (following the bed domain). Note that for the trapezoidal 2 profile, it may appear that a smaller CRVP exists near the bed inlet; however, upon closer examination, it can be observed that smaller vortices are a part of the larger vortices on either side of the bed centerline (in the axial direction). As a result of large-scale CRVP occupying almost the entirety of the bed height, the solid volume fraction was found to be concentrated near the walls for the trapezoidal inlet velocity profiles (1 and 2).

For the parabolic inlet velocity profile, a well-defined oval-shaped CRVP was observed based on the solid-phase velocity vectors. This resulted in a strong upwash of the particles in the vicinity of the channel centerline, which eventually led to an even higher concentration of the solid phase near the walls, especially in the corners near the bed inlet. Upon comparing the three inlet velocity profiles, the trapezoidal 1 profile resulted in the most uniform distribution of the solid phase in the bed domain.

6.3. Effects of Inlet Velocity Profile on the Bubble Dynamics

The effects of inlet velocity profile on the bubble dynamics were studied following the methodology described in Section 5. Figure 10, Figure 11, Figure 12 and Figure 13 show the locations of bubble coalescence and breakup (splitting) superimposed on the solid-phase vorticity and solid-phase velocity vectors.

The bubble coalescence and breakup locations were found to be at the periphery of high vorticity zones or at the edges of the vortices toward the bed centerline. As shown earlier, a large CRVP for the solid phase was observed for all velocity profiles. The bubble coalescence and breakup incidences were higher in number in the vicinity of the two interacting vortices near the bed axial centerline or in the regions with high local gradient in the vorticity values. It is noteworthy that most bubble dynamics occurred away from the wall.

For the two trapezoidal inlet profiles, due to their axially stretched nature of vortices, a small high vorticity region was observed near the entrance. This behavior was not observed for the parabolic inlet profile, which resulted in an oval-shaped CRVP for the solid phase, and as a result, the vorticity magnitudes were localized in an oval shape. Most of the bubble dynamics (coalescence or breakup) were observed in between these vortices or in the regions of high local gradient in vorticity magnitudes. Near-wall bubbles were nearly nonexistent for the parabolic inlet velocity profile.

The bubble split locations were more spread out in reference to the coalescence locations for the parabolic inlet velocity profile, with more activity toward the entrance region, when compared to the other three velocity profiles. This is attributed toward the high velocity gradient near the channel centerline for the parabolic inlet profile, which resulted in increased flow disturbance when gas interacted with the solid phase.

As noted earlier, the effect of inlet velocity profile on the bubble dynamics was studied here to gain fundamental knowledge about the bubble dynamics (coalescence and breakup) in relation to the gas-solid flow behavior, which was found to be well represented by the solid-phase volume fraction, solid-phase velocity magnitude and velocity vectors, and vorticity magnitudes. For all inlet velocity profiles, it was observed that bubble coalescence/breakup locations can be expected along the periphery of the CRVPs. For fluidized bed designers, it is recommended that TFM simulations can be used to determine the CRVP(s) to locate the highest possibility of bubble coalescence and breakup.

Apart from the coalescence and breakup locations, the bubble properties are also represented by its size and the size distribution in context with the local position in the fluidized bed. The time-averaged bubble size distribution is shown in Figure 14. The bubble size was normalized by the product of initial bed height and bed width. For all inlet velocity profiles, the bubbles grew in size as the flow moved upwards in the direction of the favorable pressure gradient. The largest bubble size was observed in the vicinity of the channel centerline in the flow direction. The bubble size was largest for the trapezoidal 1 profile, while the parabolic inlet velocity profile resulted in the smallest bubble size. Note that the occurrence of the largest bubbles was localized in the region where bubble coalescence and breakup incidences were not observed. This observation further strengthens the bubble dynamics approach as it accurately captures critical bubble phenomena.

6.4. Interphase Heat Transfer Characteristics

Interphase heat/mass transfer characteristics of fluidized beds are important to understand as it determines their effective performance, such as in fluidized bed dryers [3]. The locally resolved gas-solid flow field can be used with established correlations to obtain local heat transfer characteristics within a fluidized bed. The interphase heat transfer relationship in Nusselt number form (Equation (14)) provided by Gunn et al. [52] has been used in this study.

$$\begin{gathered} Nu=\left(7-10{ϵ}_g+5{{ϵ}_g}^2\right)\left(1+0.7R{e}_s^{0.2}{Pr}^{\left(1/3\right)}\right) \\ +\left(1.33-2.4{\mathsf{ϵ}}_g+1.2{\mathsf{ϵ}}_g^2\right)R{e}_s^{0.7}P{r}^{\left(1/3\right)} \end{gathered}$$

(14)

where

$$R{e}_s=\frac{d_p\left|\overrightarrow{v_s}-\overrightarrow{v_g}\right|{\rho }_g}{{\mathsf{\mu }}_g}$$

(15)

$$Pr=\frac{C_{pg}{\mu }_{\mathrm{g}}}{k_g}$$

(16)

Figure 15 presents the local Nusselt number distribution for the four configurations, which are superimposed with time-averaged solid-phase velocity vectors. The regions corresponding to local vortices had a lower interphase heat transfer, which is expected. The configuration of trapezoidal profile 1 has the highest near-wall interphase heat transfer, which makes it suitable for applications where heat is supplied through the walls. Further, in some fluidized bed applications, the uniformity in the interphase heat transfer within the bed volume is desirable, for example, in case of fluidized bed dryers, to achieve uniformity in grain drying. To this end, the trapezoidal profiles had more uniform interphase heat transfer statistically compared to the baseline and parabolic inlet velocity profiles.

7. Conclusions and Future Recommendations

A 2D bubbling fluidized bed filled with Geldart B particles was simulated and analyzed by conducting TFM simulations using an open-source code MFiX. The predictions were validated using the experimental data presented in Jung et al. [9]. The validated TFM model was used to study the effects of three different velocity profiles in addition to the uniform inlet velocity profile on the bed hydrodynamics and bubble dynamics. The following observations on fluidized bed hydrodynamics were made:

(i)

The time-averaged flow field reveals the presence of CRVP in the fluidized bed for all the velocity profiles. The vortex shapes were dependent on the inlet velocity profiles.

(ii)

The trapezoidal velocity profiles and the uniform velocity profile had similar axial velocity distribution at a bed height of 0.14 m, while the parabolic inlet velocity profile exhibited a marked difference in the axial velocity signature, where a peak was observed near the centerline of the bed in the axial direction.

(iii)

The solid-phase volume fraction was high near the walls for all the inlet velocity profiles, and the bed height was found to be insensitive toward the difference in inlet velocity profiles. This indicates that the bed height is simply a function of the average inlet velocity.

Regarding the bubble dynamics study, the following main observations were made:

(i)

The bubble coalescence and breakup phenomena were found to take place in between the region of two adjacent vortices. Bubble dynamics in the vicinity of the walls was nonexistent.

(ii)

The number of bubble coalescence and splitting incidences followed similar patterns and were found to be nearly the same in number, indicating that bubble sizes reached an equilibrium state, which was a phenomenon also observed by Wang et al. [20].

(iii)

The bubble coalescence and breakup dynamics incidences were related to the local change in the vorticity magnitudes in the vicinity of the axial centerline of the fluidized bed.

(iv)

The average bubble size was found to be largest for the case of trapezoidal profile 1, while the parabolic velocity profile had the smallest bubble size. The larger bubble sizes were observed in the vicinity of the bed top.

The inlet velocity profile resulted in distinct interphase heat transfer characteristics, where trapezoidal inlet velocity profiles resulted in more uniform interphase heat transfer distribution, and trapezoidal velocity profile 1 resulted in the highest near-wall interphase heat transfer. The inlet flow conditioning can be an effective flow control method in fluidized beds.