Spherical versus elongated particles – Numerical investigation of mixing characteristics in a gas fluidized bed

The possibility to offer good intermixing between particles is one of the main properties that make ﬂu- idized beds such an important industrial appliance. In this work, we use CFD-DEM simulations to compare mixing characteristics of spherical (AR-1) to elongated spherocylindrical particles (AR-4) of aspect ratio In simulation of AR-4 particles, single-particle and multi-particle correlations for hydrodynamic forces are tested. The results show that elongated particles have more vigorous intermixing and lower mixing times compared to spherical particles. Multi-particle correlations have a slight effect on particle mixing, and they increase the difference between AR-1 and AR-4 particles at higher gas velocities. Including hydrodynamic lift force and torque in the case of AR-4 particles leads to more vigorous mixing and lower mixing times. (cid:1) 2020 The Authors. Published by Elsevier Ltd. ThisisanopenaccessarticleundertheCCBYlicense(http://


Introduction
In industry many processes rely on manipulation of solid granular materials suspended by a gas, such as gasification, drying, coating, food processing and gas phase polymerization. For this kind of application fluidized beds are irreplaceable equipment, offering high levels of contact between gas and solids together with high levels of intermixing between particles. The mechanism of particle mixing in fluidized bed is closely related to the bubble formation and motion (Kunii, 2013;Rowe and Nienow, 1976;Oschmann et al., 2014).
Being such an important parameter for fluidization, the mixing characteristics of spherical particles, mostly A and B powders, have been subject of many numerical and experimental studies and are generally well understood (Rowe and Nienow, 1976;Thiel and Potter, 1978;Shi and Fan, 1984;Peeler and Huang, 1989;Shen et al., 1995;Rhodes et al., 2001;Deen et al., 2010;Banaei et al., 2018). Because many industrial processes use particles of nonspherical shape, recent years have witnessed an increase in development of numerical models for simulation of fluidized beds with non-spherical particles (Lu et al., 2015;Zhong et al., 2016;Vollmari et al., 2016;Mahajan et al., 2018a). Fluidization behavior of such particles is still not well understood, especially when it comes to fluidization in three dimensional (3D) fluidized beds. In our previous work we showed that elongated particles behave considerably different from spherical ones when suspended in gas flow  as they experience different particle-particle interactions together with additional hydrodynamic forces. Therefore findings that apply to spherical particles cannot be automatically assumed to be valid for elongated ones. So far, the only investigation that looked into mixing characteristics of non-spherical particles was done by Oschmann et al. (2014). Their work was the first that indicated a difference in mixing behavior between spherical particles and non-spherical particles of different aspect ratios. However, the exact cause of this different mixing behavior is not yet understood and needs more extensive research. Their work also did not take into account effects of lift force and hydrodynamic torque, both of which have proven to be important for accurately predicting the behavior of elongated particles in fluidized beds . Recently, Sanjeevi and Padding (2020) developed multiparticle correlations for drag and lift force and hydrodynamic torque for elongated particles of aspect ratio 4. These correlations take into account the effect of surrounding particles the so-called crowding effect on hydrodynamic forces and torque experienced by the particles. It is known that using multi-particle correlations leads to better predictions of the particles average orientation and average velocity in z-direction , however their effect on particle mixing is still not investigated.
In this work, we apply the computational fluid dynamics, coupled with a discrete element model (CFD-DEM), which was used and validated in previous works (Mahajan et al., 2018a;, to compare mixing characteristics of spherical aspect ratio 1 (AR-1) particles to elongated, spherocylindrical particles of aspect ratio 4 (AR-4) for different gas excess velocities. We also investigate the effect of multi-particle correlations for hydrodynamic forces and torque and show the importance of the effects of lift force and hydrodynamic torque on the mixing properties of elongated (AR-4) particles.

Numerical model
For this study, we used a CFD-DEM algorithm based on open source CFDEM coupling (Kloss et al., 2012). This package couples two open source softwares, OpenFOAM which uses a CFD solver to solve the fluid equations and LIGGGHTS which solves particle equations using a discrete element method (DEM). These codes have been adapted so that they can deal with spherocylinders, without relying on a multi-sphere approach. More information about the model and its validation can be found in previous works (Mahajan et al., 2018a;.

Discrete Element Model (DEM)
To simulate interactions between solid particles we have applied the discrete element model (DEM), a soft contact model first introduced by Cundall and Strack (1979) to describe interaction between granular particles. In this model trajectory of each particle is numerically integrated over time, subject to local forces and torques. The translational motion for particle i can be calculated by integrating the expression where the sum runs over all neighbours j in contact with particle i; F ij;n is the normal contact force acting on particle i due to its interaction with particle j; F ij;t is the tangential contact force acting on on particle i due to its interaction with particle j; F i;f is the total hydrodynamic force acting on the particle, F i;p represents the pressure gradient (buoyancy) force acting on the particle and F i;b is any body force acting on the particle including gravity. The rotational motion of a particle i can be solved using where I i is the particle moment of inertia tensor, x i is the angular velocity of the particle, T ij is the contact torque acting on the particle i due to its interaction with neighbouring particle j, and T i;f is the fluid-induced pitching torque. Note that for spherical particles, the contact torque T ij is only caused by tangential contact forces. However, for non-spherical particles the torque is caused by both tangential and normal contact forces. Inter-particle forces develop only when particles spatially overlap. Two adjacent spherical particles are overlapping when the distance between their centers is smaller than the sum of their radii. Calculating overlap in the case of spherocylindrical particles is more complicated as it requires the distance between their shafts to be smaller then sum of their radii. An algorithm for calculating the shortest distance between shafts is presented by Vega and Lago (1994), and more detailed information about its application in the used model can be found in Mahajan et al. (2018a), .
To calculate the normal contact force exerted on particle P 1 by particle P 2 we use a linear spring-dashpot model such that the normal contact force is given by F 12;n ¼ Àk n d n n 12 À g n v 12;n ð3Þ where k n is the normal spring constant, g n is the normal damping coefficient and v 12;n is the normal relative velocity between the particles. As shown in Eq. (3), the normal contact force is dependent on the degree of overlapping distance d n , which is calculated with the aforementioned collision detection scheme for spherocylinders.
The magnitude of the tangential contact force is calculated from the Coulomb-type friction expression In this expression t 12 is tangential unit vector defined as t 12 ¼ v 12 = v 12 j j; k t ; d t ; g t ; l and v 12;t are the tangential spring constant, tangential overlap, tangential damping coefficient, friction coefficient and tangential relative velocity respectively. d t is calculated from the time integral of the tangential relative velocity since the development of the initial particle contact and given by where t c;0 is the time of initial contact between the particles.

Computational Fluid Dynamics (CFD)
In CFD-DEM the fluid phase is treated as a continuum, and is described on the basis of the volume-averaged Navier-Stokes equations, which are discretised on a uniform grid. The equation of continuity and momentum conservation are given by where f is the fluid volume fraction, q f is the fluid density, v f is the fluid velocity, s f is the stress tensor for the fluid phase, g is gravity, and R f ;p represents the momentum exchange between the fluid and particle phase, expressed as: where p is the particle label, N p is the number of particles in the computational fluid cell, F p D is the drag force acting on particle p due to the fluid, F p L is the lift force acting on particle p due to the fluid, and V cell is the volume of the computational fluid cell. We do not consider two-way coupling of the torque since it has a negligible localized effect on the fluid.

Hydrodynamic forces
The interaction between the fluid phase and solid particles is resolved trough closures for hydrodynamic forces. The main driving force for fluidization is the drag force and it is considered for both AR-1 and AR-4 particles, while lift force and hydrodynamic torque are considered only for AR-4 particles. For calculating hydrodynamic forces in case of AR-4 particles we applied two approaches: with single particle correlations and with multiparticle correlations. Correlations applied in different cases explored in this work are listed in Table 1.
Drag force. Single-particle drag correlation In the case of AR-1 particles and AR-4 particles with single-particle correlations, the drag force experienced by a particle is calculated using Felice's (1994) modified drag force expression: where v f is the fluid velocity interpolated to the location of particle i; C D is the drag force coefficient, d p the volume-equivalent particle diameter, v i is the velocity of particle i and v is the Di Felice correction factor given by v ¼ 3:7 À 0:65 exp À 1:5 À log Re where the particle Reynolds number Re is calculated as =l f with q f being the fluid density and l f the fluid viscosity. The Di Felice drag force expression was developed to take into account the effect of surrounding particle (crowding effect) on the drag force experienced by a particle. Even though it was originally developed for spherical particles, the Di Felice expression continued to be applied in simulations of elongated particles fluidization (Oschmann et al., 2014;Vollmari et al., 2016;Mahajan et al., 2018a;Ma and Zhao, 2018;Shrestha et al., 2020) because it was until recently the only option for approximating the crowding effect. Comparison between the effects of the Di Felice approximation and a few other more recently developed models for approximating the crowding effect on drag force can be found in work by Mahajan et al. (2018a). Even though the Di Felice expression was applied, because the drag force coefficient is calculated using the single particle, Hölzer and Sommerfeld (2008) correlation, simulations with AR-4 particles done using this approach are considered as singleparticle cases. The Hölzer and Sommerfeld correlation can be applied to arbitrary shaped particles where the shape of the particle is taken into account through sphericity, and lengthwise and crosswise sphericity: where the particle Reynolds number Re p is calculated using the expression defined after the Eq. (10). The particle sphericity (U), the lengthwise sphericity (U k ), and the crosswise sphericity (U ? ) are calculated as: where V p is the volume of the particle, A p is the surface area of the particle, and A e is the cross-sectional area of the volume equivalent sphere. For a sphere U ¼ U k ¼ U ? ¼ 1.
The multi-particle correlation defined by Sanjeevi and Padding (2020) calculates the drag force experienced by a particle as: where F D;/ is the average drag (normalised by the drag on an isolated volume equivalent sphere) based on the particle orientation with respect to the fluid flow (/). Sanjeevi et al. showed that the average drag F D for different / follows a sine-square interpolation for individual particles as well as assemblies (Sanjeevi and Padding, 2020;Sanjeevi et al., 2018), where the average drag F D for any / can be calculated as: The corresponding terms are as follows: where C d;isol is calculated as proposed by Sanjeevi et al. (2018): where the coefficients (a 1 . . . a 5 ) for parallel (/ ¼ 0 o ) and perpendicular (/ ¼ 90 o ) orientation are given in Table 2.
The coefficients for Eqs. (20) and (21) for parallel and perpendicular orientation are also given in Table 2. Even though the multi-particle correlation takes into account the effect of surrounding particles on the drag force experienced by a particle, sub-grid inhomogeneities such as channeling, which are known to occur in beds of elongated particles at velocities around U mf , are not taken into account by this correlation.
In case of AR-4 particles, next to drag force, lift force and hydrodynamic torque are also considered. Their effect on the fluidization characteristics of AR-4 particles has been studied in our previous work .
Lift force. Single-particle lift correlation The magnitude of the shape induced lift force F L experienced by an isolated particle is expressed as where C L is the lift force coefficient. As the lift force is perpendicular to fluids relative velocity v 0 fi ¼ v f À v i and lies in the plane defined by the particle long axis orientation vector u i and v 0 fi , the lift force magnitude F L is multiplied by the lift force orientation vectorê L 0 which is given aŝ The resultant lift force experienced by a particle is then expressed as F L ¼ F LêL 0 . Table 1 Correlations for hydrodynamic forces applied in this study.
Sanjeevi and Padding (2020) The lift force coefficient is approximated using the correlation by Zastawny et al. (2012) Fitting coefficients used for the correlation can be found in Table 3 Multi-particle lift correlation. Sanjeevi and Padding (2020) proposed a simplified calculation for shape induced lift force at different inclination angles /, based on its relation to the normalized drag force. The average lift force F L (normalised by the drag on an isolated volume equivalent sphere) experienced in a multiparticle system at different / is calculated as: The magnitude of multi-particle lift force is calculated as: Hydrodynamic torque. Single-particle torque correlation. The magnitude of the hydrodynamic torque on an isolated AR-4 particle is calculated as where C T is the torque coefficient. The hydrodynamic torque acts perpendicularly to the plane of particle relative velocity and particle orientation vector. Hence, the torque orientation vectorê T 0 is given bŷ The resultant torque is then expressed as T p ¼ T pêT 0 . For approximating the torque coefficient on an isolated particle we have also applied a correlation derived by Zastawny et al. (2012) Table 2 Coefficients for drag force calculation as proposed by Sanjeevi et al. (2018), Sanjeevi and Padding (2020).
The fitting coefficients used for calculating lift and torque coefficient have been derived by in-house direct numerical simulations (DNS) specifically for AR-4 spherocylindrical particles. More information about the DNS simulations can be found in Sanjeevi and Padding (2017), Sanjeevi et al. (2018). The coefficients used for Zastawny et al. correlations for lift and torque are listed in Table 3 Multi-particle torque correlation. The magnitude of the multiparticle torque proposed by Sanjeevi and Padding (2020) is calculated as: where T P;/ is the average hydrodynamic torque for a multi-particle system, calculated by: Coefficients for Eq. (33) are given in Table 4.

Simulation parameters
In this study we used a square fluidized bed which is the same as used in our previous numerical  and experimental work . Column dimensions and main CFD-DEM simulation parameters are presented in Table 5. The parameters for particle properties are determined experimentally by Mahajan et al. (2018b) for particles made of alumide, a 3D printing material consisting of a mixture of nylon and aluminum dust. In this work we compare the fluidization mixing characteristics of spherical AR-1 particles to elongated, spherocylindrical AR-4 particles, using single-particle (SP) and multi-particle correlations (MP) ( Table 1). The considered particles are volume equivalent to each other as this allows us to compare beds of the same mass and same number of particles which also have relatively similar minimum fluidization velocity. The minimum fluidization velocities were determined experimentally  and particle properties are listed in Table 6. In the initial bed, particles are separated by color in two layers, each containing the same number of particles, as shown in Fig. 1. In both layers the particles have identical properties, so the color distinction was made solely for tracking purposes. Horizontal or lateral mixing can be of interest in wide fluidized beds (Oschmann Generally, in CFD-DEM simulations of spherical particles the recommended grid size is between 1.6d p and 5d p (Peng et al., 2014). In this work we applied a grid size of 2.83d p which allows us to satisfy standard practice for spherical particles and also have a cell size larger than the length of the AR-4 particle .

Mixing entropy
Literature offers a number of methods to quantify the degree of particle mixing and they all have their advantages and disadvantages (Wen et al., 2015). In fluidized bed investigations two approaches are most popular in recent research: the Lacy mixing index and mixing entropy. For this investigation we decided to focus on mixing entropy as it offers slightly smoother curves and the only parameter that needs to be considered carefully is the choice of the grid size used for measuring the particle fractions (Wen et al., 2015).
Calculating the mixing index, the domain is divided into grid cells indexed by triplets of integers (i,j,k). For each cell, the local mixing entropy is calculated using: where x 1 i; j; k ð Þ and x 2 i; j; k ð Þ are number fractions of the colored particles in the cell under investigation.
Local mixing entropies are summed in order to obtain the total normalized mixing entropy (S tot ) for each time during the simulation.
where S 0 is the maximum entropy of a randomly mixed system, which depends on the ratio of the number of the two types of particles in the system. For the 1:1 ratio used here, we have S 0 ¼ À ln 1=2 ¼ ln 2. In Eq. (35), we weigh each local mixing entropy with the number of particles N cell i; j; k ð Þ in cell (i,j,k) to properly average the total mixing entropy.
When the size of the grid cell used for measuring the fractions is too small, the measurement of the mixing entropy is hampered by statistical noise. Conversely, when the size of the grid cell is too large, the spatial distribution of the particles is ignored. For the dimensions of fluidized bed considered in this work, dividing the domain in coarse grid of 5 Â 5 Â 35 cells proved to be optimal for accurately determining the mixing entropy.
Mixing time. The mixing entropy gives valuable information on how mixing between particles progresses over time. This also enables us to estimate the mixing time for each case. The normalized mixing entropy varies between values of 0 and 1 for com-  pletely separated and perfectly mixed systems, respectively. Because in practice perfect mixing is never reached, in this work we use a value of 0.95 for the mixing entropy to define a sufficiently mixed system. The same approach was also applied by Deen et al. (2010). The mixing time t 95% is defined as the time at which the mixing entropy reaches a value of 0.95, as shown by a horizontal dashed line in Fig. 2.

Results
In this section we present results on particle mixing in terms of mixing entropy and mixing time. We compare mixing characteristics of AR-1 to AR-4 particles at 6 different gas inlet velocities (1:6U mf ; 1:7U mf ; 2U mf ; 2:3U mf ; 2:6U mf and 2:9U mf ). For AR-4 particles we look into the effect of single-particle (SP) and multiparticle (MP) correlations, as well as considering hydrodynamic lift force and torque at 4 gas inlet velocities (1:6U mf ; 1:7U mf ; 2U mf and 2:9U mf ). All simulations were run for 20 s which for this size of fluidized bed is enough to reach statistical steady-state . Fig. 2 shows mixing entropies as a function of time for AR-1 particles and AR-4 particles with single-particle and multi-particle correlations, under the influence of different gas velocities. It can be seen that varying the gas velocity between 1:6U mf and 2U mf leads to a considerable change in mixing entropy. This is particularly the case when changing from 1:6U mf to 1:7U mf . However with an increase of gas velocity above 2U mf no additional effect on the evolution of the mixing entropy can be noticed anymore. This can be observed for all considered cases. By comparing the mixing entropy curves for AR-4 particles, it can be seen that using singleparticle or multi -particle correlations for hydrodynamic forces does not have a considerable effect on particle mixing. Hovewer, the difference in mixing entropy between AR-1 and AR-4 particles is more evident.

Particle mixing
We now make a more quantitative comparison of mixing entropies for AR-1 and AR-4 (MP) particles for the two lower gas velocities 1:6U mf and 1:7U mf and the highest gas velocity 2:9U mf , as shown in Fig. 3. Solid lines represent AR-4 particles while dashed lines represent AR-1 particles. Different gas velocities are distinguished by color. From Fig. 3 clear differences between AR-1 and AR-4 particles can be observed for 1.6 and 1:7U mf where AR-4 particles show considerably higher mixing entropy values and therefore more intensive mixing. This behavior was already observed by Oschmann et al. (2014). However, at the highest considered gas velocity of 2:9U mf the difference between mixing of the two types of particles is becoming negligible.
Condensing the mixing entropy curves into one value as mixing time (t 95% ), as explained in Section 4, gives more insight in the comparison between AR-1 and AR-4 particles. Fig. 4 shows the time necessary for a fluidized bed to reach a sufficiently mixed state, with a normalized mixing entropy value of 0.95, as a function of gas inlet velocity. Here the difference between considered cases becomes even clearer. At higher gas velocities it can be noticed that even though differences in the mixing time between the two types of particles is getting lower, AR-4 particles continue to show lower mixing times. The effect of gas velocity on the mixing time is considerably higher for for AR-1 particles than for AR-4 particles in the range between 1:6U mf and 2U mf . Choosing between single-particle or multi-particle correlations for AR-4 particle has only a slight effect on the mixing time. However, it can be noticed that MP cor-relations increase the difference in mixing time between AR-1 and AR-4 particle at higher gas velocities (> 2U mf ).
The previous work by Oschmann et al. (2014) already reported less intensive mixing for spherical particles and for particles of lower aspect ratio. However the physical reason for this difference is still not clear. Figs. 5 and 6 show a visualization of fluidization of AR-1 and AR-4 particles at the lowest (1:6U mf ) and highest (2:9U mf ) gas velocities at different time instances. From Fig. 5 it is already visually clear that AR-1 particles show considerably less intensive mixing after a given amount of time. It seems that for the same gas velocity (U=U mf ), AR-4 particles have a higher bed expansion with more vigorous mixing, while AR-1 particles show distinct layers of the same particle colors. The reason for this can be that large (Geldart D type) AR-4 particles show more intensive turbulent fluidization behavior than volume equivalent AR-1 particles for the same gas velocity . In our previous experimental investigation , we have shown that while AR-1 particles remain in a constant slugging regime, AR-4 particles show more turbulent behavior and periodically switch between slugging and turbulent fluidization. As the bubbles are main carriers responsible for particle mixing in fluidized beds and the main mixing occurs in the bubble wake (Kunii, 2013;Rowe and Nienow, 1976), it is expected that fluidized beds with flat nose, raining slugs, that are characteristic for large Geldart D particles (Kunii, 2013), cannot lead to as effective particle mixing as bubbling or turbulent fluidized beds.
From Fig. 6 it can be observed that at the highest gas velocity of 2.9Umf, AR-1 particles are also in a much more turbulent fluidizing regime, more similar to AR-4 particles. Fig. 7 gives more insight in the void fraction distribution for the lower part of the bed (z < 0.3 m), where for all gas velocities we expect to be in the dense part of the fluidized bed. As void fractions between 0.4 and 0.5 respond to densely packed beds, and between 0.9 and 1.0 to cells with barely any particles in them, from Fig. 7 it can be seen that AR-1 particles generally show a much more binary distribution of void fractions (i.e. either very dense or very dilute) compared to AR-4 particles which show a much wider distributions in the middle range of void fractions. At the low gas velocity of 1:6U mf , AR-1 particles are still densely packed, with individual bubbles passing trough the bed and carrying out the mixing. While for AR-4 particles, even at such low gas velocities, the particles are more suspended in the gas flow, which can lead to faster mixing.  This is more visually evident from the snapshots in Fig. 8. It can be seen that in the case of AR-1 particles, bubbles are clearly formed and the boundary between densely packed particles and gas bubbles is clearly visible. In the case of AR-4 particles, no clear bubbles are formed and particles are generally more suspended in the gas flow.
From Fig. 7 it can be seen that with the increase of gas velocity, there is an increase in the intermediate range of void fractions for AR-1 particles, however it still remains lower than in the case of AR-4 particles. AR-1 particles also show a higher peak at void fractions in the range 0.9-1.0, which still suggests that clear bubbles (or slugs) are passing through the bed. With the increase of gas velocity, the distribution of void fractions for AR-4 particles moves to the values corresponding to more dilute cells. Overall, these observations can explain the faster mixing times we observe in the case of AR-4 particles.
In our previous work  we already started to see indications of transition to a turbulent fluidization regime for spherical particles at a gas velocity of 1:9U mf , which can explain the large drop in mixing time for AR-1 particles at 2U mf . This is consistent with our observation that above 2U mf we see negligible difference between AR-1 and AR-4 particles, as both are operating in the turbulent fluidizing regime.

Effect of hydrodynamic lift force and torque on mixing characteristics of elongated particles
So far, most numerical works have neglected the effect of hydrodynamic lift and torque on fluidizaton of non-spherical particles. However, in case of elongated particles these forces have considerable effect on the average particle orientation and velocity distribution . In this section we investigate the effect of considering lift and torque on mixing behavior of AR-4 particles. Fig. 9 compares the mixing entropy and the mixing times for cases with only multi-particle drag force considered to cases where all hydrodynamic forces (using multi-particle correlations) are considered, for different gas velocities. From Fig. 9 (a) it can be observed that cases with all hydrodynamic forces considered (full lines in Fig. 9) show higher mixing entropies than cases where only the drag force is considered (dashed lines in Fig. 9). This difference is diminishing with increasing gas velocity.
Similar behavior can be seen in the mixing time ( Fig. 9 (b)) where for 1:6U mf ; 1:7U m f and 2U mf simulations with lift and torque show lower mixing times. At the highest velocity of 2:9U mf both cases have almost identical mixing times. The fact that after 2U mf we do not see any effect of varying gas velocity (Figs. 3 and 4) or hydrodynamic forces on mixing properties indicates that we have reached a terminal mixing time and that for the considered properties of the fluidized bed (column and particle properties) it is not possible to achieve faster mixing. From the presented results it can be assumed that the terminal mixing times for AR-1 and AR-4 particles in the fluidzed bed considered in this work are around 2.6 and 2.3 s, respectively.
The effect of including lift and torque on mixing time can be attributed mainly to the lift force. Including hydrodynamic lift force leads to an increase in average particle velocity in the vertical z-direction . As can be seen from the previous Section 5.1, increasing gas velocity leads to faster and more intensive mixing.

Conclusion
In this work, we applied CFD-DEM simulations to compare mixing behavior of spherical (AR-1) particles to elongated spherocylinders (AR-4) in a fluidized bed and to investigate the effect of novel multi-particle correlations for hydrodynamic forces and torque (Sanjeevi and Padding, 2020) compared to widely used singleparticle correlations. Spherical (AR-1) particles show less vigorous mixing than AR-4 particles and have higher mixing times. The largest difference can be seen for gas velocities up to 2U mf . At higher gas velocities, even though AR-1 particles still show longer mixing times, differences can be considered negligible. This indicates that between 2U mf and 2:3U mf the fluidized bed already reaches a terminal mixing velocity and it is not possible to achieve faster mixing times for the given fluidized bed properties. More vigorous mixing for AR-4 particles can be explained by the more turbulent fluidization. These observations are characteristic of tall and narrow beds with coarse particles. Different conclusions are possible for shallower beds with lower bed height to column diameter ratios.
Applying multi-particle correlations for hydrodynamic forces and torque does not have a large effect on particle mixing times. It can be noted that the mixing times for cases with multiparticle correlations at higher gas velocities (> 2:3U mf ) have slightly lower mixing times compared to the cases with single particle correlations, however this effects can be seen as marginal.
Including hydrodynamic lift and torque has an effect on mixing properties of AR-4 particles and leads to somewhat faster mixing times. The main cause for faster mixing times can be seen in the effect of lift force which increases the particle average velocity in the vertical direction.

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