Capability of the TFM Approach to Predict Fluidization of Cohesive Powders.

The fluidization behavior of cohesive particles was investigated using an Euler-Euler approach. To do so, a two-fluid model (TFM) platform was developed to account for the cohesivity of particles. Specifically, the kinetic theory of granular flow (KTGF) was modified based on the solid rheology developed by Gu et al. J. Fluid Mech. 2019. The results of our simulations demonstrated that the modified TFM approach can successfully predict the formation of particle agglomerates and clusters in the fluidized bed, induced by the negative (tensile-dominant) pressure. The formation of such granules and clusters highly depended on the particle Bond number and the tensile pressure prefactor. To evaluate fluidization regimes, a set of simulations was conducted for a wide range of particle cohesivity (e.g., Bond number and tensile pressure prefactor) at two different fluidization numbers of 2 and 5. Our simulation results reveal the formation of four different regimes of fluidization for cohesive particles: (i) bubbling, (ii) bubbling-clustering, (iii) bubble-less fluidization, and (iv) stagnant bed. Comprehensive analysis of the shear-to-yield ratio reveals that the observed regime map is attributed to the competition between the shear stress and yield stress acting on the particles. The obtained regime map can be extended to incorporate the effect of dimensionless velocity and dimensionless diameter as a comprehensive fluidization chart for cohesive particles. Such fluidization charts can facilitate the design of fluidized beds by predicting the conditions under which the formation of particle agglomeration and clustering is likely in fluidized beds.

Introduction

Importance of Cohesive Particles in Industry

Small particles are of high interest for various industrial sectors because they offer high volumetric specific surface area, leading to higher contact/reaction rates per unit reactor/bed volume.[2] For instance, small particles serve as reacting particles in the chemical and petrochemical industry (e.g., pulverized coal combustion and gasification) and are used as solid catalysts for chemical reactions (e.g., fluid catalytic cracking or FCC).[3] In the pharmaceutical industry, typically, active pharmaceutical ingredients (APIs) are relatively fine to improve the dissolution and release rate in the body. Besides, excipients in the form of fine powders can improve the uniformity of the tablet.[4] Fine powders are also favorable in the food industry due to their higher solubility.[5] Fluidized bed technology can be employed to take advantage of such a high specific surface area. However, handling fine/small powders in a fluidized bed (FB) can face several challenges.

Most importantly, small and specifically fine powders categorized as group A and C in Geldart’s classification[6] are not easy to get fluidized due to the dominance of van der Waals force at small diameters. Powders featuring cohesive interaction forces can be categorized from mildly cohesive to highly cohesive. To quantify the level of cohesion, the Bond number (NBond) is typically used. This dimensionless number describes the ratio of cohesion forces and the gravitational force.

Before reaching the bubbling condition, mildly cohesive powder (Group A) typically experiences bubble-less expansion at which the bed gets fluidized uniformly, but no bubble is observed in the bed.[7,8] Such a behavior can be related to the dominance of yield stress over shearing the powder, which has been the topic of many research studies.[7−11]

A group of researchers evaluated the fluidization behavior of cohesive particles experimentally. For instance, Geldart,[8] as the pioneer in this field, reported that the interplay of hydrodynamic and interparticle force governs the FB characteristic, for example, bed expansion and minimum bubbling/minimum fluidization ratio.

LaMarche et al.[7] experimentally investigated the fluidization behavior of mildly cohesive particles for a range of fluidization velocities. Their experimental data suggest that the mean bed voidage is smaller at the minimum bubbling condition due to the dominance of the drag force over the cohesion force. In another study, Li et al.[11] experimentally investigated the fluidization behavior of Geldart A particles in a fluidized bed. According to their observation, particle cohesion results in the stickiness of a layer of particles on the wall. They also observed the formation of particle clusters in the form of agglomerates which are not permanent as they continuously form and break.

Such experimental observations for cohesive powder fluidization can improve our understanding of particle flow behavior, specifically fluidization regimes. However, experimental studies seem very challenging to obtain a deep insight into the main reason behind such a behavior. This is since one cannot evaluate the contribution of different interparticle forces through experimental investigations. For instance, it is known that the cohesiveness of a powder and its fine content are known to play a significant role in the rheological and fluidization behavior of Geldart’s group A powders.[12] Specifically, van der Waals cohesive forces can create agglomerates of small particles in a fluidized bed when dominating the applied shear stress.[3] Here, one can take advantage of the recently established computational models, which we will summarize next.

Mathematical and Computational Modeling

In the study of cohesive powder fluidization, isolating the effects of cohesive (van der Waals) forces from other effects, for example, electrostatics and humidity, is very difficult in experimental investigations.[3] Capable of predicting spatially distributed information, mathematical approaches enable us to evaluate the effect of individual forces on the flow behavior of cohesive powders. To put it in more detail, detailed numerical simulations are capable of predicting the distribution of interparticle forces and their contributions to the strength of the particle aggregates. Consequently, the flow behavior of such powders can be thoroughly analyzed. Various mathematical approaches have been used in the literature as detailed below.

A group of researchers employed coupled computational fluid dynamics (CFD) and −discrete element method (DEM) simulations to predict fluid (i.e., gas) flow and particle flow, respectively. In this approach, the Navier–Stokes equation is solved for fluid motion, while the motion of individual particles is simulated by solving the second Newton’s law of motion. For instance, Li et al.[11] used the CFD–DEM approach to evaluate the formation of particle agglomerate in a fluidized bed of Geldart A particles. They reported that DEM simulations fail to capture agglomerates raining down from the bubbles’ cap without accounting for the cohesive forces. Their CFD–DEM results revealed that highly cohesive particles tend to form large agglomerates in dense regions. This finding was also reported by Wu et al.[13] Their results[11] demonstrated that the formed agglomerates can be broken on the bubble’s cap. They validated the formation of agglomerates through an experimental study. However, they reported that agglomerates predicted in simulations are larger while having a smaller aspect ratio compared to the experimental observation. They concluded that considering only van der Waals forces in DEM is insufficient to predict the agglomerates’ properties accurately. They related the formation of elongated agglomerates to gas–particle interactions in the presence of cohesive forces, which requires direct numerical simulation. This means that revisiting the fluid–particle interaction is also needed when simulating a fluidized bed of cohesive particles.

In another CFD–DEM study, Liu et al.[14] reported that the defluidization curve of Geldart A particles depends on the static bed height and Young’s modulus. Their main purpose was to evaluate if the defluidization curve of cohesive powders can be predicted by simulating soft particles (i.e., those featuring low values of Young’s modulus) with the DEM. This was mainly done to reduce the computational cost of a DEM-based simulation. According to Liu et al.,[14] variation of Young’s modulus can enhance cohesive effects in DEM-based simulations simply because the particle–particle overlap is different (which is an input to the cohesion model). They also examined the effect of the static bed height to obtain a system size-independent measurement for cohesive powders. They attributed this behavior to the interaction of cohesive forces with the packing structure during defluidization and material stiffness. Fluidization/defluidization behavior of aeratable Geldart’s group A powders was also investigated by Galvin and Benyahia.[3] Their DEM simulation results demonstrated a significant effect of sliding friction in the presence of cohesion. Specifically, they reported that neglecting the effects of cohesion marginalizes the action of friction for small type A particles. Recently, Wu et al.[13] used the CFD–DEM method to connect the microscopic discrete properties of cohesive powder to the macroscopic continuum description of a cohesive fluidized bed. Their simulation results showed that the tensile pressure predicted by CFD–DEM is in accordance with the Rumpf correlation.[15]

