DPM Simulations of A-Type FCC Particles’ Fast Fluidization by Use of Structure-Dependent Nonlinear Drag Force

: Nonlinear drag force has been a research frontier in complex gas-solid systems. The literature has reported that the commonly-used drag correlations often overestimate drag force and, thus, cause unrealistic homogeneous ﬂow structures in gas-solid ﬂuidized beds of ﬁne particles. For solving this problem, the structure-dependent drag model, derived from energy-minimization multi-scale approach, is used in discrete simulations of ﬂuid catalytic cracking particles in a small riser. The gas phase is dealt with by computational ﬂuid dynamics. Particles are considered as a discrete phase and described by Newton’s second law of motion. Gas-particle phases are coupled according to Newton’s third law of motion. Simulations show that use of structure-dependent drag model results in drag reduction, the effect of which is not so apparent as that in simulations of the two ﬂuid model. The particle clustering tendency, however, is more distinct and leads to more heterogeneous ﬂow structures in riser ﬂow with a much greater amplitude of outlet solid ﬂux ﬂuctuations. Moreover, the behaviors of particle and gas back-mixing can be captured in the present simulations, which was supported by past simulations and experimental data. The simulation time resolution is discussed. The spring constant can be artiﬁcially brought down for safe setting of larger time step when modelling the collision process between ﬁne particles with a higher calculation load. To appropriately mimic the continuous decay of van der Waals force may, however, need a much smaller time step. There is also an obvious effect of space resolution on simulations. When using a grid size smaller than 3 times the particle diameter, the simulated clusters turn extraordinarily large, and the effect of gas-solid back-mixing turns insigniﬁcant.


