Non-spherical particles in a pseudo-2D fluidised bed: Modelling study

a r t i c l e i n f o


Introduction
Fluidised beds are used in biomass gasification, clean energy production and a variety of other applications such as coating, drying, granulation, food processing and gas phase polymerisation (Rapagna et al., 2000). They offer excellent solids mixing characteristics and a high specific contact area between the fluid and solid phase. As a result, fluidised beds possess high heat and mass transfer rates when compared with other modes of contacting. However, the complexity of these contactors makes their prediction and scale-up very difficult. As a consequence, numerous fundamental studies have been performed to better understand fluidisation, both experimentally and computationally. In general computational models require proper validation using experimental or analytical results. A major part of the modelling studies of fluidised beds has been carried out for spherical particles . However, particles are rarely perfectly spherical . Particles possess a variety of shapes ranging from rod-like, needle shaped to cubic and ellipsoidal. Simulations that approximate particles as spheres cannot precisely predict the behaviour of real, complex shaped particles as encountered in e.g. biomass processing. Non-spherical particles are known to produce intermittent flow and dilute packing fractions as compared to spheres. They tend to result in larger fluidised bed voidages and larger minimum fluidisation velocities due to interlocking of particles (Kodam et al., 2010). Information related to the influence of particle shape is therefore important, for example for reactor design and optimisation. Hence, there is a need to perform detailed particle-resolved simulations of fluidised beds consisting of non-spherical particles. A detailed insight into particle and gas motion in the bulk region of the bed can be obtained from simulations. The coupled Computational Fluid Dynamics (CFD) -Discrete Element Method (DEM) approach combines discrete particle tracking with continuum modelling of the fluid. This approach generates detailed transient information such as the trajectories and orientation of particles and forces and torques acting on individual particles, which is extremely difficult, if not impossible, to obtain experimentally. Such information is important to understand the fundamentals of particle-fluid interactions in dense flows. This approach has been extensively employed to simulate systems with particle-fluid interactions (Tsuji et al., 1993;Deen et al., 2007). However, this approach is still largely limited to spherical particles. Simulating non-spherical particles can be extremely complex and difficult for representing shape and detecting their contacts in the DEM system. Moreover, representation of the interaction between non-spherical particles and fluid increases complexity of data transfer between the CFD and DEM parts. At the same time, determining accurate fluid-particle interaction forces based on local orientation of particles while accounting for local voidage can be extremely challenging.
In the past 10 years, non-spherical particle fluidisation has rapidly gained interest as a field of study. Many researchers have studied non-spherical fluidised systems experimentally (Shao et al., 2013;Escudie et al., 2006;Liu et al., 2008). Even more so than experimental studies, numerical methods are employed to study non-spherical particle fluidisation. Zhong et al., investigated cylindrical particles approximated through clustered spheres with the CFD-DEM (Zhong et al., 2009). Ma et al. performed CFD-DEM simulation of rod particles in a fluidised bed using super-ellipsoids (Ma et al., 2016). They investigated the orientation effects for different aspect ratio and gas inlet velocities. Zhou et al. simulated ellipsoidal particles in a fluidised bed and demonstrated the effect of shape and particle aspect ratio on minimum fluidisation velocity of such particles (Zhou et al., 2011a,b;Gan et al., 2016). Vollmari and Kruggel-Emden conducted research on a number of irregularly shaped particles via experiments as well as simulations (Vollmari et al., 2016). Spherocylinders have been previously studied for the influence of particle aspect ratio on flow behaviour and packing (Langston et al., 2004;Zhao et al., 2012). Ren et al. performed simulations of spouting bed of rod like particles (Ren et al., 2014). Other studies of non-spherical gas-solid systems using CFD-DEM include pneumatic conveying (Hilton et al., 2010;Oschmann et al., 2014) and fixed beds .
In recent years, there is increased interest in the implementation of complex shaped particles . Lu et al. made a comprehensive summary of recent developments of the DEM, including a variety of methods used for simulating non-spherical particles (Lu et al., 2015). There is a large amount of literature on DEM methods that can realise arbitrarily-shaped particles while at the same time being fast, efficient, robust and accurate (Kodam et al., 2010;Pournin et al., 2005). A detailed overview of major DEM applications can be found in literature (Zhu et al., 2008). Among these methods to represent complex particle shapes the most popular is the glued sphere approach, whereby the particle is approximated by a certain number of spheres with overlap. A higher number of spheres results in a more accurate shape representation, but at the expense of computational time. A trade-off between robustness of the method in representing different shapes, accuracy of contact detection and overall efficiency has to be made.
An important aspect of CFD-DEM simulations is an accurate estimation of the drag force acting on the particle. A number of drag closures are available in literature for spherical particles. However, there are very few drag correlations available for nonspherical particles. Hölzer and Sommerfeld proposed a drag force model for single, complex shaped particles (Hölzer and Sommerfeld, 2008). In literature, this model is used in combination with a model of Di Felice (1994) to account for the effect of surrounding particles on the drag force. A number of researchers have followed this approach for simulating fluidised beds with CFD-DEM and demonstrated that the overall bed hydrodynamics can be fairly well predicted (Vollmari et al., 2016;Hilton et al., 2010). Zhou et al. (2011a) simulated ellipsoidal and oblate particles, Ren et al. (2014) simulated a spouted bed of corn-shaped particles, while Hilton et al. simulated a fluidised bed using superquadrics representing four types of complex shaped particles (Hilton and Cleary, 2011). A detailed summary of these works can be found in the paper of . Vollmari et al. validated the CFD-DEM model for the pressure drop and bed height. Cai et al. numerically and experimentally studied the orientations of cylindrical particles in fluidised beds (Cai et al., 2012). While there are a number of studies on computational modelling of non-spherical particle fluidisation, there are no detailed validation cases available in literature. In this study, we perform a detailed one-to-one comparison of simulations with our previously performed experiments . In a semi-2D fluidised bed filled with spherocylindrical particles, pressure drop and bed expansion is measured as a function of gas flow rate. Particle orientation, circulation patterns and local particle dynamics in the bed obtained from simulations are compared with experimental measurements obtained from Digital Image Analysis (DIA), Particle Image Velocimetry (PIV) and Particle Tracking Velocimetry (PTV). We also compare the results of stacking of particles and coordination number with the experiments. The non-spherical single particle drag model of Hölzer and Sommerfeld (2008) is compared with a DNS drag model for spherocylindrical particle developed in-house. We propose two new voidage correction models and compare results with the Di Felice model (Di Felice, 1994).

