Three-dimensional fluidized beds with rough spheres : validation of a two fluid model by magnetic particle tracking and discrete particle simulations

Two ﬂuid model simulations based on our recently introduced kinetic theory of granular ﬂow (KTGF) for rough spheres and rough walls, are validated for the ﬁrst time for full three-dimensional (3D) bub-bling ﬂuidized beds. The validation is performed by comparing with experimental data from Magnetic Particle Tracking and more detailed Discrete Particle Model simulations. The effect of adding a third dimension is investigated by comparing pseudo-2D and full 3D bubbling ﬂuidized beds containing


Introduction
Gas-solid fluidized beds are widely used in the chemical, petrochemical and metallurgical industries owing to their solids mobility and high heat and mass transfer rates resulting from intensive contact between the gas and the solid particles.To improve and scale up existing processes, it is of utmost importance to understand the hydrodynamics of such gas-solid fluidized beds.However, the collection of detailed experimental data is challenging, costly and becomes rather complicated for large scale threedimensional (3D) systems because of the lack of optical access.In addition, large scale computational simulations of fluidized beds are often limited to 2D or pseudo-2D systems because of the high computational costs associated with large 3D beds.Fortunately, with the increasing availability of computational power and more efficient numerical schemes, simulations of 3D fluidized beds have become possible.The aim of this paper is to compare prediction of a recently developed computational model for rough sphere fluidization with experiments on 3D dense bubbling fluidized beds, and to highlight differences between pseudo-2D and full 3D beds.
Eulerian-Lagrangian and fully Eulerian models are widely used to simulate gas-solid flows.In these models, the gas phase is described by the volume-averaged Navier-Stokes equations.In the Eulerian-Lagrangian Discrete Particle Model (DPM) (Tsuji et al., 1993;Hoomans et al., 1996;Xu and Yu, 1997) individual particles are tracked in the computational domain where the particle's motion is described by the Newtonian equations of motion.The DPM can account for direct particle-particle interactions in a fundamental and detailed manner.However, due to CPU constraints, only a limited number of particles (O(10 6 )) can be treated simultaneously.To reach larger scales, fully Eulerian models are preferred: the Two-Fluid Model (TFM) is better suited for the simulation of large scale gas-fluidized beds.In this approach the solid phase is treated as a second continuum, inter-penetrating with the continuous gas phase.Constitutive equations are solved using additional closure equations for the particle phase (Kuipers et al., 1992).This approach has emerged as a very promising tool because of its computational efficiency.The challenge here is to establish an accurate rheological description of the solid phase, which in most modern TFM simulations is based on kinetic theory of granular flow (KTGF).
The most widely used KTGF models (Ding and Gidaspow, 1990;Nieuwland et al., 1996) have been derived for dilute flows of slightly inelastic, frictionless spheres.In reality, however, granular materials are mostly frictional.The roughness of the granular materials has been reported to have a significant effect on stresses at least in the quasi-static regime (Sun and Sundaresan, 2011).Besides, according to Yang et al. (2017a), the particle surface friction has a strong effect on the solids flow patterns and distribution.From literature, attempts to quantify the friction effect have been somewhat limited.Yang et al. (2016a) derived a kinetic theory of granular flow for frictional spheres in dense systems which includes the effects of particle rotation and friction explicitly.Moreover, this theory has been validated by Yang et al. (2016b) for a pseudo-2D dense gas-solid bubbling fluidized bed.Experiments by Sommerfeld and Huber (1999) and numerical simulations by e.g.Lan et al. (2012) and Loha et al. (2013) have revealed the importance of wall boundary conditions (BCs) in determining the characteristics of granular flows.However, there is no consensus on the formulation of the BCs.Moreover, the physical meaning of the coefficients appearing in these formulations is not always clear.In fast granular flows, a rapid succession of almost instantaneous collisions between particles and a wall cause random fluctuations of the particle velocities, which determine the amount of momentum and fluctuation energy transferred through the walls (Louge, 1994).Yang et al. (2017b) derived new BCs for collisional granular flows of spheres at flat frictional walls.They characterized the influence of a frictional wall by physically measurable quantities: the normal and tangential restitution coefficients and a friction coefficient.Their framework accounts for both rotational and translational granular temperature.They performed simulations of a bubbling pseudo-2D fluidized bed using these new BCs, and showed that the new BCs are better capable of predicting solids axial velocity profiles and solids distribution near the walls.The most noticeable effect is a better agreement of the rotational granular temperature in comparison with that results obtained from DPM simulations.

