DEM-CFD Simulations of Fluidized Beds with Application in Mixing Dynamics †

Discrete element models of fluidized beds allow the accurate reproduction of many aspects of the flow of the fluid and solid phases on both microscopic and macroscopic scales. In the present paper, we first discuss some of the basics of the DEM-CFD model, along with an illustration of areas where it proved successful. Then, a discrete element computational code is used to analyse in detail the key mechanisms governing the mixing of solids in fluidized beds. Air fluidization of a mixture of glass ballotini and steel shots is simulated, and the steady-state concentration profiles are successfully validated with experimental data. By changing the density of the heavier component, the effect of the gas velocity in combination with the density ratio is investigated in terms of an equilibrium degree of mixing and a characteristic time to reach this condition. The maximum mixing achievable is found to depend strongly on the difference of density. For mixtures with a density ratio close to 3, full mixing is practically impossible due to the presence of a non-mixable region, rich in the heavy component, at the bottom of the bed.


Introduction
Fluidized particle bed processing is a particularly efficient technology due to the high quality of fluidsolid contacting, low fluid pressure drop and prevention of dead zones or fluid channelling.It is used in a variety of industries such as oil, petrochemicals, minerals, pharmaceuticals and food processing.Accurate models which provide detailed information on the key phenomena occurring in the bed are necessary to properly design and operate fluidized beds at the desired conditions.Computational models have proved effective in reproducing most of the features on both microscopic and macroscopic scales of complex units involving multiphase flows 1) .This is particularly true for first principles models, which typically are not directly ap-plicable to the real scale of the plant but can produce realistic data crucial for fundamental studies.The flow of particulate solid-fluid mixtures in fluidized beds is characterized by a wide multi-scale nature in both time and, especially, space.There are quick and local phenomena such as the microscopic deformation of the par ticles during collisions as well as the development of the fluid boundary layer around particles, or events occurring at intermediate time and space scales such as bubble formation and growth, up to generally slow processes such as the mixing of solids, which typically involve the whole macroscopic scale of the system.In the literature, computational models of solid-fluid systems have been developed aiming to represent one or more of the mentioned scales.They can be broadly subdivided into three categories: two-fluid, direct numerical simulations and discrete element models.
In an effort to capture the bubbling characteristics of fluidized beds by similarity to the boiling of fluids, two-fluid models (TFM) have been developed 2) .They are based upon the assumption that the fluid and the particles as a whole can both be assimilated to fluids which dynamically interpenetrate, exchanging momentum as well as heat and mass.The continuity and momentum balance equations can be derived by locally averaging the equations governing the fluid and particles motion on a scale much smaller than the macroscopic scale of the system but still larger than the particle diameter 3) .The key point underlying these so-called Eulerian-Eulerian models is the way in which solid-solid and solid-fluid interaction are taken into account.The former is by far the most complex, as Newtonian behaviour requires the definition of a solid-phase pressure and viscosity and, usually, does not adequately represent the actual rheology of granular materials.Significant results in terms of bubbling bed hydrodynamics have been achieved by describing the solid-phase properties based on the Kinetic Theory of Granular Flows (KTGF) 4) .Many of them have been summarized in a book of Gidaspow 1) .However, aspects such as multicomponent fluidization and highly polydisperse systems or fluidized beds of cohesive solids can be analysed by TFM only with extreme difficulty.At the highest level of detail, the microscopic fluid motion (Direct Numerical Simulation or DNS) around each particle can be solved by means of numerical methods on ver y fine computational grids.Then, using the calculated stress tensor acting on each particle's surface to evaluate the drag force, Newton's second law for each particle in the system can be integrated to find the granular phase dynamics.This approach can be classified as really based on first principles, and is therefore considered as a reference for the other models.The huge computational intensity of such a method for large systems involving billions of particles and evolving slowly is evident.This explains why no simulations of representative fluidized beds have been reported in the literature, although fluid-solid flows involving 1204 particles in motion 5) and the sedimentation of 10000 spheres in a liquid have been already analysed 6) .A similar degree of resolution can be achieved using an even smaller scale approach such as that of lattice-Boltzmann simulations 7) , a technique recently utilized to evaluate the drag force on clusters of spheres 8,9) .The discrete element approach represents a promising combination of the two approaches discussed above.It arises from coupling the Discrete Element Method 10) for the solid-phase hydrodynamics with a solution of the locally averaged Navier-Stokes equations for the fluid-phase flow field, similar to two-fluid models.The approach adopted for the solid phase is Lagrangian, and the scale of analysis of the fluid and solid flow field is different.This allows a considerable advantage in that assumptions on the solid-phase rheology is unnecessary, as particles move due to wellcharacterized forces.However, the fluid drag acting on the particles has to be evaluated by a function of the instantaneous particle velocity, of the local fluid velocity and, most critical, of the local voidage around the particle considered.After the pioneering work of Tsuji and co-workers 11) , significant improvement and developments were introduced by Hoomans et al. 12) , who used a hard-sphere collision model, and Xu and Yu 13) , who introduced enhancements in the fluidparticle interaction terms.Since then, a considerable number of authors demonstrated the highly promising potentiality of the model, as recently reviewed by Deen et al. 14) and Zhu et al. 15) .In the present paper, we explore the capabilities of the DEM-CFD approach in the study of the key mechanisms governing mixing and segregation in fluidized beds of two solids differing in density.A brief summary of the DEM-CFD model implemented in our code will be presented first.Then, simulations will be used to represent the behaviour of a mixture of glass and steel particles fluidized by air, validating the results in terms of steady-state concentration profiles with experimental data.Finally, the influence of the gas velocity and density ratio of the mixture on the equilibrium and dynamic degree of mixing will be investigated.

