Particle ‐ Scale Simulation of Solid Mixing Characteristics of Binary Particles in a Bubbling Fluidized Bed

: The behavior of solid mixing dynamic is of profound significance to the heat transfer and reaction efficiencies in energy engineering. In the current study, the solid mixing characteristics of binary particles in the bubbling fluidized bed are further revealed at particle ‐ scale. Specifically, the influences of gas superficial velocity, Sauter mean diameter (SMD) in the system and the range distribution of particle sizes on the performance of mixing index are quantitatively explored using a computational fluid dynamics ‐ discrete element method (CFD ‐ DEM) coupling model. The competition between solid segregation and the mixing of binary particles is deeply analyzed. There is a critical superficial velocity that maximizes the mixing index of the binary mixture in the bubbling fluidized bed. Solid mixing performs more aggressive when below the critical velocity, otherwise solid segregation overtakes mixing when above this critical velocity. Moreover, superficial velocity is a major factor affecting the mixing efficiency in the binary bubbling fluidized bed. Additionally, the mixing behavior is enhanced with the decrease of SMD while it is deteriorated in the binary system with a wide range of particle size distribution. Therefore, it is highly recommended to perform a binary particle system with smaller SMD and closer particle size distribution for the purpose of enhancing the mixing behavior. The significant understanding of mixing characteristics is expected to provide valuable references for the design, operation, and scale ‐ up of binary bubbling fluidized bed.


Introduction
Due to the advantages of operational ease and excellent gas-solid contact efficiency, the bubbling fluidized bed (BFB) has been widely adopted in the chemical engineering processes [1][2][3][4]. In practice, the BFB is commonly operated with a wide particle distribution [5][6][7], which may enhance or deteriorate heat and mass transfer inside the fluidizing apparatus [8][9][10]. However, most of the previous studies were mainly focused on the mono-disperse systems [11][12][13][14]. To unveil the relationship between the particle size distribution and physical-thermal properties of gas-solid flows, it is in urgent need to investigate the transport mechanism of binary or poly-disperse systems in fluidizing apparatuses. Nonetheless, such intrinsically chaotic and strong non-linear systems within complex multiphase flows pose great challenges for experimental measurements and semi-empirical correlations.
As the computer capacity and numerical algorithm advance, computational fluid dynamics (CFD) has become increasingly popular in the exploration of gas-solid hydrodynamic mechanisms. According to the different treatments of particles, the existing numerical methods can be generally categorized into two kinds: two-fluid model (TFM) and computational fluid dynamics-discrete element method (CFD-DEM). The solid and fluid phases in the TFM are both treated as continua while the particles in the CFD-DEM are tracked individually. Comparing with the TFM, the CFD-DEM has an intrinsic ability to accurately capturing inter-particle collisions and providing detailed solid information at particle-scale [15][16][17].
Two processes commonly appear in a binary system which includes two groups of solid particles with distinct sizes, i.e., mixing and segregation [18][19][20]. Specifically, the segregation process is preferred in metallurgical industries [21] while the high-efficiency mixing process is required in engineering fields with the aim to achieve uniform heat and mass transfer [22]. As pointed out by experimental measurements, the segregation process is dominant under low gas velocities while the mixing process is dominant under high gas velocities [23]. Lu et al. [24] experimentally and numerically investigated the minimum fluidization velocity by analyzing the pressure drop distribution in a binary mixture fluidized bed and found the mass fraction of small particles in the binary mixture plays an important role in size segregation. Zhang et al. [25] experimentally studied the mixing behavior in a binary BFB with cylindrical-shape biomass and spherical sand and observed a maximum mixing state with the increase of inlet gas velocity. Zhong et al. [26] evaluated the effect of wall boundary conditions on the performance of mixing and segregation processes in a binary BFB by the TFM approach and pointed out that the restitution coefficient exerted an insignificant effect on simulating the mixing and segregation processes. Cooper et al. [27] investigated the bed expansion and bubble behaviors in a binary fluidized bed using the TFM approach and indicated that the binary mixture mixed well in axial direction under different operating conditions. Zhou et al. [28] extended the energy minimization multi-scale (EMMS) drag model to the investigation of binary fluidized bed and found that the drag model has a significant influence on the mixing behavior of binary particles.
The behaviors of binary mixture mixing have been extensively investigated by the TFM approach, where the detailed particle-scale information is unavailable. Nevertheless, a large deviation between the TFM and CFD-DEM was presented in the simulation by Fan et al. [29]. It was found that the TFM failed to reproduce the horizontal segregation across the bed. Annaland et al. [30] found that the TFM method underestimated the segregation rates in the simulation of a binary mixture system. Additionally, they pointed out that the underestimation of the segregation rate can be improved by the discrete particle simulation. Unlike the TFM approach, the particle collisions are calculated directly with different contact models and no additional closure models are required in the CFD-DEM. Therefore, the CFD-DEM approach is adopted in the current simulation. The mixing behavior of binary particles in the fluidizing apparatus has been seldom revealed at particle-scale in the literature. Using the CFD-DEM approach, Zhang et al. [31] evaluated the segregation behavior in a binary BFB with different drag models suitable for poly-disperse systems and found that the mixing process of binary particles at high superficial velocity can be well described by various drag models. Limtrakul et al. [32] found that the density and diameter of particles will significantly affect the mixing process in the binary fluidized bed. Besides, a better mixing performance was achieved in the binary system with fewer differences in particle density and diameter. Nevertheless, the influence factors of the mixing behavior, which have guiding significances for the practical industrial operations, have not been further elucidated. Moreover, the in-depth quantitative analysis of the mixing characteristic at particle-scale is still needed. Accordingly, the current paper focuses on the mixing behavior of binary particles in a BFB at the particle-scale level. The poly-disperse drag model is used to accurately predict the dynamics (i.e., flow pattern, pressure drop, and velocity) of binary particles. The effect of gas superficial velocity on the mixing performance is comprehensively assessed. The competition mechanism between segregation and mixing under the influence of superficial velocity is revealed. Moreover, the influences of Sauter mean diameter (SMD) on the binary system and the range distribution of particle sizes are further explored.
where Vcell and np are the cell volume and particle number in the current computational cell, respectively. fd represents the drag force of particle p. The stress tensor of fluid phase in the momentum equation can be written as: The fluid volume fraction in the conservation equations is calculated as:

