# A Bubble-Based Drag Model at the Local-Grid Level for Eulerian Simulation of Bubbling Fluidized Beds

Authors: Kun Hong; Manqian Cao; Atta Ullah
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

A bubble-based drag model at the local-grid level is proposed to simulate gas-solid flows in bubbling fluidized beds of Geldart A particles. In this model, five balance equations are derived from the mass and the momentum conservation. This set of equations along with necessary correlations for bubble diameter and voidage of emulsion phase is solved to obtain seven local structural parameters (uge, upe, εe, δb, ub, db, and ab) which describe heterogeneous flows of bubbling fluidized beds. The modified drag coefficient obtained from the above-mentioned structural parameters is then incorporated into the two-fluid model to simulate the hydrodynamics of Geldart A particles in a lab-scale bubbling fluidized bed. The comparison between experimental and simulation results for the axial and radial solids concentration profiles is promising.


---

## Body

## 1. Introduction

Gas-solid fluidized bed provides thorough mixing, high heat, and mass transfers between gas and solid phases. Due to these reasons, bubbling beds are widely used in industrial applications involving coal gasification, fluid-catalytic-cracking (FCC), and mineral processing [1, 2]. Over the last few decades, in particular, bubbling fluidized bed has attracted many research studies because of its unique characteristics. In earlier works, for example, the size and the rise velocity of bubbles were studied to explore and model the bubbling phenomenon [3–7]. Recently, computational fluid dynamics (CFD) modeling has been used intensively to investigate, understand, and optimize fluidized beds [8–14].In the CFD studies, Eulerian approach or the so-called two-fluid model (TFM) has been widely used to simulate industrial-scale reactors due to the limited computational resources [15–17]. In the TFM, the gas and the solid flows are statistically averaged and treated as interpenetrating continua. However, several researchers have concluded that the coarse-grid TFM with homogeneous drag failed to predict the hydrodynamics of bubbling fluidized beds with fine particles [18–20]. One of the possible reasons is that the standard TFM adopts classical empirical formula such as Ergun equation [21] and Wen-Yu [22] correlation to close the drag force acting on the particles. These empirical correlations are derived from the experiments performed on fixed beds, homogeneous fluidized beds or sedimenting suspensions, resulting in overlooking the structural (or so-called meso-scale) effect on the interphase interaction. For simulation of bubbling fluidized beds, some researchers have argued that the standard TFM may well capture the flow behavior if the grid is fine enough to the order of 10 times particle diameter [23–25]. In recent years, therefore, subgrid modeling has emerged as a reliable tool for accurate simulations of fluidized bed [26]. Among these subgrid models, a popular method of filtered TFM has been developed to simulate gas-solid flows at the coarse-grid level [27–29]. In this approach, constitutive models of drag force and solid stress are derived through filtered data obtained from highly fine-grid simulations [28–32]. But there is still no clear agreement on the feasibility of using the fine-grid TFM to simulate heterogeneous flows with fine particles [23, 33–37].Another practical way, to account for the effects of these meso-scale structures, is to employ a modified drag model which is established on the basis of heterogeneous structures [38]. For instance, some authors provide a scaling factor to revise the homogeneous drag [39–41]. Although it yields better results, however, such a strategy is also empirical and unsystematic. Alternatively, through the energy-minimization multiscale (EMMS) analysis [42], Li and his coworkers proposed their bubble-based EMMS drag where the bubble phase was treated as the meso-scale structure [20, 43–45]. Other researchers have directly established their own drag models on the basis of the bubble-emulsion structure [46, 47]. It is to be noted that, in the above models, the global operating conditions are used to deduce the correlation between drag coefficient and voidage (or bed height), which neglect the effect of local information on the drag coefficient. Subsequently, the bubble-based EMMS model has been revised and extended to the local-grid level [48–50]. In this development, the local energy dissipation was introduced and assumed to be minimal. However, there is not sufficient evidence to suggest that this kind of local energy dissipation is bound to be minimal for gas-solid fluidization [51–53].In present work, in light of the above mentioned limitations, a novel bubble-based drag model based on local-grid information is developed to deal with the impact of heterogeneity on the drag force. Firstly, within each computational cell, heterogeneous structures are resolved into the discrete bubble phase, the continuous emulsion phase, and their interphase. Then the structural parameters within each cell are computed by solving seven equations, involving mass conservation, force balance, and reliable correlations for the bubble diameter and emulsion voidage. Next the bubble-based drag coefficient is derived from the resolved structural parameters. Finally, this new drag coefficient is incorporated into the two-fluid model to simulate the hydrodynamics of Geldart A particles in a bubbling fluidized bed. Therefore, the current work can be considered as a tuning process where subgrid corrections are incorporated based on local grid information if the system was inhomogeneous. Comparison of simulation results to available experimental data is also presented for validation of this new drag model.