Areas of application of DEM-CFD simulations
As compared to other modelling approaches derived on a larger scale, the DEM-CFD technique allows the relaxation of simplifying assumptions related to the particle-phase stress tensor.Collisions are treated on a mechanical basis and this leads to more realistic reproduction of the actual behaviour of the solids.One example of the power of this method is that there is no need to control values of the voidage in the computational cells.In fact, in TFM simulations, it is possible that unrealistic voidage values occur because they are calculated as continuous variables, so special constraining techniques have to be used to prevent excessively low values 16,17) .In contrast, in DEM this condition is virtually impossible.Another interesting feature of DEM models is the ability to provide results on the particle scale, thereby allowing careful analyses of the transient-and steady-state momentum exchange between the solids.This capability makes it appropriate for applications in the study of surface erosion for immersed objects such as, e.g.heat exchanging tube bundles 18) .The ability to track individual particles in the system also means that solids in which each particle is different from the others can be conveniently analysed.The instantaneous distribution of the components in fluidized beds of mixtures of solids can be obtained, allowing the characterization of typical mixing and segregation patterns (see, e.g.Refs. 19,20,21) )In the light of these promising results, the analysis of elutriation problems in fluidized beds of polydisperse solids as well as the characterization of the fluidization quality enhancement due to the addition of fines appear attractive fields for potential applications of the DEM technique.
Apart from purely hydrodynamic aspects, due to its relevance in industrial applications, heat exchange phenomena between fluidized beds and immersed objects or external walls has long attracted the interest of the research community, beginning with the pi-oneering work of Mickley and Fairbanks 22) .However, despite the efforts, a real comprehensive insight into the phenomena and the ability to predict the overall heat transfer coefficient are still lacking.Further insight into the microscopic heat transfer mechanisms between bubbling beds and immersed objects is likely to be gained in the near future thanks to discrete element simulations, as proved by the encouraging results found by few workers 23,24) .Useful results are also expected to be found in purely fundamental studies on gas-and liquidfluidization.It is well known that gas-fluidized beds of fine and light solids exhibit a distinction between the minimum fluidization and minimum bubbling conditions 25) .This has long been obser ved and studied, and a number of possible interpretations for it have been proposed (see, e.g.Refs. 26,27) .The key point is whether stable and homogeneously expanded beds are actually fluid dynamically stable or whether it is the onset of cohesion, known to be inversely dependent upon the particle size, that is responsible for the homogeneous fluidization.

Solid phase
Among the disadvantages of DEM-CFD simulations, the computational intensity certainly occupies one of the first positions.With the CPU speed and memory resources available today (and in the near future), simulations of large-scale industrial units are not feasible.Probably, pilot-scale fluidized beds already exceed the capability of DEM-CFD-based models.However, distributed parallel computing on grids will certainly be beneficial in this regard.

