Jiawei Liao,Litao Zhu,Zhenghong Luo*

Department of Chemical Engineering,School of Chemistry and Chemical Engineering,Shanghai Jiao Tong University,Shanghai 200240,China

Keywords:Fluidized bed reactor Bubbling gas-solid flows Mesoscale drag models Coarse-grid TFM simulations Heterogeneity analysis

ABSTRACT Mesoscale drag model is of crucial significance for the reliability and accuracy in coarse-grid Eulerian-Eulerian two-fluid model (TFM) simulations of gas-solid flow hydrodynamics in fluidized bed reactors.Although numerous mesoscale drag models have been reported in the literature,a systematic comparison of their prediction capability from the perspective of heterogeneity analysis is still lacking.In this study,in order to investigate the effect of several typical drag models on the hydrodynamic behaviors,the nonuniformity analysis and the sensitivity to material properties,extensive coarse-grid TFM simulations of a bubbling pilot-scale fluidized bed reactor are carried out.The results demonstrate that the mesoscale drag models outperform the empirical drag model in terms of nonuniformity due to the consideration of the influence of the mesoscale structures on the drag force in the bubbling region.Furthermore,the results reveal that our previously developed three-marker gradient-based drag model considering the solid concentration gradient exhibits satisfactory performance in predicting the bubbling flow hydrodynamics.Besides,the material-property-dependent drag model considering the explicit effect of material properties on drag corrections is most sensitive to the particle diameter.This work provides guideline for possible future improvements of mesoscale models to simulate gas-solid flow more accurately and universally.

1.Introduction

It has been a longstanding problem to accurately predict the gas-solid flow hydrodynamic behavior in fluidized beds due to the existence of complex heterogeneous structures over a wide range of length and time scales,such as clusters and bubbles,which significantly affect the mass and heat transfer and reactions in the region[1-6].The computational fluid dynamics(CFD)simulations can be used to explore the flow region and guide the scaleup design and optimization of reactors[7-9].In fact,reactor scaleup leads to the changes in flow structures [10].For example,researchers pointed out that the shape,diameter and rising velocity of the bubbles vary as the bed diameter[11].Moreover,in CFD simulations,the grid size has to be increased with scale-up due to the limitation of computing resources,resulting in ignoring the effect of mesoscale structures[12].Therefore,the selection of suitable numerical simulation methods is important.The common numerical simulation methods include direct numerical simulation(DNS) [13],discrete particle method (DPM) [14],and two-fluid model (TFM) [15],among which TFM is more feasible than the other methods for industrial-scale reactors due to treating two phases as interpenetrating continua [16,17].Unfortunately,the kinetic theory of granular flow (KTGF) and gas-solid interaction force model for closing the TFM equations still need further improvements.In particular,the interphase drag force model is critical for reliable modeling of industrial-scale reactors.Notably,numerous studies have shown that the interphase drag coefficient is very sensitive to the flow state,particle material properties,and inhomogeneous structure within the flow region [18-21].

The pioneering research efforts on drag force are devoted to formulating empirical drag correlations using experimental data,such as the Ergun model [22],the Wen and Yu model [23].The Gidaspow model was later proposed by combining the above two empirical models [24].Recently,the drag models correlated from the highly-resolved DNS data are gaining increasing popularities.These models are widely used in the field of fluidized beds.However,in industry-scale reactors,a very large control volume must be divided to solve governing equations of gas-solid flows by the discretization method.In recent years,many researchers have compared the results of coarse-grid simulations using the homogeneous drag model with the experimental results.Here,the homogeneous drag model denotes the classic empirical drag model or the numerically-derived drag model using the data from DNS of flow passing fixed particles.It was concluded that these homogeneous drag models overestimate the drag force in coarse-grid simulations of industrial-scale reactors [25,26].Several factors may attribute to such an overestimation: (1) These drag models are originally derived from the experiments or simulations based on nearly homogeneous systems where the flow passes the stationary particulate assembly [26,27].(2) Since particles are assumed to uniformly distribute in each computational volume within the framework of TFM,these homogeneous models do not take into account the effects of mesoscale structures and the mechanism of gas-solid interactions in each computational cell cannot be adequately captured [28-31].Many sub-grid models (SGMs) have been thus reported to solve the above issues,such as the energy minimization multi-scale (EMMS) method [12,31-33] and filtered method [21,27,34,35].