## 2. Mathematical Formulation

### 2.1. Resolution of Local Heterogeneity

For a bubbling fluidized bed, the gas-solid flow can be resolved into the bubble and the emulsion phases. The discrete bubbles are surrounded by the emulsion phase. Consequently, the local flow within each computational cell can also be divided into three subsystems: the emulsion phase, the bubble phase, and the interphase, as shown in Figure1. For simplification, it is assumed that the particles are homogeneously distributed in the emulsion phase, while no particles exist inside the bubble phase (εb=1). Then seven independent parameters are required to describe such local heterogeneity. This information includes the gas and the solid velocity in the emulsion phase uge/upe, the voidage of emulsion phase εe, the volume fraction of bubbles δb, the rise velocity of bubbles ub, the bubble diameter (size) db, and the bubble acceleration ab. As a first approximation, the emulsion phase has been treated as a pseudo-fluid with averaged density ρe, viscosity μe [54], and velocity ue, which are defined as follows: (1)ρe=ρp1-εe+ρgεe,μe=μg1+2.51-εe+10.051-εe2+0.00273e16.61-εe,ue=ρgεeuge+ρp1-εeupeρgεe+ρp1-εe.Figure 1
Resolution of local heterogeneity within each computational cell for bubbling fluidization.


### 2.2. Bubble-Based Drag Model at the Subgrid Level

For solving seven structure parameters, seven independent equations, involving mass conservation, force balance and empirical correlations, are employed as follows.Voidage inside the Emulsion. The voidage in the emulsion phase plays an important role to describe the emulsion state. The correlation of Harris et al. [55], which was derived from a lot of experiments using different solids and operating conditions, is applicable over a wide range of operating conditions. In this work, thus, the correlation of Harris et al. [55] is chosen to predict the voidage inside emulsion as follows: (2)εe=1-0.581-εg1.480.13+1-εg1.48.Mean Voidage within Each Computational Cell. According to definition of Favre averaging [56], the mean voidage is given as(3)εg=1-δbεe+δb.Mass Conservation for the Particles. The flow of particles within the grid should be equal to the particle flow rate in the emulsion phase. Then, the averaged particle velocity becomes(4)us=1-δb1-εeupe1-εg.Force Balance for Particles inside the Emulsion Phase. For the bubbling fluidization, there is no carryover of particles out of the beds. As a first approximation, the acceleration of the emulsion phase is negligible (ae=0) because of the large internal difference between the gas and the solids, as indicated in previous works [20]. Thus, the drag force can be assumed to be equal to the effective weight of the particles. If the Ergun equation [21], which is suitable for the dense emulsion, is used to calculate the drag force, then, the net force acting of the bed of particles is written as(5)1501-εe2μgεedp2+741-εeρguge-upedpuge-upeεe=1-εeρp-ρgg.Mass Conservation for the Gas. The mass flow of gas within each grid should be equal to the sum of that in the emulsion and the bubble phases. Accordingly, the averaged gas velocity within the grid is defined by(6)ug=1-δbεeuge+δbεbubεg.Bubble Diameter. Compared with ambiguous clusters or strands of particles observed in high velocity fluidization, it is relatively easy to describe bubbles. Several correlations to estimate the bubble size are available and are summarized in a recent review by Karimipour and Pugsley [57]. In our present work, the correlation of Horio and Nonaka [58] is selected owning to its broad applicability. The diameter of a bubble can be estimated as (7)db=-γM+γM2+4dbmDt2·Dt4,where γM and dbm can be expressed as7aγM=2.56×10-2Dt/gUmf,7bdbm=2.59g-0.2Ug-UmfπDt240.4.Force Balance for the Bubbles. The drag force of the emulsion phase exerted on the bubbles is equal to the effective buoyancy of bubbles in each computational cell: (8)34Cd,bδbdbρeUslip,be2=δbρe-ρgg+ab,where superficial slip velocity between the bubble and the emulsion phases Uslip,be and drag coefficient of multibubble Cd,b are, respectively, defined as8aUslip,be=1-δbub-ue,8bCd,b=1-δb-0.5Cd,b0;here8cCd,b0=2.7+24Rei,Rei>1.8,38Rei-1.5,0<Rei≤1.8,withRei=ρedbUslip,beμe.For specified local information (εg, ug and us), first this set of nonlinear equations ((2)–(8)) is solved to obtain seven structural parameters (uge, upe, εe, δb, ub, db and ab). With these structural parameters determined, then the bubble-based drag coefficient can be obtained by(9)βbubble=εgFDug-us=εgδbρe-ρgg+ab+1-δb1-εeρp-ρggug-us.