Another group of researchers tried to use microscale simulation (e.g., DEM) to develop closures for a continuum approach (e.g., the two-fluid model, TFM) to simulate cohesive particle fluidization. In the TFM approach, fluid and particles are treated as interpenetrating continua, and the momentum balance equations for the fluid and solids phase are solved. As solid particles are considered as one phase, the computational cost for such an approach is much lower when compared to the CFD–DEM approach, especially for large-scale applications. Therefore, academia and industry have a high interest in equipping the TFM approach with suitable cohesion models.

In the TFM, different approaches have been adopted to account for the effect of cohesion force. The first approach is to directly model the particle growth through, for instance, the population balance equation as suggested by Kellogg et al.[16] They developed a continuum approach for rapid cohesive-particle flow using such a population balance method. They derived a closure for the success factor of cohesive powders’ collision. Their closure relates the success factors of collision to the granular temperature, a measure of the particles’ impact velocity. They also derived a closure for the effective restitution coefficient, which accounts for the dissipation of energy due to cohesion. Their model contains several model parameters related to agglomeration’s critical velocities, which can be defined using simple discrete element method simulation.

In the second approach, cohesive forces were accounted for in a balance equation of force acting on a single solid particle, as proposed by van Wachem.[17] Based on this approach, he locally calculated the cohesive agglomerate size, which can be linked to the TFM simulation through solids phases of different diameters, representing primary particles and granules of different sizes. Therefore, the number of solids phases is defined based on the specified bin size in the particle size distribution (PSD).

In the third approach, modified excess compressibility is used to account for the effects of cohesion between particles. Ye[18] investigated the cohesion effect using discrete particle simulations. Subsequently, they modified the classical KTGF theory according to their soft-sphere discrete parcel method (DPM) simulation results. In detail, they modified the KTGF considering the modified excess compressibility, which accounts for the effects of cohesion between particles. This excess compressibility depends on the magnitude of the cohesive force and the solids volume fraction.

In the fourth approach, as the most robust approach, the rheology of the solids phase was directly modified to account for cohesion forces, as adopted first by Gidaspow and Huilin.[19] They included a cohesion pressure in their solid pressure model via the radial distribution function used in their KTGF formulation. Very recently, Gu et al.[1] used the CFD–DEM approach to derive a continuum model for the rheology of mildly cohesive powders. Specifically, they also modified the KTGF theory to account for particle cohesion. Through a CFD–DEM simulation, they modified the constitutive equations for the solids pressure, solids bulk viscosity and shear viscosity, and the rate of dissipation of pseudo-thermal energy.

Gaps in the Literature

From the abovementioned literature, one can conclude that fluidization/defluidization behavior of Geldart A group particles (i.e., mildly cohesive powder) is not well-understood due to several limitations:

Experimental evaluation of forces seems very difficult as the effect of force acting on the particle cannot be isolated. Therefore, it seems impossible to explain the reason behind the bubble-less expansion in an FB of Geldart A particles. Mathematical and computational approaches can be a promising option to tackle this issue.

Lagrangian approaches cannot be applied for large-scale FBs due to the limitation imposed by computational time. Alternatively, less computationally expensive tools such as TFM can be used.

The standard KTGF fails to predict agglomerate formation in an FB of cohesive particles as this approach has been developed for slightly inelastic particles and does not account for the cohesion force. However, the collisional behavior of particles is influenced by interparticle forces.

Some researchers coupled the TFM approach with the population balance method to account for agglomerate growth. Nonetheless, the population balance method is associated with several model parameters which are valid for the conditions of establishment only. This limits the applicability of this method. This means that the TFM approach needs to be equipped with robust rheological models for the solids phase to account for the particle cohesiveness, as recently conducted by Gu et al.[1]

The effect of the Bond number on the fluidization behavior of cohesive powders is not thoroughly investigated in the open literature.

A detailed regime map of fluidization for Geldart A particles is lacked in the literature. Depending on the level of cohesion, one can expect different regimes of fluidization for mildly cohesive powder: (i) bubble-less expansion, (ii) bubbling, (iii) bubbling–clustering, and (iv) static bed.

Goal and Outline

By considering the gaps in the literature addressed in the previous section, we distilled the following main goals of our present study:

Evaluating the capability of the modified TFM approach developed by Gu et al.[1] in predicting (i) the fluidization behavior of cohesive powders and (ii) the formation and breakage of granules. To realize this goal, Gu’s rheological model[1] was implemented in a TFM simulation platform. Subsequently, the fluidization behavior of cohesive powder was evaluated for different levels of cohesivity in terms of particle Bond number and tensile pressure prefactor (αcoh,1) in Gu’s model.[1]

Investigating the influence of particle cohesivity on the solids’ flow properties in a fluidized bed. Specifically, this is the mean solid volume fraction in the emulsion phase and the bed, solid velocity, and the solids’ tensile and compression pressure to quantify the solids’ ability to flow.

Analyzing the contribution of shear stress and yield stress in the flow properties of cohesive powders for a wide range of NBond and αcoh,1. To do so, we will evaluate the cumulative distribution of the ratio of shear and yield stress. This analysis enables us to explain the formation of different fluidization regimes for cohesive powders.

Mapping the granular flow regime for cohesive powders. To achieve this goal, a set of simulations is conducted on a wide range of particle cohesivity (characterized by the Bond number and tensile pressure prefactor) at various fluidization velocities. The obtained regime map will be connected to the cumulative distribution of the shear-to-tensile stress ratio.

Evaluating how the mixing quality depends on particle cohesiveness. To do so, the mean granular temperature and mean solids velocity variance will be computed for a studied range of cohesiveness levels. These two quantities can be associated with the level of mixing on the particle level (micromixing) and the global level (macromixing), respectively.

To achieve these goals, the open-source code MFiX (Multiphase Flow with Interphase eXchange)[20] simulation platform was extended to consider the solids rheology developed by Gu et al.[1] (see Section for the detailed description). In Section , the full description of the governing equations and the constitutive laws for the kinetic theory of granular flow will be presented. In Section , the modified solids rheology developed by Gu et al.[1] will be described. The results of our simulation will be presented in Section . Specifically, the performance of an FB of cohesive powders will be analyzed qualitatively and quantitatively in Sections and 3.2, respectively. In Section , the predicted regime map of fluidization will be presented and analyzed based on the ratio of shear stress and yield stress. Finally, we will evaluate the quality of mixing in FBs of cohesive powders in Section .

Mathematical Methodology

Hydrodynamics

As described in Section , the main goal of this study is to evaluate the capability of the TFM approach in predicting the fluidization behavior of cohesive powder. Therefore, the Eulerian–Eulerian approach was applied to simulate the gas–solid flow in the present study. In this approach, different phases are treated mathematically as interpenetrating continua. Conservation equations were derived for each phase and are linked by correlations for interphase transfer rates of momentum. To calculate solid-phase rheological properties, we employed the kinetic theory of granular flow (KTGF) in connection with a frictional stress model.[20] The complete lists of the governing equations and the constitutive laws have been summarized in Tables –3. The cohesiveness of the particles was considered through the modified KTGF model proposed by Gu et al.,[1] as described in Section .

Table 1

Momentum Equation Used in TFM Simulation

momentum conservation equations for gas and solids phases
T1.1
T1.2
interphase momentum transfer
T1.3
momentum exchange coefficient
T1.4
T1.5
T1.6
solid-phase stress tensor
T1.7
solids shear viscosity
T1.8
solids pressure
T1.9
gas-phase stress tensor
T1.10

Table 3

Constitutive Equations for the Calculation of the Solids Stress Tensor

solids viscosity
solids shear viscosity
T3.1
T3.2
T3.3
T3.4
solids bulk viscosity
T3.5
solids pressure
T3.6
T3.7
frictional stress
T3.8
solids frictional viscosity
T3.9
second invariant of the deviatoric stress tensor
T3.10
frictional solids pressure
T3.11

Constitutive Models for the Rheology of Cohesive Particles