The DEM-CFD model implemented
A Fortran90 code implementing the DEM-CFD strategy was developed in the last years 28) and details of the model are briefly summarized here.A fully 3D DEM framework, including the possibility to insert internal walls, coupled to a 2D standard CFD algorithm for the flow of a compressible or incompressible fluid through the particulate phase are implemented.A summar y of the equations used is reported in Table 1.Within the DEM context, the translational and rotational motion of each individual particle is calculated by integrating the Newton-Euler differential equations (Eqs.(T1-T2)).The total force is obtained by summation of the gravitational force, the action of the fluid accounted for through a pressure-gradient-related term and a drag force contribution Eq. (T5) 29) and, finally, the sum of the forces exerted by each of the bodies in contact with the particle considered, evaluated by a soft-sphere approach as in Eqs.(T3-T4) 30) .As is common practice, additional contributions such as added mass, Basset's history, lift forces and so on are neglected here under the assumptions that for gas fluidization, the case under examination in the present work, they are very small compared to the other terms.For the fluid-phase flow field, the spatially averaged continuity and momentum balance equations for multiphase fluid-solid systems are considered, following the approach proposed by Anderson and Jackson 2) .Since the solid phase is dealt with via DEM, only the fluid-phase part of the original study is considered, using the coupling strategy as in Eq. (T9) and discussed elsewhere 31) .As far as post-processing of the results is concerned, mixing and segregation patterns in both steady-state and transient conditions are analysed in terms of a mixing index, as introduced by Lacey 32) and defined as: where σ 2 indicates the variance of the concentration distribution in the system at the time instant considered.The subscripts S and M denote perfectly segregated and mixed conditions, respectively, whose variances can be calculated by analytical expressions.The variance σ 2 of the system is estimated by decomposing the system into rectangular areas that are small compared to the macroscopic size of the system (one-tenth of the system width by roughly one twenty-fifth of the bed height in the present study), calculating the volumetric concentration of one component (flotsam in the present study) in each area and finally evaluating the average concentration and its variance in the distribution.The procedure used allows calculation of the instantaneous mixing index even in the presence of bubbles, as each small area adapts its height to include a fixed number of particles, thus virtually eliminating void zones.

Mixing of Solids in Fluidized Beds
In multi-component fluidized beds, the degree of mixing of the components is a complex function of the particle properties and operating parameters.In order to make use of DEM-CFD simulations to improve knowledge and understanding on mixing/segregation phenomena, a system composed of 15000 monosized spheres representing a two-component mixture of glass ballotini and steel shots was simulated.To distinguish between the two components, we follow the notation introduced by Rowe et al. 33) , denoting the former as flotsam ( f ) and the latter as jetsam ( j ).
A preliminar y verification of the capabilities of DEM-CFD simulations to capture the essential features of the mixing in fluidized beds of two solids is carried out by comparing some experimental measurements with model predictions.
In order to avoid excessively large CPU times due to a fully 3D description of the system, in simulations the system is confined to a pseudo-3D geometry, i.e. spherical particles are allowed to move in a domain whose thickness equals the particle diameter, so that a limited number of particles are sufficient to represent the behaviour of the real system.However, the constraint for particles to move only along a plane determines a generally higher voidage than in real experiments, so that even key variables such as the minimum fluidization velocities can be significantly different.This prevents the possibility of directly comparing the phenomena occurring at the same absolute values of the superficial velocity, suggesting that proper scaling of variables is necessary.For this purpose, and considering that the mixing patterns in a system can be expected to be observable at velocities higher than umf,j for that system, the ratio of the actual velocity to this reference value u* = u / umf,j appears a suitably scaled variable.