Nomenclature
Even though numerical simulations can be used to generate a detailed understanding of flow structures in fluidized beds, validation of the models using advanced and detailed experiments are still crucial.Due to the opaqueness of dense 3D fluidized beds, non-optical techniques must be used like Electrical Resistance Tomography, Electrical Capacitance Tomography, Positron Emission Particle Tracking, Magnetic Resonance Imaging, and Magnetic Particle Tracking (MPT).In particular, MPT has emerged as a promising tool to investigate the hydrodynamics of fluidized beds due to its long-term stability, safety, and relatively low costs (Mohs et al., 2009).The method uses a magnetic tracer particle, which follows the bulk particle flow and whose magnetic field is continuously detected by multiple magnetic field sensors located outside the bed.This allows the reconstruction of the instantaneous position and orientation of the magnetic tracer.Based on statistical analysis of the tracer trajectory, characteristic measures of the bulk particle flow, such as the average particle velocity and particle circulation pattern, can be determined as a function of the operating conditions.The application of MPT in fluidization was initiated by Mohs et al. (2009) for the study of a spouted bed.Recently, MPT was improved by Buist et al. (2015) for use in dense granular flows prevailing in bubbling fluidized beds.
The bed geometry plays an important role with respect to the observed flow characteristics.Inspection of the open literature reveals that most combined experimental and simulation research is restricted to 2D or pseudo-2D systems.In such systems particle motion is very much confined by the front and back walls.The Two Fluid Model has been extensively used and validated in literature for such 2D systems (Cammarata et al., 2003;Lindborg et al., 2007).So far, the flow characteristics achieved for 3D fluidized beds were mostly compared with 2D or pseudo 2D systems, but the validity of such a comparison is questionable.It has been shown that 2D beds may yield significantly different results compared to 3D beds (Cranfield and Geldart, 1974).Therefore, special caution is needed when thin columns are utilized to investigate hydrodynamics of fluidized beds.Besides, there has been relatively little 3D numerical work (Wang et al., 2010;Verma et al., 2013Verma et al., , 2014;;Bakshi et al., 2015;Fede et al., 2016).
The present study focuses on validation of the KTGF model for rough spheres (Yang et al., 2016a) and the corresponding BCs for rough walls (Yang et al., 2017b), for the first time in a 3D setting.The validation will be done by comparing the model results with experimental MPT measurements in a 3D bed.It is important to make the right choice for the particles.Although the impact prop-erties of spherical particles of different materials are known (Lorenz et al., 1997), there are a number of limitations related to the choice of the magnetic tracer: the tracer necessarily has a high mass density and can interact with the bulk particles if they are paramagnetic.For this reason, for the bulk particles we used stainless steel 316 which is non-magnetisable.It also has a quite rough surface, making inclusion of friction effects essential.A systematic quantitative comparison between TFM simulations, DPM simulations and MPT experiments is carried out in both a pseudo-2D and a 3D bubbling fluidized bed with a square cross-section.We investigate the effect of different inlet gas velocities on the bed hydrodynamics, comparing time-and space-averaged quantities such as particle velocity, particle flux, particle circulation patterns and particle distribution.The aim of this comparison is to quantify the level of agreement between simulations and experiments concerning the particulate phase.Further, a careful comparison is made between the present model and the effective model by Jenkins and Zhang (2002) (which we label as "old TFM").
The paper is organized as follows.Section 2 gives a short overview of the numerical approach.The experimental setup is described in Section 3. Detailed model comparison and validation are presented in Section 4. Unexpectedly, the results obtained in the present work seem to indicate that, despite the reasonable agreement reached in the particle phase, the TFM model overpredicts the bubble size.Furthermore, results from both simulations and experiments confirm that there exists a significant difference between a 3D bed and a pseudo-2D bed.