Introduction
Nowadays, more and more industrial processes are operated in gas-solid fluidized beds, which are known to be complex systems and feature remarkable nonlinear and multi-scale flow structures [1]. Insight into these structures has an inevitable relation with the reliable design of fluidized bed reactors. A lot of experimental investigations have been carried out concerning complex phenomena, such as particle entertainment, clustering, and bubble eruption, which are still far from complete. Correlations obtained from experimental investigations and used for reactor designs are now largely uncertain. Therefore, much attention should be paid to numerical simulation, which is one of the most useful techniques for exploring the underlying dynamics in fluidized beds. The widely-used mathematical models are the two fluid model (TFM) [2] and the discrete particle model (DPM) [3][4][5] in which particles are treated as continuous and discrete phase, respectively. How to appropriately calculate gas-solid interactions, mainly drag force, is a critical problem in both models.
So far, drag force on a single particle in a boundless flow field can be better calculated by the standard drag coefficient with an error less than 5%. Drag on the particle group in a homogeneous flow field can also be better correlated by research, mostly on the basis of experimental data. Once, the drag correlations of 11 different drag models were reviewed by Esmaili and Mahinpey [6] using TFM to simulate the momentum transfer between gas-solid phases. After investigating the effect of using different drag models in DPM simulations of a fluidized bed, Li and Kuipers [7] concluded that the selected different drag models possessed a similar predictive capability, although their accuracy might differ. Link et al. [8] found that the most frequently-used drag model, say the combination of the Wen and Yu equation [9] and Ergun equation [10], produced unsatisfactory results for a distinctly heterogeneous flow in spout fluidized beds with high-velocity jets. They suggested that the remarkable heterogeneous structures be resolved from detailed lattice Boltzmann simulations. Beetstra et al. [11] also found that use of a polydispersity factor on the basis of lattice Boltzmann simulations gives better agreement with the experimental data of segregating results in their simulations of a segregating fluidized system. Helland et al. [12] adopted a first-order approximation of the porosity function in the drag correction factor to construct two drag correlations as a combined model with a frequently cited model. Their combined drag-reduced model predicted a much lower circulation rate than the classical drag model. They suggested that the more appropriate form of the combined drag law needed to be correlated with results from experiments and direct numerical simulations.
The complexity of calculating drag force lies mainly in the indeterminate effect of the heterogeneous structure on the drag coefficient. It is believed that studies on the heterogeneous structure through experiments and numerical simulations have great significance on the development of structure-dependent drag models. One research study about the spatial scale characteristics of particle clusters was carried out by Tsuji et al. [13]. Using the data from large-scale three-dimensional Lagrangian simulations, they investigated the multi-scale structure by directly performing Fourier transforms of spatial particle concentration distributions. The so-called meso-scale structure that drag force strongly depends on is quite complex because of the high degree of freedom. The behavior of the cluster cannot be completely understood so far. The energy-minimization multi-scale (EMMS) approach, proposed by Li and Kwauk [14], can quantitatively described, in a macro-scale flow field with a solid concentration, and the dynamic balance between micro-scale particles, meso-scale dense clusters and fluid. This approach can also characterize the tendency for particles to aggregate to form clusters and for fluid to pass around clusters. Yang et al. [15], Wang et al. [16] and Wang et al. [17] successively developed the original stable EMMS approach, employed this approach to calculate the drag coefficient according to EMMS structure parameters and, thus, incorporated it into TFM simulations. The advancement of an EMMS-based drag model has raised the research frontier in TFM simulations of fast fluidization.
Only rarely did researchers pay attention to drag computation in DPM with consideration of the heterogeneous structure. Noting that the velocity and the position of each particle are known in each time step, Xu et al. [18] derived a local particle configuration and introduced multi-scale properties into DPM. Their results showed that consideration of the heterogeneous structure might greatly reduce total drag force. In order to calculate Reynolds number of the clustering particle phase, Zou et al. [19] developed a cluster-based drag coefficient model using a hydrodynamic equivalent cluster diameter. Compared with the drag coefficient model proposed by Wen and Yu, the results predicted by the cluster-based drag coefficient model were in good agreement with experimental results, indicating that their drag model was suitable to describe various statuses in fast fluidized beds. As for the simulation of Type-A particle fluidized systems, the interparticle forces of adhesive effect cannot be disregarded. When interparticle distance decreases, both the possible contact force and the van der Waals force increase. Then, the accuracy estimation of these forces requires a reduction of the time step with a consequent increase in CPU time. Therefore, there are only a few studies (see Yu and Xu [20], Ye et al. [21] and Kobayashi et al. [22]) on real Group A particles using DPM with consideration of the van der Waals force because of the large calculation time and memory.
Note that those current models derived from homogeneous systems are sometimes quite questionable when modelling circulating fluidized beds, especially those of fine and light particles, which is also one of the most important theory starting points for the advancement of the EMMS-based drag model. Yang's simulations of fluid catalytic cracking (FCC) particles in a riser [15] demonstrated that the EMMS approach could resolve the heterogeneous structure and describe the dependence of the drag coefficient on the structure parameters. It was suggested by Yang et al. that the feasibility of this approach could be used as a sub-grid closure law for the drag coefficient.
This article uses Yang's drag model [15], which was derived and adapted from the structure parameters of the EMMS approach, to calculate the local drag coefficient of single particle in heterogeneous circumstances. Simulations employing this structure-dependent model are carried out on the FCC-air system in a small riser of fine particles. Large clusters, gas-solid back-mixing and a little lower, but much violently fluctuating, solid outlet flux are captured in the simulations, which is different from the common DPM simulations.

Drag Model
The original EMMS approach proposed by Li and Kwauk [14] was for describing the flow structure of a global fluidized bed system. Yang et al. extended its principle to the control volume of TFM simulations. According to the operating conditions, the global porosity ε, the gas-solid material properties, all the eight structure parameters and the drag coefficient β can be obtained by solving the revised EMMS model (see Yang et al. [15] for details). In order to facilitate the application of the structure-dependent drag coefficient to TFM, they introduced the so-called drag correction factor as where d p is the particle diameter, ρ g is the gas density, U s is the superficial slip velocity and C d is the standard drag coefficient. It is assumed that the functions of correction factor vs. porosity correlated from the whole bed can be approximately extended to each local control volume. When grid porosity ε g > 0.74, the correction factor was mildly fitted according to the data points from EMMS solutions as follows This study adopts Yang et al.'s drag coefficient and extended drag correction factor to each particle in DPM. Then, the drag on particle i is calculated as where V p is the particle volume, ε i is the local porosity, u i is the local gas velocity, v i is the particle velocity and the local drag coefficient β i is calculated as where C di is the standard drag coefficient particle i and µ g is the gas viscosity.
This model of calculating the drag force is called Model B in the present simulations. Besides, the Wen and Yu model combined with the Ergun drag model for calculating the local drag coefficient is also used as, which is called Model A in our simulations.