The rheological model developed by Gu et al.[1] for cohesive particles (due to the presence of van der Waals forces) was implemented in the MFiX source code. To do so, the conservation equation for granular energy and solid viscosity and pressure were modified, as presented in Table . In eq , Fcohmax is the maximum cohesion force, which can be calculated based on the desired Bond number and the particle weight as Fcohmax = NBompg.

Table 4

KTGF-Based Models for Cohesive Particles Developed by Gu et al[1]

solids viscosity
T4.1
T4.2
T4.3
T4.4
T4.5
solids pressure
T4.6
T4.7
Pseudo-thermal energy dissipation (interparticle collisions)
T4.8
T4.9
T4.10
T4.11
with T4.12

As indicated in eq , for εs < εc, the first term (with the prefactor of αcoh,1) has a negative sign. This means that the force is acting in the direction reverse to the solids pressure. Therefore, one can expect two types of forces between the particles: the first one is the compression pressure which is related to (i) the kinetic pressure associated with the particles’ velocity fluctuations, similar to the pressure in the kinetic theory of gas, and (ii) the force chains formed among the particles that can transmit stress over larger distances.[1] The latter requires a large particle volume fraction, as indicated in eq , for εs ≥ εa with αcoh,2 being the relevant prefactor for this term. The second one, which is negative, is the attractive pressure which is attributed to the van der Waals force and acts as a tensile pressure. As the tensile pressure increases, one should expect the formation of larger agglomerates. This means that a higher value of the parameter αcoh,1 imposes a more significant cohesive force between the particles in the opposite direction in reference to the pressure.

The tensile pressure prefactor, αcoh,1, represents the attractive pair interaction between the particles. Therefore, it is mainly active in the low solids volume fractions (i.e., εs ≤ 0.2). Specifically, the tensile pressure prefactor relates the tensile pressure in a granular material to a typical cohesion force between individual particles. According to Gu et al.,[1] low solids volume fractions do not allow the formation of force chain as the stress cannot be transmitted to a large distance (low coordination number regime). Therefore, attractive van der Waals forces decrease the solids pressure. This force features the sign opposite to the classical compression pressure.

Determining the tensile pressure prefactor is a critical factor in accurately predicting solids pressure. However, this parameter cannot be directly measured, similar to the cohesion parameters in DEM simulations. Hence, a kind of calibration is required. In the original study of Gu et al.,[21] this parameter has been estimated through DEM-based simulations of homogeneous, simple shear flows of frictional and cohesive particles. They used the Levenberg–Marquardt method[22] to estimate the value of the tensile pressure prefactor, αcoh,1. The results of their simulations revealed that αcoh,1 is a function of the interparticle friction coefficient. Thus, one can relate the tensile pressure prefactor to a typical particle roughness and the tendency of a granular material to show a significant tensile strength.

It should be noted that in the original rheological model of Gu et al.,[1] the pressure at the volume fraction close to the packed condition of jamming cannot be predicted in a physically meaningful way as the term (εs – εa)2/εc – εs will diverge at εs = εc and predict an unphysical negative pressure at εs > εc. The original model of Gu et al.[1] does not account for the existence of highly packed regions (i.e., εs > εc). In other words, their study needs to be extended to develop a pressure equation under highly packed conditions. Therefore, using the current version of their model would result in unphysical negative pressures at volume fractions exceeding the jamming point. Therefore, in the present study, the (compressive and tensile) pressure was limited to zero for εs ≥ εc to avoid unphysical behavior and convergence problems. More details can be found in Appendix A in the Supporting Information.

In Table , W(x) is the Lambert-W function, which is defined as W–1(x) = x·exp(x). For numeric purposes, and as suggested by Winitzki,[23] for x > 0, eq can be used to approximate this function. As reported by Irani et al.[24] and Gu et al.,[1]W(x)/x represents the competition between the kinetic energy provided by the shear stress and the cohesive energy.

Results and Discussion

In the present work, the main focus is given to the capability of the TFM approach in predicting the fluidization behavior of cohesive powers. Therefore, an attempt was made to perform a qualitative and quantitative analysis of solids flow properties. In detail, we first evaluate the capability of TFM in predicting the formation of particle clustering and agglomerate formation qualitatively, as detailed in Section . The predicted behavior will be analyzed based on the voidage and solids pressure distribution. Subsequently, we will quantitatively examine the effect of particle cohesiveness on solids velocity, solids volume fraction in the emulsion phase, and the bed (see Section ). Afterward, we will propose a regime map of fluidization for cohesive powders by analyzing the distribution of solids pressure, solids shear stress, and yield stress (see Section ). Finally, in the last Section , the effect of particle cohesiveness on mixing quality will be investigated.

Setup and Parameter Ranges

To evaluate the capability of TFM in predicting such behaviors, we performed a set of simulations covering a wide range of particle cohesivity at different fluidization velocities. The simulations were performed in a pseudo-2D bed from the study of Li et al.[11] with a depth of 0.032 m with appropriate wall boundary conditions to model the effect of walls on the fluidization behavior of cohesive powders. A schematic illustration of bed configuration is depicted in Figure . Conducting a mesh dependency study revealed that the mesh resolution of 4dp is enough to accurately capture particle flow behavior regarding the solids volume fraction and velocity distribution. The Johnson and Jackson boundary condition[25] was used at the wall for pseudo-thermal granular energy.

Figure 1

Schematic representation of the studied fluidized bed, the black region shows the dense particle bed at initialization.

The minimum fluidization velocity for noncohesive powders has been defined by performing a set of simulations featuring various fluidization velocities. Subsequently, the resulting bed pressure drop and voidage distribution were analyzed to determine the minimum fluidization velocity. In detail, the pressure drop remains almost constant for u ≥ umf.

The detailed properties of the particle and fluidization gas and the operating conditions are reported in Table .

Table 5

Simulation Conditions and Parameters Used in the TFM Approach

parameter	base case	studied range
bed geometry
Hbed [m]	0.13
Lbed [m]	0.0508
Wbed [m]	0.0032
H0 [m]	0.06	0.06–0.12
particle properties
ρs [kg/m3]	1500
dp [μm]	148
ew,p [−]	0.955
epp [−]	0.9
NBo [−]	960	0–3840
αcoh,1 [−]	3 × 10–3	0–10
umf [m/s]	1.73 × 10–2
gas-phase properties
ρg [kg/m3]	1.188
μg [Pa·s]	1.79 × 10–5
u [m/s]	5umf	2umf–5umf
ut [m/s]	0.649
Rep (ρgutdp/μg)	7.03
Stp (1/18ρp/ρgRep)	451
wall boundary condition
gas	no-slip
solid	partial-slip
granular temperature	Johnson and Jackson

Qualitative Analysis of the Particle Cohesion Effect on the Solid Flow Behavior

As mentioned earlier, cohesive powders exhibit different flow behaviors based on the cohesiveness level. Therefore, in this section, we will examine how the particle Bond number can influence the solid particle interactions in the form of particle clustering and agglomeration. We will also evaluate how the particle Bond number can influence the powder flow near the wall due to cohesion.

Qualitative Behavior of the Fluidized Bed of Cohesive Powders

Agglomeration Formation