Model framework
The modelling done in this work is performed using a combined CFD-DEM method. A modified version of CFDEM code is used to perform simulations. The DEM framework, based on LIGGGHTS Ò (Kloss et al., 2012), handles all particle-particle and particle-wall interactions in a Lagrangian manner. Meanwhile, the CFD code based on OpenFOAM Ò (OpenCFD, 2004) and CFDEM Ò coupling (Goniva et al., 2012) solves the gas flow through the particle bed in an Eulerian way and couples this flow to the particles. It should be noted that the gas flow is resolved on a grid larger than the particle size. Using this combined simulation method and parallel computing, it is possible to gain accurate information on both particle movement and gas flow in a lab-scale fluidised bed, while still maintaining manageable computational times. Each of these models are described in detail below.

Discrete element method
DEM is used to describe the interaction of spherocylinders with rotational and translational degrees of freedom. A soft contact model, first introduced by Cundall and Strack to describe interactions between spherical granular particles, is used in this work (Cundall and Strack, 1979). The individual particles are tracked and their trajectories are numerically integrated over time. The local contact forces and torques develop when the adjacent particles spatially overlap. Consider a spherocylinder p in a dense gasfluidised reactor. The translational motion for spherocylinder p can be calculated by integrating where F ! p;c is the contact force due to collision and F ! p;f !p is the total interaction force exerted by the fluid phase on the particles which is further explained in Section 5. Gravity is accounted for in the body force F ! p;b . The rotational motion on particle p is calculated by integrating where I ! p ; x ! p and T ! p are the moment of inertia, angular velocity and torque for particle p, respectively.