The EMMS model divides the complex gas-solid flow into the particle-rich dense phase and gas-rich dilute phase.The mesoscale heterogeneous structures are adequately solved by coupling hydrodynamic conservation equations and a stability condition[31,36].Following the EMMS method,the earlier work of Yanget al.[31] proposed the drag correction factor solely depending on the solid volume fraction and could significantly improve the performance of TFM-based macroscopic CFD simulations [37].Later study of Wanget al.[38] was extended to establish a subgrid EMMS/matrix model by the way of introducing a two-step scheme into the traditional EMMS model[32].More recent studies demonstrated that the structural effects in each grid should also be considered in calculating the mass/heat transfer and reaction rate in coarse-grid simulations [31,38-40].Moreover,the drag model has been well verified and validated in the coarse-grid simulation of CFB risers and bubbling fluidized beds [32,41,42].To sum up,most of the research on the EMMS drag model has been focused on selecting appropriate closure markers for different flow regimes to improve the accuracy.For example,Luoet al.[33] proposed a bubble-based EMMS drag model that introduced the local slip velocity as a marker on the EMMS/bubbling drag model and successfully verified the model in a bubbling and a turbulent fluidized beds.The filtered two-fluid model (fTFM) accounts for the unresolved subfilter-scale structures by performing averaging or spatial filters on the micro-scale balance equation in fine-grid simulations[43,44].Specifically,the coarse-grained closure,which is extracted from filtering the high-resolution TFM simulation resultsviacoarser filter or the Reynolds average dynamics theory model,is the key to closing the unresolved items appearing due to filtering the TFM model [35,45].Research shows that the filtered drag term is the most important constitutive term of filtered TFM equations[44,46,47].Therefore,the basic question is how to correct the drag force using the appropriate filtered quantities.There are two different ways to develop the filter drag force model.The first is to use the quantities that can be obtained directly from fTFM simulations to modify the drag force.For example,the Princeton group proposed the filtered drag force as a function of the filter size,filtered volume fraction,and filtered slip velocity,etc.[27,30,48].Although the predicted results were improved with this method in specific systems,it has been reported that the prediction is less satisfactory in other systems[44].The other method correlates the drag correction using quantities that are not directly available,such as the scalar variance of solid volume fraction,drift velocity,and the filtered kinetic energy of solid velocity fluctuations [47,49,50].

As mentioned above,there are many methods considering the influence of the sub-grid structures to correct the drag force,but most of them have been validated and applied in a lab-scale reactor with Geldart A particles [36].In fact,the sub-grid drag correction is also needed for the coarse-grid simulation of an industrial scale reactor where relatively large particles are present.Our previous studies proposed a material-property-dependent sub-grid modification model based on filtering theory [51] and an effective three-marker drag model [52],which are suitable for Geldart type A and B particles.Those two models have been well validated in the lab-scale reactors.However,the predictive ability of those models for pilot-scale bubbling beds with larger particles is still not clear.Moreover,the time average results were well evaluated by previous studies while further statistical analyses on the instantaneous results are still lack.

In the present study,the influence of various drag models on the ability to simulate the complex gas-solid flow is investigated in a pilot-scale bubble fluidized bed.Specifically,the materialproperty-dependent sub-grid modification model and effective three-marker drag model,which were developed for the fast fluidized bed risers and turbulent fluidized beds in our previous work,are extended to the bubbling fluidized bed.By analyzing the instantaneous and time-average statistical results,the effect of the drag model on the nonuniformity of gas-solid flow hydrodynamics and the sensitivity to material properties are systematically studied.

2.Numerical Simulation

2.1.TFM approach

In this work,CFD is based on the Eulerian-Eulerian TFM.In this approach,conservation equations are derived for each phase from the concept of averaging [29].Meanwhile,the KTGF is used to model the interaction forces between particles.

2.1.1.Governing equations

In a bubbling fluidized bed reactor,the mass conservation equations for phasei(i=g for gas phase,s for solid phase)can be written as:

where αi,ρiand virespectively represent the volume fraction,density and velocity for gas and solid phase.In a two-phase system,αg+αs=1.

The momentum conservation equation for gas phase and solid phase derived from Navier-Stokes equations can be written as:

