Nouman Ahmad ,Jianqiang Deng,*,Muhammad Adnan

1 Shaanxi Key Laboratory of Energy Chemical Process Intensification,Xi’an 710049,China

2 School of Chemical Engineering and Technology,Xi’an Jiaotong University,Xi’an 710049,China

Keywords:Mesoscale Bubbles Energy minimization multiscale (EMMS)Heterogeneity index Bubbling fluidized bed

ABSTRACT Mesoscale bubbles exist inherently in bubbling fluidized beds and hence should be considered in the constitutive modeling of the drag force.The energy minimization multiscale bubbling(EMMS/bubbling)drag model takes the effects of mesoscale structures (i.e.,bubbles) into the modeling of drag coefficient and thus improves the coarse-grid simulation of bubbling and turbulent fluidized beds.However,its dependence on the bubble diameter correlation has not been thoroughly investigated.The hydrodynamic disparity between homogeneous and heterogeneous fluidization is accounted for by the heterogeneity index,Hd,which can be affected by choice of bubble diameter correlation.How this choice of bubble diameter correlation influences the model prediction calls for further fundamental research.This article incorporated seven different bubble diameter correlations into EMMS/bubbling drag model and studied their effects on Hd.The performance of these correlations has been compared with the correlation used previously by EMMS/bubbling drag model.We found that some of the correlations predicted lower Hd by order of a magnitude than the correlation used by the original EMMS/bubbling drag.Based on such analysis,we proposed a modification in the EMMS drag model for bubbling and turbulent fluidized beds.A computational fluid dynamics (CFD) simulation using two-fluid model with the modified EMMS/bubbling drag model was performed for two bubbling and one turbulent fluidized beds.Voidage distribution,time averaged solid concentration and axial solid concentration profiles were studied and compared with the previous version of the EMMS/bubbling drag model and experimental data.We found that the right choice of bubble diameter correlations can significantly improve the results for CFD simulations.

1.Introduction

Almost a century has passed since the first commercial use of a fluidized bed in the 1920s (Winkler’s coal gasifier) [1].But there are still many challenges that are affecting the fluidization processes in the industry[2].Gas-solid flows are being heterogeneous in nature and manifest inhomogeneity at different time and length scales.In terms of heterogeneous structures,the fluidized beds are characterized as having discrete gas voids called“bubbles”or have strands of particles called “clusters”[3].As a common trait,both bubbles and clusters exhibit an inherent instability that can affect reactor performance[4,5].Both bubbles and clusters can be classified in terms of their size,shape,etc.,and since their sizes are larger than the micro-scale,i.e.,particle diameter,and lesser than the macro-scale,i.e.,reactor diameter,they are termed as mesoscale structures[6].These mesoscale structures are key features to fully understand the hydrodynamics of a fluidized bed reactor [7,8].

The bubble is a typical heterogeneous structure that appears in a low velocity bubbling fluidized bed.When gas velocity exceeds the minimum fluidization velocityUmf,the bed enters a bubbling regime [1].The voids containing no or very low particle content appear in the bed;these voids are referred to as “bubbles”[9].These bubbles considerably affect hydrodynamics and govern the efficiency of the gas-solid fluidized systems [10].Every essential feature of gas-solid fluidized beds like heat and mass transfer,elutriation and particle entrainment,particle and gas residence time and chemical reaction conversion rates,etc.,all are significantly affected by bubble characteristics [11].These characteristics comprise bubble velocity,shape,and size.Thus a lot of research has been dedicated to experimentally measuring and characterizing the bubble growth phenomena in fluidized beds and developing a reasonable correlation to predict the bubble size and velocity.Many correlations have appeared in the literature,but each of these correlations has a specific applicability range.

Yasui and Johanson [12] were among the pioneer researchers who studied bubble characteristics and developed a bubble size correlation based on their experimental data.They used two different fluidized bed columns with internal dia.of 10 and 15 cm.Both of the columns were employed by a light-transmitting probe and were equipped with a porous plate distributor.The particle sizes used in their experiment were in the range of 41-450 μm.They found that the bubble size near the distributor is too small to be captured by their probe,and also,the bubble velocity was constant throughout the column.Instead of excess gas velocity,they used the reduced gas velocity term defined as the ratio of superficial gas velocityUgto minimum fluidization velocityUmf.Kobayashiet al.[13] also used reduced gas velocity,correlated their experimental data,and formulated their correlation.Kato and Wen [14]established in their study that the bubble size also depends on the distributor type,and they modified Kobayashiet al.’s[13]equation for fluidized beds with the perforated distributor.Later,Geldart [15] mentioned that using reduced gas velocity is not a very good way to find bubble size because it can affect the bubble diameter calculation,and instead,the excess gas velocity is a better choice.He also gave two different diameter correlations for porous and perforated distributors.Werther and Molerus [16,17] questioned the bubble size measuring methods,and they raised the fact that the presence of big probes may change the hydrodynamics around them in the bed,i.e.,they are intrusive and may give inaccurate measurements.They used the needle-type small identical probes mounted at different locations in the bed to prove their point and compare their power spectral density signals.They found no difference in the signals and deduced that their probes do not affect the surrounding hydrodynamics.They found out that the initial bubble diameter,i.e.,db0is different for Geldart A and B particles,and hence they gave their correlation with different initial bubble sizes.Dartonet al.[18] introduced the concept of the bubble growth pattern in their studies.They proposed that as the bubble moves upwards towards the bed surface,it increases in size.Dartonet al.[18] did not consider the bubble splitting phenomena in their proposed correlation and as a result,their correlation gives a constant increasing pattern in bubble size.This limitation was overcome by Horio and Nonaka [11] and they proposed a correlation which was a balance between coalescence and splitting of bubbles.