Comparison with experimental data
Experiments were performed 34) on a mixture of glass ballotini and steel shots of comparable Sauter mean diameters (Df = 450μm and Dj = 433μm) in a 10-cm diameter column with a packed bed aspect ratio of 1.7.The axial profile of component concentration was determined by instantaneously cutting off the fluidizing gas feed, withdrawal of solids from the top of the column in horizontal layers 1.5-cm thick and subsequent measurement of the flotsam volumetric fraction in each layer.
Simulations were carried out considering a system similar to the one used in experiments, whose details are listed in Table 2.An initial packed bed condition was generated by letting particles fall randomly until complete settlement.As expected, due to the pseudo-3D geometr y, the average voidage obtained (ε0 = 0.447) is noticeably dif ferent than that measured experimentally (ε0 = 0.403).The minimum fluidization velocities of the jetsam component, calculated by using the pressure drop relationship proposed by Di Felice 29) , corresponding to these two voidage values are 68 cm s -1 and 45 cm s -1 , respectively, and will be used as references for the calculation of the dimensionless velocities u*.Simulation results and experimental data were compared at various gas superficial velocities, starting from an initial flotsam-on-jetsam segregated configuration and extracting the component distribution after steady-state conditions are attained.Fig. 1 shows the axial concentration profiles corresponding to dimensionless velocities around 1.17 and 1.31 ob-  2.
tained after about 30 s for the experiments and 25 s for the simulations.Note that significant differences exist in the actual velocities necessary to obtain similar u*, e.g.52 cm s -1 in experiments and 80 cm s -1 in simulations for the former value.Despite this, the agreement between the simulations and experimental observations may be considered as good.At the lower velocity, the model is able to capture the small amount of jetsam uniformly distributed in the flotsam (0.95 < xf < 1 in the upper region of the bed).At the higher velocity, the experimental system is vigorously fluidized, exhibiting a small jetsam-rich region at the bottom of the bed and a roughly uniform upper zone, although the flotsam fraction tends to increase slightly moving towards the top of the bed.The simulated bed shows a concentration profile that follows the same trend, though with a slightly worse agreement than the low-velocity case.
In order to check the influence of the initial configuration adopted on the steady-state degree of mixing, it is useful to compare the time evolution of the system starting from a fully segregated or fully mixed bed.However, although concentration profiles provide detailed information on the mixing and segregation patterns, they prove inconvenient to carry out analyses in transient conditions.For this purpose, an overall measure of the mixing state such as the mixing index defined in Eq. ( 1) appears more appropriate.
In Fig. 2, the time evolution of the instantaneous mixing index is as obtained by post-processing DEM-CFD results, for velocities of 80 cm s -1 (u* = 1.18) and 90 cm s -1 (u* = 1.32) and starting from the two extreme configurations.It is shown that for both values of the superficial velocity, a similar ultimate value of the mixing index is found irrespective of the initial MI.Similar indications were obtained following an analogous procedure by Feng et al. 35) , by DEM simulations and experiments.They found that the degree of mixing of two solids differing only in size eventually reaches a value of dynamical equilibrium, irrespective of the initial condition, but depending on the gas velocity.It is therefore expected that a unique equilibrium state exists at a given velocity for mixtures of solids differing in both density and size.
The presence of an equilibrium state allows measurements and simulated results obtained at different times, typically of the order of minutes for experiments and tens of seconds for simulations, to be compared as in Fig. 1, as long as transient phenomena have vanished.The overall agreement confirms that despite the unavoidable simplifying assumptions, the model is able to reproduce the key mechanisms by which mixing or segregation occurs in multi-component fluidized beds.