wherepandpsare the gas phase static pressure and solid phase pressure.gis the gravitational acceleration.β represents the drag coefficient.τgand τsrespectively denote shear stresses for gas and solid phase.

2.1.2.The kinetic theory of granular flow

The stress tensor can be derived from the Newtonian stressstrain relations due to the gas phase and the solid phase was treated as an incompressible fluid.

where I is unit stress tensor,μiand λidenote shear viscosity and bulk viscosity for gas phase and solid phase,respectively.Generally,μgis a constant or as a function of temperature,and λgis assumed to be zero[53].Furthermore,KTGF,which is based on the similarity between the flow of gas molecules and particles,is used to the closure of the solid stress term.Eventually,the solid pressure,solid shear viscosity,and solid bulk viscosity are expressed as the functions of granular temperature（Θ）[17].The variable represents the fluctuation of the particle velocity and can be defined as:

Moreover,the granular temperature can be calculated by solving the granular energy equation [55].

where the terms on the right side respectively represent energy generated by solid stress tensor and diffusion of energy,γsis denotes the collisional dissipation of energy and φgsis the energy exchange between gas and solid phases [56].The relevant closure laws are listed in Table 1.

Table 1 Closure laws used for KTGF

2.1.3.Drag models

In the momentum conservation equation for gas and solid phases,the termsdenote the interphase momentum exchange,in which only the drag force is considered since it’s the most important forces in the interphase forces [59].In this work,a homogeneous drag model and three heterogeneous filtered drag models are tested and evaluated.

The homogeneous Gidaspow drag model is a combination of the Ergun model and the Wen and Yu model,which can be expressed as [24,55]:

When αs＜0.2,the drag coefficient for gas phase and solid phase(β) can be calculated as:

The heterogeneous two-marker filtered drag model (filtered drag model) proposed by Sarkar is derived from highly resolved 3D simulations of gas-particle flows.The model was validated in a bubbling fluidized bed with Geldart A particles (78.66 μm) and more details can be foundviatheir article[60].The drag coefficient can be expressed as:

where βgs,microdenotes microscopic drag coefficient,which is described by the Wen and Yu model:

Here,the filter size Δfilteris equivalent to the volume of each CFD grid (Δgrid) in this work.νtdenotes particle terminal velocity and can be calculated by the Wen and Yu [23].

The effective three-marker drag model(gradient drag model)is applied in this study,which uses the concentration gradient as an additional marker to explain the parabolic spatial concentration distribution within a computational grid [52].Our previous work demonstrated that the gradient drag model exhibits good predictive performance of flow characteristics in a 3D turbulent fluidized bed with Geldart B particles (139 μm).Meanwhile,bubbling fluidized beds and fast fluidized beds may also be within the scope of application of this model because it is derived based on theoretical analysis.The drag coefficient can be expressed as:

whereFd（Re,αg）represents uniform drag force and can be obtained by Tanget al.[61].

The correlation coefficientHd,Gradientcan be calculated as:

where Δαgdenotes the difference existing between the neighboring grids cells and can be expressed as [52]:

In which,?αgrepresents gas phase volume fraction gradient and can be obtained directly from simulation results.

The material-property-dependent subgrid drag model (MPD model)is derived from the data generated from the periodic simulation [51].The model was validatedviathe comparison between the simulated results and experimental data in five CFB risers with the particle diameter ranging from 54 to 300 μm in our previous work.The drag coefficient can be written as:

2.2.Simulation conditions

In order to verify and compare the drag models,the CFD software Fluent (version 17.2) is used to simulate the gas-solid flow in a pilot-plant ethylene polymerization FBR in this study.The FBR was reported in detail in the earlier study by Cheet al.[56].The simplified reactor geometry is represented in Fig.1.It consists of a disengagement zone,transition zone,and a fluidized zone shown in Fig.1(a) and the specific size is shown in Fig.1(b).As reported in the literature,in order to simplify the model,the inlet gas consisting of ethylene,hydrogen,nitrogen,and 1-butylene are replaced by a gas mixture.The material properties are summarized in Table 2.In the initial state,the solid particles are patched in the bed to a 1.137 m height and then are fluidized by the gas at a specific velocity entering from the bottom.The no-slip wall boundary condition is selected for the gas phase,and the free-slip wall boundary condition is selected for the solid phase.Meanwhile,the coefficient of restitution for particle-particle collisions (es) is set to 0.9.Besides,a time step of 5×10-4s is adopted.The remaining numerical settings and operating conditions are shown in Table 2.