Nevertheless,bubble size is an important and essential aspect in the design of a bubbling fluidized bed,and it can affect the overall hydrodynamics of a fluidized bed in many ways.Karimipour and Pugsley [19] pointed out that the drawback of these correlations lies in the fact that they were established under specific operating conditions,and they are hence best suitable for those specific conditions and yet fail to give results under more generic conditions.They emphasized the fact that one correlation deduced by using a particular particle size may give a biased result for another different particle size.For example,three average particle sizes were used by Parket al.[20] in which only one size,i.e.,80 μm was belonging to the Geldart A particles category while the remaining two were of the Geldart B categoryi.e.,154 μm and 340 μm.Similarly,Horio and Nonaka[11]formulated their bubble diameter correlation based on five experimental data series of which only one data series was belonging to Geldart A classification.These correlations might not be suitable for other bigger or smaller particle sizes.Hence,a careful investigation is needed to select the bubble diameter correlation as it affects the hydrodynamics in bubbling fluidized bed like bubble volume fraction,bed expansion,bubble rise velocityetc.

Fortunately,due to the recent advancements in computational power and resources,computational fluid dynamics (CFD) has become a valuable tool in studying gas-solid flow hydrodynamics.There are two major approaches in CFD to study hydrodynamics in fluidized beds:the Euler-Euler approach and the Euler-Lagrangian approach.In the Euler-Euler approach,as the name suggests,the gas and solid phases are solved on an Eulerian framework by considering them as a continuum.The conservation equations are averaged,and some locally averaged quantities are obtained.This averaging process neglects some of the quantities at the level of element,but their consequences can appear through constitutive relations.These constitutive relations,like solid stress,drag force,etc.,need to be derived separately to close the set of equations.In the Euler-Lagrangian approach,the gas phase is treated in a similar way as in the Euler approach,and for the solid phase,the motion of each particle is tracked employing a Lagrangian grid.The Eulerian-Eulerian method (namely,two-fluid model,TFM)needs the least computational cost and is thus commonly used for industrial applications.

To describe the mesoscale structures and their effect on hydrodynamics,the energy minimization multi-scale (EMMS) drag model was introduced[6].The EMMS drag model was initially proposed for describing the steady-state structure for gas-solid fast fluidization.The EMMS drag model separated the whole system into three phases,i.e.,dense,dilute,and interphase.There are eight variables to describe the multi-scale system,i.e.,superficial gas velocityUgc,superficial solid velocityUsc,and voidage εgcfor the dense phase,Ugf,Usfand εgffor the dilute phase,and volume fraction of the dense phase (or clusters)fand the cluster diameterdclfor the interphase.However,only six equations were available;thus,one stability condition was proposed to close these equations.The stability condition was then optimized under the constraints of set of force and mass balance equations.

The EMMS-based drag modeling has been first well applied to the sub-grid level for fast fluidization by introducing the acceleration terms in the force balance equations[21-26].A heterogeneity index was introduced to quantify the difference between homogenous and heterogeneous drags [26].Later Shiet al.[27] and Honget al.[28,29]extended it to dense fluidization by replacing the clusters diameterdclwith bubble diameterdbto better characterize the mesoscale structure in bubbling fluidization.Luoet al.[30]extended the EMMS/bubbling model [28] to a sub-grid level and tried to get a mesh-independent solution.

Although the EMMS/bubbling drag model [28] has found successful applications in both bubbling and turbulent fluidized beds[27,28,30-34],relevant parameter analysis has not been fully investigated.In particular,an empirical correlation was employed to close the bubble diameter,and how this choice of bubble diameter correlation influences the model prediction calls for further fundamental research and will be presented here.There are many bubble diameter correlations available in the literature,which one is reasonable,and how to choose the suitable one will be discussed in this work.Moreover,how the choice of bubble diameter correlation can affect the overall hydrodynamics of the fluidized bed will be addressed in this work.

We aim to evaluate the performance of several widely used bubble diameter correlations by integrating them into the EMMS/-drag model[28].In the second section of this paper,we give a brief description of EMMS/bubbling drag model with emphasis on different bubble diameter correlations integrated in it.In third section,the impact of bubble diameter correlations on the heterogeneity index,Hdwill be discussed in detail.This section is important in understanding that howHdwith different bubble diameter correlations will change as a function of voidage and whether it will help capture the meso-scale bubbles in the system or not.A suitable bubble diameter correlation will then be adopted to be integrated into EMMS/bubbling drag model and the related structural parameters like volume fraction of dense phase voidage εgc,dense phasef,gas velocities in dense and dilute phase(Ugc,Ugf),particle acceleration term,etc.will be thoroughly discussed.Then in the fourth section,description for experimental data for model validation and CFD simulation model and setup will be presented.Finally,a series of simulations will be carried out with a modified version of the EMMS/bubbling drag model for two bubbling and one turbulent fluidized bed and bed’s hydrodynamics compared to experimental data will be presented in Section 5.

2.EMMS/Bubbling Model

In this section,we will briefly discuss the EMMS/bubbling model and its mathematical equations.For the detailed discussions about model derivation and solution schemes refer to the literature[28].

2.1.Multi-scale description