### 2.3. Solution Scheme

For each computational cell, the local information (εg, ug, and us) is directly extracted from the CFD solver. Once the local information has been given, as shown in Figure 2, this novel bubble-based drag model will be solved using the following scheme:(1)
Calculateεe from (2).(2)
Calculateδb from (3).(3)
Calculateupe from (4).(4)
Calculateuge from (5).(5)
Calculateub from (6).(6)
Calculatedb from (7).(7)
Calculateab from (8).(8)
Calculateβbubble from (9).Figure 2
Computation scheme of solving the bubble-based drag model.


## 2.1. Resolution of Local Heterogeneity

For a bubbling fluidized bed, the gas-solid flow can be resolved into the bubble and the emulsion phases. The discrete bubbles are surrounded by the emulsion phase. Consequently, the local flow within each computational cell can also be divided into three subsystems: the emulsion phase, the bubble phase, and the interphase, as shown in Figure1. For simplification, it is assumed that the particles are homogeneously distributed in the emulsion phase, while no particles exist inside the bubble phase (εb=1). Then seven independent parameters are required to describe such local heterogeneity. This information includes the gas and the solid velocity in the emulsion phase uge/upe, the voidage of emulsion phase εe, the volume fraction of bubbles δb, the rise velocity of bubbles ub, the bubble diameter (size) db, and the bubble acceleration ab. As a first approximation, the emulsion phase has been treated as a pseudo-fluid with averaged density ρe, viscosity μe [54], and velocity ue, which are defined as follows: (1)ρe=ρp1-εe+ρgεe,μe=μg1+2.51-εe+10.051-εe2+0.00273e16.61-εe,ue=ρgεeuge+ρp1-εeupeρgεe+ρp1-εe.Figure 1
Resolution of local heterogeneity within each computational cell for bubbling fluidization.


## 2.2. Bubble-Based Drag Model at the Subgrid Level

For solving seven structure parameters, seven independent equations, involving mass conservation, force balance and empirical correlations, are employed as follows.Voidage inside the Emulsion. The voidage in the emulsion phase plays an important role to describe the emulsion state. The correlation of Harris et al. [55], which was derived from a lot of experiments using different solids and operating conditions, is applicable over a wide range of operating conditions. In this work, thus, the correlation of Harris et al. [55] is chosen to predict the voidage inside emulsion as follows: (2)εe=1-0.581-εg1.480.13+1-εg1.48.Mean Voidage within Each Computational Cell. According to definition of Favre averaging [56], the mean voidage is given as(3)εg=1-δbεe+δb.Mass Conservation for the Particles. The flow of particles within the grid should be equal to the particle flow rate in the emulsion phase. Then, the averaged particle velocity becomes(4)us=1-δb1-εeupe1-εg.Force Balance for Particles inside the Emulsion Phase. For the bubbling fluidization, there is no carryover of particles out of the beds. As a first approximation, the acceleration of the emulsion phase is negligible (ae=0) because of the large internal difference between the gas and the solids, as indicated in previous works [20]. Thus, the drag force can be assumed to be equal to the effective weight of the particles. If the Ergun equation [21], which is suitable for the dense emulsion, is used to calculate the drag force, then, the net force acting of the bed of particles is written as(5)1501-εe2μgεedp2+741-εeρguge-upedpuge-upeεe=1-εeρp-ρgg.Mass Conservation for the Gas. The mass flow of gas within each grid should be equal to the sum of that in the emulsion and the bubble phases. Accordingly, the averaged gas velocity within the grid is defined by(6)ug=1-δbεeuge+δbεbubεg.Bubble Diameter. Compared with ambiguous clusters or strands of particles observed in high velocity fluidization, it is relatively easy to describe bubbles. Several correlations to estimate the bubble size are available and are summarized in a recent review by Karimipour and Pugsley [57]. In our present work, the correlation of Horio and Nonaka [58] is selected owning to its broad applicability. The diameter of a bubble can be estimated as (7)db=-γM+γM2+4dbmDt2·Dt4,where γM and dbm can be expressed as7aγM=2.56×10-2Dt/gUmf,7bdbm=2.59g-0.2Ug-UmfπDt240.4.Force Balance for the Bubbles. The drag force of the emulsion phase exerted on the bubbles is equal to the effective buoyancy of bubbles in each computational cell: (8)34Cd,bδbdbρeUslip,be2=δbρe-ρgg+ab,where superficial slip velocity between the bubble and the emulsion phases Uslip,be and drag coefficient of multibubble Cd,b are, respectively, defined as8aUslip,be=1-δbub-ue,8bCd,b=1-δb-0.5Cd,b0;here8cCd,b0=2.7+24Rei,Rei>1.8,38Rei-1.5,0<Rei≤1.8,withRei=ρedbUslip,beμe.For specified local information (εg, ug and us), first this set of nonlinear equations ((2)–(8)) is solved to obtain seven structural parameters (uge, upe, εe, δb, ub, db and ab). With these structural parameters determined, then the bubble-based drag coefficient can be obtained by(9)βbubble=εgFDug-us=εgδbρe-ρgg+ab+1-δb1-εeρp-ρggug-us.