Two fluid model
In the TFM both gas and particle phases are treated as fully inter-penetrating continua.The continuity equations excluding chemical reactions for both phases are given by @ @t where e; q, and v represent the local volume fraction, density and velocity, whereas the subscript k denotes the gas (k ¼ g) or solid (k ¼ s) phase.The corresponding momentum equations are given by @ e g q g v g where P is the hydrodynamic pressure, s is the shear stress tensor, and g is the gravitational acceleration.The terms e g rP g and e s rP g denote the buoyancy forces due to the gas pressure gradient.Note that closures for the solids stress tensor and the solids pressures are computed from the work of Yang et al. (2016a).Since the internal angular momentum changes of the particles are coupled to the particle flow field, the stress tensor contains anti-symmetric components which are associated with a rotational viscosity.This rotational viscosity is a consequence of the intrinsic spin of the particles: when particles collide with each other, the momentum transfer due to tangential collisions will be biased when the (average) particle velocity field is rotating, i.e. when there is a non-zero rate of rotation rv s À ðrv s Þ T .Meanwhile, the gas phase stress tensor is similar to the one for single phase fluid flow.Finally, the gassolid drag force, expressed through the interphase momentum transfer coefficient b A , plays a significant role in the fluidization process.Hence, we applied the combined Ergun (1952) and Wen and Yu (1966) where Re ¼ q g e g rjv g À v s j=l g represents the particle Reynolds number for a particle of diameter r moving through a gas of viscosity l g .
The KTGF for rough spheres by Yang et al. (2016a) leads to closure relations for the rheological properties of the particulate phase, which are explicitly expressed in terms of the friction coefficient.To describe the solids phase rheology thoroughly, an extra energy balance equation for the rotational granular temperature was derived.Thus, separate transport equations for the translational and rotational fluctuating kinetic energy need to be solved: The first term on the RHS of Eq. ( 6) is the production of translational fluctuating energy due to work done by the particle against the solids stress and pressure.The second term is the diffusion of fluctuating energy.The third term is the dissipation of fluctuating energy due to inelastic particle collisions.The last term represents the dissipation of energy by fluid drag, which also tends to reduce the particle velocity fluctuations.In the rotational energy fluctuation Eq. ( 7), only a diffusion term and a collisional source/sink term are present on the RHS.We note that a production term related to the average rotational velocity is not included here because in the development of the KTGF we assumed that the average rotational velocity of the particles is very small in the bulk, so we effectively set it to zero (the average rotational velocity near the walls can be non-zero as described next).Interested readers can refer to Yang et al. (2016a) for more details on the derivation of these closures.Also note that we have neglected the effects of rotational hydrodynamic torque for the large dense spherical particles investigated in this work.Buist et al. (2015) showed that this torque has a negligible effect on the rotational velocity profiles for spheres of the size and density used in this work.
In addition to the above bulk equations, boundary conditions (BCs) need to be specified for the gas and solid phase velocities and solid granular temperatures at the confining wall.For the gas phase, a no-slip wall boundary condition is used for all side walls (left, right, front and back side of the rectangular domain).A uniform gas velocity is specified at the bottom inlet, while atmospheric pressure (101,325 Pa) is prescribed at the top outlet.For the solid phase, a partial slip boundary condition is used for the side walls.We applied the relations derived by Yang et al. (2017b) for the solids velocity gradient and translational and rotational energy dissipation rate per unit area.In developing these boundary conditions, a distinction was made between sliding and sticking collisions, and particle rotation was included.A correlation was used, obtained from DPM simulations, that links the local solids velocity to the local average angular velocity of solid particles near the walls (Yang et al., 2017b).So contrary to the bulk regions, particles near the walls are allowed to have a non-zero average angular velocity.Depending on whether sliding or sticking motion is dominant, this leads to a sink or source term for the rotational granular temperature.
Specifically, the boundary conditions for the solids velocity at a flat frictional wall have the form where n is the direction perpendicular to the wall whereas v s;i represents the i-component of the local average solids velocity tangential to the wall.Similarly, V i is the i-component of the local average contact slip velocity.Note that the difference between the average contact slip velocity component V i and average (center-of-mass) velocity component v s;i is determined by the average angular velocity of the particles near the wall, for which we use the abovementioned correlation from Yang et al. (2017b).l w is the coefficient of friction with the wall and e w is the coefficient of normal restitution with the wall.g 0 is the radial distribution function, which we obtain from Ma and Ahmadi (1986).l t and l r are the shear and rotational contributions to the viscosity based on the KTGF of Yang et al. (2016a).A 1 and A 2 are functions of the three measurable wall collision parameters, the rotational and translational granular temperatures, and the slip velocity.They are expressed as where the abbreviations are, respectively, l w ð1þewÞ ; V t the magnitude of the tangential slip velocity, and k ¼ 5Hr 2Ht .Here e w ; b w , and l w are the normal coefficient of restitution, tangential coefficient of restitution, and friction coefficient, respectively, between a particle and the wall.
The boundary conditions for the fluxes of translational and rotational fluctuation energy at a flat frictional wall obey the following expressions: where j t and j r are the thermal conductivities from the KTGF of Yang et al. (2016a).e p is the coefficient of restitution for a binary particle-particle collision.B 1 and B 2 are expressed as We will compare the TFM simulation results based on this new rheology and new BCs with TFM simulation results based on an effective model for the particulate phase rheology (a simpler model that does not explicitly include particle rotation) (Jenkins and Zhang, 2002) which use partial slip BCs by Sinclair and Jackson (1989) with a specularity coefficient of 0.2.For all TFM simulations (both old and new) the frictional stress model by Srivastava and Sundaresan is used to take into account frictional stresses prevailing in the bulk at high solids volume fractions (Sinclair and Jackson, 1989).All simulation settings are specified in Table 1.

Discrete particle model
In the DPM, the gas phase is described in the same way as in TFM (Eqs.( 1) and ( 2)).However the solid phase is treated in much more detail.The motion of every particle is computed using the Newtonian equation of motion.The translational and rotational motion of each particle is described as: where m a ; r a , and v a are particle mass, position, and velocity, and x a , and T a are the particle moment of inertia, angular velocity, and torque around the center-of-mass of particle a.The external forces acting on particle a include gravity, buoyancy and drag forces: where V p is the particle volume.The gas velocity v g , gas fraction e g and interphase momentum transfer coefficient b A are all evaluated at the location of particle a by trilinear interpolation of the Eulerian fields.We adopt a three-dimensional soft sphere model from Cundall and Strack (1979) to calculate the total contact force on particle a.In this model, each particle in contact with particle a exerts a force which is decomposed into a normal and a tangential spring-dashpot force: where k n ; n ab ; d n ; d t ; g n ; g t ; v ab are, respectively, the spring stiffness in the normal direction, the normal unit vector, the overlap in the normal direction, the overlap vector in the tangential direction, the damping coefficients in the normal and tangential direction, and the relative velocity between the particle surfaces at the contact point.We refer to Deen et al. (2007) for details on the DPM model.Similar to the TFM model, in the present work we do not consider a hydrodynamic torque due to particle rotation in a fluid because this torque may be neglected for spheres of this size and density (Buist et al., 2015).