Mixing dynamics: effect of the gas velocity and the density ratio
To investigate the hydrodynamic effects on solids mixing, simulations were carried out changing the gas velocity from values above umf,j up to u* = 1.47 (Fig. 3), the latter corresponding to a condition where particles reach the top of the simulated domain.It can be observed that around a dimensionless velocity of about 1.25, the equilibrium state of the bed changes from predominantly segregated at low velocities to well mixed at higher values, i.e. when particles undergo vigorous bubbling.and segregated (MI ≈ 0) configuration at the two dimensionless velocities indicated.Simulation conditions are reported in Table 2.
Fig. 3 Time evolution of the mixing index at the indicated gas dimensionless velocities.Simulation conditions are reported in Table 2.
A clearer picture is obtained by looking at the instantaneous particles' positions, flotsam composition colour map and axial profile along the dimensionless height, as reported in Fig. 4 for gas velocities of 85 and 100 cm s in black.The colour maps (Figs.4c,d) show the subdivision of the particle bed into small areas containing approximately the same number of particles (60 in this case) and, on a colour scale, the flotsam concentration in each area.It is worth noting how bubbles are dealt with through higher cells.The same arrangement is utilized to calculate the variance of the distribution and, thus, the mixing index.Using the same subdivision and averaging the concentrations layer by layer, a volumetric fraction profile along the bed height (Figs.4e,f) is obtained.
With respect to the steady-state information contained in Fig. 3 (MI* = 0.62 at u* = 1.25 and MI* = 0.89 at u* = 1.47), this set of plots provides a more complete view on the state of the system.In particular, in the lowermost part of the bed, a region rich in jetsam is found for both velocities, despite the fact that the higher velocity is well above the jetsam minimum fluidization velocity.The presence of these zones determines a large deviation of the concentration profile from the value exhibited in the rest of the bed (Figs.4e,f).Consequently, the large variance of the concentration distribution leads to values of the mixing index that can hardly approach unity, especially at u* = 1.25.The glass-steel particle mixture showed that even at a velocity considerably higher than umf,j, full mixing can be extremely difficult to accomplish.This could be attributed to the high density ratio of the mixture (about 3.1) that is known to drive the system towards segregation.It therefore appears interesting to analyse the effect of the density difference in combination with that of the velocity on the equilibrium degree of mixing.For this purpose, a comprehensive set of simulations was carried out considering mixtures of particles with the same properties as those reported in Table 2 except for the density of the jetsam component, which was allowed to vary.In particular, considering a reference density of the glass ballotini, jetsam-on-flotsam density ratios ranging from 1 (degenerate case) to 3 were set by increasing the density of the other component.The corresponding minimum fluidization velocity of the jetsam, calculated at the packed bed voidage, increases accordingly from about 27 cm s -1 to about 66 cm s -1 .Table 3 reports the values of the dimensionless velocities for all cases analysed and serves as a reference for the results listed in Tables 4 and 5.
Analysis of the results allowed recognition, as a common feature, of a time evolution of the mixing index closely resembling profiles typical of an exponential dependence towards an asymptotic value.The two parameters of this type of curve are the steadystate mixing index MI* and the characteristic time τ.Thus, the simulated evolutions of the degree of mixing of the system, starting from a fully segregated configuration, were post-processed to extract these two parameters.The steady-state mixing index resulted from time averaging the values on the plateau, whereas a non-linear fit procedure was applied in each case to determine the characteristic times.
The results for MI* are collected in Table 4, from which, despite the coarse discretization of the velocity range, it can be seen that the highest values of the mixing index are achieved at velocities that increase with the density ratio.This trend can also be observed in terms of dimensionless velocity u*.For example, a mixing index of 0.83 is found at u* = 1.22 for a mixture withρj/ρf = 1.8 whereas, to achieve a similar index for a mixture with the maximum density ratio considered, a velocity u* higher than 1.44 is necessary.The characteristic times resulting from the fitting procedure are reported in Table 5.It shall be remarked that, different from MI*, the characteristic time depends on the initial condition.As already mentioned, simulations were carried out starting from an initially segregated configuration, so thatτ assumes the meaning of characteristic mixing time.For most of the systems considered, monotonically decreasing mixing times are found as the velocity is increased.
In particular, it is interesting to note how, for a given mixture, characteristic times may change by more than one order of magnitude, moving from slightly to well above umf,j.This is particularly evident for density ratios close to unity.However, some values appear in contrast to the expected decreasing trend, especially at the lower velocities.Most of these correspond to cases for which a determination factor R 2 lower than 0.95 results from the fit, indicating that these systems follow an evolution not well represented by the simplified model assumed, and interpretation of the results needs to be done with caution.
It appears evident that a normalization of the val-ues reported in Table 5 would help to understand the mixing phenomena and better characterize the dependence on the velocity.However, in contrast to what has been done for the velocity, individuation of a reference value is not trivial.Further work addressing this aspect is in progress.

Conclusions
After introducing discrete element models of fluidized beds and discussing their properties, a DEMbased code was extensively utilized to investigate the role of gas velocity and density ratio on the solids mixing in fluidized beds of two components of different density.The mixtures considered, characterized by density ratios spanning 1 to 3, were fluidized at velocities higher than the minimum fluidization velocity of the jetsam and up to about twice this value.The results, in terms of steady-state mixing index and characteristic time, showed that full mixing can be unachievable for high density ratios, due to the presence of a jetsam-rich region at the bottom of the

Fig. 1
Fig. 1 Experimental (left) and simulated (right) flotsam volumetric concentration along the dimensionless bed height.The top row corresponds to low velocity and the bottom row to high velocity, as indicated.Simulation conditions are reported inTable 2.

Fig. 2
Fig.2Evolution of the degree of mixing from initially mixed (MI ≈ 1) and segregated (MI ≈ 0) configuration at the two dimensionless velocities indicated.Simulation conditions are reported in Table2.

Table 2
System data and parameters used in the simulations

Table 3
Dimensionless velocities for the combinations of density ratios and actual velocity investigated

Table 4
Steady-state mixing index values obtained for various density ratios and gas velocities

Table 5
3haracteristic mixing times in seconds.Bold values resulted from a fit with a determination factor R 2 < 0.95 From a dynamic point of view, it was shown that, depending on the gas velocity, the time necessary to move from full segregation to the equilibrium state can vary by more than one order of magnitude.financialsupport under the Project PRIN 2005 No. 61 (Area: 09).Prof.Brunello Formisani and Dr. Rossella Girimonte are gratefully acknowledged for providing the experimental data and the valuable discussions on mixing and segregation in fluidized beds.Poisson's ratio, dimensionless ρ particle density, kg m -3 ρf fluid density, kg m -3 σ 2 variance of the concentration distribution, dimensionless ψ deviatoric stress tensor, Pa τ characteristic time, s χ exponent in Eq. (T5), dimensionless Ω control volume size, m3