## 2.3. Solution Scheme

For each computational cell, the local information (εg, ug, and us) is directly extracted from the CFD solver. Once the local information has been given, as shown in Figure 2, this novel bubble-based drag model will be solved using the following scheme:(1)
Calculateεe from (2).(2)
Calculateδb from (3).(3)
Calculateupe from (4).(4)
Calculateuge from (5).(5)
Calculateub from (6).(6)
Calculatedb from (7).(7)
Calculateab from (8).(8)
Calculateβbubble from (9).Figure 2
Computation scheme of solving the bubble-based drag model.


## 3. Simulation of Gas-Solid Bubbling Fluidization with Geldart A Particles

### 3.1. Governing Equations

In this work, the two-fluid model (TFM) is adopted to simulate bubbling fluidization. The CFD software, FLUENT® 6.3.26, was used as fluid solver to carry out all simulations. The governing equations of the two-fluid model for gas-solid flows are summarized as follows.Mass Conservation Gas phase:T1-1∂∂tεgρg+∇·εgρgug=0. Solid phase:T1-2∂∂tεsρs+∇·εsρsus=0.Momentum Conservation Gas phase:T1-3∂∂tεgρgug+∇·εgρgugug=-εg∇p+∇·τg+εgρgg-βug-us. Solid phase:T1-4∂∂tεsρsus+∇·εsρsusus=-εs∇p-∇ps+∇·τs+εsρsg+βug-us.Granular Temperature T1-5 3 2 ∂ ∂ t ε s ρ s Θ s + ∇ · ε s ρ s Θ s u s = - p s I + τ s : ∇ u s + ∇ · κ s ∇ Θ s - γ s - 3 β Θ s .Constitutive relations are as follows:Stress-strain tensor:
Gas phase:(10)τg=εgμg∇ug+∇ugT-23εgμg∇·ugΙ.Solid phase:(11)τs=εsμs∇us+∇usT+εsλs-23μs∇·usI.Solid phase pressure:(12)ps=εsρsΘs1+21+esεsg0.Solid phase viscosity:(13)μs=μs,col+μs,kin+μs,fr,μs,col=45εsρsdpg01+esΘsπ,μs,kin=10ρsdpΘsπ96εs1+esg01+45g0εs1+es2,μs,fr=pssin⁡Φ2I2D.Solid phase bulk viscosity:(14)λs=43εsρsdpg01+esΘsπ.Radial distribution function:(15)g0=1-εsεs,max1/3-1.Diffusion coefficient of granular energy:(16)κs=150ρsdpΘsπ3841+esg01+65εsg01+es2+2ρsεs2dp1+esg0Θsπ.Collision energy dissipation:(17)γs=121-es2εs2ρsg0Θs3/2dpπ.In the TFM, the continuity equations for the gas and the particle phases are given as T1-1 and T1-2. The momentum conservation for each phase is expressed by T1-3 and T1-4, where τ and β denote the stress tensor and the drag coefficient, respectively. The bubble-based drag coefficient formulated in this article can be adopted to calculate the effective drag coefficient which is incorporated into TFM through a user-defined function (UDF). However, how to calculate the effective stress tensor (or viscosity) is still a challenging problem for us. As an approximation, it is assumed that the effects of inhomogeneous structure on the stress tensor are of minor importance in comparison with the effective drag, as employed in many previous works [20, 43–45, 47]. In this article, therefore, the kinetic theory of granular flows (KTGF) is still used to close the viscosity and pressure of the solid phase, neglecting the effects of structure on the stress tensor, and the gas viscosity is assumed to be constant. The conservation equation of granular temperature is described by T1-5. In the current work, however, an algebraic approximation of KTGF was adopted to ensure quick convergence. Note that further works on the effective stress should be performed to include such effects as carried out in the filtered TFM [27–29]. How to unify this bubble-based drag and the corresponding stress into certain common knowledge is certainly an interesting issue deserving more research. A possible path is to reformulate the structure-based kinetic theory for heterogeneous gas-solid system. In the latest report [59], the EMMS group proposed their promising works on velocity distribution function as a first key step toward a generalized kinetic theory for gas-solid flow. As for the empiricism in the EMMS model, more work is required to reduce model’s dependence upon correlations. Furthermore, the structure size that is bubble diameter cannot be greater than the grid size of the simulation. Therefore, a restriction also needs to be imposed to resolve the heterogeneity smaller than the grid size. All these issues require further efforts.