Contact model
We focus on accurate contact resolution for spherocylinder particles, while keeping the computational load minimal. A spherocylindrical particle can be represented in the model via a glued sphere approach. However, that is a method based on approximations, which may introduce new errors itself (Kruggel-Emden et al., 2008). A brief summary of DEM contact detection methods from literature with associated advantages and disadvantages is given in the Table 1.
Most of these methods are applicable to variety of particle shapes. Although these methods are versatile, these implementations come at the cost of accuracy and computational time. We keep the level of complexity to a minimum and choose to perform exact analytical calculations for the contact detection of spherocylindrical particles. In order to resolve the contact between particles and between particles and walls, a linear spring-dashpot model with rolling friction is used. The simple force models like linear spring-dashpot can be substituted for more accurate force models based on Hertzian force and contact volume calculations, if bulk properties only are of interest as opposed to detailed contact information such as contact area or contact duration (Kumar et al., 2018). In the case of spherical particles, particles overlap when the distance between the particle centres is less than the sum of the particle radii. For spherocylindrical particles, the identification of contacts between particles, and the subsequent calculation of the overlap region is more complicated. Two adjacent spherocylinder particles are deemed to be overlapping once the distance between their central shafts is smaller than the sum of their radii. The only requirement for this contact is the determination of the closest distance points on the shafts of the two particles.
A modified spherocylinder contact detection algorithm originally developed by Pournin et al. for granular flows has been used in this study (Pournin et al., 2005). The shortest distance points on the shafts s ! i ; s ! j are found using an improved version of the algorithm described by Vega and Lago (1994). Fig. 1 shows an example of an overlapping contact between two spherocylinder particles P i and P j . For any particle, R is the characteristic radius or radius of the spherical part of the spherocylinder, r ! i is the centre of mass, L is the shaft length, u ! i is the orientation unit vector originating at r ! i and v ! i is the translational velocity. k and w are two arbitrary parameters for particles P i and P j which range in the interval ðÀ1; þ1Þ such that when both k and w are within range [-1,+1] the vector r ! in Eq. (5) connects the two finite rods. s ! i and s ! j are closest distance points and can be expressed as follows.  (Ren et al., 2014;Zhong et al., 2009) Easy implementation, high degree of versatility Dissipative and stiff contact, computationally expensive 3 Superquadrics or superellipsoids Delaney and Cleary, 2010) For symmetric geometric shapes Contact computations cumbersome for particles with sharp edges 4 Discrete function representation (DFR) (Williams and OConnor, 1995) More advanced version of superquadrics method w.r.t. particle shape modelling Large memory allotted for sorting particle surface points, contact calculations expensive 5 Probability based methods (Jin et al., 2011) Contact calculations simple Contacts not actual, Suitable for only for regular polygons 6 Method of potential particles (Houlsby, 2009) Versatile w.r.t. particle shapes Problem modelling tangential component during contact 7 GJK Algorithm (Wachs et al., 2012;Seelen et al., 2018) Versatile w.r.t. particle shapes Slow compared to exact analytical calculation 8 Exact analytical contact calculation (Kodam et al., 2010) Accurate and fast Specific to particular shape In other words, k and w values represent points on shafts for which the distance between two rods is minimum. For the sample contact illustrated in Fig. 1, the shortest distance between the par- where s ! i and s ! j are points on the central axes of P i and P j respectively, is given by is the vector pointing from the center of particle P i to the center of particle P j . k 0 and w 0 are the values minimizing Eq. (5) and are given by: An algorithm used for determining the closest distance points between the central axis of particles is given below: The positions ( r ! i ; r ! j ) and unit orientation vectors u ! i ; u ! j of two particles under consideration are known. The closest point between the infinite lines is calculated via Eqs. (6) and (7) After s ! i and s ! j are determined, the mid-point between these points is r ! c , the point of contact and the degree of overlap between the particles is expressed as d n . n ! ij and t ! ij are the normal and tangential unit vectors for the contact respectively. The contact force acting on a particle is the total of normal and tangential forces calculated as a result of overlap.
The normal contact force exerted on particle P i by particle P j is given by using a linear spring-dashpot model.
where k n is the normal spring constant, f n is the normal damping coefficient and v ! ij;n is the normal relative velocity between the particles at the location of the contact point. The tangential contact force is calculated from the Coulomb-type friction expression In this expression k t ; d t ; f t ; l and v ! ij;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 initial particle contact at time t c;0 . This expression represents the elastic tangential deformation of the particles since the onset of particle contact. Additional equations describing the contact model are given in Table 2.

Contact parameters
The contact parameters are set based on the assumption of a maximum 1% overlap between particles to avoid unrealistic behaviour. A set of steps is followed to determine the contact parame-

Table 2
Equations for the contact model.

# Parameter
Expression(s) ters using the equations given in Table 3. This procedure is given as follows: Set the values for coefficient of restitution (e n and e t ) obtained from experiments, assume d max % 0:01D p . Estimate maximum relative velocity dv max from characteristic particle velocities. k n;min is estimated based on the potential energy stored in a spring and the amount of kinetic energy lost. k n d 2 Estimate duration of a binary contact t c . k t is estimated with orientationally averaged moment of inertia.
Damping coefficients, g n and g t can be estimated. Orientationally averaged moment of inertia for spherocylinder is assumed to be an average over the moment of inertia along its three principle axes.
In order to accurately model the particle-particle collisions in the bed, the coefficients of restitution (normal and tangential), and coefficients of friction (static and dynamic) need to be determined. These coefficients are determined from binary collision and sliding experiments.

