Fluidization of spherocylindrical particles

. Multiphase (gas-solid) ﬂows are encountered in numerous industrial applications such as pharmaceutical, food, agricultural processing and energy generation. A coupled computational ﬂuid dynamics (CFD) and discrete element method (DEM) approach is a popular way to study such ﬂows at a particle scale. However, most of these studies deal with spherical particles while in reality, the particles are rarely spherical. The particle shape can have signiﬁcant e ﬀ ect on hydrodynamics in a ﬂuidized bed. Moreover, most studies in literature use inaccurate drag laws because accurate laws are not readily available. The drag force acting on a non-spherical particle can vary considerably with particle shape, orientation with the ﬂow, Reynolds number and packing fraction. In this work, the CFD-DEM approach is extended to model a laboratory scale ﬂuidized bed of sphe-rocylinder (rod-like) particles. These rod-like particles can be classiﬁed as Geldart D particles and have an aspect ratio of 4. Experiments are performed to study the particle ﬂow behavior in a quasi-2D ﬂuidized bed. Numerically obtained results for pressure drop and bed height are compared with experiments. The capability of CFD-DEM approach to e ﬃ ciently describe the global bed dynamics for ﬂuidized bed of rod-like particles is demonstrated.


Introduction
Non-spherical particles are frequently encountered in fluidisation processes such as in the production of biomass fuel in fluidised bed gasifiers [1].Prior to fludisation, biomass particles are milled, resulting in particle shapes such as cylinders with an aspect ratio between 3 and 5.When fluidised, these particles experience anisotropic drag and lift conditions due to the surrounding fluid [2] and anisotropic interactions with adjacent particles and boundaries.High aspect ratio particles experience much larger drag torques in comparison to spherical particles of the same volume.The altered drag and collision characteristics directly affect the voids and heterogeneous structures, as well as mass and heat transfer rates.As such it is imperative to account for orientation, alignment and wall proximity effects in a description of particle collisions and the formulation of solid-fluid drag relations.
A number of established methods are available to simulate fluidised beds comprising of spherical particles [3].At the macroscopic or industrial scale, the discrete bubble model (DBM) treats the gas or fluid phase in a Lagrangian manner while the solid particles are classified using an Eulerian approach [4].At the microscopic scale, molecular dynamics (MD) accounts for both the particles and the fluid discretely, representing a Lagrangian-Lagrangian model.Between these limiting cases, other algorithms provide a more balanced interpretation of particle and gas/fluid dynamics.For example, the two-fluid e-mail: j.t.padding@tudelft.nlmodel (TFM) approximates the gas/fluid and solid phases with generalised Navier-Stokes expressions and a description of solid phase rheological properties [3].On the other hand, a coupled computational fluid dynamics and discrete element method (CFD-DEM) approach describes the gas phase as a continuum while the particles are represented as discrete entities [5,6].The CFD-DEM approach has been extensively applied to study particle scale flows for spherical particle beds.This is an idealised view, however, and in reality particles are non-spherical in nature.In this work we present a CFD-DEM model to study a laboratory scale fluidised bed consisting of spherocylinder or rod-like particles.We implement an accurate contact detection approach for spherocylinders while employing an appropriate drag coefficient expression for non-spherical particles [7].We have also constructed an experimental apparatus to study laboratory scale fluidisation and compare experimental results on pressure drop and bed height with simulations.

Numerical Method
We have developed a CFD-DEM solver to simulate spherocylinder particles in fluidised beds that couples an Open-FOAM CFD solver and the DEM solver LIGGGHTS.The CFD-DEM model is developed on the Open Source coupling engine CFDEM, which allows for program execution for a user-prescribed number of time steps after which data is exchanged between the OpenFOAM fluid solver and the LIGGGHTS particle solver.

