CFD-DEM Simulation of Fast Fluidization of Fine Particles in a Micro Riser

: In recent years, the discrete element method (DEM) has gradually been applied to the traditional ﬂuidization simulation of ﬁne particles in a micro ﬂuidized bed (MFB). The application of DEM in the simulating fast ﬂuidization of ﬁne particles in MFB has not yet received attention. This article presents a drag model that relies on the surrounding environment of particles, namely the particle circumstance-dependent drag model or PCDD model. Fast ﬂuidization in an MFB of ﬁne particles is simulated using DEM based on the PCDD model. Simulations indicate that the local structure in an MFB exhibits particle aggregation, which is a natural property of fast ﬂuidization, forming a structure where a continuous dilute phase and dispersed concentrated phase coexist. There exists a strong effect of solid back-mixing in an MFB, leading to relatively low outlet solid ﬂux. The gas back-mixing effect is, however, not so distinct. The axial porosity shows a monotonically increasing distribution with the bed height but does not strictly follow the single exponential distribution. The solid volume fraction at the bottom of the bed is signiﬁcantly lower than the correlated value in CFB. The axial heterogeneous distribution of the cross-sectional average porosity in the lower half of the bed is also weakened. The radial porosity shows a higher distribution pattern in the central region and a lower one in the sidewall region.


Introduction
Gas-solid fluidized beds are widely used in environmental, energy, chemical and other process engineering fields [1]. Numerical simulation has become an indispensable means of studying gas-solid two-phase fluid dynamics for the cognition of a fluidized bed. The discrete element method [2][3][4][5], or DEM, providing particle-level information has been a tool needed to guide process design and operation in the field of fluidization engineering. The particle and device sizes simulated with DEM show a wide range. With the development of computer technology, high-performance algorithms and modeling capabilities, DEM has been able to simulate laboratory-scale coarse particle reactors. This method is gradually recognized by the industry and will inevitably be widely used [6]. On the other hand, with the improvement of the model and the calculation accuracy, the particle size used in DEM simulation gradually turns refined. Yu's research group once simulated and analyzed different fluidization types of particle systems with particles as small as 30 µm [7].
In order to distinguish the particles in the gas-solid system with practical significance, Geldart [8] divided the particles into four classes: A, B, C and D. This is a widely accepted particle classification method in academia. Class C particles or ultra-fine particles are usually considered difficult to be fluidized; Class D particles or ultra-coarse particles have been used in numerous DEM simulation studies but are limited to fewer operations such as spout and slugging in practice; Class B particles are fluidized materials with a relatively small size in previous DEM simulations but still classified as coarse particles; and Class A particles or fine particles have a more smaller size (30-100 µm) and density (less than 1400 kg/m 3 ). Fine particles are rarely used in DEM simulation, but they are widely used in engineering. For example, the catalytic cracking catalyst (FCC) particles, as a typical fine fluidized material, show advantages and characteristics different from other types of particles, such as a high gas-solid mass transfer rate, high bed expansion ratio and high heat transfer rate; they significantly expand before bubbling. Moreover, their cohesion is stronger, and the van der Waals force should not be ignored in many cases.
Ye et al. [9] modeled the fluidization of class A particles using DEM in 2004. The next year, Potic et al. [10] from the same study group proposed the concept of a micro fluidized bed (MFB), which is far small than an experimental circulating fluidized bed (CFB). They also studied the micro fluidized bed experimentally and theoretically. Han et al. [11] indicated that MFB has the advantages of convenience, safety, low cost and high efficiency. As an ideal isothermal reactor, it can be used as a reaction analyzer. As a small high selectivity module reactor, it can be used to produce high-value fine chemicals. Since only a few coarse particles side by side can fill the whole bed diameter, it is inevitable to abandon the use of class D particles or even class B particles. One has to choose Class A particles as fluidized materials in MFB experiments and DEM simulations. Therefore, the DEM simulation of class A particle fluidization in an MFB may become a hot research field.
Due to the lack of experimental data and simulation research, it is difficult to find reference literature. The existing DEM simulation results are mainly limited to fixed, uniform expansion, bubbling, slugging and turbulence fluidized beds [12][13][14][15][16][17]. Taking additional research in the past two years as examples, Guo et al. [18] presented a study of the solid-like and fluid-like states in the homogeneous fluidization regime of Class A particles. Li et al. [19] also performed a detailed analysis of bed hydrodynamics for polydisperse gas-solid flow of Class A particles. Both of the recent works implemented 3-D simulations. It is currently imperative to carry out the DEM simulation study of unconventional fast fluidization in MFB, even if it is a 2-D simulation study.
One difficult problem is how to reasonably simulate the flow structures, especially the macroscopic structure in an MFB. The local flow structure of particle agglomeration could be well simulated in [20]. However, the gas-solid back-mixing behavior could not be reasonably simulated, leading to the macroscopic structure seriously not insistent with the well-known core-annulus theoretical model. Can DEM reasonably simulate the fast fluidization of fine particles in MFB? Looking at the problem, this paper proposes a drag force model considering the surrounding particle environment and tries to solve the problem by reasonably calculating the drag force.