### 3.2. Computing Description and Boundary Conditions

A lab-scale bubbling fluidized bed was selected to validate this new drag model, as shown in Figure3. FCC catalyst particles (belonging to Geldart group A) were selected as fluidized material. Physical properties of gas phase and particles are summarized in Table 1. For this set-up [60], the fluidizing column is 2.464 m high with 0.267 m inner diameter, while the disengaging section has a diameter of 0.667 m and a height of 1.745 m. In all simulations, the disengaging section was neglected to save the computational cost, because it does not have any effect of bubbling behavior. The carryover of particles from the bed is very small in bubbling regime as compared to turbulent regime. Therefore, the disengaging section was replaced by recirculating the entrained particles into the bottom inlet.Table 1
Summary of the physical properties of gas and particles.
Gas density (kg/m3) 1.225 Particle density (kg/m3) 1780 Gas viscosity (Pa·s) 1.7894 × 10−5 Particle diameter (m) 6.5 × 10−5 Minimum fluidization velocity (m/s) 0.003Figure 3
Schematic diagram of simulated bubbling fluidized bed with Geldart A particles.Hexahedral mesh with a grid size of 1 cm was used in the present simulations. The mesh-independence for the current gas-solid system has already been carried out in the literature [61]. Initially, the particles were loaded as a static bed. At the start of simulation gas flow is introduced to the bottom of column uniformly. The fluidizing gas leaves from a pressure outlet at the top, where atmospheric pressure was specified. A no-slip boundary condition was described at the wall for the gas, while a partial-slip boundary for the solid phase with a specularity coefficient of 0.6 was implemented, as specified in other works [20, 45, 46, 62]. A time step of 5 × 10−4 s was used to ensure quick convergence within 30 iterations for each time step. The convergence criterion between two iterations was set to 0.001. The statistical data sampling was carried out between 20 and 30 seconds of the physical time. This data was then used to compare with experimental data reported by Zhu et al. [60]. The detailed simulation settings are listed in Table 2.Table 2
Simulation settings used in FLUENT solver.
Time Unsteady, first-order implicit Viscous Laminar Particle-particle coefficient restitution 0.9 Pressure-velocity coupling Phase-coupled SIMPLE Momentum discretization Second-order upwind Volume fraction discretization Quick Granular temperature Algebraic KTGF Granular viscosity Gidaspow Granular bulk viscosity Lun et al. Frictional viscosity Schaeffer Angle of internal friction 30° Solids pressure Lun et al. Frictional pressure Based KTGF Radial distribution Lun et al. Specularity coefficient 0.6 Close packing density 0.63 Time step 0.0005 s Max iterations per time step 30 Initial static bed height 1.2 m Initial bed packing fraction 0.6 Inlet gas velocity 0.2 m/s


### 3.3. Results and Discussion