The main difference in fluidization of noncohesive and cohesive powders is the flowability of the powder and the formation of particle clusters in the bed. Therefore, we conducted a set of simulations for different Bond numbers to see if the agglomeration formation is predicted. The results of simulations revealed that in the studied range of the Bond number, with the model parameters obtained by Gu et al.,[1] no agglomeration was observed. As described in Section , the formation of agglomerates can be closely attributed to the presence of tensile (i.e., negative part) solids pressure. This means that the tensile pressure is too low to form agglomerates in the bed. The prefactor αcoh,1 partially contributes to the tensile pressure between the particles constituting the cohesive powder. Therefore, it seems that the default value of 3 × 10–3 represents the behavior of mildly cohesive powder only. Hence, a higher value is required to see agglomerates forming in the bed. Therefore, when we increase the value of the αcoh,1 parameter, particle clustering and agglomerate formation can be predicted, as shown in Figure . This can also be supported by the CFD–DEM results of Wu et al.[13] They compared the negative solids pressure by CFD–DEM with the one from Gu’s correlation. They reported that a higher value of αcoh,1 is required to match the tensile pressure. The distribution of negative pressure is presented in the right panel in Figure , for the sake of clarity. As shown in this figure, negative pressure can be observed on the bubble caps and partially inside the bubbles.

Figure 2

Distribution of negative solid pressure (right panel) and voidage (center panel) upon the formation of agglomerates in the bubbles for NBond = 3840 and αcoh,1 = 5 × 10–2 at u = 5umf and a qualitative comparison with the data of Li et al. (the most left panel).

As seen in Figure , the agglomerates/clusters are observed on the bubble cap. This predicted behavior is in accordance with the data reported by Li et al.[11] The high negative pressure is observed on the bubble cap where particle cluster formation can emerge. However, Figure suggests that the highest negative pressure (indicated by the black color) is predicted at the interface of bubbles with the emulsion phase and specifically for the particle traveling back to the bed surface after bubble bursting in the splash zone.

Figure 3

Flow path of the agglomerates (with negative pressure) in the bubble over time for NBo = 3840 and αcoh,1 = 5 × 10–2 at u = 5umf.

It should be highlighted that Gu’s model does not predict the pressure under packed conditions above the jamming point. The reader is referred to the discussion Section for more details. We conducted simulations with an alternative model formulation (i.e., unlimited pressure for εs ≥ εc), and the results of these simulations are presented in Appendix A of the Supporting Information.

The flow path of the agglomerates is shown in Figures –5 in more detail. As seen in Figure , small particle agglomerates and large clusters can be observed in the bubble and below the bubble cap, respectively (e.g., at t = 1.52 s). These clusters agglomerate and form a larger cluster while the bubbles rise in the FB. This breaks the bubble into several smaller bubbles (e.g., see three small bubbles in the right-top half of the bed at t = 1.6 s of Figure ). Subsequently, the clusters are easily disrupted and merged again into the emulsion phase. This leads to the positive pressure in the bubble wake, as seen in Figure . On the other hand, the small agglomerates are carried by the bubble, and they can break on the bed surface during bubble bursting, where a relatively high negative pressure is predicted for the solids phase (Figure ).

Figure 5

Distribution of negative pressure for the emulsion phase (with εs ≥ 0.3) (top panel) and the distribution of voidage in the agglomerate (with negative pressure) over the time for NBo = 960 and αcoh,1 = 5 × 10–2 (bottom panel) for the face close to the wall (z = 0.001 cm) at u = 5umf.

The predicted temporary agglomerates can be considered as dynamic aggregates, as suggested by Valverde.[10] According to Valverde,[10] agglomerates can be categorized as (i) dynamic aggregates and (ii) cohesive aggregates. Dynamic aggregates can easily break by the gas flow, and their size is limited by the balance between interparticle forces and shear forces. On the other hand, the gas flow cannot disrupt cohesive aggregates.

Nonetheless, these clusters feature relatively weak bonds and can easily break (Figure ); a clearer visualization of the granule path is provided in a video (Appendix B in the Supporting Information). The agglomerates with higher negative pressure cannot be carried by bubble rise motion and fall into the emulsion phase (see Figure , the panels for t = 1.56–1.58 s). In Figure , the motion of clusters on the bubbles’ caps (denoted by negative pressure in the bottom panel) has been presented for a less cohesive powder (NBo = 960). Upon comparing with the bed of particles with a higher Bond number of 3840 (Figure ), three points can be noticed:

Figure 4

Distribution of voidage in the agglomerate over time for NBo = 960 and αcoh,1 = 5 × 10–2 in a slice located at the center of the bed at u = 5umf.

The maximum tensile pressure is lower for less cohesive powders. This results in the formation of a much smaller number of agglomerates in the bubbles for NBo = 960. As seen in Figure , no agglomerate is predicted in the bubbles in the top half of the bed. Nonetheless, a small number of agglomerates can be observed close to the distributor surface.

The clusters formed in the bubble are much smaller than the ones in NBo = 3840, which can also be related to the smaller tensile pressure for the lower Bond number.

For both Bond numbers, a relatively high negative pressure is predicted above the bubble’s cap, where particle clusters are more likely to be observed. However, a larger region is influenced by the negative pressure for NBo = 3840 over the bubble’s cap. This behavior suggests that cluster formation is more likely for more cohesive powders, as one would also expect from common sense.

To evaluate the effect of the Bond number on fluidization of cohesive powder, we carefully analyzed different numerical and experimental studies from the literature. Several points can be understood:

The value of the Bond number highly depends on how the cohesive force is calculated. For instance, in the study of LaMarche et al.,[7] the van der Waals force is defined based on an asperity radius and the height related to the particle roughness parameter. In their study, for a particle size of 70 μm, a Bond number of 7 and a maximum contact force of 18 nN were reported. However, when using the method of Gu et al.,[1] the values of 228 and 1004 nN were obtained for this powder. The main reason for such a difference is considering the small-scale roughness parameters in the study of LaMarche et al.[7] Thus, smooth particles can lead to extremely high Bond numbers.

The fluidization behavior of cohesive particles cannot be purely represented by the Bond number. As demonstrated in the present contribution, the solids pressure also plays a key role in the formation of agglomerates and particle clusters. The solids pressure is governed by the Bond number and the tensile pressure prefactor. This can also be supported by comparing the simulation results of Gu et al.[1] and Wu et al.[13] In the study of Gu et al.,[1] the Bond number is in the range of 100–4000 with a tensile pressure prefactor of 3 × 10–3. On the other hand, in the study of Wu et al.,[13] the Bond number is smaller than 20, while αcoh,1 = 682–785. This means that reducing the Bond number requires a higher tensile pressure prefactor to exhibit similar solids pressure and fluidization behavior in terms of particle agglomeration and cluster formation. This finding has also been reported by Wu et al.[13]

The study of Gu et al.[9,21] shows the dependency of αcoh,1 on the particle friction coefficient. Therefore, having considered the role of αcoh,1 alongside the Bond number in fluidization behavior, the formation of agglomerates is indirectly governed by particle friction. This can also be supported by the study of LaMarche et al.,[7] who included particle friction (through roughness) in the definition of the Bond number (see the discussion in Point i).

Particle Flow Behavior near the Wall

It would also be of interest to see how the particle cluster formation and their flow are influenced by walls. To evaluate this behavior, the contour plots for the solids negative (tensile-dominated) pressure and the voidage on the wall have been depicted in Figure . Comparing Figures and 5 (shown in a contour plot at the center and the wall, respectively) reveals the formation of larger clusters close to the wall, followed by shrinkage of the bubbles. As seen in Figure , a broader region on the wall experiences negative pressures compared to the center of the bed. In addition, the bubbles appear much more elongated close to the wall because of large clusters.

In order to evaluate the effect of the wall on the fluidized bed performance, we performed an additional simulation for the case with NBond = 3840 and αcoh,1 = 5 × 10–2 with a static bed height of 12 cm (i.e., twice the size of the base case). The same fluidization regime as the one in the base case was predicted. Moreover, the results of the performed simulation demonstrated that the solid-phase flow behavior is not affected by the static bed height. Specifically, the solids (bed-averaged) volume fraction changed by 1.09%, the solids mean velocity changed by 4.4%, and the mean granular temperature changed by 6.3%. This finding is also in accordance with the experimental data of LaMarche et al.,[7] who reported that mean bed voidage is independent of the bed size.