Coefficient of restitution
In a binary collision, based on center of mass movement, the coefficients of restitution are defined as follows (Schwager and Pöschel, 2007;Kharaz et al., 2001), where v denotes the relative velocity between the two particles.
As the coefficients of restitution depend strongly on material type, but only weakly on particle shape (Marhadi and Kinra, 2005), measurements were conducted using a flat plate and spherical particles of the Alumide material used in our fluidisation experiments. A spherical particle of volume equivalent diameter was dropped onto the flat plate from a height of 20 cm. The plate was fixed on a rotating axis to allow for angled collisions. The collision was captured using a pco.dimax HD + high-speed camera (approx. 1600 fps), as shown in Fig. 2. Particle displacement between frames was used to determine velocity. Multiple (15) perpendicular collisions were captured, as well as 5 collisions per angle from 10 to 80 degrees at increments of 10 degrees. Measurements at 0-60 degrees were used to calculate e n , 20-80 degrees were used for e t . Fig. 3 shows the results from the particle-plate collision experiments. The error bars represent standard deviation. These measurements reveal that the restitution coefficient is approximately invariant to changes in the angle of collision for the angles investigated. The mean values found are e n = 0.43 (r = 0.08) and e t = 0.76 (r = 0.10).

Coefficient of friction
The dynamic coefficient of friction can be determined from the coefficients of restitution, as described by Kharaz et al. (2001) and shown in Eq. (14), where l d is the dynamic coefficient of friction and h the angle of impact.
The static coefficient of friction is determined from a sliding experiment. A block of the alumide material is placed on a slab of the same material. The angle of the slab is then slowly increased until the block starts sliding. A schematic representation of this setup is shown in Fig. 4. The static coefficient of friction can then be calculated from the angle a at which the block starts sliding (Eq. (15), where l s is the static coefficient of friction, m is the mass of the block and g is gravitational acceleration). This experiment is repeated 50 times. 3 Normal damping coefficient (Pournin et al., 2005) f n ¼ À 2Meff tc ln en 4 Tangential spring constant (Pournin et al., 2005) p 2 þ ln ðetÞ 2 5 Tangential damping coefficient (Pournin et al., 2005)

Neighbour list
To improve the performance of the DEM model, a neighbour list is used. In the neighbour list, the indices of particles surrounding a certain particle are stored so that contact calculations are only made with reference to these neighbouring particles. This neighbour list is then re-used for multiple time steps until the maximum displace-ment of a system particle has exceeded some prescribed limit and a new list must be built.
In order to further improve efficiency, an Oriented Bounding Box (OBB) (Ericson, 2005) is used. The bounding volume used to detect particles close to each other depends on the orientation of both particles. In this study, a multi-sphere shaped bounding volume is applied. As shown in Fig. 6a, a compound shape consisting of 3 spheres tightly surrounds the spherocylinder. The distance between the satellite points and the centre of mass (COM) x sat is given by Eq. (16), which was derived from geometry. In addition to the sphere radius (Eq. (17)), a skin of radius R skin is used to guarantee the neighbour list can be used for multiple time steps. When the maximum particle displacement since the last neighbour list build surpasses R skin =2, a new list is created. A similar approach is used to build the particle-wall neighbour list as shown in Fig. 6b. Since the closest distance between a particle and a wall is always located at one of the tips of the particle, the spheres are only centred around the ends of the rod.

Flow solver
The fluid phase is described in an Eulerian manner by imposing a mesh of equal sized cells on the fluidised bed. The PISO algorithm (Versteeg and Malalasekera, 1995) is used to solve the phase continuity and momentum transport equation for incompressible, Newtonian flow. Eqs. (18) and (19) are known as Model A in literature recommended by Zhou et al. (2011a) for dense gas-solid flows. Fig. 5. Determination of the coefficient of friction, according to Kharaz et al. (2001). Fig. 6. Multi-sphere neighbourlist building method for (a) particle-particle contact and (b) particle wall-contact. xsat and R neigh are given by Eqs. (16) and (17) respectively. where is fluid volume fraction, q f is fluid density, u ! f is fluid velocity and s ! ! f is the stress tensor. f ! p!f represents the momentum exchange between the fluid and particle phases. The volumetric force acting on the fluid phase due to particles and is given as where F ! p;f !p is the force acting on an individual particle and n is the total number of particles in the CFD cell. The details of the transient solver with a PISO loop for pressure velocity coupling used in this work are given in the paper by Goniva et al. (2010).