Governing Equations for the Particle Phase
In contrast to the fluid phase, each particle in the CFD-DEM is tracked individually in the Lagrangian framework. Considering the effect of gravity, drag force, collision force, and pressure gradient force, the particle motion descript by the Newton's law is given as: where fd, fp, and fc represent the drag force, pressure gradient force, and contact force, respectively. Ip, ωp, and Tp are the particle inertia momentum, angular velocity, and torque, respectively.
The soft-sphere contact force [33], which is proved to be appropriate in the simulation of dense fluid-particle flows in our previous researches [34][35][36][37], is employed in the current study. The contact force consists of normal and tangential components. It can be formulated as: where k, δ, η, and μ stand for the spring constant, overlap, damping coefficient, and friction coefficient, respectively. As the main driving force in the dense fluid-particle flows, the drag force is of great importance to be accurately calculated. The drag force in the translational motion equation of particle is written as: where β represents the momentum exchange coefficient between fluid and solid phases. As pointed out by Beetstra et al. [38,39], the drag force in the binary system will be underestimated for the large particle but overestimated for the small particle. Thus, a modification considering the influence of poly-dispersity was proposed and suggested by Beetstra et al. [38,39]. In the current study, the correction factor is employed to investigate the hydrodynamic characteristics of a binary system [38][39][40]. 2

18
(1 ) Re , (1 ) 0.064 where dp,i is the particle diameter. Re, dSMD, and fpoly stand for the particle Reynolds number, the Sauter mean diameter of binary particles, and the correction factor for the poly-disperse system, respectively.