Bond Number Effect on the Flow Behavior

In this section, we aim on quantifying the effect of particle cohesiveness on the solid flow behavior. In detail, we will evaluate the effect of the Bond number on the flow quantities of solids phase, including solids volume fraction, solids compression pressure, and tensile pressure under the steady-state and dynamic conditions as detailed below.

Time-Averaged Voidage Distribution

The contour plot for time-averaged voidage has been depicted in Figure . According to this figure, a more uniform distribution of voidage in the bed was predicted for lower Bond numbers. This will also be supported later by the cumulative distribution of the solids volume fraction in the FB, as reported in Section . This can be attributed to the motion of particles as a cluster in the bubble at high Bond numbers (e.g., NBo = 3840), where bubble coalescence occurs at lower heights relative to the distributor. On the other hand, at low Bond numbers, bubbles prefer to flow in the off-central positions, as indicated by the regions featuring a low voidage in Figure . The time-averaged voidage distribution in this figure indicates a higher bed expansion for higher Bond numbers. However, as discerned in Figure , the mean voidage in the emulsion phase reveals that a denser emulsion phase was predicted for more cohesive powders (higher NBond). This is valid for the cases in which the bed operates in the bubbling or bubbling/clustering regimes. Another point discerned from Figure is that a dense-solid region was predicted near the walls, close to the bed surface (freeboard) for the studied range of Bond numbers. This region is expanded toward the distributor when increasing the Bond number.

Figure 6

Effect of the Bond number of the time-averaged voidage distribution in the center slice of the fluidized bed at u = 5umf.

Figure 9

Comparison of the time-averaged ratio of tensile and compression pressure for different Bond numbers in the center slice of a fluidized bed at u = 5·umf.

Time-Averaged Solid Pressure Distribution

The distribution of the solids’ negative (i.e., tensile-dominant) pressure and positive (i.e., compression-dominant) pressure is shown in Figures and 8. At constant αcoh,1, a wider region featuring negative pressure is predicted for the higher Bond numbers. In detail, at a very low Bond number (e.g., 48), negative pressure is only observed in the region immediately above the distributor surface (Figure , the most left panel). However, at a higher Bond number of 480, negative pressure was also predicted near the bed surface in addition to the distributor. The negative pressure on the bed surface is attributed to the small agglomerates residing on the bed surface (splashing zone), as shown in Figure . These agglomerates are carried over by the bubbles to the bed surface, or they may be formed when large agglomerates break on the surface through bubble bursting. Therefore, one can expect a negative time-averaged pressure on the bed surface, which is dominated by tensile pressure. As seen in Figure , an increase in the Bond number expands these two regions in the axial direction inside the dense bed. This expansion continues at higher Bond numbers (e.g., 3840) to the extent that negative (tensile-dominated) pressure can be observed everywhere in the dense bed (first right panel of Figure for compression-dominated pressure). It should be noted that as legend bars in Figure indicate, the normalized negative pressure increases with the Bond number.

Figure 7

Effect of the Bond number on the time-averaged negative solids pressure (tensile-dominant pressure) in the center slice of the fluidized bed at u = 5umf.

Figure 8

Effect of the Bond number on the time-averaged positive solids pressure (compression-dominated pressure) in the center slice of the fluidized bed at u = 5umf.

In Figure , the distribution of time-averaged positive (compression-dominated) pressure was presented. As shown in this figure, positive pressure is observed in the top-half region of the bed. The highest compression-dominated pressure is observed immediately below the top negative-pressure region. This is associated with the high solids volume fraction in the top half of the bed (see Figure ), which results in powder compression. As we approach the distributor surface, the compression-dominated pressure decreases and the tensile pressure becomes dominant.

In order to obtain an insight into the contribution of tensile pressure and compression pressure to the total pressure, we computed the ratio of time-averaged tensile pressure and the compressive part based on eq . As shown in Figure , as the particle Bond number increases, a larger area of the bed experiences . As easily discerned from this figure, at NBond = 3840, the tensile pressure is dominant in dense bed regions, except for the one close to the distributor. This means that the probability of agglomerate formation is much higher for this Bond number. The predicted distribution is also in accordance with the distribution of tensile-dominated pressure and compression-dominated pressure shown in Figures and 8, respectively.

Quantitative Analysis of the Particle Cohesion Effect on the Solids Mean Flow Properties

It is also of interest to see how the particle cohesiveness influences the mean solids phase properties in the FB quantitatively. Therefore, the solids flow properties have been quantified through (i) solids volume fraction in the emulsion phase and the bed and (ii) the mean solids velocity in the bed. To achieve this goal, we investigated the effect of the Bond number and αcoh,1 parameter on mean solids properties (e.g., domain-averaged solid volume fraction in the bed and the emulsion phase and solids velocity).

We computed the mean time-averaged solids flow properties, such as solids volume fraction in the bed and the emulsion phase (Figure ), as well as the solids velocity (Figure ). It was demonstrated that unique trends cannot be extracted for the dependency of these quantities on the Bond number and tensile pressure prefactor. This is since the bed falls into different flow regimes at different values for Bond numbers and the αcoh,1 parameter (see Section for regime mapping).

Figure 10

Effect of the Bond number and tensile stress parameter (αcoh,1) on the time-averaged mean solids volume fraction in the emulsion phase (right panel) and in the bed (left panel) in the studied fluidized bed at u = 5umf.

Figure 15

Effect of the Bond number and tensile stress parameter (αcoh,1) on the time-averaged mean solids velocity in the fluidized bed at u = 5umf.

Solids Volume Fraction Versus Bond Number and αcoh,1

As shown in Figure , at a constant Bond number, an increase in αcoh,1 reduces the solids volume fraction, εs, in the emulsion phase. For αcoh,1 = 1 × 10–2, an increase in the Bond number gives rise to the solids volume fraction of the bed and the emulsion phase. Comparing the right and left panels in Figure revealed that the predicted solids volume fraction is almost similar for the bed and the emulsion phase for NBO = 960–3840 and αcoh,1 = 1 × 10–2. The identical voidage of the emulsion phase and the bed indicates that no bubbles are formed, and the bed falls into the bubble-less expansion regime. This can also be supported by the contour plot of the mean voidage shown in Figure . The bubble-less expansion of the bed at these two Bond numbers is due to the high yield stress that is larger than the applied shear stress. Figure supports this finding: the yield stress dominates the shear stress approximately by 1 order of magnitude. This means that the fluidization gas cannot shear the particles against the action of cohesion force. As it will be later explained in Figure , the solids velocity is zero for these two conditions, demonstrating that the fluidization gas cannot carry particles.

Figure 11

Distribution of solids yield stress and shear stress for the bed experiencing bubble-less expansion (NBo = 3840 and αcoh,1 = 1 × 10–2 at u = 5umf).

The analysis of bed voidage evolution revealed that the bed initially operates in the bubbling mode (data not shown here). As time passes, the emulsion phase close to the distributor surface rapidly expands (with uniform voidage) and does not allow the formation of bubbles anymore. The height of this expanded bed increases up to the bed surface after a few seconds such that no bubble is observed in the bed, and the whole bed is operated in the bubble-less expansion regime.

However, for higher tensile pressure prefactors (i.e., αcoh,1 = 5 × 10–2 to 10–1), a higher solids volume fraction was predicted in the emulsion phase. This reduction in the emulsion phase’ voidage is since the interparticle contacts are disrupted due to bubble passage and solids circulation, as reported by Geldart.[8] It should be noted that the solids volume fraction in the emulsion phase slightly increases with the Bond number for these two values of the αcoh,1 parameter.