Gas-particle coupling
The total fluid-particle interaction force F ! p;f !p includes the drag force F ! d and the pressure gradient force F ! p;rp ¼ ÀV p rp. Here V p is the particle volume and rp is the pressure gradient in the fluid phase. Often, the pressure gradient force is accounted in the drag closure and needs to be carefully checked before use (Tang et al., 2015;Vollmari et al., 2016). The most important contribution to the particle-fluid momentum exchange is given by the drag force depending on the local relative velocity between the fluid and the particle and granular volume fraction. Additionally, for nonspherical particles the orientation of the particle needs to be accounted for to accurately calculate drag. In view of these points, we present the drag models used in this work in the two following sub-sections: Section 5.1: Single particle drag models and Section 5.2: Voidage correction models.
5.1. Single particle drag models

Hölzer and Sommerfeld drag
Hölzer and Sommerfeld (2008) derived an equation describing the drag coefficient for a single non-spherical particle in a gas flow (Eq. (21)). Here / is the particle sphericity, v ! r is the relative velocity between the particle and the gas, A e is the cross-sectional area of the volume equivalent sphere and A p;tot is the surface area of particle.
This equation incorporates the orientation of the particle in the crosswise (/ ? ) and lengthwise sphericity (/ k ), given by Eqs. (23) and (24) respectively. These are calculated based on h, the angle between the gas velocity vector and the particle orientation vector otherwise referred to as the incident angle or angle of attack.
þ 0:42 10 0:4ðÀ ln /Þ 0:2 1 / ? ð21Þ  (Mema et al., 2017) and the Sanjeevi model (Sanjeevi et al., 2018). It can be seen that the correlation by Hölzer and Sommerfeld is in good correspondence with DNS results. It can also be seen that the particle orientation has large influence on the acting drag force and therefore fluidisation of such non-spherical particles.

Sanjeevi drag
A drag closure has been developed in-house using Lattice-Boltzmann (LB) simulations for a spherocylinder rod of aspect ratio 4 (Sanjeevi et al., 2018). A number of LB simulations were performed for different angles of attack and particle Reynolds number defined as Re ¼ j v ! r jd e q f =l f . The coefficients are given in Table 4.

Voidage correction models
The drag force acting on a single particle in a gas flow is given by Eq. (27). It is crucial to take into account the effect of local voidage to determine the actual drag force acting on the particle.

Di Felice model
Di Felice (1994) developed a correlation describing the effect of local void fraction on the drag force (Eqs. (28) and (29)). Although not developed specifically for non-spherical particles, the Di Felice model is widely used to account for voidage effects (Nan et al., 2016;Zhou et al., 2011a). Fig. 7. Single particle drag coefficient as a function of Reynolds number. The black, gray and hollow symbols represent 0 , 45 and 90 angle of attack respectively.
Here is the local void fraction around the particle. Gidaspow (1994) recommends the use of the Ergun equation (Ergun, 1952) in dense regions, as it is derived for a dense, packed bed. Despite the fact that it was derived for isotropic-shaped solids, this equation also accurately describes the pressure drop over a bed of non-spherical particles, as will be shown in Section 7.2. The drag force on a particle derived from the Ergun equation is given by Eq. (30). The smallest of the Di Felice and Ergun drag forces is used.
Here / is the particle sphericity, m f is the fluid kinematic viscosity and d e is the volume equivalent particle diameter.  (Tang et al., 2015;Tenneti et al., 2011). Extending the approach of Di Felice and in the absence of multi-particle drag data for non-spherical particles, we propose using these drag correlations to account for voidage effects as follows: Re p Þ is the voidage correction factor while F ! d0 is given by drag equations evaluated for a single particle. The drag force expressions are given in Table 5. Fig. 8 shows the voidage correction factor as a function of void fraction and Reynolds number. The Tang and Tenneti models show stronger dependence on voidage compared to the model of Di Felice. This drastic change in voidage dependence can have significant impact on overall fluidisation behaviour. Further, it can be seen that Di Felice model shows hardly any dependence on particle Reynolds number above 300. Rong et al. (2013) suggested an extension for Di Felice model. The voidage correction factor (gðÞ) for the extended Di Felice equation showed relatively similar behaviour to the original Di Felice model, when compared to other correction factors presented in this work. In order to avoid the redundancy of the results, the authors decided not to include the model by Rong et al. in the current investigation.
A comparison of the overall drag coefficient for different drag model combinations is shown in Fig. 9. The comparison is done for two different particle Reynolds number and a particle oriented at 45 relative to the flow. As discussed earlier, the Tang and Tenneti models with single particle drag are more sensitive to changes in Reynolds number and thus better suited for voidage corrections than the widely-used drag model of Di Felice model. Even then, it should be noted that one common limitation of all these models is that they have been developed for (mono-disperse) spheres, which are randomly arranged in space. To summarize, there are three models available which all take into account the voidage effect but which, on the other hand, do not consider non-sphericity of particles and heterogeneities with regard to particle arrangement or coalignment patterns.