CFD
The fluid phase is described in an Eulerian manner by imposing a mesh of equal sized cells on the fluidised bed.The fluid velocity is spatially averaged at each cell.The motion of the fluid phase is calculated by solving the continuity equation (Eq.( 1)) and the momentum conservation expression (Eq.( 2)) where f is fluid volume fraction, ρ f is fluid density, u f is fluid velocity, τ f is the stress tensor, R f,p represents the momentum exchange between the fluid and particle phases and g is gravity.
As an exact drag force coefficient for the spherocylinders is unavailable, as an approximation, we use the drag force correlation formula for non-spherical particles from Hölzer and Sommerfeld [7] where Re, φ and φ ⊥ are the Reynolds number, sphericity and crosswise sphericity respectively.To account for the effects of surrounding particles, the drag force scheme of Di Felice [8] is used, which is expressed as where C D is the drag coefficient (Eq.( 3)), A ⊥ is the crosssectional area perpendicular to the flow, χ is a correction factor accounting for crowding effects, u f is the gas velocity interpolated to the location of particle i, and u i is the velocity of particle i.This expression is commonly used for spherical particles but has recently been applied to non-spherical particles [9].The correction factor χ is calculated as a function of the Reynolds number where d e is the equivalent diameter and η f is the fluid viscosity.A pressure-based solver using "Pressure Implicit Split Operator (PISO)" pressure-velocity coupling is used to solve the fluid equations, where an implicit momentum predictor followed by a series of pressure solutions and explicit velocity corrections is performed.

DEM
For the DEM component, spherocylinders are depicted as particles with rotational and translational degrees of freedom.Particle trajectories are integrated using a Velocity Verlet method based on local contact forces and torques.The motion of a particle p can be calculated from Newton's second law of motion with the expression where F p,n is the normal contact force, F p,t is the tangential contact force, F p, f is the drag force exerted by the fluid phase to the particles, F p,p denotes the pressure gradient (or buoyancy) acting on the particles.Gravity is incorporated in the body force F p,b .The angular momentum on particle p is calculated with I p dω p dt = T p where I p , ω p and T p are the moment of inertia, angular velocity and torque for particle p. Particle orientation is described by quaternions.
Figure 1: A contact between two spherocylinders.
The detection of contacts between spherocylinder particles is more involved than for the case of spherical particle contact detection.Figure 1 shows two spherocylindrical particles P 1 and P 2 experiencing a contact.For particle i, R is the characteristic radius, r i is the centre of mass, L is the shaft length, u i is the unit orientation vector, v i is the translational velocity and ω i is the angular velocity.The shortest distance between the particles is |s 2 − s 1 |, where s 2 and s 1 are points on the central axes of the respective particles.r c is the mid-point between s 2 and s 1 and δ n is the degree of overlap between the particles with t 12 the tangential unit vector and n 12 the normal unit vector.The normal contact force exerted on P 1 by P 2 is approximated as a linear dashpot given by where k n is the normal elastic constant, η n is the normal damping coefficient and v 12,n is the normal relative velocity.A Coulomb-type friction law is used for the tangential forces such that where k t is the tangential elastic constant, δ t is the tangential overlap, which is defined as the time integral of the tangential relative velocity since initial particle contact, η t is the tangential damping coefficient, μ is the friction coefficient and u 12,t is the tangential relative velocity.In simulations, the CFD time step Δt CFD = 1 x 10 -4 s, the CFD cell size is 0.015 m, the DEM time step Δt DEM = 0.1Δt CFD , μ = 0.3 and the number of particles is 32,400.
It is crucial that the fluid/solid void fraction is accurately calculated in the CDF-DEM model due to its strong influence on pressure drop in the particle bed.To enable calculation of the solid fraction contribution of any spherocylinder particle to grid cells, each particle is embedded with n sp evenly spaced satellite points throughout its volume.Each satellite point within a particle carries an equal volume weight.When the parent grid cell for a particle is identified, the particle volume is assigned to the parent cell or neighbouring cells subject to the position of the satellite points.This method is flexible and easily optimised through variation of n sp .

Setup
The focus of this study is particle behaviour near walls.Therefore we have constructed a quasi-2D fluidised bed setup of width 0.3 m, depth 0.05 m and height 1 m.The side and back walls are made of anodized aluminium while the front wall is made of 10 mm thick tempered glass.The distributor plate is made of sintered stainless steel with a pore size of 100 microns.A differential pressure sensor NXP MPX5050DP with a range from 0 to 50 kPa (0 to 72.5 psi) and 0.2 to 4.7 V output is installed just above the distributor plate to measure the pressure drop.A Dantec Dynamics FlowSense EO 16M camera is used to capture images for Digital Image Analysis (DIA).The camera has a resolution of 4872 x 3248 pixels, 4.2 frames/s at full resolution, inter-frame time of 300 ns and pixel size of 7.4 μm.Four LED lights, arranged in a 2 by 2 grid were used to illuminate the bed.The setup is shown in Fig. 2.During experiments the density of air is 1.1 kg/m 3 while the viscosity is 1.88 x 10 −5 Pa s.