As the Bond number increases, the effect of αcoh,1 on the emulsion solids volume fraction is intensified. For instance, a change in αcoh,1 marginally affects εs in the emulsion phase at NBo = 48. For αcoh,1 = 10–2, the increase in the Bond number results in higher bed expansion (lower solids volume fraction). An almost similar mean voidage was predicted for Bond numbers of 960 and 3840 due to bubble-less expansion as explained above.

For αcoh,1 = 5 × 10–2, a high difference between εs in the bed and emulsion phase demonstrated that an increase in the Bond number gives rise to the bubble volumes. This can also be supported by the voidage distribution depicted in Figure . This shows that more cohesive powders result in larger bubble formation in the fluidized bed. This behavior can be attributed to particle clustering and agglomeration. The attractive force between the particles is larger for a more cohesive powder. Therefore, the fluidization gas cannot easily shear the particles. This can also be supported by the plot of the shear rate distribution in Figure , which demonstrates that the yield stress in the emulsion phase is far larger than the shear stress.

Figure 12

Effect of the Bond number on the instantaneous voidage distribution inside the fluidized bed at t = 4 s (αcoh,1 = 5 × 10–2) at u = 5umf.

Figure 13

Distribution of the solids yield stress and the shear stress for the bed experiencing particle clustering (NBo = 3840 and αcoh,1 = 5 × 10–2 at u = 5umf) and a flow time of 1.6 s.

For a detailed investigation, the laterally time-averaged solids volume fraction was computed along the bed. To do so, the time-averaged εs was averaged over the cross-sectional area at each height. As shown in Figure , for NBo = 3840, the level of heterogeneity in εs is higher compared to that observed for lower Bond numbers. Specifically, for NBo ≤ 960, the voidage distribution is more uniform along the bed excluding the distributor surface and splashing zone. The reason is the formation of large bubbles in the bed filled with the most cohesive powders (NBo = 3840). Such a distribution results in a lower mean solids volume fraction (ε̅sbed) and consequently higher bed expansion, as supported by Figures and 12. Such a heterogeneity can also be understood by the cumulative distribution of voidage later presented in Figure .

Figure 14

Distribution of the time-averaged solids volume fraction (averaged over the cross-sectional area) along the bed for various Bond numbers and αcoh,1.

Figure 17

Cumulative distribution of the solid volume fraction for various Bond numbers and αcoh,1 at u = 5umf.

Solids Velocity Versus Bond Number and αcoh,1

The predicted flow behavior can also be explained by considering the solids mean velocity. As shown in Figure , the mean solids velocity is more influenced by the Bond number rather than αcoh,1 in the studied range of particle cohesiveness. The only exception is the cases with NBo ≥ 960 and αcoh,1 = 1 × 10–2. This is since in these cases, the bed falls into the bubble-less expansion regime, and the solids velocity drops to zero. On the other hand, the solids velocity decreases with an increase in the particle Bond number for NBO ≤ 960. This can be since the fluidization gas has less capacity to shear the particles due to the increase in the dominance of van der Waals forces with increasing Bond number. This finding highlights that the competition between the shear stress and yield stress plays a significant role in the flow behavior of the FB. As detailed in the following section, such a competition can also influence the flow regime of cohesive powders in the bed.

Predicted Fluidization Regimes

As mentioned earlier, the flow behavior of cohesive powders in a fluidized bed is mainly governed by the competition between the solids shear stress and yield stress. As the yield stress increases, higher shear stress is required to shear the particles and form bubbles in the fluidized bed. In this section, we analyze the ratio of solids shear stress and yield stress. The main objective is to quantify the fluidization regime for cohesive particles through a regime map. To do so, we will first compute the instantaneous voidage and Rτ = τshear/τyield for all cells in the simulation region. To do so, the shear stress is computed in each cell asin which the solids shear viscosity is given by eq and the shear rate, γ̇, is given byHere, is the magnitude of the rate of the deformation tensor, given by eq . The expression for the yield stress has been given in eq . Then, the cumulative distribution of Rτ was calculated for the studied range of the Bond number and αcoh,1. Afterward, the predicted regime map will be explained based on the obtained cumulative distribution of Rτ.

Contribution of Shear Stress and Yield Stress

The distribution of Rτ for a fluidization number, u/umf, of 5 has been plotted in Figure . As shown in this figure, for a constant value of αcoh,1 (e.g., αcoh,1 = 0.05), an increase in the Bond number leads to a rise in the fraction of cells with dominant yield stress (i.e., τshear/τyield < 1). This is since the yield stress directly increases with the Bond number, according to eq .

However, a different behavior was predicted for a constant Bond number over the studied range. In detail, for the constant Bond number of 48 and 480, an increase in αcoh,1 leads to a rise in the fraction of cells with dominant yield stress (τshear/τyield < 1). On the other hand, the cases with NBo = 3840 and αcoh,1 = 0.01 deviate from this trend. Therefore, such a trend for the spatial distribution of τshear/τyield suggests different flow regimes for various combinations of Bond numbers and αcoh,1.

Specifically, for the case with NBo = 3840 and αcoh,1 = 0.01, the shear stress is much smaller than the yield stress in all computational cells. This reveals the dominance of the yield stress over shear stress, meaning that the fluidization gas cannot shear the particles and form bubbles. Therefore, bubble-less expansion is expected for such cases.

To prove this claim, we analyzed the cumulative distribution of the solids volume fraction. As shown in Figure , for NBo = 3840 and αcoh,1 = 0.01, excluding the cell with εs ≤ 0.05 (as first bin), 80% of cells are occupied by 40–42% solids (εs = 0.4–0.42). This supports the uniform distribution of solids in the bed. In this manner, in all the cells, the predicted τshear is smaller than τyield, as seen in Figure . This situation surprisingly occurs for high Bond numbers with low αcoh,1. Such a behavior was also predicted at a lower fluidization velocity of u/umf = 2 and αcoh,1 = 0 for NBo ≥ 48, while at the lower Bond number of NBo = 4.8 and αcoh,1 = 0, the bubbling regime was predicted.

Figure 16

Cumulative distribution of the shear-to-yield stress ratio for various Bond numbers and αcoh,1 at u = 5umf.

Analysis of our simulation results revealed that the smaller the value of αcoh,1, the lower both the yield stress and shear stress are. However, the shear–yield stress ratio is smaller, demonstrating that the reduction in shear stress is more significant than the yield stress. It should be noted that the solids volume fraction plays a significant role in the yield stress, according to eqs and T4.5. Very low values of αcoh,1 for NBo ≥ 960 resulted in the bubble-less expansion of the bed and consequently lower bed solids volume fraction. This leads to a drop in the yield stress as the yield stress is highly influenced by εs.

Consequently, it sounds not straightforward to explain the predicted behavior by analyzing the shear-to-yield stress ratio due to the interrelation between the voidage, granular temperature, and stresses. It should be added that increasing αcoh,1 gives rise to the solids negative pressure (tensile-dominant pressure). Therefore, the particle agglomerates are more likely to be formed. Gas can bypass the agglomerates and shear the particles more easily as one can easily see in Appendix C of the Supporting Information. Our simulation results also revealed that an increase in αcoh,1 gives rise to the solids granular temperature for the cases with dominant shear stress.

Cohesive Powder Flow Regimes in BFB

The analysis of the cumulative distribution of Rτ = τshear/τyield reveals that Rτ can be considered as a good indicator for the fluidization regimes of cohesive powders. To put in more detail, the predicted cumulative distribution of Rτ can be categorized into three groups as

Fcellcumulative (Rτ < 1) ≅ 1 (e.g., NBo = 3840, αcoh,1 = 0.01): the bed experiences a bubble-less expansion at the studied fluidization velocity (e.g., 5umf). The solids yield stress dominates the shear stress. Therefore, the fluidization gas cannot fluidize the particles, resulting in the uniform distribution of solid particles in the bed. This is also supported by the cumulative distribution of the solids volume fraction shown in Figure .