Single particle impact
As mentioned before, stainless steel 316 is used for the MPT measurements.Unfortunately, the collision parameters between spheres of this material and the plexiglass wall material are still unknown.Therefore, in this work, the particle-wall contact parameters were obtained experimentally following the procedure described in Gorham and Kharaz (2000).Fig. 1 illustrates the experimental setup which is used in this work.A spherical particle of stainless steel 316 (r ¼ 3 mm) falls under the influence of gravity onto an inclined, thick plexiglass plate (thickness = 20 mm).A high speed camera (LaVision Imager pro HS) with a frequency of 4000 frames per second is used to record the movement of the particle.Knowing the frame rate of the camera and the distance between two frames gives the velocity of the particle.The normal restitution coefficient e n is calculated as the ratio of the normal  component of the rebound velocity to the normal component of the velocity before the impact.Meanwhile, we obtained the effective (center-of-mass) tangential restitution coefficient e t at different inclinations of the plate (from 10 to 60 degrees) through the ratio of rebound tangential velocity to the impact tangential velocity.This coefficient can be related to the tangential restitution coefficient of the contact patch, b, as explained in Kharaz et al. (2001).However, we note that in our fluidized bed simulations, the value of b has a negligible influence on the bed hydrodynamics according to Zhao et al. (2015).The most important contact parameters are therefore the normal restitution coefficient and the friction coefficient.
The results are shown in Fig. 2. The solid squares represent values for the normal restitution coefficient.These values vary between 0.88 and 0.93 over a range of impact angles (see Fig. 1) from 0°to 60°with an overall average of 0.91.We therefore use e n ¼ 0:91 in all our simulations.The effective tangential restitution coefficient reaches a minimum at an impact angle of 30°, and then increases towards 1.These findings are similar to the ones reported by Kharaz et al. (2001) for their study of aluminium oxide spheres impacting on a soda-lime glass plate.They derived a relation between the effective tangential restitution coefficient, the normal restitution coefficient, coefficient of friction l, and the contact angle for a rigid body in the sliding regime, which is Fig. 3 shows the plot of the effective tangential restitution coefficient e t versus ð1 þ e n Þ cot h in the sliding regime (impact angle from 30 to 60 degrees).The slope of the nearly straight line represents the value of the wall friction coefficient; l w ¼ 0:12.However, we need to take into account that during the fluidization experiments, the particles and especially the wall surfaces could have   elastic deformation, which may increase surface roughness due to particle-particle and particle-wall collisions (Sommerfeld and Huber, 1999).Thus the coefficient of friction obtained in the impact experiments represents a minimum value; we have used 0.15 for the friction coefficient between particles and the walls in our simulations.

Magnetic Particle Tracking setup
Both a 3D fluidized bed and a pseudo-2D fluidized bed are used to validate the KTGF for bulk rheology and BCs for frictional spheres.
The 3D fluidized bed has dimensions of 0.15 Â 0.15 Â 1.0 m (W Â D Â H), which is sufficiently small to fit inside the rig embedding the 3D MPT sensor array.The four side walls of the fluidized bed are made of plexiglass.The porous distributor plate is made of bronze and has an average pore size of 10 lm and a thickness of 3 mm.Fig. 4 shows a snapshot of the 3D MPT setup, with the bed located at the center.The 3D sensor array has four rings of six tri-axis AMR sensors, resulting in 72 independent signals.The diameter of each sensor array ring is 0.26 m.During our experiments on the 3D fluidized bed we have used both an older 3D sensor array and a newer more accurate 3D sensor array, both provided by Matesy GmbH.We will compare the accuracy of the two sensor arrays in this section.
The pseudo-2D bed has the same height and width as the 3D bed, but a much smaller depth of 0.015 m.The same 2D-MPT as in Buist et al. (2015) is used here for measurements in the pseudo-2D bed.The distance between the measuring domain and the sensors is maintained at less than 2 cm during our experiments.
To ensure the collection (in a statistical sense) of sufficient data, experiments are carried out for 2.5-3.0 h.The averaged bed dynamics are inferred from the motion of the tracer.The principle of the MPT measurement technique can be found in Buist et al. (2015).We follow the same procedure to filter data and post-processing steps.An overview of all settings and properties is listed in Table 2.It is known that the experimental error of the marker position increases with decreasing magnetic moment of the marker.We determined this experimental error for our experimental setup by placing 5 magnets of different magnetic moment and size in the center of the sensor array.Using sequential quadratic programming (SQP) algorithms, the measurement error was determined from the standard deviation of the data positions in x-and y-direction.Fig. 5 shows that, for the same magnetic moment, the error is generally smaller in the newer 3D MPT setup.This implies that the new setup is more accurate, especially for the magnetic marker (I = 0.014 mA 2 ) used in our experiments, where the error is less than one particle diameter.