Equation of Gas Motion
The Navier-Stokes equations for gas flow over grid k are given in Equations (6) and (7).
where p is the gas pressure, u is the gas velocity and τ g is the viscous stress tensor. The source term S p is the volume-averaged drag force calculated as where N k is the number of particles overlapped with grid k, F d i is the drag force on particle i, A is the area of particle disk, A i is the promotional disk area of particle i overlapped with grid k and V is the quasi-three-dimensional volume of the grid with a suppositional thickness of d p .
The transformation of the two-dimensional porosity ε 2D of a grid to the threedimensional porosity ε 3D is adopted, which was proposed by Hooman et al. [4] as the following.
The Navier-Stokes equations are discretized by the finite volume method based on the collocated grid, which makes complex problems easy to solve, compared with a staggered grid. The uniform velocity inlet, pressure outlet and impermeable wall are specified as the boundary conditions. The Simi-Implicit Method for Pressure-Linked Equation Revised (SIMPLER [23]) is used as the solver because generally it features less of a total computational load than the widely-used Simi-Implicit Method for Pressure-Linked Equation (SIMPLE [23]).

Equation of Particle Motion
The collisions of particles are dealt with according to the well-known soft-sphere model. So far as two-dimensional simulation is concerned, particles are disks. Then, the transitional motion of each particle i is computed according to Newton's second law of motion as the following.
where F vi is the interparticle van der Waals force, F ci is the collision force, p i is the local gas pressure and the rag force F di is computed by Model A and Model B given in Section 2.1.
The rotational motion of each particle i is computed as (11) where ω i is the particle angular velocity, I is the inertia movement of the particle and T ci is the torque due to collision. In Equation (10), the van der Waals force on particle i due to particle j is calculated as where H a is the Hamaker constant, e ij is the unit vector from particle i to particle j and d ij is the distance between particle i and particle j. The van der Waals force on particle i due to the bed wall is calculated as where e iw is the unit vector from particle i to the wall and d iw is the distance between particle i and the wall.

Simulations and Discussion
Two-dimensional DPM simulations were carried out on a miniature rectangle container with a side length of 5 mm and mainly on a miniature riser with the geometry of W × H = 2.5 mm × 40 mm, both of fine FCC particles. The former uses 10 × 10 grid number and the latter uses 10 × 40, 10 × 160 and 20 × 320 grid numbers. The last two sets of grids basically correspond to 4.6 and 2.3 times the particle diameter, respectively. Two drag models are employed for simulations of gas-solid flow in both containers, say drag Model A and drag Model B. The cut-off distance H 0 , i.e., the minimal value of d iw − d p and d ij − d p in Equations (12) and (13) for the particle-wall and particle-particle surface distance, respectively, is set equal to 0.4 nm in the present work. This critical value is the commonly-used one and is to avoid the divergent of van der Waals force when body surfaces approach infinitely. The Hamaker constant H a should be determined according to many parameters of physical and chemical properties whose accuracy values are difficult to be measured. Note that real particles absolutely have microscopic surface roughness, which can greatly reduce the van der Waals force. This work selects a fine-tuned and relatively low value of H a as 1 × 10 −21 N·m for Type-A FCC particles. In spite that, the calculated van der Waals force may be largely lower than the ideal value, normally it is still one order higher than the particle weight. Other fixed parameters are given in Table 1 where the operation condition corresponds to that in [15]. The inlet gas velocity is tuned to be 1.7 m/s and the calculated mean porosity at the bottom is about 0.9, corresponding to the superficial gas velocity of 1.52 m/s in [15]. Table 1. Fixed parameters for particle and gas.