Void fraction calculation
The drag force depends on the local void fraction at each point in the bed. To attain these values, a distributed void fraction calculation is used, meaning that the volume of a particle is assigned to not one, but multiple cells when it crosses cell boundaries. The voidage calculation is performed on cells of same size as the CFD grid. Each particle possesses 16 satellite points, placed evenly in the complete particle volume as shown in Fig. 10. Each cell containing such a satellite point is assigned a fraction of the particle volume, creating a more continuous void fraction field and preventing sudden jumps in local porosity.

Simulation setup
Simulations were performed in order to validate the CFD-DEM model for two fluidised beds (small and large), with their sizes given in Table 6. The simulations are done for three different drag models: Hölzer and Sommerfeld single particle drag with Di Felice model for voidage effects (HDF) Sanjeevi single particle drag with Tenneti model for voidage effects (STE) Sanjeevi single particle drag with Tang model for voidage effects (STA) A schematic diagram of the setup is shown in Fig. 11. The number of particles in the bed was calculated from the total bed mass measured experimentally . All model parameters are listed in Table 6. Initially, a packing was generated by dropping randomly oriented particles into the bed under the influence of gravity. The bed height was verified for a number of different random particle orientations and a de-fluidised bed. The calculated bed height was the same for all initialisation methods, thus confirming an unbiased initial configuration. Gas was then introduced uniformly from the bottom to fluidise the bed. Simulations were run with superficial gas velocities starting from 0.2 m/s with 0.2 m/s intervals. The bed was allowed to attain a quasisteady state (5 s), after which measurements were performed from 15 s. The boundary conditions used to solve the gas flow are given below.
Inlet: fixed inlet superficial gas velocity (u z ¼ U 0 ) Walls: no slip, no penetration (u x ¼ u y ¼ u z ¼ 0 m/s) Outlet: fixed outlet reference pressure (P ¼ 0 Pa)

Cell 1
Cell 2 Cell 3 Cell 4 Fig. 10. The particle volume is distributed among the cells according to the location of the 16 satellite points. Eight more points are located in the plane perpendicular to the shown cross-section.

Results and discussion
In this section, the results obtained from the simulation of the small and large fluidised beds are presented. The results of particle orientation, particle co-ordination number and particle coalignment are presented only for the small bed.

Fluidisation behaviour
The simulations are done using a wide range of superficial gas velocities from 0.2 m/s up to 2.6 m/s for the small bed and up to 3.2 m/s for the large bed. Over this range of gas velocities two regimes are observed: Packed bed: At flow rates below the minimum fluidisation velocity (U mf ), particles are stationary. The gas passes through the voids barely moving the particles. Bubbling bed: At flow rates well above the minimum fluidisation velocity, large pockets of gas move through the center of the bed and particles are thrown high in the freeboard region.
This behaviour is also observed in experiments as reported in our previous work . We do not observe a channelling regime in simulations. This is due to the inherent nature of the CFD-DEM modelling approach. In CFD-DEM simulations, the fluid cell size is usually 4 times larger than the volume equivalent particle diameter. This allows for accurate calculation of void fraction, but has the disadvantage that changes in local drag due to channel formation are not well captured. Fig. 12 shows an instantaneous snapshot of fluidised bed simulation for large particles at U 0 = 3.2 m/s.