Pseudo-2D system
From the pressure drop in the pseudo-2D bed, we found that the minimum fluidization velocity is 3.91 m/s in the experiments.The minimum fluidization velocity is 3.75 m/s in the simulations using the Ergun/Wen and Yu (Ergun, 1952;Wen and Yu, 1966) drag law.This small mismatch between the minimum fluidization velocity in the experiment and the simulation is mainly due to the particle property (e.g.particle size distribution and density) and the experimental rig.To enable a fair comparison, in our analysis we adopt the same U g /U mf in both experiment and simulation, respectively.
The time-averaged solids volume fraction e s h i and solids flux patterns e s v s h i in the DPM and TFM simulations are compared in Fig. 6.For both runs, an averaging time of 25 s is used (after allowing 10 s to reach steady state), which is sufficiently large to obtain relatively stable results independent of the exact averaging time.
The DPM simulations show dense zones of solids close to the side walls and at the bottom of the bed.This type of solids volume fraction distribution reveals that bubbles are mostly formed at the bottom and move towards the center.On the other hand, animations of the porosity patterns indicate that in TFM simulations more bubbles than in DPM are generated at the bottom of the bed and larger bubbles are formed (due to coalescence) in the center of the bed.Therefore, the TFM simulations produce slightly more dilute zones in the lower part of the bed and more dilute zones in the center, and consequently larger dense zones near the side walls in comparison with DPM simulations.Besides, the very low solids concentration at the top of the bed from both DPM and TFM simulations indicates the bursting of bubbles.With increasing superficial gas velocity, in both DPM and TFM simulations bubble coalescence is enhanced and more pronounced vertical motion of bubbles occurs, leading to a more dilute zone in the center of the bed.
The TFM and DPM simulations show very similar solids flux patterns.The particles move laterally close to the distributor, flow upwards in regions of more intense bubble activity and downwards in regions of lesser bubble activity.Consequently, a pronounced global solids circulation pattern with two symmetric vortices in the middle of the bed is formed.The height of the dense zone and vertical size of the vortices grow upon increasing the superficial gas velocity.This was also observed by Lindborg et al. (2007).and the numerical simulations at U g ¼ 1:5U mf .Overall, the particles ascend in the center and descend near the side walls due to the preferred path of the rising bubbles.Note that rough wall BCs are employed also at the bottom wall in our TFM simulations, which probably hinders particle upwards motion close to the distributor.Additionally, in the dense bottom region, long-term and multi-particle collisions are dominant (Buist et al., 2016), which are not accounted for in the TFM simulations.As a consequence, the new TFM simulations underpredict the particle velocity in the center at the lower height of 0.05 m (Fig. 7a), but produce good agreement near the wall and in the annulus.At all other heights, a good match is obtained among the new TFM, DPM simulations and   MPT experiments in the bulk.In contrast, the old TFM obtains good agreement with the MPT experiments and DPM simulations in the lower part of the bed, but underpredicts the particle axial velocity in the upper parts of the bed.In the dense wall region, DPM overestimates the downward solid velocity.This deviation between MPT experiments and DPM simulation was also reported by Buist et al. (2015) for a bubbling bed.The present TFM simulations are in excellent agreement with the experiments.
Fig. 8 shows the particle axial solid flux e s v s;z at heights obtained from the experiment and numerical simulations at U g ¼ 1:5U mf .Compared to the previous figure, this also takes into account in the time average the amount of particles that move with a certain axial velocity at a certain location and time.At all heights, an upward solids flux is obtained in the center and a downward solids flux near the side walls.Consistent with our earlier observation, the upward solids flux close to the distributor is underestimated in the new TFM simulations as a consequence of the use of rough wall BCs at the bottom.Small bubbles are generated close to the distributor and side walls, which carry particles in their wakes, producing voids filled by downward flowing particles.Due to the coalescence of these bubbles, the amount of downward solids flux increases with increasing bed height, particularly in the near wall region.As bubbles move up in the bed center, particles flow vertically upward at higher axial area (0.05-0.15 m), which indicates an increase in upward solids flux.The number of particles close to the freeboard region is so limited that a lower upward solids flux is observed at the height of 0.2 m.All results corresponds well to the results in Fig. 6.Overall, the new TFM, the old TFM and DPM simulations are all in good agreement with the experiment.
Fig. 9 is the equivalent of Fig. 7, but for a higher superficial gas velocity U g ¼ 2:0U mf .In general, at these higher velocities, the axial velocities from the new TFM simulations agree well with the results from DPM and the experiment.Unfortunately, the old TFM based on Jenkins and Zhang (2002) underpredicts the particle axial velocity in the center at all heights, and underestimates the amount of particle slip with the side walls.A similar underestimation of the particle axial velocity is reported in the work of Lu and Gidaspow (2003).Particle friction leads to the formation of more heterogeneous structures because for rough particles more energy is dissipated during particle-particle and particle-wall collisions.Jenkins and Zhang (2002) introduced an effective restitution coefficient to account for the effect of particle friction on the energy dissipation in case the coefficient of friction is small.Even though their model is simple, it fails to predict the fluidized particle behaviour for these inelastic frictional spheres.Fig. 10 is the equivalent of Fig. 8, but for the higher superficial gas velocity of U g ¼ 2:0U mf .At all heights, the new TFM simulations obtain a better match with the DPM simulations than the old TFM simulations.The MPT predictions are sometimes deviating from the DPM simulations, especially at 0.1 m and 0.15 m height.It is important to note that measurements of the bulk particle flow in MPT are based on a statistical analysis of the marker trajectory.Consequently, the MPT experiment can only provide information on the product e s h i v s;z .In contrast, e s v s;z is measured in the simulations.It is clear that these measurements are not necessarily the same if there exist (positive or negative) correlations between fluctuations in e s and fluctuations in v s;z .Thus, for this part we should focus on the comparison among simulations.Taking the detailed DPM model to provide the ground truth, the new TFM achieves a significant improvement compared to the old TFM.
In summary, from Figs. 7 until 10 it is evident that the new TFM simulations agree better with DPM simulations and with MPT experiments (when referring to the axial velocity) than the old TFM at higher superficial gas velocity.This better agreement could be due to an overall decrease in solids volume fraction at higher gas velocity, associated with more binary collisions and more chaotic flow, making the assumptions of KTGF models better applicable.This confirms that the new TFM can be used to predict the effect of particle friction on gas-fluidized bed hydrodynamics.Fig. 11 shows the time-averaged translational and rotational granular temperature from DPM and TFM simulations at U g ¼ 2:0U mf .In the DPM simulations we observe that lower values of translational and rotational granular temperature exist near the walls, and higher granular temperatures in the center of the bed.Our TFM simulations confirm this semi-quantitatively, where the agreement is better for the translational temperature than for the rotational temperature.According to the work of Yang et al. (2017b), the flux of pseudothermal energy to the wall depends on the competition between a source term due to the transfer of kinetic energy from particle flow to velocity fluctuations induced by wall roughness, and a sink term representing dissipation by inelastic particle-wall collisions.In this case, the wall serves as a sink term.Also, as can be seen from Fig. 7, in most cases, the DPM simulations predict relatively larger particle velocity gradients close to the wall.Consequently, the flow shear rate is increased near the wall.All of these contribute to strong particle fluctuation motion, which explains the reason for higher granular temperature in that region.Fig. 11 (c) and (d) show that higher rotational granular temperature is obtained from DPM simulation.Further improvement for the rotational granular temperature prediction is needed.Fig. 12 shows instantaneous snapshots of the gas volume fraction (porosity) from DPM and both new and old TFM simulations for an inlet gas velocity of 1.5U mf .As indicated before, small bubbles are formed at the bottom of the bed, grow in size in the center due to coalescence, and burst through the bed surface.After this, the bed collapses rapidly and small new bubbles re-expand the bed again (Pain et al., 2001).Because of inclusion of particle friction and rotation, our new TFM model predicts a significantly more heterogeneous and stretched structure than the old TFM model, in agreement with the DPM simulations.Moreover, the fluidization in the new TFM simulations seems to be more of the slugging kind, also in agreement with the predictions from DPM simulations.
The time-averaged bubble size and bubble count are presented in Fig. 13 as a function of height at U g ¼ 1:5U mf .Note that the equivalent bubble diameter is evaluated using the bubble area A, i.e.D e ¼ ffiffiffiffiffiffiffiffiffiffiffi ffi 4A=p p .A porosity of 0.8 is used to define the bubble boundary.Moreover, we exclude bubbles in contact with the freeboard region to avoid ambiguity.Fig. 13(a) shows that the equivalent bubble size increases with increasing bed height, in agreement with Fig. 12.As shown in Fig. 13(b), a large number of bubbles appear near the bottom, and a deceasing number of bubbles with increasing bed height reveals bubble coalescence.The old TFM simulations predict larger and fewer bubbles than the DPM simulations and new TFM simulations.This is consistent with Fig. 12, where the particle friction (present in the DPM and new TFM simulations) leads to stronger heterogeneous structures, and therefore a higher number of bubbles.Fig. 14 shows the laterally averaged solids volume fraction as a function of height from simulations and experiment at U g ¼ 1:5U mf and U g ¼ 2:0U mf .For the different superficial gas velocities qualitatively similar contour shapes are obtained from the DPM and TFM simulations.The maximum solids concentration is observed at the bottom in DPM simulations.In TFM simulations and MPT experiment the solids concentration peaks at a height between 0.1 and 0.15 m.Close to the bed bottom (height<0.1),the particle concentration in TFM simulations is lower than that in DPM simulations, while in the height range 0.1-0.22 m the reverse is true.These observations correspond well to Fig. 6.At higher gas superficial gas velocity, bubble coalescence is accelerated, thus larger bubbles are formed through the bed center and higher solids concentrations are observed at lower superficial gas velocities.
Finally, in Table 3 we compare the particulate phase energy budgets for the TFM and DPM simulations as described in Yang et al. (2016b).Detailed expressions for the mechanical energy balance can be found in Goldschmidt et al. (2004).In the present work, we distinguish between translational kinetic energy E kin , rotational granular energy E rotgran , and gravitational potential energy E pot .The translational kinetic energy is further subdivided into a convective contribution E conv and translational granular energy E gran .Generally, from Table 3 we find that the energy distribution in the new TFM is closer to DPM simulations in comparison with the energy distribution obtained from the old TFM.Clearly, most of the energy is in the form of potential energy.The contributions of rotational energy in both models is small.It is evident that the amount of convective kinetic energy for the new TFM is in better agreement with the DPM data in comparison with the convective energy contribution in the old TFM to a large extent.