Figure4 shows time sequence snapshots of solid concentration using the present novel bubble-based drag model. It can be seen that this new drag can capture the visible bubbles surrounded by continuous emulsion phase and their complex movement. The visible bubbles appearing near the bottom inlet coalesce through the middle section and erupt at the top surface of the bed. The gas appears to form a larger bubble to pass through the fluidized bed. This tendency can be attributed to the fact that the weaker interphase drag promotes the coalescence of bubbles easily, as considered in our model with consideration of meso-scale structures. In this way, the compromise of flow dominance of bubbles against emulsion can be realized. That is to say that the gas finds its way with less resistance, while the particles always seek minimum gravitational potential. This explanation is consistent with dominant mechanisms of gas-solid flows, as indicated by Li and Kwauk [42].Figure 4
Instantaneous snapshots of solid concentration in bubbling fluidized bed (blue and red denote the dilute and dense flow regions, resp.).Figure5 shows a time series plot of local solid concentration at the center and near the wall for the last 5 seconds of simulation. When the solid concentration is lower than 0.2, the gas can be considered to be in a bubble phase. At the center region, therefore, the frequent low peaks may mean that there is strong oscillation of bubbles in that region. In other words, this hints that the gas mainly flows through the center of the bed in the form of bursting bubbles, whereas the emulsion tends to accumulate near the walls.Figure 5
Time series plot of instantaneous solid concentration near the wall and the center regions at two different heights: (a)h=0.4 m and (b) h=0.8 m.
(a) (b)Figure6 further provides a time series of local particle velocity at the center and near the wall at two different heights. The particles mainly move upward at the center, while they flow down near the wall. In general, the particle velocity fluctuates around 0.5 m/s. This is consistent with the experimental data and the assumption of ae=0 in our model.Figure 6
Time series plot of instantaneous axial particle velocity near the wall and the center regions at two different heights: (a)h=0.4 m and (b) h=0.8 m.
(a) (b)Figure7 displays the comparison of the experimental data and the axial profile of time-averaged solid concentration from the simulation. It can be seen that a dilute-top and dense-bottom distribution is visible here by using the new model. In contrast, the predicted results give a qualitative agreement with experiments, while coarse-grid TFM with homogeneous drag cannot capture such distribution at all, as reported many times in literature [20, 44–46]. To avoid unnecessary repetition, here we will not present the simulation with the homogeneous drag again. For details, the interested reader is referred to the relevant literature mentioned above. The comparatively larger deviation in the axial profile at the middle may be partly explained by the uncertainty in experimental measurements. Therefore, it is necessary to perform detailed study of characteristics of this novel drag model. More investigations for the validation of current model are underway and will be presented once available.Figure 7
Comparison of experimental and simulation results for axial solid profiles.Figure8 shows detailed comparison between simulated and experimental radial distribution of solid concentration at four different heights. It can be easily observed that the radial distribution reveals higher solid concentration near the wall and lower solid concentration at the center, as shown in the cross-section contour of mean solid concentration in Figure 8. The new drag model shows good agreement with the experimental data. Hence the model is able to reasonably predict the core annulus type flow structure.Figure 8
Comparison of experimental and simulation results for radial solid concentration profiles at four heights of the bed: (a)h=0.8 m; (b) h=1.1 m; (c) h=0.4 m; (d) h=0.6 m.


## 3.1. Governing Equations