Numerical Setups
The numerical scheme of CFD-DEM has been well documented in our previous literature [34,41]. The geometry configuration refers to the experimental work conducted by Goldschmidt et al. [23]. The investigated object is a binary BFB with 0.15 m, 0.015 m, and 0.70 m in width, depth, and height, respectively. The density of solid particles is 2.5 g/cm 3 . The diameter of small particles and large ones is 1.5 mm and 2.5 mm, respectively. The total mass of the binary particles in the current cases is 326 g. The computational domain is divided into 20, 2, and 94 elements in width, depth, and height, respectively. The computational cell adopted in the current simulation is about 3-5 times the particle diameter, meeting the requirement of the conventional CFD-DEM [36,42,43]. The detailed physical parameters of small and large particles are summarized in Table 1.

Initial and Boundary Conditions
The initial and boundary conditions should be properly set for the enclosure of the algebraic equations. As illustrated in Figure 1, the small and large particles are randomly scattered in the inferior and superior part of the BFB in the bottom region, respectively. The superficial gas velocity is assigned in the distributor with a fixed value. The zero gradient boundary condition is applied to the outlet of the system. At the wall, a non-slip boundary condition is chosen for the solid phase while the gas phase is simulated with free-slip boundary conditions according to the study by Beetstra et al. [39]. Due to the bubbling fluidized regime, particles in the system will not be escaped from the top of the bed. The fluid time step is set as 1 × 10 −5 s while the DEM time step is chosen as 1 × 10 −7 according to the criterion [44]. Each case in the current study is simulated in a cluster using 16 processors with a total physical time of 20 s.

Model Validation
In order to validate the reasonability of the drag model adopted in the current study, the experimental measurement conducted by Jiang et al. [40] is used. The lab-scale BFB has the dimensions of 0.1 m, 0.34 m, and 0.014 m in width, height, and depth, respectively. Two different particles with the bed mass of 20.2 g and 13.5 g are initially packed in the fluidized bed with equal volume. The diameters and densities of investigated particles are 2.5 mm, 1450 kg/m 3 , and 1.8 mm, 985 kg/m 3 , respectively. The Sauter mean diameter in the validation simulation is approximately 2.15 mm. The computational domain is divided into 20, 100, and 3 elements in width, height, and depth, respectively. The detailed properties of gas and solid phases and numerical parameters are given in Table 2. A total physical time of 30 s is carried out for the statistical procedure in the validation. Table 2. Numerical parameters and physical properties of gas and solid phases, P1 for the white particle and P2 for the green particle.  Figure 2 displays the instantaneous distribution of solid particles in the BFB. Particles with binary sizes are respectively extracted and labeled with two colors, that is, the green and white particles respectively represent the large and small particles. The solid behaviors captured with the current simulation agree well with the experimental measurements by Jiang et al. [40]. A similar bed expansion is shown in the simulation and experiment. Good mixing between large and small particles can both be observed in the simulation and experiment. Different particles in the binary mixture system can be easily extracted owing to the flexibility of CFD-DEM simulation, while a certain system error will be introduced into the experimental extraction by the label matrix method [40]. Therefore, a slight difference in the void fraction can be seen between simulation and experiment. Nevertheless, the discrepancy is within the acceptable level.  Figure 3 shows the density distribution of particle vertical velocity and horizontal velocity with different sizes. The particle velocities are extracted at the heights ranging from 6 cm to 7 cm of the bed. The simulation results modeled with the current drag model are compared with that of the mono-disperse drag model and experimental data. It is obvious that more reasonable velocity distributions are captured when the poly-dispersity effect is considered in the drag model. The particle motions of both large and small particles perform more vigorous in the vertical direction compared with those in a horizontal direction. Moreover, similar distributions of horizontal velocity are obtained both in small and large particles while a slightly larger vertical velocity of small particles is presented in simulation and experiment. Additionally, the symmetric density distributions of particle velocity are presented both in horizontal and vertical directions.  Figure 3. Density distributions of particle vertical velocity and horizontal velocity with different particle sizes: (a) particle P1; (b) particle P2.