3D system
For the 3D bed, we found that the minimum fluidization velocity is 2.6 m/s in the experiments.The minimum fluidization velocity is 2.72 m/s in the simulations using the Ergun/Wen and Yu (Ergun, 1952;Wen and Yu, 1966) drag law.As for the 2D bed, this small difference between the minimum fluidization velocity in the experiment and the simulation exists due to the slightly polydisperse particle size distribution and small differences in the experimental rig.As a consequence of the significant friction with the front and back walls in the 2D bed, there exists a difference in the minimum fluidization velocity between a 2D and 3D bed.Here, we also adopted the same ratio U g /U mf to enable comparison with MPT experiments.
The time-averaged solids volume fraction and solids flux pattern in the DPM and TFM simulations are compared in Fig. 15.At low fluidization velocity of 1:5U mf , both models predict a rather homogeneous solids distribution, where a slightly higher solids concentration is observed in the DPM simulation.With increasing superficial gas velocity, bubble coalescence and lateral bubble movement become more pronounced in both the DPM and TFM simulations, which results in a low solids fraction near the central position at the bottom.When these results are compared with those from Fig. 6, it can be observed that particles are more homogeneously distributed in a 3D system in comparison with a pseudo-2D system.In the latter case, the frictional effects from the front and back walls substantially affect the movement of bubbles and particles.
Overall, our TFM and DPM simulations show very similar macroscopic flux patterns with an upward flux in the center and downward flux close to the side walls.At high fluidization velocity, the global solids circulation pattern becomes more pronounced.Fig. 16 shows the time-averaged particle axial flux at four crosssections above the distributor obtained from experiment and simulations, at a superficial gas velocity of U g ¼ 1:5U mf .In general, in the height range from 0.08 to 0.2 m, the DPM and (new) TFM models predict strong upward flux in the center and downward flux near the wall.As mentioned before, the rough wall BCs at the bottom influence the upward flux in the lower part of the bed, which is also observed in the bed center in Fig. 8 Fig. 17 is the equivalent of Fig. 16, but for a higher superficial gas velocity U g ¼ 2:3U mf .In the higher part of the bed (> 0:05 m), the new TFM simulations show a stronger good quantitative agreement with the DPM simulations and MPT experiments.However, the old TFM simulations underestimate the downward particle axial flux close to the wall.In the bulk, both TFM simulations predict a similar upward flux profile.At this higher gas velocity a stronger solids upward and downward flux is predicted in comparison with the one observed in Fig. 16 at the same height.In comparison with the pseudo-2D system, we found both TFM simulations obtain better agreement with the DPM simulations and MPT experiments, which implies that the particle surface friction exhibts less influence on bed hydrodynamics for the 3D system.Fig. 18 shows the time-averaged translational granular temperature from DPM and TFM simulations at U g ¼ 1:5U mf and U g ¼ 2:3U mf .On comparison with Fig. 15, it becomes clear that both DPM and TFM models produce a high granular temperature in the dilute areas with dilute particle concentration and visa versa.Furthermore, a higher fluidization velocity enhances the particle velocity fluctuations in both models.Fig. 15 shows that the TFM simulations predict a slightly more homogeneous particle distribution than the DPM simulations in the lower part of the bed, which leads to lower granular temperatures.Still, as mentioned in the work of Yang et al. (2017b), the new TFM improves the agreement with DPM simulations in simulations of rough spheres and rough walls.
Fig. 19 shows the laterally averaged solids volume fraction as a function of height from simulations and experiment at U g ¼ 1:5U mf and U g ¼ 2:3U mf .Most particles stay near the bottom of the bed due to gravity, which is in agreement with Fig. 15.At heights less than 0.12 m, the simulated particle concentration from TFM is lower than that from DPM simulation.Nevertheless, TFM simulations predict excellent agreement with the results from DPM simulation in the upper part of the bed (height > 0.12 m).It is interesting to notice that the agreement between TFM and DPM simulations becomes better in the full 3D system in comparison with the pseudo-2D system (Fig. 15).This is probably related to the fact that in the pseudo-2D system particle bridges can form between the front and back wall.The effects of such bridges are not explicitly modelled in the TFM simulations.
Fig. 20 shows snapshots of the three-dimensional bubble contours obtained from simulations at a superficial gas velocity of 1.5U mf .Bubbles are first generated close to the bottom of the bed, and move towards the bed center.Finally, these bubbles coalesce near the top of the bed.As a consequence, bubbles at the top are always larger than those close to the bottom of the bed (Pain et al., 2001).This is confirmed in Fig. 21, where we show the equivalent bubble diameter as a function of height (based on the volume of the bubble defined by a porosity contour value of 0.8).Fig. 21 reveals that the new TFM model slightly underpredicts the bubble size close to the bottom while it is reversed near the freeboard region.Nevertheless, the agreement between the TFM and DPM models is reasonable.Comparison between Figs.21 and 13 shows that there is a difference in the dependence of bubble diameter with increasing bed height between the 3D and pseudo-2D beds, especially close to the freeboard.This difference could arise from the additional freedom for bubble interaction in the 3D bed where bubbles are able to move freely in the horizontal plane instead of being restricted to essentially one direction (Geldart, 1970).