In this work, the two-fluid model (TFM) is adopted to simulate bubbling fluidization. The CFD software, FLUENT® 6.3.26, was used as fluid solver to carry out all simulations. The governing equations of the two-fluid model for gas-solid flows are summarized as follows.Mass Conservation Gas phase:T1-1∂∂tεgρg+∇·εgρgug=0. Solid phase:T1-2∂∂tεsρs+∇·εsρsus=0.Momentum Conservation Gas phase:T1-3∂∂tεgρgug+∇·εgρgugug=-εg∇p+∇·τg+εgρgg-βug-us. Solid phase:T1-4∂∂tεsρsus+∇·εsρsusus=-εs∇p-∇ps+∇·τs+εsρsg+βug-us.Granular Temperature T1-5 3 2 ∂ ∂ t ε s ρ s Θ s + ∇ · ε s ρ s Θ s u s = - p s I + τ s : ∇ u s + ∇ · κ s ∇ Θ s - γ s - 3 β Θ s .Constitutive relations are as follows:Stress-strain tensor:
Gas phase:(10)τg=εgμg∇ug+∇ugT-23εgμg∇·ugΙ.Solid phase:(11)τs=εsμs∇us+∇usT+εsλs-23μs∇·usI.Solid phase pressure:(12)ps=εsρsΘs1+21+esεsg0.Solid phase viscosity:(13)μs=μs,col+μs,kin+μs,fr,μs,col=45εsρsdpg01+esΘsπ,μs,kin=10ρsdpΘsπ96εs1+esg01+45g0εs1+es2,μs,fr=pssin⁡Φ2I2D.Solid phase bulk viscosity:(14)λs=43εsρsdpg01+esΘsπ.Radial distribution function:(15)g0=1-εsεs,max1/3-1.Diffusion coefficient of granular energy:(16)κs=150ρsdpΘsπ3841+esg01+65εsg01+es2+2ρsεs2dp1+esg0Θsπ.Collision energy dissipation:(17)γs=121-es2εs2ρsg0Θs3/2dpπ.In the TFM, the continuity equations for the gas and the particle phases are given as T1-1 and T1-2. The momentum conservation for each phase is expressed by T1-3 and T1-4, where τ and β denote the stress tensor and the drag coefficient, respectively. The bubble-based drag coefficient formulated in this article can be adopted to calculate the effective drag coefficient which is incorporated into TFM through a user-defined function (UDF). However, how to calculate the effective stress tensor (or viscosity) is still a challenging problem for us. As an approximation, it is assumed that the effects of inhomogeneous structure on the stress tensor are of minor importance in comparison with the effective drag, as employed in many previous works [20, 43–45, 47]. In this article, therefore, the kinetic theory of granular flows (KTGF) is still used to close the viscosity and pressure of the solid phase, neglecting the effects of structure on the stress tensor, and the gas viscosity is assumed to be constant. The conservation equation of granular temperature is described by T1-5. In the current work, however, an algebraic approximation of KTGF was adopted to ensure quick convergence. Note that further works on the effective stress should be performed to include such effects as carried out in the filtered TFM [27–29]. How to unify this bubble-based drag and the corresponding stress into certain common knowledge is certainly an interesting issue deserving more research. A possible path is to reformulate the structure-based kinetic theory for heterogeneous gas-solid system. In the latest report [59], the EMMS group proposed their promising works on velocity distribution function as a first key step toward a generalized kinetic theory for gas-solid flow. As for the empiricism in the EMMS model, more work is required to reduce model’s dependence upon correlations. Furthermore, the structure size that is bubble diameter cannot be greater than the grid size of the simulation. Therefore, a restriction also needs to be imposed to resolve the heterogeneity smaller than the grid size. All these issues require further efforts.


## 3.2. Computing Description and Boundary Conditions

A lab-scale bubbling fluidized bed was selected to validate this new drag model, as shown in Figure3. FCC catalyst particles (belonging to Geldart group A) were selected as fluidized material. Physical properties of gas phase and particles are summarized in Table 1. For this set-up [60], the fluidizing column is 2.464 m high with 0.267 m inner diameter, while the disengaging section has a diameter of 0.667 m and a height of 1.745 m. In all simulations, the disengaging section was neglected to save the computational cost, because it does not have any effect of bubbling behavior. The carryover of particles from the bed is very small in bubbling regime as compared to turbulent regime. Therefore, the disengaging section was replaced by recirculating the entrained particles into the bottom inlet.Table 1
Summary of the physical properties of gas and particles.
Gas density (kg/m3) 1.225 Particle density (kg/m3) 1780 Gas viscosity (Pa·s) 1.7894 × 10−5 Particle diameter (m) 6.5 × 10−5 Minimum fluidization velocity (m/s) 0.003Figure 3
Schematic diagram of simulated bubbling fluidized bed with Geldart A particles.Hexahedral mesh with a grid size of 1 cm was used in the present simulations. The mesh-independence for the current gas-solid system has already been carried out in the literature [61]. Initially, the particles were loaded as a static bed. At the start of simulation gas flow is introduced to the bottom of column uniformly. The fluidizing gas leaves from a pressure outlet at the top, where atmospheric pressure was specified. A no-slip boundary condition was described at the wall for the gas, while a partial-slip boundary for the solid phase with a specularity coefficient of 0.6 was implemented, as specified in other works [20, 45, 46, 62]. A time step of 5 × 10−4 s was used to ensure quick convergence within 30 iterations for each time step. The convergence criterion between two iterations was set to 0.001. The statistical data sampling was carried out between 20 and 30 seconds of the physical time. This data was then used to compare with experimental data reported by Zhu et al. [60]. The detailed simulation settings are listed in Table 2.Table 2
Simulation settings used in FLUENT solver.
Time Unsteady, first-order implicit Viscous Laminar Particle-particle coefficient restitution 0.9 Pressure-velocity coupling Phase-coupled SIMPLE Momentum discretization Second-order upwind Volume fraction discretization Quick Granular temperature Algebraic KTGF Granular viscosity Gidaspow Granular bulk viscosity Lun et al. Frictional viscosity Schaeffer Angle of internal friction 30° Solids pressure Lun et al. Frictional pressure Based KTGF Radial distribution Lun et al. Specularity coefficient 0.6 Close packing density 0.63 Time step 0.0005 s Max iterations per time step 30 Initial static bed height 1.2 m Initial bed packing fraction 0.6 Inlet gas velocity 0.2 m/s