Drag Reduced
In order to establish a preliminary and rough understanding on the effect of the structure-dependent drag force on the fluidization, firstly, 2000 particles in the micro rectangle container are tentatively simulated. The number of particles in and out is always balanced in the present simulations. The particles newly fed into the container do not overlap with other particles. When a particle flows out, the horizontal ordinate of the particle entering the inlet is set to the same as that at the outlet. If the particle partially overlaps with other particles, the vertical ordinate will be randomly corrected over and over. This cycle continues until the new particles do not overlap with other particles. Figure 1 shows the particle distributions simulated by the two drag models. Both models can simulate the dilute-core/dense-wall structure. There is a stronger particle aggregation effect near the wall and a more significant cavitation phenomenon in the central region in the simulations employing drag Model B. The outlet flux of solid phase G s simulated by drag Model A and Model B is shown in Figure 2. When the globally stable structure is formed, the outlet solid flux in Figure 2a fluctuates around the average value of about 100 kg/m 2 s, and that in Figure 2b fluctuates around the average value of about 70 kg/m 2 s. Although the latter is about 30% lower than the former, this difference is not so great as that reported by Yang et al., who used drag Model A and B in TFM simulations [15]. In addition, the fluctuation range of the outlet solid flux in Figure 2b is much larger than that in Figure 2a, which indicates that drag Model B predicts a much more heterogeneous flow structure and reflects the nature of dense particle clusters.   Figure 3 shows the simulated particle distributions in the riser container. In total, 8230 particles are set as randomly homogeneous in the riser at the initial moment for all the four cases. As is shown, the Wen and Yu equation combined with the Ergun equation, i.e., drag Model A, predicts globally rather than homogeneous flow structures, despite that the dilute bottom part may be affected by the inlet boundary condition, which is not entirely realistic. However, the other three groups of frames captured by the structuredependent drag force, i.e., drag Model B, show significantly heterogeneous structures with large clusters and frequent formation and breakup, and thus, an obvious evidence of fast fluidization. This indicates that the structure-dependent drag model possesses the better capability of characterizing the tendency for particles to aggregate and form clusters. Besides, it is found that the simulated clusters turn much larger with the refinement of the grid, as shown in Figure 3d. However, in Figure 3b, although the clusters form close to the wall, the sight of those extremely large clusters can hardly be observed by use of this set of the coarse grid. This is further discussed in Section 3.4.  Figure 3 shows the simulated particle distributions in the riser container. In total, 8230 particles are set as randomly homogeneous in the riser at the initial moment for all the four cases. As is shown, the Wen and Yu equation combined with the Ergun equation, i.e., drag Model A, predicts globally rather than homogeneous flow structures, despite that the dilute bottom part may be affected by the inlet boundary condition, which is not entirely realistic. However, the other three groups of frames captured by the structuredependent drag force, i.e., drag Model B, show significantly heterogeneous structures with large clusters and frequent formation and breakup, and thus, an obvious evidence of fast fluidization. This indicates that the structure-dependent drag model possesses the better capability of characterizing the tendency for particles to aggregate and form clusters. Besides, it is found that the simulated clusters turn much larger with the refinement  Figure 4 shows the simulated particles without excessive overlapping. Figure 4a shows that particles aggregate to form extremely dense clusters along the riser wall of the sectional drawing area. Figure 4b shows that high-speed particles race towards the outlet.  Figure 4 shows the simulated particles without excessive overlapping. Figure 4a shows that particles aggregate to form extremely dense clusters along the riser wall of the sectional drawing area. Figure 4b shows that high-speed particles race towards the outlet. Both regions are the most identifiable trouble spots where particles are easy to excessively overlap. Excessive overlapping in the dense region will cause nonphysical porosity lower than the minimal porosity, while excessive overlapping on the high-speed occasion will cause an unrealistic, extraordinarily higher velocity of particles. Both occasions may further cause the divergence of simulations. As can be seen from Figure 4, the simulated particles can approach close, but also do not excessively overlap. The DPM time step is determined by ∆t p ≤ π/5 m p /κ, corresponding to 7.8 × 10 −7 s for the selected material with a very small diameter of 54 µm. In this work, ∆t p = 2 × 10 −7 s, which is enough to guarantee the good performance in the simulation of particle collision. Here, the spring constant of particle κ = 50 N·m −1 , which is artificially brought down for the safe setting of relatively large time step to reduce the calculation load.   Figure 5 shows the sectional van der Waals force distribution in the simulations employing drag model B. The van der Waals force on most particles reaches extremes of either maximization or zero, and the intermediate stress values are taken just for only a small quantity of particles. The former extreme indicates the simultaneous contact between pairs of particles, while the latter indicates the surface distance is much bigger than the cut-off distance as if the action is of no force. Since the van der Waals forces are shortrange forces [21], it decays violently when the body surface distance is larger than the cutoff distance. Thus, to appropriately mimic the contentious decay process of particles, it needs very small time step, maybe smaller than the contact time step. The presentlysimulated van der Waals force distribution supports the often-used assumption that the adhesion van der Waals force could approximately be set as constant times the particle weight. For example, it was assumed that for Type-A particles, the adhesion force was constant and set to 40 times the particle weight in [22]. The presently-calculated maximal van der Waals force is 2.8 × 10 −8 N, which is about 37 times higher than the particle weight 7.5 × 10 −10 N. Sometimes there is recognizable multi-particles attraction, which can be noticed in Figure 5. Besides, the existence of recognizable adhesion stress in the dense region takes on a relative spatial continuity. The van der Waals stress on many particles also reveals the relative time continuity for particles not only in the dense phase, but also in the dilute phase, as shown in Figure 5b.  Figure 5 shows the sectional van der Waals force distribution in the simulations employing drag model B. The van der Waals force on most particles reaches extremes of either maximization or zero, and the intermediate stress values are taken just for only a small quantity of particles. The former extreme indicates the simultaneous contact between pairs of particles, while the latter indicates the surface distance is much bigger than the cut-off distance as if the action is of no force. Since the van der Waals forces are short-range forces [21], it decays violently when the body surface distance is larger than the cut-off distance. Thus, to appropriately mimic the contentious decay process of particles, it needs very small time step, maybe smaller than the contact time step. The presentlysimulated van der Waals force distribution supports the often-used assumption that the adhesion van der Waals force could approximately be set as constant times the particle weight. For example, it was assumed that for Type-A particles, the adhesion force was constant and set to 40 times the particle weight in [22]. The presently-calculated maximal van der Waals force is 2.8 × 10 −8 N, which is about 37 times higher than the particle weight 7.5 × 10 −10 N. Sometimes there is recognizable multi-particles attraction, which can be noticed in Figure 5. Besides, the existence of recognizable adhesion stress in the dense region takes on a relative spatial continuity. The van der Waals stress on many particles also reveals the relative time continuity for particles not only in the dense phase, but also in the dilute phase, as shown in Figure 5b. Beside the time resolution, the selected space resolution can also affect the simulation results. Firstly, there is a distinct effect on the simulated cluster size. That is, the smaller the selected grid size, the larger the simulated clusters, which is also shown in Figure 3 and has been discussed in the foregoing sections. Fine grid simulation can raise the resolution and provide detailed information about the fluid flow. This is not to say that the fine grid simulations are absolutely more realistic in the present study. We really find that use of fine grid smaller than 3 times the particle diameter predicts insignificant particle back-mixing. In the following section, we will expand further to discuss about this point. Figure 6 shows the distributions of particle velocity and local porosity simulated by Model B and by use of different grid resolutions. It can be noticed in Figure 6a that particles can both travel downwards and wander in a very low velocity along the wall, while in Figure 6b, there is no so distinct evidence of particle downflow. Bi et al. [24] pointed out that particle back-mixing played an important role in fast fluidization. Jin et al. [25] regarded particle back-mixing as a necessary condition for the onset of fast fluidization. In a narrow sense, back-mixing refers to the mixing of materials in a continuous process caused by movement in the opposite direction of the mainstream. The existence of such mixing influences the concentration distributions along the main flow direction. In a broad sense, near-wall particles' low-speed movement of continuous competition against Beside the time resolution, the selected space resolution can also affect the simulation results. Firstly, there is a distinct effect on the simulated cluster size. That is, the smaller the selected grid size, the larger the simulated clusters, which is also shown in Figure 3 and has been discussed in the foregoing sections. Fine grid simulation can raise the resolution and provide detailed information about the fluid flow. This is not to say that the fine grid simulations are absolutely more realistic in the present study. We really find that use of fine grid smaller than 3 times the particle diameter predicts insignificant particle back-mixing. In the following section, we will expand further to discuss about this point. Figure 6 shows the distributions of particle velocity and local porosity simulated by Model B and by use of different grid resolutions. It can be noticed in Figure 6a that particles can both travel downwards and wander in a very low velocity along the wall, while in Figure 6b, there is no so distinct evidence of particle downflow. Bi et al. [24] pointed out that particle back-mixing played an important role in fast fluidization. Jin et al. [25] regarded particle back-mixing as a necessary condition for the onset of fast fluidization. In a narrow sense, back-mixing refers to the mixing of materials in a continuous process caused by movement in the opposite direction of the mainstream. The existence of such mixing influences the concentration distributions along the main flow direction. In a broad sense, near-wall particles' low-speed movement of continuous competition against the gas flow is considered as back-mixing, which can also cause the distribution of residence time along the riser. The simulated phenomena of particle downflow and near-wall wandering are mainly caused by heterogeneous gas-solid flow structures. For example, particles aggregate to form clusters and gas flows around dense regions. Therefore, the use of a structure-dependent drag force can capture the phenomena of solid back-mixing. As should be pointed out, coarse grid simulations seem to predict a trend towards backmixing that is more pronounced, as shown in Figure 6. However, fine grid simulations show that near-wall clusters lose particles to breakup in a lower velocity; meanwhile, they further capture dispersed particles, and thus turn larger and larger. Perhaps this inconspicuous back-mixing of a narrow sense, shown in Figure 6b, explains why extremely large clusters prevail in the simulations employing fine grid simulations and drag Model B, as shown in Figure 3d. the gas flow is considered as back-mixing, which can also cause the distribution of residence time along the riser. The simulated phenomena of particle downflow and near-wall wandering are mainly caused by heterogeneous gas-solid flow structures. For example, particles aggregate to form clusters and gas flows around dense regions. Therefore, the use of a structure-dependent drag force can capture the phenomena of solid back-mixing.