Conclusions
In this work we validated a new KTGF model for rough spheres and rough walls for the first time using a full 3D system.The model is validated by comparing TFM simulations that employ this KTGF model with MPT experiments and more detailed DPM simulations of the same dense gas-solid fluidized beds.To quantify the effect of the third dimension, we investigated both a pseudo-2D and a 3D configuration.
In the pseudo-2D bed, we compared the new TFM model with a TFM model based on a simpler kinetic theory derived by Jenkins and Zhang (2002).Through comparison with the DPM simulations and MPT experiments, we conclude that the new KTGF model improves the predictions for fluidized beds of inelastic rough particles.The energy distribution from the new TFM is closer to that predicted in the DPM simulations in comparison with the TFM based on Jenkins and Zhang's model.In particular, Jenkins and Zhang's model underpredicts the particle axial velocity in the bed center and predicts a lower particle velocity close to the side walls, especially at high superficial gas velocity.Because of particle friction, larger densely packed zones are obtained from simulations using DPM and the new TFM.
In a full 3D bed, the particles are more homogeneously distributed than in a pseudo-2D system.This is attributed to frictional effects from the front and back walls, which substantially affect the movement of bubbles and particles.The present TFM model obtains excellent agreement with experiment and discrete particle simulation for the prediction of particle axial flux.
We recommend that further validation of the new TFM be carried out for Geldart A and B type particles.Additionally, since industrial fluidized beds are generally large and cylindrical in shape, validation of the new TFM model in a cylindrical bed is also necessary.