## 3.3. Results and Discussion

Figure4 shows time sequence snapshots of solid concentration using the present novel bubble-based drag model. It can be seen that this new drag can capture the visible bubbles surrounded by continuous emulsion phase and their complex movement. The visible bubbles appearing near the bottom inlet coalesce through the middle section and erupt at the top surface of the bed. The gas appears to form a larger bubble to pass through the fluidized bed. This tendency can be attributed to the fact that the weaker interphase drag promotes the coalescence of bubbles easily, as considered in our model with consideration of meso-scale structures. In this way, the compromise of flow dominance of bubbles against emulsion can be realized. That is to say that the gas finds its way with less resistance, while the particles always seek minimum gravitational potential. This explanation is consistent with dominant mechanisms of gas-solid flows, as indicated by Li and Kwauk [42].Figure 4
Instantaneous snapshots of solid concentration in bubbling fluidized bed (blue and red denote the dilute and dense flow regions, resp.).Figure5 shows a time series plot of local solid concentration at the center and near the wall for the last 5 seconds of simulation. When the solid concentration is lower than 0.2, the gas can be considered to be in a bubble phase. At the center region, therefore, the frequent low peaks may mean that there is strong oscillation of bubbles in that region. In other words, this hints that the gas mainly flows through the center of the bed in the form of bursting bubbles, whereas the emulsion tends to accumulate near the walls.Figure 5
Time series plot of instantaneous solid concentration near the wall and the center regions at two different heights: (a)h=0.4 m and (b) h=0.8 m.
(a) (b)Figure6 further provides a time series of local particle velocity at the center and near the wall at two different heights. The particles mainly move upward at the center, while they flow down near the wall. In general, the particle velocity fluctuates around 0.5 m/s. This is consistent with the experimental data and the assumption of ae=0 in our model.Figure 6
Time series plot of instantaneous axial particle velocity near the wall and the center regions at two different heights: (a)h=0.4 m and (b) h=0.8 m.
(a) (b)Figure7 displays the comparison of the experimental data and the axial profile of time-averaged solid concentration from the simulation. It can be seen that a dilute-top and dense-bottom distribution is visible here by using the new model. In contrast, the predicted results give a qualitative agreement with experiments, while coarse-grid TFM with homogeneous drag cannot capture such distribution at all, as reported many times in literature [20, 44–46]. To avoid unnecessary repetition, here we will not present the simulation with the homogeneous drag again. For details, the interested reader is referred to the relevant literature mentioned above. The comparatively larger deviation in the axial profile at the middle may be partly explained by the uncertainty in experimental measurements. Therefore, it is necessary to perform detailed study of characteristics of this novel drag model. More investigations for the validation of current model are underway and will be presented once available.Figure 7
Comparison of experimental and simulation results for axial solid profiles.Figure8 shows detailed comparison between simulated and experimental radial distribution of solid concentration at four different heights. It can be easily observed that the radial distribution reveals higher solid concentration near the wall and lower solid concentration at the center, as shown in the cross-section contour of mean solid concentration in Figure 8. The new drag model shows good agreement with the experimental data. Hence the model is able to reasonably predict the core annulus type flow structure.Figure 8
Comparison of experimental and simulation results for radial solid concentration profiles at four heights of the bed: (a)h=0.8 m; (b) h=1.1 m; (c) h=0.4 m; (d) h=0.6 m.


## 4. Conclusion

A novel bubble-based drag model for bubbling fluidized beds has been proposed in this work based on the description of heterogeneity within each computational cell. In this model, the effect of meso-scale bubbles on the drag force has been taken into account. Seven structural parameters have to be determined by solving a set of balance equations, involving mass conservation, force balance, and reliable correlations for predicting of bubble diameter and emulsion voidage. Then, a bubble-based drag coefficient has to be derived directly from local structure parameters. It is worthwhile to note that the slip velocity (or the gas and the solid velocities) plays an important role in addition to the local voidage in our model. For validation, this new drag was incorporated into the two-fluid model to simulate the hydrodynamic of Geldart A particles in a bubbling fluidized bed. The simulated results, using this new drag model, show a qualitative agreement with the experiments. More works are, however, required to verify and revise this new model.



---

*Source: 8256816-2016-10-27.xml*