Back-Mixing
As should be pointed out, coarse grid simulations seem to predict a trend towards backmixing that is more pronounced, as shown in Figure 6. However, fine grid simulations show that near-wall clusters lose particles to breakup in a lower velocity; meanwhile, they further capture dispersed particles, and thus turn larger and larger. Perhaps this inconspicuous back-mixing of a narrow sense, shown in Figure 6b, explains why extremely large clusters prevail in the simulations employing fine grid simulations and drag Model B, as shown in Figure 3d. On the other hand, it is important to find out what causes fine grid simulations to predict the inconspicuous particle back-mixing of a narrow sense. Figure 6 also shows the particle distribution as a function of local porosity. As is noticed in the figure, the use of a fine grid smaller than three times the particle diameter turns out to underestimate the local porosity, especially in the dense regions where particles aggregate and form clusters or strands. Note that porosities, whether grid porosity or local porosity, are parameters sensitive to simulations. The underestimated local porosity can cause a higher drag force On the other hand, it is important to find out what causes fine grid simulations to predict the inconspicuous particle back-mixing of a narrow sense. Figure 6 also shows the particle distribution as a function of local porosity. As is noticed in the figure, the use of a fine grid smaller than three times the particle diameter turns out to underestimate the local porosity, especially in the dense regions where particles aggregate and form clusters or strands. Note that porosities, whether grid porosity or local porosity, are parameters sensitive to simulations. The underestimated local porosity can cause a higher drag force on near-wall particles and dense clusters. This drag increase effect was also reported by Helland et al. [12], who attributed this effect inside dense clusters (with solid concentration higher than 10%) to the increased gas shear stress. Whatever the reason for the drag increase, these particles in the dense phase seem to always be suspended or even accelerated upwards. It is suggested that grid porosity should be more accurately calculated and local porosity should be calculated in more appropriate scale and by more reasonable means. In this work, two-dimensional grid porosity is precisely calculated as in one of our recent research studies [26]. However, local porosity is calculated by the area weighted average according to bilinear interpolation and is actually calculated in the fine grid scale. The more appropriate method should consider the surrounding heterogeneous particle locations and calculate in a size far larger than the fine grid [27]. That is to say, local porosity should be calculated as completely dependent on the surrounding circumstances of the object particle and in a scale that is large enough. Figure 7 shows the gas velocity distributions simulated by use of Model A and Model B. It is noticed that, in Figure 7a, gas tends to flow uniformly upwards with a slightly lower velocity close to the wall, while in Figure 7b,c, the gas flows heterogeneously and even back flows near the wall or around the clusters. This indicates that gas back-mixing can also be captured in the simulations employing the structure-dependent drag Model B. Some researchers [28,29] argued that, essentially, a plug flow existed and, thus, backmixing was impossible to be present in fast beds, compared with low-speed fluidized beds. Yang et al. [15] reported that, in their simulations of FCC-air riser flow, the solid and gas velocities in the core region were upward and much higher than those in the annulus region, while the solid and gas velocities near the wall were downward. This behavior of gas back-mixing in fast fluidization is also captured in our simulations, which occurs in fast fluidization on the occasion of high levels of solid mixing and has been found in the experimental research of Weistein et al. [30] and Li et al. [31]. Since gas back-mixing is mostly caused by particle back-mixing, here it is also noticed in Figure 6 that the fine grid simulations, corresponding to the lower degree of particle back-mixing, predict the lower degree of gas back-mixing. Figure 8 shows the variations of outlet solid flux with time. It is noticed that for both drag Model A and Model B, the simulated outlet solid flux reaches a globally steady state after 0.15 s. The average value for drag model A is about 70 kg/m 2 s, while that for drag Model B is about 60 kg/m 2 s, without such sharp wave troughs as wave crests. There is also evidence that the drag reduction of drag Model B is not so apparent in the present DPM simulations. Despite the slightly lower values of outlet solid flux, the instantaneous values predicted by the structure-dependent drag model fluctuate much more violently than those predicted by the Wen and Yu model combined with the Ergun model, indicating more disordered and unstable status in the riser due to its cluster nature. This also indicates the effectiveness of the structure-dependent drag model in improving the simulation accuracy. Although distinct heterogeneous and unstable local structures prevail in the riser, the simulated global gas-solid flow type is fast fluidization without any type of durative choking. The adopted numerical method is convergent and effective.  Figure 8 shows the variations of outlet solid flux with time. It is noticed that for both drag Model A and Model B, the simulated outlet solid flux reaches a globally steady state after 0.15 s. The average value for drag model A is about 70 kg/m 2 s, while that for drag Model B is about 60 kg/m 2 s, without such sharp wave troughs as wave crests. There is also evidence that the drag reduction of drag Model B is not so apparent in the present DPM simulations. Despite the slightly lower values of outlet solid flux, the instantaneous values predicted by the structure-dependent drag model fluctuate much more violently than those predicted by the Wen and Yu model combined with the Ergun model, indicating more disordered and unstable status in the riser due to its cluster nature. This also indicates the effectiveness of the structure-dependent drag model in improving the simulation accuracy. Although distinct heterogeneous and unstable local structures prevail in the riser, the simulated global gas-solid flow type is fast fluidization without any type of durative choking. The adopted numerical method is convergent and effective.