Fcellcumulative (Rτ < 1) ≥ 0.8 (e.g., NBo = 3840, αcoh,1 = 0.05 or NBo = 48, αcoh,1 = 3): the bed experiences the formation and breakage of particle clusters. These clusters can be carried by the bubbles and break at the bed surface or flow downward in the bubble and break. Therefore, they are not mostly permanent agglomerates considering the pressure distribution over the bubble cap. As easily discerned from Figure , the broader distribution of the solids volume fraction is predicted for this regime, which shows the highest level of heterogeneity in the FB due to the formation of large bubbles.

Fcellcumulative (Rτ < 1) ≤ 0.8 (e.g., NBo = 48, αcoh,1 = 0.0.5–0.1): the bed is operated in the bubbling regime with the marginal effect of tensile pressure. Therefore, particle clusters are rarely formed. According to Figure , the distribution of the solid volume fraction has a higher level of heterogeneity than the bubble-less expansion but is lower than the clustering bed.

It should be added that a transition zone was identified for 0.5 ≤ Fcellcumulative (Rτ < 1) ≤ 0.8. In this transition zone, the fluidization regime is mainly bubbling. However, small negative pressure values can be observed near the bubbles’ caps. Nonetheless, the tensile pressure is small compared to the stress caused by the gravitational force. This means that the formation of particle clusters and agglomerates rarely happens in this range. In other words, the higher maximum tensile pressure is expected for the cohesive bed with the higher resultant value of Fcellcumulative (Rτ < 1). To support this, the distribution of solid voidage and tensile pressure was compared, as shown in Figure , for the cases with Fcellcumulative (Rτ < 1) = 0.5 (NBond = 48, αcoh,1 = 0.05) and 0.8 (NBond = 48, αcoh,1 = 0.3). Therefore, the transition zone can be merged to Fcellcumulative (Rτ < 1) ≤ 0.8.

Figure 18

Distribution of voidage and tensile pressure for NBond = 48 and αcoh,1 = 0.05 (left panel) and αcoh,1 = 0.3 (right panel) at time = 4 s.

Therefore, four fluidization regimes for cohesive powder can be identified, as shown in Figure :

Figure 19

Regime map of fluidization for cohesive powder at u/um = 5.

Shear-dominant (bubbling); at low Bond numbers and low αcoh,1, particles are easily fluidized with the marginal effect of cohesion. In this regime, the shear stress significantly dominates the yield stress.

Yield-dominant (bubble-less expansion): at high Bond numbers and low values of αcoh,1, the dominance of the van der Waals force does not allow the fluidization gas shearing the particles, resulting in bubble-less expansion of the bed.

Intermediate (clustering–bubbling): at intermediate-to-high values of the Bond number and αcoh,1, the contribution of yield stress and shear stress is comparable. Therefore, particles can form temporary agglomerates which can easily break up with the bubble flow.

Purely yield-dominant (stationary): at high values of αcoh,1, the dominance of tensile pressure results in forming a big agglomerate in the bed. Therefore, the fluidization gas cannot fluidize the particles. It should be added that due to high tensile pressure, the solution can easily diverge.

It should be noted that one of the main purposes of the present work was to evaluate the fluidization regimes at different Bond numbers. Therefore, we have chosen a wide range of αcoh,1 values (to reflect a wide range of tensile pressures) to observe different fluidization regimes (e.g., bubbling, clustering–bubbling, bubble-less expansion, and stationary) at a fixed Bond number. However, for high Bond numbers (NBond > 960), even at αcoh,1 = 0, the bubbling regime was not predicted and the bed started to uniformly expand, as one can easily see in the predicted fluidization regime map in Figure .

As the fluidization velocity influences the flow regime in the fluidized bed,[26] we evaluated the fluidization of cohesive particles at a lower fluidization velocity, u = 2umf. As easily discerned in Figure , compared to the higher fluidization velocity of 5umf, the extent of each regime can vary in the design space of NBond – αcoh,1. Specifically, for u = 2umf the bubble-less expansion regime starts at lower Bond numbers (NBo > 48) and lower αcoh,1 (value of zero). This regime continues up to αcoh,1 = 5 × 10–2 at u = 2umf while up to αcoh,1 = 1 × 10–2 at u = 5umf. The main reason for such a behavior is that at lower fluidization velocity, the shear supported by gas cannot overcome the yield stress between the particles at constant NBo and αcoh,1. Therefore, the bed cannot get fluidized and falls into the bubble-less expansion. Hence, the bubbling regime is expected to occur at a lower Bond number. To evaluate this, a set of the simulation was performed for a lower Bond number of 4.8; the bubbling regime was predicted for αcoh,1 ≤ 5 × 10–2, which was not shown in the regime map in Figure . Therefore, one can conclude that in addition to particle cohesion properties, process parameters of fluidization gas play a significant role in the predicted regime map for cohesive powders.

Figure 20

Regime map of fluidization for cohesive powder at u/umf = 2.

Mixing Quality

In this section, the effect of particle cohesiveness was evaluated on the mixing quality at the particle level (micromixing) and global level (macromixing). Micromixing is associated with mixing due to the fluctuation of solid particles’ velocity, quantified by granular temperature.[27] On the other hand, macromixing is attributed to solids mixing due to the motion of bubbles in the fluidized bed, quantified by the solids velocity variance.[27] More information about the mixing can be found in our previous works.[27−29]

Micromixing (Particle Level)

To evaluate the level of micromixing, mean time-averaged solids granular temperature was computed for the studied range of the Bond number and αcoh,1. To do so, the mean granular temperature was weighted by the solids volume fraction in each cell aswhere NΔ is the total number of time steps and NCell stands for the total number of cells in the computational domain. The computed mean granular temperature is depicted in Figure . Generally speaking, an increase in the Bond number and αcoh,1 leads to a rise in the mean granular temperature. The only exception is when the bed falls into the bubble-less expansion regime. This finding is in accordance with our previous study[30] and the study of Luding[31] and Radl et al.[32] This revealed that particle cohesiveness increases the velocity fluctuation at the particle level. However, different behaviors were predicted for macromixing in the bubble level (see Figure ). According to Figure , for the higher Bond numbers, the effect of αcoh,1 on the micromixing becomes more significant.

Figure 21

Effect of the Bond number and tensile stress parameter (αcoh,1) on the time-averaged mean solid granular temperature in the fluidized bed at u = 5umf.

Figure 22

Effect of the Bond number on the solids axial velocity variance (left subpanel for each Bond number) and solids lateral velocity variance (right subpanel for each Bond number) in the fluidized bed for αcoh,1 at u = 5umf.

One can also analyze the predicted granular temperatures in connection with the predicted fluidization regimes. As discerned from Figure , the effect of the Bond number and αcoh,1 on the particle fluctuation intensity is marginal in the bubbling regime. Such a behavior is expected as the flow of particles is mainly dominated by the shear imposed by fluid flow. However, an increase in the Bond number or tensile pressure prefactor leads to a rise in the granular temperature in the clustering–bubbling regime. Luding[31] and Radl et al.[32] also reported this trend. A similar behavior can be observed for the bubble-less expansion regime.

Macromixing (Global Level)

In a fluidized bed, solids macromixing is mainly influenced by the motion of bubbles. Therefore, the variance of the hydrodynamic velocity is of high relevance to quantify macromixing.[33] In this regard, the solids velocity variance tensor was computed using the following equation in each computational cell[27]where and (t) represent the time-averaged and instantaneous solids velocity in the computational cell, respectively. ⊗ denotes the dyadic product. Solids mean velocity, , was weighted by the solids volume fraction and is given by