Pressure drop and bed height
The pressure drop is a key parameter for the characterisation of the fluidisation behaviour of particles. The pressure drop found from simulations of the small and large beds are shown in Figs. 13 and 14, alongside the experimental results for the same systems. Qualitatively, the pressure drop curves match their experimental counterparts very well.
In Fig. 13, the pressure drop curve for the small bed for the Hölzer-Sommerfeld drag model shows a very good prediction. On the other hand, for the drag models of Tenneti and Tang, fluidisation is achieved at gas velocities lower than in experiments. This is due to the strong dependence on voidage of these drag models as demonstrated in Section 5. In Fig. 14, the pressure drop curve for the large bed for the Hölzer-Sommerfeld drag model shows an under-prediction, especially in the packed bed regime. This results in a higher minimum fluidisation velocity than observed in experiments. In Section 5, it was shown that the Hölzer and Sommerfeld drag model does not under-predict the single particle drag coefficient. Therefore, the most probable cause for this discrepancy is the conversion from single particle drag to multiparticle drag. Contrary to the results for the small bed, for the large bed the pressure drop curves with the Tang and Tenneti drag models show a better prediction than the Hölzer-Sommerfeld model. One reason for this could be the relative accuracy of these drag models for different particle Reynolds numbers.
The bed height found from simulations of the small and the large cases are shown in Figs. 15 and 16 respectively. The measurement of bed height was achieved using a similar method as in experiments; by finding the maximum gradient in particle density. Similar to the pressure drop curves, the bed height curves qualitatively  match their experimental counterparts very well. The most probable cause of quantitative mismatch is again the inaccurate representation of the effective drag acting on the particles.
The other possible causes for the bed height mismatch are related to the contact model and the parameters, and is mainly of concern in the fluidised regime. The chosen inter-particle collision model accounts for particle penetration (overlap) based on distance but not based on the penetration (overlap) volume. This can cause inaccuracies in determination of the local voidage. Secondly, even though the contact parameters used in the model are based on experimental observations (see Section 3.2), this does not guarantee that the particle contact is captured realistically by the model. The DEM contact model, which is based on particle overlap, is different in nature from the true particle-particle and particle-wall contact, where particles deform and have an inherent surface roughness. For this reason, it is necessary to carefully tune the contact parameters (both for particle-particle and particlewall contact) to equate the simulated bed height and experimentally obtained values. Additionally, the contact stiffness has been assumed to be constant while for spherocylinder particles, varying  stiffness with transitions should be used according to the local geometry at the point of contact (Kumar et al., 2018;Kidokoro et al., 2015).

Mass flux
A comparison between mass flux and particle velocity results from experiments and simulations of the small and large bed is shown in Figs. 17 and 18 respectively. The mass flux is extracted from simulation data by multiplying the solids volume fraction in each cell with the average velocity of all particles in that cell. This is expressed by Eq. (32), where / ! m is the local mass flux, v ! the particle velocity, q p the particle density and the local void fraction.
Experiments and simulations show a good match for both large and the small bed. However, it is clear that in simulations, particles are much more mobile than in experiments, noted by the higher  mass flux in both negative and positive direction. It is likely that this discrepancy is caused by the same reasons of sufficiently accurate drag and contact parameters mentioned earlier in Section 7.2. It was shown that both particle-particle contact parameters (especially coefficients of restitution and friction) (Goldschmidt et al., 2001;Reuge et al., 2008) and particle-wall (Li et al., 2010;Ye et al., 2005) contact parameters can greatly influence the dynamics and solids mixing in fluidised beds. When particles are able to slide against each other more freely, bubbling is more vigorous, particles are thrown higher and solids mixing is faster.

Particle orientation
Figs. 19 and 20 shows the probability density function (PDF) of the particle orientation from experiments and simulations respectively for the large setup at different flow rates. The PDF is given by Eq. (33), where N p ðaÞ is the number of particles at angle a. Only particles close to and parallel to the front wall have been considered, to be comparable with the experimental results in Fig. 19. Evidently, in simulations, particle do not align themselves with the gas flow at high flow rates. Particles remain predominantly horizontal, with a small peak emerging at an angle of 0 . When the PDF is parsed for different positions in the bed (Fig. 20  bottom), it is clear that this peak is caused by interaction with the side walls. Similar results have been obtained by Oschmann et al. (2014).
The difference between particle orientation distributions may, in part, be attributed to the differing hydrodynamic forces experi-enced by particles in experiments and simulations. While spherocylindrical particles are inherently subject to drag, lift and torque in the laboratory-scale reactor, we consider only hydrodynamic drag forces in our simulations. Hydrodynamic torque is often considered to be negligible, given that gas density is low and the particle mean free path between particle collisions is short (Radl et al., 2012;Bernard et al., 2016).This assumption certainly holds for domains that have a large particle volume fraction i.e. in close proximity to the inlet of the reactor. However, in the upper domains of a reactor where particle volume fraction is comparatively lower, the trajectory of particles is likely to be more susceptible to the influence of the gas phase, and thus a particle can experience significant hydrodynamic torque and lift. As torque facilitates the rotation of a particle depending on its initial orienta-tion to an incoming flow, the inclusion of torque in the simulations may more accurately capture particle orientations as observed in experiments. Fig. 19. Probability Density Function (PDF) of the particle orientation (as observed from the front wall) in the large setup at different flow rates (top) and different positions in the bed (bottom, U 0 = 1.9 U mf ). Angles of À90 and 90 degrees correspond with particles laying down horizontally, an angle of 0 degrees corresponds with particles standing up vertically .