Following the original EMMS model[6],the overall system was also divided into dense phase (or emulsion),dilute phase (inside bubbles),and interphase(bubbles),as shown in Fig.1.The voidage inside the bubbles and dense phases are εgfand εgc,respectively,and the volume fraction of the dense phase is denoted byf(the volume fraction of bubbles is denoted by 1-f).The dense phase (denoted by subscript c) and dilute phase (denoted by subscript f)are further divided into four sub-structure continua,i.e.,dilutephase gas (gf),dilute-phase solid (sf),dense-phase gas (gc),and dense-phase solid (sc).The superficial gas velocities in the dense phase and dilute phases areUgcandUgf,respectively.The superficial solid velocities in dense and dilute phases areUscandUsf,respectively.The interphase (bubbles) rises at the speed ofUbthrough the surrounding media(the dense phase having the superficial velocityUe).The superficial gas velocityUgand the superficial solid velocityUsare the operating parameters.For the unsteadystate flow,different acceleration terms for dilute-phase and dense-phase particles (af,ac) and the acceleration term for mesoscale phase (ai) are respectively introduced in the momentum equations.

2.2.Hydrodynamic equations

The mathematical formulation of the EMMS/bubbling drag model [28] is presented below.It should be noted that only the scalars of vertical components are used here.In all,there are seven equations (Eqs.(1)-(3) and Eqs.(6)-(9)) and a stability condition(Eq.(10)) to close ten variables (Ugc,Usc,εgc,asc,Ugf,Usf,εgf,agf,f,anddb).

The momentum balance equations for the dense phase,dilute phase,and interphase are as follows:

Pressure balanced equation for different phases:

Hereaican be related with the other acceleration terms through the pressure-drop balance equation:

Mass balance equation for the gas and solid phases can be expressed as:

The mean voidage is related to the dense-phase and dilutephase voidages by

The closure of Eq.(3)i.e.meso-scale drag force requires a bubble diameter(db)correlation.As discussed in the introductory section,there exist many bubble diameter correlations in the literature.An extensive effort is thus required to determine the suitable correlation that can be used for wide flow regimes.

The correlation of bubble diameter used by the EMMS/bubbling drag model [27-29] is:

Horio and Nonaka [11]:

Fig.1.System resolution for bubbling fluidized bed reproduced from Hong et al.[28].

Besides Horio and Nonaka’s [11] correlation,we integrated the following bubble diameter correlations one by one from the literature in the EMMS/bubbling drag model [28] and studied their effect on model structural parameters.Changing the bubble diameter correlation will have an effect on meso-scale drag force and hence it will change the stability condition.The summary of these bubble diameter correlations is given in Table 1.The discussion on structural parameters will be discussed in the next section.

Table 1 Different correlations of bubble diameter and their corresponding applicable ranges

Romero and Johanson [35]:

Whitehead and Young [36]:

Geldart [15]:

Werther [37]:

Mori and Wen [38]:

Yacono [39]:

Caiet al.[40]:

The stability condition:

There are seven main equations constituting the EMMS drag model(Eqs.(1)-(3)and Eqs.(6)-(9))and 10 independent variables.To solve this optimization problem a stability condition(Nst=min)was proposed [6].In gas-solid upward flow the particles and the gas cannot dominate each other and there exists a compromise.The gas suspends the particles with least energy consumption and flows through particles with minimum possible resistance.This compromise results in a tendency of minimal energy consumption for suspending and transporting particles with respect to unit mass of particles.Later for general purposes,the stability condition is defined by a dimensionless ratio ofNst/NT[41]whereNTis the total mass-specific energy consumption.In summary,the change in bubble diameter correlation is going to affect the meso-scale drag force which consequently will affect the stability condition.

The above model equations are assumed to satisfy the stability condition in terms of the mass-specific energy consumption rate for suspending particles.

The model can be solved using each of the above-mentioned bubble diameter correlations integrated with the EMMS/bubbling drag model.A total of 8 solutions could be obtained following the solution scheme of Honget al.[28].The flow sheet of the numerical scheme of the above model is given in Fig.2.

3.Model Evaluation with Different Bubble Diameter Correlations

3.1.Heterogeneity index Hd

Wen and Yu [41] proposed an empirical correlation that is based on the voidage function (εg),which reveals how the drag force is affected by the bed voidage in homogenous fluidized beds.The Wen and Yu [41] coefficient can be expressed as:

wheredsis the solid particle diameter,ρgis the gas density,ugand usare the velocities of gas and solid phases,respectively.The drag coefficient of the single-particleCd0is defined as

whereReis the Reynolds numberi.e.Re=and μgis the dynamic viscosity of the gas.It should be noted that the Wen and Yu [41] drag correlation was derived from the experimental data of almost homogenous fluidized bed systems,and hence they should not be applied directly to study the hydrodynamics of heterogeneous fluidized beds.

Wang and Li [26] proposed aHdto quantify the effects of heterogeneous structures on the drag coefficient.Hdcan be defined as a ratio of mean drag force in a heterogeneous system to that in a homogeneous system when the voidage and slip velocity are the same for both cases.

Here βeis the structure-dependent drag coefficient and is defined as [28]:

We incorporated different bubble diameter correlations into the EMMS/bubbling model,and then the structure-dependent drag coefficient was calculated from Eq.(15).TheHdwas calculated from Eq.(14),and the effect of different bubble diameter correlations on it is presented in the next section.

3.2.Effects of bubble size correlation on Hd

The effect of different bubble diameter correlations mentioned in Table 1 on the heterogeneity indexHdobtained by coupling them with the EMMS/bubbling model [28] will be discussed here.The material properties used for this analysis are shown in Table 2.

Fig.2.The flow sheet about the numerical scheme of the EMMS/bubbling drag model.