Summary and Prospect
The commonly-used drag correlations are not enough to accurately compute the ga solid interaction in the discrete simulation of gas-solid fluidized beds. For solving th problem, the structure-dependent drag model is used, which was derived from EMM

Summary and Prospect
The commonly-used drag correlations are not enough to accurately compute the gas-solid interaction in the discrete simulation of gas-solid fluidized beds. For solving this problem, the structure-dependent drag model is used, which was derived from EMMS approach by Yang et al. Simulations show that use of the structure-dependent drag model in DPM results in the drag being reduced to some extent. This drag reduce effect in DPM simulations of riser flow of real, fine FCC particles is not so apparent as that in TFM simulations. The particle clustering tendency, however, is far more distinct. Use of the structure-dependent drag model can reasonably predict more heterogeneous flow structures in riser flow with a much greater amplitude of outlet solid flux waves. Moreover, the behavior of particle and gas back-mixing can be successfully captured in the present simulations, which was reported in other simulations and supported by data presented by experimental researchers. The simulation resolution is also discussed. Although one can artificially bring down the spring constant to achieve a safe setting of larger time step, especially in high calculation load of fine particles simulation, to appropriately mimic the contentious decay process requires a very small time step, maybe smaller than the contact time step. There is also an obvious effect of space resolution on simulations. When using grid size smaller than 3 times the particle diameter, the simulated cluster turns extraordinary big and the effect of gas-solid back-mixing turns insignificant.
How to improve the accuracy of fine grid simulations is of great interest. Fine grid DPM simulations are advantageous in resolving the characteristic structures of clusters and bubbles and can also directly provide more detailed information of fluid motion. Grid porosity or local porosity are parameters that are both sensitive to simulations. The underestimated local porosity in ultra-fine grid simulations can cause a higher drag force on near-wall and clustering particles, thus leading to insignificant gas-solid back-mixing of a narrow sense. It is suggested that two-dimensional grid porosity should be precisely calculated and local porosity should be calculated completely according to the surrounding circumstances of the object particle and in a scale larger than the fine grid size.

Data Availability Statement:
The data presented in this study are available upon request from the corresponding author.

Acknowledgments:
The authors acknowledge the financial supported by the National Natural Science Foundation of China (61962051).

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
A particle disk area, m 2 A i particle disk area overlapped with grid, m 2 C d standard drag coefficient C di standard drag coefficient for particle i d i diameter of particle i, m d ij distance between particle i and j, m F force, N G s