Particle coalignment
Particle coalignment has been studied only for particles close to, and parallel to the front wall. Fig. 21 shows the result of the orientation autocorrelation in the large bed from experiments and the simulations. The orientation autocorrelation value drops at a distance much less the particle length. Comparison shows a good quantitative agreement between simulation and experiment. However, the curves for experiments show little difference as a function of distance from the side walls. Contrary to this, in simulations although very small, a clear distinction is see with respect to proximity to side walls, although this effect is very small. As discussed earlier, this might be due to the simplicity of our model i.e. neglecting hydrodynamic torque effects.

Particle coordination number
Since in simulations there is direct access to particle positions and velocities, it is straightforward to compare simulation data with results from experimentally determined particle tracking velocimetry measurements. Fig. 22 shows an example snapshot of the large bed with particles colored according to their coordination number (CN). Fig. 24 shows the average particle coordination number in a way analogous to the experimental results in Fig. 23. As in the experiments, the height of the bed is divided into sections with height equal to half the bed width. In order to obtain the 2D CN that is measured using PTV, the 3D CN is divided by the number of particle layers in the cut-off radius.
From comparison of the figures, it is clear that both experimental and simulation results are in good agreement. In the lower section of the bed velocities are small and particles are densely packed, as also shown by Fig. 24 (bottom) in the previous section. Higher in the bed, particles mainly move upwards in dense structures and rain down in dilute regions. The same observations were made in experiments (Fig. 23).

Conclusions and outlook
In this study, the fluidisation behaviour of rod-like particles has been investigated numerically with the CFD-DEM approach and compared with experiments. A description of the model was given, focussing on the extensions developed in this study. An exact representation of particle shape through analytical calculation is used instead of a comparatively slow and inaccurate multi-sphere approach. The fluid forces acting on particles are calculated using particle orientation, shape and local void fraction. The bed behaviour is investigated at different gas inlet flow rates.
Comparison of experimental and numerical pressure drop and bed height results has been carried out for three different drag models in two differently sized beds. The results show good qualitative agreement. However, when compared quantitatively, the results show that the particle drag force is under-predicted in case of the combined Hölzer and Di Felice drag model while overpredicted in the other two cases. This is associated with the voidage correction term for the drag coefficient, which is based on data for spherical particles in absence of this information for sphe- rocylindrical particles. A combination of single particle drag for orientation effects and multi-particle drag correlation of spherical particles for voidage correction clearly is not the most accurate approach. For more accurate results (also in the channelling regime), a multi-particle drag closure is needed. This drag closure should be derived from direct numerical simulations of same particle shape and size, in this case, spherocylinder particles aspect ratio 4 (Sanjeevi et al., 2018). More importantly, such simulations should consider voidage, relative particle spacing, mutual particle orientation and particle Reynolds number as parameters. In addition, it could also take into account channeling within a single CFD cell.
Comparison of simulations with PIV and PTV experimental results shows that qualitatively the solids circulation pattern is well captured by the model. Comparison of the average particle orientation at different flow rates shows the importance of the hydrodynamic torque. In experiments, particles align themselves along the flow direction at high gas velocities. In simulations, this effect is not observed as particles remain predominately horizontally oriented. It was also shown that this observation is not an effect of particle-wall collisions but rather depends on the particle-gas interaction. In literature, the hydrodynamic torque is generally regarded as of very little influence, as gas viscosity is low and the mean free path between particle collisions is very short. However, our results suggest that the incorporation of hydrodynamic torque is necessary for accurate modelling of nonspherical particles.
The effects of lift forces and hydrodynamic torque are excluded in this work. The recent work by Mema et al. (2017) has shown that these forces cannot be neglected for rod-like particles. Inclusion of these additional forces should be considered for better prediction of non-spherical fluidised bed hydrodynamics. Furthermore, with this validation study, this CFD-DEM model can further be used to study the rheological behaviour of non-spherical particle suspensions. The DEM model can be used to measure and study the developed particle stresses. With a comprehensive study considering all the Fig. 22. An example snapshot of particle coordination number, obtained from simulations in the large bed at U0 = 3.2 m/s. The colors indicate the coordination number for the respective particles, blue: particles with 9 or more neighbours, green: 5 to 8 neighbours and red: less than 5 neighbours. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) relevant parameters encountered for fluidised bed, a stress closure can be developed, which can be used in modelling particle stresses in a more coarse-grained model for industrial scale simulations.