Parameters
In order to further verify the accuracy of the poly-disperse drag model in the current work, the mixing index is compared between different numerical models and experiment results. For the purpose of quantitative description, the Lacey mixing index (MI) [45] is adopted to quantify the mixing degree of the whole bed, which is defined as: where, the equivalent number of particles ne in the binary system [40,46] is employed. c represents the averaged concentration of the smaller particles. The variance S is given as: where ci is the concentration of smaller particles in the sample i. ki demonstrates a weighting factor. Obviously, the mixing index is in the range of 0 and 1, and the mixing degree at each instant can be descript by this value. Specifically, the MI of 1 and 0 stands for the perfectly mixed state and the completely segregated state at the initial time, respectively. Figure 4 illustrates the evolutionary profiles of the mixing index between simulations and experiments. In harmony with the experimental measurement, the mixing index in the current simulation is calculated from the initial time. It can be found that the mixing index in the BFB gradually increases at the initial time period. At the time of 2 s, the mixing index achieves an equilibrium state and finally fluctuates around a fixed value. The profile obtained by the polydisperse drag model is comparable to that of the experiment. Moreover, the relatively large deviations between the prediction with the mono-disperse drag model and experiment are presented in the current study. Observing from the instantaneous distributions of solid mixing behavior at 2 s, the poly-disperse drag model well predicts the mixing state between large and small particles while the traditional mono-disperse drag model obviously underestimates it, as the latter model underestimates the solid motion of large particles and overestimates that of the small particles. Based on the qualitative and quantitative analysis, the proposed drag model shows responsibility in the simulation of poly-disperse particle systems.

General Flow Pattern in the Binary Bubbling Fluidized Bed
Exploring the flow pattern of binary particles in the BFB is essential for not only the understanding of mixing behavior but also the design and optimization of the binary system. Figure 5 displays the evolution of solid motion in the binary BFB over time. After introducing the gas flow from the bottom distributor, the bed height is apparently increased with the rapid expansion of mixture particles. Meanwhile, owing to the formation of void structure between the interface of small particles and large particles, the small particles in the bottom region penetrates into the large particles in the upper layer. As the bubble rises and eventually breaks at the bed surface, the large particles are thrown to both sides of the left and right walls. After a period of bubble disturbance from formation to collapse and the influence of solid backmixing along the wall region, the binary system is considered to reach a stable state at 5 s. Comparatively, a better mixing performance can be found in the upper region of the BFB because of the strong restriction of solid motion in the bottom region. Nevertheless, a good mixing state between small and large particles can be observed and the classical bubble phenomenon in the binary BFB is also captured.