In this work,four grid sizes are used to analyze the grid sensitivity.For comparing the drag models,four models are chosen:Gidaspow model,gradient model,MPD model,and filtered model,respectively.In addition,the particle diameter is also investigated for comparison and analysis of drag models.In summary,a total of 19 cases are designed to simulate the pilot-plant ethylene polymerization FBR,as listed in Table 3.From the curves of solid volume fraction (SVF) and pressure drop fluctuating with flow time shown in Fig.2,it can be concluded that the case in this study requires approximately 4 s to reach a pseudo steady state.Thus,the simulation results obtained from flow time of 10 to 20 s is used for all the remaining analyses.

Table 2 Model settings and parameters [56]

Table 3 Modeling scheme

Fig.1.The pilot-plant ethylene polymerization FBR [56].(a) Schematic 3D geometry.(b) Schematic details.

3.Results and Discussion

3.1.Grid sensitivity analysis and model validation

Before comparing the models,a mesh sensitivity analysis is first performed to exclude the effect of the grid size.It is worth notingthat the grid sensitivity analysis for the sub-grid drag models has been well investigated by open reports[25,51].Moreover,because the homogeneous Gidaspow drag model is very sensitive to the variations of mesh size,the following grid analysis is based on the heterogeneous gradient drag model.Four grid resolutions of Δ=80ds,61ds,45dsand 35dsare used in this study.The timeaveraged pressure drop simulated by the gradient model with various grid resolutions are shown in Fig.3.It is found that reasonable results in good agreement with experimental data can be obtained based on various grid resolutions.Meantime,the figure indicates that as the grid size decreases,the predicted results are less dependent on the grid size and the result predicted by Δ=45dsis very close to that by Δ=35ds.

Fig.2.Cross-section averaged solid volume fraction and pressure drop at the height of 1.0 m in the FBR.

Fig.3.Time-averaged pressure drop predicted by Gradient model with different grid resolutions.

To quantitatively analyze the effect of grid size on flow parameters,the grid convergence index(GCI)is introduced,which can be expressed as [52]:

where ther,p,andfm,fnrespectively represent grid refinement ratios,discretization order,and the results computed by grid resolutions with fine and coarse.Because the bed expansion height(H) and the average solids concentration () in the fluidized zone are the key parameters for bubbling fluidized bed.Therefore,we use theHandas examples to calculate the GCI.

Table 4 shows the GCI calculated from Eq.(41).The table reveals that the GCI34calculated byHandat grid resolutions of Δ=35dsand Δ=45dsis far lower than the others.The result qualitatively indicates that there is no significant difference between the simulation results using grid resolutions of Δ=35dsand Δ=45ds.Considering that the bubbling and clustering structures captured by a grid resolution of Δ=35dsare finer,this grid resolution is used in the following simulations.

Table 4 The bed expansion height and average solids volume fraction in the fluidized zone and GCI calculated by gradient model with different grid numbers

3.2.Heterogeneity analysis

3.2.1.Heterogeneous flow hydrodynamics

To compare and analyze the influence of drag models on hydrodynamic behavior,Fig.4 provides the transient and time-averaged snapshots for solid volume fraction predicted by four models.It can be clearly seen that the typical phenomenon of bubbling fluidization is evident in all simulations,but there are qualitative differences in the distribution of particles.Compared with the homogeneous Gidaspow drag model,simulations by the three heterogeneous drag models predict much denser flow regions and lower bed expansion heights due to the particle aggregation at the bottom of the bed.Moreover,that indicates that the heterogeneous drag corrections are more evident in some regions and the similar phenomenon is also reported in the literature [51,52,60].Meanwhile,the dense clustering structures predicted by three heterogeneous drag models is more significant than the homogeneous drag model,among which the filtered model is the most prominent.In fact,these mesoscale structures that continuously form and breakup are due to compromise in competition between the phases of gases and particles.Furthermore,the drag reduction is associated with the occurrence of these structures.In Fig.4,the size,number,and shape of the bubbles and clusters predicted by different drag models are different.These mesoscale structures have a profound impact on the accurate prediction of hydrodynamic behavior.Therefore,the choice of the drag model is very important.