Fig. 2 .
Fig. 2. Variation of coefficients of restitution with different impact angles for stainless steel 316 beads impacting on a plexiglass plate.

Fig. 3 .
Fig. 3. Measurements of effective (center-of-mass) tangential coefficient of restitution et versus ð1 þ enÞ cot h for stainless steel 316 beads impacting on a plexiglass plate.

Fig. 4 .
Fig. 4. Snapshot of the 3D fluidized bed in the MPT setup with 24 triaxis AMR sensors, arranged on 4 rings (at different vertical heights z) with 6 sensors each.

Fig. 5 .
Fig. 5. Experimental error on marker position as a function of magnetic moment for the old setup (black squares) and for the new setup (red circles).(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 9 .
Fig. 9. Time-averaged (10-35 s) particle axial velocity as a function of lateral (x) position, obtained from MPT experiments (black squares) and simulations (DPM: red circles; new TFM: blue up triangles; old TFM: magenta down triangles) of the pseudo-2D bed, at heights 0.05, 0.1, 0.15 and 0.2 m, for a superficial gas velocity Ug ¼ 2:0U mf .(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 10 .
Fig. 10.Time-averaged (10-35 s) axial solid flux as a function of lateral (x) position, obtained from MPT experiments (black squares) and simulations (DPM: red circles; new TFM: blue up triangles; old TFM: magenta down triangles) of the pseudo-2D bed, at heights 0.05, 0.1, 0.15 and 0.2 m, for a superficial gas velocity Ug ¼ 2:0U mf .(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 7
Fig.12.Instantaneous snapshots of the gas volume fraction (porosity) obtained from simulations for a superficial gas velocity of 1.5 U mf in the pseudo-2D bed.

Fig. 14 .
Fig. 14.Profiles of the time-averaged (10-35 s) and laterally averaged solids volume fraction versus height from simulations and experiment at two different superficial gas velocities in the pseudo-2D bed.

Fig. 15 .
Fig. 15.Time-averaged (10-35 s) solids flux patterns (arrows) and solids volume fraction distribution (colors, colorbar is for solids volume fraction) obtained from the central vertical plane of DPM simulations (left) and TFM simulations (right) of the 3D bed, at superficial gas velocities Ug ¼ 1:5U mf (top row) and Ug ¼ 2:3U mf (bottom row).(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 16 .
Fig. 16.Time-averaged (10-35 s) axial solid flux in the central vertical plane as a function of lateral (x) positions, obtained from MPT experiments (black squares) and simulations (DPM: red circles; new TFM: blue up triangles) of the 3D bed, at heights 0.1 m, 0.15 m, 0.2 m and 0.25 m, for a superficial gas velocity Ug ¼ 1:5U mf .(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 17
Fig. 17.Time-averaged (10-35 s) axial solid flux in the central vertical plane as a function of lateral (x) positions, obtained from MPT experiments (black squares) and simulations (DPM: red circles; new TFM: blue up triangles) of the 3D bed, at heights 0.05 m, 0.1 m, 0.15 m, 0.2 m and 0.25 m, for a superficial gas velocity Ug ¼ 2:3U mf .(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 17.Time-averaged (10-35 s) axial solid flux in the central vertical plane as a function of lateral (x) positions, obtained from MPT experiments (black squares) and simulations (DPM: red circles; new TFM: blue up triangles) of the 3D bed, at heights 0.05 m, 0.1 m, 0.15 m, 0.2 m and 0.25 m, for a superficial gas velocity Ug ¼ 2:3U mf .(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
(b).Nevertheless, the present two-fluid model simulations are in good agreement with the DPM simulation results at the other three heights.At the top of the bed, according to Laverman et al. (2008), the uncertainty in experimental measurement might be relatively high because the solids particles rain down from the top of the bubbles where the data accuracy is low.Comparison among Fig. 8(b)-(d) and Fig. 16(a)-(c) indicates that the magnitude of solids flux in the 3D system is smaller in comparison with the pseudo-2D system.

Fig. 19 .
Fig. 19.Profiles of the time-averaged (10-35 s) and laterally averaged solids volume fraction versus height from simulations and experiment at two different superficial gas velocities in the 3D bed.

Fig. 20 .
Fig. 20.Snapshots of three-dimensional bubble contours obtained from simulations with inlet gas velocity of 1.5 U mf .A porosity of 0.8 is used to define the bubble boundary.Top row: DPM simulations; bottom row: TFM simulations.

Table 1
Simulation settings.Note that in this work, x; y, and z coordinates are taken in the width, depth and height direction, respectively.

Table 3
Time averaged result of contributions to the mechanical energy for particles in the pseudo-2D bed (10-35 s).The total energy Etot consists of translational kinetic energy E kin , gravitational potential energy Epot and rotational granular energy Erotgran.The translational kinetic energy is further subdivided into a convective energy Econv , based on cellaveraged velocities, and a granular energy Egran, based on the particle fluctuating velocities within a cell.