In Figure , the distribution of velocity variance was depicted for various Bond numbers and a constant tensile pressure prefactor of αcoh,1 = 5 × 10–2. This figure suggested that the lowest velocity variance can be predicted close to the wall. This is because the bubbles mainly flow toward the central bed region after forming on the distributor surface. Therefore, the bubble flow is experienced in a lower extent by the particles moving along the wall (bubble-off position). The time-averaged voidage distribution can also support this. As seen in Figure , in the regions close to the walls, the time-averaged voidage corresponds to the jamming condition, which locally reduces particle movability. Consequently, one should expect lower velocity variance in this region. Another point discerned from Figure is that the increase in the Bond number broadens the region with the low-velocity variance. This is associated with the voidage distribution plotted in Figure . As the Bond number increases, this “jamming” region expands toward the center and the distributor.

On the other hand, for αcoh,1 = 1 × 10–2, a different trend was predicted. Specifically, an increase in the Bond number reduces the velocity variance in both axial and radial directions in such a way that for the cases bubble-less fluidization is experienced (e.g., NBo ≥ 960), the variances fall to zero. As discerned from Figure , the axial velocity variance decreases with the Bond number for NBo ≤ 960. However, higher velocity variance was predicted for NBo = 3840. Such a behavior can be supported by the instantaneous voidage distribution in the bed. As shown in Figure , relatively large bubbles are formed in the bed for this Bond number. Such large bubbles result in higher fluctuation in the particle axial velocity over time. Consequently, higher axial velocity variance can be expected. For a higher value of αcoh,1, this behavior is observed at a lower Bond number of 960, as depicted in Figure . This means that the axial velocity variance hits a minimum at NBo = 480.

Figure 23

Effect of the Bond number and tensile stress parameter (αcoh,1) on the time-averaged mean solids axial (left panel) and radial (right panel) velocity variances in the fluidized bed at u = 5umf.

According to Figure , high Bond numbers (NBond ≥ 960) and high tensile pressure prefactors (αcoh,1 ≥ 5 × 10–2) lead to a rise in the predicted mean velocity variances. This is to some degree counterintuitive as one would expect that these velocity fluctuations, and consequently the global mixing rate, reduce upon an increase in the particle level of cohesion. However, the mixing dynamics of wet granular matter are highly complex, as shown in many studies, e.g., Radl et al.[11] The reason is that comparably large agglomerates when sheared at a given shear rate lead to a larger reference velocity simply because of dimensionality reasons. Also, kinetic theory of granular flow would predict increased granular temperature (and hence mixing) in the case of larger particles (or agglomerates) when considering a fixed shear rate.

The obtained trend of increased velocity variance at high Bond numbers can be explained to some degree by the contour plot of solids velocity variance in Figure and time-averaged voidage distribution in Figure . As easily discerned in Figure , the axial solids velocity variance can be divided into two relatively uniform regions featuring (i) high-variance region in the central region of the bed and (ii) low-velocity variance developing from the wall to the central region of the fluidized bed. Considering the voidage contour plot in Figure , one can conclude that the high velocity variance is predicted in the region where voidage is relatively high due to the motion of bubbles. This can also be supported by the instantaneous voidage distribution in Figures and 3. The reader is referred to Appendix C of the Supporting Information to see the motion of bubbles over time for a very cohesive bed. As we were unable to identify individual agglomerates in our present study, the exact mechanism that leads to increased velocity fluctuations at high Bond numbers could not be investigated in our present study.

Conclusions

The present study focused on the capability of the TFM approach in predicting the fluidization behavior of cohesive particles. To do so, the modified solids rheology for cohesive powders developed by Gu et al.[1] was implemented in the MFiX platform. The main advantage of the implemented model of Gu et al.[1] is that it is based on particle-level simulation of cohesive powder flow. This means that we do not need any empirical constitutive equation (e.g., agglomeration kernel in the population balance model) to predict particle agglomeration if the grid resolution is sufficiently fine. Specifically, this is since particle cohesion and agglomeration are already considered via tensile stress (i.e., a negative “tensile pressure”) in the solids rheology model of Gu.[1]

Our simulation results demonstrated that the TFM approach successfully predicted cohesive powder fluidization for the first time. The model was qualitatively validated against the data reported by Li et al.[11] in terms of (i) the formation of particle clusters in the bubble and (ii) the formation of a solids tensile stress at the bubble cap. Most importantly, the flow of powders as clusters and agglomerates was predicted. Moreover, the distribution of negative pressure was presented, which can be attributed to the formation of tensile stress due to particle cohesion. Furthermore, our simulation results revealed that the particle Bond number and the tensile pressure prefactor critically affect the level of particle clustering in the fluidized bed. To thoroughly evaluate this, we performed a set of simulations for a wide range of particle cohesion levels to obtain a regime map of cohesive fluidization. Analyzing the voidage and pressure contour plot revealed that four different particle flow regimes are developed in the Bond number–tensile pressure prefactor space. The predicted flow regimes were successfully attributed to the competition between the acting shear stress and the yield stress of the powder that tends to form larger particle assemblies (i.e., agglomerates). Considering the cumulative distribution of τshear/τyield, the four observed regimes can be categorized as follows:

Bubbling regime: the shear stress dominates the yield stress, resulting in no (or no significant) particle agglomeration.

Bubbling–clustering regime: the shear stress and the yield stress are of comparable magnitude such that the formation of particle clusters and agglomerates is expected mainly in bubbles.

Bubble-less expansion: the yield stress dominates the applied shear stress such that the fluidization gas cannot shear the particles. Hence, bubble formation is suppressed, and the bed becomes uniformly expanded.

Stationary bed: the yield stress is much larger than the applied shear stress such that the bed cannot get fluidized, and large agglomerates can be formed in the bed. Gas bypasses the particles through channeling. The computational cost for simulating such highly cohesive powders (e.g., high values of NBo or high values of αcoh,1) is significantly larger than that for powders in the other regimes. Simulation can even become unstable due to the very small time step that is necessary to handle the extremely high negative pressures.

The solids rheology model developed by Gu et al.[1] shows significant progress in the application of the TFM for cohesive fluidized beds. However, their model needs to be revisited for volume fractions beyond the jamming point. Their original model suggests high negative pressure at εs ≥ εc, which is unphysical. In our present study, the pressure due to cohesion was set to zero for densely packed situations, that is, above the jamming point.

Apart from this, the analysis conducted in our present work should be extended in future studies. Most importantly, the regime map of fluidization for cohesive powders was obtained at two different fluidization velocities at constant particle and gas properties. However, this regime map should be extended to consider the effect of the dimensionless diameter (d*) and the dimensionless velocity (u*) similar to Grace’s fluidization chart.[26] Furthermore, quantitative validation of the model requires (i) characterization of a set of powders with different cohesion levels to obtain accurate values of model parameters (e.g., the tensile and compression pressure prefactors) and subsequently (ii) a set of experiments using these powders in a fluidized bed.

Finally, predicting the agglomerate size would require a deeper analysis, which necessitates a detailed analysis of solids tensile pressure, as well as shear stress and yield stress. Specifically, the gradients of these quantities need to be tracked over time to evaluate if the particles have enduring contacts. The boundary between a specific agglomerate and the fluid could then be demarcated in a next step. One agglomerate can only be larger than one computational cell. Therefore, to compute the agglomerate size, one needs to perform an image analysis of solids volume fraction, solids tensile pressure, and yield stress to extract the agglomerate size. As the main focus of our present work is demonstrating the capability of the TFM approach in predicting agglomeration phenomena of a cohesive powder, the analysis of the agglomerate size distribution was postponed to a future study.

Table 2

Transport Equation for Granular Kinetic Energy

T2.1

Pseudo-thermal conductivity of solids phase

T2.2

T2.3

T2.4

Pseudo-thermal energy generation due to the interaction of solids phase with gas

T2.5

Pseudo-thermal energy dissipation due to inelastic particle–particle collision

T2.6