To further analyze the differences in drag model predictions,Fig.5 shows a quantitative comparison of the time-averaged axial distributions of solids simulated by the Gidaspow model,gradient model,MPD model,and filtered model.Most of the solid particles are concentrated in the fluidized zone of the pilot reactor,in which the average solid phase concentration and bed expansion height(threshold value of solid phase volume fraction αs=0.003) predicted by the four models are shown in Table 5.It is clear that the average solid phase concentration predicted by the filtered model in the fluidized zone is much higher than the other model.This means that the coarse drag force predicted by the filtered model is smaller than that by the other models,resulting in an excessive clustering of solids in the dense fluidization zone.This result may be attributed by the fact that the data for deriving the filtered model are derived from fine grid simulations based on Geldart A particles.Although diverse variables are rendered in a dimensionless form in the construction process of the filtered model in order to make it more applicable for elaborate flow patterns,this may still be insufficient to establish a generalizable model dependent on material properties.Moreover,the gradient model and the MPD model are close and slightly larger than the Gidaspow drag model.In addition,it can be observed from Fig.5 that the solid phase concentration distributions predicted by the four models as a whole show an ‘‘S-shaped” distribution with a dense bottom and a dilute top.Moreover,the trend of solid concentrations predicted by different models in the dense region at the inlet is analogous to the average solid concentration in the vertical region mentioned above,but the prediction by the gradient model is significantly denser than that by the MPD model.This implies that the drag force predicted by the filtered model is lower than that by the other three models.Furthermore,the gradient model and the MPD model are generally similar but there are differences of solids concentrations between the dense region and the dilute region.Besides,Fig.5 displays that the average solid concentration in the fluidized zone predicted by the four models decreases with the increase of height.Especially,the filtered model has the fastest decreasing rate,followed by the gradient and MPD models.However,the average solid concentration in the fluidized zone predicted by the Gidaspow model remains relatively stable.This corresponds to the degree of uniformity predicted by the different models.

Table 5 The bed expansion height and average solids volume fraction in the fluidized zone calculated by different models

Fig.4.The transient and time-averaged (TA) snapshots for solid volume fraction (SVF): (a) gradient model,(b) Gidaspow model,(c) MPD model,(d) filtered model.

Fig.5.The time-averaged axial solid volume fraction distributions predicted by four models.

Fig.6 gives the time-averaged radial distributions of the solid volume fraction and gas-solid slip velocity simulated by the four drag models at different heights.Fig.6(a)-(c) presents the radial solid distributions at bed heights ofH=0.5,1.0,and 1.5 m.The figure indicates clearly that the typical non-uniform distribution of solid volume fractions predicted by the four models is captured,in which the solid particles tend to be more concentrated near the wall due to the influence of the wall,i.e.,the core-annular structure.Meanwhile,the solid concentration predicted by the Gidaspow drag model shows a flat-to-increase steady trend with the increasing radial position (r/R),while those by the heteroge-neous drag models tend to decrease at first and then increase,which is more pronounced atH=0.5 m.This means that the heterogeneity predicted by the heterogeneous drag models is reflected in both the axial and radial positions.It is worth noting that although the Filtered model predicts a larger mean solid concentration than the other models in the previous description,its trend in the radial spatial position is about the same as that by the gradient and MPD models.In addition,it can be seen from Fig.6(d)-(f) that the gas-solid slip velocities at different heights predicted by the models all decrease near the bounding walls except for those predicted by the gradient model.This is because the continuous aggregation of the downward particulate flow near the wall takes place due to the influence of the solid circulation mode and the wall.Since the concentration gradient at the wall region is significantly larger than that in the core region,and it is the key marker for the gradient drag model,the drag force at the wall region predicted by this model is smaller than that of other models,resulting in an increase in the slip velocity between the gas phase and solid phase.

3.2.2.Heterogeneity analysis of the gas-solid flow

The previous analysis is aimed at the comparison of the gas and solid flow behavior simulated by different models.In order to further analyze the rationality of the heterogeneous drag correction model,this section focuses on the evaluation of the flow nonuniformity predicted by the drag model.The following analysis is carried out at various heights of 0-1.8 m,where the flow highly develops,in order to exclude the influence of reactor boundary conditions on hydrodynamic analysis.