Influence of Superficial Velocity on the Solid Mixing
As pointed by some researchers [47,48], elevating the superficial gas velocity can promote solid mixing in the mono-disperse fluidized system. In the current study, the effect of superficial gas velocity on the solid mixing characteristic of the binary BFB is comprehensively investigated. The superficial gas velocity changes from 1.15 m/s to 3.0 m/s while the physical properties of the small and large particle are the same as the description in Table 1.
The pressure characteristic is of great importance to the evaluation of gas-solid behaviors in the fluidization [49,50]. Figure 6 displays the influence of superficial gas velocity on the instantaneous pressure fluctuation and the associated power expressed as mean square amplitude (MSA). It is found that the pressure drop under different conditions randomly fluctuates around a fixed value. The amplitude of pressure drop is gradually strengthened with the increase of superficial gas velocity. The oscillation of pressure drop mainly depends on the disturbance caused by the formation and collapse of bubbles in the BFB. Observing from the power spectra analysis, a peak MSA value under a similar frequency is performed in each case. Generally, a higher MSA value appears at a larger superficial gas velocity.  Figure 7 illustrates the time evolutionary profiles of the mixing index in the binary BFB with different superficial gas velocities. The mixing index gradually increases from zero and then oscillates around a fixed value, which respectively corresponds to the completely separated state at the initial moment and the equilibrium state of a binary system. The larger the superficial gas velocity is, the shorter the time is required for the system to reach an Over the current range of superficial gas velocity from 1.15 m/s to 1.6 m/s, a better mixing behavior is obtained in the system with a larger superficial gas velocity. Owing to the more chaotic and unsteady fluidization in the system with a larger superficial gas velocity [48], the oscillation of MI is a bit more vigorous with the increase of superficial gas velocity. To further explore the influence of superficial gas velocity on the mixing performance in the binary system, a wider range of superficial gas velocities are comprehensively investigated. Furthermore, the independent study of statistical time is conducted in the current study. Figure  8 shows the time-averaged mixing index with different statistical time and various superficial gas velocities ranging from 1.15 m/s to 3 m/s. The results of the initial 5 s are discarded for avoiding the startup effect. In the case of 1.15 m/s, a big deviation can be found with a statistical time of the last 15 s and those of 5 s and 10 s, which indicates that it requires a long time to reach the equilibrium state as illustrated in Figure 7. Beyond that, the discrepancy among the results with different statistical time is fairly insignificant. Therefore, the statistical procedure conducted over the last 10 s of the simulation is reliable in the present study. Interestingly, the solid mixing in the binary BFB is firstly enhanced then deteriorated with the increase of superficial gas velocity. This phenomenon is different from previous observation as documented in the mono-disperse system [15,48].  Qualitative and quantitative analysis about the influence of superficial velocity on the mixing behavior is comprehensively discussed in the current study. Figure 9 displays the instantaneous structures of bubbles in the BFB under the superficial gas velocity of 1.5 m/s and 3.0 m/s, which represent two different representative operating conditions. Owing to the introduction of gas flow, small bubbles are randomly formed at the bottom region and coalesce into larger bubbles in the rising process, and finally break up at the bed surface. In the case with lower superficial gas velocity, the distribution of bubbles is more uniform. At a larger superficial gas velocity, an extremely large bubble can be observed in the bed. As pointed by previous researchers [51], decreasing the size of bubbles is beneficial for the gas-solid interactions especially in the solid mixing. Therefore, the large bubble observed in the system with a larger superficial gas velocity can provide a good descriptive metric for the deterioration of mixing characteristics. As shown in Figure 10, the vertical velocity of the small particle and large particle are quantitatively analyzed in the current study. Compared with the profiles in the middle and upper regions, the solid vertical velocity is more moderate in the lower layer. Owing to the solid back-mixing phenomenon in the upper region, the upward solid motion is more drastic in the central region. In the case with low superficial gas velocity, similar distributions of a vertical velocity of small and large particles can be observed in the central and upper regions while the solid vertical velocity of the small particle is a little smaller than that of the large particle in the lower layer. It is considered that the deviation of solid motion between small and large particles in the bottom region will enhance the disturbance of a binary mixture, thus promoting the mixing process in the binary system. As for the situation with a higher superficial gas velocity, the small particle preforms a larger velocity in the upper layer than that of the large particle. It is considered that this deviation in the upper region will cause a segregation phenomenon between small and large particles, which demonstrates the deterioration of solid mixing under high superficial gas velocities. Generally, the mixing characteristic in the binary BFB is determined by the competition between the solid segregation and mixing. Concluding from Figure 8, there is a critical superficial gas velocity that maximizes the mixing index of the binary mixture in the BFB. Mixing beats segregation and segregation overtakes mixing under the operating conditions below and above the specific critical velocity, respectively.