Particles
The spherocylinder particles have been 3D printed using the selective laser sintering technique, a process that is cheap and highly precise.Made from a composite of aluminium and nylon dust, the particles have high strength and low density.Each particle has a diameter of 3 mm and length of 12 mm and have a density of 1395 kg/m 3 .Thus they can be classified as Geldart D type particles Prior to the experiments, the particles were tested for attrition.At rest the bed height is approximately 0.3 m.The total weight of the bed in experiment and simulation is the same.to extreme Geldart D particles i.e. with spherical diameter ≥ 3 mm.Notably the variation of ΔP here differs from spherical particles , where the pressure drop curve is smoother, since we observe the development of an initial plateau in the channeling regime followed by a step up to a second plateau in the BF regime.Our observation is in qualitative agreement with previous studies on rod-like particles [10].In addition we observe a hysteretic response in ΔP.We estimate the minimum fluidisation velocity U m f = 1.70 m/s, which marks the intersection of the maximum pressure drop or plateau in the BF regime with the ΔP curve for decreasing U 0 in the PC regime (Fig. 5).

Pressure Drop and Bed Height
Fig. 6 shows the variation of the bed height h bed with U 0 from experiments and simulations.The bed height is extracted from experimental images using a maximum gradient routine implemented using MATLAB and the Image Processing Toolbox TM .Over all regimes we find that h bed calculated from simulations is less than the experimental value.This difference in the trends might be attributed to a lower value of μ that is used in the simulations in comparison to the coefficient of friction of the aluminium/nylon composite particles in experiments.

Conclusions and Outlook
We have presented results on the variation of pressure drop and bed height from CFD-DEM simulations and experiments of the fluidisation of spherocylinder particles.With varying flow velocity U 0 , we have identified four distinct flow regimes.For low U 0 , we have a packed bed regime while for high U 0 we observe the onset of intense bubbling dynamics.For intermediate flow velocity, we observe passive channelling, which is succeeded by active channelling with increasing U 0 (Fig. 4).Simulation results recover to good agreement both the variation in pressure drop and bed height from experiments.
As part of the future work, we will include drag force closures specific to spherocylinder particles.These closures can be defined using DNS simulations for particles oriented at different angles, under the effect of walls, at different Reynolds number and particle volume fraction.We will determine the exact coefficient of friction and restitution for the particles by experimentally studying the interaction of a single particle with differing surfaces made from the particle composite material, glass and aluminium.Given its importance for calculation of the pressure drop, we will also revise the approach for estimating the fluid/solid void fraction.For comparative purposes, we will also study fluidized spherical particles both experimentally and numerically.

Figure 3 :
Figure 3: Variation of pressure drop with fluid velocity U 0 .Error bars represent the standard deviation.

Figure 4 :
Figure 4: Examples of the four flow regimes observed in experiments (as indicated in Fig.3).

Fig. 3
Fig.3shows the pressure drop ΔP for varying fluid velocity U 0 from experiments and simulations.In the experiments, pressure measurements were made at 5 m 3 /h intervals in gas flowrate by sampling at 100 Hz for 100 seconds for increasing and decreasing flowrates as indicated by the arrows and the trends in Fig.3.Four flow regimes have been identified, as indicated on Fig.3; namely packed bed (PB), passive channeling (PC), active channeling (AC) and bubbling fluidisation (BF).Snapshots of the bed in each regime are shown in Fig.4.In the case of spherical particles with equivalent particle volume, the channeling regimes are not observed.This behaviour is specific

Figure 5 :
Figure 5: Close-up of pressure drop curve showing calculation U m f = 1.70 m/s from experiments.

Figure 6 :
Figure 6: Variation of bed height (h bed ) with U 0 from experiments and simulations.