Fig.6.Time-averaged radial solid distributions and slip velocities simulated by the four drag models at different heights:(a)-(c)radial solid distributions at H=0.5,1.0,1.5 m heights,(d)-(f) radial slip velocities at H=0.5,1.0,1.5 m.

Fig.7.Probability distribution functions(PDF)of the solid phase concentrations predicted by different drag models:(a)PDF of the time-averaged solid phase concentrations,(b)-(f) PDF of the solid phase concentrations at flow time 16-20 s.

Fig.7 shows the probability distribution functions (PDF) of the solid phase concentrations predicted by different drag models.It is worth noting that Fig.7(a) is obtained from time-average data while Fig.7(b)-(f) are obtained from instantaneous data at flow time of 16-20 s.In Fig.7(a),the PDFs derived from the timeaveraged results predicted by the gradient model,MPD model,and Gidaspow model have an obvious peak,while that by the filtered model has two similar peaks.Moreover,there are obvious fluctuations near the high concentration of the solid phase in PDF curves.On the one hand,the peak predicted by the Gidaspow model is higher and narrower than the remaining models,and it occurs in the case of the lowest solid concentration.Moreover,the peak predicted by the filtered model is evidently lower and broader.This indicates that the average solid concentration predicted in the fluidized zone by the Gidaspow model is lower than those by the other models,while its maximum value is predicted by the filtered model.This observation is consistent with the previous analysis.Meantime,the narrow and high peaks indicate that more free-moving particles in the core region in the reactor are predicted by the Gidaspow model,and they are difficult to form clusters.On the other hand,by comparing the above analysis,it can be concluded that the fluctuation in the high-concentration area is related to the particle aggregation near the wall.This suggests that particle clusters tend to form in the bounding wall region,and the clustering degree predicted by the Filtered model is highest,followed by the gradient and MPD models.In Fig.7(b)-(f),it can be seen that the curves of PDF derived from instantaneous data are quite different from those derived from timeaveraging data.First,the maximum solid concentration in PDF derived from the instantaneous data is 0.63,which is clearly far greater than the maximum solid concentration derived from the time-averaged data.Second,the solid concentration probability densities derived from instantaneous data are much smaller than those derived from time-averages data.Finally,the peaks in the PDF derived from instantaneous data are wider than those derived from time-averaging data.This means that the particulate flow in the fluidization region is far more intense and complex.In this situation,the accurate conclusion cannot be solely drawn from the time-averaged results,especially for the non-uniform phenomena that exist in the flow region.For the gradient and MPD models,it can be seen from the Fig.7(b)-(f) that the probability of the high solid concentration region predicted by the gradient model is greater than that of the MPD model.However,the width and position of PDF peaks by the gradient model are similar to those by the MPD model.This indicates that the maximum solid holdup of the clusters predicted by the gradient model is larger than that of the MPD model,and more clusters with high solid concentration are predicted,but the overall heterogeneity is similar.Notably,it can be seen from Fig.7(b)-(f) that the filtered model predicts a much denser bed than the other models.This suggests that the solids in the fluidization region predicted by the filtered model are more likely to form clusters,and it reflects a higher degree of inhomogeneity.

Fig.8 shows probability distribution functions of the slip velocity between the gas and solid phases for different drag models.Similar to Fig.7,the figure also indicates that the curves of PDF derived from instantaneous data and those from the timeaveraged data are different.Fig.8 reveals that the peaks predicted by the Gidaspow model correspond to the lowest slip velocities and those by the Filtered model are the highest,which is similar to the trend of the solid concentration probability distribution function.This suggests that the inhomogeneous phenomenon in the flow region is not only reflected in the particle distribution but also affects the slip velocity between the two phases.Besides,it can be seen that the peak of PDF predicted by the gradient model is higher and narrower than that of the MPD model in the probability density plot of slip velocity.That is,the slip velocity in the fluidized zone predicted by the gradient model is more uniform than that by the MPD model.This may be because the MPD model is derived from the gas-solid slip velocity in the initial homogeneous state under different materials,it is more sensitive to the slip velocity than the three-marker gradient model.