Influence of SMD on the Solid Mixing
As commonly assumed in the previous researchers [52][53][54][55][56], the SMD is adopted in the simulation when dealing with the poly-disperse systems. Therefore, the influence of the mean diameter assumption in the simulation of the poly-disperse system is comprehensively explored in the current section. Leaving the physical properties (i.e., diameter and density of small and large particles) untouched, the number of small and large particles in different cases is various while the equivalent number of particles in each system is fixed. Table 3 gives detailed information of the simulation settings in the present section. Table 3. Details of the gas-solid properties in the current simulation, P1 for the small particle and P2 for the large particle.  Figure 11 illustrates the time evolutionary profiles of the solid mixing index in the BFB with different SMD. Comparting the equilibrium state under various conditions, the smaller the SMD of the binary mixture is, the better mixing behavior can be achieved in the system. This finding is consistent with the practical engineering applications (i.e., coal combustion, gasification, and chemical looping combustion) where small particles perform a better mixing and heat transfer performance [57][58][59]. Moreover, the mixing behavior is considerably deteriorated in the system with the SMD of 2.23 mm while the oscillation in this case is relatively more drastic compared with other conditions. In addition to the evaluation of the mixing index, the mixing efficiency is also compared under different conditions. Observing from the evolutions, the slopes of mixing index curves at the initial time are quite similar, which indicates that the SMD in the binary system has a little effect on the mixing efficiency. Owing to the intense interphase momentum exchange and inter-particle interaction, the disturbance characteristics of small and large particles are quantitatively analyzed by the granular temperature in the current study. Figure 12 exhibits the time-averaged granular temperature of binary solid mixture with various SMD at the height of 0.05 m and 0.1 m. It is found that the horizontal distributions of granular temperature at different height are quite distinct. Specifically, parabolic shape distributions of granular temperature are observed in the lower region while the profiles in the upper layer are flatter. Moreover, the disturbance of small particles in each case is definitely more drastic than that of large particles. Comparing with the solid disturbance in the near-wall regions, a greater granular temperature magnitude appears in the central part.  For further understanding of the influence of SMD on the mixing behaviors in the binary BFB, the granular temperature difference between small and large particles under different conditions is comprehensively investigated in Figure 13. The granular temperature differences distribute parabolic profile and 'M' shape at the height of 0.05 m and 0.1 m, respectively. Moreover, the granular temperature difference is relatively greater in the upper layer. The granular temperature difference between small and large particles gradually increases with the increase of SMD in the binary system both in the lower and upper layers. A small deviation of the granular temperature indicates the fluctuation of small particles and large particles is closer,

Influence of Particle Size Distribution on the Solid Mixing
In the current section, the influence of particle size distribution on the solid mixing characteristic is comprehensively revealed. Three ranges of the particle size distribution are investigated while the equivalent number of particles and SMD in the binary system are maintained in each case. Detailed properties of the binary mixture are summarized in Table 4. Table 4. Details of the gas-solid properties in the current simulation, P1 for the small particle and P2 for the large particle. 1.88 SMD in case 5 (mm) 1.88 SMD in case 6 (mm) 1.88 Gas velocity (m/s) 1.5 Figure 14 illustrates the time-averaged mixing index in the binary system with different ranges of particle size distribution. Comparing with the system with a wide range of particle size distribution, the mixing index is relatively higher in the system with a narrower particle size distribution. With the decrease of particle size distribution range, the mixing behavior in the binary system is gradually enhanced. A similar phenomenon can be concluded in the previous researches [32,40].

Conclusions
The mixing characteristics of binary BFB are numerically investigated in the current study. The poly-disperse drag model used in this paper can accurately predict the properties of a binary system. The influence of superficial gas velocity, SMD in the system, and particle size distribution on the behavior of mixing index in the binary BFB are comprehensively revealed. The results show that the mixing characteristic in the binary BFB is firstly enhanced then deteriorated with the increase of superficial gas velocity. Generally, the mixing characteristic in the binary BFB is determined by the competition between the solid segregation and mixing. A critical superficial gas velocity exists, which corresponds to the optimal mixing of binary particles. Solid mixing is dominant below the critical velocity while solid segregation is dominant above the critical velocity. The mixing process is enhanced with the decrease of SMD while it is deteriorated with a wide range of particle size distribution. Therefore, it is highly recommended to perform a binary particle system with smaller SMD and closer particle size distribution for the purpose of enhancing the mixing behavior. Superficial gas velocity is a major factor affecting the mixing efficiency, while the SMD and particle size distribution in the binary system have little influence on the mixing efficiency.
Author Contributions: Conceptualization, J.L. and K.L.; data curation, J.L.; methodology, J.L. and S.W.; writing-original draft preparation, J.L.; writing-review and editing, S.W. and L.S.; supervision, K.L.; funding acquisition, K.L. and J.F. All authors have read and agreed to the published version of the manuscript.