Drag Model
The drag model formally adopts the drag formula of a single particle in the particle group where d p is the particle diameter, C di is the apparent drag coefficient of the particle group around particle i, ρ g is the gas density, u i is the local gas velocity, v i is the local gas velocity and ε i is the local porosity. A gird-independent calculation method for ε i is given in the following. A neighborhood of the target particle is firstly determined to consider the influence of each surrounding particle in the neighborhood on the target particle. The influence of each surrounding particle is distinguished according to the different distances from the surrounding particle. This grid-independent calculation method is similar to the kernel approximation method in smooth particle dynamics (SPH) [21]. In this way, Xu et al. [22] calculated ε i as where N i is the number of particles in the neighborhood, h is the characteristic radius or smooth length of the neighborhood, r i − r j is the distance between particle i and j and W(d, h) is the normalized kernel function which is calculated as However, the anomalies close to the boundary have been addressed for the above Equations (2) and (3). The local porosity close to the boundary is effectively improved by adding boundary virtual particles in a previous study [23]. Figure 1 shows the diagram of the neighborhood, real particles and virtual particles in the kernel approximation method.  (2) where i N is the number of particles in the neighborhood, h is the characteristic radius or smooth length of the neighborhood, j i r r − is the distance between particle i and j and ) , ( h d W is the normalized kernel function which is calculated as However, the anomalies close to the boundary have been addressed for the above Equations (2) and (3). The local porosity close to the boundary is effectively improved by adding boundary virtual particles in a previous study [23]. Figure 1 shows the diagram of the neighborhood, real particles and virtual particles in the kernel approximation method. In ref. [23], h is set to be 2.5 d p . In fact, the accurate tune is difficult even in SPH simulations. On the other hand, in Equation (2), the contribution of each particle j for the local solid volume fraction of particle i is W( r i − r j , h) · πd p 2 /6. However, this contribution is formally not three-dimensional but two-dimensional. Moreover, the contribution should not be fixed for different particle systems. Thirdly, the basis for the correction is that the global mean values of local porosity cannot deviate from the global mean values of grid porosity. This problem is solved by introducing a solid volume fraction multiplier as In order to avoid parameter testing, λ is determined by the following equation as where ε ts is the particle filling ratio of the bed, i.e., the total solid volume fraction and N is the total number of real particles that are randomly distributed throughout the bed. Combined with the apparent drag coefficient proposed by Wen and Yu [24] and the single-particle drag coefficient proposed by Schiller and Naumann [25], the apparent drag coefficient is calculated as where ε i is calculated according to Equation (4). Thus, the particle circumstance-dependent drag (PCDD) model is constructed according to Equations (1) and (6).

Simulation Method
The Navier-Stokes equations of the gas phase motion are expressed as Equations (7) and (8), where ε g is the mean porosity, p is the pressure, u is the gas velocity, t is the time, τ g is the stress tensor and S p is the momentum exchange source term. S p is calculated as where N k is the number of particles that overlap with grid k, A is the particle disc area, A i is the overlap area of particle i with grid k and V is the grid volume, which is calculated as if the grid had a thickness of d p .
The exact area fraction model is used [26] to calculate the 2-D porosity, and the following conversion formula [3] is used to convert ε 2D to the 3-D porosity ε 3D as The finite volume method is used for discretizing the Navier-Stokes equations, using the consistent velocity inlet, pressure outlet and wall impenetrable conditions as boundary conditions and using the SIMPLER method [27] to solve the discretized algebraic equations.
The particle translation is described according to the Newton's Second Law as where F ci is the sum of contact force on particle i, F vi is the sum of van der Waals force on particle i, V i is the volume of the particle, p i is the local pressure. The soft sphere method is used to handle the particle-particle and particle-wall collision. The mass of the wall is infinite. The drag force F di is calculated using the PCDD model. The van der Waals force between particles 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 of the wall applied to particle i is calculated as where e iw is the unit vector from the particle to the wall, and d iw is the distance between the particle and the wall. Moreover, to calculate the non-contact van der Waals force according to Equations (12) and (13), one needs to set the cutoff distance H 0 , i.e., the lower bound of d ij − d p and d iw − d p . The current simulations take H 0 = 0.4 nm. The particle rotation is described as where ω i is the particle angular velocity, I is the inertia of the particle and T ci is the torque of the collision.
In the present two-dimensional simulations, the size of the micro riser is D × H = 2.5 mm × 40 mm. The number of grids discretizing the corresponding flow field is 10 × 160. The maximal particle number used is 11,452. Initially, 8230 real particles are set static and randomly distributed throughout the whole bed, which determines λ = 0.948 according to Equation (5). The particle feed way follows the import-and-export balance. Table 1 shows the other fixed parameters of the gas-solid phases in the simulations. Table 1. Fixed parameters for particle and gas.

Particle Agglomeration
The transient particle and porosity distribution are given in Figure 2. The coexisting structures of the dilute phase as disperse particles and the dense phase as granular aggregates are consistently presented in the MFB. Kuipers et al. [28] recommended using the porosity of 0.85 as a standard to distinguish the emulsion phase and the bubble phase in a bubbling fluidized bed. Figure 2 also uses this standard to distinguish the dense phase and the dilute phase. As can be seen from the porosity distribution, the dilute phase is continuous while the dense phase is dispersed. In different areas of the bed, the boundaries of the two phases are sometimes blurred and sometimes clear. With the change of time and space, clusters show the dynamic behaviors of continuous formation and fragmentation, and they vary in shape and size. porosity of 0.85 as a standard to distinguish the emulsion phase and the bubble phase in a bubbling fluidized bed. Figure 2 also uses this standard to distinguish the dense phase and the dilute phase. As can be seen from the porosity distribution, the dilute phase is continuous while the dense phase is dispersed. In different areas of the bed, the boundaries of the two phases are sometimes blurred and sometimes clear. With the change of time and space, clusters show the dynamic behaviors of continuous formation and fragmentation, and they vary in shape and size.   Figure 3 presents the instantaneous distribution of the particle velocity within the 0.01 m height range in the bed. As shown in the figure, there is a serious back-mixing effect of the particles. The large cluster forms near the wall and then moves down along the wall. The radial velocity of the particles inside the cluster is small, and the axial velocity is large. Where the cluster passes, a strong tail vortex is formed to attract the tail particles to follow down. At the same time, a part of the peripheral particles begins to break away from the cluster due to the large radial velocity. Eventually, the cluster disintegrates due to the penetration of the fluid. At the bottom of the bed, the newly entering particles have not been fully accelerated, and the upper particles are in frequent contact. This results in a low radial and axial velocity of the particles in this area and also in a disorder orientation. In addition to the concentrated phase area, the particle velocity near the wall is relatively small, while the velocity in the central area is relatively large. Moreover, for the particles in the central area, whether they are the dispersed particles or particle groups, the axial velocity is mainly upward. As also indicated in Figure 3, the positive and negative velocity might be detected in almost all radial and axial positions. been fully accelerated, and the upper particles are in frequent contact. This results in a low radial and axial velocity of the particles in this area and also in a disorder orientation. In addition to the concentrated phase area, the particle velocity near the wall is relatively small, while the velocity in the central area is relatively large. Moreover, for the particles in the central area, whether they are the dispersed particles or particle groups, the axial velocity is mainly upward. As also indicated in Figure 3, the positive and negative velocity might be detected in almost all radial and axial positions.  Figure 4 shows the instantaneous distribution of the gas velocity in the 0.01 m height range, which corresponds to time in Figure 3. The gas tends to bypass the large cluster, around which the gas has a high radial velocity. In the interior of the cluster, both the radial and axial gas velocity are greatly reduced. Although the orientation of the radial  Figure 4 shows the instantaneous distribution of the gas velocity in the 0.01 m height range, which corresponds to time in Figure 3. The gas tends to bypass the large cluster, around which the gas has a high radial velocity. In the interior of the cluster, both the radial and axial gas velocity are greatly reduced. Although the orientation of the radial velocity becomes chaotic, the larger part of the axial velocity is negative. Only in this area, the gas shows a sign of back-mixing. This proves that the back-mixing of the gas phase is mainly caused by the back-mixing of particles, i.e., the solid phase flow dominates. Except for the dense phase area, the gas velocity tends to be evenly distributed. Gas flow is mainly axially upward, and the radial velocity is small. That is to say, the dilute phase area is generally dominated by the upward gas flow. Thus, most of the dispersed particles can be continuously transported upward. the gas shows a sign of back-mixing. This proves that the back-mixing of the gas phase is mainly caused by the back-mixing of particles, i.e., the solid phase flow dominates. Except for the dense phase area, the gas velocity tends to be evenly distributed. Gas flow is mainly axially upward, and the radial velocity is small. That is to say, the dilute phase area is generally dominated by the upward gas flow. Thus, most of the dispersed particles can be continuously transported upward.  Figure 5 shows the outlet solid flux over time. After about 0.2 s, the two-phase flow within the MFB reaches a relatively stable state. In the relatively stable state, the outlet solid flux still has large fluctuations. The average value is about 7.8 kg/(m2 ·s), and it is much smaller than in the literature [20], which is 90~110 kg/(m2 ·s). This significant difference indicates that the overall drag force calculated with the PCDD model is largely reduced in the current simulations.

Gas-Solid Back-Mixing
The saturated entrainment rate of the particle * s G can be calculated using the following correlation formula of Bai and Kato [29] as  Figure 5 shows the outlet solid flux over time. After about 0.2 s, the two-phase flow within the MFB reaches a relatively stable state. In the relatively stable state, the outlet solid flux still has large fluctuations. The average value is about 7.8 kg/(m 2 ·s), and it is much smaller than in the literature [20], which is 90~110 kg/(m 2 ·s). This significant difference indicates that the overall drag force calculated with the PCDD model is largely reduced in the current simulations.  The saturated entrainment rate of the particle G * s can be calculated using the following correlation formula of Bai and Kato [29] as where Fr is the Frodes number and Ar is the Archimedes number. When the operating gas velocity is 1.52 m/s, the correlated G * s is 16.18 kg/(m 2 ·s). The experimental result given by Li and Kwauk [30] and the simulation result presented by Yang et al. [31] are both approximately 14.3 kg/(m 2 ·s), indicating that the correlation formula has good predictive power. The correlation value at the current gas operating velocity 1.7 m/s is 19.9 kg/(mG * s ·S). Therefore, due to the failure to reasonably simulate the back-mixing behavior of particles, the outlet solid flux is overestimated in [20], which is even much higher than the saturation entrainment rate of particles. The currently selected gas-solid properties are basically the same as those in Yang et al. In this work the operating gas velocity is slightly higher, and the particle filling ratio is slightly lower. However, the simulated average outlet solid flux is much lower than the simulated and correlated values in CFB. The authors speculate that the most direct reason might be the geometry effect. In MFB compared with CFB, the bed diameter is very small, and thus for the wall, the relative area of contact with particles significantly enlarges [11]. The strong wall friction force causes the particles near the wall to be unable to be transported upward, which through particle collisions, gradually passed to the surrounding particles and even the central area. And the gas phase needs to consume a lot of energy or pressure drop for suspended particles, which seriously hinders the transport of the particles. Figure 6 shows the distribution of the mean porosity at the bottom and top sections of the bed. It shows that the top of the bed is a dilute phase area, and the fluctuation range of the mean porosity is small or even close to 1. At the bottom of the bed, the porosity is greatly reduced, and the fluctuation range is large. The lowest value can descend to below 0.8, which indicates that the dense phase can continuously form in different radial positions of the bottom.  Figure 6 shows the distribution of the mean porosity at the bottom and top sections of the bed. It shows that the top of the bed is a dilute phase area, and the fluctuation range of the mean porosity is small or even close to 1. At the bottom of the bed, the porosity is greatly reduced, and the fluctuation range is large. The lowest value can descend to below 0.8, which indicates that the dense phase can continuously form in different radial positions of the bottom.    In CFB, the axial porosity distribution of fast fluidized bed is usually a monotonically increasing exponential distribution unless the outlet solid flux reaches the saturated entrainment rate of the particles. The current simulated axial distribution, although a monotonically increasing function of height, does not strictly obey the simple exponential distribution. It appears to have an inflection point between 0.015125 m and 0.025125 m. Because the parameters of the simple exponential function given by the vast majority of experimental studies are not uniform, it is difficult to compare the current simulated distribution of axial porosity with the results of the literature. Consider the solid phase volume fraction at the bottom of the bed ε sd . As G s < G * s ; the correlation formula of Bai and Kato [29] is used as

Axial Structure
where u t represents the terminal velocity of the particles. In the moderate Reynolds number condition, the formula of Lewis et al. [32] can be used to calculate u t . From (16) the correlated ε sd is 0.22, while the presently simulated value is 0.136. The simulated value is far below the correlated value. The possible reasons for such a deviation include the following. The two-dimensional simulation cannot reach the reality of the threedimensional simulation; the simple set import and export boundary conditions cannot simulate the real complex import and export effects. Additionally, it should be more inclined to think that the relative area of wall-particle contacts increase significantly in the MFB, resulting in a large decrease of the outlet solid flux. It is not difficult to analyze according to Formula (16), and ε sd decreases with the decrease in the outlet solid flux. Due to the low level of ε sd , the heterogeneous distribution of cross-sectional porosity in the bottom area of the bed is weakened. This may be the main reason why the simulated axial voidage does not strictly obey the simple exponential distribution. where t u represents the terminal velocity of the particles. In the moderate Reynolds number condition, the formula of Lewis et al. [32] can be used to calculate t u . From (16) the correlated sd ε is 0.22, while the presently simulated value is 0.136. The simulated value is far below the correlated value. The possible reasons for such a deviation include the following. The two-dimensional simulation cannot reach the reality of the three-dimensional simulation; the simple set import and export boundary conditions cannot simulate the real complex import and export effects. Additionally, it should be more inclined to think that the relative area of wall-particle contacts increase significantly in the MFB, resulting in a large decrease of the outlet solid flux. It is not difficult to analyze according to Formula (16), and sd

Radial Structure
The radial flow law of gas-solid phases in the fast fluidized bed is extremely important. The results of extensive experimental analysis show that as long as the cross-sectional porosity is determined, the radial distribution of the porosity is only a function of the radial position, i.e., the dimensionless radius, regardless of the operating conditions. The cross-sectional porosity at the height 0.015125 m is 0.912. The correlation formula of Patience and Chaouki [33] can be used as the following

Radial Structure
The radial flow law of gas-solid phases in the fast fluidized bed is extremely important. The results of extensive experimental analysis show that as long as the cross-sectional porosity is determined, the radial distribution of the porosity is only a function of the radial position, i.e., the dimensionless radius, regardless of the operating conditions. The crosssectional porosity at the height 0.015125 m is 0.912. The correlation formula of Patience and Chaouki [33] can be used as the following The correlation result of the radial porosity in CFB is then obtained. Figure 8 shows the radial distribution of the porosity at the height 0.015125 m. r can be the dimensionless radius of the bed, and only the results from the right of the MFB can be selected. The simulation results in Figure 8 and the correlation results both show the high porosity in the central and the low porosity in close the wall, which follows the typical core-annulus distribution in the trend. The difference between the correlated and the simulated lies in the relatively low porosity in the central area and the relatively high porosity in the close wall in the simulated results, indicating a weak core-annulus effect. This significant difference lies in the amplification effect of wall friction in the MFB. The narrow bed diameter itself is easy to cause material slug, which can easily fill most of the bed diameter. As has also been stated, there exists the trace of clusters in the central area of the bed. Therefore, the heterogeneity of the radial distribution is also weakened.
The correlation result of the radial porosity in CFB is then obtained. Figure 8 shows the radial distribution of the porosity at the height 0.015125 m. r can be the dimensionless radius of the bed, and only the results from the right of the MFB can be selected. The simulation results in Figure 8 and the correlation results both show the high porosity in the central and the low porosity in close the wall, which follows the typical core-annulus distribution in the trend. The difference between the correlated and the simulated lies in the relatively low porosity in the central area and the relatively high porosity in the close wall in the simulated results, indicating a weak core-annulus effect. This significant difference lies in the amplification effect of wall friction in the MFB. The narrow bed diameter itself is easy to cause material slug, which can easily fill most of the bed diameter. As has also been stated, there exists the trace of clusters in the central area of the bed. Therefore, the heterogeneity of the radial distribution is also weakened.

Conclusions
In this paper, we give a drag model considering the particle environment, which is heterogeneous and complex. This model is implanted into CFD-DEM for simulation of

Conclusions
In this paper, we give a drag model considering the particle environment, which is heterogeneous and complex. This model is implanted into CFD-DEM for simulation of fast fluidization in the MFB of fine particles. Simulations show that gas-solid two-phase flow in the MFB satisfies the common law of fast fluidization to some extent. Because the geometric size of the container is much smaller than that of a CFB, the relative area of the wall-particle contact is large. The wall friction factor becomes very significant, and thus, simulations exhibit a special rule of a large difference with those in CFB. The specific conclusions are drawn as follows: (1) The local structure in the MFB satisfies the natural property of fast fluidized particle agglomeration, forming a disperse dilute phase and continuous dense phase; (2) There is serious gas-solid back-mixing in the MFB, and the dense phase is the main area of gas-solid back-mixing. The wall friction factor aggravates the particle remixing effect, resulting in a relatively low outlet solid flux; (3) The axial porosity presents an increasing distribution with the bed height but does not strictly satisfy the monotonic exponential distribution. The solid volume fraction at the bottom of the bed is much lower than the correlated results for a CFB; (4) The radial porosity exhibits a weak core-annulus structure with a higher central region and lower side region. Compared with the correlated results for a CFB, the central region in the MFB is relatively dense, while the side region is relatively dilute; (5) All distinct variation results of the present simulation can be successfully explained using the relatively strong friction factor. Thus, the present drag model and simulation are both validated and effective, at least in the sense of a qualitative trend.

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

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.