To quantitatively compare the nonuniformity of gas-solid flows simulated by different models,the nonuniformity index (NI) proposed by Zhu and Manyele [62] is introduced.The NI is defined as the ratio of the standard deviation of the solid concentration in the flow region to its maximum standard deviation,expressed as:

where the standard deviation of the solid concentration in the flow region σ（αs）and the maximum standard deviation σmax（αs）can be calculated by.

Fig.8.PDF of the slip velocity predicted by different drag models: (a) PDF of the time-averaged slip velocity,(b)-(f) PDF of the slip velocity at flow time 16-20 s.

Fig.9.The global nonuniformity index (NI) calculated by different models at different flow times.

Fig.9 shows the NI of the four models for different flow time at a height of 0.0-1.6 m.The figure plots that the NIs calculated by different models range from 0.3 to 0.5 over time.Furthermore,it can be found that the NI manifests different trends over time.Fig.10 shows the transient snapshots for solid volume fraction at flow time of 10,18 and 20 s corresponding to the instants of maximum NI predicted by different models,respectively.The figure also shows that there are more bubble and cluster structures in the snapshot predicted by the gradient model fort=10 s corresponding to the maximum NI.Meanwhile,the bubble and cluster structures predicted by the filtered model are the most abundant at a flow time of 18 s,followed by the gradient and MPD model,while the structures predicted by the Gidaspow model are relatively unclear.This trend is consistent with that described in Fig.9.It indicates that structures such as bubbles and clusters have a profound effect on nonuniformity.Finally,the time-averaged NI and its standard deviation for different models are shown in Fig.11.The figure demonstrates that the nonuniformity predicted by the filtered model is larger than those by the other models.More specifically,that by the gradient model is slightly greater than that by the MPD model while that by the Gidaspow model is lower than those by the remaining models.In addition,the bubbling and clustering structures predicted by the gradient model are more remarkable than those by the other models.As a result,the standard deviation of NI predicted by this model is larger than that of other models.

Fig.10.The transient snapshots for solid volume fraction(SVF)predicted by different models:(a)analyzed region,(b)-(e)flow time at 10 s,(f)-(i)flow time at 18 s,(j)-(m)flow time at 20 s.

Fig.11.The time-averaged and standard-deviation of NI for the different models.

Fig.12.The transient snapshots for solid volume fraction predicted by different models: (a)-(d) particle diameter (D) 250 μm,(e)-(h) particle diameter 446 μm,(i)-(l)particle diameter 590 μm,(m)-(p) particle diameter 830 μm.

Fig.13.The global nonuniformity index (NI) calculated by different models with different particle diameters at flow time 20 s.

4.Conclusions

In this study,the CFD simulation is carried out to investigate the effect of several typical drag models on hydrodynamic behaviors,inhomogeneity of flow region,and the sensitivity to material properties.By analyzing the results of time-averaged,instantaneous,and statistical results from coarse-grid simulations,such as the spatial distribution of solids,bed expansion height,gas-solid slip velocity,averaged solids concentration and inhomogeneity in the fluidized zone,the following major conclusions can be drawn:

(1) All of the evaluated models can simulate reasonable bubbling flow phenomena such as the ‘‘S-shaped” distribution and ‘‘core-annular” structure.

(2) The gradient model,MPD model,and filtered model appear to be more reasonable than the homogeneous Gidaspow model in terms of the drag force prediction and nonuniformity due to the consideration of the influence of mesoscale structures on the drag force in the flow region.

(3) Since the gradient model uses the concentration gradient as a closure marker,the correction mainly works in the areas with large concentration gradients such as near the wall region and the reactor inlet region where the predicted clusters and bubbles are obvious.

(4) It is revealed that an increase in the particle diameter (Geldart B to Geldart D)will lead to a whole increase in nonuniformity in the bubbling fluidized bed and this phenomenon was most pronounced for the MPD model.

In summary,it is necessary to use the heterogeneous drag model for coarse-grid simulations of industrial scale reactors.However,it is also important to ensure the reasonableness of the mesoscale structure predicted by the heterogeneous drag model.In addition,the important role of the particle diameter in predictions of flow hydrodynamics is emphasized in this study.This factor may be considered for possible future improvements of the mesoscale drag model.

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

The authors thank the computing resources provided by the π 2.0 cluster supported by the Center for High Performance Computing at Shanghai Jiao Tong University.

Supplementary Material

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