Since the previous versions of the EMMS/ bubbling drag model[27,28]used the bubble diameter correlation of Horio and Nonaka[11],so we call it our base case and compare the heterogeneity obtained with other diameter correlations(Table 1)with this base case.It should be noted that some bubble diameter correlations showed a largerHdthan the base case,and some showed a lowerHdthan the base case.For clarity,we showed the comparison in two figures so that the results would not be merged into each other.

Table 2 Summary of material properties in experimental data of Zhu et al.[42]

Fig.3(a) below shows the comparison ofHdobtained from the EMMS/bubbling drag model [28] coupled with bubble diametercorrelation of Horio and Nonaka [11] which is our base case,withHdobtained from the bubble diameter correlations of Romero and Johanson [35],Geldart [15],Mori and Wen [38],and Caiet al.[40].These above-mentioned correlations showed the largerHdas compared to the base case.

The heterogeneity in the case of Romero and Johanson [35]dropped until the voidage was increased up to 0.59 approximately,and then it started to increase again till voidage reached 0.98.After this point,theHdstarted to drop again,which is strange because the bed must have reached the homogenous state.As compared to the base case,the correlation of Romero and Johanson [35]showed a minimal drop inHd,showing a quick transition between heterogeneous and homogenous behavior.The correlation of Geldart [15] showed a continual drop inHduntil around voidage of 0.55 and then started to increase gradually but failed to reach unity and instead started to decline again at a voidage of around 0.96.Quantitatively the correlations of Mori and Wen [38],and Caiet al.[40]showed the drop inHdwithin most of the voidage range like that of the base case.In comparison to the correlation of Horio and Nonaka [11],these correlations failed to converge to unity at the later end of voidage.

Fig.3(b)shows the comparison ofHdobtained from the EMMS/bubbling drag model [28] coupled with bubble diameter correlation of Horio and Nonaka [11] (base case),withHdobtained from the bubble diameter correlations of Yacono [39],Whitehead and Young[36]and Werther[37].These above-mentioned correlations showed the lowerHdas compared to the base case.The correlation of Yacono [39] showed a drop inHduntil the voidage of around 0.47 and started to increase again with the increase in voidage until 0.99.The correlations of Whitehead and Young [36] and Werther [37] almost showed the same quantitative trend,i.e.,a lowerHdby order of magnitude as compared to the base case.In their case,Hddropped until at the voidage of around 0.43 and then started to increase again towards unity.The correlation of Werther[37]converged to unity at the later end of the voidage,but the correlation of Whitehead and Young [36] failed to converge after around voidage of 0.99.This problem of convergence at the later end of the voidage can be solved by lettingHd=1 at some higher values of voidage (εg）for example 0.9991 or above in CFD simulations.Nevertheless,the correlations of Whitehead and Young [36]and Werther[37]showed a lowerHdas compared to the base case reflecting that the formation of heterogeneous structures results in the reduction of drag coefficient.

The above analysis showed that the bubble diameter correlations have a profound effect on predicting the heterogeneity index.Thus,we propose a modification in EMMS/bubbling drag by integrating it with the correlations like that of Whitehead and Young[36] and Werther [37] to predict a betterHd.In the next section,we will perform some structural parameter analysis on the modified EMMS/bubbling drag model.Since the correlations of Whitehead and Young [36] and Werther [37] showed an almost identical trend,here we chose the correlation of Werther [37] for further analysis.

3.3.Discussion on structural parameters for the modified EMMS/bubbling drag model

After obtaining the model solution for the modified EMMS/bubbling drag model,the related structural parameters like volume fraction of dense phase voidage εgc,dense phasef,gas velocities in dense and dilute phase (Ugc,Ugf),particle acceleration term,etc.are very important.We will examine the variation trend of these parameters with overall voidage εg.Here,εgis the overall voidage of unit volume and within the range (εmf,εmax).

3.3.1.The volume fraction of the dense phase and the dense phase voidage

Fig.4(a)shows the variation of the dense phase voidage εgcwith the mean voidage εg.With the increase of εg,initially,the εgcshowed a constant value that was equal to the minimum fluidization voidage εmf(here is 0.4).This means that the bed was not homogeneously expanded at the beginning of the fluidization,and instead,void or bubbles may form.After around the mean voidage value of 0.44,the dense phase voidage started to increase until catching up with the mean voidage at the dilute end,indicating that more and more gas gradually entered the emulsion phase.The change in the trend of the volume fraction of the dense phase is associated with the dense phase voidage.This trend is pictured in Fig.4(b),where it can be seen that as the voidage εgincreased,the dense phase volumefdecreased consequently.After around a voidage of 0.45,thefstarted to level off,but εgckept on increasing to maintain the bed expansion.

3.3.2.Gas velocities in dense and dilute phases

Fig.5 shows the variation of the superficial gas velocities in dense and dilute phases with mean voidage.Fig.5(a) is for the gas velocity in the dilute phase and showed that the gas velocity decreased with the increase of the mean voidage because more gas entered the emulsion phase.Fig.5(b) shows that gas velocity in the dense phase declined at the beginning of the fluidization;meanwhile,the dilute phase’s gas velocity also decreased,indicating that more bubbles were formed (the volume fraction of dilute phase 1-fincreases).Above the voidage of about 0.44,the gas velocity in the dense phase started to increase with the mean voidage.Generally,we can know from the figure that mostly the gas was bypassing the dense phase and making room for itself to make its way through minimum resistance hence achieving a compromise of flow dominance.This analysis is consistent with previous studies [21,26].

3.3.3.Acceleration term

Fig.6 shows the variation of particle acceleration term in the dense phase with the mean voidage.It can be seen that the particle acceleration decreased gradually as the mean voidage increased.This means that gas velocity was sufficient to push the particles in the upward direction towards higher acceleration.The acceleration term in the dilute phase is assumed to be equal to gravitational acceleration,and as a result,the drag force exerted by dilute phase solid particles can be neglected.It should be noted that Wang and Li [26] pointed out the acceleration terms do not equate to the real particle accelerations.Their values are indicators of the non-equilibrium states where a high value of acceleration means a higher non-equilibrium system.

3.3.4.Variation of Nst/NT with voidage

Li and Kwauk [6] have explained that in gas-solid concurrentup flow in a riser,there exists a compromise between the gas and particles.The gas tends to flow in the upward direction with the least resistance,and particles tend to arrange themselves with minimal potential energy.As a result of this compromising arrangement,the energy consumed to suspend and transport the particlesNst,i.e.,Eq.(10) (stability condition),is minimized.Later,Ge and Li[43]mentioned that it is better to characterize the stability condition as a ratio ofNstandNT(Nst/NT) whereNTis the total mass-specific energy consumption,i.e.,Eq.(11).

Fig.3.Comparison of heterogeneity indexes functions with respect to voidage using the EMMS/bubbling model coupled with different bubble diameter correlations.

Fig.4.(a) Voidage and (b) volume fraction of the dense phase vs.the mean voidage.

Fig.6.Variation of the dense phase particle acceleration with the mean voidage.

Fig.7 gives the variation of energy consumption ratio (Nst/NT)with the mean voidage εg.Firstly as the voidage increased,theNst/NTdropped until around the voidage εg=0.44.After this point,theNst/NTgradually increased with an increase in voidage.This variation is possibly due to the heterogeneous nature of the flow.The sudden drop in the energy consumption until around voidage εg=0.44 is because,at first the system showed heterogeneity.After a further increase of voidage εgbeyond 0.44,the energy consumption ratioNst/NTstarted to increase,showing that the flow heterogeneity was gradually reduced with the increase of the voidage.

Fig.7.Variation of energy consumption ratio with the mean voidage.

4.Validation and Numerical Simulation

4.1.Experimental data

The experimental setups for validation that we chose for this study are given in Fig.8 below.

Fig.8.Schematics of the experimental setup for bubbling fluidized beds of(a) Zhu et al.[42],(b) Dubrawski et al.[44],and (c)turbulent fluidized bed of Venderbosch [45].

For bubbling fluidized beds,we used the experimental data of Zhuet al.[42] and Dubrawskiet al.[44].In their experimental setup,Zhuet al.[42] used a Plexiglas column with an internal diameter of 0.267 m,and the height between the top and distributor is 2.464 m.They used an aluminum perforated plate distributor to feed the air.An orifice plate was also employed to measure the gas flow rate.The particles were fluidized upward with a gas velocity of 0.2 m·s-1.They used Geldart A particles(FCC particles)with a Sauter mean diameter of 65 μm and density of 1780 kg·m-3.Dubrawskiet al.’s [44] experimental setup consisted of a traveling fluidized bed column with inner cyclones integrated into it.They used particles with a Sauter mean diameter of 103 μm and density of 1560 kg·m-3.

For turbulent fluidized bed,we used the experimental data of Venderbosch [45].In their experimental set up they used a Pyrex glass apparatus.The fluidizing column was 0.75 m long with an internal diameter of 0.05 m.The column’s top section had an increased cross-sectional area that allowed the entrained particles to fall back into the bed.They used Geldart A particles with a Sauter mean diameter of 90 μm and a density of 1375 kg·m-3.The other relevant properties and experimental operating conditions are given in Table 3.

Table 3 Summary of material properties in experimental data

4.2.Model selection

All the numerical simulations were carried out by using the continuum model,also called the two-fluid model available in commercial software Ansys Fluent?14.5.The air was the primary phase,and the solid particles were the secondary phase.The governing and constitutive equations are summarized in Table 4.Related model selection will be given in the next section.

4.3.Simulation setup

The 2D schematics of the three selected fluidized beds are shown in Fig.9.As a previous study [46] showed that both 2D and 3D simulations could capture typical qualitative flow features in fluidized beds,we restricted our study to 2D simulations for saving computational cost.

Fig.9.2D geometry of bubbling fluidized bed of (a) Zhu et al.[42],(b) Dubrawski et al.[44],and (c) turbulent fluidized bed of Venderbosch [45].

For the bubbling fluidized bed of Zhuet al.[42],the fluidized column was 2.464 m in height and 0.267 m in inner diameter.The disengaging section was neglected to save computational time.Particles with minimum fluidization voidage were initially packed up to an initial height of 1.2 m.The gas (0.2 m·s-1) uniformly flowed into the bottom of the bed and left from the top outlet,where atmospheric pressure was prescribed.For the bubbling fluidized bed of Dubrawskiet al.[44],the fluidized column was 0.96 m in height and 0.133 m in inner diameter.The disengaging section was 0.59 m in height above the fluidizing column.Particles were initially packed up to an initial height of 0.8 m.The gas velocity was set at 0.4 m·s-1according to the experimental conditions.

For the case of a turbulent fluidized bed of Venderbosch [45],the main section of the fluidized column was 0.75 m in height and 0.05 m in diameter,above which an expanded bed section(ID 0.1 m,0.3 m long)was also included to keep particles from serious entrainment.The initial height of packed particles with a solid concentration of 0.5625 was set to be at 0.2 m.

For all the cases,the SIMPLE algorithm was used to achieve the pressure-velocity coupling.All of the modified EMMS drag models with different bubble diameter correlations were integrated into the solver through user-defined functions (UDFs) subroutines using C language.The fitting functions of EMMS drag (integrated with different bubble diameter correlations)of all the cases applied in this paper are listed in Supplementary Material.The kinetic theory of granular flows(KTGF)[47]was used to close the solid pressure and viscosities,thus neglecting the structure-dependent stress and pressure.The algebraic form of the granular temperature equation was used,which has been found to save computational time while also maintaining a reasonable precision[23,32,48].The kinetic theory of granular flow (KTGF) [47] was applied to close the particle stress terms.The no-slip boundary condition and partial-slip boundary conditions were applied for the gas and solid phases.The solid specularity coefficient and restitution coefficient of particle collision for both cases were set to be 0.6 and 0.9,respectively [27].All simulations were performed for 20s [28,29],of which,first 10 s were to reach the steady-state,and the data statistics were performed on the later 10 s.Relevant simulation settings are shown in Table 5.

Table 4 Summary of governing and constitutive equations

Table 5 Simulation settings

4.4.Grid size selection

To achieve the best computational cost and numerical accuracy,grid size selection is critical in CFD simulations.Square grids were generated uniformly by using Gambit?2.4,and we used the structured computational grid size of about 5 mm for bubbling fluidized bed.This grid size has been reported to be mesh-independent by previous researches [49,50].

The independent mesh study was then performed for a turbulent fluidized bed.Three different grid sizes (4,5,8 mm) were tested.This study aimed to identify the largest possible grid size with no significant influence on the numerical results.All the simulations for this study were carried out with the EMMS/bubbling drag model [28] integrated with the Werther’s bubble diameter correlation.The data statistics for each simulation case were performed after the initial 10 s of the simulation to ensure statistical steady-state behaviour inside the bed.The instantaneous solid concentration and time-averaged axial solids concentration profiles predicted using three different grid resolutions are shown in Fig.10.All selected grid sizes gave nearly the same solid volume concentration along with the bed height.Thus,keeping in mind the computational cost and numerical accuracy,in the further sections,we will use the 4 mm grid size for simulations of a turbulent fluidized bed.

Fig.10.Instantaneous solid concentration contours (left) and time-averaged axial solids concentration profiles (right) predicted from different grids.

5.Results and Discussion

5.1.Bubbling fluidized bed (Zhu et al.[42])

5.1.1.Instantaneous voidage distribution

Fig.11 shows the instantaneous voidage distributions predicted by the original EMMS/bubbling drag model (i.e.,base case with Horio and Nonaka’s[11]correlation)and the modified EMMS/bubbling drag model coupled with different bubble diameter correlations.The voidage distributions plotted in Fig.11 are for those bubble diameter correlations that predicted a larger heterogeneity index than the base case(see Fig.3(a)).It can be seen that the correlations of Romero and Johansson[35]and Geldart[15]over predicted the bed expansion height and showed more uniform voidage distributions without clear separation of the bubble and emulsion phases as compared to the base case.The reason for this discrepancy (with the base case) is that of the larger gas-particle interaction force of these correlations,which is also evident from the heterogeneity index curves (Fig.3(a)).The correlations of Caiet al.[40]and Mori and Wen[38]improved the distributions compared with Romero and Johansson [35] and Geldart [15],as their heterogeneity index curves were closer to the base case.However,the bed expansion height predicted from these correlations was slightly higher than the base case because of their higher drag force.

Fig.11.Instantaneous voidage distribution for bubbling fluidized of Zhu et al.[42]using the EMMS/bubbling drag model coupled with bubble diameter correlations that showed larger Hd than the base case.

Fig.12 shows the comparison of instantaneous voidage distributions between the base case and the correlations of Yacono[39],Werther [37],and Whitehead and Young [36].These correlations gave a lower heterogeneity index than the base case (see Fig.3(b)).Thus,these three correlation’s predicted results showed more heterogeneous distributions with a clear separation of the bubble and emulsion phases,even better than that of the base case.The heterogeneity index calculated from these correlations is lower than that of the base case,hence reflecting a lower drag force.Here the modified EMMS model integrated with the correlation of Werther [37] gave the lowest bed height,and mesoscale structures in the forms of bubbles or voids are well captured.

Fig.12.Instantaneous voidage distribution for bubbling fluidized of Zhu et al. [42]using the EMMS/bubbling drag model coupled with bubble diameter correlations that showed lower Hd than the base case.

5.1.2.Instantaneous and time-averaged solid concentration distributions

Fig.13 shows the comparison of instantaneous and timeaveraged solid concentration distributions between the base case and the correlations of Romero and Johanson [35],Geldart [15],Mori and Wen [38],and Caiet al.[40].The correlation of Romero and Johanson [35] could not predict the typical characteristics of the bubbling bed (i.e.,the dense region at the bottom and freeboard region at the top of the bed) and showed uniform distributions in the whole domain because of its larger gas-particle interaction force.Meanwhile,the distributions predicted from the correlations of Geldart [15],Mori and Wen [38],and Caiet al.[40]were improved to some extent as compared with the Romero and Johanson[35]because of lower gas-solid interaction force (as evidenced from the heterogeneity index curve).However,the bed expansion height predicted from these correlations was still larger than that of the base case.

Fig.13.Instantaneous and time-averaged distribution of solid concentration of the bubbling fluidized bed of Zhu et al.[42]by using the EMMS/bubbling drag model coupled with bubble diameter correlations that showed larger Hd than the base case.

Fig.14 represents the comparison of instantaneous and timeaveraged solid concentration distributions obtained from the base case and the correlations of Yacono [39],Whitehead and Young[36],and Werther [37].The solid concentration distributions predicted from these correlations were similar to the base casei.e.,a dense concentration at the bottom,dilute concentration at the top,dense concentration near the walls,and dilute in the center.But a typical core annulus structure was better captured by the correlation of Whitehead and Young[36]and Werther[37]as compared to the base case(i.e.,Horio and Nonaka’s[11]correlation).It can be seen that the concentration in the case of Whitehead and Young[36]and Werther[37]was denser at the wall(as compared to the base case) and dilute at the center.Meanwhile,the bed expansion height reported from the correlations of Whitehead and Young [36] and Werther [37] was also lower than the base case.

Fig.14.Instantaneous and time-averaged distribution of solid concentration of the bubbling fluidized bed of Zhu et al.[42]by using the EMMS/bubbling drag model coupled with bubble diameter correlations that showed lower Hd than the base case.

5.1.3.Time-averaged axial solid concentration profiles

Fig.15 shows the comparison of time-averaged axial distribution of solid concentration in bubbling fluidized bed obtained from the base case and the correlations of Romero and Johanson [35],Geldart [15],Mori and Wen [38],and Caiet al.[40].The crosssectionally averaged axial solid concentration predicted by using Romero and Johanson [35] and Geldart [15] was around 0.25 and 0.27,respectively,which was much less than that predicted by the base case and experimental data.Mori and Wen [38] ’s correlation failed to predict the solid concentration in the lower bed than the base case but gave a better prediction at the upper part of the bed.The correlation of Caiet al.[40] gave a better overall prediction,but the bed height is slightly higher than the base case.

Fig.15.Time-averaged axial solids concentration profile for bubbling fluidized bed of Zhu et al. [42] predicted from the EMMS/bubbling drag model coupled with bubble diameter correlations that showed higher Hd than the base case.

Fig.16 gives the comparison of time-averaged axial distribution of solid concentration in bubbling fluidized bed obtained from the base case and the correlations of Yacono [39],Whitehead and Young[36],and Werther[37].With Yaccono[39],a lower concentration was again well captured as compared to the base case and a slightly larger bed expansion height was found.The correlations of Whitehead and Young [36] and Werther [37] better predicted the solid concentration in the lower bed and agreed well with the experimental data as compared to the base case.The relative error is calculated between experimental and simulated results and with base case is 3.8% while with the modified version is 2.3%approximately.

Fig.16.Time-averaged axial solids concentration profile for bubbling fluidized bed of Zhu et al.[42] predicted from EMMS/bubbling drag model coupled with bubble diameter correlations that showed lower Hd than the base case.

5.2.Bubbling fluidized bed (Dubrawski et al.[44])

In this section,we only used the correlation of Werther[37]and compared it with the base case.Fig.17 represents the comparison of instantaneous and time-averaged solid concentration distributions obtained from the base case and the correlation of Werther[37].Again a typical core annulus structure was better captured by the correlation of Werther [37] as compared to the base case(i.e.,Horio and Nonaka’s [11] correlation).It can be seen that the concentration in the case of Werther [37] was denser at the wall(as compared to the base case)and dilute at the center.Meanwhile,the bed expansion height reported from Werther [37] was also lower than the base case.

Fig.17.Instantaneous and time-averaged distribution of solid concentration of the bubbling fluidized bed of Dubrawski et al. [44] by using the EMMS/bubbling drag model coupled with bubble diameter correlation of Werther [37].

Fig.18 shows time-averaged axial distribution of solid concentration in bubbling fluidized bed obtained from the base case and the correlation of Werther [37].Werther’s [37] correlations better predicted the solid concentration in the lower bed and agree well with the experimental data compared to the base case.The relative error is calculated between experimental and simulated results and with base case is 17% while with the modified version is 10%approximately.

Fig.18.Time-averaged axial solids concentration profile for bubbling fluidized bed of Dubrawski et al. [44] predicted from EMMS/bubbling drag model coupled with Werther’s bubble diameter correlation [37].

5.3.Turbulent fluidized bed

This section used only those bubble diameter correlations with the EMMS/bubbling drag model that showed the lower heterogeneity index than the base case.Fig.19 shows the comparison of heterogeneity index (Hd) obtained from the EMMS/bubbling model [28] coupled with bubble diameter correlation of Horio and Nonaka [11] (base case),withHdobtained from the bubble diameter correlations of Yacono [39],Whitehead and Young [36],and Werther [37].These above-mentioned correlations showed the lowerHdas compared to the base case.

5.3.1.Instantaneous voidage distributions

Fig.19.Comparison of heterogeneity indexes functions with respect to voidage using the EMMS/bubbling drag model coupled with bubble diameter correlations that showed lower Hd than the base case.The model input parameters are based on the experiment of Venderbosch [45] (Table 3).

The instantaneous voidage distributions of the turbulent fluidized bed predicted from the base case and three different bubble diameter correlations(i.e.,Yacono[39],Whitehead and Young[36],and Werther [37]) are shown in Fig.20.Similar to the base case,voidage distributions predicted from different bubble diameter correlations can reveal the turbulent fluidized bed’s essential characteristics,i.e.,low voidage with the presence of clusters in the dense bottom part and high voidage with no particles at the top of the bed.The void area occupies more space in the bed and the gas becomes the continuous phase,suspending dense particle clusters.The major difference in the voidage distribution contour plots predicted from different bubble diameter correlations can be observed from their bed expansion heights.The correlation of Yacono [39] predicted bed expansion height that was closer to the base case,while the distributions predicted from thecorrelations of Whitehead and Young [36] and Werther [37] revealed lower bed expansion heights.The reason for such behavior is that the heterogeneity index of Werther [37] and Whitehead-Young[36] correlations were much lower than that of the base.

Fig.20.Instantaneous voidage distribution for turbulent fluidized of Venderbosch[45] using the EMMS/bubbling drag model coupled with bubble diameter correlations that showed lower Hd than the base case.

5.3.2.Solid concentration distributions

The instantaneous and time-averaged snapshots of solid concentration distributions predicted from the base case and three different bubble diameter correlations are shown in Fig.21.The distribution of solid particles predicted from different bubble diameter correlations in the axial direction showeda dense macroscopic region at the bottom of the bed and a dilute region at the top.In the radial direction,the distributions predicted from different bubble diameter correlations revealed microscopic characteristics of the turbulent bed,i.e.,the dense region with a high concentration of solid particles near the walls and a dilute region low concentration in the core,making a famous core-annulus structure.Although our results predicted from different bubble diameter correlations show the characteristics of the turbulent bed,however,some differences can still be seen from the solid concentration snapshots.With Yacono’s [39] correlations,the solid concentration distributions showed a larger bed expansion height with some solid particle entrainment in the expanded bed section(also noticed in the base case).

Fig.21.Instantaneous and time-averaged distribution of solid concentration of the turbulent fluidized bed of Venderbosch [45] by using the EMMS/bubbling drag model coupled with bubble diameter correlations that showed lower Hd than the base case.

Meanwhile,the distributions predicted from correlations of Whitehead and Young [36] and Werther [37] showed lower bed expansion height with a clear separation of the dense bottom and diluted top regions in the main fluidization column.The solid concentration distributions predicted in the present work from the modified EMMS/Whitehead and Young [36] and EMMS/Werther[37] are more realistic than those obtained from the base case and EMMS/Yacono [39].This is again attributed to the lower heterogeneity index predicted by these two correlations than the base case.

5.3.3.Time-averaged axial solid concentration profiles

This sectionstudied the time-averaged axial solid concentration profiles obtained from the base case and with different bubble diameter correlations and compared them with the experimental measurements.The time-averaged axial solid concentration profiles predicted from the different bubble diameter correlations are shown in Fig.22.The correlations of Yacono [39] revealed the axial solid concentration profile close to the base case,which wasunder-predicted in the dense bottom section and overpredicted in the middle transition section and diluted top region.In comparison to the base case,the axial profiles predicted from the correlations of Whitehead and Young [36] and Werther [37]showed significantimprovement in the dense,dilute,and transitional regions.From the axial solid concentration profiles,we can conclude that correlation of Whitehead and Young [36] and Werther[37]are more suitable than the other tested bubble diameter correlations for the applications of turbulent fluidized beds.The relative error is calculated between experimental and simulated results and with base case is 15% while with the EMMS/Werther is approximately 7%.

Fig.22.Time-averaged axial solids concentration profile for turbulent fluidized of Venderbosch [45] predicted from EMMS/bubbling drag model coupled with bubble diameter correlations that showed lower Hd than the base case.

6.Summary

A modification in the EMMS/bubbling drag model is proposed and investigated thoroughly with parameter analysis and model validation.A parametric analysis of the original EMMS/bubbling model emphasizing the bubble-diameter correlation was first conducted.The original EMMS/bubbling model was coupled with different bubble diameter correlations,and their effects on heterogeneity indexHdwere then studied.Werther’s correlation was found to be the most suitable one to predict the heterogeneous nature of the gas-solid flow for bubbling and turbulent fluidized beds.Then some structural parameters for the modified EMMS/bubbling model (EMMS/Werther),including dense phase voidage εgc,the volume fraction of dense phasef,gas velocities in dense phase (Ugc) and dilute phase (Ugf),and energy consumption ratio(Nst/NT),etc.were obtained to analyze their variation with the mean bed voidage.Then this modified drag model combined with the TFM was examined through simulating two bubbling and one turbulent fluidized beds.Compared to the previous EMMS/bubbling model,the modified EMMS/bubbling drag model shows better qualitative and quantitative agreement with the experimental data.

Nomenclature

aacceleration term,m·s-2

Cdeffective drag coefficient for a particle

Dtbed diameter,m

dbbubble diameter,m

dsparticle diameter,m

Fdragforce,N

fvolume fraction of dense phase

ggravitational acceleration,m·s-2

Hdheterogeneity index

Nstmass-specific energy consumption for suspending and transporting particles,W·kg-1

NTtotal mass-specific energy consumption,W·kg-1

ReReynolds number

Usuperficial velocity,m·s-1

Ufccomplete fluidization velocity,m·s-1

Umfminimum fluidization velocity,m·s-1

Uslipsuperficial slip velocity,m·s-1

ureal velocity,m·s-1

β drag coefficient,kg·m-3·s-1

ε volume fraction

μ viscosity,Pa·s

ρ density,kg·m-3

Subscripts

b bubble

c dense phase

f dilute phase

g gasphase

i interphase

s solid phase

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This work is financially supported by the National Natural Science Foundation of China (21978227).

Supplementary material

Supplementary data to this article can be found online at https://doi.org/10.1016/j.cjche.2021.10.006.